Next Article in Journal
Assessment of Easily Accessible Spectroscopic Techniques Coupled with Multivariate Analysis for the Qualitative Characterization and Differentiation of Earth Pigments of Various Provenance
Next Article in Special Issue
Improvement of Carbon Dioxide Sequestration of Anorthite through Bacterial: Release of Calcium and Destruction of Crystal Structure
Previous Article in Journal
The Efficiency of Black Mass Preparation by Discharge and Alkaline Leaching for LIB Recycling
Previous Article in Special Issue
Study on the Alteration of Pore Parameters of Shale with Different Natural Fractures under Supercritical Carbon Dioxide Seepage
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Effects of Charged Solute-Solvent Interaction on Reservoir Temperature during Subsurface CO2 Injection

by
Christopher Paolini
Department of Electrical and Computer Engineering, San Diego State University, San Diego, CA 92182, USA
Minerals 2022, 12(6), 752; https://doi.org/10.3390/min12060752
Submission received: 16 April 2022 / Revised: 8 June 2022 / Accepted: 10 June 2022 / Published: 14 June 2022

Abstract

:
A short-term side-effect of CO2 injection is a developing low-pH front that forms ahead of the bulk water injectant, due to differences in solute diffusivity. Observations of downhole well temperature show a reduction in aqueous-phase temperature with the arrival of a low-pH front, followed by a gradual rise in temperature upon the arrival of a high concentration of bicarbonate ion. In this work, we model aqueous-phase transient heat advection and diffusion, with the volumetric energy generation rate computed from solute-solvent interaction using the Helgeson–Kirkham–Flowers (HKF) model, which is based on the Born Solvation model, for computing specific molar heat capacity and the enthalpy of charged electrolytes. A computed injectant water temperature profile is shown to agree with the actual bottom hole sampled temperature acquired from sensors. The modeling of aqueous-phase temperature during subsurface injection simulation is important for the accurate modeling of mineral dissolution and precipitation because forward dissolution rates are governed by a temperature-dependent Arrhenius model.

1. Introduction

The first geologic carbon sequestration test in the United States was the Frio Brine Pilot experiment funded by the Department of Energy (DOE) National Energy Technology Laboratory (NETL) under the leadership of the Bureau of Economic Geology (BEG) at the Jackson School of Geosciences, The University of Texas at Austin, with major collaboration from the GEO-SEQ project, a national lab consortium led by Lawrence Berkeley National Laboratory (LBNL) [1]. In this test, 1600 metric tons of CO2 was injected into a sandstone layer residing in the Frio Formation, located east of Houston, Texas. An injection well was drilled to 1750 m KB (5755 ft KB) and perforated in the upper part of the upper Frio C sandstone region. An existing well 100 ft (30 m) updip was used as an observation well for measuring temperature, pressure, and solute concentration. Table 1 shows the depth of each region. The examination of water samples taken from the observation well over the span of one month, in which injection took place over a 10-day period, revealed a decrease in pH (from 6.5 to 5.7) occurred well before the arrival of fluid with increased alkalinity (from 100 to 3000 mg/L as HCO3 [2]. A major concern of this finding was the potential for the rapid dissolution of carbonate minerals that encounter the acidic front, which could lead to undesirable increased permeability in rock seals or well cements and the subsequent potential leakage of CO2 and brine. On the other hand, an advancing low pH front may increase the wettability of mineral surfaces, which would improve CO2 sequestration efficiency with respect to the capillarity of the moving effluent, as well as aid in enhanced oil recovery operations. Previous research by Pruess and Müller [3] showed, through 1D numerical simulation, that the injection of CO2 into saline aquifers may cause formation dry-out and the precipitation of NaCl(s) near an injection well and that such precipitation would reduce formation porosity, permeability, and injectivity, leading to a pressure build up. While Pruess and Müller’s simulations were isothermal, they concluded that nonisothermal effects could result in greater effects on reducing formation porosity, permeability, and injectivity if CO2(g) is injected at a temperature different from the formation aqueous brine phase. In subsequent research, Pruess derived [4] an analytical expression for finding the saturation of NaCl(s) (i.e., the volume fraction of solid NaCl precipitate) in the dry-out region.

1.1. Separation Distance between the Diffusive Acidic and Alkaline Bulk Effluent Fronts

1D numerical simulations were conducted to examine the distance between the developing acidic front and the alkaline front as a function of reservoir temperature, seepage velocity, and time [5]. Figure 1 shows the separation distance as a function of seepage velocity and reservoir temperature after 5 years computed using a 1D reaction transport modeling (RTM) application [6] configured for a sandstone lithology similar to the high-permeability brine-bearing sandstone of the Frio Formation. A substantial contributing factor to the development of a low pH front is the order of magnitude difference in solute diffusivity (see Table 2) between hydron H+ (1.4770 × 10−4 cm2 s−1) and bicarbonate ion HCO3 (2.1560 × 10−5 cm2 s−1), as approximated using the model of Boudreau in Equation (1) [7], with the parameters Tc and Tf representing the solutes used in a 1D reaction transport application to model the diffusion and sweep front distance given in Figure 1. A superimposed plot of the observation-well-sampled temperature and pressure measurements [8] for the month of October 2004 is shown in Figure 2. The vertical green line in the figure denotes the arrival of an introduced gas-phase fluorescent tracer to indicate CO2 breakthrough. The three pressure drop shown in the figure identifies the pump shutoff. The data show an initial, pre-injection reservoir temperature of 58.4 °C and a final post-injection temperature of 58.8 °C at the final 30-day sampling, with a remarkable drop in temperature to 57.6 °C at breakthrough. One may suspect the sharp temperature drop is a result of Joule–Thomson cooling, but the Joule–Thomson coefficient μJT has been shown to be positive for pure CO2 at several reservoir representative temperatures across a pressure range of 1 to 500 bar [9],
D = 10 6 ( T c + T f T ) ,   [ c m 2 s 1 ]
μ J T = ( T P ) H Δ T Δ P
And the corresponding drop in temperature occurs with an increase in pressure at breakthrough.

1.2. Contribution of Soulte-Solute Interaction on Temperature

In this paper we show, using numerical 3D reaction transport modeling, how the interaction of charged electrolytes with the formation water solvent produces a temperature profile in agreement with sampled wellbore measurements. The results of numerical simulation also show the arrival of the low pH diffusion front results in relative cooling, while the later arrival of the high alkalinity sweep front results in a steady temperature increase and the maintenance of a relatively heated plume that approaches 58.8 °C for 17 days post-injection up to the 30-day final measurement, due to a calculated aqueous-phase thermal diffusivity α of approximately 4 × 10−7 m2 s−1. To model the temperature change resulting from the charged solute-solvent interaction, we implemented the Helgeson–Kirkham–Flowers (HKF) model [10], which is based on the Born solvation model [11,12] that governs the aqueous-phase thermodynamics of the ion–solvent interaction, within an existing RTM application [6] that was extended to 3D with new numerical models, developed by the author, for poroelasticity, thermal transport, and porous media flow. The modeling of the aqueous-phase temperature during the simulation of the subsurface injection is important for the accurate modeling of mineral dissolution and precipitation because forward dissolution rates are governed by the temperature-dependent Arrhenius model.

1.3. Organization of This Article

The remainder of this article is organized as follows. Nomenclature provides a table of definitions of symbols used throughout this document. Section 2.1, The Born Solvation Model, provides a theoretical background of the model used to find the change in free energy due to solvation. Section 2.2, Aqueous Phase Thermal Transport, discusses the construction of the source term used in an advection–diffusion model to simulate thermal transport and a discussion of the finite element method (FEM) used to implement the model. Section 2.4, Porous Media Flow, discusses coupled models used to simulate fluid velocity and density in brine saturated sandstone. Section 2.7, Poroelastic Model, discusses coupled models used to simulate shale and sandstone strain, stress, and porosity. Section 2.9, Aqueous Phase Mass Transport, discusses the advection–diffusion reaction model to simulate solute (mass) transport and the kinetic mechanism used to model mineral dissolution and precipitation. Section 2.11, Frio Formation Simulated Configuration, presents the design of a new software application developed to simulate subsurface CO2 injection and how the application was configured to simulate CO2 injection into the Frio Formation. Section 3, Results and Discussion, presents results of simulating CO2 injection into the Frio Formation with respect to a developing temperature gradient between the acidic and alkaline fronts. Section 4, Conclusions, summarizes the major findings of this paper.

2. Methodology

2.1. The Born Solvation Model

The electric field of a charged electrolyte (ion) solute will interact with the dipoles that form a water solvent. The permittivity (dielectric constant) ε of the water solvent determines how the ion’s applied electric field affects charge migration (electric flux) and dipole reorientation in the solvent. The Born model provides an approximation of the standard state molar Gibbs free energy of solvation, Δ G ¯ s which works by establishing a polarization charge in a solvent upon the introduction of an ion. The permittivity ε is a measure of capacitance that is encountered by an ion when forming an electric field in a medium, such as the reservoir formation water, and describes the amount of charge needed to generate one unit of electric flux in a given medium. In our RTM application, calculations for the relative permittivity of H2O use the power series Equation (3) and parameters (Table 3) generated by Johnson and Norton [13],
ε r = k = 1 5 c k D w k 1 ,   where   c 1 = 1 , c 2 = a 1 T n , c 3 = a 2 T n + a 3 + a 4 T n , c 4 = a 5 T n + a 6 T n + a 7 T n 2 , c 5 = a 8 T n 2 + a 9 T n + a 10
An ion’s charge will yield more electric flux in a solvent with low permittivity than in a solvent with high permittivity. From Equation (3), the relative permittivity of H2O decreases with increasing temperature. At lower temperatures, H2O dipoles can become more aligned and establish a polarization charge, while at higher temperatures the dipoles are less able to align due to random molecular motion from increased kinetic energy, which results in more flux or charge migration through the solvent [14]. From Coulomb’s Law, the force on a test charge Q due to a single point charge q (at rest) a distance r apart is
F = 1 4 π ε r ε 0 q Q r 2 r ^ , N   m 2 C 2 C 2 m 2 = N r
From (4), the electric potential from a point charge q at a distance r from q is given by
ϕ ( r ) = q 4 π ε r ε 0 r , N   m 2 C 2 C m = N   m C = J C = V
And from (5), it follows that the electric potential from a uniform set of Z (i.e., valence) elementary charges e (coulombs) at a distance r from Ze is
ϕ ( r ) = Z e 4 π ε r ε 0 r ,   V
The energy U of an electric field E from the successive introduction of point charges in a medium with absolute permittivity ε can be found by first assuming the energy from introducing the first point charge q1 into an empty medium will give Uq1 = 0. From (5), the energies Uq1, q2, Uq1, q2, q3, and Uq1, q2, q3,…, qn from introducing the second, third, and nth point charges, respectively, are
U q 2 = 1 4 π ε q 2 q 1 r 1 , 2 , U q 3 = 1 4 π ε q 3 ( q 1 r 1 , 3 + q 2 r 2 , 3 ) , , U q n = 1 4 π ε 1 2 i = 1 n j i n q i q j r i , j
In the last term of (7), a division by 2 is required since the product charges qi qj and qj qi are each counted [15]. Factoring the potential function (5) from (7), we have
U q n = 1 2 i = 1 n q i ( j i n 1 4 π ε q j r i , j ) = 1 2 i = 1 n q i ϕ ( r ¯ i )
For a charge density ρ, integrating over an arbitrary volume V with permittivity ε, we have
U = 1 2 V ρ ϕ d V = 1 2 V ( ε E ) ϕ d V = ε 2 [ V E ϕ d V + S ϕ E n ^ d S ]
where Gauss’s law ε ( · E ) = ρ is employed in the second term and the product rule · ( ϕ E ) = ϕ ( · E ) + E · ϕ and divergence theorem are employed in the third term. From the fundamental theorem of calculus for line integrals, ∇φ = −E, and observing as r → where the electric potential at infinity is assumed to be zero [15],
S = 4 π r 2 r 2 ,   E = 1 4 π ε q r 2 r ^ 1 r 2 ,   ϕ ( r ) = q 4 π ε r ε 0 r 1 r S ϕ E n d S 1 r 0
Equation (9) reduces to
U = ε 2 V | E | 2 d V = ε 0 ε r 2 V | 1 4 π ε 0 ε r q r 2 r ^ | 2 d V = ε 0 ε r 2 q 2 16 π 2 ( ε 0 ε r ) 2 V 1 r 4 d V
For the energy of uniformly a charged spherical shell (ion) of total charge q and radius r. Substituting dV = dA dr = r2sinθdθdφ dr and q = Ze into (11) and integrating over the volume outside of the ion’s surface from r = ri to r = , where ri is the radius of the ion, we obtain the Born expression [12] for the energy of an ion with charge Ze in a medium of permittivity ε,
1 2 ε r ε 0 Z 2 e 2 16 π 2 ( ε r ε 0 ) 2 V r 2 1 r 6 d r d A = 1 2 ε r ε 0 Z 2 e 2 16 π 2 ( ε r ε 0 ) 2 V r 2 1 r 6 d r r 2 s i n   θ d θ d ϕ = 1 2 ε r ε 0 Z 2 e 2 16 π 2 ( ε r ε 0 ) 2 4 π r i 1 r 2 d r = Z 2 e 2 8 π ε [ 1 1 r i ] = Z 2 e 2 8 π ε r i = k e Z 2 e 2 2 ε r i
Denoting U0 as the energy of the electric field of an ion in a vacuum and U as the energy of the electric field of an ion in a medium (solvent) of relative permittivity εr, the energy difference U-U0 provides the work done introducing an ion into a solvent, or equivalently, the Gibbs free energy of solvation,
Δ G s = U U 0 = k e Z 2 e 2 2 ε r r i k e Z 2 e 2 2 r i = k e Z 2 e 2 2 r i ( 1 ε r 1 )
Multiplying (13) by the Avogadro constant NA and defining the absolute Born coefficient ω j a b s , the absolute standard molal Gibbs free energy of solvation of the jth ion in a solvent is given by [14,16]
Δ G ¯ s , j , a b s = N A Z j 2 e 2 8 π ε 0 r j ( 1 ε r 1 ) = ω j a b s ( 1 ε r 1 )
Substituting the effective ionic radius, r e , j = r j , T r , P r + | Z j | g ( T , P ) , defined by Tanger and Helgeson [17], where rj is the crystallographic ionic radius (~0.3Å to 2Å) and g is a solvent function at reference temperature and pressure with units of Å obtained by regression of and heat capacity and volume data for NaClaq [14,17,18]. To derive the thermodynamic properties, such as specific heat capacity and specific enthalpy, of charged solutes, which are needed in our thermal transport model, Tanger and Helgeson define a conventional electrostatic Born coefficient ωj. The conventional coefficient is the difference between the absolute Born coefficient for ion j, which represents Δ G ¯ s , j °   of ion j, and Δ G ¯ s , h s ° which is the Gibbs free energy required to form an inner hydration sphere of H2O dipoles encompassing ion j. The hydration sphere Gibbs free energy is the valence of ion j, Zj,, multiplied by Δ G ¯ s , H + °   , the solvation energy for a proton. The conventional Born coefficient defined by Tanger and Helgeson is thus ω j = ω j a b s Z j ω H + a b s , and the relative standard state molar Gibbs free energy of solvation of ion j becomes
Δ G ¯ s , j = ω j ( 1 ε r 1 )
For εr > 1 in (15), Δ G ¯ s , j < 0, which indicates the movement of an ion from a less polar solvent, such as a vacuum, to H2O, is thermodynamically favorable. In addition, the magnitude of Δ G ¯ s , j °   decreases as the ionic radius rj increases since the voltage φ(r) at the ion’s outer shell, defined by Equation (6), decreases with the distance from the ion’s center where all charge is assumed (by the Born model) to reside [19]. From (15), the standard state molar entropy of solvation is obtained through the relation
( Δ G ¯ s T ) P = Δ S ¯ s = ( ω j ( 1 ε r ( T ) 1 ) ) / T | P = ω j ( 1 ε r 2 ( T ) ( ε r T ) P ) ( 1 ε r ( T ) 1 ) ( ω j T ) P
The standard state molar enthalpy of solvation is then found by substituting (16) into Δ G ¯ s ° = Δ H ¯ s ° T Δ S ¯ s ° ,
Δ H ¯ s = Δ G ¯ s + T Δ S ¯ s = ω ( 1 ε r 1 ) + ω T Y T ( 1 ε r 1 ) ( ω T ) P
where the partial derivatives of relative permittivity εr at constant pressure are denoted as functions Z, Y, and X [16],
Z = ( 1 ε r ) , Y = ( Z T ) P = ( 1 ε r 2 ) ( ε r T ) P [ 1 K ] , X = ( Y T ) P = ( 1 ε r 2 ) ( 2 ε r T 2 ) P ( 2 ε r Y 2 ) [ 1 K 2 ]
For Δ H ¯ s ° < 0, the ion–solvent interaction is exothermic, and there will be a temperature increase in the formation water. The first and second partial derivatives of relative permittivity εr from Equation (3) at constant pressure are given by
( ε r T ) P = j = 0 4 D w j [ ( d c d T ) j + 1 j α c j + 1 ] ( 2 ε r T 2 ) P = j = 0 4 D w j [ ( d 2 c d T 2 ) j + 1 j ( α ( d c d T ) j + 1 + c j + 1 d α d T ) j α ( ( d c d T ) j + 1 j α c j + 1 ) ]
The first and second partial derivatives of the conventional Born coefficient ωj are given by [16].
( ω T ) = η X 1 ( g T ) p [ J mol   K ] ,   X 1 = a b s ( z 3 ) ( r i ) 2 [ z ( 3.082 + g ) 2 ] ( 2 ω T 2 ) = 2 η X 2 ( g T ) 2 p η X 1 ( 2 g T 2 ) p [ J mol   K 2 ] ,   X 2 = a b s ( z 4 ) ( r i ) 3 [ z ( 3.082 + g ) 3 ]
From (17), H.C. Helgeson, D.H. Kirkham, and G.C. Flowers (Helgeson et al., 1981), derived a revised model for specific molar heat capacity c ¯ P
c ¯ P = c 1 + c 2 ( T Θ ) 2 2 T ( T Θ ) 3 [ a 3 ( P P r ) + a 4 ln ( Ψ + P Ψ + P r ) ] nonsolvation   contribution + ω T X + 2 T Y ( ω T ) P T ( 1 ε r 1 ) ( 2 ω T 2 ) P solvation   contribution
and specific molar enthalpy h ¯ [14,16]
h ¯ = c 1 ( T T r ) c 2 [ ( 1 T Θ ) ( 1 T r Θ ) ] + a 1 ( P P r ) + a 2 ln ( Ψ P Ψ P r ) + ( 2 T Θ ( T Θ ) 2 ) [ a 3 ( P P r ) + a 4 ln ( Ψ + P Ψ + P r ) ] + ω ( 1 ε 1 ) + ω T Y T ( 1 ε 1 ) ( ω T ) P ω Pr , T r ( 1 ε Pr , T r 1 ) ω Pr , T r T r Y Pr , T r
For an aqueous electrolyte species j, in (22), coefficients ai and ci are species dependent non-solvent regression parameters ([20], p. 92).

2.2. Aqueous Phase Thermal Transport

The aqueous-phase volumetric energy generation rate ST can be derived by summing, over all solutes j for n solutes, the rate of change of molar concentration of species j multiplied by that species’ specific molar enthalpy given by (22),
S T = j = 1 n d c j d t h ¯ j [ J m 3 s ]
The transient heat advection–diffusion reaction model is then
( ρ c p T ) t + ( ρ c p T u ) = ( k T ) + S T
where u is fluid velocity, k = k(T) is the thermal conductivity of H2O, and ρ = ρ(T) is the density of H2O.

2.3. Transient Finite Element Method Formulation

A Crank–Nicolson discretization scheme is used with time step Δt and implemented using the deal.II general-purpose, object-oriented, finite element library [21]. A weak formulation with weight function v and solution wn at time tn is given by
1 Δ t m ( w n + 1 , v ) + 1 2 a ( t n + 1 , w n + 1 , v ) = 1 Δ t m ( w n , v ) + 1 2 a ( t n , w n , v ) + 1 2 ( F ( t n ; v ) + F ( t n + 1 ; v ) )
Functionals with thermal diffusivity α and a solute-solvent reaction thermal rate ST/ρcp with units K s−1 are
m ( w , v ) = Ω w v d V   a ( t n , w , v ) = Ω α w n v + u ( t n ) ( w n ) v + u ( t n ) w n v d V   F ( t , v ) = Ω ( S T / ρ c p ) v d V
For Nh node points, define the function u from Lagrange basis polynomials φ
u = j = 1 N h u j ϕ j
Let UnRNh be the vector of coefficients uj for the finite element solution at time step n. The aqueous-phase temperature at the n + 1th time step is solved via the implicit system
( M + Δ t 2 A ( t n + 1 ) ) U n + 1 = ( M + Δ t 2 A ( t n ) ) U n + Δ t 2 ( f n + 1 + f n )
where
M i j = m ( ϕ i , ϕ j )   A i j ( t n ) = a ( t n ; ϕ i , ϕ j ) f n ( i ) = F ( t n ; ϕ i )

2.4. Porous Media Flow

A Darcy Law formulation for fluid flux is implemented to model aqueous-phase fluid injection,
v res = K μ ρ 0 ( p res ρ r e s g z )
Equation (30) is solved simultaneously with the continuity Equation [22]
t ( ρ 0 ϕ * ) + v res = q res

2.5. Mixed Finite Element Method Formulation

We use a mixed finite element method to simultaneously calculate the reservoir fluid velocity vector field vres, together with a scalar reservoir pressure field p and test functions ψ and χ, respectively. A weak formulation is constructed by multiplying (30) by a vector valued function ψ and (31) by a scalar valued function χ and then integrating the domain Ω,
Ω ψ T ( v r e s + K ρ 0 μ ( p r e s ρ r e s g z ^ ) ) d V = 0
Ω χ ( t ( ρ 0 ϕ * ) + v r e s q r e s ) d V = 0
For a given pair of test functions, (ψ, χ), (32) and (33) states that the weighted residual (error) of PDEs (30) and (31), when applied to solution functions vres and p, will be zero. Equations (32) and (33) are a weak form of (30) and (31), respectively, because these weak-form equations may also be true for other test functions that would then approximate the so-called strong solutions that solve the original, strong form of the PDEs. The weak form equations can be rewritten with additive terms separated,
( ψ , v r e s ) Ω + ( ψ , K ρ 0 μ p r e s ) Ω ( ψ , K ρ 0 μ ρ r e s g z ^ ) Ω = 0
( χ , t ( ρ 0 ϕ * ) ) Ω + ( χ , v r e s ) Ω ( χ , q r e s ) Ω = 0
The pressure gradient term ∇pres in (34) can be eliminated by applying the product rule and the Divergence Theorem,
( p r e s ψ ) = p r e s ( ψ ) + ψ p r e s
Applying the Divergence Theorem introduces an integral term in terms of the divergence of test function ψ, and a boundary integral in terms of the pressure values and the normal component of the test function, ψ · n. Normalizing (34) by multiplying each term by μ K−1 ρ0−1 gives
( ψ , μ K 1 ρ 0 v r e s ) Ω + ψ n , p r e s Ω ( ψ , p r e s ) Ω ( ψ , ρ r e s g z ^ ) Ω = 0
To obtain a unique solution, a constraint is placed on the pressure function. In our model, the far field pressure on the boundary furthest from the injection location and above the caprock layer is given a Dirichlet boundary condition set to equal the initial pressure. The assumption is that the far field pressure is beyond the range of influence of the fluid injection. We call the portion of the boundary where this condition is set ΓD, and we split the boundary integral as
ψ n , p r e s Ω = ψ n , p r e s Ω Γ d ψ n , p 0 Γ d
Separating the known terms from the unknown terms in Equation (37) and grouping the known terms on the right-hand side, we obtain
( ψ , μ K 1 ρ 0 v r e s ) Ω + ψ n , p r e s Ω Γ d ( ψ , p r e s ) Ω = ( ψ , ρ r e s g z ^ ) Ω + ψ n , p 0 Γ d , ψ V h
( ϕ * ρ 0 n , χ ) δ t ( v r e s , χ ) + ( ϕ * ρ 0 n 1 + δ t q r e s , χ ) = 0 , χ W h
where in (39) and (40), we find fields (vres, p) ∈Vh × Wh where Vh and Wh are vector and scalar components of a mixed finite element space on a quadrilateral/hexahedral mesh of the domain Ω. A mixed formulation in terms of two variables, vres and p, is locally conservative. The inner product (v, w) is defined as
( v , w ) Ω = Ω w T v   d V
Note that definition (41) applies to scalar-valued functions by noting that the transpose operator is the identity for scalars. Raviart–Thomas elements are used for fluid flux and discontinuous Lagrange polynomials for pressure. This choice of basis elements satisfies the La-dyshenskaya, Babuska, and Brezzi (LBB) condition for the mixed Laplace formulation given by (41). Rothe’s method of first discretizing in time and then in space is used, with a theta method for the time discretization. We choose Rothe’s method because our RTM application utilizes the coupling of the physical and chemical domains that can lead to changes in flow-related domain parameters and boundary conditions over time. For instance, changes in permeability are incorporated by changing the permeability tensor on a time-step by time-step basis without any changes to the computational framework. The theta method is a semi-implicit method that allows for adjusting the discretization’s dependence on the forward and backward timestep through parameter θ ∈ [0, 1]. When θ is chosen as 0.0, 0.5, or 1.0, the method gives the forward Euler, Crank–Nicolson, or backward Euler methods, respectively. We use the Crank–Nicolson method, though slight deviations from θ = 0.5 can introduce some numerical diffusion, if needed, to smooth out oscillations in time, while sacrificing the second-order accuracy of the Crank–Nicolson method. Designating the chosen time step as Δt, the implementation of Rothe’s method with a theta scheme for time discretization on Equation (40) is
( χ , ρ 0 ϕ * ) Ω n ( χ , ρ 0 ϕ * ) Ω n 1 + θ Δ t ( χ , v r e s ) Ω n + ( 1 θ ) Δ t ( χ , v r e s ) Ω n 1 θ Δ t ( χ , q r e s ) Ω n ( 1 θ ) Δ t ( χ , q r e s ) Ω n 1 = 0

2.6. Fluid Density Equation of State

Fluid density is modeled with a slightly compressible linearized equation of state [22]
ρ r e s ( p ) = ρ 0 ( 1 + c p )
Expanding (43) and expressing in terms of inner products, we obtain
( χ , ρ r e s ϕ * ) Ω n = ( χ , ρ 0 ( 1 + c p ) ϕ * ) Ω n = ( χ , ρ 0 ϕ * ) Ω n + ( χ , ρ 0 c p ϕ * ) Ω n

2.7. Poroelastic Model

The simultaneous solution of reservoir rock strain ϵ(u), stress σ(u), effective stress σpor(u,p), and displacement u is solved using a finite element formulation according to the system,
ϵ ( u ) = 1 2 ( u + u T ) σ ( u ) = λ ( u ) I + 2 G ϵ ( u ) σ p o r ( u , p ) = σ ( u ) α B W p I σ p o r ( u , p ) = f
In (45), rock strain ϵ(u) is the symmetric gradient of displacement. Rock stress σ(u) accounts for internal rock stresses not considering the external pressure field p or body force f. Under a quasi-static assumption, the body force due to gravity, f, of the rock mass is balanced by the divergence of the reservoir stress rank−2 tensor σpor(u,p).

2.8. Reservoir Porosity Model

Effective porosity ϕ in the reservoir is modelled as a function of the rock strain ϵ(u) and pore pressure p [22],
ϕ * = ϕ 0 + α B W ( ϵ v ϵ v 0 ) + 1 M B ( p p 0 )

2.9. Aqueous Phase Mass Transport

An advection–diffusion reaction mass transport model (Park, March 2014) that conserves mass at the elemental scale is used to model solute transport in the pore water,
e β t = ϕ j = 1 N j ν β j [ ( D j c j ) ( c j u ) ] γ = 1 M γ ν β γ A γ G γ
The evolution of chemical elemental mass depends on mass-transfer from diffusive and advective forces as well as the precipitation and dissolution of minerals governed by kinetic reaction rates. The source term in (47) accounts for the net rate of the increase or decrease of a mineral in a fluid element due to chemical kinetics.

2.10. Mineral Kinetics

The kinetic mechanism governing mineral dissolution and precipitation is given by a mineral reaction rate model for each mineral γ,
G γ = k γ , d i s s ( j = 1 , ν j γ < 0 N j c j ν j γ 1 K γ j = 1 , ν j γ > 0 N j c j ν j γ )
where the forward reaction rate for mineral dissolution is approximated using an Arrhenius model,
k γ , d i s s = A γ exp ( E a R T )
Equilibrium reactions for minerals k-feldspar, calcite, dolomite, quartz, halite, and illite, respectively, which are modeled in the simulated Frio Formation, are given by,
K A l S i 3 O 8 + 2 H 2 O   K + + O H + 3 S i O 2 , a q + A l ( O H ) 3 , a q
C a C O 3   C a 2 + + C O 3 2
C a M g ( C O 3 ) 2   C a 2 + + M g 2 + + 2 C O 3 2
S i O 2 ( s ) S i O 2 ( a q )
N a C l ( s ) N a + + C l
( 0.01 C a , 0.13 N a , 0.53 K ) ( 1.27 A l , 0.44 M g , 0.36 F e ) 2 ( 3.55 S i , 0.45 A l ) 4 O 10 [ ( O H ) 2 , ( H 2 O ) ] + 12 H 2 O + 22 H + 0.01 C a 2 + + 0.13 N a + + 0.53 K + + 1.27 A l ( O H ) 3 , a q + 0.36 F e 3 + + 0.44 M g 2 + + 3.55 S i O 2 , a q

2.11. Frio Formation Simulated Configuration

2.11.1. Computational Domain

The Frio Formation physical domain is modeled as a 22-m-thick slab from z = −1517.00 m to z = −1539.00 m with a 500 m2 area. Figure 3 shows a plot in ParaView of the computational mesh. The y-z plane at x = 0 is the injection plane where a perforation zone from z = −1527.50 m to z = −1533.75 m models an injection well centered at y = 250 m. The injection wellbore diameter was configured to be 5.5 inches. A Neumann boundary condition is imposed on the y-z plane at x = 0, such that the temperature gradient is 0 to implement a symmetry assumption on that wall. All other walls have a Dirichlet boundary condition of T = Ti imposed to implement a far-field constant reservoir temperature assumption. The initial reservoir temperature was configured to be 58.4 °C, based on the observation well measurement data shown in Figure 2.
The overburden saturated rock density was configured to ρS = 2260 kg/m3 [23]. The far field pressure on the boundary furthest from the injection perforation zone and above the shale (caprock) layer is given a Dirichlet boundary condition set equal to the initial reservoir pressure Pi = 148 bar. A homogeneous Neumann boundary condition ( P n = 0 ) is imposed for pressure along the plane of symmetry.
The reservoir water velocity is configured to be initially quiescent. At simulation time t0, the injectant mass flow rate is defined to be 3 kg/s, which remains until t = 7 days have elapsed, after which the injectant mass flow rate is defined to be 0 kg/s to simulate pump shutoff. From t = 7 days to t = 30 days, the water–rock interaction in the reservoir chemically equilibrates.
Multiple grid configurations consisting of an initial control-volume count of 100 × 100 × 16, 100 × 100 × 32, 200 × 200 × 16, 40 × 40 × 16, 40 × 40 × 32, 40 × 40 × 64, 80 × 80 × 16, 20 × 20 × 08, 20 × 20 × 16, and 20 × 20 × 32 were used. Simulations were executed with deal.II adaptive mesh refinement (AMR) enabled.

2.11.2. Lithologic Configuration

The Frio Formation lithology was modeled as a five-layer shale-sandstone system according to the depths and layers given in Table 1. The initial time t = t0 compositions of the shale and sandstone are given in Table 4 in terms of solid-phase volume fraction and grain radii. Minerals with an initial volume fraction of 0 are minerals not present at the start of the simulation but are allowed to form through precipitation according to the kinetic mechanism specified by Equations (50) through (55). Poroelastic rock properties defined in the simulation are given in Table 5. The initial and boundary concentrations of solutes residing in the pore space formation water are given in Table 6.

3. Results and Discussion

Figure 4 shows a contour plot of pH and [ H C O 3 ] on a slice along the perforation zone after 4 days of injection at a mass rate of 3 kg/s, with a background color mapped to temperature. The bicarbonate ion concentration front identifies the bulk effluent, indicated by thick purple to gold contour lines (a.k.a. “plasma colormap”). The initial reservoir water has a much lower bicarbonate concentration of 0.002 M. From the figure, one can see an acidic front develops ahead of the effluent indicated by black to white contour lines. A lower temperature fluid region develops in the lower pH (acidic) region, indicated by a darker blue background. In our simulation, we see the same pH drop before the arrival of a high concentration of bicarbonate ion 30 m from the injection zone.
Figure 5 shows snapshots of a simulation where an aqueous solution with CO2,aq mass flow rate of 3 kg s−1 is injected for 7 days into the lower sandstone. The upper row shows aqueous-phase temperature using a colormap in the range 58.4 °C (blue) to 60.0 °C (red) at selected times in units of days. The bottom row shows the source term volumetric energy generation rate ST computed using Equation (23) divided by volumetric heat capacity ρcp to give a temperature rate in units K s−1. As expected, ST vanishes (=0) after 7 days when injection terminates.

4. Conclusions

A short-term side-effect of subsurface CO2 injection in brine-bearing sandstones is a developing low-pH front that forms ahead of the bulk water injectant, due to differences in solute diffusivity. The rapid dissolution of CaCO3(s) and other carbonate minerals resulting from low-pH reservoir water can have significant environmental implications in geologic carbon sequestration. The dissolution of carbonate minerals can create unwanted pathways that will increase permeability and lead to a possible loss in CO2 containment. Additionally, the formation of H2CO3 can be vulnerable to wellbore integrity through the corrosion of structures made from calcium-silica-hydrate compounds in Portland cements used to seal the annular volume between a casing and borehole wall. Temperature, pressure, and solute concentration measurements taken from an observation well during the Frio Pilot CO2 sequestration experiment showed a reduction in aqueous-phase temperature with the arrival of the low-pH front, followed by a gradual rise in temperature upon the arrival of a high concentration of bicarbonate ion. We present a possible contributing factor to the variation in temperature between these two fronts by modeling aqueous-phase transient heat advection and diffusion, with the volumetric energy generation rate computed from solute-solvent interaction using the Helgeson–Kirkham–Flowers (HKF) model, which is based on the Born solvation model, for computing specific molar heat capacity and the enthalpy of charged electrolytes. We show how a simulated injectant water temperature profile is shown to agree with sampled temperature time-series data acquired from bottom hole observation well sensors. From Figure 4, we see our simulations reveal a lower temperature fluid region develop within an acidic region ahead of the bulk effluent. This acidic region is a diffusion front that is lower in temperature than the surrounding formation water being displaced and followed by a higher temperature sweep (effluent or injectant) front. The figure shows a plot of temperature vs. time at a single point on the slice, 30 m from the injection source, corresponding to the experimental observation well in the Frio experiment. From this plot, we see agreement between the simulated and experimental temperature profile.
The computation of aqueous-phase temperature during subsurface injection simulation is important for the accurate modeling of mineral kinetic mechanisms, as forward dissolution rates are governed by an Arrhenius model. In conclusion, the contribution of aqueous electrolytic reactions could play a role in the experimentally observed drop in temperature that corresponds to the arrival of a low-pH fluid region that forms ahead of the bulk effluent injectant, identified by a high concentration of bicarbonate ion ( H C O 3 ) that forms from dissociated aqueous CO2.

Funding

This work was supported by a grant from the US Department of Energy (DOE), National Energy Technology Laboratory (NETL), titled "Recovery Act: Web-based CO2 Subsurface Modeling", DOE Award Number DE-FE0002069, Geologic Sequestration Training and Research Funding Opportunity Number DE-FOA-0000032, Simulation and Risk Assessment.

Conflicts of Interest

The author declares no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Nomenclature

Aγmineral specific surface area, cm2 cm−3 (i.e., cm−1)
Djdiffusion coefficient of molecular species j, cm2 s−1
Dwdimensionless density of water = ρw g−1 cm3
Eelectric (vector) field, N C−1 = V·m−1
EaArrhenius model activation energy, J mol−1
Fforce (vector) field, Nr
Kpermeability of the porous medium, cm2
Kγequilibrium constant for mineral γ dissolution reaction
KBmeasured depth below the Kelly bushing
Gsecond Lame parameter (shear modulus), gram-force cm−2
Gγmineral reaction rate, mol cm−2 s−1
Mmolar concentration, mol L−1
MBBiot modulus, gram-force cm−2
Mwjmolecular weight of species j, g mol−1
Qtest charge, C (Coulombs)
Paqueous-phase pressure, bar
Piinitial reservoir aqueous-phase pressure, bar
Praqueous-phase reference pressure, 1 bar
Rgas constant, 8.31446261815324 J mol−1 K−1
STaqueous-phase volumetric energy generation rate, J m−3 s−1
THKF model aqueous-phase temperature, K
Tireservoir aqueous-phase initial temperature, 58.40 °C
TrHKF model aqueous-phase reference temperature, 228K
Uenergy of an electric field E, J (Joules)
Zion valence, integer
cjconcentration of species j, mol L−1
cpfluid compressibility constant, 5.8 × 10−10 Pa−1
C P , J specific heat capacity of ion j, J g−1 K−1
C ¯ P , J specific molar heat capacity of ion j, J mol−1 K−1
eelementary charge, 1.60217662 × 10−19 C
eβconcentration of atoms of element β, mol cm−3
fbody force, gram-force
ggravitational acceleration, cm s−2
h j specific enthalpy of ion j, J g−1
h ¯ j specific molar enthalpy of ion j, J mol−1
jaqueous-phase molecular species index
kthermal conductivity of H2O, W m−1 K−1
keCoulomb constant, 8.9875517873681764 ×109 N m2 C−2
kγ,dissforward reaction rate for mineral dissolution
qpoint charge, C
qresreservoir fluid source density rate, g cm−3 s−1
preservoir fluid pressure, gram-force cm−2
p0initial reservoir fluid pressure, gram-force cm−2
urock displacement or fluid velocity vector field, cm s−1
vresinjectant fluid mass flux into the reservoir, g cm−2 s−1
zgravitational direction unit vector, cm
Δ G ¯ s , h s ° standard state molar Gibbs free energy of formation of an inner hydration sphere, J mol−1
Δ G ¯ s , j ° standard state molar Gibbs free energy of solvation of ion j, J mol−1
Δ G ¯ s ° standard state molar Gibbs free energy of solvation, J mol−1
Δ S ¯ s , j ° standard state molar entropy of solvation, J mol−1 K−1
αaqueous-phase thermal diffusivity, m2 s−1
αBWBiot–Willis coefficient, dimensionless
βelement index
εrrelative permittivity of a medium (e.g., H2O), dimensionless = ε εo−1
εabsolute permittivity of a medium, F m−1
εopermittivity of a vacuum, 8.85 × 10−12 F m−1
ϵVreservoir volumetric strain, dimensionless
ϵ v 0 initial reservoir volumetric strain, dimensionless
γmineral index
λfirst Lame parameter, gram-force cm−2
φreservoir porosity, dimensionless
φ*effective reservoir porosity, dimensionless
φ(r)electric potential at a distance r, V (Volts)
ρdensity of water, g m−3
ρ0uncompressed density of reservoir fluid, g cm−3
ρrescompressed density of reservoir fluid, g cm−3
ρsoverburden saturated rock density, kg m−3
μdynamic viscosity of reservoir fluid, g cm−1 s−1
μJTJoule–Thomson coefficient, dimensionless
νβjelemental stoichiometric conversion factor equal to the number of atoms of element β in molecular solute species j
νβγelemental stoichiometric conversion factor equal to the number of atoms of element β in mineral γ
νmolecular stoichiometric conversion factor equal to the number of molecules of species j in mineral γ
ωjconventional electrostatic Born coefficient ωj
ω j a b s absolute standard molal Gibbs free energy of solvation of the jth ion, J mol−1
ω H + a b s absolute standard molal Gibbs free energy of the H+ ion, 0.5387 J mol−1
ΓDDirichlet boundary condition
ΓNNeumann boundary condition
ΨH2O solvent singularity pressure, 2600 bar
ΘH2O singularity temperature, 228 K

References

  1. Hovorka, S.D.; Sakurai, S.; Kharaka, Y.K.; Nance, H.S.; Doughty, C.; Benson, S.M.; Freifeld, B.M.; Trautz, R.C.; Phelps, T.; Daley, T.M. Monitoring CO2 storage in brine formations: Lessons learned from the Frio field test one year post injection. In Proceedings of the 2006 UIC Conference of the Groundwater Protection Council, San Antonio, TX, USA, 23–25 January 2006; Gulf Coast Carbon Center: Washington, DC, USA, 2006. [Google Scholar]
  2. Kharaka, Y.K.; Cole, D.R.; Hovorka, S.D.; Gunter, W.D.; Knauss, K.G.; Freifeld, B.M. Gas-water-rock interactions in Frio Formation following CO2 injection: Implications for the storage of greenhouse gases in sedimentary basins. Geology 2006, 34, 577–580. [Google Scholar] [CrossRef]
  3. Pruess, K.; Müller, N. Formation dry-out from CO2 injection into saline aquifers: 1. Effects of solids precipitation and their mitigation. Water Resour. Res. 2009, 45, W03402. [Google Scholar] [CrossRef] [Green Version]
  4. Pruess, K. Formation dry-out from CO2 injection into saline aquifers: 2. Analytical model for salt precipitation. Water Resour. Res. 2009, 45, W03403. [Google Scholar] [CrossRef] [Green Version]
  5. Paolini, C.; Park, A.; Castillo, J. An investigation of the variation in the sweep and diffusion front displacement as a function of reservoir temperature and seepage velocity with implications in CO2 sequestration. In Proceedings of the 9th Annual International Energy Conversion Engineering Conference, San Diego, CA, USA, 31 July–3 August 2011; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2011. [Google Scholar]
  6. Park, A.J. Water-rock interaction and reactive-transport modeling using elemental mass-balance approach: I. The methodology. Am. J. Sci. 2014, 314, 785–804. [Google Scholar] [CrossRef]
  7. Boudreau, B.P. Diagenetic Models and Their Implementation: Modelling Transport and Reactions in Aquatic Sediments; Springer: Berlin/Heidelberg, Germany, 1997. [Google Scholar]
  8. Frio Brine Project Data, Frio Brine Pilot Study, Gulf Coast Carbon Center, Bureau of Economic Geology, Jackson School of Geosciences; The University of Texas at Austin, Regional Carbon Sequestration Partnership: Austin, TX, USA, 2018. Available online: https://edx.netl.doe.gov/dataset/frio-brine-project-data (accessed on 15 April 2022). [CrossRef]
  9. Carneiro, J.N.E.; Pasqualette, M.A.; Reyes, J.F.R.; Krogh, E.; Johansen, S.T.; Ciambelli, J.R.P.; Rodrigues, H.T.; Fonseca, R. Numerical Simulations of High CO2 Content Flows in Production Wells, Flowlines and Risers. In Proceedings of the OTC Brasil. Offshore Technology Conference, Rio de Janeiro, Brazil, 4–7 May 2015; p. 12. [Google Scholar]
  10. Helgeson, H.C.; Kirkham, D.H.; Flowers, G.C. Theoretical prediction of the thermodynamic behavior of aqueous electrolytes by high pressures and temperatures; IV, Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600 degrees C and 5kb. Am. J. Sci. 1981, 281, 1249–1516. [Google Scholar]
  11. Bjerrum, N. Neuere Anschauungen über Elektrolyte. Ber. Der Dtsch. Chem. Ges. A B Ser. 1929, 62, 1091–1103. [Google Scholar] [CrossRef]
  12. Born, M. Volumen und Hydratationswärme der Ionen. Z. Für Phys. 1920, 1, 45–48. [Google Scholar] [CrossRef]
  13. Johnson, J.W.; Norton, D. Critical phenomena in hydrothermal systems: State, thermodynamic, electrostatic, and transport properties of H2O in the critical region. Am. J. Sci. 1991, 291, 541–648. [Google Scholar] [CrossRef]
  14. Anderson, G.M. Thermodynamics of Natural Systems, 2nd ed.; Cambridge University Press: Cambridge, UK, 2005. [Google Scholar]
  15. Griffiths, D.J. Introduction to Electrodynamics, 4th ed.; Cambridge University Press: Cambridge, UK, 2014. [Google Scholar]
  16. Marini, L. Geological Sequestration of Carbon Dioxide, Volume 11 Thermodynamics, Kinetics, and Reaction Path Modeling, 1st ed.; Elsevier Science: Amsterdam, The Netherlands, 2006. [Google Scholar]
  17. Tanger, J.C.; Helgeson, H.C. Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures; revised equations of state for the standard partial molal properties of ions and electrolytes. Am. J. Sci. 1988, 288, 19–98. [Google Scholar] [CrossRef]
  18. Shock, E.L.; Oelkers, E.H.; Johnson, J.W.; Sverjensky, D.A.; Helgeson, H.C. Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures. Effective electrostatic radii, dissociation constants and standard partial molal properties to 1000 °C and 5 kbar. J. Chem. Soc. Faraday Trans. 1992, 88, 803–826. [Google Scholar] [CrossRef]
  19. Cherie. Thermodynamic Treatment of Ion-Solvent-Interactions: The Born Model. 2014. Available online: https://www.scribd.com/document/184980407/Born-Solvation-from-Cherie (accessed on 15 April 2022).
  20. Johnson, J.W.; Oelkers, E.H.; Helgeson, H.C. SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar and 0 to 1000 °C. Comput. Geosci. 1992, 18, 899–947. [Google Scholar] [CrossRef]
  21. Bangerth, W.; Hartmann, R.; Kanschat, G. deal. II—A general-purpose object-oriented finite element library. ACM Trans. Math. Softw. 2007, 33, 24. [Google Scholar] [CrossRef]
  22. Ganis, B.; Mear, M.E.; Sakhaee-Pour, A.; Wheeler, M.F.; Wick, T. Modeling fluid injection in fractures with a reservoir simulator coupled to a boundary element method. Comput. Geosci. 2014, 18, 613–624. [Google Scholar] [CrossRef]
  23. Huang, Z.Q.; Winterfeld, P.H.; Xiong, Y.; Wu, Y.S.; Yao, J. Parallel simulation of fully-coupled thermal-hydro-mechanical processes in CO2 leakage through fluid-driven fracture zones. Int. J. Greenh. Gas Control. 2015, 34, 39–51. [Google Scholar] [CrossRef]
  24. Gregory, A.R.; Kendall, K.K.; Lawal, A.S. Study Effects of Geopressured-Geothermal Subsurface Environment on Elastic Properties of Texas Gulf Coast Sandstones and Shales Using Well Logs, Core Data, and Velocity Surveys; Bureau of Economic Geology: Austin, TX, USA, 1980. [Google Scholar]
  25. Molina, O.; Vilarrasa, V.; Zeidouni, M. Geologic Carbon Storage for Shale Gas Recovery. Energy Procedia 2017, 114, 5748–5760. [Google Scholar] [CrossRef] [Green Version]
  26. Bruno, M.S.; Nakagawa, F.M. Pore pressure influence on tensile fracture propagation in sedimentary rock. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1991, 28, 261–273. [Google Scholar] [CrossRef]
Figure 1. The separation distance between the low pH diffusion front and the high alkalinity sweep front as a function of seepage velocity, reservoir temperature, and time.
Figure 1. The separation distance between the low pH diffusion front and the high alkalinity sweep front as a function of seepage velocity, reservoir temperature, and time.
Minerals 12 00752 g001
Figure 2. The observation well bottom hole pressure and temperature measurements during the month of October 2004. The injection period is represented by the yellow highlight. The breakthrough time is represented by the green vertical line.
Figure 2. The observation well bottom hole pressure and temperature measurements during the month of October 2004. The injection period is represented by the yellow highlight. The breakthrough time is represented by the green vertical line.
Minerals 12 00752 g002
Figure 3. A mesh representing a 500 square-meter area of the Frio Formation. A Neumann boundary condition is imposed on the y-z plane at x = 0, such that the temperature gradient is 0 to implement a symmetry assumption on that wall. All other walls have a Dirichlet boundary condition of T = Ti imposed to implement a far-field constant reservoir temperature assumption. The simulated injection wellbore perforation interval was configured to −1527.5 m to −1533.75 m with respect to the z-axis at a y-axis location of 500 m. The injection fluid and reservoir initial temperature was configured to be Ti = 58.4 C. The injection mass flow rate was configured to be 3.0 kg/s.
Figure 3. A mesh representing a 500 square-meter area of the Frio Formation. A Neumann boundary condition is imposed on the y-z plane at x = 0, such that the temperature gradient is 0 to implement a symmetry assumption on that wall. All other walls have a Dirichlet boundary condition of T = Ti imposed to implement a far-field constant reservoir temperature assumption. The simulated injection wellbore perforation interval was configured to −1527.5 m to −1533.75 m with respect to the z-axis at a y-axis location of 500 m. The injection fluid and reservoir initial temperature was configured to be Ti = 58.4 C. The injection mass flow rate was configured to be 3.0 kg/s.
Minerals 12 00752 g003
Figure 4. The contours of pH and [ H C O 3 ] on a slice along the perforation zone after 4 days of injection at a mass rate of 3 kg/s. The background is color-mapped to the formation water temperature.
Figure 4. The contours of pH and [ H C O 3 ] on a slice along the perforation zone after 4 days of injection at a mass rate of 3 kg/s. The background is color-mapped to the formation water temperature.
Minerals 12 00752 g004
Figure 5. The simulated 30-day temperature profile at grid point (x = 30 m, y = 250 m, and z = −1530.625 m), corresponding to the location where the Frio Pilot observation well temperature sensor was placed. In the left plot, the y-axis is temperature in °C and the x-axis is time in days. The right image is a plot of temperature on the x-y plane at z = −1530.625 m. In the left plot, one can see the early arrival of a relatively lower temperature, low-pH diffusion front followed by a relatively higher temperature, high-alkalinity sweep front.
Figure 5. The simulated 30-day temperature profile at grid point (x = 30 m, y = 250 m, and z = −1530.625 m), corresponding to the location where the Frio Pilot observation well temperature sensor was placed. In the left plot, the y-axis is temperature in °C and the x-axis is time in days. The right image is a plot of temperature on the x-y plane at z = −1530.625 m. In the left plot, one can see the early arrival of a relatively lower temperature, low-pH diffusion front followed by a relatively higher temperature, high-alkalinity sweep front.
Minerals 12 00752 g005
Table 1. The shale, sandstone, and perforation regions of the injection and observation wells.
Table 1. The shale, sandstone, and perforation regions of the injection and observation wells.
ZoneInjection Well (ft, m)Observation Well (ft, m)
Top “A” sandstone
Top “B” Shale
Top “B” Sandstone−4966, −1514−4938, −1505
Top “C” shale−4978, −1517−4950, −1509
Top “C” sandstone−5034, −1534−5000, −1524
Top Perforation −5014, −1528
Bottom Perforation −5034, −1534
Bottom “C” sandstone
Table 2. The model parameters and diffusion coefficient values D computed at T = 65 °C using the model of Boudreau [7].
Table 2. The model parameters and diffusion coefficient values D computed at T = 65 °C using the model of Boudreau [7].
H+OHSiO2(aq)K+CO2(aq)HCO3CO3−2Na+Al(OH)3(aq)Ca+2Fe+2Mg+2
Tc54.4025.905.009.555.505.064.336.064.463.603.313.43
Tf1.5551.0940.5000.4090.3250.2750.1990.2970.2430.1790.1500.144
D1.6 × 10−49.7 × 10−53.8 × 10−53.6 × 10−52.7 × 10−52.3 × 10−51.7 × 10−52.5 × 10−52.0 × 10−51.5 × 10−51.3 × 10−51.3 × 10−5
Table 3. The regression parameters for Equation (3) from Johnson and Norton [13].
Table 3. The regression parameters for Equation (3) from Johnson and Norton [13].
a1a2a3a4a5
0.1470333593 × 1020.2128462733 × 103−0.1154445173 × 1030.1955210915 × 102−0.8330347980 × 102
a6a7a8a9a10
0.3213240048 × 102−0.6694098645 × 101−0.3786202045 × 1020.6887359646 × 102−0.2729401652 × 102
Table 4. The initial and boundary formation lithology and porosity.
Table 4. The initial and boundary formation lithology and porosity.
MineralVolume FractionGrain Radius, mm
Shale
Calcite, CaCO30.200.001
Quartz, SiO20.280.02
Halite, NaCl0.000.01
K-Feldspar, KAlSi3O80.010.03
Illite, (0.65 K, 0.08 Na) (2.27 Al, 0.14 Fe, 0.2 Mg)3.41 SiO10(OH)20.410.0001
Initial shale porosity φ = 10%
Sandstone
Calcite, CaCO30.000.001
Dolomite, CaMg(CO3)20.000.020
Quartz, SiO20.550.02
Halite, NaCl0.000.01
K-Feldspar, KAlSi3O80.130.03
Initial sandstone porosity φ = 32%
Table 5. The poroelastic property definitions.
Table 5. The poroelastic property definitions.
PropertyShaleSandstone
Bulk modulus, GPa [24]10.18.83
Shear modulus, GPa [24]11.09.64
Biot Modulus, GPa8.8212.56
Biot–Willis coefficient, αBW0.350.79
Poisson’s ratio [25]0.430.32
Tensile Strength, Pa [26]8.25 × 1068.45 × 106
Table 6. The initial and boundary formation water composition.
Table 6. The initial and boundary formation water composition.
SoluteConcentration, M
Initial Formation Water
OH1.0 × 10−7
Ca2+2.5 × 10−3
Al(OH)3,aq1.7 × 10−5
K+5.0 × 10−5
SiO2,aq1.0 × 10−3
CO2,aq2.0 × 10−3
Na+4.0 × 10−1
Fe2+1.0 × 10−4
Mg2+1.0 × 10−5
Cl4.0 × 10−1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Paolini, C. Effects of Charged Solute-Solvent Interaction on Reservoir Temperature during Subsurface CO2 Injection. Minerals 2022, 12, 752. https://doi.org/10.3390/min12060752

AMA Style

Paolini C. Effects of Charged Solute-Solvent Interaction on Reservoir Temperature during Subsurface CO2 Injection. Minerals. 2022; 12(6):752. https://doi.org/10.3390/min12060752

Chicago/Turabian Style

Paolini, Christopher. 2022. "Effects of Charged Solute-Solvent Interaction on Reservoir Temperature during Subsurface CO2 Injection" Minerals 12, no. 6: 752. https://doi.org/10.3390/min12060752

APA Style

Paolini, C. (2022). Effects of Charged Solute-Solvent Interaction on Reservoir Temperature during Subsurface CO2 Injection. Minerals, 12(6), 752. https://doi.org/10.3390/min12060752

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop