Introduction

Due to increasing human impact on the environment, more contaminants are reaching aquifers and causing deterioration of groundwater quality (Kleczkowski 1984; Appelo and Postma 1999; Dragon 2006; Bear and Cheng 2010; Witczak et al. 2013). It is necessary not only to investigate the processes of contaminant transport, but also forecast the fate of contaminants in the subsurface (Johnson et al. 2003). Forecasting contaminant transport in aquifers requires the formulation of an adequate mathematical model and identification of the model parameters. The mathematical description of such transport has been evolving, starting from the 1970s. Initially, sorption isotherm models were analyzed, which lead to the development of experimental methods focusing on these isotherms (Limousin et al. 2007). Nowadays, reactive transport models are essential to investigate related physical, chemical and biological processes (Lichtner et al. 1996). These models are associated with constitutive laws and assume conservation of energy, momentum and solute mass (Steefel et al. 2005). Lichtner et al. (1996) highlights that the continuum hypothesis underlies theoretical considerations concerning transport in aquifers. It is assumed that a physical system can be approximated by a mathematical model of continuous media and the properties of the system change smoothly enough to use calculus and related numerical tools to describe transport processes. Such an assumption makes it possible to analyze the dynamics and kinetics of the system using partial differential equations.

Transport processes can be analyzed on various time and spatial scales. The mathematical models of dynamically occurring processes can be verified during column experiments, using columns of various sizes (Liu et al. 2008). On the other hand, processes occurring on the geological time-scale, especially diffusion processes, are more difficult to verify (Steefel and Lichtner 1994). The upscaling of results from the column scale to the local or regional scale in the field is a separate problem. Bear and Cheng (2010) presented a synthetic overview of water balance modelling, groundwater flow and contaminant transport.

In order to describe the transport of pollutants in a sand deposit, it is necessary to recognize the hydrogeochemical conditions. The hydraulic conductivity and the hydraulic gradient are of most importance for the transport of pollutants in the sand deposit (Zhu 2019). The geological properties are determined by the lithology of the deposit and the content of clay fraction and organic materials. The chemical composition of the water as well as the type and concentration of the pollutant are also important factors (Appelo and Postma 1999; Fetter 2001).

The analysis of the hydrogeochemical conditions alone may not provide unambiguous information about the processes occurring in the case of transporting a specific pollutant in the investigated sand deposit (Dai et al. 2012). For the recognition of processes occurring during a column experiment, it may be helpful to analyze the recorded breakthrough curve of the reactive tracer and compare it with the recorded breakthrough curve of the conservative tracer.

Transport analyses are most often conducted in laboratory conditions, less often in the field. Laboratory column experiments are performed using conservative and reactive tracers (Fetter 1999; Małecki et al. 2006; Liu et al. 2008; Lin et al. 2019). As a result, breakthrough curves for both types of tracers are obtained. The selection of the transport model requires choosing an appropriate mathematical description of the processes taking place, particularly an adequate sorption model for the deposit, which will generate breakthrough curves for conservative and reactive tracers with a high goodness of fit. The current practice in the selection of the mathematical model of sorption relies on choosing a certain group of models—hereinafter this procedure will be called qualitative preselection. Next, the effectiveness of these models is compared in terms of accordance between the model predictions and the observation results, which is the basis for the proper selection. Methods of model preselection used so far are based on qualitative criteria such as the identification of physical and chemical mechanisms of sorption for the investigated deposit and the migrating substances (an adequacy criterion) and on the basis of the simplicity criterion (Limousin et al. 2007).

The proper selection involves assessing the goodness of fit between calculated and experimental breakthrough curves through measures or statistical tests (Saiers and Hornberger 1996). Koeppenkastropt and De Cario (1993) mention several adequate statistical indicators such as the standard deviation, coefficient of determination, correlation coefficient, and an index called the model selection criterion (MSC). The last one, besides characterizing the goodness of fit, also represents the information contained in a set of parameters of each model, and takes into account the ratio of the number of model parameters to the number of experimental data.

When selecting the best model, the parameters of all preselected models should be identified; furthermore, for each preselected model the goodness of fit should be determined. With this approach, selecting the appropriate model is laborious and time consuming, especially when dealing with poorly examined aquifers where the transport of unknown contaminants takes place and the existing physical, chemical and biological parameters are not recognized. Despite an extensive study, no other methods of selecting the adequate transport model have been found in the literature available to the authors. The complexity of transport processes and the variety of the proposed models make the selection of the transport model difficult, thus justifying the search for tools supporting the selection process.

The aim of the presented research was to develop an algorithm which, on the basis of quantitative measures (descriptors) of breakthrough curves, enables a significant narrowing down of the number of model groups that should undergo proper selection. The authors focused on analyzing impulse breakthrough curves for a reactive tracer in relation to the curves for a conservative tracer.

This new approach to the issue of preselection relies on the use of descriptors for both breakthrough curves, therefore it was called ‘quantitative preselection’. The development of the algorithm was based on a parametric analysis of breakthrough curves. The variability range of transport parameters adopted in the analysis corresponded to the variability of these parameters in natural hydrogeological conditions. In the initial stage of research, six descriptors which quantify changes in the shape of breakthrough curves were analyzed. The quantitative preselection algorithm of the transport model based on breakthrough curve descriptors was initially limited to a set of five simple sorption models (Okońska et al. 2019b) or six hybrid models (Okońska et al. 2019a). Ultimately, it turned out that for an effective quantitative preselection of a model from a set of five simple models and six hybrid models, it is sufficient to analyze three descriptors.

The algorithm proposed in this article facilitates and accelerates the selection of the best mathematical model of transport processes taking place during the column experiment. The algorithm also objectifies the model selection thanks to the use of breakthrough curve descriptors without the need to analyze the molecular mechanisms of sorption. It should be noted that the proper selection of a transport model, apart from statistical tests which determine theoretical and experimental curve fitting, must also include the expertise of a hydrogeologist, who always makes the final decision related to the model selection. The usefulness of the quantitative preselection algorithm of the transport model is illustrated on the example of a water intake located in the Prosna River valley, Poland.

Methods

Column experiment

The column experiment involves the flow of water and solutes through a column filled with sand (Domenico and Schwartz 1998; Liu et al. 2017). First, water flows in the experimental column until there is hydrochemical stabilization of the sand deposit, which can be manifested by the stabilization of electrical conductivity at the outflow; after that, water with tracers is injected. Usually, the simultaneous injection of a conservative tracer (only transported by means of advection, diffusion and dispersion) and a reactive tracer (additionally subject to sorption processes) is executed as an impulse with the duration of tin (Fig. 1; see Appendix for a summary of the nomenclature).

Fig. 1
figure 1

Schematic diagram of a column experiment

Full size image

The injection of tracers is described by the injection curve Cin(t), which shows the changes in tracer concentration over time at the input of the column (Fig. 1). Curves depicting the change of conservative and reactive tracer concentrations over time Cout(t), referred to as breakthrough curves, are recorded at the column output.

The dimensions of columns are significant in column experiments. It is important that the relation of the diameter D to the length L of the column provides a representative sample and minimal wall effects. Liu et al. (2008) conducted experiments using columns with the following dimensions: D = 3.4 cm, L = 10 cm to D = 15 cm, L = 80 cm. Analyzing the experiences of other authors and taking into account previous research conducted by the authors, columns with intermediate dimensions were used. In the presented example, a D = 9 cm, L = 50 cm column was used, which can be considered sufficiently representative for the investigation of transport processes. Moreover, the structure of the sand deposit is also important. It should be noted that analyzing undisturbed sand samples may seem unrealistic. On the other hand, a uniformly grained deposit prepared by selecting appropriate grains hardly reflects natural conditions; therefore, in the presented research, deposits containing grains of a natural range of sizes were tested.

Mathematical models

The transport of the conservative tracer in a porous medium involves processes such as advection, diffusion and dispersion. In addition, the reactive tracer is also subject to equilibrium or nonequilibrium sorption processes. Generally, one-dimensional (1D) transport of the tracer subject to transport processes is represented by the following equation (Luckner and Schestakow 1991; Weber Jr et al. 1991; Lichtner 1996):

$$ \frac{\partial C}{\partial t}=-\frac{k\cdotp i}{n_{\mathrm{e}}}\frac{\partial C}{\partial x}+\left({D}_{\mathrm{M}}+\alpha \frac{k\cdotp i}{n_{\mathrm{e}}}\right)\frac{\partial^2C}{{\partial x}^2}-\left(1-n\right)\frac{\rho_{\mathrm{s}}}{n_{\mathrm{e}}}\frac{\partial s}{\partial t} $$
(1)

where the following terms are as defined:

C:

Tracer concentration in the liquid phase [M L−3]

DM:

Coefficient of diffusion molecular [L2 T−1]

i:

Hydraulic gradient [L L−1]

k:

Hydraulic conductivity [L T−1]

n:

Total porosity [−]

ne:

Effective porosity [−]

s:

Tracer concentration in the solid phase (total sorbed-phase concentration) [M M−1]

t:

Time [T]

x:

Distance [L]

α:

Longitudinal dispersivity [L]

ρs:

Density of the porous medium [M L−3]

This equation does not involve complex chemical reactions with the soil, biodegradation, chemi- and electro-osmotic processes or capillary effects.

In the case of the conservative tracer, the last component of Eq. (1) is left out, as sorption does not occur (model A-D). In the case of the reactive tracer, the last component of Eq. (1) refers to sorption processes, whose mathematical descriptions may display various levels of complexity (Kaczmarek and Kazimierska-Drobny 2007). Simple transport models with equilibrium sorption include so-called sorption isotherms:

  • Linear Henry model (H model)

    $$ s={K}_{\mathrm{H}}C $$
    (2)

    where KH is the Henry distribution coefficient [L3 M−1]

  • Nonlinear Freundlich model (F model)

    $$ s={K}_{\mathrm{F}}{C}^{n_{\mathrm{F}}} $$
    (3)

    where KF is the Freundlich sorption coefficient [L3 M−1] and nF is the Freundlich sorption exponent [−].

  • Nonlinear Langmuir model (L model)

    $$ s=\frac{\alpha_{\mathrm{L}}{\beta}_{\mathrm{L}}C}{1+{\alpha}_{\mathrm{L}}C} $$
    (4)

    where αL is the Langmuir constant [L3 M−1] and βL is the total sorption capacity of the solid phase [M M−1].

The simple kinetic transport models with non-equilibrium sorption are:

  • Linear irreversible sorption model (I model)

    $$ \frac{\partial s}{\partial t}={k}_1C $$
    (5)

    where k1 is the irreversible sorption rate coefficient [L3 M−1 T−1]

  • Linear reversible sorption model (R model)

    $$ \frac{\partial s}{\partial t}={k}_2C-{k}_3s $$
    (6)

    where k2 is the first reversible sorption rate coefficient [L3 M−1 T−1] and k3 is the second reversible sorption rate coefficient [T−1].

The mathematical model describing 1D advection-dispersion transport with sorption is a system of two equations, i.e., the transport Eq. (1) and the sorption Eqs. (2)–(6) (Marciniak et al. 2009). In the case of equilibrium sorption, these would be Eq. (1) and one of the following: Eq. (2), Eq. (3) or Eq. (4). For nonequilibrium sorption, these would be Eqs. (1) and (5) or Eq. (6).

Sorption processes occurring during the transport of the reactive tracer through porous media in natural geological settings are often more complex than it appears from the simple mathematical description of equilibrium or nonequilibrium sorption (Liu et al. 1991; Johnson et al. 2003; Maraqa 2007; Yu et al. 2019). In the macroscopic sense, equilibrium sorption may occur on some parts of the surface of soil grains while nonequilibrium sorption may occur on the other part of the soil grain surface (Van Genuchten and Wierenga 1976). Tracer concentration in the solid phase throughout equilibrium sorption was denoted as se, and tracer concentration in the solid phase throughout nonequilibrium sorption as sn. In that case, instead of Eq. (1) combined with one of the simple sorption models, i.e. Eqs (2)–(6), a system of Eqs. (7), called the hybrid model, should be considered:

$$ \Big\{{\displaystyle \begin{array}{c}\frac{\partial C}{\partial t}=-\frac{k\cdotp i}{n_{\mathrm{e}}}\frac{\partial C}{\partial x}+\left({D}_{\mathrm{M}}+\alpha \frac{k\cdotp i}{n_{\mathrm{e}}}\right)\frac{\partial^2C}{\partial {x}^2}-\left(1-n\right)\frac{\rho_{\mathrm{s}}}{n_{\mathrm{e}}}\frac{\partial s}{\partial t}\\ {}s={s}_{\mathrm{e}}+{s}_{\mathrm{n}}\end{array}} $$
(7)

where the equilibrium component se and nonequilibrium component sn are described by one of the six hybrid models:

  • Henry and linear irreversible sorption model (H-I model)

    $$ \Big\{{\displaystyle \begin{array}{c}{s}_{\mathrm{e}}={K}_{\mathrm{H}}C\\ {}\frac{\partial {s}_{\mathrm{n}}}{\partial t}={k}_1C\end{array}} $$
    (8)

  • Freundlich and linear irreversible sorption model (F-I model)

    $$ \Big\{{\displaystyle \begin{array}{c}{s}_{\mathrm{e}}={K}_{\mathrm{F}}{C}^{n_{\mathrm{F}}}\\ {}\frac{\partial {s}_{\mathrm{n}}}{\partial t}={k}_1C\end{array}} $$
    (9)

  • Langmuir and linear irreversible sorption model (L-I model)

    $$ \Big\{{\displaystyle \begin{array}{c}{s}_{\mathrm{e}}=\frac{\alpha_{\mathrm{L}}{\beta}_{\mathrm{L}}C}{1+{\alpha}_{\mathrm{L}}C}\\ {}\frac{\partial {s}_{\mathrm{n}}}{\partial t}={k}_1C\end{array}} $$
    (10)

  • Henry and linear reversible sorption model (H-R model)

    $$ \Big\{{\displaystyle \begin{array}{c}{s}_{\mathrm{e}}={K}_{\mathrm{H}}C\\ {}\frac{\partial {s}_{\mathrm{n}}}{\partial t}={k}_2C-{k}_3{s}_{\mathrm{n}}\end{array}} $$
    (11)

  • Freundlich and linear reversible sorption model (F-R model)

    $$ \Big\{{\displaystyle \begin{array}{c}{s}_{\mathrm{e}}={K}_{\mathrm{F}}{C}^{n_{\mathrm{F}}}\\ {}\frac{\partial {s}_{\mathrm{n}}}{\partial t}={k}_2C-{k}_3{s}_{\mathrm{n}}\end{array}} $$
    (12)

  • Langmuir and linear reversible sorption model (L-R model)

    $$ \Big\{{\displaystyle \begin{array}{c}{s}_{\mathrm{e}}=\frac{\alpha_{\mathrm{L}}{\beta}_{\mathrm{L}}C}{1+{\alpha}_{\mathrm{L}}C}\\ {}\frac{\partial {s}_{\mathrm{n}}}{\partial t}={k}_2C-{k}_3{s}_{\mathrm{n}}\end{array}} $$
    (13)

To solve the transport equations, the authors adopted the impulse injection, for which appropriate boundary conditions were defined. The initial condition referred to the zero concentration of the substance in the whole analyzed volume at the initial time t = 0:

$$ C\left(x,t=0\right)=0\ \mathrm{for}\ x\ge 0 $$
(14)

At the column input, Dirichlet’s boundary condition was imposed in the form of an impulse injection with amplitude C0 and duration tin:

$$ \Big\{{\displaystyle \begin{array}{c}C\left(x=0,t\right)={C}_0\kern2.25em \mathrm{for}\kern2em 0<t\le {t}_{\mathrm{in}}\\ {}C\left(x=0,t\right)=0\kern2.75em \mathrm{for}\kern2.25em t>{t}_{\mathrm{in}}\end{array}} $$
(15)

At the column output, Neumann’s boundary condition was imposed in the following form:

$$ {\left.\frac{\partial C\left(x,t\right)}{\partial x}\right|}_{x=L}=0\ \mathrm{for}\ t\ge 0 $$
(16)

Descriptors

In order to investigate the impact of the changes in model parameter values on the shape of the breakthrough curves, the authors conducted a sensitivity analysis for all eleven models of pollutant transport in a saturated porous medium (Okońska et al. 2019a, b). The analysis was conducted for 12 parameters of water flow and tracer transport in the filtration column. The basic values of each parameter, as well as the range of their variability from the lower to the higher values were determined on the basis of a literature review (Table 1). The lower and higher values were chosen to correspond to the actual variability of these parameters occurring during the transport of both conservative and reactive tracers in natural hydrogeological settings.

Table 1 Variability range of hydrogeological parameters assumed for the transport processes, based on Okońska et al. (2019b)
Full size table

Six descriptors were used in the sensitivity analysis stage to quantify the breakthrough curves. When looking for a solution to the problem of transport model preselection, it was found that the number of descriptors can be limited to three. For the quantitative characterization of the shape changes of the breakthrough curves, the authors used the following three descriptors:

  • Relative tracer recovery ε [%]:

    $$ \varepsilon =\frac{100}{C_0{t}_{\mathrm{in}}}\underset{t_1}{\overset{t_2}{\int }}C(t) dt $$
    (17)

    where time t1 was determined from the breakthrough curve C(t) as the last moment at which C(t1) < 0.01·C0, whereas time t2 was determined from the breakthrough curve C(t) as the first moment at which C(t2) < 0.01·C0

  • Time of maximum tracer concentration at the output tmax [min]:

    $$ {\left.{t}_{\mathrm{max}}=t\right|}_{C={C}_{\mathrm{max}}} $$
    (18)

    where Cmax is the maximum tracer concentration at the output

  • Spread of breakthrough curve σ [min]:

    $$ \sigma =\sqrt{\frac{\int {\left(t-{t}_{\mathrm{max}}\right)}^2C(t) dt}{\int C(t) dt}} $$
    (19)

The meaning of all symbols present in definitions of descriptors for the conservative and reactive tracers are illustrated in Fig. 2. Detailed results showing the influence of model parameters on breakthrough curves, as well as descriptor values for the conservative and reactive tracers, were presented in Okońska et al. 2019a, b.

Fig. 2
figure 2

Symbols used in the description of descriptors: a the conservative tracer breakthrough curve; b the reactive tracer breakthrough curve

Full size image

Relative changes in the values of the three descriptors defined previously, determined for breakthrough curves of the reactive tracer and the conservative tracer, were calculated as percentage deviations δD:

$$ \delta D=\frac{D_{\mathrm{r}}-{D}_{\mathrm{c}}}{D_{\mathrm{c}}}100\% $$
(20)

where D is one of the descriptors ε, tmax and σ are defined by Eqs. (17)–(19), Dr is the descriptor value for the breakthrough curve of the reactive tracer, Dc is the descriptor value for the breakthrough curve of the conservative tracer.

Results

Algorithm of sorption model preselection

The results of the parametric analysis in the form of calculated percentage deviations δD are presented in tables as bar graphs for simple models (Fig. 3) and for hybrid models (Figs. 4 and 5). The third columns of the tables contain the lower, basic and higher parameter values. The percentage deviations of the descriptors in the range from −100 to 300% are depicted in columns 4, 5 and 6. Table rows include parameter values that correspond with particular models in the range provided in Table 1.

Fig. 3
figure 3

The variability range of descriptor deviations for simple models

Full size image
Fig. 4
figure 4

The variability range of descriptors deviations for hybrid models with irreversible sorption

Full size image
Fig. 5
figure 5

The variability range of descriptor deviations for hybrid models with reversible sorption

Full size image

For each model the descriptors whose percentage deviations in the conducted sensitivity analysis for all the considered parameters or combinations of sorption parameters exceed 5% highlighted. The assumed value of deviation results from the estimation of measurement uncertainties related with the execution of column experiments. Then, the baseline (minimum) δDb values were determined for percentage deviations of descriptors: δεb = 5%, δtmax b = 15% and δσb = 15%. The percentage deviation thresholds of the individual descriptors must be greater than the uncertainties associated with the measurement of the breakthrough curves. The proposed algorithm verifies if critical descriptor values were exceeded (in the following order: δεb, δtmax b and δσb). Since the values of descriptors and their deviations may be determined from the experimental breakthrough curves, it was assumed that the algorithm of model classification obtained as part of the sensitivity analysis can be applied to preselect the models of sorption for deposits and tracers occurring in natural settings. The resulting flowchart of the model preselection algorithm is presented in Fig. 6.

Fig. 6
figure 6

The algorithm of transport model selection based on the experimentally recorded breakthrough curve of the tracer. See text for symbol explanations

Full size image

Example of column experiment and algorithm application

The municipal water intake is located in the Holocene river terrace of the Prosna River (Poland) in the town of Tursko and supplies water to the city of Pleszew (Fig. 7). The aquifer has a limited spatial extent. The lithology is dominated by river and glacial sediments of varying thickness ranging from a few to even 60 m. The near-surface zone, 5–10 m thick, is dominated by Holocene peat, silt and river sands of various grain size, dispersed with organic matter. In the deeper part of the aquifer, Pleistocene sands and gravels are present. The aquifer is unconfined and highly susceptible to pollution from the surface.

Fig. 7
figure 7

Groundwater intake in the Prosna River valey

Full size image

The water intake consists of three pumping wells and six observation wells. Under natural conditions, groundwater flows from the south and west towards the Prosna River. Pumping intensified the seepage from local streams and drainage ditches. One of the drainage ditches is used to transport water from the nearby sewage treatment plant in the town of Gołuchów. The seepage of partially treated sewage from the drainage ditch into the aquifer deteriorates the quality of groundwater (Dragon et al. 2016).

While drilling the P3 observation well (Fig. 7), a sample of medium-grained sand was taken from a depth of 7–8 m, from the zone where contaminants seeping from the drainage ditch are transported. This sample was used to perform a column experiment in order to determine the parameters of water flow and transport of selected ions. The column experiment involved recording the breakthrough curve for chloride ions (Cl), treated as a conservative tracer, and breakthrough curves for selected reactive tracers, i.e., lithium ions (Li), zinc ions (Zn), and copper ions (Cu) migrating through a sample of medium-grained sand (Fig. 8). The concentration of chloride ions in the water sample was determined using the ion chromatography (IC) method. The detection limit for those ions was 0.011 mg dm−3. The concentration of lithium ions, zinc and copper was determined by means of flame atomic absorption spectrometry (FAAS). The FAAS method is widely used in research on the transport of pollutants in surface waters (Walna and Siepak 2012) and in groundwater (Ciążela et al. 2018). The detection limit for ions concentration was, respectively, 0.001 (Li), 0.002 (Zn) and 0.001 (Cu) mg dm−3.

Fig. 8
figure 8

Injection curve and experimental breakthrough curves for tracers: C/C0 – tracer concentration in relation to its concentration in the injected solution

Full size image

The uniformity coefficient of the sand collected from the aquifer at the Tursko water intake was U = 2.12 (d60 = 0.55 mm; d10 = 0.26 mm). The coefficient U is defined as: U = d60/d10, where d10 is the grain diameter [mm] corresponding to 10% of the sample passing by weight (which means that 10% of the particles are smaller than the diameter d10) and d60 is the grain diameter [mm] at 60% passing (Holtz and Kovacs 1981). Total organic carbon (TOC) content was 0.021%. Trace amounts (<0.01%) of the clay fraction were detected. Clay minerals identified in the sand sample were smectite, kaolinite and illite. The experimental column was filled with sand, taking care of saturation and the even distribution of the grains. For 2 days, filtration in column with water taken from the P3 observation well was carried out in order to obtain hydrochemical stabilization of the sand deposit. The hydraulic conductivity of the examined deposit was determined using the Darcy formula, by measuring the water flow rates several times. The total porosity was determined using the volumetric method. At the input of the column, the solution was injected as an impulse for 115 min. At the output of the column, samples of water were collected in chemically neutral vessels for 24 h. In this way, the breakthrough curves for each analyzed ion were obtained. These curves are shown in Fig. 8, as well as in Fig. 9a–d More details about the methodology of implementing the column experiment can be found in the works of Marciniak et al. (2009) and Okońska et al. (2017).

Fig. 9
figure 9

Fit of numerical breakthrough curves to experimental curves for the following ions: a chloride; b lithium; c copper; d zinc

Full size image

To calculate the values of descriptors which characterize the shape of experimental breakthrough curves, the authors used the computation code Descriptors2.m written in the MATLAB environment (Okońska et al. 2019c). Then, percentage deviations of the descriptors were determined (Table 2). Next, in accordance with the proposed algorithm (Fig. 6, δεb = 5%, δtmax b = 15% and δσb = 15%), mathematical models which best described tracer transport through the deposit were chosen for each analyzed ion. The previously analyzed baseline values of descriptor deviations δDt were taken into consideration. An example of using the proposed algorithm to select the model of copper ion transport is shown in Fig. 10.

Table 2 Values and percentage deviations of descriptors and selected transport models of specific ions
Full size table
Fig. 10
figure 10

An example of Cu ion transport model selection

Full size image

The parameters of transport were identified using the Id_Test_Kol4.m code in the MATLAB environment. First, transport parameters were determined for the advection-dispersion model based on the breakthrough curve of chloride ions (Table 3). Then, the parameters that govern the transport of reactive ions were calculated for mathematical models chosen by the algorithm. These parameters are presented in Table 4. To identify the parameters, the authors used the numerical solutions of transport equations Eqs. (1)–(13) for the assumed boundary conditions (Eqs. 14–16) with the use of optimization methods implemented in the MATLAB environment (Okońska et al. 2017).

Table 3 Identified parameters of the advection-dispersion model
Full size table
Table 4 Values of parameters determined for the analyzed transport. The models with the best fit are highlighted with an asterisk (*); the model parameters indicated by the proposed algorithm (based on MATLAB) are italic
Full size table

Matching of the numerically calculated breakthrough curves with the experimental curves for the models selected by the algorithm is presented in Fig. 9. The goodness-of-fit between numerical and experimental curves was established by calculating the root mean square error (RMSE) and the correlation coefficient r (Małecki et al. 2017).

Discussion

In order to check whether the proposed algorithm gives the best fitting transport model, the authors used the code Id_Test_Kol4.m to determine the parameter values of all sorption models considered in section ‘Mathematical models’ for the three reactive ions: lithium, zinc and copper. The five simple and the six hybrid mathematical models of sorption processes were tested. The results obtained for all the models are presented in Table 4. In the table, the models with the best fit, i.e., those with the lowest value of RMSE and the highest correlation coefficient r value, are highlighted with an asterisk. In addition, the model parameters indicated by the proposed algorithm are italic.

The conducted research resulted in the following findings:

  • For lithium (Li) ions the preselection algorithm proposed simple H and F models. The best fit between numerical and experimental breakthrough curves was obtained for the F, F-I and F-R models. The analysis of the identified parameter values shows that the differences between the breakthrough curves obtained from the F, F-I and F-R models are insignificant. In the hybrid models, the values of parameters k1, k2 and k3 are very low. This shows that equilibrium sorption is dominant in the process of lithium ion transport; therefore, the F model describes the process of lithium ions sorption with sufficient accuracy.

  • For copper (Cu) ions and zinc (Zn) ions, the algorithm proposed a simple I model and hybrid H-I and F-I models. For these ions, both components of equilibrium and nonequilibrium sorption are significant. Goodness-of-fit indices, RMSE and r, prefer the F-I and F-R models. However, the dominance of irreversible sorption is indicated by descriptor ε, which suggests that the recovery rate of copper and zinc ions is 22.3 and 32.4%, respectively.

Many authors have conducted research on the transport of Cu and Zn ions in groundwater. Krogulec and Zabłocki (2015) and Małecki et al. (2017) conducted research in the Kampinos National Park in Poland. They showed that the mobility of copper and zinc in near-surface groundwater is a function of the lithology of the hypergenic zone. Both ions sorbed irreversibly, the process being more intense in the case of copper ions. Karathanasis (1999) conducted column experiments with soil colloids and showed that the presence of colloids typically enhanced metal transport by 5- to 50-fold over that of the control treatments, with Zn being consistently more mobile than Cu. Arias et al. (2005) studied the adsorption and desorption of Cu and Zn in the surface layers of acidic soils. Both for copper and zinc ions, better fit was obtained in the case of the Freundlich equation than the Langmuir equation. Arias et al. (2005) only focused on simple sorption models; however, they found higher adsorption than desorption both in arable soils and forest soils, which would suggest the conditions of irreversible sorption and confirm that the selection of F-I model was right.

Studies by various authors confirm the occurrence of irreversible sorption of Cu and Zn ions, which was also noticed in the area of the Tursko water intake. Therefore, the F-I model can be considered as aptly indicated by the proposed preselection algorithm.

Summary and conclusions

In order to investigate transport in a sand deposit during a column experiment, a qualitative preselection of sorption models was conducted, taking into consideration sorption mechanisms and model simplicity. As a result, a set of 11 models was obtained. For this set, a new algorithm that implements quantitative preselection has been proposed, which, in the authors’ opinion, should precede and facilitate the proper selection of the best transport model.

A previously performed detailed parametric analysis and partial preselection algorithms that incorporate descriptors characterizing the geometry of breakthrough curves for conservative and reactive tracers served as the foundations of the proposed algorithm. Initially, a preselection algorithm was proposed for the five simple models (Okońska et al. 2019a), and then for the six hybrid models (Okońska et al. 2019b). In this study, the sorption model preselection algorithm was generalized with the joint consideration of simple and hybrid models using three breakthrough curve descriptors.

The starting points for the application of the algorithm are the experimentally recorded breakthrough curves of the conservative tracer and the reactive tracer. Next, the values of three descriptors D and percentage deviations of descriptors δD of the breakthrough curve of the reactive tracer in relation to the breakthrough curve of the conservative tracer should be calculated in accordance with the equations provided in this paper. Having the results of these calculations, the procedure is shown in Fig. 6. This algorithm will allow selection of the best fitting transport model. The usability of the proposed algorithm was confirmed by examples, one of which is presented in the article.

It must be stressed that the parametric analysis used as the basis for developing the proposed algorithm was conducted for a range of flow and transport parameters established by studying relevant existing literature. The issue that is worth exploring in further studies is the extension of the range of transport parameters adopted for the parametric analysis. This may result in the modification of the proposed baseline values of deviations of particular descriptors; furthermore, the conducted analysis refers to homogeneous media. A high diversification of the deposit grain sizes may be a factor disturbing the effectiveness of the algorithm. Broader application of the algorithm should show the influence of factors such as: the credibility of data derived from hydrogeological and hydrochemical exploration, the degree of hydrochemical complication of the transport process, soil gradation or the height of the deposit in the laboratory column.

It was demonstrated that the proposed algorithm can facilitate and accelerate the preselection of the best fitting transport model. The process of the proper model selection should be finalized by making the ultimate decision concerning the model choice. Such decision should always be made by the hydrogeologist interpreting the results of column experiments, who needs to take into consideration the hydrogeochemical determinants of contaminant transport in an investigated aquifer.

In the authors’ opinion, the developed methodology of creating the quantitative preselection algorithm can be used for other sets of models (sorption models, but not only), other materials (with a different composition, soil structure or the type of liquid) and a different scale of transport processes. A broader implementation of the algorithm could improve the reliability of its predictions. Therefore, the authors provide the Descriptors2.m code in the electronic supplementary material (ESM), which enables the calculation of breakthrough curve descriptors.