Skip to main content
  • Research Article
  • Open access
  • Published:

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.


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.

Fig. 1
figure 1

a The numbers of Marcellus Shale water release accidents in Pennsylvania from 2005 to June 8, 2015, with 78% of spills occurred in Northeastern PA. Red spot indicated the location of Bradford County. 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

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 water–rock 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 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).

Table 1 Mineral composition and flow velocity in the natural waters

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.

Table 2 Composition of natural waters and Marcellus Shale water (mg/L)

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 m3 with the median value being 0.144 m3 [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 m3 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. The spill rates are varied to examine the importance of release characteristics.

We define the dilution factor (DF):

$$ DF = \frac{{Q_{MSW} + Q_{NW} }}{{Q_{MSW} }} $$

where QMSW and QNW are the volumetric flow rates (m3/s) of MSW and the receiving natural waters, respectively. The QNW values are calculated as the product of flow velocity (m/day) and cross-sectional area of 1 m2 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,50,51,51]. The governing mass conservation equation for a chemical component i that participates in ion exchange reactions can be written as follows:

$$ \frac{{\partial (\phi C_{i} )}}{\partial t} + \rho \frac{\partial S}{\partial t} = \nabla \cdot \{ \phi {\mathbf{D}}_{{\mathbf{i}}} \nabla (C_{i} ) - \phi {\mathbf{u}}C_{i} \} + \sum\limits_{r = 1}^{Nr} {v_{ir} R_{r} } $$

Here \(\phi\) is porosity, Ci is total concentration (mol/m3 pore volume) of i, t is time (s), D i is diffusion/dispersion tensor (m2/s), u is flow velocity (m/s), Nr is total number of kinetic reactions that involve species i, vir is stoichiometric coefficient of species i associated with reaction r, Rr is the rate of chemical reactions in which the species i is involved (mol/m3/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 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 \( \rho \frac{\partial S}{\partial t} \) represents mass exchange with solid phase through ion exchange, with ρ being solid bulk density (g/m3 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).

Table 3 Reaction network, Reaction thermodynamics, and kinetics for mineral–water interactions

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 100–102 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 4 Simulation scenarios for Marcellus Shale water release

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]:

$$ R_{Ca} = kA\left( {1 - \frac{IAP}{{K_{eq} }}} \right) $$

Here RCa is the rate of calcite dissolution (mol/s), A is the reactive surface area (m2). The ion activity product (IAP) is defined as \( a_{{{\text{Ca}}^{2 + } }} a_{{{\text{CO}}_{3}^{{2{ - }}} }} \), and Keq is the equilibrium constant. The IAP/Keq measures the distance from equilibrium. If IAP is lower than Keq, the water is under saturated and calcite dissolves; if IAP is higher than Keq, 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]:

$$ uBX_{v} \left( s \right) + vA^{u + } \left( {aq} \right) \Leftrightarrow vAX_{u} \left( s \right) + uB^{v + } \left( {aq} \right) $$

Here (aq) and (s) refer to aqueous and exchanged phases, respectively; X denotes negatively charged exchange sites occupied by cations Au+ and Bv+ 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, Cmax, 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 m3/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). The concentrations increase upon release and return to background concentration when the release stops. Their concentrations in the MSW are 185.00 and 29,252.00 mg/L, respectively. With the dilution factor of 21.85, the calculated Br and Cl concentrations during release are 8.82 and 1404.00 mg/L, respectively, approximating their MSW concentrations divided by the dilution factor plus the background concentration.

Fig. 2
figure 2

Evolution at the release point for Br under four scenarios. All four color lines overlap. The grey shaded zone represents the release period. Due to its non-reactive nature, the inclusion of different processes does not affect their evolution

Na, Ca, and Mg

The Na concentration ([Na]) is the highest (13,200.00 mg/L) among the three species in the MSW. The Na exchanges with presorbed Ca and Mg at the solid concentrations of 2.72 × 10−7 and 2.55 × 10−8 mol/g, respectively. The Ca therefore depends on the mixing with the ground water, dissolution and precipitation of calcite, and ion exchange. In the MIX case, Ca behaves similarly to Cl. The [Ca] in the MIX + DISS/PPT case is lower than that in the MIX case because of calcite precipitation, as indicated by the positive calcite rate in Fig. 3g. In the MIX + IEX case, the [Ca] increases sharply upon release, which echoes the fast Ca decrease on the surface in 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 indicate the dominance of ion exchange and relatively minor role of calcite dissolution/precipitation when both processes coexist. Compared to the MIX + DISS/PPT case, the increase in [Ca] in the MIX + DISS/PPT + IEX case also leads to much higher calcite precipitation rate during release (Fig. 3a, g).

Fig. 3
figure 3

Evolution at the release point for a Ca (mg/L) in logarithmic scale, b Ca on exchange sites (mol/g solid), c Mg (mg/L) in logarithmic scale, d Mg on exchange sites (mol/g solid), e Na (mg/L) in logarithmic scale, f Na on exchange sites (mol/g solid), g calcite reaction rate (mol/m2/s) (negative indicates dissolution and positive values indicate precipitation), and h pH. Grey line overlaps with the black line

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.

Ba and Sr

Barium and strontium exchange with presorbed cations Ca and Mg, leading to decreased aqueous [Ba] and [Sr], and increased aqueous [Ca] and [Mg] in MIX + DISS/PPT + IEX (Fig. 4). After the release stops, Ba and Sr slowly desorb over a longer period of time. Although not shown here, barite and celestite precipitate in negligible rates, indicating the dominant role of ion exchange.

Fig. 4
figure 4

Evolution at the release point for a Ba in water (mg/L), b Ba on surface (mol/g solid), c Sr (mg/L), d Sr on surface (mol/g solid). Ion exchange controls concentrations of these species while mineral dissolution and precipitation play a minor role

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 m3/s.


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 migrates out of the domain until the system returns to its background.

Fig. 5
figure 5

Spatio-temporal evolution of Br concentration in the sandstone aquifer in the MIX + DISS/PPT + IEX case on days 11, 25, 27 and 160. Release starts on day 10 and ends on day 25. The other tracer Cl behaves the same as Br

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.

Fig. 6
figure 6

Spatio-temporal profiles of major species in the sandstone aquifer under the MIX + DISS/PPT + IEX scenario on days 11, 25, 27 and 160. Left column is for aqueous concentrations (mg/L); right column is for concentrations on solid surface (mol/g solid). Rows from the top to bottom: Ca (a, b), Mg (c, d), Na (e, f), Ba (g, h), and Sr (i, j)

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” Cmax 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 Cmax 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), Cmax 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 m3 however at different release rates. The “High” release rate is 1.11 × 10−7 m3/s for 15 days, the same as the case in the Sect. 4.1. The “Medium” rate is 5.55 × 10−8 m3/s for 30 days. The “Low” release rate is 1.11 × 10−8 m3/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, Cmax 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.

Fig. 7
figure 7

Profiles of Br, Ca, Ca on solid surface, Na, Na on solid surface in the sandstone aquifer during release (left column) and after release (right column) under the three release cases. The High, Medium, and Low release rates are 1.11 × 10−7 m3/s for 15 days, 5.55 × 10−8 m3/s for 30 days, and 1.11 × 10−8 m3/s for 150 days, respectively. The “during release” curves are on day 10 after the release starts. The “after release” curves are on day 5 after release stops

Effect of receiving water bodies

The river has the highest flow velocity (2.76 × 104 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 m3/s for 15 days. The dilution factors for the three receiving natural waters are 21.85, 42.70, and 2.87 × 106, 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.

Fig. 8
figure 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 m3/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

Impacts of the release incidents

Values of Cmax and τrecovery quantify the impacts and time scales of release accidents. The numerical experiments indicate that Cmax of Ca, Na and Cl are 2 orders of magnitude higher than Ba, Sr and Mg (Fig. 9). The river has the lowest Cmax compared to the groundwater aquifers with all Cmax values below the drinking water standard [85,86,87]. In the SG aquifers, all species behave as if they are non-reactive with their Cmax values proportional to their concentrations in the original MSW. Only the Cmax 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 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.

Fig. 9
figure 9

The memory index of natural waters: Cmax 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 Cmax 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


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 spatio-temporal 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 × 106). 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.

Table 5 Cases with contamination detected during direct discharge of MSW into rivers [14, 15, 89]


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 Cmax 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, 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.


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 Cmax, 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 Cmax and τrecovery of all species. The dilution factor determines Cmax 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.


  1. Chapman EC, Capo RC, Stewart BW, Kirby CS, Hammack RW, Schroeder KT et al (2012) geochemical and strontium isotope characterization of produced waters from Marcellus Shale natural gas extraction. Environ Sci Technol 46(6):3545–3553. doi:10.1021/Es204005g

    Article  Google Scholar 

  2. Acharya HR, Henderson C, Matis H, Kommepalli H, Moore B, Wang H (2011) Cost effective recovery of low-TDS frac flowback water for re-use. GE Global Research, Niskayuna. Report No.: Contract No.: DE-FE0000784

  3. Olmstead SM, Muehlenbachs LA, Shih JS, Chu ZY, Krupnick AJ (2013) Shale gas development impacts on surface water quality in Pennsylvania. P Natl Acad Sci USA 110(13):4962–4967. doi:10.1073/pnas.1213871110

    Article  Google Scholar 

  4. Haluszczak LO, Rose AW, Kump LR (2013) Geochemical evaluation of flowback brine from Marcellus gas wells in Pennsylvania, USA. Appl Geochem 28:55–61

    Article  Google Scholar 

  5. Vidic RD, Brantley SL, Vandenbossche JM, Yoxtheimer D, Abad JD (2013) Impact of Shale gas development on regional water quality. Science. 340(6134):1235009. doi:10.1126/science.1235009

    Article  Google Scholar 

  6. Myers T (2012) Potential contaminant pathways from hydraulically fractured Shale to aquifers. Ground Water. doi:10.1111/j.1745-6584.2012.00933.x

    Google Scholar 

  7. Rozell DJ, Reaven SJ (2012) Water pollution risk associated with natural gas extraction from the Marcellus Shale. Risk Anal 32(8):1382–1393

    Article  Google Scholar 

  8. Haag WR, Hoigne J (1983) Ozonation of bromide-containing waters: kinetics of formation of hypobromous acid and bromate. Environ Sci Technol 17(5):261–267

    Article  Google Scholar 

  9. Brenniman G, Kojola W, Levy P, Carnow B, Namekata T (1981) High barium levels in public drinking water and its association with elevated blood pressure. Arch Environ Health 36(1):28–32

    Article  Google Scholar 

  10. Judd RM, Levy BI (1991) Effects of barium-induced cardiac contraction on large-and small-vessel intramyocardial blood volume. Circ Res 68(1):217–225

    Article  Google Scholar 

  11. Kiviat E (2013) Risks to biodiversity from hydraulic fracturing for natural gas in the Marcellus and Utica shales. Ann N Y Acad Sci 1286(1):1–14

    Article  Google Scholar 

  12. Mastrocicco M, Prommer H, Pasti L, Palpacelli S, Colombani N (2011) Evaluation of saline tracer performance during electrical conductivity groundwater monitoring. J Contam Hydrol 123(3):157–166

    Article  Google Scholar 

  13. Brantley SL, Yoxtheimer D, Arjmand S, Grieve P, Vidic R, Pollak J et al (2014) Water resource impacts during unconventional shale gas development: the Pennsylvania experience. Int J Coal Geol 126(1):17

    Google Scholar 

  14. Ferrar KJ, Michanowicz DR, Christen CL, Mulcahy N, Malone SL, Sharma RK (2013) Assessment of effluent contaminants from three facilities discharging Marcellus Shale wastewater to surface waters in Pennsylvania. Environ Sci Technol 47(7):3472–3481. doi:10.1021/Es301411q

    Article  Google Scholar 

  15. Warner NR, Christie CA, Jackson RB, Vengosh A (2013) Impacts of shale gas wastewater disposal on water quality in western Pennsylvania. Environ Sci Technol 47(20):11849–11857. doi:10.1021/es402165b

    Article  Google Scholar 

  16. Padep. DEP Office of oil and gas management compliance report. Accessed 15 Aug 2015

  17. Osborn SG, Vengosh A, Warner NR, Jackson RB (2011) Methane contamination of drinking water accompanying gas-well drilling and hydraulic fracturing. Proc Natl Acad Sci USA 108(20):8172–8176. doi:10.1073/pnas.1100682108

    Article  Google Scholar 

  18. Warner NR, Jackson RB, Darrah TH, Osborn SG, Down A, Zhao KG et al (2012) Geochemical evidence for possible natural migration of Marcellus formation brine to shallow aquifers in Pennsylvania. Proc Natl Acad Sci USA 109(30):11961–11966. doi:10.1073/pnas.1121181109

    Article  Google Scholar 

  19. Llewellyn GT, Dorman F, Westland J, Yoxtheimer D, Grieve P, Sowers T et al (2015) Evaluating a groundwater supply contamination incident attributed to Marcellus Shale gas development. Proc Natl Acad Sci 112(20):6325–6330

    Article  Google Scholar 

  20. Sang W, Stoof CR, Zhang W, Morales VL, Gao B, Kay RW et al (2014) Effect of hydrofracking fluid on colloid transport in the unsaturated zone. Environ Sci Technol 48:8266–8274

    Article  Google Scholar 

  21. Bearup LA, Navarre-Sitchler AK, Maxwell RM, McCray JE (2012) Kinetic metal release from competing processes in aquifers. Environ Sci Technol 46(12):6539–6547. doi:10.1021/es203586y

    Article  Google Scholar 

  22. Bertsch PM, Seaman JC (1999) Characterization of complex mineral assemblages: implications for contaminant transport and environmental remediation. Proc Natl Acad Sci USA 96(7):3350–3357

    Article  Google Scholar 

  23. Howarth RW, Ingraffea A, Engelder T (2011) Natural gas: should fracking stop? Nature 477(7364):271–275. doi:10.1038/477271a

    Article  Google Scholar 

  24. Steefel CI (2009) CrunchFlow software for modeling multicomponent reactive flow and transport user’s manual. Lawrence Berkeley, National Laboratory, Berkeley

    Google Scholar 

  25. Wen H, Li L, Crandall D, Hakala JA (2016) Where lower calcite abundance creates more alteration: enhanced rock matrix diffusivity induced by preferential carbonate dissolution. Energy Fuels 30(5):4197–4208. doi:10.1021/acs.energyfuels.5b02932

    Article  Google Scholar 

  26. Li L, Maher K, Navarre-Sitchler A, Druhan J, Meile C, Lawrence C et al (2016) Expanding the role of reactive transport models in critical zone processes. Earth Sci Rev. doi:10.1016/j.earscirev.2016.09.001

    Google Scholar 

  27. Swistock B (2007) A quick guide to groundwater in pennsylvania. The Pennsylvania State University, College of Agriculture Science

  28. Xu TF, Apps JA, Pruess K (2004) Numerical simulation of CO2 disposal by mineral trapping in deep aquifers. Appl Geochem 19(6):917–936

    Article  Google Scholar 

  29. Rogers RJ (1989) Geochemical comparison of ground water in areas of New England, New York, and Pennsylvania. Ground Water 27(5):690–712

    Article  Google Scholar 

  30. Denny CS, Lyford WH, Goodlett JC (1963) Surficial geology and soils of the Elmira-Williamsport region, New York and Pennsylvania, with a section on forest regions and great soil groups. USGS, Reston

    Google Scholar 

  31. Fulton J (1878) Report of progress in Bradford and Tioga counties, Harrisburg

  32. Fisher RS, Stueber AM (1976) Strontium isotopes in selected streams within Susquehanna river basin. Water Resour Res 12(5):1061–1068. doi:10.1029/Wr012i005p01061

    Article  Google Scholar 

  33. Glass HD, Potter PE, Siever R (1956) Clay mineralogy of some basal Pennsylvanian sandstones, clays, and shales: geological notes. AAPG Bull 40(4):750–754

    Google Scholar 

  34. Pirc S (1979) Uranium and other elements in the Catskill formation of east-central Pennsylvania. Pennsylvania State University

  35. Zhou Q, Birkholzer J, Rutqvist J, Tsang C-F (eds) (2007) Sensitivity study of CO2 storage capacity in brine aquifers with closed boundaries: dependence on hydrogeologic properties. In: Sixth annual conference on carbon capture and sequestration

  36. Heath R (1983) Basic ground-water hydrology: US geological survey water-supply paper 2220. US Geological Survey, Washington, DC

    Google Scholar 

  37. Domenico PA, Schwartz FW (1998) Physical and chemical hydrogeology. Wiley, New York

    Google Scholar 

  38. Williams JH, Taylor LE, Low DJ (1998) Hydrogeology and groundwater quality of the glaciated valleys of Bradford, Tioga, and Potter Counties, Pennsylvania. Pennsylvania Geologican Survey, United States Geological Survey

  39. Trapp H, Horn MA (1997) Ground water atlas of the United States: Delaware, Maryland, New Jersey, North Carolina, Pennsylvania, Virginia, West Virginia. Segment 11: the Survey

  40. Oram B (2012) Pennsylvania groundwater quality: your private well: what do the results mean? In: Halsor S, Redmond B (eds). B. F. Environmental Consultants Inc., Dallas

  41. Reed LA, Stuckey MH (2002) Prediction of velocities for a range of streamflow conditions in Pennsylvania. US Department of the Interior, US Geological Survey

  42. Schulze K, Hunger M, Döll P (2005) Simulating river flow velocity on global scale. Adv Geosci 5:133–136

    Article  Google Scholar 

  43. Watkins DM, Cornuet TS (2012) Evaluation of geology and water well data associated with the EPA hydraulic fracturing retrospective case study Bradford county, Pennsylvania. WESTON, Oklahoma city

    Google Scholar 

  44. Boyer EW, Swistock BR, Clark J, Madden M, Rizzo DE (2012) The impact of Marcellus gas drilling on rural drinking water supplies

  45. America’s Natural Gas Alliance A (2013) Bradford and Susquehanna county, Pennsylvania retrospective case study characterization report

  46. Hayes T (2009) Sampling and analysis of water streams associated with the development of Marcellus Shale gas (Final). Report to Marcellus Shale coalition. Gas Technology Institute, Des Plaines

  47. Gradient (2013) National human health risk evaluation for hydraulic fracturing fluid additives. Halliburton Energy Services, Inc., Houston

    Google Scholar 

  48. Amos RT, Mayer KU (2006) Investigating ebullition in a sand column using dissolved gas analysis and reactive transport modeling. Environ Sci Technol 40(17):5361–5367

    Article  Google Scholar 

  49. Bao C, Wu H, Li L, Newcomer D, Long PE, Williams KH (2014) Uranium bioreduction rates across scales: biogeochemical hot moments and hot spots during a biostimulation experiment at rifle, Colorado. Environ Sci Technol 48(17):10116–10127. doi:10.1021/es501060d

    Article  Google Scholar 

  50. Salehikhoo F, Li L, Brantley SL (2013) Magnesite dissolution rates at different spatial scales: the role of mineral spatial distribution and flow velocity. Geochim Cosmochim Acta 108:91–106. doi:10.1016/j.gca.2013.01.010

    Article  Google Scholar 

  51. Zheng LG, Apps JA, Zhang YQ, Xu TF, Birkholzer JT (2009) On mobilization of lead and arsenic in groundwater in response to CO2 leakage from deep geological storage. Chem Geol 268(3–4):281–297. doi:10.1016/j.chemgeo.2009.09.007

    Article  Google Scholar 

  52. Gelhar LW, Welty C, Rehfeldt KR (1992) A critical review of data on field-scale dispersion in aquifers. Water Resour Res 28(7):1955–1974

    Article  Google Scholar 

  53. Valocchi AJ, Street RL, Roberts PV (1981) Transport of ion-exchanging solutes in groundwater: chromatographic theory and field simulation. Water Resour Res 17(5):1517–1527

    Article  Google Scholar 

  54. Wollast R, Chou L (1988) Rate control of weathering of silicate minerals at room temperature and pressure. In: Physical and chemical weathering in geochemical cycles. Series C: Mathematical and Physical Sciences, vol 251. Springer Netherlands, Dordrecht, pp 11–32

  55. White AF, Brantley SL (2003) The effect of time on the weathering of silicate minerals: why do weathering rates differ in the laboratory and field? Chem Geol 202(3–4):479–506

    Article  Google Scholar 

  56. Brandt F, Bosbach D, Krawczyk-Barsch E, Arnold T, Bernhard G (2003) Chlorite dissolution in the acid pH-range: a combined microscopic and macroscopic approach. Geochim Cosmochim Acta 67(8):1451–1461. doi:10.1016/S0016-7037(02)01293-0

    Article  Google Scholar 

  57. Chakraborty S, Wolthers M, Chatterjee D, Charlet L (2007) Adsorption of arsenite and arsenate onto muscovite and biotite mica. J Colloid Interface Sci 309(2):392–401. doi:10.1016/j.jcis.2006.10.014

    Article  Google Scholar 

  58. Carrollwebb SA, Walther JV (1988) A surface complex-reaction model for the pH-dependence of corundum and kaolinite dissolution rates. Geochim Cosmochim Acta 52(11):2609–2623

    Article  Google Scholar 

  59. Perng YS, Wang EIC, Lu CC, Kuo LS (2006) Performance of sericite as a pigment for coating paper. Appita J 59(5):378

    Google Scholar 

  60. Sherman LA, Barak P (2000) Solubility and dissolution kinetics of dolomite in Ca–Mg–HCO3/CO3 solutions at 25 °C and 0.1 MPa carbon dioxide. Soil Sci Soc Am J 64(6):1959–1968

    Article  Google Scholar 

  61. Tomson MB, Fu G, Watson MA, Kan AT (2003) Mechanisms of mineral scale inhibition. Soc Pet Eng 18(3):575–582

    Google Scholar 

  62. Brandt F, Bosbach D (2001) Bassanite (CaSO4·0.5H2O) dissolution and gypsum (CaSO4·2H2O) precipitation in the presence of cellulose ethers. J Cryst Growth 233(4):837–845

    Article  Google Scholar 

  63. Setoudeh N, Zamani MAA, Welham NJ (2011) Carbothermic reduction of mechanically activated mixtures of celestite and carbon. World Acad Sci Eng Technol 50:526–529

    Google Scholar 

  64. Russell R, Rinehart D, Smith H, Peterson R (2009) Development and characterization of gibbsite component simulant. Pacific Northwest National Laboratory, Richland. Contract No.: Pnnl-18013

  65. Wolery TJ, Jackson KJ, Bourcier WL, Bruton CJ, Viani BE, Knauss KG et al (1990) Current status of the EQ3/6 software package for geochemical modeling. ACS Symp Ser 416:13

    Google Scholar 

  66. Meunier A (2005) Clays. Springer, Berlin

    Google Scholar 

  67. Li L, Steefel CI, Kowalsky MB, Englert A, Hubbard SS (2010) Effects of physical and geochemical heterogeneities on mineral transformation and biomass accumulation during biostimulation experiments at Rifle, Colorado. J Contam Hydrol 112(1–4):45–63. doi:10.1016/j.jconhyd.2009.10.006

    Article  Google Scholar 

  68. Appelo CAJ, Postma D (1993) Geochemistry, groundwater and pollution. A. A. Balkema publishers, Leiden

    Google Scholar 

  69. Palandri JL, Kharaka YK (2004) A compilation of rate parameters of water–mineral interaction kinetics for application to geochemical modeling. US Department of the Interior, Washington

  70. Brantley SL, Kubicki JD, White AF (2008) Kinetics of water–rock interaction. Springer, Berlin

    Book  Google Scholar 

  71. Li L, Salehikhoo F, Brantley SL, Heidari P (2014) Spatial zonation limits magnesite dissolution in porous media. Geochim Cosmochim Acta 126:555–573. doi:10.1016/j.gca.2013.10.051

    Article  Google Scholar 

  72. Dagan G, Fiori A, Jankovic I (2013) Upscaling of flow in heterogeneous porous formations: critical examination and issues of principle. Adv Water Resour 51:67–85. doi:10.1016/j.advwatres.2011.12.017

    Article  Google Scholar 

  73. Navarre-Sitchler A, Brantley S (2007) Basalt weathering across scales. Earth Planet Sci Lett 261:321–334. doi:10.1016/j.epsl.2007.07.010

    Article  Google Scholar 

  74. Li L, Peters CA, Celia MA (2006) Upscaling geochemical reaction rates using pore-scale network modeling. Adv Water Resour 29(9):1351–1370

    Article  Google Scholar 

  75. NJDEP (2013) Development of a dilution-attenuation factor for the impact to ground water pathway. In: Protection NJDoE (ed). NJDEP, Trenton

  76. USEPA (1996) Soil screening guidance: user’s guide. In: Agency USEP (ed) Office of solid waste and emergency response. United States Environmental Protection Agency, Washington DC, pp 89

  77. Lasaga A (1998) Kinetic theory in the earth sciences. Princeton Univ Press, Princeton

    Book  Google Scholar 

  78. Wolery TJ, Jackson KJ, Bourcier WL, Bruton CJ, Viani BE, Knauss KG et al (eds) (1990) Current status of the EQ3/6 software package for geochemical modeling. In: ACS symposium series. Oxford University Press, Oxford

  79. Vanselow AP (1932) The utilization of the base exchange reaction for the determination of activity coefficients in mixed electrolytes. J Am Chem Soc 54(4):5

    Article  Google Scholar 

  80. Appelo CAJ, Willemsen A (1987) Geochemical calculations and observations on salt water intrusions, I. A combined geochemical/minxing cell model. J Hydrol 94:313–330

    Article  Google Scholar 

  81. Bethke CM (2008) Geochemical and biogeochemical reaction modeling. Cambridge University Press, Cambridge

    Google Scholar 

  82. Zabochnicka-Świątek M, Stępniak L, Stanczyk-Mazanek E (2010) Hazard evaluation of the inorganic contaminants (B, Ba, Sr, Cu, Zn) expansion under various pH values in the quaternary aquifer. Pol J Environ Stud 2:255–260

    Google Scholar 

  83. Brady PV (1996) The physics and chemistry of mineral surfaces. CRC Press, Boca Raton

    Google Scholar 

  84. Zhang P-C, Brady PV, Arthur SE, Zhou W-Q, Sawyer D, Hesterberg DA (2001) Adsorption of barium (II) on montmorillonite: an EXAFS study. Colloids Surf A 190(3):239–249

    Article  Google Scholar 

  85. Vrionis HA, Anderson RT, Ortiz-Bernad I, O’Neill KR, Resch CT, Peacock AD et al (2005) Microbiological and geochemical heterogeneity in an in situ uranium bioremediation field site. Appl Environ Microbiol 71(10):6308–6318

    Article  Google Scholar 

  86. Edition F (2011) Guidelines for drinking-water quality. WHO Chron 38:104–108

    Google Scholar 

  87. Canada H (2014) Guidelines for Canadian drinking water quality—summary table. In: Branch HEaCS (ed) Ottawa

  88. Digiulio DC, Wilkin RT, Miller C, Oberley G (2011) Investigation of ground water contamination near Pavillion, Wyoming. Office of Research and Development, National Risk Management Research Laboratory, Ada, United States Environmental Protection Agency

  89. Kyshakevych RG, Prellwitz HS (2001) Monongahela and Youghiogheny rivers, pools 2 and 3 riverbank geology, conditions, and access reports. Carnegie Mellon University, Pittsburgh

    Google Scholar 

  90. Li L, Steefel CI, Williams KH, Wilkins MJ, Hubbard SS (2009) Mineral transformation and biomass accumulation associated with uranium bioremediation at rifle, Colorado. Environ Sci Technol 43(14):5429–5435. doi:10.1021/es900016v

    Article  Google Scholar 

  91. Haggerty R, Gorelick SM (1995) Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. Water Resour Res 31(10):2383–2400

    Article  Google Scholar 

  92. Heidari P, Li L (2014) Solute transport in low-heterogeneity sandboxes: the role of correlation length and permeability variance. Water Resour Res 50(10):8240–8264

    Article  Google Scholar 

  93. Dentz M, Cortis A, Scher H, Berkowitz B (2004) Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport. Adv Water Resour 27(2):155–173

    Article  Google Scholar 

  94. Sudicky EA, Illman WA, Goltz IK, Adams JJ, McLaren RG (2010) Heterogeneity in hydraulic conductivity and its role on the macroscale transport of a solute plume: From measurements to a practical application of stochastic flow and transport theory. Water Resour Res. doi:10.1029/2008wr007558

    Google Scholar 

  95. Luo J, Cirpka OA (2008) Traveltime-based descriptions of transport and mixing in heterogeneous domains. Water Resour Res. doi:10.1029/2007wr006035

    Google Scholar 

  96. Wang L, Li L (2015) Illite spatial distribution patterns dictate Cr(VI) sorption macrocapacity and macrokinetics. Environ Sci Technol 49(3):1374–1383

    Article  Google Scholar 

  97. Salehikhoo F, Li L (2015) The role of magnesite spatial distribution patterns in determining dissolution rates: when do they matter? Geochim Cosmochim Acta 155:107–121. doi:10.1016/j.gca.2015.01.035

    Article  Google Scholar 

Download references

Authors’ contributions

ZC carried out the reactive transport numerical simulation, analyzed the data and drafted the manuscript. LL, corresponding author, initiated the idea, helped analyze simulation data, and revised multiple versions of the draft. Both authors read and approved the final manuscript.


This work is funded by DOE National Energy and Technology Laboratory (NETL). We acknowledge Susan L. Brantley for providing information on accidental Marcellus Shale flow back/produced water spill. We thank Reviewers for their meticulous and insightful reviews that has helped improve the manuscript.

Competing interests

The authors declare no financial and non-financial competing interests.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Li Li.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Cai, Z., Li, L. How long do natural waters “remember” release incidents of Marcellus Shale waters: a first order approximation using reactive transport modeling. Geochem Trans 17, 6 (2016).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: