Virtual experiments to assess opportunities and pitfalls of CSIA in physical-chemical heterogeneous aquifers

Degradation of chlorinated ethenes (CEs) in low conductivity layers of aquifers reduces pollution plume tailing and accelerates natural attenuation timeframes. The degradation pathways involved are often different from those in the higher conductive layers and might go undetected when only highly conductive layers are targeted in site assessments. Reactive transport model simulations (PHT3D in FloPy) were executed to assess the performance of dual carbon and chlorine compound specific stable isotope analysis (CSIA) in degradation pathway identification and quantification in a coupled physical-chemical heterogeneous virtual aquifer. Degradation rate constants were assumed correlated to the hydraulic conductivity: positively for oxidative transformation (higher oxygen availability in coarser sands) and negatively for chemical reduction (higher content of reducing solids in finer sediments). Predicted carbon isotope ratios were highly heterogeneous. They generally increased downgradient of the pollution source but the large variation across depth illustrates that monotonously increasing isotope ratios downgradient, as were associated with the oxidative component, are not necessarily a common situation when degradation is favorable in low conductivity layers. Dual carbonchlorine CSIA performed well in assessing the occurrence of the spatially separated degradation pathways and the overall degradation, provided appropriate enrichment factors were known and sufficiently different. However, pumping to obtain groundwater samples especially from longer well screens causes a bias towards overestimation of the contribution of oxidative transformation associated with the higher conductive zones. As degradation was less intense in these highly conductive zones under the simulated conditions, overall degradation was underestimated. In contrast, in the usual case of limited CSIA data, dual CSIA plots may rather indicate dominance of chemical reduction, while oxidative transformation could go unnoticed, despite being an equally important degradation pathway.


Introduction
The potential degradation of chlorinated ethenes (CEs) through both biotic and abiotic redox-sensitive pathways in groundwater (Vogel et al., 1987) opened opportunities for monitored natural attenuation (MNA) as an alternative to invasive, energy-intensive technologies for the remediation of polluted sites. Until recently, MNA was focused on biotic reductive dechlorination (BRD) of TCE. Recent studies suspect that chloroethene degradation in the presence of iron(II) minerals might have been underestimated at field sites, including at sites where degradation occurs through simultaneous abiotic and biotic degradation (Brown et al., 2007). One of the reasons would be the general lack of detection of degradation products, either due to being indiscernible from the background concentrations or simply because they are not measured (He et al., 2009). For instance, when trichloroethylene (TCE) is in contact with the iron hydroxides green rust or magnetite, the daughter products acetylene, ethene, and ethane are likely to degrade further in aquifers (Berns et al., 2019;Han et al., 2012;Jeong and Kim, 2007;Lee and Batchelor, 2002;Liang et al., 2009). For similar reasons, the aerobic cometabolism of TCE, which ultimately produces CO 2 and Cl, could have been overlooked.
Compound specific stable isotope analysis (CSIA) has been demonstrated to be a valuable tool for determining CE degradation pathways and extent of degradation (Bloom et al., 2000;Clark et al., 2016;Palau et al., 2014;Sherwood Lollar et al., 1999). For instance, during TCE degradation, the bond containing a heavy carbon isotope is cleaved slightly slower than the bond containing a light carbon isotope. The resulting carbon isotope fractionation depends on the extent of TCE degradation and on the fractionation potential of the reaction, characterized by the process-specific enrichment factor ε (Hunkeler et al., 2002). Both degradation by iron(II) minerals and the aerobic cometabolism of TCE are associated with a large range large carbon isotope enrichment factors (respectively from −9‰ to −33‰ for chemical reduction (Liang and Dong, 2007;Zwank et al., 2005) and from −1‰ to −20‰ for aerobic cometabolism (Barth et al., 2002;Chu et al., 2004) and could be detected and quantified at field sites using CSIA.
When two degradation pathways are co-occurring at a site, the CSIA data of a unique element does not allow a proper attribution of degradation to the one or the other pathway (Höyng et al., 2015). Dual carbonchlorine CSIA (C-Cl CSIA) is specifically promising for the study of coexisting degradation pathways of CEs at field sites, as the relative extent of C and Cl isotope fractionation depends on the reaction mechanism. Where several processes were likely to occur, C-Cl CSIA slopes have helped in selecting the main degradation pathway (Badin et al., 2016;Clark et al., reaction pathway are known and the relative chlorine to carbon isotope enrichments (εCl/εC) are sufficiently different, the respective participation of each of the two different degradation pathways to overall degradation as well as the extent of overall degradation can be determined based on a C-Cl CSIA dataset ( van Breukelen, 2007).
Physical and chemical heterogeneities of aquifers directly impact reactive transport of pollutants. Physical heterogeneities induce preferential flow and therefore mixing of reactive components. Geochemical heterogeneity induces zones of either abundance or scarcity of reactive minerals and sedimentary organic matter. It is likely that zones of low conductivity will be more reductive, while conductive zones, where oxygen can be transported faster than it is consumed, will be more oxic. Additionally, different geological layers might present different mineralogy, specific surface areas, or host different microbial groups, inducing an additional correlation between hydraulic conductivity and reactivity. For example, at a field site sampled with high resolution, the reduction and therefore the immobilization of uranium was higher in low permeability zones than in permeable sediments, resulting from both abiotic and biotic processes (Janot et al., 2015). However, hydraulic conductivity is not the only parameter expected to be correlated with reactivity. For instance, the ferrous mineral content seemed to be an important parameter to explain the difference of reactivity of two clayey soils in the abiotic degradation of chlorinated ethenes (Schaefer et al., 2017). In both low-conductivity zones and in conductive zones, availability of substrate and bacterial activity plays a role in the biotic reduction or oxidation of pollutants (Berns et al., 2019;McMahon, 2001;Wanner et al., 2018a;Yan et al., 2016). With time, reduced conditions develop in chlorinated solvent plumes (Christensen et al., 2000). Such additional conditions influence reactivity while being potentially uncorrelated to the hydraulic conductivity, and were addressed in previous modeling studies (Atchley et al., 2014;Jang et al., 2017;Sanz-Prat et al., 2016). Despite those additional complexities, heterogeneity of hydraulic conductivity in aquifers plays a major role in degradation (Jang et al., 2017).
Well-characterized aquifers are rare, and aquifers are often considered as 'black boxes' in which the degradation and transport patterns of CEs must be determined through spatially and temporally sparse datasets. In this context, synthetic data sets created using models combining reactive transport with isotope fractionation are employed for advising practitioners on the optimal CSIA sampling strategies (Höyng et al., 2015;Xu et al., 2017). The effects of hydraulic conductivity heterogeneity on resulting C-CSIA patterns of toluene degrading with various electron-acceptors was assessed before by reactive transport modeling of aquifer analogs (Höyng et al., 2015). Using such a similar approach, the present study investigates whether dual C-Cl CSIA could detect TCE degradation through two different reaction pathways when those pathways occur dominantly in different conductivity zones. Low-permeability sediments such as clay lenses are more likely to present reductive conditions susceptible to host biotic reductive dechlorination of CEs (Damgaard et al., 2013;Takeuchi et al., 2011;Wanner et al., 2018b), and might offer large reactive iron mineral contents and surfaces (Schaefer et al., 2017). High-conductivity zones as permeable sands are more likely to be aerobic and therefore can host TCE biotic oxidation, either directly or through aerobic cometabolism. Applying a simple correlation between the hydraulic conductivity and the two degradation rate constants of TCE (aerobic transformation and chemical reduction), alike the approach of (Cunningham and Fadel, 2007), we investigate the resulting spatially heterogeneous carbon and chlorine isotope ratio patterns of TCE and the implications for polluted site investigation. The main objectives of this study are (1) to understand how spatially heterogeneous CSIA patterns evolve from physicalchemical heterogeneous conditions, (2) to assess the performance of dual C-Cl CSIA in detecting the occurrence of both degradation pathways and in calculating overall degradation under these conditions, and (3) to assess the effect of pumping (preferentially from the more permeable zones) on degradation assessment.

General setup of model
The multicomponent transport and degradation is modeled with PHT3D (Prommer et al., 2001), which combines the 3-D multicomponent transport model MT3DMS with the geochemical model PHREEQC-2 (Parkhurst and Appelo, 1999). PHT3D was accessed through the python interface FloPy (Bakker et al., 2016). The thirdorder total-variation-diminishing (TVD) solution scheme was used for the advective term, which tends to minimize numerical dispersion and artificial oscillation. Model run time on a computer with an Intel CPU at 2.5 GHz and 4 GB of RAM was about 9400 min and the produced data file size was of about 5 GB.
A TCE source with a thickness of 4 m is set at the upgradient boundary of the aquifer. The source includes a non-degradative tracer with otherwise similar characteristics to TCE. The degradation rate constants of aerobic cometabolism and chemical reduction of TCE are either positively (oxidation) or negatively (chemical reduction) correlated to the hydraulic conductivity (K). The first-order degradation rate constants range between 0 and 10 y −1 (discussed in detail later).
The simulated 2D aquifer domain is 100 m long and 10 m thick, with a homogeneous rectangular discretization of 0.5 m in the x (spreading) direction and 0.05 in the z (depth) direction. For simplicity, the porosity was set homogeneously to 0.35 across the aquifer. This value is typical for unconsolidated sands. The dry bulk density was set to 1.855 kg·L −1 throughout the aquifer system.
The heterogeneous hydraulic conductivity field was obtained using the field generator available in ModFlow (Chiang, 2013), which is based on Mejia's algorithm, which was previously used for hydrogeological modeling studies (Zammouri and Ribeiro, 2017). In our model, the mean of the log-normally distributed K was set to 3.3·10 −4 m·s −1 , with a standard deviation for ln(K) of 1.5, in line with prior modeling studies (Atchley et al., 2014;Uçankuş and Ünlü, 2008) and field data (Gelhar and Welty, 1992). Heterogeneity is supposedly lesser in the horizontal than in the vertical direction in sediment layers, which is obtain in the generated K field by setting correlation lengths to 7 m longitudinally and 1 m vertically. Those values are in of the order of magnitude of published correlation lengths at the as relatively homogeneous considered the Borden site (0.91-8.3 m in the horizontal direction, 0.07-0.34 m in the vertical direction) (Maghrebi et al., 2015).
Those input parameters yield a medium heterogeneous aquifer (see Fig. 1) with hydraulic conductivities corresponding to a sandy aquifer (Gelhar and Welty, 1992) and are in line with prior modeling studies (Atchley et al., 2014;Uçankuş and Ünlü, 2008). The simulated conductivity field is presented in Fig. 2a and the distribution of K is presented in Fig. 1. The flow direction is from left to right, with an imposed average hydraulic gradient of 0.0033 m·m −1 calculated to achieve an average groundwater velocity of 30 m·y −1 based on the geometric mean of the simulated hydraulic conductivity of 3.2·10 3 m·y −1 . The simulation duration is of 6 years, at which steady-state is reached. The longitudinal and vertical dispersivities were set to 1 m and 5 mm, respectively. The selected grid spacing (mentioned above) and time step (1.1 days) ensure limited numerical dispersion. TCE aqueous diffusion coefficient D w was set to 5.1·10 −5 m 2 ·d −1 , calculated after (Worch, 1993). To account for tortuosity, the effective diffusion coefficient was calculated as D eff = D w × τ app , where τ app is the apparent tortuosity, which was taken equal to the porosity therefore τ app = 0.35. This value is in the low range for τ app values (Carey et al., 2016). Diffusion effects on TCE isotope fractionation in a heterogeneous aquifer are likely to be more important at the plume fringe in case of "plume-fringe degradation" (van Breukelen et al., 2005) whereas in the present model only "core degradation" is simulated. Additionally, such isotope effects are best captured by short wells screens (Thouement et al., 2019;van Breukelen et al., 2005). In this study the effects of diffusion were expected to go undetected as sampling is simulated through long screens in the core of the plume. Diffusion related isotope effects were therefore not simulated. Diffusion was expected to have negligible effects on TCE isotope fractionation under the simulated conditions. Similarly, sorption is likely to influence isotope ratios only on the spreading front of pollution plumes and was not further simulated.

Correlation between hydraulic conductivity and degradation rates
The modeling study relies on the hypothesis is that the spatial distribution of degradation pathways is connected to the hydraulic conductivity in the aquifer. Recent detailed studies of subsurface sediments in the Netherlands show that clays present a generally higher potential reactivity than sands, with notably higher contents of reactive iron, pyrite, and organic matter in clays versus sands (Griffioen et al., 2016;Griffioen et al., 2012). The more reduced clay zones are therefore likely to provide the necessary reactive iron minerals such as pyrite required for abiotic reductive dechlorination (Lee and Batchelor, 2002;Liang et al., 2009;Weerasooriya and Dharmasena, 2001). Biotic reductive dechlorination of CEs has often been detected in reductive high permeability zones. When occurring in clay, biotic degradation of CEs might profit from elevated contents of organic matter which serve as substrate to the bacterial population (Wanner et al., 2018b), to the point that biotic degradation is favorable in the low conductive zones compared to the conductive zones (Wanner et al., 2016). Finally, TCE aerobic (co)metabolism, for instance on methane (Fogel et al., 1986), requires the presence of oxygen but might benefit from the vicinity of low-permeability layers provided they act as substrate source for the microbial population in the otherwise less reactive sand zones (McMahon, 2001). We assume that more permeable parts of the aquifer tend to contain higher levels of oxygen and thus may entail higher rates of aerobic (co)metabolism. In this paper, the first-order degradation rate constants were therefore correlated to the spatially varying hydraulic conductivity following (Cunningham and Fadel, 2007). Despite the simplification of a perfect correlation, we expect that the employed correlation is sufficient for addressing our research objectives.
The first-order rate constant of TCE oxidation (OX), k OX , positively correlated with the hydraulic conductivity, is such as ln(k OX ) = f(ln(K)). The first-order rate constant of chemical reduction (CR), k CR , negatively correlated to the hydraulic conductivity, is such as ln(k CR ) = f(− ln (K)). For the rate constants to evolve between 0 y −1 and 10 y −1 , they were then defined at the aquifer scale such as lnk OX = a OX ⋅ ln (K) + b OX and lnk CR = a CR ⋅ ln (K) + b CR , with each pair of a and b values being calculated as follows: Fig. 1. On the left axis, the distribution of the hydraulic conductivity K (grey histogram) and the proportion of the oxidative, k OX , to chemical reductive, k CR , first-order rate constants (black). On the right axis, the change of degradation rate constants for TCE chemical reduction (blue) and TCE oxidation (red) against the decimal logarithm of K (K in m·s −1 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) where min ln(K) , min -ln (K) , and max ln(K) , max -ln(K) , represent respectively the minimums and the maximums of the values of the ln(K) permeability field and of -ln(K). The resulting distribution of the rate constants k OX and k CR are presented in Fig. 2b and c, respectively, and in relation with K in Fig. 1. Across the simulated aquifer, the average and median of the pathways rate constants are of the same order of magnitude (not shown). TCE chemical reduction is slightly more favorable than oxidation due to the skewness of the K values distribution (Fig. 1). Although the model assumptions based on (Cunningham and Fadel, 2007) lead to co-existence of reduction and oxidation, always one pathway (strongly) dominates. Only at the median permeability they occur at equal but very low rates. Note that we simulated a medium heterogeneous sand aquifer which not specifically included clay lenses but finer grained and more impermeable sections assumed to have a similar reactivity for CR as clay/silts as known from the literature. As a heterogeneous sandy aquifer was modeled, it was chosen not to increase the standard deviation and/ or to lower the mean of K to produce zones with K values closer to clay hydraulic conductivity, as this would also have yielded more zones with K values representative of gravel. The outcomes of this study are not expected to change by this assumption.

Simulation of isotope fractionation
The TCE source is modeled as a constant TCE concentration at the upgradient aquifer boundary. The initial TCE concentration is set to about half the solubility, i.e., 600 mg·L −1 . The initial TCE carbon and chlorine isotope ratios are −30‰ and + 2.5‰, respectively, which corresponds to the isotope signature of industrially produced TCE (Shouakar-Stash et al., 2003). The initial concentration of each carbon and chlorine isotope is calculated following (Hunkeler et al., 2009) as presented in the Appendix in section A1. Concentrations below the usual isotope quantification limit (1 μg/L) are not reported.
First-order degradation kinetics was applied to both reaction pathways. Both carbon and chlorine isotope fractionation were modeled independently through the "bulk isotope" method. For each isotopic element (C and Cl), both the light ( l TCE) and heavy ( h TCE) isotopes in the TCE molecule were simulated. The difference in degradation rate constants associated with the heavy and the light isotopes, h k bulk and l k bulk , respectively, is described by the kinetic isotope fractionation factor α k = h k bulk / l k bulk , where α k = ε + 1. Isotope fractionation is incorporated in the reactive transport model by setting the degradation rate constants of the heavy and the light isotopes, h r and l r, respectively, as follows Van Breukelen et al. (2005): where r is the overall first-order degradation rate of TCE associated with a specific process: where k is the first-order TCE degradation rate constant of a specific process.
The magnitude of degradation-induced isotope effects is reported as the bulk effect, i.e., the average of all reactive and non-reactive positions of a molecule. Bulk isotope fractionation is sufficient for simulating isotope enrichment of parent compounds (Badin et al., 2018;Cretnik et al., 2014;Hunkeler et al., 2002;van Breukelen et al., 2017van Breukelen et al., , 2005. Bulk isotope enrichment factors of TCE CR were set to εC CR = −15‰ and εCl CR = −2.5‰ for carbon and chlorine, respectively (Audí-Miró et al., 2013), which corresponds to degradation by zero-valent iron. For TCE OX, bulk isotope enrichment factors were set to εC OX = −11.5‰ and εCl OX = +0.3‰ (Gafni et al., 2018), which corresponds to TCE aerobic cometabolism by the toluene degrader Pseudomonas Putida F1. These values correspond to εCl/εC slopes of +0.16 for TCE CR and − 0.03 for TCE OX. Different εCl/εC slopes associated with each pathway as selected here are necessary to enable pathway identification with dual C-Cl CSIA ( van Breukelen, 2007).
The slightly lower degradation rates of the heavy isotopes results in the total summed rate of all isotopes to be somewhat lower than the intended degradation rate of the molecule. As the reaction networks for carbon and chlorine isotopes are independent, the predicted C-Cl CSIA slopes gradually deviate from reality over time. Following Van Breukelen et al. (2017) a correction factor was applied to each isotope reaction rate in order to align the two independent isotope networks with the overall reaction progress (Appendix Section A2).

Virtual groundwater samples
The TCE concentrations and isotope ratios were calculated for virtual groundwater samples obtained from various virtual monitoring well networks. The concentration of TCE or of the tracer in each virtual sample, C (sample) , were calculated as the pumping flow-rate-weighted mean concentrations of the model grid cells representing the immediate surroundings around the well screen (Höyng et al., 2015). It was assumed that the pumping local flow rate was proportional to the local hydraulic conductivity. This method takes into account the influence of the local hydraulic conductivity on the groundwater flux into the well screen during pumping: where C i is the concentration of the compound and K i the hydraulic conductivity of a grid cell i. The isotope ratio of element E (in this study, carbon or chlorine) in the sample is weighted according to the respective concentrations of TCE in each grid cell: where δ (sample) and δ i are the heavy to light isotope ratio of the element E, of the sample and of a cell i, respectively, expressed in the δ notation δ (sample) = (R (sample) − R 0 )/R 0 , and where R (sample) and R 0 are the ratios of the heavy to light isotope concentrations of an element E in a sample and of the international standard, respectively.

Degradation estimates and rate ratio
Degradation estimates were obtained using the previously developed modified Rayleigh equation for assessing degradation through two different processes using dual isotope analysis ( van Breukelen, 2007). With this concept, degradation can be attributed to the one or the other pathway, or to a combination of those, based on the dual C-Cl CSIA data, provided that their respective carbon and chlorine isotope enrichment factors εC and εCl are known. The proportion of the influence of the one pathway over the other is quantified through the rate ratio F as defined in Van Breukelen (2007).
The derivation of the modified Rayleigh equation is briefly summarized below. We refer to Van Breukelen (2007) for the details. Based on the commonly used Rayleigh equation, the isotopic shift of an element E, Δ, is calculated as follows: where δ source and δ (sample) are the initial isotope ratio of the element E in the source and a downgradient sample, respectively, ε is the enrichment factor for the element E for a given reaction pathway, and f is the fraction of remaining (non-degraded) TCE. The ratio of the apparent enrichment factors Φ is estimated based on the ratios of the isotopic shifts: For our study, the rate ratio F of our two competing pathways is calculated as: where εC CR , εC OX , εCl CR , and εCl OX are the carbon and chlorine isotope enrichment factors of the chemical reduction and the aerobic cometabolism pathways, respectively. The distribution of the two degradation pathways can therefore be directly obtained from the relative isotopic shifts and from the enrichment factors of the two pathways, provided they are known. The apparent carbon isotope enrichment factor εC app is calculated as follows: The modified Rayleigh equation to quantify the remaining TCE fraction f using dual C-Cl CSIA is then:  Fig. 3 illustrates the impact of the coupled physical and chemical heterogeneity on the spatial distribution of the TCE concentrations, the carbon and chlorine isotope ratios, and the degree of overall degradation within the plume at steady-state. TCE concentrations are the highest in the high-conductivity zones which induce preferential flow (Fig. 3b). TCE carbon and chlorine isotope ratios show limited enrichment in the high-conductivity zones compared to the low-conductivity zones due to relatively lower residence times (Fig. 3d, e). Both carbon and chlorine isotopes display a net mean increase with distance but with a large variation with depth (Fig. 4a). The observed patterns differ from a prior study of heterogeneous toluene degradation where degradation and consequently isotope enrichment was concentrated at the fringes of the plume (Höyng et al., 2015). Similar results were obtained in the case co-existence of reductive and oxidative processes was not allowed (i.e. by putting the rate of the minor process to zero; results not shown).

General patterns
Two additional simulations in which either TCE chemical reduction (CR) or TCE oxidation (OX) occurs reveal a major difference in associated isotope patterns (Fig. 4bc). While the standard deviation of δ 13 C values at a certain depth profile for TCE OX is smaller than 0.5‰ regardless of the travel distance (Fig. 4b), for TCE CR the standard deviation is large and varies between 2‰ and 10‰ (Fig. 4c).
Interestingly, despite the heterogeneity of TCE OX degradation rate constants, isotope patterns for TCE OX (Fig. 4b) are similar to the monotonous enrichment expected in case of uniform degradation in the aquifer (Hunkeler et al., 1999). The contrasting large depth variation of isotope ratios associated with TCE CR is sustained by the inverse correlation of k CR with the hydraulic conductivity. The implications are broader as this pattern might be observed when any degradation pathway yielding large isotope enrichment (for instance biotic reductive dechlorination) is favorable in the lesser permeability layers. As degradation mostly occurs in the low conductivity zones, those zones will show high isotope enrichments especially when conductivity is relatively low and contaminants reside there for a long time. In contrast, limited isotope enrichment is observed in the connecting high permeability zones where degradation by CR is limited. The correlation between the chemical and the physical heterogeneity plays thus a determining role in the distribution of isotope enrichment.
Using the synthetic datasets, the true state of TCE degradation, D True , is calculated at the point scale. D True shows a pattern similar to the pattern of carbon and chlorine isotope ratios (Fig. 3c). In the conductive zones of the aquifer, D True increases gradually with distance from the source, reaching up to 40% degradation, while low-conductivity zones are hotspots of intense degradation with D True varying from 60% to 80%.
The share of TCE degraded through oxidation versus reduction is estimated using the rate ratio, F, calculated using the isotope ratios δ 13 C and δ 37 Cl for each point of the model through eq. (10) to (12). The rate ratio F varies between 0 and 1, with 0 indicating that TCE degradation is entirely attributed to CR and 1 that TCE is entirely degraded through OX. At the source zone, the limited carbon isotope enrichment (< 0.5‰) did not allow for calculating F. Despite zones of intense TCE oxidation, F reaches at most 0.73 in the simulated domain. F approaches 0 only in the upgradient zone of the aquifer where the aquifer is relatively impermeable and CR dominates. Downgradient, F bears the influence of the alternation of CR and OX across the aquifer and presents no extreme values as F varies approximately between 0.2 and 0.7. F as calculated from the isotope ratios poorly compares to the real local rate ratio k OX /(k OX + k CR ) which was set as input (Fig. 5b). F is therefore not an indicator of the local conditions of degradation but rather an integrator of the degradation processes acting on the pollutant from the source zone until the sampling location.

Sampling effect
Field investigation often relies on water samples obtained through pumping from (long) well screens. This approach might induce an overestimation of the importance of the degradation pathway (dominantly) occurring in the high-permeability zones as relatively more water is pumped from these zones compared to the low-permeability zones. Under the assumption that the obtained synthetic dataset is representative of realistic transport and degradation behavior of TCE, both the degradation extent and the contribution of each degradation pathway to TCE degradation were evaluated from synthetic sample CSIA data.
Dual C-Cl CSIA data plots are commonly used to determine degradation pathways likely to occur at field sites. The proportion of isotope enrichment of dual CSIA data is compared to compiled laboratory enrichment factors in order to identify which pathway(s) may be taking place. Dual C-Cl CSIA data are plotted for each of the 6 monitoring networks (each consisting of 5 wells placed at equal depth and with 20 m spacing) differing in their start position from the source (10 m or 20 m) and well screen elevation (3 m, 5 m, or 7 m) (Fig. 6). The slopes of δ 37 Cl/δ 13 C linear regressions at each of the monitoring networks (+0.04 to +0.12) are always larger than the theoretical εCl/εC slope of TCE OX (−0.03). A first result is that independently from the positioning of the monitoring network, the sampled C-Cl CSIA data are not only influenced by the process occurring in the high-but also by the process in the low-conductivity zones, and this despite the preferential flow through the high-conductivity zones. For practitioners, this implies that the presence of a degradation process occurring in low-permeability zones might be revealed by a more or less randomly set monitoring network on the condition that the long contact time in low permeability zones does not lead to complete degradation.
A potential drawback for the assessment of degradation is that the small positive (inverse) chlorine isotope enrichment factor (+0.3‰) associated with TCE OX is overridden by the normal chlorine isotope enrichment factor (−2.5‰) induced by TCE CR (Fig. 6 and Fig. 4a). Each set of CSIA data displays strong chlorine isotope enrichment which could result in TCE OX going undetected. For instance, the slope of the network presented in Fig. 5d (+0.12 ± 0.04) is not significantly different from the εCl/εC slope for TCE chemical reduction (+0.16).
The slopes values are in the range of recently published TCE εCl/εC slopes during biotic reductive dechlorination (BRD) obtained after the microorganisms were precultivated on lesser chlorinated CEs in a laboratory experiment (range:~0.05-0.11, published as εC/εCl slopes of 9.0 ± 1.1 to 18.2 ± 4.3 (Lihl et al., 2019)). Those growth conditions and therefore the small slopes are not expected to be common at field site. However, when the observed slopes are lower than the previously published sets of εCl/εC values both for TCE biotic and abiotic reductive degradation (Fig. 6), co-occurrence of oxidation and chemical reduction pathways may be suspected. The present simulation shows that the resulting εCl/εC slope is in between the slopes of the reductive and the oxidative pathway. BRD as met so far at field sites has a much larger εCl/εC slope than CR. A similar scenario with BRD in the low conductivity layers could produce δ 37 Cl/δ 13 C slopes in the range of published εCl/εC slopes of CR or the lower end of BRD. In that case, dual C-Cl CSIA would not be suitable for pathway distinction. The formation of the metabolites DCE and VC might inform that biotic degradation occurs unless they are consumed by oxidation. Additionally, this underlines that the εCl/εC slopes of the degradation pathways must Fig. 3. Cross section showing a case of combined TCE oxidation (in high-permeability layers) and chemical reduction (in low-permeability layers) at steady state of (a) K, (b) TCE concentration, (c) the true TCE degradation extent D True , and of the isotopic shifts (d) Δ 13 C and (e) Δ 37 Cl. The TCE source is indicated by the thick black line at the upgradient (left) boundary. The 6 different regularly spaced monitoring networks consisting of 5 wells of 1 m long screens placed at equal depth (3 m, 5 m, or 7 m) and with 20 m spacing differing in their start position from the source (10 m (black rectangles) or 20 m (grey rectangles)). At 50 m and 90 m from the source the grey contours correspond to two 4-m-long wells with average screen elevation of 3 m and 7 m. be sufficiently different to be able to use C-Cl CSIA data for pathway distinction as pointed out in Van Breukelen (2007). Currently, natural attenuation through aerobic cometabolism of TCE is considered to be limited, but could have been overlooked due to the lack of specific degradation products. In the line of recent research to determine lines of evidence of TCE aerobic cometabolism and rate constants characterization (Wilson et al., 2019), dual C/Cl CSIA might be an additional line of evidence for judging of the relevance of TCE aerobic degradation at a field site. Moreover, this approach could be applied to other compounds which can degrade through both aerobic and anaerobic conditions such as BTEX compounds or the TCE daughter compounds DCE and VC.
In order to illustrate the potential of dual CSIA data to inform on the importance of the two pathways occurring at the field site, the calculated rate ratio, F, as obtained from simulated flow-rate-weighted well data, F (sample) , is compared to the average of F in the pore water along the well screen, F (core) . The values are compared at two selected distances, 50 and 90 m, and for two well screen sizes, 1 m and 4 m (Fig. 7c, 7d). When sampling through long well screens, the importance of TCE oxidation in the aquifer was overestimated as result of sampling (F (sample) > F (core) ) as groundwater is preferentially abstracted from the more permeable layers. Overestimation was the most at the lower plume fringe in zones of lowconductivity (bottom long wells, Fig. 3), and was the least at the top wells screens, which are set in the core of the TCE plume and in high-conductivity zones (top long wells, Fig. 3) (Fig. 7c, 7d). Despite the overestimation, the F (sample) values (< 1) always indicate some co-occurrence of TCE CR (Fig. 5, 7c, 7d). For short wells, F (sample) matches well with F (core) . Interestingly, some limited underestimation (F (sample) < F (core) ) can be noticed for the top short well screen at 90 m from the source (Fig. 7d), likely due to TCE back diffusing from a low-conductivity zone a few meters upgradient from the screen (Fig. 3). In this modeling study, the virtual samples correspond to optimal sampling conditions. Under real conditions, insufficient pumping rates or durations might lead to the abstracted water to originate even more from the high-conductivity zones (McMillan et al., 2014), in which case overestimation of the importance of degradation processes associated with high-conductivity zones when sampling across low-conductivity zones is likely the rule.
The CSIA-based extent of degradation at point scale in the aquifer, D Rayleigh , is calculated from applying eq. (9) to (12) to the local δ 13 C and δ 37 Cl values. When plotting D Rayleigh against depth, the D Rayleigh curves slightly underestimate D True (Fig. 7e, 7f). The underestimation is a result of dispersion which fades the isotope ratio gradients resulting from degradation (Abe and Hunkeler, 2006;van Breukelen and Prommer, 2008). For the presented cross-sections, the underestimation of D True is in the range of the previously estimated underestimations at field sites and of less than 5% (Abe and Hunkeler, 2006). The limited difference between D Rayleigh and D True suggests that flow segregation as described   Thouement and B.M. Van Breukelen Journal of Contaminant Hydrology 231 (2020) 103638 in (Kopinke et al., 2005) had little impact on the results as flow segregation limits the applicability of CSIA for degradation estimation. A similar simulation with larger reactive zones or rate constants is likely to exemplify flow segregation in which case CSIA-based degradation estimates would underestimate degradation (Kopinke et al., 2005). The overestimation of the occurrence of TCE OX in the sample illustrated by F (sample) > F (core) (Fig. 7e, 7h) suggested that the εC app calculated for the groundwater sample εC app(sample) would be biased towards εC OX . As |εC OX | < |εC RD |, εC app(sample) is expected smaller in absolute value than εC app calculated from the average pore water isotope values (εC app(core) ). Such underestimation of |εC app | would overestimate degradation extent (D Rayleigh(sample) > D True(core) ), of which the opposite is observed as also the difference in isotope enrichment between the groundwater sample and the pore water needs to be taken into account. With the example of the deepest long well at 90 m from the source, which shows the most contrast between the sample and the core degradation values (Fig. 7 h), the sampled isotope values (δ 13 C (sample) = −23.1‰ (Fig. 7c), δ 37 Cl (sample) = 3.1‰) yield a εC app (sample) of −13.3‰. In comparison, the apparent enrichment value for pore water is indeed slightly larger, with εC app(core) = −13.6‰. However, the carbon isotope ratio in the sample (−23.1‰) is notably less enriched than in the core (δ 13 C (core) = −21.4‰). When applying the smaller εC app(sample) to the less enriched δ 13 C (sample) , the resulting D Rayleigh(sample) is underestimating degradation as occurred in the pore water (D (core) = 52%, D Rayleigh(sample) = 41%). The underestimation reflects the lesser degradation in the highly conductive zones (degradation of 30% to 40%) compared to the lesser conductive zones (up to 70%) as shown by the depth profiles at 80 m from the source (Fig. 7h, the conductivity field is given in Fig. 3a). If degradation would be more intense in the highly conductive zones compared to the lesser conductive zones, pore scale degradation would be overestimated. Finally, the slight underestimation of D True(sample) by D Rayleigh(sample) is of 3%, in the range of the underestimation seen for the depth profiles D True and D Rayleigh . The underestimation is likely due to dispersion. As shown by the smaller difference between D Rayleigh(sample) and D True(core) for smaller well screens or for the top long well (Fig. 7h), the ability to describe degradation in the aquifer depends strongly on the sample method and the lithology crossed by the well screen.

Conclusion and implications
This simulation study illustrates that reductive degradation of CEs in low-permeability zones might be detected using dual carbon (C) and chlorine (Cl) CSIA, also when TCE oxidation occurs in highly conductive zones. Simulated degradation extent and isotope enrichments were heterogeneous, with the lesser conductive zone forming pockets of more intense degradation compared to the more conductive zones. When heterogeneous aquifers are sampled, large variations of C CSIA values across depth and non-linear increase with distance could indicate that degradation is favorable in the less permeable layers. The δ 37 Cl/δ 13 C slopes of the various virtual monitoring networks were intermediate between the published sets of εCl/εC slopes for TCE Fig. 6. Dual element isotope plot (δ 37 Cl (sample) vs. δ 13 C (sample) ) representing simulated well data at steady-state for the 6 virtual monitoring networks consisting of 5 wells each, with the corresponding calculated linear regression (thick black line, slopes of the linear regression (slp.) indicated in the plot). Slopes of the colored lines correspond to selected sets of εCl/εC for biotic reductive dechlorination in the literature (thin green lines) and to the input parameters (thick red line: (1) TCE OX, −0.03 (Gafni et al., 2018); thick blue line: (2) TCE CR, +0.16 (Audí-Miró et al., 2013)). Relative εCl/εC for TCE biotic reductive dechlorination (BRD) (3) Dehalococcoides: +0.21 (Kuder et al., 2013), (4) Geobacter lovely strain SZ: 0.3 (Cretnik et al., 2014), (5) mixed bacterial culture: +0.37 (Wiegert et al., 2013). The error bars on the well samples correspond to realistic uncertainties for field data: ± 0.5‰ for carbon and ± 1‰ for chlorine. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) chemical reductive dechlorination and TCE oxidation. When the εCl/εC slopes of the pathways occurring at the site are known and sufficiently different, the contrast with the sampled data could point towards the presence of two degradation pathways that are simultaneously occurring. The results could be extended to other degradation pathways, for instance in the case of favorable biotic reductive dechlorination in the low permeability zones instead of chemical reduction. Biotic reduction of TCE is associated with larger εCl/εC compared to chemical reduction. The resulting δ 37 Cl/δ 13 C from the combination of oxidative and reductive biotic pathways could match previously published εCl/εC slopes for chemical reduction and may prevent pathway distinction.
Additionally, the effect of sampling bias should not be neglected. The preferential abstraction from the more conductive layers led to the overestimation of TCE oxidation to overall TCE degradation in the aquifer. Despite this wrong attribution of degradation, dual CSIA from the sampled data provides a good estimate of the degradation state of the sampled TCE provided εC and εCl of both occurring pathways are known. Limiting sampling across layers of difference conductivities is critical to obtain a good estimate of degradation at the pore scale. Degradation might otherwise be over-or under-estimated, depending on whether the conductive layers are more or less reactive than the lesser conductive zones. In practice, the bias induced either by the sampling method or by a limited amount of CSIA data should be taken into account when designing sampling campaigns or interpreting data in view of determining to which extent various degradation pathways are taking place.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. (e,f) CSIA-based rate ratio F, and (g,h) TCE extent of degradation: true (D True ) and CSIA-based (D Rayleigh ). For all plots, the value/parameter is given across depth at the point scale (thin black lines, blue dashed line for D Rayleigh ) and for flow-rate-weighted groundwater samples (dashed red vertical lines). Average pore water values along well screens are represented with thick black vertical lines for (e,f) F (F (core) ) and (g, h) D True (D True(core) ). Values from flow-rate-weighted groundwater samples are represented with dashed red vertical lines for (a, b) TCE concentration, (c,d) F, and (e,f) D Rayleigh . The true extent of degradation at the groundwater sample level D True(sample) is indicated by grey vertical lines in panels e and f. Results are shown at 50 m (left panel) and 90 m (right panel) downgradient from the source. At both locations, the shallow wells have been placed in more permeable zones than the deep wells. The positions of the shallow and deep wells are shown in Fig. 3.