Modeling 3D-CSIA data_ Carbon, chlorine, and hydrogen isotope fractionation during reductive dechlorination of TCE to ethene

Reactive transport modeling of multi-element, compound-specific isotope analysis (CSIA) data has great potential to quantify sequential microbial reductive dechlorination (SRD) and alternative pathways such as oxidation, in support of remediation of chlorinated solvents in groundwater. As a key step towards this goal, a model was developed that simulates simultaneous carbon, chlorine, and hydrogen isotope fractionation during SRD of trichloroethene, via cis-1,2-dichloroethene (and trans-DCE as minor pathway), and vinyl chloride to ethene, following Monod kinetics. A simple correction term for individual isotope/isotopologue rates avoided multi-element isotopologue modeling. The model was successfully validated with data from a mixed culture Dehalococcoides microcosm. Simulation of Cl-CSIA required incorporation of secondary kinetic isotope effects (SKIEs). Assuming a limited degree of intramolecular heterogeneity of δCl in TCE decreased the magnitudes of SKIEs required at the non-reacting Cl positions, without compromising the goodness of model fit, whereas a good fit of a model involving intramolecular CeCl bond competition required an unlikely degree of intramolecular heterogeneity. Simulation of H-CSIA required SKIEs in H atoms originally present in the reacting compounds, especially for TCE, together with imprints of strongly depleted δH during protonation in the products. Scenario modeling illustrates the potential of H-CSIA for source apportionment.


Introduction
At many contaminated sites, monitored natural attenuation (MNA) of chlorinated ethenes is the preferred and cost-effective remediation approach (Meckenstock et al., 2015). Microbial sequential reductive dechlorination (SRD) of chlorinated ethenes is usually the main transformation process in MNA. The dechlorination proceeds from the primary contaminants tetrachloroethene (PCE) and/or trichloroethene (TCE), to daughter products cis-and trans-1,2-dichloroethene (cDCE and tDCE), 1,1-dichloroethene (1,1-DCE), vinyl chloride (VC), and finally to nontoxic ethene (ETH). Degradation may, however, also occur via alternative transformation pathways such as (an)aerobic oxidation (Bradley, 2011;Bradley and Chapelle, 2011;Chu et al., 2004;Pooley et al., 2009) and chemical reduction (Damgaard et al., 2013;Darlington et al., 2013;Ferrey et al., 2004;Lee and Batchelor, 2002) of lower and higher chlorinated ethenes, respectively. At field sites, the efficacy of SRD can be verified by quantitation of the degradation products that are pathway-specific. However, assessment of the alternative pathways of chlorinated ethene destruction is more difficult, since the degradation products (Cl − , CO 2 ) blend with the natural background levels. Moreover, degradation of SRD products can lead to an underestimation of the SRD performance as it seems that parent compounds are not well degraded because product concentrations are low. Consequently, less sustainable remedies, such as pump and treat, may be instituted or continued unnecessarily.
Compound-specific stable isotope analysis (CSIA) has been applied in contaminant studies, to detect and identify degradation processes (Hunkeler et al., 2008). One complication in interpretation of CSIA data is potential occurrence of sequential/parallel transformations, e.g., SRD followed by aerobic degradation of SRD products (Arp et al., 2001;Mundle et al., 2012). Multi-dimensional CSIA, i.e., CSIA of more than one element (C, Cl, H), as started with 2-D carbon and hydrogen CSIA of MTBE to discriminate anaerobic from aerobic transformation, holds particular promise also for detecting individual degradation pathways of chlorinated solvents. In addition to a growing number of reports on combined C and Cl isotope effects in various chlorinated solvent degradation systems (Abe et al., 2009;Audi-Miro et al., 2013;Cretnik et al., 2013;Wiegert et al., 2012) several CSIA applications focus specifically on discrimination of such alternative reaction mechanisms for chlorinated ethenes (Badin et al., 2016;Badin et al., 2014;Dogan-Subasi et al., 2017).
Studies of Cl isotope fractionation in reactions of organochlorine compounds focus on the so-called primary isotope effects, i.e., isotope fractionation of the atoms positioned at the reacting molecular bonds (Paneth, 2006). In accordance with that consensus, the existing models of Cl isotope fractionation were developed for isotope fractionation involving only the primary effects Jin et al., 2013). However, it was also suggested that reactions of chlorinated hydrocarbons (Palau et al., 2014a) may result in fractionation of nonreactive Cl atoms (secondary isotope effects) or possibly combine both primary and secondary effects (Abe et al., 2009;Cretnik et al., 2014;Palau et al., 2014a). Höhener (2016) recently extended the analytical BIOCHLOR model with calculation of carbon and chlorine stable isotope ratios in chloroethenes and accounting for secondary chlorine isotope effects in the TCE to cDCE step. Likewise, to address the potential contributions from Cl isotope effects occurring at non-reactive positions, the present numerical model considers isotope effects occurring at all Cl positions within the reacting molecules.
Recently,  presented a dataset on SRD of TCE to ETH by a Dehalococcoides culture, including 3D-CSIA (C, Cl, H) results, at high temporal resolution. This dataset allows for further validation and improvement of the chlorine isotopologue fractionation model developed by Hunkeler et al. (2009) and the development and validation of a hydrogen isotope fractionation model. The objectives of the current study were (i) to extend the current chlorine isotope fractionation model  as adopted in subsequent studies (Höhener, 2016;Jin et al., 2013;Palau et al., 2014b;Wiegert et al., 2012) for isotope effects occurring at multiple Cl positions of TCE (Höhener, 2016) and of DCE, and for intramolecular heterogeneity in δ 37 Cl of the source compound (TCE); (ii) to develop a completely novel model of hydrogen isotope fractionation during SRD; and (iii) to validate the models with the experimental data . The validated models requiring Monod kinetics are expected to form a template for CSIA interpretation of SRD of halogenated hydrocarbons and contribute towards CSIA-based support of remediation of chlorinated solvent groundwater pollution at field sites.

Experimental data for model validation
Model validation used data from a microcosm experiment on dechlorination of TCE by a Dehalococcoides (Dhc) culture. A detailed description of that experiment is available elsewhere . In summary, microcosms were set up with Bio-Dechlor Inoculum (BDI) culture (Amos et al., 2008) a consortium of Dhc strains that is capable of complete dechlorination of PCE via TCE, DCE, and VC to ETH. The microcosms were amended with TCE and lactate electron donor. Concentrations and C, Cl, and H isotope ratios were determined for TCE and the aforementioned reaction products over the course of degradation.

Nomenclature of the chlorine and hydrogen isotope effects
Parameters describing the magnitude of isotope effects use the rate constants (k) for heavy vs. light isotope species, where: lightk/ heavy k = 1/α; ε = α − 1. The fractionation factor (α) and enrichment factor (ε) can describe position-specific or "bulk" effects. The latter, indicated by the "bulk" subscript (e.g., ε bulk ) are averaged for all reactive and non-reactive positions of a molecule. KIE (kinetic isotope effect) is the ratio of light k and heavy k for a specific molecular position and a specific transformation process. In a primary KIE, the isotope substituted-position is at the reaction center. A secondary KIE (SKIE) is the effect of isotope substitution at a position remote from the reaction center (Elsner et al., 2005).
In the discussion, we refer to position-specific isotope effects, to reflect the different fates of individual atoms of the reacting compounds. We refer to ε RP (at "reactive position", i.e., the Cl position undergoing dechlorination) or as ε NRP (at "non-reactive position", i.e., the Cl position not undergoing dechlorination). Isotope effects for the Cl atoms remote from the dechlorination site (e.g., α and βt for TCE in Fig. 1a) are by definition secondary KIEs (SKIEs). Observable isotope effects for the Cl atom undergoing dechlorination may be in fact primary or secondary, depending on the dechlorination process involved. The latter may occur if the initial transformation step does not involve CeCl cleavage. It was previously proposed that the initial step in abiotic dechlorination of VC and cDCE by cobalamin is nucleophilic addition of cobalamin and the chloroethene species (Glod et al., 1997). A similar process was postulated for biodegradation of the same species by Dhc (Abe et al., 2009;. It is also possible that the Fig. 1. (a) A scheme illustrating chlorine isotope effects that occur during SRD of TCE to ETH. Note that the isotope effects at the non-reacting positions (NRP) are defined with respect to the location of the reacting position (RP): α is the germinal position, β is/are the vicinal position(s); c (cis) and t (trans) represent the stereochemical position relative to the reacting position. In TCE transformation, SKIEs (ε NRP ) occur of type α and βt, while a KIE or SKIE (ε RP ) occurs at the dechlorination site (see text for additional explanations). (b) A scheme of hydrogen isotope fractionation in SRD of TCE to ETH, combining SKIEs (in TCE transformation, type βc) and the effects of inserting H atoms by protonation, ε Hprotonation (see below and also Eqs. 5 and 6). For simplicity, the reactions involving tDCE are omitted. nucleophilic addition mechanism may apply to biodegradation of TCE . In such reactions, the only observable isotope effects may in fact be secondary KIEs of the addition step.
For simplicity, the model treats the ε RP as a primary KIE (see Supplementary Table S2), and, in the narrative, we will refer to the effect at the reactive position as KIEs and the effects at the non-reactive positions as SKIEs; however, the model structure can equally accommodate reaction scenarios with SKIEs at all Cl positions of the parent compound and no observable primary KIEs.
Recently, Cretnik et al. (2014) proposed an alternative model of TCE dechlorination, where cDCE resulted from conversion of tDCE product. In that pathway, both germinal Cl atoms would contribute to cDCE. That mechanism will be discussed in detail in Supplementary data Section S4.
In the case of H isotope fractionation, the change of H isotope ratios of TCE over the progress of reaction is only controlled by the (secondary KIE) ε NRP at the single H atom of that compound. In the remaining reactions, the net changes of isotope ratio of the daughter compounds combine ε NRP (H atoms transferred from the parent compounds) and the isotope composition of the H atom added in protonation (ε Hprotonation ; Fig. 1b). Protonation is discussed more extensively in Section 3.4.1.

Reaction kinetics
The model was developed with the PHREEQC code (Parkhurst and Appelo, 1999). SRD of TCE to ETH, via DCE and VC, was simulated. cis-DCE was the main DCE isomer produced in the microcosm, but minor quantities of trans-DCE and 1,1-DCE were detected . For model simplicity, the sum of the latter two DCE isomers was explicitly simulated as trans-DCE. Two minor pathways were, therefore, added to the model: TCE to trans-DCE and trans-DCE to VC. Monod kinetics was applied without microbial growth (Bekins et al., 1998) for all reactions: where Rate m is the reaction rate of molecule m (− Rate m is production rate of its daughter product), v m is the substrate utilization rate (M·L − 3 T − 1 ), C m is the concentration of the molecule (M·L − 3 ), and K s is the half-saturation constant (M·L − 3 ). Individual lag periods, i.e., the period before the reaction in question began (T), were used for all reactions. The selection of kinetic parameter values is discussed in Section 3.1. A time step of 1 h was taken.

Simulation of carbon isotope fractionation
The "bulk isotope" method was applied (Van Breukelen et al., 2005). In other words, for each compound, the light and the heavy isotope species were defined reflecting the compound's concentration fraction of the light and heavy isotopes, respectively. The concentration of an isotope species was taken to be equal to its fraction multiplied by the compound's concentration. Reaction rates were as follows: where Rate L and Rate H are the rates of the light and heavy isotopes, respectively, C L and C H are the concentrations of the light and heavy isotopes, respectively, and ε bulk is the bulk kinetic isotope enrichment factor of the reaction step. Isotope ratios were calculated from the simulated concentrations of the light and heavy isotopes.

Simulation of chlorine isotope fractionation
The isotopologue approach was applied ( ); see Fig. 2) which considers all isotopologues in the reaction network, and, for TCE, also all Cl isotopomers (i.e., isotopologues with same number of heavy isotopes but located at different positions). This model developed by Hunkeler et al. (2009) was extended in the current study to also account for isotope effects at the positions away from the dechlorination center (SKIEs). The model also addresses the possibility of intramolecular δ 37 Cl heterogeneity of TCE and the potential of intramolecular CeCl bond competition (IBC) in TCE degradation ( Fig. 2; Supplementary Section S4).
The initial Cl isotopologue/isotopomer concentrations assuming either intramolecular homogeneity or heterogeneity in δ 37 Cl of initial TCE were calculated as described in Supplementary Section S2. Isotope ratios were calculated from the simulated concentrations of the isotopologues/isotopomers. The reaction rate for each isotopologue/isotopomer (Rate mi ) was obtained by: Fig. 2. Schematic of the reaction network on chlorine isotope fractionation following the isotopologue approach showing all chlorine isotopologues and isotopomers during reductive dechlorination of TCE. The PCE isotopologues are shown for sake of completeness. Blue solid arrows indicate (the approximate chance of) reaction pathways without reaction rates being affected by 37 Cl at reacting positions (but with potential effect of 37 Cl at non-reacting positions), whereas red dashed arrows indicate (the approximate chance of) reaction pathways with rates being affected by 37 Cl at the reacting position (together with potential additional effects of 37 Cl at non-reacting positions). The two boxes show the TCE isomers containing 1 and 2 heavy chlorine isotopes, respectively. Note only the chlorine atoms in yellow boxes react in the formation of cDCE. The chlorine atoms in green boxes and the green dashed arrows indicate a mechanism postulated by Cretnik et al. (2014) were both geminal chlorine substituents of TCE are reactive in the formation of cDCE (intramolecular CeCl bond competition (IBC); see section S4). Potential secondary isotope effects are not shown in this schematic but given in Table  S2.
where Rate m (Eq. (1)) is the degradation rate of the corresponding parent compound, C i (M·L − 3 ) is the concentration of the isotopologue/ isotopomer of interest, n i is the number of Cl atoms at reactive positions (undergoing CeCl cleavage) in the isotopologue/isotopomer, H i is the number of heavy chlorine isotopes at reactive positions in isotopologue/isotopomer i, and Π(ε (N)RPi + 1) or Πα (N)RPi is the multiplication of the inverse (α) of the applicable primary KIEs and SKIEs during transformation of isotopologue/isotopomer i to either one or two daughter isotopologues. A single daughter isotopologue is formed if where n m is the total number of chlorine atoms in the molecule (Fig. 2). A single daughter isotopologue is thus formed during transformation of entirely heavy or entirely light isotopologues and for all TCE isotopomers. The first part of the last term between brackets of Eq. (4) describes the chance that a Ce 35 Cl bond of an isotopologue is broken and the potentially secondary KIEs that apply. The second part of this term describes the chance that a Ce 37 Cl bond is broken and all primary and/or secondary KIEs that apply. The following conversions were used: α (N)RP − 1 = ε (N)RP; α (N)RP = 1/(S)KIE; (S)KIE = 35 k/ 37 k; and ε bulk ≈ Σε (N)RP /(n m /x) (Elsner et al., 2005), where x is the number of these atoms in reactive positions (1 for all reactions). Note, the sum of ε (N)RPs is nearly equal to their multiplication product (in α (N)RPs equivalents; Table S2).

Simulation of hydrogen isotope fractionation
The bulk hydrogen isotope ratio of TCE is affected only by the SKIE at the single H atom present (Fig. 1b). The bulk δ 2 H of the daughter products is affected predominantly by the isotope signatures of the H atoms incorporated during dechlorination/protonation ( Fig. 1b) (Ertl et al., 1998;Shouakar-Stash et al., 2003). Hydrogen isotope ratios were simulated with an extended "bulk isotope" method (see Section 2.3.2). To simulate δ 2 H of a daughter product, the model considered (i) isotope fractionation of the H atoms transferred from the parent to daughter product and the subsequent daughter product (Eqs. (2)-(3), where ε bulk results exclusively of SKIEs (ε bulk = ε NRP(MEAN) ), and (ii) the rates of the light, Rate 1H , and the heavy, Rate 2H , H isotopes replacing the Cl atom of the parent compound, i.e., through protonation, at each dechlorination step calculated as the total rate multiplied by the light and heavy H isotopic abundance, respectively: where Rate m (Eq. (1)) is the degradation rate of the corresponding parent compound, the terms between parentheses in Eqs. (5)-(6) are the isotopic abundances of light and heavy hydrogen isotopes, respectively, ϵ Hprotonation is the overall hydrogen isotopic enrichment factor expressed with respect to δ 2 H water and associated with this reaction step, and VSMOW is the ratio of 2 H/ 1 H of the international standard for the H isotopic composition of water. Note that the values of ϵ Hprotonation and δ 2 H water as input for Eqs. (5) and (6) are not converted to permil, following the IUPAC recommendation for isotope ratio notation (Coplen, 2011 (1)-(6) describe the degradation rates of the individual isotopes/isotopologues as part of the molecule. Note that the sum of their rates based on Eqs. (1)-(6) is somewhat lower than the intended rate of the molecule as described by Monod kinetics in Eq. (1). Because the heavy isotopes/isotopologues react slightly slower than the light isotopes due to isotope fractionation, the total summed rate of all isotopes/ isotopologues in the molecule based on Eqs. (1)-(6) is lower than the rate given by Eq. (1). Furthermore, the overall reaction progress varies slightly among the independent (C, Cl, H) isotope reaction networks since the heavy isotope abundances and fractionation effects vary for the different elements and reaction pathways. Therefore, the predicted isotope ratios of different elements at a certain moment in time are not ideally aligned with the overall reaction progress. This effect is irrelevant in assessment of individual isotope networks, but for a multidimensional isotope analysis, predicted slopes (e.g., Δ 37 Cl/Δ 13 C) may deviate somewhat from reality under extreme conditions (Jin et al., 2013), although this was not relevant for the current simulation.
The three independent isotope networks were aligned with the overall reaction progress with the following simple correction term applied to each isotope/isotopologue reaction rate: where A H and A L are the abundances of the heavy and light isotope, respectively, and α is the bulk isotope fractionation factor in the case of C, and the multiplication of α (N)RPs for Cl and H. Values of correction terms consequently varied per isotope pair and compound. This correction term is in principle similar as the corrections done by Hunkeler et al. (2009) for first-order chlorine isotope fractionation. Note that as the heavy isotopes react a factor α slower, the abundance of the heavy isotopes multiplied with α together with the abundance of the light isotopes (multiplied with 1) (= A H × α + A L ) represents the extent that the summed isotope/isotopologue network reacts slower than aimed with Eq. (1). Therefore, multiplying the individual isotope/isotopologue rate equations with the inverse of (A H × α + A L ) makes the sum of the individual isotope/isotopologue rate equations to become identical to the intended Monod kinetics rate as defined in Eq. (1). With this correction, compound concentrations (sum of their isotopes/isotopologues) were identical among C, H, and Cl isotope networks over the entire course of reaction. Therefore, simulation of multielement isotopologues (e.g., 12 C-13 C-37 Cl-37 Cl-1 H-2 H as one of the 27 multi-element C/Cl/H DCE isotopologues) as performed for combined CeCl isotopologues of chlorinated ethenes by Jin et al. (2013) albeit an elegant and sophisticated approach, was not required, avoiding the dramatic increase in the number of isotopologues and reactions to be modeled.

Metrics of model fit
The goodness of model fit was evaluated following, e.g. Karlsen et al. (2012), by means of computing the root mean squared weighted error (RMSWE): where w i , obs i , and sim i are the weights, observed, and simulated values, respectively, of observation i, and n obs is the total number of observations. Weights are included so that different groups of data (concentrations, C-/Cl-/H-CSIA data) can be directly compared (dimensionless) based on the estimated observation error. Mean of the error is used so that datasets with different amount of observations can be compared. The weights are calculated for each observation as: where the accuracy is taken as the measurement error of the CSIA data (C-CSIA: 0.5‰; Cl-CSIA: 0.8‰, TCE and VC; 1.0‰, cDCE; H-CSIA: 20‰) or the observed concentration value multiplied by the measurement error (CEs: 7%; ETH: 15%). Thus observations with higher accuracy get a higher weight and vice versa. Note that the value 1.96 in Eq. (9) is the approximate value of the 97.5 percentile point of the normal distribution, i.e., the 95% confidence interval lies within roughly 1.96 standard deviations of the mean. Consequently, a RMSWE of 1.96 means that the model on average deviates from the observations with a value equal to the accuracy of the observations.

Reaction kinetics
Using the model described, the substrate utilization rate, v m , the half-saturation constant, K s , and the lag period were determined for each pathway of the SRD reaction network by manually fitting the model to the concentration data. Table 1 lists all kinetic parameter values. A reasonably good agreement was obtained for all concentration observations (Fig. 3a). As metrics of goodness of model fit, the root mean squared weighted error (RMSWE) was computed for all simulations and for each parameter (see Table S7).
The transformation of TCE to DCE started almost immediately. The overall linear concentration decline reflects zero-order kinetics, which follows from Monod kinetics with a low K s relative to the TCE concentration. However, the initial decay rate was very slow and could not be simulated with this kinetic model. The lag period had to be set longer (2.4 days) than actually observed. Consequently, the first δ 37 Cl-cDCE observation after 1.1 days was not simulated. TCE dechlorination produced mostly cDCE but minor quantities of tDCE and 1,1-DCE were determined and were part of the total mass balance. Their summed molar concentration ranged between 4 and 9% with respect to cDCE. Due to their relative recalcitrance, they became proportionally more important when cDCE was almost gone, but then their concentrations were already below 1 μM. The model accounted for the small effect of the transient presence of these minor DCE isomers under the simplification that only tDCE occurred. The model did, however, not describe the shape of the cDCE concentration peak around day 5 well. This may relate to the slight mass balance variations observed during this period. Moreover, co-production of cDCE, tDCE, and 1,1-DCE greatly increased the number of model variables that required calibration for this reaction step.
Conversion to VC was complete after 8 days. Reductive dechlorination of VC to ETH was slow in the subsequent period but high rates occurred after 41 days. The model fitted VC and ETH concentrations well. The relatively high K s compared to the previous reaction steps (200 μM versus 8-33 μM; Table 1) resulted in pseudo first-order degradation kinetics of VC. The concentration and C isotope mass balances indicated that no further conversion of ETH occurred.

Carbon isotope fractionation
The model fitted C-CSIA of all compounds very well (Fig. 3b). The C isotope enrichment factors of the TCE to DCE transformation (ϵC TCE → cDCE and ϵC TCE → tDCE , − 16.4 ± 0.4‰) and of the VC to ETH step (ϵC VC → ETH , − 26.7 ± 1.9‰) were taken equal to the Rayleigh equation derived values . Enrichment of the cDCE to VC step (ϵC cDCE → VC , −26.8‰) was derived by  from the difference between δ 13 C-cDCE and the initial δ 13 C-VC following Hunkeler et al. (1999), whereas enrichment of the tDCE to VC step (ϵC tDCE → VC , −30.3‰) was taken similar as Hunkeler et al. (2002). The C isotope enrichment factors were consistent with previous results for Dhc and for biotic RD in general , and the observed trends of evolution of δ 13 C over time, for TCE and the daughter products were typical for SRD (Van Breukelen et al., 2005).

. Secondary KIEs in chlorine isotope fractionation
In the early stage of transformation, the Cl isotope ratios showed clear offsets of δ 37 Cl among the chlorinated ethenes (Fig. 3c). Instantaneously produced cDCE and VC were thus more depleted in δ 37 Cl than their precursors. In SRD, as illustrated in Fig. 1a, Cl atoms at the CeCl bonds undergoing dechlorination are split off the ethene skeleton of the reacting molecule, while the non-reactive Cl atoms are transferred to the daughter compounds. Therefore, if the bulk chlorine isotope effect consists exclusively of a primary KIE (at the reactive position), then isotope fractionation does not affect the Cl atoms transferred to the daughter compounds. However, in the case where the bulk chlorine isotope effect also includes contributions from SKIEs, the δ 37 Cl value of the daughter chlorinated ethene compound should differ from that of the precursor, and the difference should be equal the average SKIEs of the reaction, as postulated by Hunkeler et al. (2009). Indeed, inclusion of SKIEs in the model was required to adequately simulate the evolution of δ 37 Cl values of all reaction products during SRD, as explained below. An alternative explanation for the TCE to cDCE δ 37 Cl offset might be intramolecular heterogeneity (IH) in δ 37 Cl of initial TCE (discussed in Section S4).
The position-specific isotope effects and bulk enrichment factors of all reaction steps are listed in Table 2 (determined experimentally in the microcosm study by  with the exception of the α and βc isotope effect in the transformation of TCE to tDCE). The isotope effects are explained as follows. First, in the transformation of TCE to cDCE, the mean of the SKIEs (type α and βt; ϵ NRP(MEAN) , −3.3‰) was calculated from the difference between δ 37 Cl-TCE and initial δ 37 Cl-cDCE observed ; since only the average of the α and βt effects could be determined experimentally, they were considered as equal in the model. This was done for the sake of model simplicity even though SKIEs in the β position tend to be smaller than those in the α position (Elsner, 2010;Elsner et al., 2005); the βc effect (in the TCE to tDCE step) was assumed to be of equal magnitude (ϵ NRP(βc) = −3.3‰). Second, the reactive position effect in the TCE to cDCE step (ϵ RP , − 4.2‰) followed from 3 × ϵ (bulk) − 2 × ϵ NRP(MEAN) , where ε bulk was − 3.6 ± 0.3‰ . Third, the SKIE in the cDCE to VC step (ε NRP(βc) = − 1.7‰) followed from the difference between δ 37 Cl-cDCE and initial δ 37 Cl-VC observed . Fourth, the reactive position effect in the cDCE to VC step (ϵ RP , − 4.5‰) followed from − 2 × (δ 37 Cl-VC final − δ 37 Cl-TCE initial ) + ϵ NRP(βc) (cf. Eq. S1 in ). The βt SKIE (ϵ NRP(βt) ) of the tDCE to VC step was assumed equal to the βc SKIE (ϵ NRP(βc) ) of the cDCE to VC step. Fifth, as only a ϵ RP occurred in the VC to ETH step, its value followed directly from fitting the Rayleigh equation to the observations (ϵ RP = ϵ bulk , − 2.7 ± 0.4‰; ).
The extended model using parameters derived from direct observations fitted the microcosm Cl-CSIA data very well. Note that the addition of SKIEs to the model of Hunkeler et al. (2009) implied that isotope fractionation occurred for nearly all chlorine isotopologue reactions except for those that only contained light chlorine atoms (Table  S2). This confirmed the conclusion from  where the SKIEs (in a nucleophilic addition reaction) were postulated for all SRD reactions.

Effect of intramolecular heterogeneity on chlorine isotope patterns
Intramolecular heterogeneity (IH) in δ 37 Cl of TCE could be an alternative to SKIEs in explaining the initial offset between δ 37 Cl-TCE and δ 37 Cl-cDCE. In manufacturing of TCE, chlorinated organic compounds of dissimilar chlorine isotope signatures may be combined, resulting in different position-specific chlorine isotope ratios . In the event that the non-reactive chlorine positions in the TCE to cDCE reaction have a lower average δ 37 Cl versus the δ 37 Cl-TCE, cDCE produced at the outset of the transformation will be 37 Cl-depleted relatively to the TCE. Note that IH cannot account for the offset between δ 37 Cl-cDCE and δ 37 Cl-VC since the transformation of DCE to VC is not Cl position-specific. Fig. 4 (and also Fig. S1) present various simulations to illustrate the effect of IH and the assumed absence or presence of SKIEs. Table 3 (and  also Table S3) presents an overview of the main characteristics of these simulations. A baseline model incorporating only primary KIEs (Table 3) does not describe the observations well, and results in δ 37 Cl values of daughter products always equal to or higher than the initial value of TCE (Fig. 4a). Note that the two non-reactive chlorine atoms of TCE become part of cDCE and, therefore in the absence of SKIEs, instantaneously produced cDCE always has the same δ 37 Cl as source TCE. Similarly, unless the SKIE is present, the initially produced VC cannot have a lower δ 37 Cl than the cDCE precursor.
The model including IH, but no SKIEs, reproduces the first few δ 37 Cl-cDCE observations but does not fit the last two δ 37 Cl-cDCE data points and all δ 37 Cl-VC data (Fig. 4b). Furthermore, this model requires the assumption of a wide difference in the isotopic ratios for different positions, i.e., a TCE molecule in which the reacting position had a δ 37 Cl of +9.8‰ and the two non-reacting positions had δ 37 Cl values of − 0.1‰ ( Fig. 4b; Table 3). As discussed by  such a large IH is unlikely because the required δ 37 Cl value of +9.8‰ at the reactive position is unreasonably large in comparison with the isotope ratios of chloride evaporites used as industrial chlorine sources and the range of δ 37 Cl reported for synthetic organochlorine compounds (− 5 to Table 2 Calibrated isotope enrichment factors (‰).

Reaction
Carbon Chlorine Hydrogen . c Note only the average of ϵ NRP(α) and ϵ NRP (βt) can be determined and should equal ϵ NRP(MEAN) . d Taken as the observed difference between the parent and initial daughter isotope ratio. e Carbon and chlorine isotope effects were taken the same as for the dominant TCE to cDCE and cDCE to VC steps, respectively. However, the ϵC for the latter step was taken equal to the one found by Hunkeler et al. (2002). f ϵ RP of cDCE to VC follows from −2 × (δ 37 Cl-VC final -δ 37 Cl-TCE initial ) + ϵ NRP(βc) (cf. Eq S1 in  +6‰; (Hoefs, 2009)). The occurrence of SKIEs in RD of TCE and of cDCE was required to simulate the δ 37 Cl-cDCE and δ 37 Cl-VC observations (Fig. 4c-d). The final calibrated model assumed absence of IH for the sake of model simplicity ( Fig. 4c; Table 3), whereas the alternative model assumed a mild degree of IH and as result a larger difference between the reacting and non-reacting position effects in RD of TCE and of cDCE was required ( Fig. 4d; Table 3). The goodness of fit of these two models was similar in terms of RMSWEs (final versus alternative: 2.4 vs. 2.0 (δ 37 Cl-cDCE), 0.6 vs. 0.9 (δ 37 Cl-VC); Table S7). Therefore, absence of IH cannot be ascertained but the occurrence of SKIEs is required to provide a good model fit, particularly with respect to the initial depletion of VC relative to cDCE. Noteworthy, the alternative model assuming a mild degree of IH led to a perhaps more plausible larger difference in the deduced magnitudes of the isotope effects at the reacting and non-reacting positions.

Evaluation of intramolecular CeCl bond competition
We tested whether the mechanism of intramolecular CeCl bond competition (IBC) postulated by Cretnik et al. (2014) is consistent with our experimental data set. The model was extended as explained in the Supplementary Section S4 to describe the reaction of TCE to trans-DCE with selective interconversion of trans-DCE to cDCE. We assumed the rates of this route and of the normal TCE to cDCE route as equal. Cretnik et al. proposed IBC to explain an unexpectedly high level of position specific Cl isotope effects without invoking significant SKIEs at those positions. We observed that fitting the IBC mechanism to our data is indeed possible, but only for improbably high extent of IH (Supplementary Section S4). Our alternative model Fig. 4d (cf. Table 3) enables reduction of the ε NRP (from −3.3‰ down to −1.7‰) under the assumption of half to a quarter of the degree of IH needed for the IBC model (see Supplementary Section S4, Table S5), and therefore seems more probable.

Effects of protonation
As SRD of TCE progressed, a 'stepped' H isotope fractionation pattern was observed as the most striking feature, whereby each subsequent daughter product was more depleted in 2 H that its precursor, with ETH reaching δ 2 H values around − 270‰, whereas initial δ 2 H-TCE was + 530‰ (Fig. 3d). Furthermore, in the case of TCE and less for cDCE, δ 2 H continued to decline as the reaction progressed, and related to inverse secondary isotope effects (discussed later). A good agreement was achieved between the modeled and measured δ 2 H patterns (Fig. 3d).
The stepped decrease in δ 2 H values observed in the order TCE, cDCE, VC, and ETH is the result of the addition of a H atom at each dechlorination step that, on average, is strongly depleted in 2 H relative to both water (−42‰) and the original TCE (+530‰) (Kuder et al.,Fig. 4. Modeled isotope fractionation patterns for TCE, cDCE, and VC for the first 10 days (without further VC conversion) obtained by four different model assumptions: (a) primary KIEs only, with ϵ TCE,bulk = − 3.6‰ and ϵ cDCE,bulk = − 3.1‰ (both the same for tDCE production and consumption, respectively); (b) intramolecular heterogeneity (IH) of TCE in δ 37 Cl, with for TCE a δ 37 Cl of + 9.8‰ at the reacting position, and − 0.1‰ at the two nonreacting positions, and primary KIEs only as in model a; (c) the final calibrated model, including isotope effects at non-reacting positions but without IH (see text); and (d) the alternative model, with mild degree of IH (δ 37 Cl of + 6.4‰ at the reacting position, and + 1.6‰ at the two non-reacting positions) and larger differences between the reacting and non-reacting position effects of TCE (ε RP = − 7.4‰, ε NRP(MEAN) = − 1.7‰) and of cDCE (ε RP = −6.9‰, ε NRP = − 1.7‰).

Table 3
Model parameter values (‰) of simulations shown in Fig. 4.

Simulation Explanation
Initial δ H-TCE was  + 530‰), respectively, with sucrose and cane-sugar as electron-donor (δ 2 H-H 2 O was −60‰). The current study showed that the overall fractionation effects with respect to water, ε Hprotonation , increased with the extent of dechlorination (TCE → DCEs: − 170‰; DCEs → VC: -580‰; VC → ETH: -740‰). These values were obtained by manual calibration of the model and were similar as those obtained by  except for the TCE to DCE step (TCE → DCEs: −130‰; DCEs → VC: -590‰; VC → ETH: -750‰) who applied the following equation where δ 2 H addition is the δ 2 H of the newly added hydrogen atom in a reaction step, n is the number of hydrogen atoms in a given daughter product, "bulk" refers to the average δ 2 H of the daughter and parent compounds, e.g., cDCE and TCE. Subsequently, ϵ Hprotonation follows from δ 2 H addition -δ 2 H water . We ascribe the relatively large underestimation by Eq. (10) for ε Hprotonation of the TCE to DCE step to the large inverse hydrogen isotope effect of TCE SRD resulting in a much reduced offset between δ 2 H-TCE and δ 2 H-cDCE. Whereas Eq. (10) ignores this additional isotope effect, the model accounts for it and enables accurate quantification of protonation effects.
The δ 2 H of H added during protonation is the outcome of a complex series of fractionation processes. Dhc species used in the experiment require molecular hydrogen as the immediate electron donor . Consequently, the hydrogen atom replacing the Cl atom during RD derives from H 2 produced by the fermentation of lactate. The δ 2 H of produced H 2 should, therefore, depend on the δ 2 H of lactate and isotope fractionation effects during fermentation, both of which were unknown. In isotopic equilibrium, δ 2 H of produced H 2 is strongly depleted with respect to δ 2 H-H 2 O. Its equilibrium value in the current experiment would be −757‰ (α = (D/H) H2O /(D/H) H2 = 3.95 at 20°C (Horibe and Craig, 1995). In biological systems, isotopic equilibration is fast, provided hydrogenases are present (Campbell et al., 2009;Valentine et al., 2004). As an illustration, results from experiments in which D. autotrophicum grew on formate and where δ 2 H-H 2 O and δ 2 H-formate were varied independently showed that δ 2 H of the cells fatty acids was entirely controlled by δ 2 H-H 2 O (Campbell et al., 2009). The strongly depleted δ 2 H of the H atom inserted into SRD product (ε Hprotonation ) appear to mainly reflect the depleted δ 2 H of the H 2 . Therefore, we normalized ε Hprotonation with respect to ambient water.
While ε Hprotonation is likely a composite parameter (combining the depleted δ 2 H of the H 2 and possibly kinetic effects or effects associated with hydrogen insertion), model results showed we could simplify this complex series of fractionation steps with a single and constant processspecific overall 'protonation' fractionation effect, ε Hprotonation .

Secondary kinetic isotope effects in hydrogen isotope fractionation
Imprinted on the main pattern of decreasing δ 2 H of subsequent reaction products, a clear and surprisingly large inverse isotope effect (+ 34 ± 11‰)  was observed to occur during the transformation of TCE decreasing instead of increasing its δ 2 H (Fig. 3d). Since primary KIEs do not occur for H during SRD, the inverse isotope effect must therefore be an SKIE of type βc that occurs for the TCE to cDCE step (Fig. 1b). Inverse effects are typical for H atoms bonded to the C adjacent to the site of nucleophilic addition (Elsner et al., 2005), as in the proposed TCE reaction with cobalamin for the current experiment . SKIEs of the other dechlorination steps could be kept to zero to obtain a good model fit (Fig. S4: simulation S4a).
However, assuming occurrence of SKIEs in the cDCE to VC step (ε NRP(MEAN) = +10‰, average of SKIEs types α and βt) improved the fit for the last δ 2 H-cDCE observation (RMSWE 0.9 vs. 1.4) without changing the results for VC and ETH much (Fig. S4: simulation S4b). For the sake of consistency, the final calibrated model (Figs. 3d;S4) applied the same values for these three SKIEs for each reaction step (see Table S6). For example, an average SKIE (ε NRP(MEAN) = +14.7‰) of types α, βt (from cDCE), and βc (from TCE) was applied for the VC to ETH step. It must be stated that considering the analytical error of H-CSIA of ± 20‰ , it cannot be determined with certainty whether the SKIEs related to the cDCE and VC dechlorination steps truly deviate from zero. Note that the drop in δ 2 H-VC values around day 5 was captured well by the model (Fig. S4).
It remains to be seen how well the model is able to reproduce data from field sites where ε Hprotonation values may be less consistent than observed during the present lab experiment. Relatively stable and depleted δ 2 H-cDCE field values (− 211 ± 20‰, n = 10, one outlier excluded) pointing to RD of PCE seem promising in that respect (Audi-Miro et al., 2015). Further experimental studies are needed to test how ε Hprotonation varies with microbial culture, reaction rate, and temperature in order to obtain further mechanistic understanding of the magnitudes of these prime parameters affecting the δ 2 H offsets of chlorinated solvents and reaction products.

Exploring potential H-CSIA patterns in aquifers by means of scenario modeling
In order to explore the potential use of H-CSIA in source apportionment of TCE versus PCE source zones, we extended the H-CSIA model with the PCE to TCE step to assess the δ 2 H values of TCE and daughter products in scenarios of pure and mixed PCE and TCE sources (see Supplementary Fig. S5).

Model extension
The PCE to TCE step only involves the simulation of protonation of TCE. PCE was added as a molecule to the model and its degradation rate linked to the production and thus protonation rate of TCE. The extended PHREEQC model (see Fig. S5) was used in 1-D advection/dispersion transport mode during complete reductive dechlorination. The summed concentration of PCE and TCE in the source was 1 mmol/L in all simulations. The groundwater flow velocity was 20 m per year, the longitudinal dispersivity coefficient was 1 m, and the duration of simulations was 15 years. Table 4 shows the applied input parameter values selected for the degradation and H isotope fractionation processes. First-order kinetics was assumed with a set of degradation rate constants in line with previous modeling studies to SRD (Hunkeler et al., 2009;Van Breukelen et al., 2005). The values on hydrogen isotope fractionation were adopted from simulation S4a; Thus ϵ NRP(MEAN) values were zero except for TCE (+ 34‰). The value of ϵ Hprotonation was not known for the PCE to TCE step and was taken equal to the TCE to DCE step. δ 2 H-H 2 O was taken as −42‰. δ 2 H-TCE was taken as + 500‰, within the range (+ 467‰ to + 682‰) of published values for manufactured TCE .  Fig. 5 presents the 1-D model simulation results representing concentrations and H isotope patterns for complete dechlorination of PCE and/or TCE to ETH along the simulated flow path. Note that the sum of the CEs and ETH concentrations decreases beyond 250 m downgradient because of longitudinal dispersion with the displaced clean background groundwater.

Model results
The parent compound TCE becomes depleted in δ 2 H (Fig. 5: left: solid lines) during reductive dechlorination due to the inverse isotope fractionation as observed for this reaction step in this study. Daughter products are increasingly depleted the less they are chlorinated. During protonation strongly depleted hydrogen atoms replace the Cl atoms resulting in strong depletion of the final metabolite, ETH.
Note PCE does not contain H atoms and consequently H isotope ratios are not shown for PCE. In this case (Fig. 5: middle), strongly depleted δ 2 H-TCE is produced, about 550-700‰ more depleted than the source TCE of the previous simulation ( Fig. 5: left). Correspondingly, the other daughter products are also considerably more depleted than in the TCE as parent compound scenario. Note that the difference in δ 2 H between the two simulations (PCE versus TCE) decreases the less chlorinated the compound is. In this scenario of a PCE source, δ 2 H-DCE exceeds δ 2 H-TCE because of (i) the inverse H isotope effect during the TCE to DCE step; and (ii) the isotope fractionation effects associated with protonation are assumed similar for both the PCE to TCE and the TCE to DCE steps. As a result, both H atoms added during protonation in the sequential steps PCE to DCE are equally depleted and their δ 2 H will, on average, increase in the TCE to DCE step related to the inverse fractionation effect. It might be that the fractionation factor related to protonation is in fact different and possibly smaller than assumed for the PCE to TCE step. In that case, δ 2 H-DCE and δ 2 H-TCE will be more similar along the flow path. Note that δ 2 H-TCE equals δ 2 H-DCE in absence of any hydrogen SKIE in SRD (Fig. 5: middle: thin black lines). The simulated δ 2 H-DCE values agree well with the δ 2 H-cDCE values observed at a field site (− 211 ± 20‰, n = 10, one outlier excluded) where cDCE was the reaction product of PCE RD (Audi-Miro et al., 2015). We used the model to evaluate the impact of the ϵ NRP(MEAN) value of TCE (+ 34‰, solid lines; 0‰, thin black lines; − 34‰, dashed lines) on the predicted δ 2 H values of cDCE, VC, and ETH. In the case of a normal isotope fractionation effect in TCE RD (Fig. 5: left and middle: dashed lines), δ 2 H of cDCE/VC/ETH slightly increase instead of decrease downgradient. Initially produced TCE daughter products are more depleted but become enriched as also TCE becomes more enriched over distance/time (note result not shown in Fig. 5: left). Interestingly, the final δ 2 H values of intermediate daughter products cDCE and VC become more enriched compared to the case of an inverse isotope effect in TCE RD (Fig. 5: left and middle: dashed versus solid lines), whereas final δ 2 H-ETH remains identical irrespective of the SKIE values of the preceding steps because of reasons of isotope mass balance. Thus, provided source TCE is strongly enriched, H isotope analysis could be useful to distinct among source TCE and TCE produced through PCE reductive dechlorination. Besides δ 2 H-TCE, also δ 2 H of lower chlorinated daughter products and ETH could be informative about their source compound (PCE or TCE) as their δ 2 H is strongly different between the two scenarios (about 200‰ or more, which is ten times or more the uncertainty of H-CSIA (±20‰)).
Finally, Fig. 5 (right) shows the simulation results of a mixed PCE/ TCE source (1:1 molar ratio). PCE reductive dechlorination produces strongly depleted δ 2 H-TCE which mixes with the pool of strongly enriched source TCE. As a result δ 2 H-TCE values decrease rapidly away from the source and become intermediate of the pure PCE and TCE source scenario values. Likewise, δ 2 H values intermediate of these pure source scenarios are simulated for the degradation products. Note that in the case of pure sources (PCE or TCE) and SRD as sole reaction pathway, the δ 2 H values of the chlorinated ethenes and ETH are relatively constant in the flow direction. However, in the case of a mixed PCE/TCE source, strong decreases in δ 2 H are simulated especially near the source area ( Fig. 5: right). Summarizing, the model scenarios indicate that PCE daughter products are considerably more depleted than those produced from a pure (high δ 2 H) TCE source. Predicted δ 2 H values of specific daughter products remain fairly constant along flow suggesting the potential of H-CSIA for source apportionment.

Conclusions and outlook
The developed numerical model serves as a template model to interpret C, H, and Cl CSIA data in SRD of halogenated hydrocarbons in Fig. 5. Results of model scenario transport simulations on hydrogen isotope fractionation following first-order sequential reductive dechlorination of chlorinated ethenes: source is TCE (left), source is PCE (middle), source is PCE/ TCE (50%/50%) (right). Fractionation effects were as in simulation S4a: ϵ NRP of TCE = + 34‰ (solid lines), assuming instead a normal isotope effect of TCE RD of equal magnitude: ϵ NRP of TCE = − 34‰ (dashed lines), and assuming all ϵ NRP(MEAN) values are equal to zero (thin black lines). general, with the aim of investigating (S)KIEs, intramolecular halogen isotope ratio heterogeneity, and protonation effects. Furthermore, this model has great promise in application to CSIA-based MNA of sites polluted with chlorinated solvents to demonstrate and clarify the mechanisms of contaminant destruction. Extending the model to include the PCE to TCE step in SRD and alternative one-step degradation pathways such as chemical reduction and biological oxidation is straightforward. 3-D simulations using this PHREEQC model with either PHAST (Parkhurst et al., 2010) or PHT3D (Prommer and Post, 2010) are possible (Kuder et al., 2014) but would require long calculation times. Alternatively, the recently developed analytical BIO-CHLOR-ISO model (Höhener, 2016) enables rapid 3-D simulations of concentrations and carbon and chlorine isotope ratios but has limitations: it cannot cope with heterogeneous conditions, Monod kinetics, multiple DCE isomers, and does not consider the possibility of intramolecular heterogeneity of δ 37 Cl in TCE.
The main limitations of model application are besides uncertainty on intramolecular heterogeneity of TCE, probably the current limited sets of fractionation factors, particularly those of hydrogen, which are still uncertain under field conditions for SRD and completely unknown for oxidation. The model may, however, obtain such fractionation factors via model calibration provided the level of field site complexity is low and data coverage is high. Finally, the modeling of H-CSIA data may improve source apportionment of daughter products deriving from TCE or PCE since those derived of PCE should be considerably be more depleted, provided that the source TCE is strongly enriched as reported for the majority of modern TCE products Shouakar-Stash et al., 2003).