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.
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:
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;
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.
Natural hillslope description
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.
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
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.
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
Effect of rate of exponential decline in porosity with depth on subsurface sorption
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.
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
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.
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).
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.
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.