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
m/d, where refers to the soil depth. In addition, the best exponential fit function to porosity-depth measurements of the mineral soil was .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.
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
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.
Random walk particle tracking and kinetic sorption
RESULTS
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%.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.
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
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.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.
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.
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.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.
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
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.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.
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.
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 | ||||
| 3 | 2 | 1 | 0 | |
| 0.0 | 8.14 | 8.81 | 12.80 | 14.81 |
| 0.1 | 7.9 | 8.66 | 11.4 | 13 |
| 1 | 7.4 | 8.0 | 11 | 12.33 |
| 10 | 6.4 | 6.5 | 9.08 | 9.82 |
| Relative change | −21% | −26% | −29% | −34% |
| D0 | ||||
| 3 | 2 | 1 | 0 | |
| 0.0 | 0.97 | 1.10 | 1.68 | 2.25 |
| 0.1 | 0.97 | 1.13 | 1.72 | 2.30 |
| 1 | 0.97 | 1.15 | 1.76 | 2.34 |
| 10 | 1 | 1.16 | 1.80 | 2.45 |
| Relative change | 3% | 5% | 7% | 8% |
| Nsb | ||||
| 3 | 2 | 1 | 0 | |
| 0.0 | 8.14 | 8.81 | 12.80 | 14.81 |
| 0.1 | 7.9 | 8.66 | 11.4 | 13 |
| 1 | 7.4 | 8.0 | 11 | 12.33 |
| 10 | 6.4 | 6.5 | 9.08 | 9.82 |
| Relative change | −21% | −26% | −29% | −34% |
| D0 | ||||
| 3 | 2 | 1 | 0 | |
| 0.0 | 0.97 | 1.10 | 1.68 | 2.25 |
| 0.1 | 0.97 | 1.13 | 1.72 | 2.30 |
| 1 | 0.97 | 1.15 | 1.76 | 2.34 |
| 10 | 1 | 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
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 .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.
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.
