How long do natural waters “remember” release incidents of Marcellus Shale waters: a first order approximation using reactive transport modeling

Natural gas production from the Marcellus Shale formation has significantly changed energy landscape in recent years. Accidental release, including spills, leakage, and seepage of the Marcellus Shale flow back and produced waters can impose risks on natural water resources. With many competing processes during the reactive transport of chemical species, it is not clear what processes are dominant and govern the impacts of accidental release of Marcellus Shale waters (MSW) into natural waters. Here we carry out numerical experiments to explore this largely unexploited aspect using cations from MSW as tracers with a focus on abiotic interactions between cations released from MSW and natural water systems. Reactive transport models were set up using characteristics of natural water systems (aquifers and rivers) in Bradford County, Pennsylvania. Results show that in clay-rich sandstone aquifers, ion exchange plays a key role in determining the maximum concentration and the time scale of released cations in receiving natural waters. In contrast, mineral dissolution and precipitation play a relatively minor role. The relative time scales of recovery τrr, a dimensionless number defined as the ratio of the time needed to return to background concentrations over the residence time of natural waters, vary between 5 and 10 for Na, Ca, and Mg, and between 10 and 20 for Sr and Ba. In rivers and sand and gravel aquifers with negligible clay, τrr values are close to 1 because cations are flushed out at approximately one residence time. These values can be used as first order estimates of time scales of released MSW in natural water systems. This work emphasizes the importance of clay content and suggests that it is more likely to detect contamination in clay-rich geological formations. This work highlights the use of reactive transport modeling in understanding natural attenuation, guiding monitoring, and predicting impacts of contamination for risk assessment.


Background
The development of unconventional natural gas in the Marcellus Shale formation has grown rapidly in recent years. Significant concerns arise in parallel due to their possible impacts on water resources. Here Marcellus Shale waters (MSW) are defined as waters from gas wells including both flowback and produced waters. Marcellus Shale waters are typically characterized by high total dissolved solids (TDS, usually >200,000.00 mg/L), elevated concentrations of Br, Cl, major cations (Na, Ca, Mg, K), as well as Ba and Sr, often accompanied by natural occurring radioactive materials [1][2][3][4]. Accidental release of MSWs has been reported to occur through impoundments, drilling site discharge, spills, among others [5][6][7]. Although these major ions are of less environmental concern than toxic metals, their high concentrations can still pose adverse effects on human health. For example, Br may produce bromate through ozonation, a human carcinogen [8]. High Ba concentration can cause muscle weakness and affects blood pressure, nervous and circulatory system [9,10]. Their release can deteriorate water quality and aquatic ecological systems [6]. In 2013, four northeastern amphibian species have been recorded to be adversely affected by 50-1000 mg/L chloride, suggesting small accidental releases can impede breeding habitats [11]. They are also important indicators of fracking fluid, flowback and produced water, and brine contamination in aquifers or rivers [5,12,13].
Recent evidence highlighted the risk of MSW leakage into natural waters. Direct discharge of MSW into surface waters has been frequently reported [14,15]. In Pennsylvania, a total of 229 spills occurred from 2005 to 2015 [16], as illustrated in Fig. 1a. High concentrations of methane, saline brine [17,18] and 2-n-butoxyethanol (often used in the fracking fluids) [19] were found in drinking groundwater aquifers in Pennsylvania, indicating potential leakage associated with Marcellus Shale gas development. The discharge of MSW has been found to increase downstream Br and Cl concentrations by more than three orders of magnitude [14,15]. Ferrar et al. [14] found Ba and Sr surpassed the US maximum concentration level (MCL) after a deliberate MSW discharge. Sang et al. [20] reported 32-36% of heavy metals associated with colloids mobilized by flowback water flush.
These studies raise questions regarding the impacts of release incidents. How long and to what extent do natural waters (rivers and aquifers) "remember" the release of MSW? In other words, how long do MSW stay in natural waters? The ultimate transport and fate of released chemical species can be affected by many processes (Fig. 1b). Mixing of different waters occur immediately upon release, which means that the relative magnitude of water release rate and the background flow rate in the receiving waters can play a significant role in determining their concentrations [21]. Cations in the Marcellus Shale waters can participate in multiple water-rock interactions, including mineral dissolution and precipitation, ion exchange, and surface complexation when clay minerals are abundant. The geochemical conditions of receiving aquifers, therefore, can be important in determining dominant reactions, natural attenuation potential, and impacts of accidental releases [22]. There has been a significant lack of key measures that quantify and predict reactive transport and fate of chemical species from MSW.
The objective of this study is to (1) understand key processes that govern the reactive transport and fate of major cations from MSW; and to (2) to quantify time scales and magnitude of the release impacts on water quality under various release and receiving water conditions. It is important to note that here we focus on abiotic interactions, instead of microbe-mediated biodegradation of organic contaminants. Heavy metals are not included in this study as they deserve a separate study. The insights learned here can facilitate fundamental understanding of natural attenuation and assess environmental risks. Simulations were done under conditions relevant to natural waters in Bradford County in the Pennsylvania, where local residential concerns on water quality arise in parallel with the large number of drilled wells [23]. We use the multicomponent reactive transport model CrunchFlow [24], which solves conservation equations with respect to mass, momentum, and energy. It has been extensively used to understand and predict reactive transport of contaminants, and waterrock interaction in porous media [25,26]. To the best of our knowledge, this work is among the early studies that use reactive transport modeling tools to understand the impacts of Marcellus Shale waters in natural water systems.

Problem setup
As shown in Fig. 1b, MSWs are introduced into homogeneous and isotropic natural water systems including ground water in sandstone (S) aquifers and sand and gravel (SG) aquifers and surface water. This represents a base case scenario with major focus on the coupling of transport and geochemical reactions without considering spatial heterogeneities. The interactions between chemical species in MSWs and sediment (typically <2 vol%) in rivers are assumed negligible. The yellow numbers are the numbers of spills. b A schematic diagram of 1-dimensional modeling setup. We assume a release point where the Marcellus Shale waters are introduced into the surface water (river) or groundwater (aquifers). The release can occur through spills, discharge, leakage, seepage, among others The S aquifers and SG aquifers were chosen as representative aquifers because they dominate in Northeastern Pennsylvania [27]. They differ in mineralogical compositions, with the S aquifers containing much more clay. We chose a branch of the Susquehanna River to represent the river. The release characteristics of MSWs, including release rates, time duration, and therefore total volumes, can vary significantly. All these factors can influence the impacts of accidental release on natural water compositions.

Properties of natural waters and MSWs Natural water systems
We used the characteristics of a sandstone aquifer with dominant clay mineral of 21.7% in the Catskill Formation in Bradford County, PA. The S aquifer has a groundwater velocity of 0.20 m/day and is predominantly a low-rank graywacke with major minerals being quartz, mica (represented by muscovite) and other clays, and trace amount of carbonate (mostly calcite) [28]. In contrast, the Sand and Gravel aquifer has a groundwater velocity of 0.40 m/day and a lower clay amount than that of the S aquifer [29,30]. For rivers we choose conditions relevant to the Susquehanna River segment in Bradford County, PA [31], considering 2% (v/v) of suspended sediments [32]. The major difference between the surface and subsurface water systems are the orders of magnitude higher flow rates and the negligible presence of solid phases compared to the aquifers (Table 1).

Water composition
The three natural waters differ in their chemical composition [38] ( Table 2). The surface water has higher concentrations of sulfate and cations including iron, potassium, and zinc, while the ground waters are richer in calcium, magnesium, and sodium. The major difference between the surface and subsurface water systems are the orders of magnitude higher flow rates and the negligible presence of solid phases compared to aquifers. The MSW composition was chosen to be in the low concentration level of produced and flowback waters.

Characteristics of Marcellus Shale water release incident
A total of 9179 unconventional wells were installed in the Marcellus Shale formation in Pennsylvania from 2005 to 2015 [16]. A total of 229 spill accidents have occurred, dictating a spill possibility of 2.40% per well in average. The spill volumes varied from 0.003 to about 11.35 m 3 with the median value being 0.144 m 3 [47]. With the same spill volume, a release can occur at small rates for a long duration or high rates for a short time frame. The MSWs reached groundwater by seeping into groundwater aquifers, which is a relatively slow process. Here we assume a net water volume of 0.144 m 3 reaching natural waters; the actual spill water can be much larger as the vadose zone tends to trap a large percent of spilled water [47]. Here we do not explicit consider vadose zone processes.

Table 1 Mineral composition and flow velocity in the natural waters
a Four secondary minerals, including gypsum, celestite, barite, and gibbsite, are initially assigned with a volume fraction of 10 −10 for precipitation in simulated natural water domain [29,33] b Porosity and flow velocity are within the typical range for S aquifers in this area [34,35] c Porosity and flow velocity are within the typical range for SG aquifers [36,37]   The spill rates are varied to examine the importance of release characteristics. We define the dilution factor (DF): where Q MSW and Q NW are the volumetric flow rates (m 3 /s) of MSW and the receiving natural waters, respectively. The Q NW values are calculated as the product of flow velocity (m/day) and cross-sectional area of 1 m 2 in the model. As such, we focus on understanding processes at the immediate vicinity of the leakage point and flow path. The DF quantifies the extent of dilution upon release into natural waters. A high DF value means that the released MSW is quickly diluted by the fast background natural waters. It is important to note here that fluid injection into an aquifer typically only causes limited mixing at the fringes. Here by assuming well-mixed intruding fluid and background water at the injection point, we can use this as a rough estimation of the relatively magnitude of the injection fluid rate versus the background fluid rate.

Reactive transport modeling
Upon accidental release into natural water systems, the chemical species in the MSWs interact with natural waters and solid phases. Major processes include mixing, transport, and various types of water-rock interactions.

Reactive transport equations
Reactive transport models (RTM) have been extensively used to understand complex interactions among physical, chemical, and biological processes in porous media [48][49][50][51]. The governing mass conservation equation for a chemical component i that participates in ion exchange reactions can be written as follows: Here φ is porosity, C i is total concentration (mol/m 3 pore volume) of i, t is time (s), D i is diffusion/dispersion tensor (m 2 /s), u is flow velocity (m/s), N r is total number of kinetic reactions that involve species i, v ir is stoichiometric coefficient of species i associated with reaction r, R r is the rate of chemical reactions in which the species i is involved (mol/m 3 /s). The diffusion/dispersion coefficients and flow velocities are set constant with a dispersivity of 1.0 cm [52]. Here kinetic reactions include mineral dissolution and precipitation. Ion exchange and aqueous complexation are considered as fast and are equilibrium-controlled. This equation implies that v ir R r mass change rate of species i depends on diffusion/dispersion represented by the first term in the right hand side (rhs), advection described by the second term in the rhs, and reaction described by the third term. The term ρ ∂S ∂t represents mass exchange with solid phase through ion exchange, with ρ being solid bulk density (g/m 3 pore volume), and S being solid phase concentration of i (mol/g). This term is essentially a storage term taking into account mass accumulation of i on the solid phase [53]. The geochemical system here includes 18 chemical components ( Table 2) and 14 kinetic mineral reactions ( Table 3).
The RTM was implemented within a 10 m one-dimensional domain with 100 grid cells and a fixed spatial discretization of 0.1 m. The spatial discretization was chosen as the lowest one that results in the same modeling output as those from spatial resolutions higher than 0.1 m. The injection point is the first grid cell. As such, we are simulating the first 10 m immediately down gradient of an injection point. We choose not to do numerical experiments in a large spatial domain of kilometers because the goal here is to understand dominant geochemical processes that govern natural attenuation of Marcellus Shale waters. A domain of 10 m is sufficient for such purpose. As will be discussed later, the dimensionless time derived from this work is not confined to the physical length of simulated domain. Running simulation at large spatial scales however presents additional challenges, largely because reaction parameters in literature are mostly measured in relatively small scale laboratory systems at the spatial scale of 10 0 -10 2 cm. Reaction parameters, in particular reaction kinetic constants and effective surface areas, are often orders of magnitude lower in large scale heterogeneous systems [49,55,71,72]. Running simulations at the scale of kilometers therefore requires overcoming upscaling of reaction processes, which has been a long-standing and unresolved puzzle [73,74].
We examined five cases with different types of water systems and release characteristics ( Table 4). The three release rates were determined by using reported dilution factors in literature [47,75,76] and Eq. (1).

Table 3 Reaction network, Reaction thermodynamics, and kinetics for mineral-water interactions
a SSA values are from the laboratory studies in the literature which are generally observed to be faster than those from the fields [55,70]  NaX ⇔ Na + + X − 5.0 × 10 −5 eq/g 3.0 × 10 −5 eq/g 0.00

Mineral dissolution and precipitation
Mineral reactions are listed in Table 3 with their equilibrium constants and reaction kinetics. In the systems in this paper, most waters are at close to neutral conditions so we only use rate laws based on neutral mechanisms and follow the classical transition-state-theory-based (TST) rate law [77]: Here R Ca is the rate of calcite dissolution (mol/s), A is the reactive surface area (m 2 ). The ion activity product (IAP) is defined as a Ca 2+ a CO 2− 3 , and K eq is the equilibrium constant. The IAP/K eq measures the distance from equilibrium. If IAP is lower than K eq , the water is under saturated and calcite dissolves; if IAP is higher than K eq , the system is over saturated and calcite precipitates. The equilibrium constants are from the standard EQ 3/6 geochemical database [78].

Ion exchange
Ion exchange is represented as follows [79,80]: Here (aq) and (s) refer to aqueous and exchanged phases, respectively; X − denotes negatively charged exchange sites occupied by cations A u+ and B v+ of charge u and v for A and B, respectively. Ion exchange reactions are commonly calculated through the Vanselow convention using cation mole fractions on the exchange sites [81]. The overall cation exchange capacity was calculated based on volume fraction and surface area of clay minerals including muscovite, illite, kaolinite, clinochlore-14A and daphnite-14A. The selectivity coefficients in Table 3 indicate cation affinity to solid surface. The species Ba and Sr have higher affinity than Ca and Mg, which in turn have higher affinity than Na and K. This means that under similar concentration conditions, Ba and Sr tend to be exchanged onto clay surface first before Ca, Mg, and K. The very high Na concentration in Marcellus Shale waters also induces the exchange of Na onto solid surface compared to Ca and Mg.

Quantification of release impacts
We define several terms to quantify release impacts on natural water composition. The maximum concentration in receiving waters during release, C max , quantifies the magnitude of the impacts. The residence time τ r is calculated by the domain length divided by the natural water flow velocity; it quantifies the time scale at which a non-reactive species stays in the domain of interest. The recovery time, τ recovery , is the time scale for each species to return to within 5% difference from its original concentration. Because different species involve different types of water-rock interactions (e.g., mineral precipitation versus ion exchange), this time scale can vary drastically among species. The relative recovery time τ rr is defined as the ratio of τ recovery over τ r . The τ rr quantifies the time duration (relative to residence times) that the released chemical species still remain in the simulation domain. All these terms are calculated based on modeling observations from output of numerical experiments. Note that τ rr is a dimensionless number and its value is not constrained to the length or time scale of the calculation domain. The τ rr is similar to the concept of effective retardation coefficient and is species specific [53]. For instance, the retardation factors of Ba and Sr are 111 and 60, respectively under neutral condition [82]. The cations generally follow the retardation sequence of Mg<Ca<Sr<Ba [83,84]. A higher affinity to solid surface leads to a larger retardation and therefore a slower movement, longer residence time and ultimately longer τ rr and memory.

Results and discussion
"Controlling processes in the S aquifer" section focuses on understanding processes that control transport and fate of major species in the S aquifer. "Effect of release characteristics in the S aquifer" section assesses the role of release rates. "Effect of receiving water bodies" section compares reactive transport of major species under different release rates under different natural water conditions.

Controlling processes in the S aquifer
Here we examine the spatio-temporal evolution of major species after release into the S aquifer under four scenarios of increasing process complexity: a case including only mixing without any reactions ("MIX"), a case including mixing and mineral dissolution/precipitation ("MIX + DISS/PPT"), a case with mixing and ion exchange without mineral dissolution/precipitation ("MIX + IEX"), and a case including mixing, mineral dissolution/precipitation, and ion exchange ("MIX + DISS/ PPT + IEX"). The release occurred from day 10 to day 25 at the rate of 1.11 × 10 −7 m 3 /s in all cases. Before the release accident, initial water-rock equilibria are established by continuously injecting natural waters into the simulated domain until their compositions are stabilized.

Temporal evolution at the release point
Br and Cl The breakthrough of Br and Cl in the four scenarios are the same due to their non-reactive nature (Fig. 2) Fig. 3b and Na increase on the solid phase in Fig. 3f. This indicates that the quick increase is caused by the ion exchange between Na and the presorbed Ca. This desorbed Ca leads to calcite precipitation with sharply increasing rates during release (Fig. 3g positive calcite rates), which decreases aqueous Ca significantly and cause calcite dissolution afterwards (Fig. 3g negative calcite rates). At the time when release stops, the precipitation even draws Ca concentration to below background concentration.
The system eventually relaxes back to the original state. Despite the differences in MIX + IEX and MIX + DISS/ PPT + IEX cases, similar [Ca] in these two cases indi- in the MIX + DISS/PPT + IEX case also leads to much higher calcite precipitation rate during release (Fig. 3a, g). Similar to Ca, Mg also participates in mineral dissolution/precipitation (clinochlore-14A and dolomite) and ion exchange (Table 3). Compared to Ca, its concentration is about an order of magnitude lower in both background and MSW. Its evolution at the release point resembles that of Ca (Fig. 3c). Although not shown here, dolomite is close to equilibrium while the dissolution rate of clinochlore-14A is in the order of 10 −10 mol/s. Comparison between the 4 cases show that Mg behaves similarly to Ca and is primarily controlled by ion exchange. The highly elevated Na in MSW leads to massive exchange on the solid surface. After the release stops, Na slowly desorbs, resulting in a long tail for over more than 150 days (Fig. 3e, f ). Conversely, Ca and Mg sorb back to the solid (Fig. 3b, d), which results in lower aqueous Ca and Mg concentration when compared to the background concentration and calcite dissolution, as indicated by the negative calcite rates in the MIX + DISS/ PPT + IEX case after the release. They eventually return to background concentrations after continuous groundwater flushing and reach equilibrium again.
The original pH values are 7.40 and 6.90 in the S aquifer and MSW, respectively. Values of pH drop upon release in all cases (Fig. 3h). In the MIX and MIX + IEX cases, pH drops slightly and returns immediately to its background when release stops. In the other two cases that involve mineral dissolution and precipitation, pH drops much more significantly during the release, primarily due to calcite precipitation. In the MIX + DISS/PPT + IEX case, because the ion exchange kicks out sorbed Ca and increased aqueous [Ca], the higher calcite precipitation rates lead to more significant pH decrease (Fig. 3g). The calcite dissolution leads to pH increase for an extended period of time until Ca dominates the solid surface again. In general, the pH curve mirrors the shape of calcite rate. The pH values relax back to its background immediately after the release in all cases except the MIX + DISS/ PPT + IEX case where pH is mostly controlled by calcite dissolution and precipitation reactions.  (Fig. 4). After the release stops, Ba and Sr slowly desorb over a longer period of time.

Ba and Sr
Although not shown here, barite and celestite precipitate in negligible rates, indicating the dominant role of ion exchange.

Spatio-temporal evolution in the MIX + DISS/PPT + IEX case
Here we examine the spatio-temporal evolution of major species in the MIX + DISS/PPT + IEX case where all processes are included. The release occurs between day 10 and 25 at the rate 1.11 × 10 −7 m 3 /s. Tracers During release, the [Br] amd [Cl] in the down gradient rapidly increase (Fig. 5). After the release, Cl returns to background concentration starting from the release point. The high concentration zone gradually  Reactive species The five major cations can be categorized into two groups (Fig. 6). Group I includes Ca and Mg (top two rows), both of which are in the MSW and are originally on exchange sites. They are mobilized through ion exchange with cations in the MSW, primarily Na, Ba, and Sr. During release, their aqueous concentration peaks in some zone while the corresponding solid concentration show "valley" of low concentrations. The peaks expand over time during the release. At the end of the incident, their aqueous concentrations are lower than the background concentrations due to their exchange back to the surface. Correspondingly, their solid phase low concentration valleys migrate down gradient slowly over a much longer time scale, long after the release stops on day 25. The depletion zone also becomes wider and shallower due to dispersion as they migrate out of the system.
The Group II species consist of Na, Ba, and Sr, which are abundant in the MSW and exchange with solid surface upon release, displacing Ca and Mg. During release they all show highest aqueous and solid concentrations at the release point, while quickly decrease down gradient. Both peak aqueous and solid concentrations increase over time during release. After the release stops, these cations on the exchange sites gradually become remobilized back into the aqueous phase through ion exchange. Compared to Group I species, Group II species show peaks in both aqueous and solid phases that migrate at similar rates down gradient. The concentration peaks become wider and shallower over time.

Quantification of memory indexes from spatio-temporal profiles
The "maximum concentration" C max and the "recovery time" τ recovery can be calculated from the spatial profiles discussed above (Figs. 5, 6). These two measures differ significantly from one species to another. The C max of tracers (Br and Cl) are controlled by the mixing process. After release the system returns to their background after approximately one residence time. For Group I species (Ca and Mg), C max values are higher than those estimated by their dilution factor because they are mobilized from the solid surface during release. For Group II species (Na, Ba, and Sr), their peak concentrations equal to or are lower than those predicted by dilution factor because they exchange onto solid surface. The memory or the time scales of the reactive species are dictated by their affinity to the surface. On day 160, the peak for Na has disappeared, indicating its migration out of the system. In contrast, the peaks of Ba and Sr are approximately at 8 m at that time. As indicated in the ion exchange coefficients in Table 3, the affinity to the surface is Ba/Sr>Ca/Mg>Na. The Ba and Sr therefore migrate out of the system much slower. The τ rr values are 6. 79, 9.25, 9.38, 20.09, 18.76 for Na, Ca, Mg, Ba, Sr, respectively. This means that it takes 6.79 residence times to flush out Na, 9.25/9.38 residence time for Ca/Mg, and 20.09/18.76 for Ba/Sr, which are consistent with their affinity to the solid surface. This gradient of time scales consistent with their gradient of the affinity to the solid surface is similar to the chromatographic effects in literature [53].

Effect of release characteristics in the S aquifer
Three cases were compared with the same release volume of 0.144 m 3 however at different release rates. The "High" release rate is 1.11 × 10 −7 m 3 /s for 15 days, the  Ca (a, b), Mg (c, d), Na (e, f), Ba (g, h), and Sr (i, j) same as the case in the Sect. 4.1. The "Medium" rate is 5.55 × 10 −8 m 3 /s for 30 days. The "Low" release rate is 1.11 × 10 −8 m 3 /s for 150 days ( Table 4). The corresponding dilution factors are 21.85, 42.70, and 209.54, respectively. Figure 7 shows the spatio-temporal evolution for Br (tracer), Ca (Group I), and Na (Group II). In general, the higher release rate, the higher impact on the water chemistry. For the tracers, C max are essentially the MSW concentrations divided by the corresponding dilution factor in each case. For the reactive species, the low release rate leads to much lower aqueous and/or solid concentrations than in the high rate case. In addition, it takes shorter time to flush out Na in the low release rate case and therefore the system recovers sooner.

Effect of receiving water bodies
The river has the highest flow velocity (2.76 × 10 4 m/ day) compared to the S aquifer (0.20 m/day) and SG aquifers (0.40 m/day). The S aquifer has 21.7 vol% of clay content compared to 0.9% in the SG aquifers and zero in the river. The release occurred at the same high rate of 1.11 × 10 −7 m 3 /s for 15 days. The dilution factors for the three receiving natural waters are 21.85, 42.70, and 2.87 × 10 6 , for S aquifer, SG aquifers, and river, respectively. Figure 8 shows the effects of receiving water characteristics on the reactive transport of major species. With orders of magnitude higher flow velocity, the river has no memory of MSW-all concentrations are at the background concentration. The MSW however leaves their footprint on the ground water aquifers. The [Br] during the release is lower in the SG aquifers than in the S aquifer due to the higher flow velocity in the SG aquifers. Note that the background [Br] in the two aquifers are also different, with lower [Br] in the SG aquifers. The reactive species behave similarly to the tracers in the SG aquifers because of the low clay content and the lack of ion exchange. The higher dilution factor in the SG aquifers lead to a concentration about half of the maximum [Na] in the S aquifer at the release point, while in the down gradient [Na] is higher in the SG aquifers because negligible ion exchange occurs compared to that in the S aquifer. The [Na] and [Ca] return to the background concentration much faster in the SG aquifers than in the S aquifer.

Impacts of the release incidents
Values of C max and τ recovery quantify the impacts and time scales of release accidents. The numerical experiments indicate that C max of Ca, Na and Cl are 2 orders of magnitude higher than Ba, Sr and Mg (Fig. 9). The river has the lowest C max compared to the groundwater aquifers with all C max values below the drinking water standard [85][86][87]. In the SG aquifers, all species behave as if they are non-reactive with their C max values proportional to their concentrations in the original MSW. Only the C max of Cl and Na exceed the drinking water standard. In the S aquifer, however, almost all species exceed drinking water standards in the High and Medium release rate cases. In the Low release rate case, only Br and Ba exceed the drinking water standard. The τ recovery values vary by orders of magnitude and depend on specific characteristics of natural waters, release incidents, and individual The "during release" curves are on day 10 after the release starts. The "after release" curves are on day 5 after release stops species (Fig. 9). The S aquifer remembers the incident longer compared to SG aquifers due to the lower flow velocity and higher clay content. The Low release rate case to recover much fast back to the background concentration than the Medium and High cases. Their corresponding relative time scales, τ rr , however, vary only from 1.0 to a maximum of about 20 (Fig. 9b). In fact, under all conditions where chemical concentrations are controlled by the mixing process, values of τ rr are close to 1. This includes non-reactive species in all natural waters at all release rates, and reactive species in natural waters with negligible clay content (rivers and SG aquifers). Only in S aquifer with abundant clay, τ rr values depend on cation affinity to solid surface with τ rr between 5 and 10 for Na, Ca, and Mg, and 15-20 for Sr and Ba.

Environmental implications
Despite the fact that multiple minerals are involved in dissolution and precipitation, these reactions play a relatively small role compared to ion exchange. This highlights the importance of clay content in determining the time scales and impact of incidental release on natural waters.
The results have interesting implications in understanding reactive transport, monitoring, and detection of contaminants from Marcellus Shale waters in natural water systems. In a controversial example involving unconventional gas wells in a sandstone formation near Pavillion, Wyoming, EPA detected contamination in shallow monitoring wells from 2010 to 2011 [88]. Synthetic organic compounds used in hydraulic fracturing fluids were detected in monitoring wells; [Cl] and [K] were found more than one order of magnitude higher in a monitoring well than the background concentrations. In a second time sampling in the same wells in April and May 2012, some previously detected compounds (e.g., xylenes, toluene) were not found and a number of other compounds have lower concentrations than the previous analysis. As shown in the spatiotemporal figures, there are only certain "time windows" that the signature of MSW can be observed in a specific location, which indicates the ephemeral nature of contamination events and the transient and elusive contamination plumes. This imposes significant challenges to monitoring and detection of groundwater contamination [13].
There have been several cases that discharged MSW were detected in rivers. For example, at the discharge point, [Cl] and [Br] were 6000-and 12,000-fold higher than that in the stream background, both exceeding the drinking water standard [15]. This is a case where the MSW was discharged into river with large volume and therefore the dilution factor of 739 is more than three orders of magnitude lower than that in our model (2.87 × 10 6 ). During dry season, low flow rates in rivers lead to lower DF [7], which also increase the possibility of contamination detection. Table 5 shows a few cases where elevated chloride concentrations were reported when MSWs were discharged into river. The DF values in these cases, estimated as the ratio between the flow rates of the river and the discharge rate, vary between 510 and 1246. These values are 3-4 orders of magnitude lower than the DF value in the incidental release case in this work. Fig. 8 Profiles of Br, Ca, Ca on solid surface, Na, Na on solid surface during release (left column) and after release (right column) in the sandstone aquifer, sand and gravel aquifer, and river, respectively. The release rate is 1.11 × 10 −7 m 3 /s for 15 days. The "during release" is on day 10 after the release starts. The "after release" is on day 5 after release stops

Limitations
This study is for the specific hydrological and geochemical conditions in Northeastern Pennsylvania in homogeneous systems of one-dimensional 10 m immediately down gradient of the release point where the impacts on natural waters are most significant. This is different from three dimensional natural water systems in reality that have larger dispersion and spreading. As such, the calculation here likely overestimates C max and τ recovery and therefore represents the worst case scenario. However, the quantitative term defined here, especially the relative recovery time τ rr , is dimensionless and is not restricted to the length scale of the simulation domain. For example, if the estimated τ rr for a particular species is 5.0, it means that the time needed for recovery is five times of the water residence time. This estimation can be used for systems of different lengths and of flow velocities, because residence times scale with length and flow velocity. As such, τ rr provides the approximation needed for estimating memory or time scales of release incidents in natural waters. In addition, as long as geochemical conditions remain relatively similar, the dominant water-rock interactions are similar.
Here we mainly focus on water-rock interactions of major cations in the MSW without considering redox reactions and biodegradation of organic contaminants that can be present in MSWs. If organic contaminants are present and used by microbe as carbon source, Fig. 9 The memory index of natural waters: C max and a τ recovery and b τ rr of major species in the river (filled squares), SG aquifer (filled triangles), and S aquifer with high release (filled circles), medium release (crossed circles), and low release rates (open circles). Both are calculated from the modeling output of spatio-temporal concentration evolution. The C max is determined as the maximum aqueous concentration during release. The τ recovery is the time scale for each species to return to within 5% difference from its background concentrations in natural waters. The relative recovery time τ rr , calculated as the ratio of τ recovery over τ r , is a measure of the time scale that natural waters remember the incident relative to their residence time. Each species is represented by one color, with dashed line of the same color being their drinking water standard. In S aquifer with abundant clay, τ rr values depend on cation affinity to solid surface with τ rr between 5 and 10 for Na, Ca, and Mg, and 15-20 for Sr and Ba biodegradation reactions will transform organic contaminants into dissolved inorganic carbon, which can increase the concentrations of bicarbonate significantly. Under such circumstances, carbonate precipitation may play a much more significant role, as indicated in literature [49,67,90]. This study also considers homogeneous systems. Natural groundwater aquifers are typically layered with heterogeneous distribution of hydrological and geochemical properties [91,92]. Such spatial heterogeneities have long been reported to cause order-of-magnitude longer tail for non-reactive tracers [93][94][95] and lower reaction rates [25,96,97]. The specific characteristics of different natural water systems, including spatial distribution of clay lenses and layers, therefore, will play a significant role in determining the recovery time of natural waters from incidental release.
In addition, we did not consider the vadose zone processes. Vadose zone processes will affect how much spill volume and chemicals will get into aquifers. However, the major reactive transport processes in natural waters will remain the same and the time scales that the released chemicals remain in aquifer will still be determined by their affinity to the solid surface-this aspect is not going to change whether we include vadose zone processes or not.

Conclusions
Recent studies on MSWs have mostly focused on evidence linking altered water composition to possible release of Marcellus Shale waters. Process-based understanding and quantification on reactive transport of accidentally released chemicals, however, are largely lacking.
Here we use major cations as tracers of release events and reactive transport numerical experiments to illustrate key processes that determine the impacts of accidental release.
The magnitude of the impacts is quantified by C max , the maximum observed concentration during release, while the time scale of the impact, τ recovery , the required time duration to recover to within (100 ± 5%) of its back ground concentration. We also define a dimensionless number τ rr that is the relative ratio of the τ recovery compared to the residence time of natural waters τ r . Our results show that in rivers and SG aquifers with negligible clay content, mixing process controls C max and τ recovery of all species. The dilution factor determines C max while τ recovery approximates the residence time. In clay-rich natural water systems, ion exchange plays a dominant role compared to mineral dissolution and precipitation. The S aquifers with abundant clay selectively remember Sr and Ba for 10-20 residence times due to their higher affinity to clay surface compared to 5-10 residence times for Na, Ca, and Mg. This highlights the importance of clay content in both monitoring and natural attenuation of chemicals from Marcellus Shale waters. This suggests that under otherwise similar conditions, it is more likely to detect contamination in clay-rich geological formations because it takes longer for chemicals to return to its original state in these formations.
This work highlights the usefulness of reactive transport modeling in process understanding and in guiding sampling and monitoring in natural water systems. Findings from this work facilitates prediction of contaminant transport and fate, quantifies impacts of released MSWs in natural waters, and provides insights on risk assessment and strategies for sustainable shale gas development.