A model-based assessment of the potential use of compound-specific stable isotope analysis in river monitoring of diffuse pesticide pollution

Compound-specific stable isotope analysis (CSIA) has, in combination with model-assisted interpretation, proven to be a valuable approach to quantify the extent of organic contaminant degradation in groundwater systems. CSIA data may also provide insights into the origin and transformation of diffuse pollutants, such as pesticides and nitrate, at the catchment scale. While CSIA methods for pesticides have increasingly become available, they have not yet been deployed to interpret isotope data of pesticides in surface water. We applied a coupled subsurface-surface reactive transport model (HydroGeoSphere) at the hillslope scale to investigate the usefulness of CSIA in the assessment of pesticide degradation. We simulated the transport and transformation of a pesticide in a hypothetical but realistic twodimensional hillslope transect. The steady-state model results illustrate a strong increase of isotope ratios at the hillslope outlet, which resulted from degradation and long travel times through the hillslope during average hydrological conditions. In contrast, following an extreme rainfall event that induced overland flow, the simulated isotope ratios dropped to the values of soil water in the pesticide application area. These results suggest that CSIA can help to identify rainfallrunoff events that entail significant pesticide transport to the stream via surface runoff. Simulations with daily rainfall and evapotranspiration data and one pesticide application per year resulted in small seasonal variations of concentrations and isotope ratios at the hillslope outlet, which fell within the uncertainty range of current CSIA methods. This implies a good reliability of in-stream isotope data in the absence of transport via surface runoff or other fast transport routes, since the time of measurement appears to be of minor importance for the assessment of pesticide degradation. The analysis of simulated isotope ratios also allowed quantification of the contribution of two different reaction pathways (aerobic and anaerobic) to overall degradation, which gave further insight into the transport routes in the modelled system. The simulations supported the use of the commonly applied Rayleigh equation for the interpretation of CSIA data, since this led to an underestimation of the real extent of degradation of less than 12 % at the hillslope outlet. Overall, this study emphasizes the applicability and usefulness of CSIA in the assessment of diffuse river pollution, and represents a first step towards a theoretical framework for the interpretation of CSIA data in agricultural catchments.


Introduction
Modern agriculture makes use of a variety of pesticides to increase crop yield and reduce pests and the growth of weeds.As a result, pesticides have become ubiquitous organic contaminants in agricultural catchments.Diffuse pollution by pesticides can pose a risk for the terrestrial and aquatic environment, and human health.Pesticide residuals and their metabolites have been found in groundwater and surface water and affect drinking water quality (Donald et al., 2007; Published by Copernicus Publications on behalf of the European Geosciences Union.Kolpin et al., 1998;Kjaer et al., 2005).It is therefore vital to assess the fate of diffuse pollutants, and to identify major pesticide sources in a catchment.

S. R. Lutz et al.: Potential use of CSIA in river monitoring
After application, pesticides are subject to various transfer, transformation, and transport processes (Gavrilescu, 2005;Flury, 1996).Important transfer processes are volatilization (the transfer of compounds from the solid or liquid phase to the gas phase), and sorption (the transfer from the liquid to the solid phase of the soil matrix).Pesticide molecules that have not been volatilized can undergo transformation processes on the ground surface and in the soil that lead to destruction of the compound.The transformation of pesticides is either ascribed to microbial activity or abiotic processes.Microbial degradation in the soil matrix can occur under aerobic or anaerobic conditions.The transport of pesticides in the aqueous phase towards surface water bodies can occur via surface runoff or subsurface flow.If the pesticide is not directly removed by surface runoff, it leaches into the subsurface.Drainage systems, preferential flow (e.g. in earthworm burrows and cracks) and subsurface storm flow can thereupon cause fast transport to surface water (Gavrilescu, 2005;Müller et al., 2003;Leu et al., 2004b).In contrast, pesticide leaching to groundwater represents a slow subsurface transport mechanism (Holvoet et al., 2007;Flury, 1996).
Previous research has shown that rainfall events increase the risk of surface water contamination by pesticides (Taghavi et al., 2011;Müller et al., 2003).The highest pesticide loads happen in response to rainfall events shortly after pesticide application, and are attributed to surface runoff and preferential flow (Leu et al., 2004a).A secondary contribution to overall pesticide export is ascribed to the baseflow component (Louchart et al., 2001;Squillace and Thurmanz, 1992).
Monitoring of pesticide concentrations in rivers enables us to assess the extent of diffuse pollution at the catchment scale.However, concentration data do not provide clear evidence of degradation processes, since a reduction in concentration might, for example, indicate degradation of the contaminant, changes in the application pattern and amount, or dilution (Battaglin and Goolsby, 1999;Schreglmann et al., 2013).In contrast to degradation, non-destructive processes such as dilution do not reduce the contaminant loads to subsurface and surface water.In this context, compound-specific isotope analysis (CSIA) has emerged as a valuable tool for the analysis of the transformation of organic contaminants.
CSIA is the measurement of the isotopic composition, i.e. the ratio between the abundance of a heavy and a light stable isotope of an element in a compound.This isotope ratio tends to increase during degradation since chemical bonds that contain a heavy isotope are less amenable to degradation processes than those that solely consist of light isotopes (Meckenstock et al., 2004;Elsner, 2010).This phenomenon, called isotope fractionation, does not occur during dilution (van Breukelen, 2007a), and it is only relevant under specific conditions for other physical non-destructive processes such as diffusion (van Breukelen and Rolle, 2012) and sorption (van Breukelen and Prommer, 2008).In contrast to concentration monitoring, CSIA thus serves as a direct indicator of the transformation of a contaminant and, by measuring multiple elements, even allows for the determination of the degradation mechanism that leads to the transformation (Elsner, 2010).
Reactive transport models that incorporate isotope fractionation effects (isotope fractionation reactive transport models, IF-RTMs) have become a popular method to model and interpret CSIA data from point-source pollutants in groundwater systems (Atteia et al., 2008;D'Affonseca et al., 2011;Pooley et al., 2009;Prommer et al., 2009;van Breukelen et al., 2005;Wanner et al., 2012b).In particular, assumptions concerning reaction kinetics can be validated by comparing model results to measured concentration and CSIA data (D'Affonseca et al., 2011;Pooley et al., 2009;Prommer et al., 2009;Atteia et al., 2008).IF-RTMs have also proven useful to study the transformation of the agrochemical pollutant nitrate (Chen and MacQuarrie, 2004;Green et al., 2010).As opposed to the simulation of concentrations of a degrading contaminant, which solely yields a certain concentration reduction between the emission source and a downgradient point, these studies have demonstrated that IF-RTMs allow for the quantification of in situ degradation and the distinction between destructive and non-destructive processes.Moreover, IF-RTMs permit the assessment of the performance of the Rayleigh equation, which is the mathematical basis of CSIA.The Rayleigh equation tends to underestimate the extent of biodegradation in real flow systems because of the attenuation of isotopic enrichment due to dispersion or mixing processes (Abe and Hunkeler, 2006;van Breukelen and Prommer, 2008).
The aforementioned IF-RTM studies have solely focussed on groundwater systems.In contrast, the current study presents the first incorporation of isotope fractionation effects into a reactive transport model of a coupled surfacesubsurface system, which involves more complex transport routes.The aim of the modelling study was to determine whether CSIA measured in surface water can help to identify transport routes and quantify the extent of degradation of diffuse agricultural pollutants in subsurface-surface systems.We modelled a hypothetical but realistic situation of pesticide application, transport, and degradation, including isotope fractionation at the hillslope scale.To this end, we applied a distributed physically based model, which allowed for the simulation of hydrological processes in different flow domains in a detailed and spatially explicit way (Kampf and Burges, 2007) and the description of pesticide fluxes and concentrations (Thorsen et al., 1996).This study thus ties in with the virtual experiment approach in hillslope hydrology of Weiler and McDonnell (2004), Hopp et al. (2009), and Mirus et al. (2011), and additionally considers reactive transport and isotope fractionation.We opted for a hillslope transect because hillslopes are a fundamental landscape element and form the basic hydrological unit of catchments.Hence, understanding processes at the hillslope scale is an important and relevant first step towards interpreting CSIA at the catchment scale.By means of scenario modelling, we examined the evolution of isotope ratios under average hydrological conditions, in response to an extreme rainfall event, and under transient, daily varying hydrological conditions.In order to advance the transition from these virtual experiments to real applications, we give, based on the simulation results for this hillslope model, general guidelines for the monitoring and interpretation of isotope data in the context of diffuse pollution.

Model code description: HydroGeoSphere
We simulated isotope fractionation during the transport and transformation of a hypothetical pesticide with the software HydroGeoSphere (HGS) (Brunner and Simmons, 2012;Therrien et al., 2010).HGS is a fully coupled, subsurfacesurface flow and solute transport model.Surface flow is described by the diffusion-wave approximation of the Saint Venant equation, while the variably saturated form of the Richards equation is used for the subsurface.The Newton-Raphson technique is implemented to solve the non-linear equations of variably saturated flow.Solute transport is simulated by solving the advection-dispersion equation; the degradation of the solute is modelled with first-order kinetics.Furthermore, it is possible to include interception, transpiration and evaporation processes.Applications of HGS include large watershed modelling (Goderniaux et al., 2009;Li et al., 2008), simulations of aquifer-river interactions in hypothetical model domains (Doble et al., 2012;McCallum et al., 2010), and the analysis of contaminant transport in groundwater (Rivett et al., 2006;Sudicky et al., 2010).

Hillslope geometry and model grid
The model domain consists of a two-dimensional hillslope that is 200 m long and stretches 15 m in the vertical direction (Fig. 1).Conceptually, the hillslope represents part of a headwater catchment.The hillslope outlet corresponds to a river monitoring point.Inflow of water and pesticides from upstream parts was not considered in this model.The hillslope is convex with an average gradient of 5 %, which is comparable to the average slope in previously studied agricultural catchments (Doppler et al., 2012;Leu et al., 2004a).A river bank was incorporated as a vertical drop of 2 m over the last 5 m in the x direction, which increased the average slope to 5.8 %.The model domain does not represent, nor was calibrated for a specific field site, but is supposed to resemble a realistic agricultural hillslope.
Figure 1a shows the model mesh and the three subsurface zones (see Sect. 2.3) in the hillslope.The subsurface comprises a total of 43 layers.The vertical grid spacing is a few centimetres in the topsoil and increases with depth, which results in a maximum cell height of more than one meter at the bottom of the model domain.The mesh nodes are 0.5 m apart in the horizontal direction, except for the zone near the river, which has a horizontal discretization of 0.25 m (Fig. 1b).

Hydraulic properties and flow simulation
We wanted to model a layered system with variable hydraulic properties, as many sites are characterized by relatively shallow soils on fractured bedrock.The hillslope model therefore consists of three zones with different properties: the topsoil extends 0.3 m below the ground surface, the subsoil is located between 0.3 and 2 m below the ground surface, and the remaining part of the subsurface represents the bedrock.A similar layering was used for the numerical simulation of runoff generation mechanisms in different catchments by Mirus et al. (2011).The saturated hydraulic conductivity for the three subsurface zones was set as follows: 1.0 m d −1 in the topsoil, 0.5 m d −1 in the subsoil, and 0.1 m d −1 in the bedrock (Table 1).These values represent the frequently observed decline in hydraulic conductivity with depth and are comparable to those used by Christiansen et al. (2004) to simulate pesticide transport in a small agricultural catchment in Denmark.The hydraulic parameters for the soil were obtained from an experimental site at Canadian Forces Base Borden (north of Toronto, Ontario, Canada) because HGS has been validated against irrigation experiments at this site (Abdul, 1985;Therrien et al., 2010).These simulations therefore provided verified HGS parameters for the porosity, the residual saturation, and the parameters α and β of the Van Genuchten model (which determine the saturation-pressure relation).
Previous HGS simulations of the Borden site used uniform soil parameters and only considered the upper 4 m of the profile (Therrien et al., 2010), but we preferred to explicitly represent the bedrock.The bedrock parameters were therefore taken from simulations of the forested Coos Bay site in Oregon, USA (Mirus et al., 2011).They represent weathered bedrock and fractured sandstone.The specific storage coefficient was kept at its default value (1 × 10 −4 m −1 ) for all three zones.HGS explicitly represents the overland flow domain as a layer on top of the subsurface domain.The nodes of the overland flow domain coincide with the top elements of the subsurface domain.The overland flow parameters rill storage height and coupling length were both set to 0.01 m based on previous studies using HGS (Goderniaux et al., 2009;Therrien et al., 2010).The rill storage height defines the depth of depressions on the ground surface, which inhibit the generation of surface runoff.The coupling length determines the degree of continuity in pressure heads between the surface and the subsurface domain and thus defines the exchange flux between these domains (Verbist et al., 2012).A smaller coupling length results in increased coupling of the overland flow domain to the subsurface domain.
The lateral and bottom boundaries of the hillslope were set to zero-flux boundaries such that the water could leave the hillslope only via the overland flow domain.This was realized by applying a critical depth boundary at the hillslope outlet, which yields a time varying flux that is given by Manning's equation (Therrien et al., 2010).The hillslope outlet was represented by two nodes in the surface layer, since the hillslope extends one cell into the y direction.Rainfall was applied at a spatially uniform rate on the entire model domain.
Rainfall and potential evapotranspiration data were supplied as model input.The actual evapotranspiration in HGS depends on the moisture content and several vegetationrelated parameters.The difference between precipitation and actual evapotranspiration is then added to or abstracted from the model domain, depending on which component is dominant.The parameters for evapotranspiration were kept at their default values, except for the leaf area index, which was set to the value for temperate and tropical crops given in the HGS manual (leaf area index of 4.2), and the interception by canopy, which was assumed to be negligible.

Reactive solute transport
We used HGS to simulate the reactive transport of a soluble, non-volatile, and non-sorbing hypothetical pesticide.These properties are characteristic of the widely used compounds MCPA, bentazone, metam-sodium and clopyralid (University of Hertfordshire, 2013), as well as nitrate.We therefore neglected volatilization and sorption processes.The pesticide was applied to the cells at the ground surface between x coordinates 50 m and 130 m (Fig. 1a).Depending on the simulated scenario, the solute transport boundary condition was either set to a specified concentration or a specified mass flux (see Sect. 2.6).The transport parameters were chosen as follows: an aqueous diffusion coefficient of 7.8 × 10 −5 m 2 d −1 , which is the value that was assumed for diffusion of the herbicides metolachlor and alachlor in groundwater by Lee and Benson (2004); and the default values for dispersivities in HGS of 1.0 m in the longitudinal direction and 0.1 m in the vertical transverse direction.
The hypothetical pesticide was assumed to degrade via two different pathways: an aerobic reaction in the topsoil and an anaerobic reaction in the subsoil and bedrock.A pesticide half-life of 51 days, corresponding to an overall degradation rate constant of 5 per year, was chosen for the aerobic reaction (Table 2), which is within the range of reported half-lives of widely used pesticides (e.g.linuron and clopyralid; University of Hertfordshire, 2013).The half-life a Overall degradation rate constant (Eq.3); b enrichment factor for carbon isotope fractionation; c enrichment factor for hydrogen isotope fractionation.
for the anaerobic reaction was assumed to be much longer (1265 days; rate constant of 0.2 yr −1 ).
The model yielded concentrations for the entire model domain.Additionally, for each time step and solute, the concentration at the hillslope outlet was obtained by tracking the solute and water fluxes through the boundary notes.These concentrations were verified against the model-based concentrations in the overland flow domain at the boundary notes.
The simulations were run in finite difference mode with upstream weighting.The adaptive time stepping scheme was used; the maximum time step was set to one day for the steady-state simulations and 0.01 days for the transient simulation (see Sect. 2.6).A maximum concentration change of 0.1 relative to the source concentration was allowed in each time step for the steady-state simulation; this value had to be increased to 500 for the transient simulation because of the temporary high mass load.

Simulation of isotope fractionation effects
Isotope fractionation effects were included in the model by simulating the concentrations of the light and heavy isotopes of the pesticide.Since the assumption of two degradation pathways requires the analysis of two isotopic elements (van Breukelen, 2007b), this resulted in the simulation of four solutes.The two-dimensional isotope analysis was implemented by simulating carbon and hydrogen isotopes because carbon represents the element with the most isotope data available for pesticides (e.g.Annable et al., 2007;Badea et al., 2009;Meyer et al., 2009;Penning et al., 2010), and the cleavage of chemical bonds in pesticide molecules with hydrogen atoms induces a significantly stronger isotope fractionation effect than, for example, with nitrogen atoms (Hartenbach et al., 2008;Penning et al., 2010).
Following Hunkeler et al. (2009), van Breukelen et al. (2005), and van Breukelen and Prommer (2008), the differential equations for the degradation of the light carbon and hydrogen isotopes were simulated as a first-order reaction with the degradation rate constant L k: with L C being the concentration of the light carbon or hydrogen isotopes, respectively.The reaction kinetics of the heavy isotopes, involving a different degradation rate constant H k, were specified as where H C is the concentration of the heavy carbon and hydrogen isotopes, respectively.The kinetic isotopic fractionation factor, α, represents the ratio between the degradation rate constants of the heavy and light isotopes (i.e.H k / L k).It thus determines the strength of the isotope fractionation effect for a specific reaction.Since it typically has a value close to one, it is reported in per mil (‰) as the kinetic isotopic enrichment factor ε (‰) = (α − 1) • 1000.The degradation rate constant of the light isotopes ( L k) was set to the overall degradation rate constant (k eff , Table 2).Given α = H k / L k and L k = k eff , the degradation of the heavy carbon and hydrogen isotopes was thus simulated as . The enrichment factors (ε) for aerobic and anaerobic degradation for the carbon isotopes were chosen to be representative for fractionation effects during biotic pesticide degradation (Meyer et al., 2009;Penning et al., 2010).Data about fractionation effects for hydrogen isotopes are much scarcer; the enrichment factors for hydrogen were therefore set to values following the general trend of stronger enrichment in hydrogen compared to carbon isotopes (Hunkeler and Elsner, 2009).
In a system with two transformation pathways, the ratio between the enrichment factors of two isotopic elements, for example of hydrogen and carbon (ε H /ε C ), can be indicative of a specific degradation mechanism (Meyer et al., 2009;Fischer et al., 2008).The enrichment factors for the two reaction pathways were therefore chosen such that they yield distinct ε H / ε C -ratios (Table 2).This is in agreement with considerably diverging ε H / ε C -ratios that were, for example, observed for aerobic and anaerobic biodegradation of MTBE (Zwank et al., 2005).
Combining and integrating Eqs. ( 1) and ( 2) leads to the Rayleigh equation, which allows for the quantification of in situ degradation on the basis of isotope ratios.Its simplified form can be expressed as where R represents the ratio between the abundance of a heavy and a light isotope of an element in a compound for a sample (R sample ) and at the emission source (R source ), respectively, and f denotes the non-degraded fraction of the compound in the sample with respect to the emission source.
To facilitate inter-sample comparison, the isotope ratio of a sample (R sample ) is expressed in the δ-notation, which is the relative difference of R sample from a standard ratio R standard S. R. Lutz et al.: Potential use of CSIA in river monitoring (Schmidt and Jochmann, 2012): where the δ-value or isotopic signature is commonly reported in per mil (‰).The standard ratios that were used in this study are the international standards for carbon and hydrogen isotope ratios, i.e.Vienna Pee Dee Belemnite (VPDB; R standard = 0.011237) and Vienna Standard Mean Ocean Water (VSMOW; R standard = 1.5575 × 10 −4 ), respectively.The δ 13 C-value at the source was fixed at −30 ‰, which corresponds to typical values for carbon isotopes in pesticides (Annable et al., 2007;Kawashima and Katayama, 2010;Badea et al., 2009).Since no typical source values were available for hydrogen isotopes in pesticides, a value of −100 ‰ was chosen, which is consistent with δ 2 H-values reported for other organic contaminants (Mancini et al., 2008;Wang et al., 2004).These initial δ values were used to determine the ratio of the concentrations of light and heavy isotopes at the pesticide source (Eq.5).The simulated concentrations of the light and heavy isotopes were used to calculate the δ values for the entire hillslope domain and the hillslope outlet.

Simulated scenarios
In a preliminary model run, a recharge rate of 250 mm yr −1 was applied to the whole surface domain in order to achieve a steady-state flow field.This value was considered representative for average hydrological conditions for the Netherlands and northern Germany (Otto, 2001;Querner, 2000) and did not cause any surface saturation, except for a few nodes at the hillslope outlet that represent the river bank and bottom.The distribution of the hydraulic head values at the end of this simulation was used as the initial condition for subsequent simulations with solute transport.We simulated three different scenarios with solute transport: scenario 1 represents steady-state conditions in order to analyze the pattern of isotopic enrichment of the pesticide during transport and degradation under average hydrological conditions; scenario 2 focuses on the response of concentrations and isotope ratios to an extreme rainfall event to determine the effect of surface runoff; and scenario 3 incorporates periods of baseflow conditions and extreme rainfall events to study transient pesticide concentrations and isotope ratios in the course of the year.

Scenario 1: steady-state flow conditions
The first scenario was designed to mimic diffuse pollution under average hydrological conditions.The emission source was implemented as a specified concentration boundary with a constant relative concentration of C 0 = 1.0 for the sum of the light and heavy isotopes of each isotopic element.In addition, one conservative tracer with C 0 = 1.0 was applied across the entire surface of the model domain and another tracer (also with C 0 = 1.0) at the application area only, to allow for the calculation of the mean travel time of groundwater.This scenario was run until the concentrations at the hillslope outlet reached steady state.

Scenario 2: extreme rainfall event
As a hillslope system is in reality exposed to varying hydrological conditions, it was subsequently studied how concentrations and isotope ratios responded to a single extreme rainfall that leads to surface runoff.To facilitate the occurrence of surface runoff, the coupling length for this scenario was increased from 0.1 to 0.8 m.Rainfall was applied with a uniform intensity of 60 mm h −1 for 30 min.Based on a rainfall depth-duration-frequency curve (Overeem et al., 2009), the return period for such an event in the Netherlands is more than 58 yr.For the rest of the simulation time, the recharge rate was held constant at the same value as for the steadystate scenario (scenario 1).Solute transport was initialized with the concentration results from the steady-state simulation; the concentration in the source area was kept at a constant value of C 0 = 1.0 for the total concentration of each isotopic element.

Scenario 3: transient simulation for water flow and solute input
Finally, a time-varying hydrological system was simulated by applying the daily rainfall and evapotranspiration data of 2010 from the meteorological station "Twente" (Database of the Royal Netherlands Meteorological Institute, Koninklijk Nederlands Meteorologisch Instituut KNMI) for 20 yr.This time series (Fig. 2a) includes an extraordinary rainfall event at the end of August with a maximum intensity of 27.6 mm h −1 and a total rainfall amount of 106.4 mm; the rainfall event lasted for 20.8 h.Based on Overeem et al. (2009), the return period of this event exceeds 260 yr.The use of these data thus allowed for the simulation of a broad range of possible rainfall intensities and runoff responses for mid-European climates; it resulted in a net groundwater recharge of approximately 360 mm yr −1 (Fig. 2b).The annual rainfall and evapotranspiration data set was repeated for a total of 20 simulation years.This facilitated the interpretation of the results as opposed to an actual data set for 20 consecutive years, which would have produced a different hydrological and isotopic response for every simulation year.Since pesticides are not applied at a constant rate throughout the year, a more realistic setup comprises short emission pulses.To this end, the pesticide was applied once a year on a dry day in spring (day 100; 11 April) using a specified relative mass flux of 1 for the sum of light and heavy isotopes of each isotopic element for every grid cell at the pesticide source.The specified mass flux was chosen as the boundary condition for this scenario because the pesticide application is supposed to occur on a day without precipitation.If the boundary condition had been set to a specified concentration, the pesticide would not have been transferred to the subsurface because of the lack of infiltration on dry days.The model domain did not contain any pesticide at the beginning of the simulation.In order to assure the presence of pesticide in the system, which is required for the simulation of degradation, the model run was started on the first application day.

Post-model calculations
For the steady-state conditions, the cumulative transit time distribution at the hillslope outlet was derived from the concentration of the conservative tracer that was applied to the entire surface of the model domain using the relation between the transit time distribution and conservative solute breakthrough described by Duffy and Lee (1992) and further explored by Eberts et al. (2012).The mean travel time (MTT) of the system was calculated by integrating the transit time distribution.The MTT for the pesticide at steady state was obtained from the concentration at the outlet of the conservative tracer that was only applied at the pesticide application area.
Two-dimensional CSIA permits the assessment of the relative contribution of each transformation pathway to the overall degradation in a system of two competing reactions.For the simulated system, which involves an aerobic and an anaerobic reaction affecting carbon and hydrogen isotope ratios, this contribution was calculated according to van Breukelen (2007b): where the subscripts C and H indicate the enrichment factors for carbon and hydrogen isotopes, the subscripts aerobic and anaerobic denote the respective pathway, and is the ratio of the isotopic shifts 2 H and 13 C for carbon and hydrogen: = 2 H/ 13 C.The isotopic shift, i.e. the change in the δ value, is defined as where δ sample and δ source are the δ values for the carbon or hydrogen isotopes at a point in the model domain and at the pollution source, respectively.The non-degraded fraction of the compound at a point in the model domain, f deg , was obtained by substituting F (Eq. 6) into the Rayleigh equation for two-dimensional CSIA (van Breukelen, 2007b): The fraction f deg yields the extent of degradation based on the Rayleigh equation approach: The exact value of the fraction of non-degraded pesticide, f deg,m , was calculated from the simulated concentrations of the pesticide and the conservative tracer, and was used to quantify the underestimation of the extent of degradation that would result from the application of the Rayleigh equation to the simulated isotope data (Abe and Hunkeler, 2006).Given that the rate constant, k, of a first-order reaction satisfies −kt = ln(f ), this underestimation, θ, was calculated as (van Breukelen and Prommer, 2008) The simulation results were also used to derive the residual fractions of the pesticide that remain at the hillslope outlet  8) to the simulated isotope ratios at the hillslope outlet.The dilution factor f dil was then obtained as where f tot is the ratio between the concentration at the hillslope outlet and the source concentration.

Scenario 1: steady-state flow conditions
The constant recharge rate in scenario 1 produced a steadystate flow and transport regime with a mean travel time (MTT) of 6.7 yr for the groundwater and of 5.0 yr for the pesticide.At the hillslope outlet, this resulted in a steady-state concentration of 0.09 for the degrading pesticide and 0.41 for the conservative tracer that was applied at the pesticide application area relative to the initial concentration of 1.0 at the pollution source.Pesticide concentrations were lower than concentrations of the conservative tracer in the entire model domain.They were highest below the application area and decreased with depth and distance from the source (Fig. 3a).
In contrast to the shallow subsurface beneath the application area, the shallow soil layers at the hillslope bottom were characterized by low solute concentrations.
Figure 3b shows the carbon isotope ratios in the hillslope domain, which became progressively enriched with distance from the source area.With a shift of up to 4 ‰ from the initial value of −30 ‰, the shallow subsurface at the lower hillslope was characterized by the largest enrichment in δ 13 C. Hydrogen isotopes were enriched by about 20 ‰ in the shallow soil layers close to the hillslope outlet (Fig. 3c).In contrast to carbon isotopes, the strongest fractionation effects occurred in the deeper bedrock, which showed an enrichment of up to 25 ‰.
Degradation during transport induced isotope fractionation, which resulted in a steady-state isotope ratio of −26.7 ‰ for δ 13 C and −80.8 ‰ for δ 2 H at the hillslope outlet.According to the two-dimensional Rayleigh equation (Eq.8), this corresponds to an extent of degradation of 73 % at the hillslope outlet.It follows from Eq. ( 11) that the nondegraded fraction was approximately 27 % (f deg = 0.27 in Eq. 11), while the residual fraction after the effect of dilution was 34 % (f dil = 0.34 in Eq. 11).Dilution therefore contributed slightly less to the concentration decrease than degradation.The extent of degradation, derived from the simulated isotope values (Eq.9), increased with depth and distance from the pollution source (Fig. 3d).With a value of 80 %, it was highest in the deep bedrock and in the shallow soil layers at the footslope.The θ value, which compares the extent of degradation given by the two-dimensional Rayleigh equation to the model-based extent of degradation (Eq.10), was largest below the application area (up to 25 %).Apart from this area, the underestimation was < 20 %; it was approximately 5 % in the deep bedrock (Fig. 3f) and 11.5 % at the hillslope outlet.
Figure 3e shows the relative contribution F of the aerobic reaction to the overall degradation of the pesticide according to Eq. ( 6).F decreased rapidly with depth from a fraction of more than 0.9 in the topsoil below the pesticide application area to less than 0.2 in the deeper bedrock.At the hillslope outlet, it reached a value of 0.39.

Scenario 2: extreme rainfall event
Figure 4a illustrates the pesticide concentrations and mass flux at the hillslope outlet in response to the extreme rainfall event.It can be seen that with the onset of rainfall, concentrations first dropped and subsequently reached a distinct maximum.In contrast, the pesticide mass flux shows two peaks, with the second peak being much more pronounced than the first.In response to the rainfall event, the isotope signatures at the hillslope outlet dropped from their steady-state values of δ 13 C = −26.7 ‰ and δ 2 H = −80.8‰ to δ 13 C = −30.0‰ and δ 2 H = −100.0‰ (Fig. 4b).The latter values are the characteristic isotopic signatures of the pollution source.The minima in the isotope ratios coincided with the maximum in pesticide concentrations and the second peak in pesticide mass flux (Fig. 4a).

Scenario 3: transient simulation for water flow and solute input
The changing hydrological conditions and the annual pesticide application in scenario 3 induced a transient concentration response in the subsurface.We consider the first fifteen years of the transient simulation as spin-up period and focus on the results of the last five simulation years.Figure 5 shows a comparison of the pesticide concentration at two different times: after 19 yr (panel a) and after 19 yr and 50 days (panel b), thus 50 days after the last pesticide application.It illustrates that the pesticide plume from the new application mixed with the residual plume from the previous applications.
The inclusion of evapotranspiration into the simulation resulted in several periods of negative net infiltration, which followed the time of pesticide application (shaded areas in Fig. 6).Hereafter, this phase will be denoted as the dry period.The dry period was characterized by an increase in pes- ticide concentrations till the extreme rainfall event in August (Fig. 6b).Two other concentration peaks occurred later in each simulation year (beginning of November and end of February), while concentrations sharply dropped following increased rainfall in November.After the third concentration peak in February, the pesticide concentrations decreased again until the next pesticide application.In contrast to the steady-state scenario, the pesticide concentrations at the hillslope outlet thus varied significantly and followed the same seasonal pattern every year.In particular, some of the days with high rainfall induced distinct minima in the pesticide concentrations.
Figure 6c shows the δ 13 C-and δ 2 H-values, which alternated between periods of enrichment and depletion.Moreover, despite the 15 yr spin-up period, the isotope ratios still increased slightly.The hydrogen isotope ratios declined to their minimum in the middle of the wet season and reached the highest value at the end of the dry season (shaded area).In contrast, the peak in carbon isotope ratios occurred at the onset of the dry season.While the sharp drops in pesticide concentrations correlate with the minima in δ 2 H, they correspond to maxima in δ 13 C.
The contribution of the aerobic degradation to overall degradation (F ) mirrored the seasonal pattern of the δ 13 Cvalues: the maximum occurred at the beginning of the dry season (Fig. 6d).In contrast, the extent of degradation (B) increased with increasing concentrations during the dry season and reached its maximum at the end of the dry season before the extreme rainfall event in August and subsequent precipitation.It thus shows a similar response as the δ 2 Hvalues (Fig. 6c).In addition, the parameters B and F reflect the response of the isotope ratios to rainfall events: while the distinct peaks in δ 13 C-values correspond to maxima in F , the drops in δ 2 H-values coincided with minima in B.
The underestimation of the extent of degradation (θ ) that results from the use of the two-dimensional Rayleigh equation approach increased during the wet period and with subsequent simulation years; it reached a value of 9.2 % at the hillslope outlet at the end of the transient simulation (Fig. 6e).In addition, θ shows distinct peaks that are correlated with the concentration minima during rainfall events.With respect to the entire model domain, the maximum in θ occurred in the shallow subsurface below the application area in response to every pesticide input.Subsequently, this zone of the highest θ moved downgradient with the plume center (not shown).

Pesticide movement through the hillslope
Scenario 1 yielded an average travel time of five years for the pesticide, which implies a relatively long response time to the emission of diffuse pollutants.This is consistent with elevated atrazine concentrations in a spring in a small catchment several years after the use of the herbicide had been abandoned (Gutierrez and Baran, 2009).Similarly, the pesticide reached the hillslope outlet approximately 1.4 yr after the first application in the spin-up period of the transient scenario (scenario 3).This, in turn, highlights the potential time lag between the first use and the first detection of a newly introduced pesticide in stream flow.Analogous to the continuous detection of an abandoned compound, the extent of contamination caused by a new compound might manifest itself only after several consecutive application years.Moreover, the isotope ratios and derived parameters in the transient scenario did not reach an oscillatory steady state at the end of the 20 yr simulation.Correspondingly, the extent of degradation (B) still increased over the course of the last five simulation years (Fig. 6d).This is another indicator of the long response time of the modelled hydrological system to diffuse pollution.
At steady state, the pesticide concentrations in the shallow soil layers at the footslope were low (Fig. 3a), which resulted from their position above the water table during steady-state conditions: the transfer of pesticides into the unsaturated shallow soil close to the outlet could only occur via dispersion and diffusion from deeper layers, whereas the transport through the saturated zone was driven by faster advective transport.The part of the pesticide that was transported to the shallow subsurface at the footslope had, therefore, undergone more degradation than the pesticide that was present directly below the water table (Fig. 3d).Correspondingly, the shallow subsurface soil at the footslope also shows more enriched isotope ratios relative to other parts of the model domain (Fig. 3b and c).
The first two periods of increasing concentrations in each simulation year in transient scenario 3 ended abruptly in response to rainfall events (Fig. 6b).These rainfall events thus led to a pronounced dilution effect in the pesticide concentrations.In contrast, the third concentration peak in every simulation year did not precede a large rainfall event, but was related to a transition to wetter conditions (Fig. 6a), which induced a stronger dilution of pesticide concentrations.In summary, concentrations seem to be rather driven by the hydrological conditions than by the application pattern of the pesticide.However, the application pattern might have influenced the second phase of rising concentrations following the extreme rainfall event because the pesticide reached the hillslope outlet 1.4 yr after the first application.The second increase thus coincided with the arrival of the pesticide plume from the previous year.Consequently, the hydrological conditions immediately affected the concentration pattern, whereas the response to the emission of the contaminant into the system was much more delayed and subdued.
The results of the steady-state simulation illustrate how the relative contribution of the aerobic degradation pathway (F ) rapidly decreased with depth (Fig. 3e).This is due to the fact that aerobic conditions solely occurred in the topsoil.Accordingly, flow lines that display long travel times through the bedrock exhibit a significant extent of degradation under anaerobic conditions.However, even directly below the application area, F did not attain its maximum value of one because of the effect of vertical dispersion that caused pesticide transport from the anaerobic subsoil into topsoil layers.F reached a steady-state value of 39.0 % at the hillslope outlet, although the topsoil accounts for less than 3 % of the total subsurface volume.This mainly resulted from the higher degradation rate in the topsoil (5 per year) than in the subsoil and bedrock (0.2 per year), which disproportionately increased the contribution of the aerobic reaction to overall degradation.Furthermore, the pesticide traversed the topsoil (i.e.aerobic zone) during infiltration and exfiltration, which enhanced the imprint of the aerobic degradation pathway.Correspondingly, the relatively small influence of the anaerobic transformation indicates that a considerable part of the pesticide was transported via shallow flow pathways, which did not traverse the deeper bedrock and thus had a relatively short residence time in the anaerobic zone.The significant contribution of shallow flow pathways to hillslope discharge resulted, in turn, from the decrease in saturated hydraulic conductivity with depth.Hence, the chosen hydraulic properties also influenced the relative contribution of the aerobic reaction pathway to overall degradation.

Carbon and hydrogen isotope ratios
In general, isotope ratios increased with travel distance, as more degradation leads to more pronounced isotope fractionation.This becomes apparent when considering the progressive enrichment in carbon and hydrogen isotope ratios with distance from the source area during steady state (Fig. 3b  and c).However, carbon and hydrogen isotopes show a different strength of enrichment in the deep bedrock, which can be ascribed to the choice of enrichment factors for the aerobic and the anaerobic zones.In contrast to the hydrogen isotopes, the isotopic enrichment of carbon was stronger for the aerobic reaction pathway in the topsoil than for the anaerobic reaction pathway in the subsoil and bedrock.Hence, the pesticide that crossed the topsoil for the second time during exfiltration shows the most pronounced enrichment in δ 13 C.In addition, this re-entrance into the topsoil was driven by a dispersive flux, which is, as explained before, associated with a higher degree of degradation than the advective flux.In contrast, the most enriched δ 2 H-values occurred in the deep bedrock, since long travel times through the subsurface induced, together with the high hydrogen enrichment factor under anaerobic conditions, the strongest fractionation effect for the hydrogen isotopes.Corresponding to these different areas of largest isotope fractionation for carbon and hydrogen, the two-dimensional Rayleigh equation (Eq.8) yielded a high degree of degradation for the deep bedrock and the shallow soil layers at the footslope (Fig. 3d).
In scenario 3, the overall extent of degradation increased with longer residence times of the pesticide in the hillslope (i.e. with pesticide transport via deeper flow pathways).Therefore, B increased under drier conditions and was highest at the end of the dry season (Fig. 6d).The concentration increase during the dry season (Fig. 6b) thus resulted from less dilution due to little rainfall, and not from less degradation.Similar to scenario 1, the choice of enrichment factors also explains the differences in the seasonal patterns of the hydrogen and carbon ratios in scenario 3. Since the enrichment for hydrogen was assumed to be stronger in the anaerobic zone than in the aerobic zone, the increase in δ 2 H during the dry season reflects the increasing contribution of deeper groundwater.This agrees with the shape of the curve of the extent of degradation (Fig. 6d), which resembles more the pattern of δ 2 H than the pattern of δ 13 C.In general, an increasing contribution of deeper flow pathways, which implies longer travel times, also caused a higher degree of fractionation for the carbon isotopes.Nevertheless, a smaller relative contribution of the aerobic reaction pathway can account for a less pronounced overall fractionation effect for carbon because of the higher enrichment factor in the topsoil.Accordingly, δ 13 C-values increased during the wet season as a result of more degradation in the topsoil (Fig. 6d), which demonstrates that shallow flow through the topsoil was more pronounced under wet conditions.This finding agrees with Brown et al. (1999) and Rozemeijer and Broers (2007), who emphasized the increasing importance of shallow subsurface flow during rain events in headwater catchments and lowlands, respectively.
The simulation of isotope ratios allowed for quantification of the extent of pesticide degradation (B).It was thus possible to compare the concentration reduction by dilution to the concentration reduction by transformation processes.In addition, the Rayleigh equation approach yielded the distribution between the aerobic and anaerobic reaction pathway (F ), which suggests that isotope data of two elements can, under favorable conditions, allow for the distinction between two reaction mechanisms.Moreover, F was used for the analysis of the relative importance of shallow and deeper flow pathways to pesticide fluxes to the stream.In summary, these analyses support the additional benefit of CSIA, as pesticide concentration data alone would not have allowed for the derivation of B and F in a real hillslope system.

Responses to rainfall events
The response of concentrations and isotope ratios to the extreme rainfall event in scenario 2 highlights the advantages of combined concentration and CSIA data.Judging by the peak in pesticide mass flux, the concentration minimum at the onset of the rainfall event has to be ascribed to initial dilution by rainwater (Fig. 4a).The subsequent concentration peak occurred shortly after the rain had ceased.Given that the mass flux shows a second maximum, this concentration peak stemmed from increased pesticide transport.However, it cannot be concluded from the concentration and mass flux data which mechanism led to increased pesticide transport.In contrast, the analysis of isotope ratios at the hillslope outlet reveals the underlying mechanism.During the rainfall event, the isotope ratios remained at the constant pre-event level despite the occurrence of surface runoff.This indicates that the surface runoff reaching the hillslope outlet did not contain any pesticide yet, which is due to the uphill location of the pollution source and the pesticide-free area at the lower hills-lope.Subsequently, the isotope ratios decreased to the source values of −30 ‰ and −100 ‰, respectively (Fig. 4b).This resulted from the arrival of contaminated surface runoff at the hillslope outlet, which occurred, by coincidence, at the end of rainfall.As the time lag between the release from the pollution source and the arrival at the hillslope outlet was too short (less than 30 min) to allow for significant degradation, the pesticide did not undergo any detectable fractionation.It follows that transport via surface runoff dominated the overall pesticide flux to the stream at that time.This scenario thus shows how the combined analysis of concentrations and isotope data can reveal the occurrence of surface runoff.
The first peak in mass flux prior to the arrival of contaminated surface runoff suggests that the rainfall event initially caused enhanced pesticide transport via exfiltration of shallow groundwater.This illustrates the concept of emission of "old" contaminant residues with pre-event water at the onset of rain (Burt and Pinay, 2005).Discharge of pre-event water can also be inferred from the analysis of the parameter F , which indicates a slightly enhanced contribution of topsoil degradation after the onset of rainfall and prior to the discharge of contaminated surface runoff (not shown).
Surface runoff in scenario 2 occurred as infiltration excess overland flow, which can be an important mechanism of pesticide transport from agricultural land (Doppler et al., 2012).The infiltration excess overland flow was generated by increasing the coupling length from 0.1 m to 0.8 m, which decreased the coupling between overland flow domain and subsurface domain.This can be considered an analogue for the decreased infiltration capacity of the soil as a result of surface sealing during high-intensity rainfalls.The same approach was applied in the HGS simulation by Verbist et al. (2012) to generate infiltration excess overland flow.A layer of low hydraulic conductivity at the surface would have had the same effect, but this would have required a modification of the model domain.Alternatively, the rainfall intensity could have been increased, but without changing the coupling length, a sufficiently long period of overland flow could have only been achieved with an unrealistically high rainfall amount (data not shown).
With a return period of nearly 60 yr, the simulated rainfall event in scenario 2 corresponds to an exceptional event for mid-European climate.This extreme rainfall event was required to generate continuous overland flow from the area of the pesticide application to the hillslope outlet and thereby produce a discernible response in the isotope ratios in the stream.Nonetheless, scenario 2 can be considered representative of situations where the surface is sealed or has been disturbed (see above), or where a pesticide product is spilled on impermeable areas such as paved farmyard.In this case, subsequent surface runoff from these contaminated areas may lead to concentration peaks during relatively dry periods (Holvoet et al., 2007;Kreuger, 1998).Similarly, preferential flow to drainage systems represents a fast transport route of pesticides that can lead to elevated pesticide concentrations in the absence of surface runoff at the application area (David et al., 2003).Pesticide transport via drain flow and surface runoff from impermeable areas can thus occur in response to rainfall of lower intensity.The response in isotope ratios would be comparable to the one simulated in scenario 2, and has already been observed for oxygen isotope ratios of nitrate: Ging et al. (1996) found elevated δ 18 O-values in a storm sewer, which were indicative of an atmospheric nitrate source, and ascribed them to surface runoff from impervious areas.Correspondingly, CSIA data that is measured during rain events could also be used for forensic source determination, i.e. for the distinction between different pesticide products and application areas that are characterized by different isotopic compositions.
The sharp drops in pesticide concentrations in response to large rain events in scenario 3 (Fig. 6b) resulted from dilution.Moreover, these drops also occurred on days with smaller events, provided that the rainfall (in combination with previous wet conditions) was sufficient to saturate the footslope.This had a similar dilution effect as high intensity rainfalls.The concentration minima correlate with minima in hydrogen isotope ratios, as high intensity rainfall and footslope saturation were accompanied by the discharge of more shallow groundwater, which had a shorter residence time and thus shows less hydrogen isotope fractionation.In contrast, the concentration minima correspond to peaks in the carbon isotope ratios, since they are associated with a relatively large contribution of topsoil degradation and thus more pronounced carbon isotope fractionation.Consequently, analogous to the seasonal pattern of isotope ratios, the different strength of isotope fractionation for the aerobic and the anaerobic reaction pathway resulted in the opposite behavior of carbon and hydrogen isotope ratios during single rainfall events.This demonstrates how CSIA can give insights into transport routes of a contaminant in a hydrological system, provided that reaction pathways can be attributed to different zones.
The rainfall events in scenario 3 that led to footslope saturation were accompanied by minima in the extent of degradation (Fig. 6d), since they induced discharge of more recently applied and thus less degraded pesticide.This agrees with previous findings about nitrate contamination where rainfall events led to a fast mobilization of soil nitrate with a relatively depleted δ 15 N-value, while increasing denitrification in soil between rainfall events resulted in progressively enriched δ 15 N-values (Kellman and Hillaire-Marcel, 2003).In contrast to the simulation of the extreme rainfall event in scenario 2, the isotope ratios did not drop to the initial values of the pesticide application area in scenario 3, which implies the absence of pesticide transport via surface runoff.This can be explained by the time lag of 137 days between the pesticide application and the extreme rainfall event in August, which allowed for infiltration and degradation of the pesticide prior to intense rainfall as opposed to the simultaneous application of rain and pesticide in scenario 2.

Validity of model assumptions
As shown in Sect.4.1, the relative importance of the aerobic and anaerobic reaction pathways (F ) allowed for the analysis of the role of shallow and deeper flow pathways for pesticide transport.This analysis was facilitated by the model assumption of mutually exclusive reaction pathways with distinct enrichment factors that were active in specific hillslope layers.We considered this a reasonable assumption for aerobic and anaerobic degradation.In reality, however, reaction mechanisms might not be spatially exclusive, or their spatial distribution might be unknown.The analysis of transport routes based on F might therefore be restricted to cases where the reaction mechanisms occur in specific subsurface zones (e.g., anaerobic beneath the aerobic zone), and would be more complex for competing degradation pathways that occur simultaneously in space, entail a similar extent of isotope fractionation, or whose spatial distribution is unknown.
The simulations did not explicitly account for pesticide sorption, volatilization, or preferential flow pathways.Volatilization was not considered, as it would mainly lead to a decrease in the pesticide mass load at the pollution source.Moreover, as rainfall occurred soon after pesticide application in the simulations, we assumed a rapid mobilization of the applied compound and thus minor losses due to volatilization for the modelled system.This implies negligible volatilization-induced isotope fractionation under comparable hydrological conditions.However, depending on the properties of the pesticide, soil parameters, and meteorological conditions, volatilization losses can be up to 50 % of the applied amount for some pesticides (van den Berg et al., 1999).If isotope fractionation at the pollution source due to volatilization is not considered in this case, CSIA might yield an inaccurate assessment of the extent of degradation.
Sorption to the solid phase can be relevant for pesticide transport in two ways: as sorption to suspended matter (e.g. in surface runoff), and as sorption to the soil matrix.Calculations under the assumption of equilibrium sorption revealed that transport with suspended matter in surface runoff is only relevant for highly sorbing compounds (log K OC > 4; not shown).We therefore decided to disregard this transport route in our simulations.Simulation of sorption to the soil matrix would have had a retardation effect on pesticide transport to the hillslope outlet.If sorption had been included, it would have been possible to also simulate sorption-induced isotope fractionation in the soil matrix.This can result in an enrichment in heavy isotopes at the front and a depletion at the tail of a migrating groundwater plume (Kopinke et al., 2005;van Breukelen and Prommer, 2008).If isotopic enrichment is only attributed to degradation, CSIA can thus yield an overestimation of degradation for the plume front.However, we did not explicitly simulate sorption-related isotope fractionation, as this is only relevant in non-stationary plumes when degradation is slow (van Breukelen and Prommer, 2008), and would thus be insignificant in the modelled hillslope system.We also disregarded diffusion-induced isotope fractionation effects, as their importance for the simulated system should be minor given the widespread nature of the emission and the relatively large spatial scope (van Breukelen and Rolle, 2012).
Similar to other virtual experiments with physically based models (e.g.Hopp andMcDonnell, 2009, 2011;James et al., 2010;Mirus et al., 2011;Mirus and Loague, 2013), we did not include preferential flow in the simulations.We anticipate that vertical preferential flow leads to a faster transition of water through the soil, and thus a decrease in the extent of degradation in the soil layers.This will, in turn, increase pesticide concentrations and decrease the isotope ratios at the hillslope outlet.Furthermore, vertical preferential flow would result in a lower contribution of the aerobic reaction pathway to overall degradation.The simulation of lateral preferential flow in the soil would decrease flow to the bedrock and result in a faster subsurface flow response to rainfall.This would, similar to vertical preferential flow, lead to less degradation and isotopic enrichment.However, the effect of lateral preferential flow pathways on the relative contribution of aerobic and anaerobic degradation to overall degradation, and thus on carbon and hydrogen isotope ratios at the hillslope outlet, depends on their location in the soil.
Even if the influence of preferential flow on concentrations and isotope fractionation is small, it may still affect the accuracy of the CSIA method, as it would cause enhanced mixing between recently applied pesticide and partially degraded pesticide in the soil.This might result in a stronger attenuation of apparent isotopic enrichment, and thus amplify the underestimation of degradation by the Rayleigh equation approach (Kopinke et al., 2005).
CSIA at downstream river monitoring points might be affected by polluted influent water from upstream parts, which could mask isotope fractionation effects in hillslope discharge at the stream monitoring point and, therefore, bias the quantification of pesticide degradation.Similarly, dry conditions can cause infiltration of polluted river water into the streambank or aquifer, which would perturb the pattern of gradual enrichment with increasing travel time in the hillslope.Hence, the simulation results primarily apply to gaining streams in headwater catchments or river sections that do not show mixing with residues of the same pesticide from upstream sources.

Implications for the applicability of CSIA to assess pesticide transport and transformation
In order to examine the applicability of compound-specific isotope analysis in the context of diffuse agricultural pollutants, the concentration decrease in scenario 1 needs to be compared to the difference between soil water concentrations and detection limits of CSIA.Initial soil water concentrations of up to a few milligrams per liter have been reported for pesticides at realistic application rates (Liu et al., 2012), while carbon isotopes of organic contaminants, including the pesticide atrazine and its metabolite desethylatrazine, could be measured at concentrations of around 100 ng L −1 (Jochmann et al., 2006;Schreglmann et al., 2013).Consequently, the difference between pore water concentrations and the detection limit for carbon isotope analysis can be four orders of magnitude.By way of comparison, the simulation yielded a concentration reduction by a factor 10 between the application area and the hillslope outlet.Therefore, the model results indicate that, given appropriate sampling and preconcentration techniques, low environmental concentrations of pesticides would not impede CSIA of diffuse river pollutants.
The magnitude of the simulated enrichment between the source area and the hillslope outlet in scenario 1 exceeds the uncertainty range of CSIA.Isotope ratios increased by 4 ‰ for δ 13 C and by 20 ‰ for δ 2 H, while the instrumental uncertainty is about 0.5 ‰ for carbon and 5 ‰ for hydrogen isotope analysis (Sherwood Lollar et al., 2007).This indicates that the enrichment in the isotopic composition at the hillslope outlet relative to the source values is detectable, which supports the use of CSIA in the analysis of diffuse river pollution under average hydrological conditions.We suggest that isotope ratios during baseflow conditions indicate the maximum potential for degradation, as they are not influenced by fast pesticide transport routes such as surface runoff (see scenario 2).It should be noted, however, that the magnitude of isotope enrichment depends, among others, on the site-specific travel times and isotope fractionation effects.
In contrast to the response to the extreme event in scenario 2, the isotope ratios in scenario 3 display only small variations, which did not exceed 0.3 ‰ for δ 13 C and 1.4 ‰ for δ 2 H during each simulation year.Hence, these small fluctuations would not be detectable based on the current accuracy of carbon and hydrogen CSIA.In the absence of overland flow, the assessment of the extent of in situ degradation on the basis of CSIA would thus yield the same results throughout the year.This suggests that grab samples for CSIA measurements are sufficient for a representative assessment of pesticide transformation between the pollution source and river monitoring point, except during events that result in direct transport via overland flow.This might, however, not apply to systems that show a larger temporal variation in the relative contribution of baseflow versus stormflow to streamflow, which could, for example, result from the activation of preferential flow pathways during rainfall events or a larger difference between the permeabilities of the soil and bedrock.A more pronounced variation in in-stream CSIA might also arise from pesticide use closer to the river than for the modelled hillslope, as rainfall events would then be more likely to lead to transport of recently applied pesticide into the river, while the dry season would be characterized by transport of strongly degraded pesticide via baseflow.Moreover, catchments with well-drained soils, long mean travel times and a high groundwater contribution to streamflow tend to dampen solute input signals to a much larger extent than catchments that exhibit more responsive soils and short mean transit times (Hrachowitz et al., 2013).Consequently, the latter systems might show a significant seasonal variability in isotope data that reflects the short emission pulses of pesticide application and the fast subsurface transport mechanisms in response to rainfall events.The determination of factors that would cause larger seasonal variabilities in instream CSIA was not the objective of this study, but should be addressed in future research.
During extreme rainfall events, only a fine temporal sampling resolution can allow for the detection of transient surface runoff, as the latter might only occur for a short time.According to the comparable response in concentrations and isotope ratios in scenario 2, the detection of surface runoff by CSIA requires the same temporal resolution as concentration measurements.Analysis of isotope ratios in addition to concentrations during rain events requires the use of automatic sampling devices, which would also allow capturing surface runoff from potentially contaminated areas.If a concentration peak in the monitoring data occurs, CSIA could then additionally be performed in order to detect isotope ratios that correspond to typical values for pesticide products.This analysis can, however, only allow for the detection of contamination via surface runoff if the isotopic signature of the diffuse pollution source is known and distinguishable from isotope ratios that are associated with pesticide transport via subsurface flow pathways.
The model results for the isotope ratios allowed for a detailed analysis of the seasonal pattern of the extent of degradation (B), the relative contribution of the aerobic reaction pathway (F ), and the underestimation by the use of the Rayleigh equation approach (θ).In the case of measured CSIA data, these parameters might, however, appear constant over time because of minor seasonal fluctuations in isotope ratios within analytical uncertainties.Nonetheless, even a constant value of these parameters can facilitate the analysis of the underlying transport and transformation mechanisms in the studied flow system.Since CSIA represents a unique method for the determination of these parameters, this highlights an additional benefit of isotope analysis in the context of diffuse river pollutants.
The extent of degradation during pesticide transport through the hillslope was determined by applying the Rayleigh equation to the simulated isotope ratios.Nevertheless, the Rayleigh equation is in principle only applicable to closed and fully mixed systems (e.g. to degradation experiments in microcosms) (van Breukelen and Prommer, 2008).Hydrological systems such as the modelled hillslope are, however, open systems, which display a variety of transport and transfer processes.For example, dispersion leads to different transport routes with varying flow velocities to the same measurement point.Therefore, while some molecules might have been subject to strong isotope fractionation during the transport along a flow pathway, others might have travelled much faster via a different pathway.The latter molecules exhibit less fractionation because of less exposure to degradation processes (Aravena and Hunkeler, 2009).In other words, physical heterogeneities of the aquifer and hydrodynamic dispersion can result in the attenuation of isotope fractionation effects that occurred in a part of the sample mixture only (van Breukelen and Prommer, 2008).This masking effect might also appear in areas where flow pathways from different layers mix (Kopinke et al., 2005).Consequently, if the Rayleigh equation approach is applied in such open flow systems, it tends to result in an underestimation of the extent of in situ degradation (Abe and Hunkeler, 2006;van Breukelen and Prommer, 2008).For the simulated hillslope system, the extent of degradation at steady state was significantly underestimated for the topsoil below the application area.This is due to the high degradation rate under aerobic conditions, which induced larger concentration gradients and thus larger dispersive fluxes than the slower reaction rate under anaerobic conditions.Large dispersive fluxes cause, in turn, an attenuation of isotopic shifts and, therefore, result in a more significant underestimation of the actual extent of degradation (van Breukelen and Prommer, 2008).Correspondingly, the zone of the strongest underestimation in scenario 3 was located at the plume center, as the pesticide plume is associated with the strongest concentration gradients in the system (not shown).
Although transient hydrological conditions are likely to cause enhanced mixing of different flow pathways, the degree of underestimation at the hillslope outlet in the steadystate simulation was slightly larger than for the transient scenario.This might result from a greater importance of deeper flow pathways to pesticide transport at steady state, which agrees with a smaller contribution of the aerobic pathway at steady state than in the transient scenario.Pathways with a long travel time thus controlled the isotope ratios at the hillslope outlet to a larger degree in scenario 1 than in scenario 3. Consequently, the contribution of pesticide that displays advanced degradation and thus strong isotope fractionation to overall pesticide export was more significant at steady state.This, in turn, favored the masking of strong isotope fractionation in scenario 1, which resulted in a higher θ for scenario 1 compared to scenario 3. Scenario 3 shows a higher θ during the wet period than during the dry period (Fig. 6e), which illustrates the effect of mixing of different flow pathways on the accuracy of the Rayleigh equation: the growing relative contribution of more shallow flow pathways during wet conditions resulted in increased discharge of more recently applied pesticide, which led to an enhanced masking of the isotopic enrichment in deeper groundwater.The same mechanism caused the distinct peaks in θ following rainfall events.
The underestimation of degradation due to the use of the Rayleigh equation is typically below 5 % for groundwater plumes (Abe and Hunkeler, 2006).However, it has been shown that the underestimation can exceed 50 % at fringes of pollution plumes, especially for high degradation rates and at large distances from the contaminant source (van Breukelen and Prommer, 2008;van Breukelen and Rolle, 2012).In view of the large underestimation for groundwater systems, the maximum θ value at the hillslope outlet of 11.5 % and 10.0 % in scenario 1 and scenario 3, respectively, can be considered negligible.Therefore, the simulation results suggest that CSIA yields a good assessment of in situ degradation, not only for aquifer systems, but also for more complex subsurface-surface systems that are subject to mixing of different flow pathways.

Conclusions
In this paper, we present a model study of compound-specific stable isotope analysis (CSIA) in the context of diffuse pollution.The objective was to examine whether CSIA qualifies as a feasible and expedient technique for the analysis of transport pathways and the assessment of the extent of degradation of diffuse pollutants.We simulated reactive solute transport and isotope fractionation effects for a hypothetical hillslope.The model results support the usefulness of CSIA data in this context: the simulated isotope data allowed for the quantification of the extent of in situ degradation and the relative contribution of two competing pathways to overall degradation, which would not have been possible on the basis of simulated concentration data only.
The two-dimensional Rayleigh equation provided a reliable estimate of the overall extent of degradation under transient conditions.In particular, the inherent underestimation of the Rayleigh equation approach was small, considering the high degree of mixing of groundwater flow pathways from different depths at the hillslope toe.The attenuation of isotope signals, which partly results from this mixing, did not exceed the degree of attenuation reported in previous studies of groundwater pollution plumes.
The simulation of an extreme rainfall event illustrated how isotope data can, as opposed to concentration data alone, reveal the occurrence of surface runoff and thereby indicate the fast transport of a diffuse pollutant to a river.In this way, CSIA might allow for the distinction between pollution via surface runoff or direct spillage, and via groundwater exfiltration solely.However, the simulation results also showed that surface flow might only be discernible in CSIA data for a very short period, which requires the use of automated sampling procedures during large rainfall events.
The simulation of transient hydrological conditions resulted in small seasonal variations in isotope ratios and derived parameters, which would not be detectable in CSIA data because they would fall within the uncertainty range of current analytical methods (with the exception of pesticide transport via surface runoff in response to rain events).For systems with a larger seasonal variation in isotope ratios (e.g.resulting from shallower and more permeable soils or the activation of preferential flow pathways in response to rain events), CSIA could yield a time-dependent estimate of the extent of degradation.In the case of a system with low seasonal variability such as the modelled hillslope, however, CSIA would give a stable result throughout the year, regardless of the temporal sampling resolution.This, in turn, supports the feasibility of CSIA in the analysis of pollutants in stream flow.Provided that the relevant underlying degradation mechanisms and associated isotope fractionation factors are known, CSIA thus offers a unique tool for the assessment of pesticide transformation, and, if the spatial distribution of the degradation mechanisms is known, even for a qualitative description of the interplay between transport via shallow and deep flow pathways.
Future modelling studies might extend this study to a three-dimensional catchment and incorporate in-stream degradation, or include several pollution sources to test CSIA as a tool for source identification and apportionment.Furthermore, transport via suspended matter should be considered for highly sorptive pesticides and erosion-prone sites.In addition to modelling studies, it is crucial to further test CSIA in experimental studies, especially in view of the low environmental concentrations.These studies have not been conducted yet, but the model results advocate the applicability and advantages of CSIA.In conclusion, this study emphasized the potential benefits of CSIA in the characterization of diffuse river pollution.

Fig. 1 .
Fig. 1.Model domain with the pesticide application area (red) and three subsurface zones: topsoil (yellow), subsoil (green) and bedrock (blue) (a).The discretization is finer in the upper soil and close to the hillslope outlet (b).The hillslope outlet is represented by two boundary nodes (red circle).

Fig. 2 .
Fig. 2. Daily rainfall and evapotranspiration data for the meteorological station Twente (a); and net infiltration as the difference between precipitation and modelled actual evapotranspiration (b).The data set shown in panel a was repeated for 20 yr in the transient simulation (scenario 3).The pesticide application (11 April; day 100) is marked by the red dashed vertical line.The simulation was started on this day to ensure the presence of pesticide in the hillslope system.

Fig. 3 .
Fig. 3. Steady-state results for the pesticide concentration (a); δ 13 C-values (‰) of the pesticide (b); δ 2 H-values (‰) of the pesticide (c); extent of degradation (%, d); relative contribution of the aerobic reaction to overall pesticide degradation (e); and underestimation (%) of the true extent of degradation when the Rayleigh equation is used (f).Areas with a concentration reduction of more than three orders of magnitude relative to the source are blanked because CSIA would not be possible due to detection limits.The black arrows indicate streamlines of the steady-state flow field.The purple line shows the position of the water table.Vertical exaggeration is five times.

Fig. 4 .
Fig. 4. Response to the extreme rainfall event at the hillslope outlet immediately before and after the rainfall event: pesticide concentration and mass flux (a, note the logarithmic scale for the mass flux); and carbon and hydrogen isotope ratios of the pesticide (b).The timing of the rainfall event is shaded in grey.

Fig. 5 .
Fig. 5. Pesticide plume before and after the last application: concentrations after 19 yr (a) and after 19 yr and 50 days (b).Areas with a concentration reduction of more than three orders of magnitude relative to a concentration of 1.0 are blanked because CSIA would not be possible due to detection limits.The black line indicates the position of the water table.

Fig. 6 .
Fig. 6. Results of the last five years of the transient simulation at the hillslope outlet (after 15 yr of spin-up): application days and net infiltration (a); pesticide concentrations (b); δ 13 C and δ 2 H isotope ratios (c); extent of degradation and relative contribution of the aerobic reaction pathway to overall degradation based on the twodimensional Rayleigh equation approach (d); and underestimation of the extent of degradation resulting from the use of the Rayleigh equation (e).Periods with a negative net infiltration (dry periods) are shaded in grey.

Table 1 .
Hydraulic properties of the subsurface domain.

Table 2 .
Parameters for degradation and isotope fractionation.