Do CSIA data from aquifers inform on natural degradation of chlorinated ethenes in aquitards?

Back-diffusion of chlorinated ethenes (CEs) from low-permeability layers (LPLs) causes contaminant persistence long after the primary spill zones have disappeared. Naturally occurring degradation in LPLs lowers remediation time frames, but its assessment through sediment sampling is prohibitive in conventional remediation projects. Scenario simulations were performed with a reactive transport model (PHT3D in FloPy) accounting for isotope effects associated with degradation, sorption, and diffusion, to evaluate the potential of CSIA data from aquifers in assessing degradation in aquitards. The model simulated a trichloroethylene (TCE) DNAPL and its pollution plume within an aquifer-aquitard-aquifer system. Sequential reductive dechlorination to ethene and sorption were uniform in the aquitard and did not occur in the aquifer. After 10years of loading the aquitard through diffusion from the plume, subsequent source removal triggered release of TCE by back-diffusion. In the upper aquifer, during the loading phase, δ 13 C-TCE was slightly enriched (up to 2‰) due to diffusion effects stimulated by degradation in the aquitard. In the upper aquifer, during the release phase, (i) source removal triggered a huge δ 13 C increase especially for higher CEs, (ii) moreover, downstream decreasing isotope ratios (caused by downgradient later onset of the release phase) with temporal increasing isotope ratios reflect aquitard de- gradation (as opposed to downstream increasing and temporally constant isotope ratios in reactive aquifers), and (iii) the carbon isotope mass balance (CIMB) enriched up to 4‰ as lower CEs (more depleted, less sorbing) have been transported deeper into the aquitard. Thus, enriched CIMB does not indicate oxidative transformation in this system. The CIMB enrichment enhanced with more sorption and lower aquitard thickness. Thin aquitards are quicker flushed from lower CEs leading to faster CIMB enrichment over time. CIMB enrichment is smaller or nearly absent when daughter products accumulate. Aquifer CSIA patterns indicative of aquitard degradation were similar in case of linear decreasing rate constants but contrasted with previous simulations assuming a thin bioactive zone. The Rayleigh equation systematically underestimates the extent of TCE degradation in aquifer samples especially during the loading phase and for conditions leading to long remediation time frames (low groundwater flow velocity, thicker aquitards, strong sorption in the aquitard). The Rayleigh equation provides a good and useful picture on aquitard degradation during the release phase throughout the sensitivity analysis. This modelling study provides a framework on how aquifer CSIA data can inform on the occurrence of aquitard degradation and its pitfalls.


Introduction
Chlorinated ethenes (CEs) are widely spread toxic groundwater contaminants, due to their inadequate disposal and the spilling of the commercially widely used compounds perchloroethylene (PCE) and trichloroethylene (TCE). Their potential complete reductive dechlorination to ethene in anaerobic aquifers opened opportunities for monitored natural attenuation (MNA) as an alternative to invasive, energy-intensive technologies for polluted site remediation. MNA might be specifically relevant for the monitoring of back-diffusion from lowpermeability layers, which is increasingly recognized as an important cause for increased CEs source longevity (Chapman and Parker, 2005;Liu and Ball, 2002;Seyedabbasi et al., 2012). Back-diffusion corresponds to the mobilization of pollutants from silt or clay layers to transmissive zones (Chapman and Parker, 2005;Johnson and Pankow, for contaminants due to their diffusive and sorptive properties (Liu and Ball, 2002). Even a thin contaminated clay layer might maintain CEs concentrations in the aquifer above maximum concentration limits (MCLs) for decades after source removal (Parker et al., 2008;Yang et al., 2015). Usual remediation technologies are not applicable to lowpermeability layers notably due to their inherent low transmissivity (Khan et al., 2004). Fortunately, recent studies show that CEs biotic reductive dechlorination can take place in low-permeability layers despite the small pore size (Chambon et al., 2010;Damgaard et al., 2013aDamgaard et al., , 2013bManoli et al., 2012;Takeuchi et al., 2011;Verce et al., 2015). At some sites low-permeability layers present better conditions for reductive degradation (RD) than the transmissive zone (Parker et al., 2008;Petitta et al., 2013;Wanner et al., 2016). Consequently, natural attenuation in the low-permeability layers may naturally mitigate 'back-diffusion' and thereby limit long remediation time frames (Carey et al., 2015;Damgaard et al., 2013a;Sale et al., 2008;Schaefer et al., 2013;Wanner et al., 2018aWanner et al., , 2018bWest and Kueper, 2010). MNA would benefit from tools able to detect and quantify the RD potential of the low-permeability layers. However, the study of RD in low-permeability layers requires coring them, which is usually prohibitive in regular site assessment (Chapman and Parker, 2005;Filippini et al., 2016;Parker et al., 2004). Compound specific isotope analysis (CSIA) of CEs, which is increasingly employed for MNA at field sites, could potentially help degradation assessment in the low-permeability layers as monitored from the transmissive zone as shown in a recent study simulating CSIA data at the Borden site (Wanner et al., 2018b).
As a result of kinetic effects, CEs degradation leads to enrichment in heavy carbon isotopes ( 13 C) of the degrading compound (Hunkeler et al., 2011b(Hunkeler et al., , 1999. The degree of carbon isotope enrichment is influenced by both the advancement of the reaction and the kinetic isotope fractionation factor, characteristic of the degradation reaction, and therefore informs on the extent of degradation (Abe and Hunkeler, 2006;Hunkeler et al., 1999;Van Breukelen, 2007). Reaction products are initially depleted in 13 C, while the weighted sum of the carbon isotope ratios of the CEs (carbon isotope mass balance, CIMB) is theoretically constant during reductive dechlorination in batch experiments (Aeppli et al., 2010;Bloom et al., 2000;Hunkeler et al., 1999). The combination of these CSIA tools proved itself powerful for detection and quantification of CEs degradation in the subsurface (Höhener, 2016;Höhener et al., 2015;Hunkeler et al., 2011a;Van Breukelen et al., 2005;Wiegert et al., 2013Wiegert et al., , 2012. One potential complication, however, is that several non-destructive processes might also lead to isotope fractionation as shown in small scale lab studies and aquifer scale model simulations. Those processes include highly relevant processes in low-permeability layers, i.e., sorption; with preferential sorption for isotopologues with light carbon isotopes ( 12 C) (Höhener and Atteia, 2010;Kopinke et al., 2005;Qiu et al., 2013;Van Breukelen and Prommer, 2008;Wanner et al., 2017) and diffusion, with slightly smaller diffusion coefficients of the isotopologues containing 13 C, as diffusion is influenced among others by molecular weight (Jin et al., 2014;LaBolle et al., 2008;Rolle et al., 2010;Wanner and Hunkeler, 2015;Xu et al., 2017). Associated isotope fractionation effects were detected in low-permeability layers. Diffusion reportedly leads to depletion in 13 C with depth during short diffusion periods in clay (Wanner and Hunkeler, 2015). Sorption enriches chlorinated hydrocarbons in 13 C with depth (Wanner et al., 2017). While the effects are small compared to RD, sorption-induced isotope fractionation for instance could lead to significant enrichment in the aquitard (Wanner et al., 2017).
Both numerical and analytical modelling were proven useful to the study of volatile compounds or tracer concentrations in the subsurface during back-diffusion processes (Adamson et al., 2016(Adamson et al., , 2015Liu and Ball, 2002;Parker et al., 2008;Yang et al., 2015). Modelling studies were shown essential for investigating aquifer parameters influencing TCE plume tailing (Carey et al., 2015;Maghrebi et al., 2015Maghrebi et al., , 2014Yang et al., 2017).
The isotope ratios of volatile organic compounds in the transmissive zone in the vicinity of a non-reactive aquitard were studied through such RTMs. Diffusion towards low-permeability layers was shown to have a potential for carbon isotope enrichment within the aquifer (LaBolle et al., 2008), the latter being likely attenuated and negligible compared to potential degradation-induced enrichment that could occur in the transmissive zone, also after back-diffusion (Xu et al., 2017). However, processes leading to stronger isotope fractionation in the aquitard such as degradation could potentially lead to detectable enrichment at the aquifer level (Wanner et al., 2018b;Xu et al., 2017). Moreover, diffusion of daughter products from the low-permeability layers to the aquifer was shown to contribute to plume tailing (Rasa et al., 2011;Wanner et al., 2018b). It is relevant to investigate whether any CSIA data of TCE or of its degradation products sampled in the aquifer could be used for degradation assessment despite the low concentration levels. Recently, Wanner et al. (2018b) demonstrated with a modelling study that TCE to cDCE degradation in shallow bioactive zones of aquitards could cause unique CSIA patterns in above lying aquifers.
This study extends to the aforementioned work and sets out to illustrate to which extent spatial and temporal carbon CSIA patterns could be informative of CEs reductive degradation in aquitards despite the impacts of non-destructive processes on both the concentration and isotope ratios. To this end, we apply virtual experiments, i.e., run reactive transport models fed with sets of representative model parameters for common field conditions. For example, the hydraulic parameters chosen for the model were based on previous hydrological models predicting plume tailings for different types of low-permeability layers. Following up to the work of Wanner et al. (2018b), we simulated the full sequence of reductive dechlorination; assessed the effect of thinner aquifers and deeper bioactive zones; and evaluated the performance of the Rayleigh equation and CIMB as common CSIA data interpretation tools. The main objectives of this modelling study are (1) development of a model that simulates carbon CSIA and concentration patterns in an aquifer-aquitard-aquifer system before and after source removal, and (2) evaluating the potential of CSIA data to detect reductive dechlorination and the extent of TCE degradation as occurring in the aquitard while setting a sensitivity analysis on model parameters. With the sensitivity analysis we demonstrate the effect of critical model parameters like aquitard thickness and degree of sorption on simulation results. The performance of the Rayleigh-equation to determine TCE degradation and the use of the carbon isotope mass balance under these specific conditions is also assessed. Finally, we provide recommendations to practitioners on how CSIA data from aquifers can be used and interpreted for the detection and quantification of reductive dechlorination in low-permeability layers.   .A.A. Thouement, et al. Journal of Contaminant Hydrology 226 (2019) 103520 scenario consisting of a 100 m long and 10 m high cross-sectional domain. A 1 m thick clay aquitard separates an 8 m thick upper sand aquifer from a lower aquifer, of which only the upper 1 m is modelled. The model simulates two temporal phases: a loading phase, during which dissolution of a TCE DNAPL source present at the bottom of the upper aquifer leads to TCE advection through the aquifer following groundwater flow and TCE diffusion into the clay aquitard; followed by a release phase, which focusses on back-diffusion of CEs from the clay aquitard to the aquifer after the DNAPL is removed after 10 years (Adamson et al., 2016). A natural hydraulic gradient of 0.0042 m·m −1 was imposed through the boundary conditions, leading together with the selected permeability and porosity (see later) to a groundwater flow velocity of 0.12 m·d −1 in the sand aquifers. For simplicity, sorption of CEs is assumed to occur only in the aquitard. Likewise, sequential reductive dechlorination (from TCE to DCE, VC, and finally ethene) is only assumed to occur with a uniform rate constant across the aquitard (discussed further in Section 2.4). Localised reductive dechlorination (RD) in the aquitard in an otherwise unreactive aquifer is likely as the aquifer, unlike the aquitard, can be too oxic to host reductive processes, as it was previously observed at the Borden site (Wanner et al., 2016). Aquitards usually contain higher contents of soil organic matter, that can potentially act as electron donor in reductive dechlorination (Wanner et al., 2016). The model is described in detail in the following.

Modelling approach
The multicomponent reactive transport code PHT3D (Prommer et al., 2001) was applied using the FloPy Python package (Bakker et al., 2016) as interface. PHT3D allows combining the geochemical model PHREEQC-2 (Parkhurst and Appelo, 1999) with the 3-D multicomponent transport model MT3DMS (Zheng et al., 2012). The implicit TVD method was used. Similarly as in prior publications Van Breukelen et al., 2017, PHREEQC-2 was employed for modelling carbon isotope fractionation during sequential degradation of chlorinated ethenes. MT3DMS has been applied before to simulate back-diffusion and has been validated with analytical solutions and experimental data (Chapman et al., 2012). Each carbon isotopologue was simulated with different values of effective diffusion coefficients (D e ), solid-phase organic carbon-water partitioning (sorption) coefficients (K oc ), and first-order degradation rate constants (k). The selection procedure of the aforementioned model parameters is elaborated below.

General model settings
The model grid for the base case consists of 200 columns, 200 rows, and 1 layer. Grid discretization is of 0.5 m in the horizontal direction (columns) and of 0.05 m in the vertical direction (rows) (Fig. 1). Results obtained with grids of higher and lower resolutions are presented for comparison in Fig. S2 in the Appendix. For simplicity, the porosity was set to 0.35 for both the aquitard and aquifers. This value is typical for unconsolidated sands albeit at the lower end of clay layers. The dry bulk density was set to 1.855 kg·L −1 throughout the system. The aquifer and aquitard were given a hydraulic conductivity of 3650 m·y −1 and 1 m·y −1 , respectively. The low hydraulic conductivity of the clay aquitard leads to negligible flow. Longitudinal and vertical transverse dispersivities were set to 1.0 m and 1.5 mm, respectively. These values are in the range of previously reported dispersivities values of plumes in sandy aquifers (Garabedian et al., 1991;Klenk and Grathwohl, 2002). The selected grid spacing (above) and time step (1 day) insure satisfying grid Peclet and Courant criteria.
Aqueous diffusion coefficient (D w ) values were calculated as 5.9·10 −10 m 2 ·s −1 , 6.9·10 −10 m 2 ·s −1 , 9.7·10 −10 m 2 ·s −1 , and 1.3·10 −9 m 2 ·s −1 for TCE, DCE, VC, and ethene, respectively, as elaborated in Section 2.5. The fraction of organic content f oc of the clay aquitard was set to 0.0033 (0.33%), a representative value for fluvial sediments (Griffioen et al., 2016). Together with the selected values of other model parameter a retardation factor of 3.5 for TCE was obtained similar to prior back-diffusion studies (Carey et al., 2015;Parker et al., 2008). The calculation of the retardation factors is elaborated in Section 2.6. During the simulations linear sorption is assumed. The high CEs concentrations following the dissolution of the DNAPL could saturate the sorption capacity, in which case Freundlich or Langmuir sorption isoterms would provide more realistic results. The implications of the model simplifications are described further in the discussion.

Pollutant source and degradation
A TCE DNAPL source was simulated as a constant concentration boundary condition at the lower 2 m of the upper aquifer for 10 years (Fig. 1). Such source positioning was previously used in the study of back-diffusion of both NAPL and DNAPL (Xu et al., 2017). The TCE source concentration (~0.6 g·L −1 ) was about half the TCE solubility (1.28 g·L −1 ) and~6·10 5 times larger than the typical carbon CSIA quantification limit (1 μg·L −1 ). Concentrations and isotope ratios are shown only for concentrations above this quantification limit.
The source was given a constant isotope composition δ 13 C-TCE Source of −30‰. Setting a constant isotope ratio for a dissolving DNAPL is acceptable as the model does not simulate degradation in the sand aquifer and therefore no weathering of the DNAPL (Hunkeler et al., 2004;Hwang et al., 2013). The three carbon isotopologues of each CE and ethene are simulated as done before for chlorine isotopologues ). The TCE carbon isotopologue concentrations in the source were calculated from the TCE concentration and isotope ratio of the source, as described in Appendix Section S3. After a simulation time of 10 years, the loading phase was stopped by setting the concentration of each CE at the constant concentration boundary condition to 0 g·L −1 .
Biological reductive dechlorination in the aquitard degrades sequentially TCE to DCE, VC, and finally ethene. The model assumes that only the major cis-DCE isomer is formed. First-order degradation rate constants for biotic reductive dechlorination are selected in the line with the values employed in (Van Breukelen et al., 2005), i.e., 1 y −1 , 0.2 y −1 , and 0.2 y −1 for TCE, DCE, and VC, respectively (Table S4 in Appendix). Although the high source concentration applied calls for Monod-kinetics (Bekins et al., 1997), first-order kinetics often used at field sites was chosen for simplicity, and also further additional complexity such as microbial growth, inhibition, and donor limitation (Chambon et al., 2013) was ignored. Degradation rate constants were assumed to remain constant with depth in the base case model as opposed to two recent studies where bioactivity was assumed to be (i) retrained a centimeters thin bioactive zone at the top of a thick aquitard as modelled by exponentially declining rate constants with depth  (Wanner et al., 2018b), or (ii) declining gradually with depth as modelled by nearly linearly decreasing rate constants till about 1 m depth (Wanner et al., 2018a). In those previous studies, likely the microorganisms distribution, lead to decline of the degradation rate constants with depth. In the one study, the decreasing availability of nutrients and electron donors, diffusing from the aquifer to the aquitard, likely influenced the degradation rate constants (Wanner et al., 2018b). In the other study, the limited presence of microorganisms responsible for degradation of the chlorohydrocarbons of concern deeper than 0.35 cm together with the near linearly decreasing degradation rate constants below this depth could indicate that transport of microorganisms remains a potential limitation for biotic degradation across thick aquitards (Wanner et al., 2018a). For thin aquitards (< 1 m) with rich organic carbon content, the contact with aquifers both above and below is likely to allow a more homogeneous bacterial distribution than for the aquitards described in those prior studies. Moreover, several bacterial strains including Dehaloccocoides spp., a strain able of complete TCE to ethene degradation, were present across the studied depth of~1.1 m (Wanner et al., 2018a). Unfortunately, at this stage studies on degradation and bacterial activity in aquitards are still limited. The choice of low uniform degradation rate constants in the aquitard is therefore a simplification which could be refined later once bacterial degradation in aquitard is better understood.
Each isotopologue degrades with a slightly different rate constant, leading to isotope fractionation of the CEs. For carbon, this fractionation is usually modelled as a bulk isotope fractionation, i.e., simplifying the model to the evolution of the compound's fraction of the light and the heavy isotope . For this isotopologue model, we express the respective degradation rates of each isotopologue employing the bulk kinetic carbon isotope enrichment factor ε Degradation : where LL C, LH C, and HH C are the isotopologue concentrations (moles) of the carbon isotopologues containing two light (LL) carbon isotopes, one light (L) and one heavy (H) isotope (LH), and two heavy isotopes (HH), respectively; and k (y −1 ) the first-order degradation rate constant. Carbon isotope fractionation of chlorinated ethenes during biotic reductive dechlorination has been intensively studied. For this model, average values of bulk carbon isotope enrichment factors ε Degradation obtained with microcosm experiments were used for the base case scenario, i.e., −12‰, −21‰, and − 23‰ (Table S3), during biotic reductive dechlorination of TCE, DCE, and VC, respectively. The variability of enrichment factors reflects the diversity of degradative microorganisms. For instance, ε Degradation during TCE reductive dechlorination was found to vary from −3.3‰ to 18.9‰ (Cichocka et al., 2008). We restrained from assessing this variability on the model outcome, testing only a single lower value for each ε Degradation in an additional simulation in Section 3.5.

Diffusion and diffusion-induced isotope fractionation
Aqueous diffusion coefficients of each CE isotopologue were acquired from the aqueous diffusion coefficients of the CEs. Each CE's aqueous diffusion coefficient D w was first calculated through (Worch, 1993): where T is the water temperature (°K), set for the model to 283.15 K (10°C), η T the water's viscosity at T (Pa·s), and M CE is the CE's weighted molar mass using isotopologues abundances at 0‰ (Table S5 in Appendix). Each isotopologue was then assigned a specific diffusion coefficient to model diffusion-induced isotope fractionation. Recent work showed that for a CE with N isotopologues, isotopologue-specific aqueous diffusion coefficients (D w,n ) of the isotopologue n are defined by (Jin et al., 2014): where A i and m i are the relative abundance and the mass of the isotopologue i, respectively, and β is an empirical coefficient determined experimentally for each compound. The abundances at 0‰ and masses of the isotopologues, given in Table S5 in Appendix, were then employed to calculate the diffusion coefficients of each isotopologue. The diffusion coefficient of the compound will be slightly modified by the change in its isotopologue composition. Measurements obtained with a modified Stroke's diffusion cell (Wanner and Hunkeler, 2015) yielded β values of 0.029 ± 0.008 and 0.024 ± 0.002 for carbon and chlorine isotope fractionation of TCE, respectively. Through the modelling of diffusion experiments in gel diffusion tubes a slightly higher β value for chlorine isotope enrichment of TCE was obtained (0.043 ± 0.008; while of DCE = 0.088 ± 0.015 (Jin et al., 2014)). The difference between TCE and DCE might be explained by the role of the hydration shell in the extent of diffusioninduced isotope fractionation (Jin et al., 2014).
In the present model, we used the β values reported by Wanner et al. (Wanner and Hunkeler, 2015). Considering the similarity of β values for carbon and chlorine isotope fractionation of TCE, we assumed that the same applies to DCE. Therefore, we used β = 0.088 for the carbon isotope effects of DCE. No experimental β values for VC and ethene are available. For those compounds, we used the value of 0.088, given that DCE in the closest available analogue for these compounds (Jin et al., 2014). Resulting diffusion-induced carbon isotope enrichment factors ε Diffusion are −0.2‰, −0.9‰, −1.4‰, and − 3.1‰ for TCE, DCE, VC, and ethene, respectively (Table S4). Only 1-D w,LH /D w,LL is given here for simplicity, where D w,LH and D w,LL are the aqueous diffusion coefficients for the isotopologue containing one heavy and one light and two light isotopes, respectively (Table S5 in Appendix lists the isotopologuespecific aqueous diffusion coefficients).
In porous media, molecules diffuse around the sediment grains of the rock matrix. Therefore, the effective diffusion coefficient D e is smaller than D w in order to account for the added travel distance. D e was taken as D e = D w ⋅ τ app (7), where D w is defined for each isotopologue using Eq. 4, and τ app is the apparent diffusion tortuosity coefficient (also referred to as the tortuosity factor) (Carey et al., 2016). The coefficient τ app is often written as a function of the accessible porosity for simplicity (Boving and Grathwohl, 2001). Recent work suggests that the coefficient τ app is correlated to the permeability, not to the porosity of a media, with an average τ app of 0.33 for low-permeability media (Carey et al., 2016). Fitting of TCE diffusion profiles in clay suggested much lower τ app values (0.13, (Wanner and Hunkeler, 2015); 0.03, (Wanner et al., 2016)). In the present study, τ app is taken equal to the porosity across the model (0.35) in the base case scenario.

Sorption and sorption-induced isotope fractionation
The retardation factor, R, represents the degree of retardation of a compound with respect to the local groundwater flow velocity and is calculated as follows: Thouement, et al. Journal of Contaminant Hydrology 226 (2019) 103520 where ρ b is the bulk density across the aquifer-aquitard-aquifer system (kg·L −1 ), K oc is the solid-phase organic carbon-water partition coefficient (L·kg −1 ), and n the porosity. K oc values are usually derived from its correlation with the octanol-water partitioning coefficient, K ow (Lu et al., 2011;Nguyen et al., 2005 (Table S4) (Valsaraj et al., 1999).
The present simulation sets f oc to 0.3% in the aquitard so TCE retardation is of 3.5 as in prior modelling studies of aquifer-aquitard systems (see Table S1 in Appendix). The selected f oc value results in retardation factors, R, of 2.3, 1.8, and 1.5 for DCE, VC, and ethene, respectively, after applying Eq. 6 and Eq. 7.
Sorption-induced isotope fractionation results from the slightly stronger sorption to organic matter of molecules with light isotopes than those with heavy isotopes (Kopinke et al., 2005). The fractionation factor for linear equilibrium sorption, α Sorption, has been defined as the ratio of the K oc of the heavy ( H K oc ) and the one of the light isotope ( L K oc ) (Höhener and Yu, 2012;Kopinke et al., 2005;Van Breukelen and Prommer, 2008;Wanner et al., 2017). The resulting sorption-induced isotope fractionation is included in the model by attributing different K oc values to each isotopologue (See Table S5 in Appendix). The calculation of isotopologue-specific K oc values was inspired by (Höhener and Yu, 2012). We extended their approach for comparing data of differently deuterated compounds to carbon isotopologues, for which no specific relation was yet developed. The resulting relation is: where LH α Sorption and HH α Sorption are the sorption-induced fractionation factors of the isotopologue containing one light (L) and one heavy (H), and two heavy (HH) carbon isotopes, respectively. The K oc values of the two isotopologues containing at least one heavy carbon atom can be calculated from: where LL K oc , LH K oc, and HH K oc are the K oc values of the isotopologue containing two light carbon isotopes, one light and one heavy, and two heavy, respectively. Isotopologue-specific K oc values were estimated using the fractionation factor α Sorption taken from the literature (see below) and the reported K oc values. LL K oc ≈ K oc was assumed, which is an acceptable approximation considering the large variation of K oc . An example of the calculation of the isotopologue-specific K oc is presented in the Appendix Section S5, and the values are reported in Table S5. Simulations of sorption-induced isotope fractionation during advection with an analytical solution (Van Breukelen and Prommer, 2008) yielded identical results for the isotope and this isotopologue approach (see Fig.  S1 in Appendix).

Post-modelling calculations 2.7.1. Isotope ratio
The heavy to light carbon isotope ratios of a sample R H/L are calculated from the concentrations of the isotopologues as follow: The isotope ratios are reported in the usual delta notation δ 13 C. The δ 13 C is defined as δ 13 C = R R/L /R VPDB − 1 where R VPDB is the carbon isotope ratio of an international standard. Isotope ratios are calculated only when the CEs concentration is above isotope quantification limit, taken as 1 μg/L for all CEs in this study.

Carbon isotope mass balance
The carbon isotope mass balance (CIMB) is the concentrationweighted average isotope ratio of the CEs and ethene combined: where δ 13 C i and C i are carbon isotope ratios and concentrations, respectively, of the CEs and ethene. CIMB is not calculated when any of the CEs concentrations are below the minimum quantification limit for δ 13 C (taken as 1 μg·L −1 for all CEs for simplicity). The concept of CIMB has shown useful in other studies (Aeppli et al., 2010;Amaral et al., 2011;Hunkeler et al., 1999) as its enrichment points to, for example, oxidative transformation, whereas reductive dechlorination yields constant CIMB. The CIMB can also be calculated excluding ethene, which δ 13 C is not always available, in which case constant CIMB would be a signal of a lack of VC transformation to ethene, for instance.

Performance of the Rayleigh equation to estimate TCE degradation
The Rayleigh equation is commonly applied to estimate the extent of biodegradation of a compound based on its degree of isotope enrichment and the applicable isotope enrichment factor (Hunkeler et al., 2002 where D Rayleigh is the Rayleigh-based estimate of TCE degradation, and ε Degradation is the bulk carbon enrichment factor during TCE degradation. A D Rayleigh of 0.1 signifies that 10% of the initial TCE was degraded. The true extent of biodegradation as occurred in the model, D Sample , is calculated from the simulated TCE concentration and the concentration of a non-degrading tracer with otherwise similar characteristics to TCE. D Sample is calculated for a cell i of the non-sorbing aquifer as follows: where [TCE] i and [Tracer] i are the TCE and tracer concentrations, respectively, in cell i. The Rayleigh equation systematically underestimates the true extent of degradation (Abe and Hunkeler, 2006;Höyng et al., 2015;Van Breukelen and Prommer, 2008) but could under some conditions also overestimate if δ 13 C enrichment resulting from non-degradative processes such as sorption at plume fronts or diffusion in plume cores would be attributed to degradation (Höhener and Yu, 2012;Kopinke et al., 2005;Van Breukelen and Prommer, 2008;Van Breukelen and Rolle, 2012). In order to decrease below the MCL of 5 μg/L, the initial TCE concentration of 0.6 g/L must be reduced by 99.999%. At such a high number of decimals, the ratio of the extents of degradation D Rayleigh /D Sample would be a poor quantifier of the underestimation at overall high extents of degradation. Using instead the ratio of the rate constants as a measure of the underestimation as done in previous studies (Lutz et al., 2013;Van Breukelen and Prommer, 2008;Van Breukelen and Rolle, 2012) circumvents this issue as it equals the ratio of the logarithms of the residual fractions, given that the rate constant k of a first-order reaction satisfies −k ⋅ t = ln (1 − D) where D is the extent of degradation at the time t. In the present study, the deviation between the true extent of degradation and the Rayleigh-based estimate is addressed through the following measure θ, which quantifies the underestimation of the true first-order degradation rate constant k Sample by the Rayleigh-based estimate of the first-order degradation rate constant k Rayleigh . The values of the rates themselves are not the subject of this paper.
A θ of 1 signifies that both rate constants are equal, i.e., no underestimation of D Sample by D Rayleigh ; a θ of 0.8 indicates that the estimate k Rayleigh is 1.25 times smaller than the true k Sample , whereas a θ of 0.1 means that k Rayleigh is 10 times smaller than the true k Sample , i.e., D Rayleigh strongly underestimates D Sample . θ is calculated only for significant δ 13 C-TCE enrichment, i.e., the isotopic shift Δ 13 C between δ 13 C-TCE and δ 13 C-TCE Source , is beyond 1‰. A value below a threshold of 2‰ cannot be taken with certainty as an indication of degradation considering the analytical uncertainty of carbon CSIA and site heterogeneity (Hunkeler et al., 2008); however, in practice, enrichment beyond 1‰ is considered to be an indication of degradation.
Two additional concepts are defined for the sake of this study to estimate the true degradation effect at a larger scale than the scale of a sample. The total mass decrease of TCE versus the one of the tracer as occurred in the model is calculated at a certain distance downgradient for (1) the entire depth of the aquitard, D Aquitard , and (2) the entire depth of the aquifer-aquitard-aquifer system, D System . D Aquitard and D System are therefore the ratios of TCE and the tracer's mass across a column of the model as in Eq. 16: where N is the number of rows in the aquitard or in the aquifer column, respectively; [TCE] i and [Tracer] i the concentration for this column at the row i (mg·L −1 ); f oc is the local fraction of organic carbon (only different between the aquifer and the aquitard); η is the porosity (dimentionless); K oc,TCE the solid-phase organic carbon-water partition coefficients of TCE (identical for the tracer) (L·kg −1 ); and ρ b is the bulk density, V i the the cells volume (m 3 ), the length and width of the model cells are always the same but the height of the cells varies accross depth. D Aquifer and D Aquitard are 0% when there is no degradation and 100% when degradation is maximal.

Complete model including CEs degradation: loading phase
The base case scenario is employed to simulate CEs concentrations and isotope ratios following flow, diffusion, sorption, and degradation including associated carbon isotope effects. This section focusses on the loading phase, until the plume is removed. During the loading phase, the dissolving DNAPL forms a plume with high TCE concentration in the aquifer. The concentration gradient between the aquifer and the aquitard drives TCE diffusion into the initially pristine aquitard. In the aquitard, transport is governed by diffusion only; advective transport does not occur. At the end of the loading phase in year 10, daughter products show the highest concentration a few centimeters (2 cm to 15 cm) below the interface (Fig. 2a), resulting from their sequential production in the aquitard (close up of the aquitard available in the Appendix Fig. S4). The centers of the CEs masses are spatially separated, with ethene present the deepest, followed by VC, DCE, and TCE ( Fig. 2a and b). CEs sequential reductive dechlorination, sorption, and diffusion create these patterns in concert. The more dechlorinated the compound, the later it is formed in the degradation sequence, the less it sorbs, and the faster it diffuses. For instance, ethene is the end product of reductive dechlorination, it sorbs least, and diffuses fastest compared to the other compounds. Whereas the peak concentration of TCE in the aquitard decreases downgradient, those of the daughter products remain remarkably similar (Fig. 2a, plain vs. dashed lines). Daughter products which back-diffused to upgradient parts of the aquifer are transporteddowngradient where they again diffuse into the aquitard. Consequently, daughter product concentration peaks remain high right below the aquifer-aquitard interface throughout the aquitard, also away from the source (dashed line Fig. 2a).
Simulations allow comparing the base case scenario DySy, which includes both the diffusion-and sorption induced isotope fractionation effects, to a model disregarding both those effects (DnSn), as well as two intermediary scenarios (DySn and DnSy, see description of these simulation codes in caption of Fig. 2). Details of additional simulations in the absence of degradation are presented in the Appendix Section S7. For all scenarios including TCE degradation, δ 13 C-TCE enriches with depth in the aquitard with as much as +45‰ at 30 cm depth in year 10 (Fig. 2c). This enrichment is strongly larger than the negligible δ 13 C-

Fig. 2. Depth profiles of concentrations (a), molar fractions (b), and carbon
isotope ratios (c) of CEs and ethene (TCE, red; DCE, blue; VC, cyan; ethene, green; CIMB, black) at the end of the loading phase (year 10), near the aquiferaquitard interface, 2 m downgradient of the source (except indicated otherwise on (a)). (c) The base case simulation DySy is compared to simulations DySn, DnSy, and DnSn: DySy is the base case scenario, with diffusion-induced (Dy) and sorption-induced (Sy) isotope fractionation; DySn includes diffusion-induced isotope fractionation (Dy) but omits sorption-induced isotope fractionation (Sn); DnSy includes sorption-induced isotope fractionation (Sy) but omits diffusion-induced isotope fractionation (Dn). DnSn does not include any isotope fractionation induced by physical processes. The vertical dashed black line (c) represents δ 13 C-TCE Source . The shadded grey zones indicates the aquitard. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) TCE enrichment predicted in the absence of reductive dechlorination (< 0.5‰, Fig. S3c), resulting solely from the combined effects of sorption-and diffusion induced isotope fractionation. Each lesser chlorinated ethene presents a more depleted δ 13 C than its precursor, which is the usual pattern for reductive dechlorination Slater et al., 2000;Van Breukelen et al., 2005). For the base case scenario DySy, δ 13 C of daughter products tend to curve with depth opposed to the straight enrichment of δ 13 C-TCE. After initial enrichments compared to their starting values, δ 13 C level off, and even show a declining trend deep in the aquitard (Fig. 2c). When DCE or VC eventually travels ahead of its precursor, the net δ 13 C enrichment becomes smaller. Eventually, δ 13 C decreases with depth below the δ 13 C peaks for both the models including diffusion-induced isotope fractionation (Fig. 2c: DySy and DySn). The ethene δ 13 C pattern is different as ethene does not degrade; δ 13 C-ethene will not further enrich once it is away from its precursor VC. With the selected ε Diffusion , ethene shows stronger depletion in the base case scenario compared to the simulation omitting diffusion-induced isotope fractionation (Fig. 2c:~− 7‰ for the base case model DySy vs. DnSy at 90 cm depth). Despite a f oc of 0.3%, ε Sorption itself leads to a small enrichment for all CEs and is significant (> 2‰) only in the deepest parts of the aquitard (Fig. 2c: DySy vs. DySn and DnSy vs. DnSn). In the aquifer, the upwards decline of δ 13 C of DCE, VC, and ethene by a few ‰ (Fig. 2c) is caused by the faster diffusion of the light isotopes following the upward concentration gradients.
The CIMB increases just below the aquifer-aquitard interface where enriched TCE and DCE are dominant. The CIMB increase results from the faster downward diffusion of the depleted daughter products (Fig. 2c). Beneath, the CIMB strongly decreases with depth until being equal to δ 13 C-Ethene once ethene is the dominant compound (depletion by about −25‰ at −0.9 m, Fig. 2c). The occurrence of depleted CIMB was described before in simulations of CEs transport and degradation through an aquifer in presence of sorption (Van Breukelen et al., 2005). The different transport velocities of the CEs is the main reason for downwards CIMB decrease, further enhanced by diffusion-induced isotope fractionation (Fig. 2c: DnSy and DnSn). The most depleted CIMB predicted in this simulation (−55‰) is close to the range of the lowest reported isotope signatures of PCE and TCE at a field site (−67.9‰ and − 70‰ respectively (Nijenhuis et al., 2013)). As CEs isotope profiles in aquitards are employed for retracing CEs source history (Adamson et al., 2015), such low CIMB in a sorbing aquitard should not be confused for the prior presence of another strongly depleted TCE spill.
Concentration and δ 13 C patterns of CEs across the aquifer at the end of the loading phase (year 10) are presented in Fig. 3. Daughter products accumulate in the aquifer where their concentrations increase with distance from the source (Fig. 3a). At the end of the loading phase, daughter product concentrations at the bottom of the aquifer are already several orders of magnitude above the MCLs (MCL = 70 μg·L −1 for cis-DCE, 2 μg·L −1 for VC) (Fig. 3a). At year 10, although no degradation occurs in the aquifer, δ 13 C-TCE downgradient enrichment just above the aquifer-aquitard interface is strongly larger than in absence of degradation (isotope enrichment Δ 13 C~5‰ vs.~0.05‰; Fig. 3c and S3 g, respectively). TCE degradation in the aquitard enhances the concentration gradient at the aquifer/aquitard interface. The resulting faster diffusion to the aquitard leads to a stronger isotope enrichment in the aquifer. The same process leads to the increase of δ 13 C-DCE and δ 13 C-VC with distance from the source (up to~5‰ for DCE, Fig. 3c). However, δ 13 C-Ethene decreases with distance as ethene accumulates and as degradation is overall less advanced downgradient from the source (Fig. 3c).
Until the release phase, daughter product concentrations are much lower than TCE concentrations in the aquifer and δ 13 C-TCE is enriched. Consequently, the CIMB is dictated by δ 13 C-TCE at the bottom of the aquifer at the end of the loading phase, and therefore enriches with distance (Fig. 3e).

Complete model including CEs degradation: release phase
After source removal, the main plume is flushed away from the 100 m long stretch of aquifer by upgradient uncontaminated groundwater within about 2.7 years (the groundwater flow velocity is of 44 m·y −1 ). The TCE concentration gradient is reversed and TCE diffuses now back from the aquitard to the aquifer. Prior to source removal, the high TCE concentration levels forced the δ 13 C of CEs to low values. Five years after source removal, at year 15, δ 13 C values of all compounds are considerably higher than during the loading phase and are above their minimum detection limit of 1 μg/L. Interestingly, δ 13 C now strongly decreases downgradient (Fig. 3d). The downgradient decline is as much as −10‰ for DCE and VC, and near −30‰ for TCE (Fig. 3d). The ongoing uniform degradation in the aquitard explains δ 13 C-TCE enrichment over time in the entire system, while the overall downgradient δ 13 C-TCE depletion relates to the delay in flushing of the downgradient parts of the aquifer (in Appendix; and discussed further in Section 3.3).
The temporal evolution of the CEs and ethene are evaluated near the aquitard, 95 m downgradient of the source. The TCE concentration drops below the MCL within 10 years after the source is removed (Fig. 4). Concentrations of daughter products near the aquitard also decline, but slower, and the accumulating product of CEs sequential degradation, ethene, resides in the aquifer nearly as long as the nondegrading tracer (Fig. 4a). VC concentration away from the source is still above its MCL at year 40, i.e., 30 years after source removal, although TCE declined below its MCL long before, at year 18 (Fig. 4a). Simulations of VC or DCE stall, presented in the Appendix, show that their concentrations would also remain nearly as high as that of the tracer in the aquifer (Fig. S7 and S8). Despite degradation occurs in the aquitard and not in the aquifer, the potentially high concentrations and the longevity of the degradation products in the aquifer should be considered in the study of back-diffusion.
Shortly after the cleanup, CIMB significantly increases by~3‰ and remains high until the end of the simulation (Fig. 4b), also for different well screen lengths and positions (Fig. S5 in the Appendix). The CIMB increase points to the loss of depleted compounds from the aquifer during the release phase. The more dechlorinated and more depleted compounds have been transported deeper into the aquitard (see Section 3.1) and do not diffuse back to the aquifer as fast as the more enriched and chlorinated compounds. The difference in CIMB, ΔCIMB, was evaluated for several reaction chains, with DCE or VC as end product instead of ethene. ΔCIMB is smaller, but still > 1.5‰ after source removal when VC is the end product (Fig. S7). ΔCIMB decreases in a few years to near 0‰ when DCE is the end product (Fig. S8). At field sites, CIMB is expected to increase only in presence of degradation processes other than reductive dechlorination. For instance, DCE and VC are more readily oxidizable than TCE and could degrade in the aquifer provided the right conditions are met (Coleman et al., 2002a(Coleman et al., , 2002b. Therefore, such significant CIMB increase could be confused for alternative degradation pathways for the CEs.

Distinct CSIA patterns in the aquifer in presence of aquitard degradation
With the goal to assess spatial CSIA patterns characteristic of either aquitard RD or aquifer RD, the base case simulation on aquitard RD was compared with another simulation where RD only occurred in the aquifer (Fig. 5). Aquifer RD causes monotonously increasing CSIA patterns of all compounds in downgradient direction during both the loading and release phase (Fig. 5a,b) similarly as known for aquifer RD without aquitard interactions (Van Breukelen et al., 2005). Strikingly, downgradient CSIA patterns and their magnitude remain almost unaffected after source removal and remain constant in time long after source removal (results not shown). During the loading phase, the CIMB is not increasing downgradient, but it slightly increases by about 0.5‰ during the release phase. This minor CIMB increase after source  removal can be explained by the lower back diffusion rates of the more dechlorinated and more depleted daughter products which have penetrated deeper into the aquitard as they diffuse faster and sorb less.
These patterns are in strong contrast with those of uniform aquitard RD as discussed shortly before, where carbon isotope ratios only limitedly increase downgradient during the loading phase, whereas most strikingly, after source removal huge isotope enrichment is simulated of all compounds as the aquitard serves as a reservoir of degraded CEs. As in both cases, concentration drop during the release phase, isotope ratios is more applicable to discriminate between both hypotheses. Furthermore, the CSIA patterns decreasing downgradient (Fig. 5b), and downgradient CIMB increase is substantial especially during the release phase with about 4‰ (Fig. 5c,d). A simulation assuming linear decreasing degradation rate constants dropping to zero at the base of the aquitard, showed similar CSIA patterns as the base case simulation (Fig.  S6). Note that in the case bioactivity is limited to the upper centimeters of aquitards further isotope enrichment over time after the initial jump will not occur, while isotope ratios will also not decrease downgradient (Wanner et al., 2018b) because the residence time of back-diffusing TCE in the bioactive zone is invariant. Interestingly, CSIA patterns may thus inform on the depth of bioactive zones in aquitards, or in general the depth of reactivity as aquitards with potentially elevated contents of iron-sulfur minerals could also host chemical reductive dechlorination as CE attenuation process.
In the case RD occurs within the entire system (both in the aquitard and the aquifer(s)), the CSIA patterns become a mixture of the unique features of aquitard and aquifer RD. During the loading phase, downgradient CSIA patterns increase for the higher chlorinated CEs alike aquifer RD but are more stable for lower chlorinated CEs and ethene similar to aquitard RD. Whereas during the release phase aquifer RD shows increasing CSIA patterns, and aquitard RD shows strongly elevated but decreasing CSIA patterns, RD occurring within the entire system shows elevated but relatively constant CSIA patterns downgradient at least under the simulated conditions of spatially homogenous first-order rate constants. Even under those conditions aquitard RD is still noticeable as CSIA patterns of daughter products are much more enriched than in the absence of aquitard RD, notably near the source. The CIMB increase is in between those of purely aquitard and aquifer RD except near the source, where the increase is large which differs strongly from the aquifer RD pattern. Concluding, aquitard RD can still be clearly detected by aquifer CSIA in case also aquifer RD occurs because of (i) the significant CIMB increase, and especially (ii) the strong jump of CSIA ratios after source removal also near the source, and (iii) the absence of clear increasing and possibly decreasing CSIA trends downgradient.

Performance of the Rayleigh equation in assessing the extent of degradation
At contaminated sites with a predominant role of aquitard RD, CSIA of samples from conventional wells situated in the aquifer, but close to the aquifer-aquitard interface could detect and potentially be used to quantify the extent of degradation. The base case scenario shows that aquitard RD leaves an imprint on δ 13 C-TCE in the aquifer. Additionally, daughter products are detected in the aquifer throughout the simulation. This section evaluates whether a CSIA calculation tool often used for biodegradation assessment at field sites, the Rayleigh equation, provides trustworthy information on the occurrence and degree of TCE degradation in the aquitard when based on CSIA samples from the aquifer.

Performance of the Rayleigh equation for local degradation assessment
The direct application of the Rayleigh equation to δ 13 C-TCE data obtained from a polluted site is in theory sufficient to obtain a conservative estimate of TCE degradation as occurred in the samples, D Rayleigh . The underestimation by the Rayleigh equation of the true extent of degradation for the simulated conditions can be assessed by comparing D Rayleigh to the true extent of degradation D Sample as both are obtained from the synthetic dataset (Lutz et al., 2013). Prior modelling studies pointed out that D Rayleigh systematically underestimates the true extent of degradation under various conditions (Abe and Hunkeler, 2006;Höyng et al., 2015;Lutz et al., 2017;Lutz and van Breukelen, 2014;Van Breukelen and Prommer, 2008;Van Breukelen and Rolle, 2012). The underestimation was then the result of the attenuation of degradation-induced δ 13 C enrichment through mixing and dispersion of pollution plumes in aquifers. At high extents of degradation, for instance above 90%, the ratio of the extents of degradation D Rayleigh and D Sample would be a poor indicator of the underestimation. The ratio of the Rayleigh-equation estimate to the true first-order degradation rate constant, θ, better qualifies the underestimation also at high extent of degradation (see Section 2.7). Fig. 6 illustrates the variation of θ in the upper aquifer at various elevations and distances from the source, and screen lengths. Sampling from long well screens or near the source lowers the chance to detect significantly enriched δ 13 C-TCE, either because of the high TCE concentration during the loading phase, or due to strong dilution during the release phase. Because θ is only calculated for Δ 13 C > 1‰, the longest time series of θ can be computed for short screens, away from the source, and near the aquitard (Fig. 6). When calculated, θ varies between 0.5 and 0.9, i.e., the Rayleigh equation underestimates the real first-order degradation rate constant with a factor 2 to 1.1, respectively.
Away from the source, θ is initially low and progressively decreases from~0.65 to~0.5 during the loading phase (Fig. 6a, cyan lines). In this initial stage, the TCE concentration and its isotope ratio have plateaued, while the tracer concentration keeps rising (Fig. 4a,b), indicating that the true extent of degradation, D Sample , increases. Consequently, the gap increases between the CSIA-based degradation estimate and the true extent of degradation, causing θ to drop.
Within 10 years after source removal, θ increases to high values, meaning that the Rayleigh-equation only slightly underestimates degradation as occurred in samples and a close estimate is obtained. TCE back-diffusing from the aquitard where it has been partially degraded is now the main source of TCE in the aquifer, and the δ 13 C-TCE enrichment corresponds better to the actual degradation at the level of the aquifer sample.

Can mass destruction in the aquitard be estimated with CSIA data?
It remains the question whether locally obtained D Rayleigh values from the aquifer can be used in assessing the overall TCE mass decrease as occurred in the aquitard which would be highly valuable for estimating remediation time frames. For this analysis, we employ two additional concepts: the total mass decrease of TCE versus the one of the tracer as occurred in the model at a certain distance downgradient of the source for (1) the entire depth of the aquitard, D Aquitard , and (2) the entire depth of the aquifer-aquitard-aquifer system, D System . These mass decreases are calculated from the ratio of TCE mass versus the mass of the non-degradative tracer with similar properties as TCE following Eq. 16. D System and D Aquitard are compared at both 2 m and 95 m from the source to D Rayleigh and to the true extent of degradation D Sample , both obtained at the bottom of the aquifer (Fig. 7).
During the loading phase and near the source, both the lack of detectable δ 13 C-TCE enrichment and the high TCE concentration levels prevent D Rayleigh and D Sample to depart from 0% (Fig. 7, black and cyan lines, respectively). However, the large D System and D Aquitard values of 15% and~70%, respectively ( Fig. 7a red and blue dashed lines, respectively) point to TCE mass decrease in both the aquitard and the aquifer.
During the loading phase and away from the source (95 m), D Rayleigh increases soon after the arrival of the plume but underestimates D Sample considerably at year 10 ( Fig. 7b). D Rayleigh even more strongly underestimates D Aquitard . At year 10, D Rayleigh is in the range of D System but those similar values are coincidental as D System depends among others on the thickness of the plume. Although δ 13 C-TCE enrichment is detectable in the aquifer already in the loading phase, D Rayleigh values may give pessimistic views on the actual high degree of mass destruction in the aquitard, D Aquitard .
Once the source is removed after 10 years, D Aquitard quickly increases both near and away from the source as a result of the continuous TCE degradation in the aquitard without further admixing of non-degraded TCE from the source. The poorly degraded TCE mass in the aquifer is flushed by clean water, therefore the mass remaining in the system is mostly contained in the reactive aquitard and D System converges towards D Aquitard . Both D System and D Aquitard are large due to the ongoing degradation in the aquitard. Both D Sample and D Rayleigh in the aquifer are now influenced only by the back-diffusing TCE and rapidly increase similarly to D Aquitard . Once the source is removed, degradation assessment through CSIA is improved due to the enriched isotope signature of the back-diffusing TCE, which interestingly provides an insight to the state of degradation in the aquitard. Quickly after source removal, TCE mass effectively decreased in the aquifer system and D Rayleigh provides a conservative estimate of the state of degradation.

Sensitivity analysis of degradation assessment
The results of a sensitity analysis of key model parameters on CIMB and θ values are shown and discussed in detail in the Appendix S13. The main outcomes are reported here. With respect to CIMB, the most influential parameters appear to be the aquitard f oc (Fig. S10c) and aquitard thickness (Fig. S10e). CIMB is sensitive to the spatial separation of the depleted daughter products from their more enriched precursors as discussed before. Especially during the release phase, a higher aquitard f oc therefore yields a larger CIMB enrichment. Noteworthy, also without sorption, the CIMB is predicted to enrich in later stages of the release phase. Probably this is caused by the loss of depleted daughter products which have diffused through the aquitard and reached the lower aquifer. With respect to aquitard thickness, a thinner aquitard results probably in an increased mass loss of depleted daughter products to the underlying aquifer while the more enriched precursors sorb stronger and are still present in thinner aquitards during the back-diffusion phase causing a strong CIMB enrichment. Preferential mass loss of depleted daughter products to the underlying aquifer is thus another phenomenon contributing to high CIMB enrichments that could be observed in overlying aquifers. The CIMB of samples taken near aquitards is in most cases not a reliable indicator of the occurrence of alternative degradation pathways besides reductive dechlorination like oxidation. The sensitivity analysis further shows that Rayleigh-based estimates of degradation are always conservative. Moreover, we show that both δ 13 C-TCE enrichment prior to source remediation and a sudden change in δ 13 C values once the source is remediated are to be expected under many different field conditions. The Rayleigh equation systematically underestimates the extent of TCE degradation in aquifer samples especially during the loading phase and for conditions leading to long remediation time frames (low groundwater flow velocity, strong sorption in the aquitard, thicker aquitards).

Model limitations
The simulations required simplifications which influence the results.
The key assumptions concern mostly reaction, i.e., the simulation of linear sorption, and of first-order and uniform degradation kinetics in the aquitard. Diffusion-related assumptions are also discussed. Their implications on the main results are described in this section.
As presented in the method section, the model simulates linear sorption of chlorinated ethenes in aquitards. Sorption is site-specific, and while CEs sorption was near linear in natural clayey tills (Lu et al., 2011), TCE sorption was best represented with Freundlich isotherms in several clay-rich aquitards in a prior study (Allen-King et al., 1996). Freundlich or Langmuir sorption behavior implies that high CE concentrations compared to the sorption site concentration would overwhelm the sorption capacity of the aquitard. The δ 13 C patterns of the CEs would therefore tend towards the patterns observed in the absence of sorption as in Fig. 9c and 10c. A non-linear sorption behavior would mainly limit the CIMB enrichment induced by the different sorption potentials of the CEs.
First-order kinetics were assumed for simplifying the transport model. The high source concentration calls for Monod kinetics, in which case degradation at the high CE concentration levels would comparatively be slower than in our simulation. In our simulations, slower degradation rates lead to better Rayleigh-based degradation assessment and to lesser increase of the CIMB during the release phase; Fig. 7. CSIA-based estimate of degradation (D Rayleigh , black solid line), true extent of degradation (D Sample , cyan dashed line), determined for short well screens at the bottom of the aquifer at two distances from the source (2 m, (a), and 95 m, (b)); true extents of TCE mass reduction at the respective distances across depth of the aquifer-aquitard-aquifer system (D System , red dashed line), and across depth of the aquitard (D Aquitard , blue dashed line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) therefore, our simulations provide a higher range of the overall impact of the aquitard on the CSIA data. Our assumption of biodegradation occurring throughout the aquitard (either uniform, or non-uniform as linear decreasing rate constants) yielded different CSIA patterns in the aquifer than under the assumption of a thin bioactive zone at the top of the aquitard except for the jump in isotope ratios after source release. Note that resulting CSIA patterns could be similar in case of abiotic RD. CEs can degrade abiotically in contact with ferrous minerals in aquitards made of clayey soils or rock matrix (Schaefer, 2016;Schaefer et al., 2017). While several parameters such as the magnetite susceptibility or the contact surface might influence degradation rates (Schaefer, 2016;Schaefer et al., 2017), those parameters are not necessarily related to the aquitard depth and degradation rate constants might be more uniform across such an aquitard. Clearly more research is needed to the depth distribution of bioactivity and chemical reactivity in aquitards both for interpretation of CSIA data and for assessing remediation time frames.
The aquifer could also be sufficiently oxic to oxidize the daughter products diffusing out of the aquitard. The detection of enrichment of δ 13 C-TCE in the absence of daughter compounds could be confused for TCE oxidation in the aquifer. Provided chlorine isotope enrichment factors for TCE cometabolism are available, the study of combined carbon and chlorine CSIA data could prevent this uncertainty (Abe et al., 2009;Wiegert et al., 2012). For instance, as the slopes of DCE carbon and chlorine CSIA data differ for biotic and abiotic reductive degradation, combined carbon and chlorine CSIA was employed to determine DCE degradation pathways at a field site (Audí-Miró et al., 2015). The study of carbon and chlorine isotope ratios is therefore promising for pathway distinction in this case.
Finally, the ε Diffusion values of DCE, VC, and ethene were estimated based on the literature. The variations of ε Diffusion had little impact on both the CIMB and θ. Provided they are in the selected range, an improved estimate of ε Diffusion values will not change the implications of the simulations. Additionally, a lower τ app value as suggested by Wanner and coworkers for the aquitard would reduce the distance travelled in the aquitard by the CEs within the simulation time and reduce the diffusion-induced isotope effect.

Summary and implications
Scenario simulations performed with a reactive transport model illustrate the potential of CSIA data to detect and quantify naturally occurring reductive dechlorination of chlorinated ethenes (CEs) in an aquifer-aquitard-aquifer system. The model input parameters were based on an extensive literature review. During presence of a TCE DNAPL source, TCE sequential reductive dechlorination to ethene in the underlying reactive aquitard enhances TCE diffusion from the aquifer to the aquitard. A resulting δ 13 C-TCE enrichment of a few permil (Δ 13 C > 2‰) is detectable in the aquifer near the aquitard, while in the absence of degradation, δ 13 C-TCE enrichment induced by TCE diffusion to the aquitard is negligible. This signifies that before source removal, large δ 13 C-TCE enrichment is unlikely to be attributable to degradation only in the aquitard, but would correspond to degradation at the aquifer level.
After source removal, the diffusion of all CEs out of the aquitard is accompanied by a quick and steady increase of the δ 13 C of all CEs and ethene with time, and a decreasing δ 13 C of all CEs with distance from the source. These characteristic CSIA patterns during the loading and releases phases contrast strongly to those associated with reductive dechlorination occurring in the aquifer only. These characteristic aquitard degradation patterns remain visible when degradation co-occurs in the aquifer and could serve as a line of evidence for degradation in the aquitard. The sudden jump in isotope ratios in the aquifer after source removal was recognizable in most of the simulations during the sensitivity analysis, suggesting that practitioners could use this line of evidence at many sites for effective DNAPL source removal (Wanner et al., 2018b).
Our simulated CSIA patterns differ from those found in the modelling study of Wanner et al. (2018b) where it was assumed that major TCE degradation only occurred in centimeters to decimeters thick bioactive zones at the top of aquitards. Under those conditions δ 13 C-TCE did not further increase over time after the initial jump that occurred after source removal and δ 13 C-TCE also slightly increased instead of decreased downgradient. These aquifer CSIA patterns are characteristic for aquitards where bioactivity drops to zero already at shallow depth. In a more recent study it was shown that bioactivity can spread deeper into an aquitard and rate constants decrease near linearly below a certain depth (Wanner et al., 2018b). We demonstrate that also in the case of non-uniform degradation similar CSIA patterns are expected as in the case of uniform degradation provided that the aquitard remains sufficiently bioactive (or reactive in the case of abiotic degradation) until CEs spreading depths. CSIA patterns during the release phase are thus also informative over aquitard reactivity with depth. Note that our predicted characteristic spatiotemporal CSIA patterns assuming relatively steady reactivity with depth could also be representative for abiotic chemical reductive dechlorination in aquitards where chemical reactivity as content of reactive minerals is likely to be relatively uniformly distributed.
The carbon isotope mass balance (CIMB) is influenced by the spatial separation of the CEs and ethene during transport. During the release phase, the CIMB increases in the upper aquifer as a result of the relatively limited back-diffusion of depleted daughter products. While this is the case also when degradation occurs only in the aquifer, the CIMB increase is stronger when degradation occurs only in the aquitard. The CIMB increase also occurs but to a lesser extent when VC is the end product of degradation, as such it only occurs for a short period when DCE stalls. The CIMB increase is strongest for thinner sorbing aquitards. Unfortunately for degradation monitoring, such increase could be confused with the degradation of the CEs through other pathways. Overestimating the daughter products' degradation is problematic as daughter products in the aquitard can act as a new contaminant source.
During the loading phase, application of the Rayleigh equation to δ 13 C-TCE data sampled in the aquifer strongly underestimates TCE degradation in the aquitard, provided the aquifer is unreactive. Within only a few years after source removal the precision of Rayleigh-based degradation quantification is improved to reasonable levels of underestimation (< factor 1.25), and seems to well represent the extent of mass destruction in the aquitard at this stage. Unfortunately, conditions increasing remediation timeframes such as a higher f oc in the aquitard, thicker aquitards, or a slower groundwater velocity increase degradation underestimation, especially before source remediation. Sorptionand diffusion-induced isotope fractionations have limited impact on TCE degradation assessment.
Scenario modelling provides insights into the potential of aquifer CSIA data for monitoring the reactivity of an aquitard. The study of a simplified conceptual site could be extended to more complex aquifer geometry, for instance multiple reactive clay beds of different thickness and lengths. Further work, including experimental work, should address the potential of chlorine and hydrogen isotope fractionation in identifying aquitard degradation and also to distinct between biological and abiotic degradation that might co-occur in aquitards (Wanner et al., 2018a). Further experimental research to both biological and chemical aquitard reactivity with depth is highly needed as this largely impacts resulting CSIA patterns in the release phase. Vincent Post for providing support with the Python interface FloPy of PHT3D, and Guillaume Chereau for providing technical assistance in high performance computing.