Earth and Planetary Science Letters 450 (2016) 1–9 Contents lists available at ScienceDirect Earth and Planetary Science Letters www.elsevier.com/locate/epsl Noble gas fractionation during subsurface gas migration Kiran J. Sathaye a , Toti E. Larson a , Marc A. Hesse a,b,∗ a b Department of Geological Sciences, University of Texas at Austin, USA Institute for Computational Science and Engineering, University of Texas at Austin, USA a r t i c l e i n f o Article history: Received 9 October 2015 Received in revised form 8 May 2016 Accepted 22 May 2016 Available online xxxx Editor: M. Bickle Keywords: noble gases shale gas two-phase flow migration–fractionation geochemical tracers geological carbon storage a b s t r a c t Environmental monitoring of shale gas production and geological carbon dioxide (CO2 ) storage requires identification of subsurface gas sources. Noble gases provide a powerful tool to distinguish different sources if the modifications of the gas composition during transport can be accounted for. Despite the recognition of compositional changes due to gas migration in the subsurface, the interpretation of geochemical data relies largely on zero-dimensional mixing and fractionation models. Here we present two-phase flow column experiments that demonstrate these changes. Water containing a dissolved noble gas is displaced by gas comprised of CO2 and argon. We observe a characteristic pattern of initial coenrichment of noble gases from both phases in banks at the gas front, followed by a depletion of the dissolved noble gas. The enrichment of the co-injected noble gas is due to the dissolution of the more soluble major gas component, while the enrichment of the dissolved noble gas is due to stripping from the groundwater. These processes amount to chromatographic separations that occur during twophase flow and can be predicted by the theory of gas injection. This theory provides a mechanistic basis for noble gas fractionation during gas migration and improves our ability to identify subsurface gas sources after post-genetic modification. Finally, we show that compositional changes due to two-phase flow can qualitatively explain the spatial compositional trends observed within the Bravo Dome natural CO2 reservoir and some regional compositional trends observed in drinking water wells overlying the Marcellus and Barnett shale regions. In both cases, only the migration of a gas with constant source composition is required, rather than multi-stage mixing and fractionation models previously proposed.  2016 Elsevier B.V. All rights reserved. 1. Introduction Noble gases are used as tracers in a variety of subsurface fluid flow processes due to their non-reactivity, low natural concentrations, and distinct isotopic signatures (Ballentine et al., 2002). Noble gas isotopes have been used to trace the fluid sources and migration pathways of hydrocarbons in sedimentary basins (Battani et al., 2000; Prinzhofer et al., 2010; Hunt et al., 2012), to quantify trapping processes in natural CO2 reservoirs (Gilfillan et al., 2008, 2009; Sathaye et al., 2014) and to identify fugitive natural gas in shallow groundwater (Jackson et al., 2013; Darrah et al., 2014, 2015). Shale gas now provides the largest share of US natural gas production (Kerr, 2010; Hughes, 2013) and the potential contamination of shallow groundwater by fugitive gases is a major public * Corresponding author at: The University of Texas at Austin, Jackson School of Geosciences, Department of Geological Sciences, 2305 Speedway, Stop C1160, Austin, TX 78712-1692, USA. E-mail address: mhesse@jsg.utexas.edu (M.A. Hesse). http://dx.doi.org/10.1016/j.epsl.2016.05.034 0012-821X/ 2016 Elsevier B.V. All rights reserved. concern (Osborn et al., 2011; Vidic et al., 2013; Jackson et al., 2013; Brantley et al., 2014). The identification of fugitive gases in drinking water aquifers is challenging because most aquifers overlying shale gas resources already contain dissolved thermogenic methane due to regional groundwater flow (Warner et al., 2012; Darrah et al., 2014; Moritz et al., 2015). A promising method to distinguish fugitive gases from natural gas migration is the use of noble gases to constrain both the migration distance and source (Darrah et al., 2014). Further development of this method requires a mechanistic understanding of the compositional changes that occur during gas migration. A second emerging energy technology is the mitigation of CO2 emissions by carbon capture and storage in deep saline aquifers (Orr, 2009; Szulczewski et al., 2012). This technology builds on decades of industry experience with gas flooding for enhanced oil recovery (Orr and Taber, 1984; Johns and Dindoruk, 2013). Several pilot projects show that CO2 injection and monitoring are feasible (Michael et al., 2010). Nonetheless, the long-term security of geological CO2 storage and the possibility of leakage are major concerns that have to be addressed (Gasda et al., 2004; Trautz et al., 2013). Natural CO2 reservoirs provide important 2 K.J. Sathaye et al. / Earth and Planetary Science Letters 450 (2016) 1–9 analogs for the fate of CO2 in the subsurface and demonstrate that very large accumulations of CO2 can be stored over millennial timescales (Sathaye et al., 2014). The interpretation of natural analogues is largely based on noble gases (Gilfillan et al., 2008, 2009) and therefore requires an understanding of their behavior during CO2 migration in the subsurface. Similarly, noble gases may provide important information about the pore space accessed and the dissolution that has occurred during the initial stages of a CO2 storage project (Györe et al., 2015). Despite the recognized importance of two-phase flow (Ricchiuto and Schoell, 1988; Hunt et al., 2012: Darrah et al., 2014, 2015), noble gases are commonly interpreted using zero-dimensional mixing and fractionation models (Gilfillan et al., 2008, 2009; Darrah et al., 2014). Currently the relationship between these models and the underlying physical processes that control the evolution of gas composition is not clear. The aim of this manuscript is therefore to provide an understanding of the behavior of noble gases in twophase flow. To this end, Section 2 presents well-characterized two-phase displacement experiments that demonstrate the compositional changes due to partitioning of noble gases between the phases. In Section 3 these observations are first discussed using mass-balance arguments, before we outline the application of gas injection theory to noble gas fractionations in Section 4. This theory predicts the observed compositional trends and provides a framework for understanding post-genetic changes in gas composition. Finally, the experimental observations and theoretical predictions are compared with published observations in natural CO2 reservoirs and groundwater overlying shale gas resources in Section 5. Fig. 1. Experimental elution curves from a two-phase displacement of neonsaturated water by CO2 –argon mixture. Gas volume flow rates at the column outlet are shown as function of pore volumes injected: A) Neon, B) Argon, and C) CO2 . Subplots A and B show the arrival of two highly enriched noble gas banks ahead of the main CO2 front. 2. Two-phase gas fractionation experiments Previous observations from shale gas production and geological CO2 storage show that gas migration leads to the enrichment of both noble gases initially dissolved in the brine and those that are introduced with the gas phase (Gilfillan et al., 2008, 2009; Darrah et al., 2014, 2015). To investigate the mechanism that leads to the simultaneous enrichment of both co-injected and dissolved atmospheric noble gases, we performed four-component two-phase column displacement experiments. Gas with a CO2 /argon ratio of 50 was injected with a constant volume flow rate. The gas displaced the water from the top of the column towards the bottom to obtain a gravity stabilized displacement. The effluent liquid and gas were collected in helium-filled glass vials and subsequently analyzed in a gas chromatograph. For a detailed description of the experimental method see supplementary materials Section 2. Two experiments with different initial dissolved gases were conducted. Experiment 1 was performed in a 122 cm long and 1 cm inner diameter glass column packed with glass beads of 150–212 micron diameter. The column was initially filled with water saturated with 2 Atm of neon. Experiment 2 was performed in a steel column of length 100 cm and 1 cm inner diameter, filled with quartz sand with a mean grain diameter of 100 microns. The column was initially filled with water saturated with CH4 at 2 Atm. We chose CH4 as the dissolved gas in the second experiment, because its higher solubility leads to a larger chromatographic separation. Figs. 1 and 2 show the results of the four-component gas displacement experiments. To account for the finite sampling intervals, gas composition is presented in units of volume flow rate as a function of pore volumes injected. Gas breakthrough occurs before one pore volume has been displaced, indicating the presence of pore water behind the gas front (Fig. 2A). Typical residual water saturations in the column at the end of the experiment were 40%. Two distinct noble gas banks are seen in the first experiment. Fig. 1A shows the effluent flow rate of neon, which was initially Fig. 2. Experimental elution curves from a two-phase displacement of CH4 -saturated water by CO2 –argon mixture. Volume flow rates at the column outlet are shown as function of pore volumes injected: A) Water, B) Argon and CH4 , and C) CO2 . Subplot B shows the arrival of two highly enriched gas banks ahead of the main CO2 front. present as a dissolved component in the water. The first gas eluting from the column is highly enriched in neon that has been stripped from the residual water remaining behind the gas front. Shortly thereafter, the injected gas arrives and is enriched in argon, due to the dissolution of the highly soluble CO2 (Fig. 1B and 1C). Finally, the retarded CO2 front breaks through and the CO2 /Argon ratio reverts to the injected ratio of 50. In the second experiment, water with dissolved CH4 was displaced with the same CO2 /argon mixture (Fig. 2). The initial mass of dissolved gas, in experiment 2 is greater due to the higher solubility of CH4 relative to neon. The stripping of this dissolved gas leads to the formation of a larger and more distinct bank (Fig. 2B) that coincides with a decrease in the volume flow rate of water (Fig. 2A). The argon bank is clearly separated from the CH4 bank and precedes the CO2 front, similar to experiment 1. After the CO2 breakthrough, the gas effluent composition matches the injected gas and CH4 has been completely stripped from the water remaining in the column. Note, the small CH4 concentrations ahead of the gas front are due to exsolution of dissolved CH4 into the sampling 3 K.J. Sathaye et al. / Earth and Planetary Science Letters 450 (2016) 1–9 head space, giving an indication of the initial dissolved concentration. This illustrates the reduction of CH4 concentration in the residual water behind the gas front, relative to the initial dissolved concentration. Our experiments show two types of behavior that are compatible with the overarching trends that emerge from previous work on noble gas fractionations associated with shale gas production and geological CO2 storage: First, co-injected noble gas isotopes that are less soluble than the main gas component become increasingly enriched as migration proceeds. This is due to the preferential dissolution of the more soluble major gas components, an effect that is pronounced in CO2 -rich gases. Second, atmospheric noble gases are stripped from air-saturated groundwater and enriched above atmospheric values as migration proceeds. This combination of processes enriches the front of the migrating gas in both co-injected and atmospheric noble gases, but depletes atmospheric noble gases behind the gas front. 3. Mass-balance interpretation of experimental results The quantitative interpretation of these experiments requires a coupling of the dynamics of two-phase flow with the partitioning of the components between the differentially moving phases. Before addressing this coupling in Section 4 we will first show that the basic features can be understood by considering the mass balance in a plug-flow. We consider a piston like displacement with a constant residual saturation, s wr . The length of the column is ℓ and its pore volume is V p . Gas of uniform composition is injected at the constant pressure, p g , and constant volumetric rate, Q . Assuming an ideal gas, the molar injection rate into the column is given by p g Q /( R T ), where R is the gas constant and T is the temperature of the injected gas. Due to the low viscosity of the gas most of the pressure drop occurs in the water ahead of the gas front and we will assume that the gas pressure is constant at the injected value throughout the column. In both experiments gas breakthrough, tb , is observed before one pore-volume has been injected, tb < V p / Q (Figs. 1 and 2). This is a result of the residual water saturation, s wr , which reduces the pore space accessible to the gas to V g = V p (1 − s wr ). In the absence of dissolution, mass balance requires that gas breakthrough occurs at tb∗ = V g / Q after nb = p g V g /( R T ) moles have been injected (Fig. 3A). Typical residual saturations in our experiments are 0.4, suggesting that an immiscible gas should break through after 0.6 pore volumes. However, the observed gas breakthrough in the experiment occurs at pore volumes between 0.70 and 0.75. This delayed arrival of the gas is due to dissolution of the injected CO2 (Fig. 3B). The fraction of the gas that remains in the gas phase is given by Fi = 1 − s wr 1 + ( R T H i − 1)s wr , (1) where H i is the Henry’s law solubility constant for species i. If nb moles of a pure gas are injected, its position is given by xi = Fi ℓ. The breakthrough will occur after tb = tb∗ /Fi pore-volumes have been injected. In the case of pure CO2 , breakthrough is predicted to occur after 0.93 pore-volumes have been injected. If argon is co-injected, it will advance further than the more soluble CO2 , according to Equation (1). In a piston-like displacement with constant gas saturation behind the front, the argon must accumulate in a bank just ahead of the CO2 to maintain constant gas pressure (Fig. 3C). The width of this argon bank is xAr = p Ar / p g (FAr − FCO2 )ℓ, where p Ar is the partial pressure of argon in the injected gas and p g is its total pressure. The front of the injected gas after injection of nb moles is located at xinj = xCO2 + xAr . The breakthrough is expected at tb = tb∗ ℓ/xinj . Fig. 3. Plug-flow mass balance models with increasing number of components. All figures show injection of a gas from left to right into a water-filled porous medium. Displacement is assumed to be a plug flow with residual water saturation, s wr = 0.4. A) Injection of an immiscible gas into pure water. B) Injection of pure CO2 into pure water. Dissolution reduces the length of the gas plume to xCO2 /ℓ = 0.63. C) Injection of a 80% CO2 and 20% argon into pure water. Accumulation of low-solubility argon increases the length of the gas plume to xinj /ℓ = 0.72. D) Injection of 80% CO2 and 20% argon into neon-saturated water. Neon degassed from the residual water forms a bank at the front of the gas plume, increasing the plume length to xg /ℓ = 0.73. The lengths of all gas plumes are scaled to the column length ℓ. The Henry’s law coefficients at 25 ◦ C are H Ar = 1.4 · 10−3 , H Ne = 4.5 · 10−4 , and H CO2 = 0.034 mol/(L · atm). For the case of injection of CO2 and Ar at a ratio of 4:1 the breakthrough occurs at approximately 0.84 pore volumes injected. Finally, if the initial water contains a dissolved noble gas such as neon, it will be stripped from the residual water and accumulate in a bank of pure neon at the front of the migrating gas (Fig. 3D). The width of this neon bank is given by xNe = xinj s wr H Ne p Ne (1 − s wr ) p g /( R T ) + s wr H Ne ( p g − p Ne ) , (2) where p Ne is the partial pressure of neon that the water has been equilibrated with. The position of the gas front after the injection of nb moles is therefore x g = xinj + xNe and the breakthrough occurs after tb = tb∗ ℓ/x g moles have been injected. This leads to breakthrough of the gas after the injection of 0.83 pore volumes. The discussion above used a injected gas that contained more argon than the experiments to make the argon bank visible in Fig. 3. With the experimental values we expect neon, argon and CO2 breakthrough after 0.93, 0.92, and 0.91 pore volumes injected. The CO2 breakthrough is 0.16 pore volumes later than that observed experimentally in Fig. 1. Some of this discrepancy may be due to the simplifying assumption of a plug flow. However, we show in Section 4.6 that the gas saturation changes very little after breakthrough. Therefore, the late breakthrough suggests too much CO2 dissolves in the mass-balance calculation. Given that equation (1) is independent of the gas pressure and the Henry’s law solubility constant for CO2 is well known, it is likely that the CO2 has not accessed the entire residual water at breakthrough. This suggests that aqueous CO2 diffusion is not fast enough to transport CO2 from the gas–water interface into the entire residual water volume during the two hours required for CO2 breakthrough. Using equation (1) and a CO2 breakthrough of 0.75, we estimate that the effective residual saturation to be s∗wr = 0.25. 4. Theory of gas injection The analysis of the gas fronts presented in Section 3 assumes a plug flow with constant saturation and is based on additional implicit assumptions such as the ordering of the fronts and the 4 K.J. Sathaye et al. / Earth and Planetary Science Letters 450 (2016) 1–9 formation of pure noble gas banks. These are motivated by experimental observations but difficult to justify otherwise. Analytic solutions to the governing conservations laws for two-phase flow with partitioning of components provide these answers. Such a theory of gas injection has been developed in petroleum engineering and successfully predicts compositional changes in laboratory experiments with fluid–fluid displacements (Welge et al., 1961; Wachmann, 1964; Dumore et al., 1984; Laforce et al., 2010). Insights derived from this theory have been applied to design gasinjection operations in the oil and gas industry (Johns and Dindoruk, 2013). However, it is currently not recognized that the theory of gas injection also provides a framework to understand the noble gas fractionations observed above. Due to the non-linearity of two-phase flow and the coupling of multiple partial differential equations, the full theory is rather complex (Lake, 1989; Orr, 2007). Here, we can only outline the development and illustrate it by constructing the particular analytical solution for the experiments discussed above. In Section 4.1 we review the basic characteristics of immiscible two-phase displacements that form the foundation for the theory of gas injection. Sections 4.2 to 4.4 review the basic framework of the theory of gas injection, before we develop the analytic solution relevant to our experiments in Section 4.5. Fig. 4. Solution for immiscible gas injection into water. A) Relative permeabilities to gas and water, krg , kr w . B) Fractional flow of gas phase, f g , with M = 10 and s wr = 0.2. The shock saturation for initial gas saturation of zero is shown by the tangent point to the fractional flow curve. C) Derivative of the fractional flow function with respect to gas saturation df g /ds g . D) Gas saturation profile as a function of speed, assuming q/φ = 1. E) Gas saturation profiles at various times. The profile at a given time is obtained by multiplying the propagation speed by t. F) Outlet gas (s g ) and water (s w ) saturation as a function of pore volumes injected (PVI). PVI is proportional to time for an injection with constant flow rate. Therefore, at the end of the column, x = 1, PVI is equal to the inverse of speed. 4.1. Saturation profiles in two-phase displacements Fig. 2A shows a typical water effluent history for a two-phase displacement experiment. Before the injected gas has reached the end of the column, only water is eluted. The injected gas breaks through before one pore-volume has been injected, indicating that significant amounts of water remain in the column. The flux of water drops significantly during breakthrough and afterwards the flux declines slowly towards zero. Below we show that these basic features of the effluent history are captured by Buckley– Leverett theory for incompressible and immiscible two-phase flow in a porous medium (Buckley and Leverett, 1942; Lake, 1989; Orr, 2007). During two-phase flow, gas and water only occupy a fraction of the pore-space, called their saturation 0 ≤ s w , s g ≤ 1, where s w + s g = 1. The relative permeability of each phase, krp , is proportional to its saturation (Wyckoff et al., 1961; Leverett, 1939) (Fig. 4A). For the primary drainage, i.e., the injection of a gas into a porous medium fully saturated with water, the relative permeabilities are given by ∗ krg = krg  s g − s gc 1 − s wr − s gc n g and kr w =  1 − s g − s wr 1 − s wr − s gc n w , (3) ∗ is the end-point relative permeability of the gas, and where krg n g and n w are the relative permeability exponents of the gas and water phases, respectively. In the column flow experiments considered here the residual water saturation is typically 0.4. The critical gas saturation s gc is the minimum gas saturation required for the gas phase to flow. Consider the flow of two immiscible and incompressible phases in a homogeneous porous column in the absence of gravity. The mass balance of the gas is given by φ ∂ f gq ∂ sg + = 0 or ∂t ∂x ∂ sg q d f g ∂ sg + = 0, ∂t φ ds g ∂ x (4) where q is the volumetric flux, f g is the fractional flow of the gas, and φ is the porosity (Leverett, 1941). The fractional flow of gas is controlled by the relative permeabilities and viscosities of the two fluids, f g = krg /(krg + Mkr w ), where M = μ g /μ w is the viscosity ratio. The fractional flow curve shown in Fig. 4B illustrates that the fraction of fluid transport in the gas phase is not linearly proportional to the gas saturation. At very low gas saturations, the low relative permeability prevents efficient transport in the gas, so that f g < s g . As the gas saturation increases, transport in the gas becomes dominant, f g ≫ s g . This is due to the low viscosity of the gas relative to water, M ≪ 1. Equation (4) is an advection equation for s g . The speed at q df which each gas saturation propagates is given by λ(s g ) = φ ds g , g where q/φ is the interstitial fluid velocity in the saturated porous medium ahead of the gas front. Equation (4) is non-linear because λ depends on s g as shown in Fig. 4C. The location of each saturation along the column at any particular time, t, is given by x = λ(s g )t. This leads to the multi-valued saturation profile shown in Fig. 4D, because λ is a non-monotonic function of s g . To avoid this, a shock front must exist and can be found by the tangent construction shown in Fig. 4B (Buckley and Leverett, 1942; Lake, 1989; LeVeque, 1992). This shock propagates at a speed equal to the slope of the tangent line. For injection into a fully saturated column the shock speed is given by (ŝ g ) = f g (ŝ g )/ŝ g , were ŝ g is the shock saturation. The shock saturation is determined by the tangency condition, i.e., the requirement that the propagation speed is equal to the shock speed, λ(ŝ g ) = (ŝ g ). This construction gives rise to a gas saturation effluent profile where an initial rapid change in saturation is followed by a gradual increase (Fig. 4D and E). The arrival of the shock in gas saturation leads to an initial rapid drop in the flux of water eluted that is followed by a gradual decline. This matches the patterns of eluted water observed in our experiments (Fig. 2A). These basic features are retained even if the gas is compressible and in the presence of weak gravitational forces. The basic concepts outlined in this section can also be applied to a two-phase displacement with partitioning of components. In this case, multiple coupled advection equations must be solved. For the conditions relevant to our experiments, this leads to the formation of multiple shock fronts which bound the two noble gas banks that are observed. The main difficulty in the multi-component case is the identification of the correct solution path in composition space as outlined below. 5 K.J. Sathaye et al. / Earth and Planetary Science Letters 450 (2016) 1–9 4.2. Governing equations for multi-component two-phase flow In the multi-component gas–liquid system the bulk composition of the pore fluids is given by C = φ(1 − s g )cl + φ s g c g , (5) where φ is the porosity, s g is the gas saturation and the compositions within the liquid and gas phases are given by cl and c g , respectively. In this notation, bold symbols indicate a vector of compositions. Individual components within this vector will be indicated by an additional subscript, so that cl,Ar , c g ,Ar and C Ar are the concentrations of argon in the liquid, gas and mixture, respectively. At local chemical equilibrium, the partitioning of the components between the phases is governed by a set of partition coefficients, such that cl,i = H i c g ,i , where H i is Henry’s law coefficient of component i. If C and the H i values are known, the gas saturation, s g , and the phase compositions, cl and c g , can be determined (Whitson and Michelsen, 1989). In the limit of negligible hydrodynamic dispersion and capillary diffusion, local chemical equilibrium, incompressible phases, and a homogeneous porous medium, the change in C due to two-phase flow with partitioning between the phases is governed by a system of coupled non-linear hyperbolic partial differential equations ∂ C ∂ F(C) + = 0, ∂t ∂x (6a) where F(C) is a flux vector that only depends on C, so that the system is quasi-linear (LeVeque, 1992). This system of conservation equations is the extension of Equation (4) to the multi-component case. The entries of F are the fluxes of the individual chemical components and F is not to be confused with the flux vector in a multi-dimensional flow problem that indicates the direction of flow. Figs. 5A and 5B show that these fluxes are piece-wise defined (Orr, 2007): linear in the single-phase region and non-convex in the two-phase region. In the two-phase region the fluxes of the components are closely related to the fractional flow function discussed in Section 4.1. While the theory is developed for an arbitrary number of components (Orr, 2007), a four-component system is sufficient for the discussion of the experimental data presented above and the field data discussed below. This allows the graphical presentation of the analysis in the three-dimensional tetrahedron representing composition space, shown in Fig. 5C. In all experiments discussed in Section 2, the injected gas comprises a major component, G, and a co-injected noble gas, g. The initial liquid comprises the major component, L, i.e., water, and a dissolved gas, l. The initial condition representing the injection of a two-component gas into a liquid containing a dissolved gas can be written as C=  Cinj = (C G , C g , 0, 0), x < 0, Cini = (0, 0, C L , C l ), x > 0. (6b) Given this particular initial condition, the analytic solution is restricted to two of the triangular surfaces of the tetrahedron. The analysis can be presented in a two-dimensional diamond, as shown in Fig. 5C. For the theoretical example considered below, the Henry’s law coefficients are ordered as follows: H g < H l < H G ≪ H L , so that the minor gases are less soluble than the main gas component. Fig. 5. A) Flux function for the main liquid component L along the injection tie line is shown in red. The injection tie line connects C1 and C3 in subfigure C. The tangent constructions for waves W1 and W2 are shown by black lines and the corresponding tangent points are C2 and C3 . B) Flux functions for the g–L crossover and l–L initial tie lines are shown in green and blue. The chords representing waves W3 and W4 are shown in black. C) The 3D figure in the upper left shows the four-component tetrahedron representing composition space. The large 2D diamond shows the g–G–L and g–L–l faces of the tetrahedron. The red and blue shaded regions represent the single-phase gas and liquid regions, respectively. The tie lines connecting the equilibrium liquid and gas compositions in the white two-phase region are shown as dashed black lines and form a set of composition paths. The other set of composition paths is shown as full black lines. The solution path for initial condition (6b) and solution structure (9) is shown as a green line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) the eigenvalues of the Jacobian matrix, ∇c F, are not ordered and therefore system (6a) is not strictly hyperbolic. This introduces an additional transitional wave and associated intermediate state (Isaacson et al., 1990). The general solution of the four-component two-phase flow problem has the following structure W1 W2 W3 W4 Cinj −−→ C2 −−→ C3 −−→ C4 −−→ Cini (7) and comprises four waves, W p , that separate three regions of constant composition given by the intermediate states, C2 , C3 , C4 . From a geochemical and observational perspective, the composition of the intermediate states is of greatest interest. Due to the non-linear couplings in (6a), the compositions of these intermediate states are not bounded by the initial and injected compositions. This allows noble gas fractionations that evolve towards compositions that are not simple mixtures between the initial and injected compositions. 4.4. Solution construction in composition space Similar to incompressible two-phase flow discussed in Section 4.1, the solution of (6) is self-similar and the concentration profiles simply stretch over time in the direction of displacement (Orr, 2007). This allows the reduction to a non-linear eigenvalue problem 4.3. Solution structure (∇c F − λ p I)r p = 0, The solution of a genuinely non-linear system of three strictly hyperbolic equations comprises three waves or fronts that separate regions of constant composition (LeVeque, 1992). However, The eigenvectors, r p = dC/dη , determine directions in composition space that allow all components to travel with the same speed, given by λ p (Rhee et al., 2001). These eigenvalues are analogous p ∈ [1, 2, 3]. (8) 6 K.J. Sathaye et al. / Earth and Planetary Science Letters 450 (2016) 1–9 to the propagation speed of the saturation front, λ, in Section 4.1. Here η is a arc-length variable in composition space. The solution is constructed in composition space where the eigenvectors of (8) define the grid of composition paths shown in Fig. 5C. Starting with Cinj , the solution follows paths with increasing speeds towards Cini . The intermediate states, C2 , C3 and C4 , are located at the points where the solution switches paths. If λ p increases along the path, the associated wave is a continuous rarefaction, R p , that spreads with time. In contrast, a decrease of λ p along the path leads to the formation of a discontinuous self-sharpening shock, S p , that propagates with a speed given by the jump condition  p = [F] p /[C] p , where the brackets indicate the jump in the quantity across the shock. Finally, the non-convex two-phase flux function, shown in Fig. 5A, allows the formation of composite waves RS p , SR p , or SRS p . The Buckley–Leverett solution for immiscible two-phase flow discussed in Section 4.1 is such a composite wave, namely a rarefaction-shock, RS . Once the composition path of the solution is known, the concentration profiles are obtained from the variation of λ p along the path and the jump condition for the shocks. For details of the solution construction, see SI Appendix Section 1 and the following references (Johansen et al., 2005; Orr, 2007). 4.5. Solution for noble gas fractionation during gas migration The initial condition (6b) corresponds to the migration of a gas phase containing a noble gas though groundwater containing a different dissolved (noble) gas. Fig. 5C shows the composition path of the solution and Fig. 6A the gas saturation profile. The wave structure of this solution is given by SR1 S2 S3 S4 Cinj −−−→ C2 − −→ C3 −−→ C4 −−→ Cini . Fig. 6. Self-similar saturation and composition profiles for initial condition (6b) and solution structure (9) are given as function of propagation speed. A) Gas saturation along the direction of displacement. The speeds of the discontinuous waves are given by the chords in Fig. 5A and 5B. The speed along the continuous portion between △ and ✷ is given by dF L /dC L in Fig. 1A. Gas B) and liquid C) phase concentration of components along the direction of displacement. The three regions of constant composition correspond to the three intermediate states. Two banks highly enriched in l and g form at the front of the displacement. No gas is present in the white area corresponding to the initial liquid, ahead of W4 . Note that g is present in the liquid phase to the left of W2 , but its concentration is too low to be visible. In contrast, the dissolved gas l is only present in the liquid to the right of W3 . (9) The gas saturation shows a typical Buckley–Leverett profile with gradual increase in gas saturation along SR1 , as discussed in Section 4.1. The slow increase in gas saturation is due to the low viscosity of the gas relative to the displaced liquid. However, ahead of the main saturation shock, S2 , is a region of lower saturation. The concentration profile in Fig. 6B shows that the main gas component, G, is absent because it has dissolved into the residual water behind the main shock. The low-saturation region is divided into two compositionally distinct segments by S3 . In the section directly ahead of the main shock, the composition is given by C3 and is highly enriched in the injected noble gas, g, after the dissolution of G. The front of the low-saturation region is highly enriched in the atmospheric noble gas, l, and its composition is given by C4 . This noble gas has been stripped from the groundwater inside the gas plume and is concentrated at the front of the migrating gas to values much larger than the initial air saturated values (Fig. 6C). 4.6. Comparison with neon–argon–CO2 experiment The previous section presents a four-component gas displacement with idealized coefficients to illustrate the main features of the solution. In particular, the solubilities had to be exaggerated to make the single-phase regions visible in Fig. 5. The concentration profiles are also shown as a function of velocity, similar to Fig. 4D. In contrast, the experimental observations show effluent concentration as a function of volume injected. To facilitate comparison between theory and experiments, Fig. 7 shows theoretical elution curves using the appropriate solubilities of neon, argon, and CO2 for the experiment shown in Fig. 1. Similar to the experimental results, this solution shows a small bank of neon at first gas arrival, followed by a larger bank of argon, and finally the arrival of the injected CO2 –argon gas mixture. Matching the arrival times of the different fronts requires a nearly constant Fig. 7. Analytical gas injection solution approximating Experiment 1. A) Relative permeability curves for liquid and gas phases. The liquid phase relative permeability (blue) has an exponent n w = 1.05 and residual water saturation of 0.4. Gas phase relative permeability (red) has an exponent n g = 3 with a critical gas saturation of ∗ ) of 0.2. B) Fractional flow of the gas 0.2 and end point relative permeability (krg phase as a function of gas saturation. The gas:liquid viscosity ratio is 60. C) Theoretical effluent curve using realistic solubilities of neon, argon, and CO2 at standard pressure and temperature. The shaded areas show the composition of the effluent gas, while the blue curve represents the saturation of the liquid phase at the column outlet. Similar to the mass balance model, this creates a narrow neon bank, followed by a bank of pure argon, before the system reverts to the injected CO2 – argon mixture. After the end of the pure argon bank, the liquid phase saturation s L slowly decreases towards the residual saturation of 0.4. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) and high gas saturation within the gas plume. In the two-phase flow model, this is achieved by reducing the mobility of the gas at low saturations. The combination of a nearly constant gas saturation within the plume and noble gas banks of pure composition K.J. Sathaye et al. / Earth and Planetary Science Letters 450 (2016) 1–9 makes the simple mass balance model of Section 3 a useful approximation. Similar to the mass balance model, the gas injection solution approximation to the experiment arrives later than the experimental gas arrival, due to incomplete CO2 dissolution into the residual water during the timescales of the experiment. The theoretical solution neglects the effects of dispersion and therefore predicts sharp fronts that separate pure noble gas banks. In the experiments, both hydrodynamic dispersion and the discrete sampling intervals contribute to smearing of the compositional fronts. The theoretical model predicts banks of similar concentration and differing width, while the experiment shows overlapping banks of similar width and differing concentration. This indicates that the length-scale of smearing is comparable to the width of the noble gas banks. The maximum concentration of argon and neon within the banks does not reach their theoretical value due to this smearing effect. The ratio of the maximum observed concentrations in the banks is comparable to the ratio of the noble gas bank widths in the theory. 5. Field observations of noble gas fractionation in two-phase flow Here we introduce three field data sets from previous publications. These studies illustrate the results of noble gas fractionation during two-phase flow of natural gas through groundwater. The first data set contains samples of gas composition from a continuous gas emplacement process at Bravo Dome, NM. The second and third data sets contain dissolved gas measurements from shallow groundwater in the Marcellus and Barnett shales. The observed distributions in these groundwater areas are likely the result of multiple gas migration events. As such, the Bravo Dome case is conceptually closer to the experimental results of Section 2. 5.1. Bravo Dome natural CO2 reservoir Bravo Dome is a large natural CO2 reservoir that contains approximately 1.3 Gt of almost pure magmatic CO2 trapped beneath a regional evaporite seal (Broadhead, 1990; Sathaye et al., 2014). In map view the flow was approximately radial away from the pressure maximum in the west (Fig. 8A). After entering the reservoir, the gas flowed nearly horizontally beneath the cap rock for approximately 60 km (Fig. 8C). Noble gases in Bravo Dome have been studied extensively and provide constraints on the dynamics of the reservoir (Gilfillan et al., 2008, 2009). The concentrations of atmospheric 20 Ne are lower than expected from degassing and require that CO2 entered the reservoir as a free gas phase that stripped 36 Ar, 20 Ne and other dissolved gases from the reservoir brine. The extreme fractionation of 20 Ne/36 Ar is interpreted as evidence for re-dissolution of noble gases stripped from the brine. The noble gas dynamics are therefore interpreted as a two-stage process of initial stripping and later re-dissolution (Gilfillan et al., 2008). Recent mass balance estimates show that significant dissolution must have occurred during CO2 emplacement (Sathaye et al., 2016, 2014). The effective residual water saturation in the field is similar to that of the laboratory experiments, approximately 40%, and should cause fractionation of noble gases during the emplacement. Bravo Dome provides an opportunity to test the theoretical predictions against field data, because of the simple filling history, the long migration distance, and the large effective residual brine saturation. In the context of the theory introduced above, the main gas component, G, corresponds to CO2 , the noble gas isotope introduced with the gas, g, is 3 He and the dissolved noble gas isotope, l, is 20 Ne. Figs. 8B and 8D show concentration profiles of previously published 3 He and 20 Ne data in the Bravo Dome gas as a function 7 Fig. 8. A) Bottom hole pressure at the onset of commercial CO2 extraction from the Bravo Dome reservoir interpolated from 373 wells (Sathaye et al., 2014). The high pressure zone and noble gas isotopes suggest that the CO2 source is in the west. The markers indicate locations of wells that were sampled for noble gases (Cassidy, 2005). B) 3 He and 20 Ne concentrations from the Bravo Dome natural CO2 field as a function of distance from entry point shown in subfigure A. Both isotopes are enriched in banks at the gas front. of approximate migration distance (Gilfillan et al., 2008). Close to the source, 3 He is constant at the injected value. At the front of the migrating gas the dissolution of CO2 has increased the concentration of 3 He. Close to the source of the gas, 20 Ne concentrations are significantly reduced relative to the atmospheric levels ahead of the gas. The 20 Ne stripped from the groundwater in the interior of the migrating gas plume has accumulated at the front. Here, its concentration is significantly higher than the initial atmospheric value in the brine leading to re-dissolution. These field data are consistent with the experimental observations (Figs. 1B and 2B), and theoretical predictions of noble gas enrichments (Fig. 6B). The model of combined multi-phase flow and partitioning shows that stripping and re-dissolution can occur simultaneously in different parts of the migrating gas. Monitoring of the produced gas during the Cranfield CO2 enhanced oil recovery (CO2 -EOR) project in Mississippi showed similar enrichments of 4 He, 3 He and 20 Ne in the initial gas after breakthrough. This illustrates that the formation of noble gas banks also occurs on engineering timescales and may therefore provide and important monitoring tool for gas injection projects (Györe et al., 2015). 5.2. Dissolved gases in groundwater above the Marcellus and Barnett shales In regions where shallow groundwater contains dissolved thermogenic methane, dissolved noble gases have been invoked as a tool to differentiate natural methane migration from fugitive gases associated with shale gas development (Darrah et al., 2014). In this case, the thermogenic source gas characteristics are identical, but the natural methane migration is though to occur over much longer distances. Therefore, an understanding of compositional changes during migration is necessary to discriminate these two processes. In this case, gas composition has to be inferred from the dissolved gases measured in the groundwater. Due to near instantaneous equilibration, the dissolved gas concentrations mirror the gas phase compositions. During gas migration from the reservoir to the shallow groundwater, more soluble gas components such as CH4 dissolve into the water so that less soluble components, such as 4 He, become enriched in the gas. Simultaneously, atmospheric noble gases initially dissolved in the groundwater, such as 20 Ne, are stripped from the water and accumulate at the front of the gas plume. This process is 8 K.J. Sathaye et al. / Earth and Planetary Science Letters 450 (2016) 1–9 Fig. 9. A) Gas phase compositions for experimental two-phase displacement of neon-saturated water by a CO2 –argon mixture shown Fig. 1. Similar to the Marcellus field data, dispersion causes the two noble gas banks to overlap, leaving the front of the plume simultaneously enriched in argon (g) and neon (l). B) Marcellus and C) Barnett groundwater noble gas concentrations (Darrah et al., 2014), normalized to concentration in air-saturated water. The solid black line connects estimated dissolved gases in water at equilibrium with Marcellus reservoir gas (Hunt et al., 2012), the theoretical noble gas bank, and air-saturated water. The noble gas bank represents the saturation limit of dissolved 4 He and 20 Ne in shallow groundwater. Note that the source noble gas composition for the Barnett shale is unknown, and is assumed to be identical to the Marcellus gas. analogous to the theoretical displacement considered above, where CH4 acts as G, 20 Ne acts as l, and 4 He as g. While theory predicts chromatographic separation between 4 He and 20 Ne banks (Fig. 6B and 6C), our two-phase flow experiments show that hydrodynamic dispersion will cause the noble gas banks to overlap (Fig. 1). This leads to a simultaneous enrichment of both 4 He and 20 Ne at the gas front (Fig. 9A). This co-enrichment process leads to noble gas concentrations which cannot be explained using simple mixing models between air-saturated water and thermogenic natural gas. Marcellus groundwater samples show similar co-enrichment of 4 He, 20 Ne and CH4 (Fig. 9B). This suggests that the groundwater contacted methane that had migrated a substantial distance, indicating natural migration (Darrah et al., 2014). In contrast, some samples from the Barnett shale lack this co-enrichment suggesting a shorter migration pathway with less noble gas exchange. The samples with the highest CH4 concentrations are depleted in 20 Ne and enriched in 4 He relative to air-saturated groundwater. Previous work has interpreted depletion of dissolved atmospheric noble gases as a signal of fugitive gas contamination (Darrah et al., 2014). The counter-intuitive observation that natural gas migration can either enrich or deplete dissolved atmospheric noble gases is explained by the experiments and theory presented here. Gas migration simultaneously enriches 20 Ne (l) at the front of the gas plume and depletes it in the interior (Figs. 6B and 1B). The groundwater samples in the Barnett shale (Fig. 9C) with depleted 20 Ne represent the stripped interior of the gas plume. In the case of fugitive gases associated with faulty well construction, the migration distance is reduced and the accumulation of dissolved atmospheric noble gases during migration to the shallow groundwater is limited. Upon arrival in the shallow groundwater, the unaltered methane is still depleted in atmospheric noble gases relative to air saturated groundwater, thereby stripping the dissolved noble gases from the shallow groundwater. 6. Conclusion We combine experiments and theory to explain noble gas fractionation observed in the field. These separations are a result of two-phase flow in the subsurface. Understanding these changes is required to fingerprint subsurface gas sources. This understanding is essential for environmental monitoring of shale gas production and geological CO2 storage, as well as hydrocarbon exploration. Two-phase flow experiments and gas injection theory show that dissolved gases are stripped from groundwater within the migrating gas and enriched in a bank at the front. Similarly, insoluble components co-injected with the migrating gas are also enriched at the front, due to the preferential dissolution of the more soluble main gas component. This leads to a characteristic pattern of co-enrichment of insoluble gas components at the front of migrat- ing gas, followed by depletion of initially dissolved gases from the groundwater. This pattern has been observed in concentrations of atmospheric and co-injected noble gases in natural CO2 reservoirs (Gilfillan et al., 2009; Cassidy, 2005), CO2 -EOR projects (Györe et al., 2015), conventional gas fields (Hunt et al., 2012), and in shallow groundwater overlying shale gas resources (Darrah et al., 2014). These compositional changes provide a powerful tool to monitor fugitive gases during shale gas development and to estimate dissolution trapping during geological CO2 storage. Gas injection theory provides a framework for noble gas fractionation during subsurface gas migration and hence the interpretation of noble gas observations in the field. Acknowledgements This work was supported as part of the Center for Frontiers of Subsurface Energy Security, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award # DE-SC0001114. K.J.S. acknowledges a doctoral fellowship through the Statoil Fellows Program at The University of Texas at Austin. The authors thank Esben Pedersen for help with laboratory experiments. Appendix A. Supplementary material Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.epsl.2016.05.034. References Ballentine, C.J., Burgess, R., Marty, B., 2002. Tracing fluid origin, transport and interaction in the crust. Rev. Mineral. Geochem. 47 (1), 539–614. Battani, A., Sarda, P., Prinzhofer, A., 2000. Basin scale natural gas source, migration and trapping traced by noble gases and major elements: the Pakistan Indus basin. Earth Planet. Sci. Lett. 181 (1–2), 229–249. Brantley, S.L., Yoxtheimer, D., Arjmand, S., Grieve, P., Vidic, R., Pollak, J., Llewellyn, G.T., Abad, J., Simon, C., 2014. Water resource impacts during unconventional shale gas development: the Pennsylvania experience. Int. J. Coal Geol. 126 (November 2013), 140–156. Broadhead, R.F., 1990. Bravo Dome carbon dioxide gas field. In: AAPG Treatise of Petroleum Geology, Atlas of Oil and Gas Fields, Structural Traps I, pp. 213–232. Buckley, S.E., Leverett, M.C., 1942. Mechanism of fluid displacements in sands. Trans. AIME 146 (1337), 107–117. Cassidy, M., 2005. Occurrence and origin of free carbon dioxide gas deposits in the earth crust. PhD thesis. University of Houston. Darrah, T.H., Jackson, R.B., Vengosh, A., Warner, N.R., Poreda, R.J., 2015. Noble gases: a new technique for fugitive gas investigation in groundwater. Groundwater 53 (1), 23–28. Darrah, T.H., Vengosh, A., Jackson, R.B., Warner, N.R., Poreda, R.J., 2014. Noble gases identify the mechanisms of fugitive gas contamination in drinking-water wells overlying the Marcellus and Barnett Shales. Proc. Natl. Acad. Sci. USA 111 (39), 14076–14081. K.J. Sathaye et al. / Earth and Planetary Science Letters 450 (2016) 1–9 Dumore, J.M., Hagoort, J., Risseeuw, A.S., 1984. An analytical model for onedimensional, three-component condensing and vaporizing gas drives. Soc. Pet. Eng. J. 24 (2), 169–179. Gasda, S.E., Bachu, S., Celia, M.A., 2004. Spatial characterization of the location of potentially leaky wells penetrating a deep saline aquifer in a mature sedimentary basin. Environ. Geol. 46 (6–7), 707–720. Gilfillan, S.M.V., Ballentine, C.J., Holland, G., Blagburn, D., Sherwood-Lollar, B., Stevens, S., Schoell, M., Cassidy, M., 2008. The noble gas geochemistry of natural CO2 gas reservoirs from the Colorado Plateau and Rocky Mountain provinces, USA. Geochim. Cosmochim. Acta 72 (4), 1174–1198. Gilfillan, S.M.V., Sherwood Lollar, B., Holland, G., Blagburn, D., Stevens, S., Schoell, M., Cassidy, M., Ding, Z., Zhou, Z., Lacrampe-Couloume, G., Ballentine, C.J., 2009. Solubility trapping in formation water as dominant CO2 sink in natural gas fields. Nature 458 (7238), 614–618. Györe, D., Stuart, F.M., Gilfillan, S.M.V., Waldron, S., 2015. Tracing injected CO2 in the Cranfield enhanced oil recovery field (MS, USA) using He, Ne and Ar isotopes. Int. J. Greenh. Gas Control 42, 554–561. Hughes, J.D., 2013. Energy: a reality check on the shale revolution. Nature 494 (7437), 307–308. Hunt, A.G., Darrah, T.H., Poreda, R.J., 2012. Determining the source and genetic fingerprint of natural gases using noble gas geochemistry: a northern Appalachian basin case study. AAPG Bull. 96 (10), 1785–1811. Isaacson, E.L., Marchesin, D., Plohr, B.J., 1990. Transitional waves for conservation laws. SIAM J. Math. Anal. 21 (4), 837–866. Jackson, R.B., Vengosh, A., Darrah, T.H., Warner, N.R., Down, A., Poreda, R.J., Osborn, S.G., Zhao, K., Karr, J.D., 2013. Increased stray gas abundance in a subset of drinking water wells near Marcellus shale gas extraction. Proc. Natl. Acad. Sci. USA 110 (28), 11250–11255. Johansen, T., Wang, Y., Orr, F.M., Dindoruk, B., 2005. Four-component gas/oil displacements in one dimension: Part I: global triangular structure. Transp. Porous Media 61 (1), 59–76. Johns, R.T., Dindoruk, B., 2013. Gas flooding. In: Sheng, J. (Ed.), Enhanced Oil Recovery Field Case Studies, 1st edition. Gulf Professional Publishing, pp. 1–30 (Chapter 1). Kerr, R.A., 2010. Natural gas from shale bursts onto the scene. Science 328, 1624–1626. Laforce, T., Cinar, Y., Johns, R.T., Orr Jr., F.M., 2010. Experimental confirmation of analytical composition routes in three-phase partially miscible flow. Soc. Pet. Eng. J. (March), 160–170. Lake, L.W., 1989. Enhanced Oil Recovery. Prentice Hall. LeVeque, R.J., 1992. Numerical Methods for Conservation Laws. Birkhäuser. Leverett, M.C., 1939. Flow of oil–water mixtures through unconsolidated sands. Trans. AIME 132 (1), 149–171. Leverett, M.C., 1941. Capillary behavior in porous solids. Trans. AIME 142 (1), 152–169. Michael, K., Golab, A., Shulakova, V., Ennis-King, J., Allinson, G., Sharma, S., Aiken, T., 2010. Geological storage of CO2 in saline aquifers. A review of the experience from existing storage operations. Int. J. Greenh. Gas Control 4 (4), 659–667. 9 Moritz, A., Hélie, J.-F., Pinti, D.L., Larocque, M., Barnetche, D., Retailleau, S., Lefebvre, R., Gélinas, Y., 2015. Methane baseline concentrations and sources in shallow aquifers from the shale gas-prone region of the St. Lawrence Lowlands (Quebec, Canada). Environ. Sci. Technol. 49 (7), 4765–4771. Orr, F.M., 2007. Theory of Gas Injection Processes. Tie-Line Publications, Copenhagen. Orr, F.M., Taber, J.J., 1984. Use of carbon dioxide in enhanced oil recovery. Science (New York, N.Y.) 224 (4649), 563–569. Orr Jr., F.M., 2009. CO2 capture and storage: are we ready? Energy Environ. Sci. 2 (5), 449–458. Osborn, S.G., Vengosh, A., Warner, N.R., Jackson, R.B., 2011. Methane contamination of drinking water accompanying gas-well drilling and hydraulic fracturing. Proc. Natl. Acad. Sci. USA 108 (20), 8172–8176. Prinzhofer, A., Dos Santos Neto, E.V., Battani, A., 2010. Coupled use of carbon isotopes and noble gas isotopes in the Potiguar basin (Brazil): fluids migration and mantle influence. Mar. Pet. Geol. 27 (6), 1273–1284. Rhee, H.-R., Aris, R., Amundson, N.R., 2001. First-Order Partial Differential Equations, Volume 2: Theory and Applications of Hyperbolic Systems of Quasilinear Equations. Dover Books on Mathematics. Ricchiuto, T., Schoell, M., 1988. Origin of natural gases in the Apulian basin in south Italy: a case history of mixing of gases of deep and shallow origin. Org. Geochem. 13 (1–3), 311–318. Sathaye, K.J., Hesse, M.A., Cassidy, M., Stockli, D.F., 2014. Constraints on the magnitude and rate of CO2 dissolution at Bravo Dome natural gas field. Proc. Natl. Acad. Sci. USA 111 (43), 15332–15337. Sathaye, K.J., Smye, A.J., Jordan, J.S., Hesse, M.A., 2016. Noble gases preserve history of retentive continental crust in the Bravo Dome natural CO2 field, New Mexico. Earth Planet. Sci. Lett. 443, 32–40. Szulczewski, M.L., MacMinn, C.W., Herzog, H.J., Juanes, R., 2012. Lifetime of carbon capture and storage as a climate-change mitigation technology. Proc. Natl. Acad. Sci. USA 109 (14), 5185–5189. Trautz, R.C., Pugh, J.D., Varadharajan, C., Zheng, L., Bianchi, M., Nico, P.S., Spycher, N.F., Newell, D.L., Esposito, R.A., Wu, Y., Dafflon, B., Hubbard, S.S., Birkholzer, J.T., 2013. Effect of dissolved CO2 on a shallow groundwater system: a controlled release field experiment. Environ. Sci. Technol. 47 (1), 298–305. Vidic, R.D., Brantley, S.L., Vandenbossche, J.M., Yoxtheimer, D., Abad, J.D., 2013. Impact of shale gas development on regional water quality. Science 340 (6134), 1235009. Wachmann, C., 1964. A mathematical theory for the displacement of oil and water by alcohol. Soc. Pet. Eng. J. 231, 250–266. Warner, N.R., Jackson, R.B., Darrah, T.H., Osborn, S.G., Down, A., Zhao, K., White, A., Vengosh, A., 2012. From the cover: geochemical evidence for possible natural migration of Marcellus Formation brine to shallow aquifers in Pennsylvania. Proc. Natl. Acad. Sci. USA 109 (30), 11961–11966. Welge, H.J., Johnson, E.F., Ewing Jr., S.P., Brinkman, F.H., 1961. The linear displacement of oil from porous media by enriched gas. J. Pet. Technol. 13 (8). Whitson, C.H., Michelsen, M.L., 1989. The negative flash. Fluid Phase Equilib. 53, 51–71. Wyckoff, R.D., Botset, H.D., Ewing Jr., S.P., Brinkman, F.H., 1961. The flow of gas– liquid mixtures through unconsolidated sands. J. Appl. Phys. 7 (921), 325–345.