Subsurface transport of a sorbing contaminant is poorly understood and characterized. Here, a new semi-analytical saturated–unsaturated flow and transport model is coupled to a kinetic sorption algorithm to assess the impact of changes in the subsurface permeability architecture and flow rate on sorption characteristics. The model outputs reveal the pronounced effect of the rate of vertical decline in Ks on the frequency of occurrence and spatial distribution of subsurface sorption as well as the timing and rate of sorbing contaminants discharged into stream. Sorption potential is weakened with infiltration rate. The impact of infiltration rate on the decline in sorption potential becomes more accentuated as the degree of subsurface vertical heterogeneity in saturated hydraulic conductivity increases. Porosity pattern also impacts sorption characteristics; but its effects highly depend upon the degree of vertical heterogeneity in Ks. The results and methodology presented in this paper have potential implications for assessing water quality in integrated groundwater–surface water systems as well as designing remediation systems.

INTRODUCTION

Sorption to soil minerals is one of the most influential factors on the transport and fate of sorbing contaminants (e.g., phosphorus and pesticides) and thus water quality and pollution in integrated groundwater–surface water systems. The sorption of a contaminant to soil minerals changes the chemical properties of the soil solid phase and reduces the contaminant's level in the adjacent water column, and thus decreases its exposure and transport to other parts of the ecosystem. Contaminant sorption to soil can also affect the base flow chemistry of regional surface water bodies.

Subsurface hydrology and physical heterogeneity can impact sorption kinetics (Elfeki et al. 1997; Stenemo & Jarvis 2007; Fuchs et al. 2009). Certain local subsurface conditions can enhance the subsurface transport of flow and sorbing contaminants from the ground surface to surface water bodies (e.g., Ameli et al. 2015). Many field-based research studies depict the influence of soil heterogeneity creating preferential and non-preferential pathways, with lower sorption and a higher rate of phosphorus transport in preferential pathways compared to non-preferential pathways (e.g., Carlyle & Hill 2001; McCarty & Angier 2001; Fuchs et al. 2009). A laboratory experiment by Fuchs et al. (2009) also suggested that phosphorus sorption increases as subsurface flow velocity decreases.

Linking subsurface hydrology and heterogeneity to the concentration of sorbing contaminants in the riparian zone and stream is a big challenge in critical zone science (Cai et al. 2016; Surinaidu 2016). This is particularly important as understanding of sorption mechanism and consequences in different environments is essential for developing solutions to critical contamination issues, such as those associated with a global increase in phosphorus use reported by Tilman et al. (2001) and Vitousek et al. (2009). A basic but important step in this undertaking is the exploration of how sorption varies with subsurface permeability architecture and subsurface flow rate. Despite the emergent impacts of subsurface permeability pattern and flow velocity on the fate and transport of sorbing contaminants, to the author's best knowledge, there is currently no research study that adequately and ‘systematically’ assesses these impacts.

Typically, soil column laboratory leaching experiments have been used to explore some of the factors governing transport of sorbing contaminants as well as to characterize the sorption parameters (Dusek et al. 2010; Kulluru et al. 2010). However, a realistic implementation of boundary conditions (e.g., recharge rate (Leterme et al. 2012)) as well as inclusion of heterogeneity in soil permeability and thus assessment of their impacts on sorption kinetics can be challenging using these laboratory experiments. Dusek et al. (2015) suggested that the sorption parameters obtained using soil column experiments can be considerably different from the ones obtained in field leaching experiments; they attributed this discrepancy to the lack of realistic characterization of soil heterogeneity in laboratory soil column experiments.

So what is the way forward for understanding the links between subsurface permeability pattern, subsurface flow rate and sorption kinetics? Grid-based (discrete) numerical flow and advection–dispersion–reaction (ADR) contaminant transport models are potentially useful tools with which to explore controls on sorption characteristics. However, most numerical models used to characterize sorption behaviour have been restricted to 1-D homogenous unsaturated conditions (e.g., Dusek et al. 2011, 2015). Grid-based numerical subsurface flow and ADR contaminant transport models face high computational costs and numerical difficulties, such as numerical instability, numerical dispersion and artificial oscillation as the degree of subsurface heterogeneity increases (Salamon et al. 2006a, 2006b; Cherry & Clarke 2007; Starn et al. 2012; Ameli & Abedini 2016). This is particularly true when discrete numerical models are confronted with a defining feature of many shallow glacial till unconfined aquifers, where rapid exponential vertical changes in saturated hydraulic conductivity (ks) are the norm (e.g., Nyberg 1995; Seibert et al. 2011; Grip 2015).

Grid-free (continuous) analytical models can be robust alternatives to grid-based numerical schemes for simulation of subsurface flow and reactive transport (Destouni 1991). Such models, however, are often restricted to geometrically regular systems where the interactions between saturated and unsaturated zones, and between unconfined aquifers and surface water bodies are neglected or overly simplified (e.g., Wang et al. 2011). Recently, Ameli & Craig (2014) and Ameli et al. (2013) have relaxed the constraints of traditional analytical methods by enhancing these schemes with a simple numerical least square algorithm; the resulting semi-analytical groundwater–surface water interaction models were used to simulate 2-D and 3-D subsurface flow and flow pathlines in multi-layer unconfined aquifers. These coupled saturated–unsaturated solutions satisfy exactly the saturated and unsaturated (linearized) governing equations. These approaches more recently have extended to explicitly and systematically account for different rates of exponential decline in ks and porosity with soil depth (Ameli et al. 2016), a characteristic feature of many till-mantled unconfined environments in the boreal/hemiboreal landscape of Canada, Fennoscandia and many other places. This grid-free scheme provides continuous maps of head, velocity and dispersion tensors in the entire hillslope in response to arbitrary infiltration functions, without any interpolation processes. This allows for an explicit and efficient incorporation of the non-uniform Lagrangian random walk particle tracking (RWPT) transport method coupled with kinetic sorption algorithm. Thus, this model can be an effective tool for quickly and systematically posing hypotheses about how different degrees of vertical exponential decline in ks and porosity in interaction with other factors may be influencing sorption characteristic.

Here, first the semi-analytical model offered by Ameli et al. (2016) is calibrated (using non-linear parameter estimation software (PEST)) and tested against the observed hillslope hydrodynamics at an extensively studied till-mantled hillslope at the Västrabäcken subcatchment in Sweden (Laudon et al. 2013). The developed model, then, is coupled with a non-uniform RWPT algorithm, which incorporates a kinetic sorption model (Michalak & Kitanidis 2000), to theoretically assess the effect of subsurface permeability architecture and flow rate on sorption characteristics.

This paper specifically intends to:

  • (1)

    calibrate and assess the performance of the semi-analytical approach presented by Ameli et al. (2016) in simulating groundwater depth–discharge relationship along the S-transect hillslope in the Västrabäcken subcatchment;

  • (2)

    explore how subsurface heterogeneity in vertical pattern in ks and porosity, and longitudinal mechanical dispersion interact with infiltration rate and together influence the spatio-temporal pattern of sorption within the riparian zone as well as discharge loads of sorbing contaminant into stream.

Although the model developed here is based on the characteristic features of till-mantled systems, the conclusions on the impacts of permeability pattern and flow velocity on sorption characteristics can be extended to many other environments; indeed, the approach presented here can ‘systematically’ assess sorption behaviour for various permeability architectures by changing the rate of exponential decline in ks and porosity with soil depth.

METHOD

Natural hillslope description

The S-transect hillslope is located on the Västrabäcken subcatchment in northern Sweden (Figure 1). Daily stream discharge was taken from continuous measurements at a V-notch weir located half a kilometer downstream of the S-transect (Figure 1(a)). The 10th percentile (low flow), median, average and 90th percentile (high flow) daily discharge flow during the study period (from 01-01-2013 to 10-10-2013) were 0.2, 0.50, 0.80 and 1.80 mm/d, respectively. Daily groundwater depth measurements were collected with pressure transducers at the groundwater well located 16 m from the stream. The ks-depth relationship was measured in the subcatchment using the permeameter method (Bishop 1991). The best fit exponential function to the observed ks-depth data was m/d, where refers to the soil depth. In addition, the best exponential fit function to porosity-depth measurements of the mineral soil was .
Figure 1

Layout of the S-transect: (a) study location; (b) plan view of Västrabäcken subcatchment; (c) plan view of S-transect along with the location of the discharge measurement site and groundwater measurement well. (d) 2D cross section of the S-transect with a length of L = 140 m. The topographic surface was generated from a 5 m LiDAR DEM. The hillslope is characterized by an exponential decline in saturated hydraulic conductivity and porosity with soil depth (i.e., and ). [LT−1] and are the saturated hydraulic conductivity and porosity along the topographic surface (zt). α and ɳ are the parameters defining the exponential relationship between saturated hydraulic conductivity with soil depth (i.e., z–zt) and porosity with soil depth (i.e., z–zt), respectively. The circle shows the stream location. Figure developed originally by Thomas Grabs, Uppsala University.

Figure 1

Layout of the S-transect: (a) study location; (b) plan view of Västrabäcken subcatchment; (c) plan view of S-transect along with the location of the discharge measurement site and groundwater measurement well. (d) 2D cross section of the S-transect with a length of L = 140 m. The topographic surface was generated from a 5 m LiDAR DEM. The hillslope is characterized by an exponential decline in saturated hydraulic conductivity and porosity with soil depth (i.e., and ). [LT−1] and are the saturated hydraulic conductivity and porosity along the topographic surface (zt). α and ɳ are the parameters defining the exponential relationship between saturated hydraulic conductivity with soil depth (i.e., z–zt) and porosity with soil depth (i.e., z–zt), respectively. The circle shows the stream location. Figure developed originally by Thomas Grabs, Uppsala University.

Hypothetical hillslope description

We used virtual experiment framework (Weiler & McDonnell 2004) to explore how vertical exponential decline in ks and porosity, longitudinal mechanical dispersion and infiltration rate interact and together influence the sorption characteristic within the riparian zone as well as discharged loads of sorbing contaminant into stream. To do that four hypothetical rates (α) of exponential decline in ks with soil depth including α = 3 (extreme heterogeneous case with rapidly declining ks with depth), α = 2, α = 1, α = 0 (homogenous case) were considered. To focus only on the impact of ks vertical heterogeneity pattern, an identical average ks throughout the hillslope was assumed for four cases. Thus, while a ks0 [LT−1] (the saturated hydraulic conductivity along the topographic surface) of 100 m/d was assumed for the heterogeneous case with α = 3, of 85, 55 and 10 m/d were considered for the other cases with α = 2, α = 1 and α = 0, respectively. Furthermore, three hypothetical rates of exponential decline in porosity with soil depth including , , were considered with an identical (the porosity along the topographic surface) for three cases. Four values of uniformly distributed longitudinal dispersivity of 10, 1, 0.1 and 0 cm were also assumed. Note that for all examples solved in this paper the transverse dispersivity was assumed as . Four values of infiltration rates equal to the measured low flow (0.2), median flow (0.5), average flow (0.8) and high flow (1.8) at the S-transect were also used for these virtual experiments. The selected values of controlling factors used in the virtual experiments were from the range observed in till environments with the exception of zero values for , ɳ and which represent pure homogenous media.

Mathematical formulation

The governing equation of the steady-state saturated moisture movement in an unconfined aquifer with an exponentially declining relationship between saturated hydraulic conductivity (ks) and depth, can be represented in terms of the saturated discharge potential function as: 
formula
1
where is the parameter defining the exponential relationship between ks and depth and is the total hydraulic head. A series solution to the above governing equation can then be obtained as follows: 
formula
2
 
formula
where N is the total number of terms in the series solution, n denotes the coefficient index and , are the unknown series coefficients corresponding to the nth coefficient index.
In the unsaturated zone, the linearized form of steady-state Richard's equation in an unconfined aquifer with an exponential decline in ks with soil depth can be represented in terms of terms of Kirchhoff potential as follows: 
formula
3
The general series solution to Equation (3) can also be obtained as: 
formula
4
 
formula
where [L] represents pressure (suction) head, [L−1] is the sorptive number of Gardner model (Gardner 1958) and [L] is air entry pressure. Note that for all examples solved in this paper the Gardner sorptive number of = 1 m−1 and an air entry pressure of = 0.05 m were used. In addition, M denotes the total number of series terms, m represents the coefficient index and , are the unknown series coefficients associated with the mth coefficient index. Hereafter, and denote saturated and unsaturated properties/variables.

To complete the potential solutions (Equations (2) and (4)), the unknown coefficients (, and , ) are calculated by enforcing the boundary conditions. The a priori unknown free boundary between saturated and unsaturated domains (top of the capillary fringe interface, ), is also obtained through a robust iterative scheme. The water table elevation, , is then located as an interface with zero pressure head. Readers are referred to Ameli et al. (2016) for the details of the calculation of potential solutions, linearization of the unsaturated governing equation, implementation of boundary conditions and the iterative scheme used to determine the location of the free boundary condition.

Darcy fluxes in the saturated zone and both directions ( and ) and Darcy–Buckingham fluxes in the unsaturated zone and both directions ( and quz (x, z)) are then calculated as: 
formula
5a
 
formula
5b
A continuous map of pore water velocity in the entire domain then can be obtained as: 
formula
6a
 
formula
6b
where the saturated moisture content is equal to the porosity and is obtained as a function of soil depth . The unsaturated moisture content is also obtained based on both the suction pressure head and soil depth at each location .

Random walk particle tracking and kinetic sorption

A sorbing contaminant moves through the subsurface at a retarded flow velocity. The sorbing contaminant velocity, , can be estimated as pore water velocity divided by a dimensionless retardation factor (R) as (Michalak & Kitanidis 2000): 
formula
7a
 
formula
7b
 
formula
7c
where [ML−3] is soil bulk density and [L3M−1] is sorption distribution coefficient. The value represents the partitioning potential of a contaminant between sorbed and aqueous phases. This empirical coefficient depends upon geological conditions in the subsurface. Thus, varies with depth due to an exponential decline in with depth considered in this paper. It is known that there is a negative correlation between and (Tompson 1993; Qin et al. 2013). Indeed, sorption can be assumed to mechanistically relate to grain surface area. Therefore, a larger specific surface area leads to more enhanced sorptive capacity and a smaller hydraulic conductivity. In a manner similar to Tompson (1993) and Qin et al. (2013), is approximated as: 
formula
8
Parameter a2 is equal to −0.5 if the sorption property of the aquifer has perfect negative correlation with the grain size (Garabedian 1987; Tompson 1993). Similar to Qin et al. (2013) parameter a1 is also assumed as zero to ensure the negative relationship between and .
The calculated continuous maps of sorbing contaminant velocity (Equation (7)) are then used to track contaminant particles from the topographic surface to the surface water course using a non-uniform RWPT scheme. The non-uniform random walk step of a contaminant particle is given by: 
formula
9a
 
formula
9b
 
formula
and [L2T−1] are longitudinal and transverse hydrodynamic dispersion coefficients, and and [L] are longitudinal and transverse dispersivities of the porous medium, respectively. and are random numbers drawn from normal distributions with zero mean and unit variance for each particle and each time step . The asterisk denotes the correction of contaminant velocity via the implementation of non-uniformity in flow within the RWPT method. The corrected velocities at both directions are: 
formula
10a
 
formula
10b
In the above equation the tensor of dispersion is given by: 
formula
11a
 
formula
11b
 
formula
11c
where
Various methods to incorporate a kinetic sorption algorithm into the RWPT are available in the literature. Readers are referred to Valocchi & Quinodoz (1989) and Salamon et al. (2006a) for general reviews of these schemes. In this paper, for the transport of a linearly sorbing contaminant undergoing kinetic sorption, the method of moments (Michalak & Kitanidis 2000; Salamon et al. 2006b) is used. This method assumes that spatial moments of particles are identical to spatial moments of concentrations, provided that the number of particles is high enough. The normalized zeroth spatial moment of contaminant concentration in aqueous and sorbed phases represents the distribution of mass in these phases. Therefore, similar to Michalak & Kitanidis (2000) and Salamon et al. (2006b) the normalized zeroth spatial moment can be used as a phase transition probability function as follows: 
formula
12a
 
formula
12b
where is mass transfer coefficient [T−1], and is the probability that a particle starts in the aqueous phase and ends in the aqueous phase. Similarly, refers to the probability of a particle starting in the sorbed phase and ending in the sorbed phase. These transition probability functions can then be used to determine if a particle is in the aqueous phase and thus allowed to move via advection and dispersion or in the sorbed phase and thus not allowed to move after a time step . Time steps must be specified such that they are small enough to ensure that the chance of having more than one phase change is negligible in each time step. For each time step and each particle, a simple Bernoulli trial is used and randomly drawn numbers, from normal distributions with zero mean and unit variance, are compared to the corresponding probability. All particles are initially released in the mobile phase at the evenly spaced locations along the topographic surface. Additionally, all particles are assumed to have uniform mass. Mass arrival is measured at the surface water course located in the vicinity of the aquifer.

RESULTS

First, the semi-analytical series solution approach offered by Ameli et al. (2016) was calibrated using non-linear PEST software for ks0 and to match the simulated and observed relationship between groundwater depth and stream discharge (infiltration rate) at the annual median discharge (0.5 mm/d, the upward triangle in Figure 2) and at the groundwater well located 16 m from the stream (Figure 1). The model was then tested against the other discharge rates between low flow and high flow (the downward triangles in Figure 2). Three robust statistical Wilcoxon rank sum test, F-test and the Kolmogorov–Smirnov test were used to assess whether or not the median, variance and the cumulative distribution function of the observed and simulated groundwater depths (at specific stream discharges) are statistically similar. The decision to reject or accept the null hypotheses (similarity between observed and simulated) was based on comparing the p-values with the significant level, specified here at 5%.
Figure 2

Simulated (triangles) and observed (circles) groundwater depth at specific stream discharges. The groundwater depth was measured continuously at the groundwater well located 16 m from the stream within the S-transect hillslope and stream discharge was continuously measured at a V-notch weir located half a kilometer downstream of the S-transect (Figure 1). These data were collected from 01-01-2013 to 10-10-2013. The black dashed line (GW depth = 0.55–0.2*Log(q)) represents the best fit to the observed groundwater depth–discharge relationship with a R2 of 0.40. A statistical t-test of the slope suggested that the slope is statistically significant (p < <0.001) and groundwater depth is highly related to Log(q). A bootstrap analysis (based on 1,000 bootstrap samples) of the standard error of the coefficients of the best fit curve revealed the standard errors of 0.013 and 0.017 for intercept and slope, respectively. The upward triangle represents the simulated groundwater depth under the annual median stream discharge (0.5 mm/d) in the calibration phase. The downward triangles represent the simulated relationship between groundwater depth and discharge in the validation phase. The horizontal axis is in logarithmic scale.

Figure 2

Simulated (triangles) and observed (circles) groundwater depth at specific stream discharges. The groundwater depth was measured continuously at the groundwater well located 16 m from the stream within the S-transect hillslope and stream discharge was continuously measured at a V-notch weir located half a kilometer downstream of the S-transect (Figure 1). These data were collected from 01-01-2013 to 10-10-2013. The black dashed line (GW depth = 0.55–0.2*Log(q)) represents the best fit to the observed groundwater depth–discharge relationship with a R2 of 0.40. A statistical t-test of the slope suggested that the slope is statistically significant (p < <0.001) and groundwater depth is highly related to Log(q). A bootstrap analysis (based on 1,000 bootstrap samples) of the standard error of the coefficients of the best fit curve revealed the standard errors of 0.013 and 0.017 for intercept and slope, respectively. The upward triangle represents the simulated groundwater depth under the annual median stream discharge (0.5 mm/d) in the calibration phase. The downward triangles represent the simulated relationship between groundwater depth and discharge in the validation phase. The horizontal axis is in logarithmic scale.

The subsurface flow model was then coupled with a non-uniform RWPT scheme and kinetic sorption algorithm to theoretically investigate how subsurface permeability architecture interacts with flow rate and together influence sorption characteristics including sorbing contaminant pathways, transit time, spatial distribution as well as total number of sorption occurred within the hillslope; the latter is represented here as average number of sorption per contaminant pathline (Nsb). For each example, a total of 560 evenly spaced particles (i.e., 12.5 cm intervals) were released from the topographic surface and used within the non-uniform RWPT algorithm.

Model development

Figure 2 shows the observed groundwater depth and stream discharge as well as the best fit curve to their relationship. A t-test analysis of the slope of the best fit curve depicted a statistically significant relationship between groundwater depth in the well and stream discharge. A negligible bootstrapped standard error of the coefficients of the best fit curve further supports the established relationship between observed groundwater depth and stream discharge. The calibrated parameters were = 79 m/d and = 2.6 1/m which are consistent with the observed values ( = 84 m/d and = 2.46) at the S-transect. The simulated groundwater depth–discharge relationship in the validation phase was statistically similar to the observed one. The non-parametric Wilcoxon rank sum test, F-test and Kolmogorov–Smirnov test showed that median, variance and cumulative distribution of groundwater depths did not differ statistically between observed and simulated results at the significant level of 0.05; the corresponding p-values were 0.45, 0.98 and 0.33, respectively.

Effect of rate of exponential decline in with depth on subsurface sorption

Given an identical average throughout the hillslope, the more rapidly declines (larger ), the more superficial the contaminant circulation becomes (Figure 3). Thus the pathlines pass through more permeable porous media (both ks and porosity are higher in shallower portions). This together with the fact that the contaminant pathlines become (on average) shorter as increases lead to a decrease in the average number of sorption per contaminant pathline (Nsb). As changes from 0 to 3, Nsb decreases from 9.82 to 6.04.
Figure 3

Layout of sorbing contaminant pathline for different rates of exponential decline in with soil depth (α): (a) α = 3, (b) α = 2, (c) α = 1, (d) α = 0. Additionally, ɳ= 0.26, = 0.49 (measured values at S-transect), R= 1.8 × 10−4 (observed annual high flow), , , = 1.6 and = 0.5 were considered for this set of virtual experiment. Nsb refers to the average number of sorption per contaminant pathline. The pathlines of only 70 particles and their corresponding sorption out of a total of 560 particles used in the RWPT scheme are shown in this figure.

Figure 3

Layout of sorbing contaminant pathline for different rates of exponential decline in with soil depth (α): (a) α = 3, (b) α = 2, (c) α = 1, (d) α = 0. Additionally, ɳ= 0.26, = 0.49 (measured values at S-transect), R= 1.8 × 10−4 (observed annual high flow), , , = 1.6 and = 0.5 were considered for this set of virtual experiment. Nsb refers to the average number of sorption per contaminant pathline. The pathlines of only 70 particles and their corresponding sorption out of a total of 560 particles used in the RWPT scheme are shown in this figure.

The distribution of transit time of sorbing contaminants discharged into the river and spatial distribution of sorption within the hillslope are also influenced by the strength of decline in with soil depth (Figure 4). The proportion of contaminants discharged into the stream in the first year is significantly higher from the hillslope with a heterogeneous pattern compared to the homogenous hillslope; while in the second and third years this trend becomes opposite (Figure 4(a)). Furthermore, sorption occurrences in the deep portions of the aquifer relative to the shallow portions are weakened (Figure 4(b)). The average of depth of sorption (D0) also decreases from 2.38 to 1 m as α increases.
Figure 4

Impact of the rate of exponential decline in with soil depth on the distribution of contaminant transit time discharged into the stream and sorption spatial distribution within the hillslope. (a) Total mass fraction arrived in the stream at each year, (b) probability density function of depth of sorption. Nsb refers to the average number of sorption per pathline and D0 refers to the average of depth of sorption.

Figure 4

Impact of the rate of exponential decline in with soil depth on the distribution of contaminant transit time discharged into the stream and sorption spatial distribution within the hillslope. (a) Total mass fraction arrived in the stream at each year, (b) probability density function of depth of sorption. Nsb refers to the average number of sorption per pathline and D0 refers to the average of depth of sorption.

Effect of rate of exponential decline in porosity with depth on subsurface sorption

Different values equal to 0.50, 0.25 and 0 were used to assess the effect of rate of exponential decline in porosity with depth on subsurface transport of a sorbing contaminant. As increases, the vertical heterogeneity in porosity increases; represents a hillslope with homogenous porosity. Four different values of equal to 3, 2, 1 and 0 were also considered to investigate the impact of porosity pattern on sorption potential for various rates of exponential decline in with soil depth. For all different values (rate of exponential decline in ), an increase in the rate of exponential decline in porosity with depth increases the average number of sorption per pathline (Nsb) as well as the average depth of sorption (D0) (Figure 5). This can be attributed to a lower porosity and therefore larger retardation factor (Equation (7c)) in the deep portions of the aquifer compared to shallow parts as increases. However, the impacts of vertical heterogeneity in porosity on (1) the proportion of contaminant discharged into the stream in the first year, (2) the average number of sorption per pathline (Nsb) and (3) D0 are weakened as the strength of exponential decline in with soil depth increases. Indeed, porosity pattern becomes a more important control on subsurface sorption when the rate of vertical heterogeneity in ks decreases.
Figure 5

Impact of rate of exponential decline in porosity with soil depth on the distribution of contaminant transit time discharged into the stream and sorption spatial distribution within the hillslope for various rates of exponential decline in with soil depth . (a) Total mass fraction arrived in the stream at each year, (b) probability density function of depth of sorption. Nsb refers to the average number of sorption per pathline and D0 refers to the average of depth of sorption. In addition, for this set of virtual experiment, = 0.49, R= 1.8 × 10−4, = 10 cm, m = 0.1 cm, = 1.6 , = 0.5 were considered.

Figure 5

Impact of rate of exponential decline in porosity with soil depth on the distribution of contaminant transit time discharged into the stream and sorption spatial distribution within the hillslope for various rates of exponential decline in with soil depth . (a) Total mass fraction arrived in the stream at each year, (b) probability density function of depth of sorption. Nsb refers to the average number of sorption per pathline and D0 refers to the average of depth of sorption. In addition, for this set of virtual experiment, = 0.49, R= 1.8 × 10−4, = 10 cm, m = 0.1 cm, = 1.6 , = 0.5 were considered.

Effect of longitudinal dispersivity on subsurface sorption

Four different values of longitudinal dispersivity equal to 10, 1, 0.1 and 0 cm as well as different rates of exponential decline in with soil depth equal to 3, 2, 1 and 0 were considered to investigate the impact of pore scale heterogeneity on sorption potential for various degrees of macro-scale vertical heterogeneity in (Table 1). An increase in longitudinal dispersivity decreases the average number of sorption per pathline (Nsb) and increases the average depth of sorption for all different values. The impacts of longitudinal dispersivity on frequency and spatial distribution of sorption are enhanced as decreases. The relative change in the number of sorption that occurred within the hillslope (by increasing longitudinal dispersivity from 0 to 10 cm) increases from 21% to 34% as value decreases from 3 to 0.

Table 1

Impact of change in longitudinal dispersivity on the average number of sorption per pathline (Nsb) and the average depth of sorption (D0) for various rates of exponential decline in with soil depth

Nsb     
(1/m) 3 2 1 0 
(cm)     
0.0 8.14 8.81 12.80 14.81 
0.1 7.9 8.66 11.4 13 
7.4 8.0 11 12.33 
10 6.4 6.5 9.08 9.82 
Relative change −21% −26% −29% −34% 
D0     
(1/m) 3 2 1 0 
(cm)     
0.0 0.97 1.10 1.68 2.25 
0.1 0.97 1.13 1.72 2.30 
0.97 1.15 1.76 2.34 
10 1.16 1.80 2.45 
Relative change 3% 5% 7% 8% 
Nsb     
(1/m) 3 2 1 0 
(cm)     
0.0 8.14 8.81 12.80 14.81 
0.1 7.9 8.66 11.4 13 
7.4 8.0 11 12.33 
10 6.4 6.5 9.08 9.82 
Relative change −21% −26% −29% −34% 
D0     
(1/m) 3 2 1 0 
(cm)     
0.0 0.97 1.10 1.68 2.25 
0.1 0.97 1.13 1.72 2.30 
0.97 1.15 1.76 2.34 
10 1.16 1.80 2.45 
Relative change 3% 5% 7% 8% 

For this set of virtual experiments, ɳ=0.26, = 0.49, R= 1.8 × 10−4, = 1.6 , = 0.5 were considered. The last row refers to the relative change in the number of sorption that occurred within the hillslope by increasing the strength of longitudinal dispersivity

Effect of infiltration rate (R) on subsurface sorption

An increase in R decreases the average number of sorption per pathline (Nsb) as well as the average depth of sorption for all different values (Figure 6). Furthermore, the proportion of the arrival of sorbing contaminants within the first year increases as infiltration rate increases; while this trend is almost opposite in the following years. These impacts become most accentuated as the rate of exponential decline in with soil depth increases. As infiltration rate increases, the water table reaches the shallower portions of the hillslope and thus the shallow flow and contaminant circulation are enhanced (as shown in Figure 3, this behaviour becomes more pronounced with an increase in ). As the contaminant circulation becomes more superficial, the contaminant pathlines from land surface to the stream become shorter in length (thus the contact time between aqueous and solid phases decreases), and also contaminant pathlines traverse more permeable portions of the hillslope (since both and porosity become higher in shallower parts). These facts describe why an increase in infiltration rate weakens the sorption potential with a more pronounced impact for the hillslope with a higher rate of vertical decline in .
Figure 6

Impact of infiltration rate (R) on the distribution of contaminant transit time discharged into the stream and sorption spatial distribution within the hillslope for various rates of exponential decline in with soil depth . (a) Total mass fraction arrived in the stream at each year, (b) probability density function of depth of sorption. Nsb refers to the average number of sorption per pathline and D0 refers to the average of depth of sorption. Additionally, ɳ=0.26, = 0.49 (measured values at S-transect), = 10 cm, = 0.1 cm, = 1.6 and = 0.5 were considered for this set of virtual experiment.

Figure 6

Impact of infiltration rate (R) on the distribution of contaminant transit time discharged into the stream and sorption spatial distribution within the hillslope for various rates of exponential decline in with soil depth . (a) Total mass fraction arrived in the stream at each year, (b) probability density function of depth of sorption. Nsb refers to the average number of sorption per pathline and D0 refers to the average of depth of sorption. Additionally, ɳ=0.26, = 0.49 (measured values at S-transect), = 10 cm, = 0.1 cm, = 1.6 and = 0.5 were considered for this set of virtual experiment.

DISCUSSION

Subsurface sorption and transport of sorbing contaminants can be influenced by permeability architecture and subsurface flow velocity (e.g., Carlyle & Hill 2001; McCarty & Angier 2001; Fuchs et al. 2009). These factors have been difficult to assess as existing grid-based numerical subsurface flow and reactive transport models have had difficulties in incorporating heterogeneity in permeability, particularly in situations where permeability changes rapidly with soil depth (Salamon et al. 2006a, 2006b; Starn et al. 2012). Such heterogeneity in permeability is a characteristic feature of shallow till soils where porosity can vertically decline with an exponential rate up to 0.5 and saturated hydraulic conductivity may decline with soil depth with an exponential rate up to 4 (Grip 2015).

Recently, Ameli et al. (2016) proposed a new grid-free semi-analytical saturated–unsaturated subsurface flow and transport approach with an ability to efficiently account for rapid permeability changes with soil depth. The current study was designed to augment this semi-analytical model with a kinetic sorption algorithm to explore how subsurface vertical heterogeneity in ks and porosity interact with subsurface flow rate and together influence the frequency and spatial distribution of sorption in the riparian zone as well as the timing and rate of sorbing contaminant load discharged into the stream. Of course, the lack of validation of the reactive transport algorithm used here reduces what can be tested in this regard; but the efficiency of this integrated saturated–unsaturated flow and transport method in emulating various systematic patterns of vertical exponential decline in ks and porosity, as well as an efficient implementation of non-uniform RWPT reactive transport algorithm can help to systematically explore the impacts of permeability architecture and flow rate on subsurface sorption characteristics, and might promote new research questions where this approach could serve as a tool for rapid assessment.

Results suggest that as the strength of exponential decline in with soil depth increases (i.e., more rapidly declines with depth), deep contaminant circulation further away from the surface watercourse reduces and shallow contaminant circulation is enhanced. The total number of sorption (or average number of sorption per pathline) also decreases within the hillslope. More importantly, the distribution of depth of sorption favours shallower depths, with a higher sorption concentration close to the ground surface (Figures 3 and 4(b)). This, together with a lack of deep contaminant percolation, implies that the hillslope's capacity to retain sorbing contaminant reduces as vertical heterogeneity in increases, since the shallower soil depths may rapidly reach their sorption saturation level due to higher exposure to sorbing contaminants. An increase in the rate of vertical exponential decline in also increases the proportion of sorbing contaminant discharged into stream within the first year. The results also show that for a given porosity along the topographic surface, the more rapidly porosity exponentially declines with depth, the more sorption occurs within the subsurface. This effect is compounded as the rate of decline with soil depth decreases. These results are in good agreement with experimental studies which have shown that subsurface physical macro-heterogeneity can affect sorption potential (Carlyle & Hill 2001; McCarty & Angier 2001; Fuchs et al. 2009).

The results also suggest that increasing the infiltration rate (i.e., the velocity at which fluid enters the subsurface) decreases the sorption potential in the critical zone and increases the proportion of sorbing contaminant discharged into stream within the first year. This finding corroborates well with the results obtained based on small-scale laboratory flow-cell experiments performed by Fuchs et al. (2009), where they showed that an increase in subsurface flow velocity can weaken phosphorus sorption to soil minerals. The impact of infiltration rate on sorption potential in the critical zone and stream concentration of sorbing contaminants are enhanced as the rate of vertical decline increases.

Future research

The results and modelling approach presented in this paper can provide helpful understanding into the controls on subsurface transport of sorbing contaminants such as phosphorus and pesticides. Knowledge of vertical distribution of sorption and sorbing contaminant pathways within the critical zone as well as timing and rate of sorbing contaminant discharge into stream has potential implications for designing environmental remediation systems and engineered management practices in landscapes with non-point source phosphorus surplus. Coupling the integrated flow and reactive transport approach presented here with the field and laboratory measurements of sorbing contaminant concentration in the future can (1) help in a more certain estimation of poorly known sorption empirical parameters such as distribution and mass transfer coefficients and retardation factor (as has been done in Dusek et al. (2015)) and (2) aid in a robust assessment of the validity of different equilibrium and (non)linear kinetic models for sorption in various environments (e.g., Langmuir 1918; Crittenden et al. 1986; Pignatello & Xing 1995; Matott & Rabideau 2010; Singh et al. 2014).

CONCLUSION

The new grid-free semi-analytical integrated subsurface flow and reactive transport approach was used to investigate the impacts of vertical changes in permeability architecture as well as subsurface flow rate on sorption. The results suggested that changes in the rate of macro-scale vertical heterogeneity in subsurface permeability including exponential decline in saturated hydraulic conductivity with soil depth can significantly impact subsurface sorption frequency as well as sorbing contaminant pathline and arrival time into stream. Changes in the vertical pattern of porosity, changes in longitudinal dispersivity and infiltration rate can also affect the sorption characteristic. Their impacts, however, could be enhanced or weakened as macro-scale heterogeneity in saturated hydraulic conductivity varied. The findings are of practical relevance for developing management and conservation practices within the context of water quality and water pollution of integrated groundwater–surface water systems. Future research will be needed to couple the current approach with long-term sorbing contaminant concentration measurements to assess the validity of widely used kinetic sorption models and their corresponding equations.

ACKNOWLEDGEMENTS

Dr Jeffrey McDonnell, Dr Kevin Bishop, Dr Thomas Grabs, Dr Irena Creed, Mark Ranjram and Nino Amvrosiadi are thanked for their constructive support and motivations throughout the process. I thank SLU's Krycklan Catchment Study and the Swedish Infrastructure for Ecosystem Science Svartberget Field Station for collecting the hydrometric data.

REFERENCES

REFERENCES
Bishop
K. H.
1991
Episodic Increases in Stream Acidity, Catchment Flow Pathways and Hydrograph Separation
.
University of Cambridge
,
Cambridge
, pp.
246
.
Cherry
G. S.
Clarke
J. S.
2007
Simulation and particle-tracking analysis of selected ground-water pumping scenarios at Vogtle Electric Generation Plant, Burke County, Georgia. 2331-1258, Geological Survey (US)
.
Crittenden
J. C.
Hutzler
N. J.
Geyer
D. G.
Oravitz
J. L.
Friedman
G.
1986
Transport of organic compounds with saturated groundwater flow: model development and parameter sensitivity
.
Water Resources Research
22
(
3
),
271
284
.
Dusek
J.
Sanda
M.
Loo
B.
Ray
C.
2010
Field leaching of pesticides at five test sites in Hawaii: study description and results
.
Pest Management Science
66
(
6
),
596
611
.
Dusek
J.
Dohnal
M.
Vogel
T.
Ray
C.
2011
Field leaching of pesticides at five test sites in Hawaii: modeling flow and transport
.
Pest Management Science
67
(
12
),
1571
1582
.
Elfeki
A. M.
Uffink
G. J.
Barends
F. B.
1997
Groundwater Contaminant Transport: Impact of Heterogenous Characterization: A New View on Dispersion
.
CRC Press
,
Boca Raton, FL
,
USA
.
Fuchs
J. W.
Fox
G. A.
Storm
D. E.
Penn
C. J.
Brown
G. O.
2009
Subsurface transport of phosphorus in riparian floodplains: influence of preferential flow paths
.
Journal of Environmental Quality
38
(
2
),
473
484
.
Garabedian
S. P.
1987
Large-Scale Dispersive Transport in Aquifers: Field Experiments and Reactive Transport Theory
.
PhD Thesis
,
Massachusetts Institute of Technology
.
Kulluru
P.
Das
B.
Panda
R.
2010
Evaluation of sorption and leaching potential of malathion and atrazine in agricultural soils of India
.
International Journal of Environmental Research
4
(
1
),
75
.
Langmuir
I.
1918
The adsorption of gases on plane surfaces of glass, mica and platinum
.
Journal of the American Chemical Society
40
(
9
),
1361
1403
.
Laudon
H.
Taberman
I.
Agren
A.
Futter
M.
Ottosson-Lofenius
M.
Bishop
K.
2013
The Krycklan Catchment Study – a flagship infrastructure for hydrology, biogeochemistry, and climate research in the boreal landscape
.
Water Resources Research
49
(
10
),
7154
7158
.
Leterme
B.
Mallants
D.
Jacques
D.
2012
Sensitivity of groundwater recharge using climatic analogues and HYDRUS-1D
.
Hydrology and Earth System Sciences
16
(
8
),
2485
2497
.
McCarty
G.
Angier
J.
2001
Impact of preferential flow pathways on ability of riparian wetlands to mitigate agricultural pollution. Preferential flow: water movement and chemical transport in the environment
. In:
Proceedings of the 2nd International Symposium
,
January 3–5, 2001
.
American Society of Agricultural Engineers
,
Ala Moana Hotel, Honolulu, Hawaii, USA
, pp.
53
56
.
Michalak
A. M.
Kitanidis
P. K.
2000
Macroscopic behavior and random-walk particle tracking of kinetically sorbing solutes
.
Water Resources Research
36
(
8
),
2133
2146
.
Pignatello
J. J.
Xing
B.
1995
Mechanisms of slow sorption of organic chemicals to natural particles
.
Environmental Science & Technology
30
(
1
),
1
11
.
Salamon
P.
Fernàndez-Garcia
D.
Gómez-Hernández
J. J.
2006a
A review and numerical assessment of the random walk particle tracking method
.
Journal of Contaminant Hydrology
87
(
3
),
277
305
.
Salamon
P.
Fernàndez-Garcia
D.
Gómez-Hernández
J.
2006b
Modeling mass transfer processes using random walk particle tracking
.
Water Resources Research
42
(
11
).
Seibert
J.
Bishop
K.
Nyberg
L.
Rodhe
A.
2011
Water storage in a till catchment. I: distributed modelling and relationship to runoff
.
Hydrological Processes
25
(
25
),
3937
3949
.
Singh
A.
Allen-King
R. M.
Rabideau
A. J.
2014
Groundwater transport modeling with nonlinear sorption and intraparticle diffusion
.
Advances in Water Resources
70
,
12
23
.
Starn
J. J.
Bagtzoglou
A. C.
Robbins
G. A.
2012
Methods for simulating solute breakthrough curves in pumping groundwater wells
.
Computers & Geosciences
48
,
244
255
.
Tilman
D.
Fargione
J.
Wolff
B.
D'Antonio
C.
Dobson
A.
Howarth
R.
Schindler
D.
Schlesinger
W. H.
Simberloff
D.
Swackhamer
D.
2001
Forecasting agriculturally driven global environmental change
.
Science
292
(
5515
),
281
284
.
Valocchi
A. J.
Quinodoz
H. A.
1989
Application of the random walk method to simulate the transport of kinetically adsorbing solutes
.
Groundwater Contamination
185
,
35
42
.
Vitousek
P. M.
Naylor
R.
Crews
T.
David
M. B.
Drinkwater
L. E.
Holland
E.
Johnes
P. J.
Katzenberger
J.
Martinelli
L. A.
Matson
P. A.
Nziguheba
G.
Ojima
D.
Palm
C. A.
Robertson
G. P.
Sanchez
P. A.
Townsend
A. R.
Zhang
F. S.
2009
Nutrient imbalances in agricultural development
.
Science
324
(
5934
),
1519
1520
.