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 *K _{s}* 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

*K*. 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.

_{s}## 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 (*k _{s}*) 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 *k _{s}* 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

*k*and porosity in interaction with other factors may be influencing sorption characteristic.

_{s}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

*k*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._{s}

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 *k _{s}* and porosity with soil depth.

## METHOD

### Natural hillslope description

*k*-depth relationship was measured in the subcatchment using the permeameter method (Bishop 1991). The best fit exponential function to the observed

_{s}*k*-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 .

_{s}### Hypothetical hillslope description

We used virtual experiment framework (Weiler & McDonnell 2004) to explore how vertical exponential decline in *k _{s}* 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

*k*with soil depth including

_{s}*α*= 3 (extreme heterogeneous case with rapidly declining

*k*with depth),

_{s}*α*= 2,

*α*= 1,

*α*= 0 (homogenous case) were considered. To focus only on the impact of

*k*vertical heterogeneity pattern, an identical average

_{s}*k*throughout the hillslope was assumed for four cases. Thus, while a

_{s}*k*

_{s}_{0}[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

*k*) and depth, can be represented in terms of the saturated discharge potential function as: where is the parameter defining the exponential relationship between

_{s}*k*and depth and is the total hydraulic head. A series solution to the above governing equation can then be obtained as follows: where

_{s}*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

*n*coefficient index.

^{th}*k*with soil depth can be represented in terms of terms of Kirchhoff potential as follows: The general series solution to Equation (3) can also be obtained as: where [L] represents pressure (suction) head, [L

_{s}^{−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

*m*coefficient index. Hereafter, and denote saturated and unsaturated properties/variables.

^{th}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.

*q*(

_{uz}*x*,

*z*)) are then calculated as: A continuous map of pore water velocity in the entire domain then can be obtained as: 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

*R*) as (Michalak & Kitanidis 2000): where [ML

^{−3}] is soil bulk density and [L

^{3}M

^{−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: Parameter

*a*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

_{2}*et al.*(2013) parameter

*a*is also assumed as zero to ensure the negative relationship between and .

_{1}^{2}T

^{−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: In the above equation the tensor of dispersion is given by: where

*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: 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

*et al.*(2016) was calibrated using non-linear PEST software for

*k*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

_{s0}*p*-values with the significant level, specified here at 5%.

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 (*N _{sb}*). 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

*k*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 (

_{s}*N*). As changes from 0 to 3,

_{sb}*N*decreases from 9.82 to 6.04.

_{sb}*D*) also decreases from 2.38 to 1 m as

_{0}*α*increases.

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

*N*) as well as the average depth of sorption (

_{sb}*D*) (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 (

_{0}*N*) and (3)

_{sb}*D*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

_{0}*k*decreases.

_{s}### 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 (*N _{sb}*) 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.

N_{sb} | ||||

(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 |

1 | 7.4 | 8.0 | 11 | 12.33 |

10 | 6.4 | 6.5 | 9.08 | 9.82 |

Relative change | −21% | −26% | −29% | −34% |

D_{0} | ||||

(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 |

1 | 0.97 | 1.15 | 1.76 | 2.34 |

10 | 1 | 1.16 | 1.80 | 2.45 |

Relative change | 3% | 5% | 7% | 8% |

N_{sb} | ||||

(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 |

1 | 7.4 | 8.0 | 11 | 12.33 |

10 | 6.4 | 6.5 | 9.08 | 9.82 |

Relative change | −21% | −26% | −29% | −34% |

D_{0} | ||||

(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 |

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

*R*decreases the average number of sorption per pathline (

*N*) 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 .

_{sb}## 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 *k _{s}* 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

*k*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.

_{s}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.