Next Article in Journal
The Geological Structure and Tectonic Complexity of Northern Thessaly That Hosted the March 2021 Seismic Crisis
Next Article in Special Issue
A Neural Network Model for Estimation of Failure Stresses and Strains in Cohesive Soils
Previous Article in Journal
Stress–Strain–Time Description and Analysis of Frozen–Thawed Silty Clay under Low Stress Level
Previous Article in Special Issue
Development of Soil Moisture Content and Soil Matric Suction Model Based on Field Instrumentation and Electrical Resistivity Imaging (ERI) for Highway Slopes Constructed on High Expansive Clay Soil
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Finite Element Simulations of Fluids Leakage through the Faulted Reservoir

by
Mohammad Nurul Islam
Transportation & Infrastructure, WSP-USA, Morristown, NJ 07960, USA
Geotechnics 2022, 2(4), 908-934; https://doi.org/10.3390/geotechnics2040043
Submission received: 5 September 2022 / Revised: 17 October 2022 / Accepted: 18 October 2022 / Published: 29 October 2022
(This article belongs to the Special Issue Recent Advances in Geotechnical Engineering)

Abstract

:
Carbon dioxide (CO2) capture and storage (CCS) in geological formation as a supercritical fluid is a viable option to reduce anthropogenic greenhouse gas emissions. Due to the density difference between CO2 and formation fluid, CO2 shows a buoyant tendency. Thereby, if CO2 migrates towards the fault in a compromised faulted reservoir, it may escape the storage reservoir. Therefore, it is essential to predict fluids leakage through the faulted reservoir into the aquifer, associated pressure development, and fluids properties over time to assess associated risk and quantification of leakage. We present finite element simulations of miscible fluids flow through the faulted reservoir to elucidate this behavior. There are very few attempts to model multicomponent fluids non-isothermal model during phase change including the Equation of State (EoS) which we addressed by coupling the mass balance equation of fluids, the fractional mass transport, and the energy balance equation. To obtain fluids mixture thermo-physical properties, we used the Peng-Robinson EoS. For validation of the coupled formulation, we compared the simulation results with Ketzin Pilot project field monitoring data, which shows good agreement. A faulted reservoir comprised of five layers is used to investigate fluids leakage through a compromised reservoir. These layers are a CO2 storage reservoir, overlain by alternating caprocks and aquifers. We also considered three different CO2 injection rates to study the injection rate effect to assess the pressure buildup during injection process. We present the thermal effect by comparing the isothermal and the non-isothermal conditions. For the latter case, we assumed three different thermal gradients. Additionally, to assess the fault aperture effect, we studied three different apertures. We observed that developed pressure and fluids properties have effects on injection rates, temperature gradient, and fault aperture. Additionally, such responses in the near-field and the far-field from the injection well are critical to assess the risk, which we discussed in this paper.

1. Introduction

After the industrial revolution in 1750, the amount of greenhouse gases (GHG) increases that trap and holds energy from solar radiation in the earth’s atmosphere in the form of heat, have contributed to a gradual rise in mean atmospheric temperature. Among these greenhouse gasses, carbon dioxide (CO2) emitted from large-scale combustion of fossil fuels is seen as a primary contributor to this effect (Bachu [1]; Volk [2]; EPA [3]; Metz et al. [4]). After 1990, GHG in the United States increased by 3.5% while the contribution of CO2 was 82.2%, and CO2 concentration is increasing [5]. In this regard, Eggleton [6] reported that the global amount of CO2 increased by 45% after the industrial revolution. Bruckner et al. [7] presented that the estimated CO2 emission from the energy section in 2010 was 14.4 GtCO2/year and is projected to increase by 25 to 30 GtCO2/year by 2050. To offset anticipated future increases in emissions, the capture of CO2 emissions from large anthropogenic point sources, utilization, and long-term, safe geologic storage (CCUS) is being actively considered as a means to manage and reduce atmospheric CO2 concentrations.
The essential steps in CCUS are separation of CO2, compression to a dense-phase fluid, transportation from the location of separation to the point of disposal, and injection into deep geologic formations for long term storage. Several candidate geologic storage targets are described in the scientific and industry literature (see also Nordbotten and Celia [8]). These include depleted oil and gas reservoirs, coal seams, and deep saline formations as the most prospective storage options, with basalt formations and unconventional organic-rich shale formations have also been considered (see also NETL [9]).
Viable storage targets will have a storage interval with sufficient injectivity and capacity to store large volumes of CO2, an injection depth sufficient to ensure that CO2 is stored in dense fluid phase (liquid or supercritical fluid; nominally 800 m) (Rutqvist and Tsang [10]; Benson and Cole [11]), and a laterally extensive caprock sequence above the storage interval capable of preventing unwanted vertical fluid migration. Bouyance of the stored CO2 relative to the connate fluid (e.g., brine) and elevated pressure of the storage reservoir in response to large-scale CO2 injection are driving forces for potential unwanted migration of CO2 and brine into overlying receptors of concern (groundwater formations and the atmosphere). Geomechanical, thermal, and geochemical phenomena also influence storage system performance.
To capture the full complexity of the dynamics of geologic carbon storage systems requires coupled simulation of fluid flow, energy, and fractional mass transport. This is the motivation for model formulation in this study. Presence of a fault through an impermeable caprock presents a potential pathway for fluid migration that could potentially impact the containment of CO2 [12]. Consideration of the potential of faults to leak fluids flow behavior through the faulted reservoir during and after the CO2 injection in the storage reservoir considering the thermodynamics of fluids, which is the point of interest to study the fluids migration behavior.
The finite element (FE) modeling of multicomponent fluids (e.g., gas and liquid) mixture properties (e.g., density, viscosity) change is challenging during change of pressure and temperature. For example, fluids temperature and pressure in a CO2 storage tank is different than the storage reservoir. Therefore, the bottom hole pressure suddenly increases due to the fluid volume expansion for temperature variation. Additionally, the accuracy of results and numerical convergence conditions in FE solutions are important (Helmig [13]). A similar intricacy also involved in enhanced oil and gas recovery, bio-mechanics and hazardous waste disposal projects, where fluids flow through a porous media and non-isothermal condition exist in the phase transition process. For such a case, it is important to introduce the equations of state (EoS) in the governing process of the model formulation. The isothermal study of fluids and the non-isothermal study of a single fluid are well documented in the literature (Bear and Bachmat [14]). However, there are very few attempts to model multicomponent fluids non-isothermal model during phase change including the EoS. The reduced order models (ROM) are most popular to address a case-specific problem (Vasylkivska et al. [15]). Among others, Rutqvist and his co-workers [10,16] also coupled FLAC and TOUGH for CO2 sequestration applications. Rutqvist [16] reported the method as a ‘pragmatic approach’. In this regard, Beck et al. [17] presented that node sharing in finite difference approach to achieve a solution of two different codes may results accuracy issues during incremental time steps. There are few literatures where coupled FE solutions are available showing an industry scale verification representing the accuracy of the results. To address such a challenge, we presented a couple FE model formulation, validation and application in this paper.
In this paper, we present a finite element simulation that couples the fluid mass balance equation, the energy balance equation, and the fractional mass transport. We also use the Helmholtz free energy (Wagner and Pruß [18]) for two fluids: carbon dioxide and water, and the Peng-Robinson Equation of State (EoS) (Peng and Robinson [19]) to describe the thermo-physical properties of fluid mixtures. A comparison of the field monitoring measured pressure response of the Ketzin CCUS project (see Kempka and Kühn [20]) with simulated pressure using a version of the model developed herein tuned to Ketzin project geologic and injection data from literature and TOUGH (Pruess [21]) provides validation of our modeling approach.
We present a coupled finite element model to address fluids flow through the faulted reservoir. The study domain is comprised of a CO2 storage reservoir, aquifers, and caprocks. To ensure the supercritical CO2, we also assumed that the model domain’s top surface is located 1000 m below the ground surface and considered miscibility of aqueous and CO2-rich phases. A constant fault aperture is assumed and potential fault propagation is not considered assuming the parallel plate theory [22]. The model captures pressure evolution with CO2 injection and distribution of fluids mole fraction (e.g., CO2, H2O) (see Table A1, Figure A1 and Figure A2). To reduce the computational time, we considered a simplified geometry and 300 days of simulated site performance.
The following sections present detail of the numerical model, finite element implementation, benchmark validation, and application of the model to capture fluid flow behavior through the faulted reservoir.

2. Governing Equations

The composite medium in the model represents fluid-filled porous rock, so is comprised of fluids phase (f) and solid phase (s). The fluid phase is consisting water (H2O) and carbon dioxide (CO2), with pore spaces in the solid-phase represented as macropore (inter-aggregate, e.g., fracture) and micropore (intra-aggregate, e.g., solid matrix) following the dual porosity conceptualization presented by Wilson and Aifantis [23]. The mobility of fluids depends on their physical disposition relative to those pores (i.e., fluids in macropores are mobile or relatively mobile while fluids in micropores are immobile). We obtain such a dual porosity concept following the macroscopic continuum method and using the Representative Elementary Volume (REV) and satisfying the inequality (see also Bear and Bachmat [14]). When fluids are composed of more than one species and their flow accompanies its transport, the governing equations require both the mass balance equation and the fractional mass transport equation (Hassanizadeh [24]; Bear and Bachmat [14]). In the following section, we discuss both of them and the energy balance.
We assumed that (i) the porous media is non-deformable considering solid domain and fluid pressure, (ii) fluid constituents are miscible and satisfy the local thermodynamics equilibrium condition, (iii) fluid phases obey the Darcy’s law which complies with the Reynolds number and (iv) fluid composition is governed with respect to its mass fraction (Bear and Bachmat [14]).
Considering the above assumptions and Hassanizadeh [24]; Bear and Bachmat [14], the balance equations are present as follows:

2.1. Kinematics

Truesdell and Noll [25] presented suppositions which, also known as the “metaphysical principles” to describe the continuum mechanics of a single constituent to many constituents, i.e., mixtures. Atkin and Craine [26]; Bowen [27] also illustrated reviews of the continuum theory of mixtures. The mixture can be homogeneous (e.g., miscible fluids, like soluble salt in water) or heterogeneous (e.g., immiscible fluids or fluid-solid, like, water-sand) which are also state-dependent, and one may obtain exceptions by changing the physical conditions, e.g., temperature. CO2 injection in a brine or water saturated storage reservoir presents a case in which (i) fluids may be miscible or immiscible depending on the state variables and (ii) the body may be deformable or non-deformable considering fluid flow in a porous media and associated initial and boundary conditions. In this paper, we limited our discussion to the miscible case and relative rigid body, while in Islam et al. [28], we presented a coupled formulation for the immiscible case and deformable body. Moreover, in Islam and Gnanendran [29] and Islam et al. [30], we also discussed elasto-viscoplastic solutions of the material body considering the creep and relaxation based on the two-phases theory, i.e., solid and fluid couple.
In this section, we defined “point” and “particle” as a specific location in a particular space and material point, i.e., a volumetric element), respectively. Assuming a domain, Ω d ; d = 1 ,   2 ,   3 composed of N continuous bodies, while at any reference time ( t = 0 ) , each continuous body comprises particles and their constituents. We present the configuration of a continuum body, where the position vector at t = 0 and t = t are X = X k I ^ k and x = x i e ^ i . We present the motion of a deformed body in terms of the Lagrangian formulation and Eulerian formulation as
x α = χ α ( X a , t )
X α = χ α 1 ( x α , t )
where χ α is the deformation function and α = 1   to   N , which represents the number of constituents. Equation (1) represents that at any time t = t , the place occupied by X α . We also assumed that χ α and χ α 1 continuously differentiable and sufficient number of derivatives exist (see Gurtin [31]).
The velocity and acceleration of α t h constituent at time t = t are obtained as
v α = x ˙ α = χ α t ( X α , t )
a α = x ¨ α = 2 χ α t 2 ( X α , t )
Assuming, ρ ( x , t ) and ρ α ( x , t ) are the total mass density of the mixture and the density of the a t h constituent or bulk density, respectively, so that
ρ ( x , t ) = a = 1 N ρ α ( x , t )
In Equation (5), ρ α can be written as
ρ α ( x , t ) = d m α d V
where d m α and d V are the instantaneous mass of the α t h constituent and volume occupied by constituents. Moreover, the mole fraction for each component in the fluid ( w α f ) can be defined as
w α f = ρ α ρ
α = 1 N w α f = 1
It is to note that relation between the true density ( γ α ( x , t ) ) and the bulk density ( ρ α ( x , t ) ) are obtained as (Bowen [27])
n α ( x , t ) = ρ α ( x , t ) γ α ( x , t )
where n ( x , t ) represent the volume fraction or porosity of the α t h constituent. Additionally, we also assumed that
α = 1 N n α ( x , t ) = 1
Again, the mean velocity of any mixture ( v ) is presented as
v = 1 ρ α = 1 N ρ α v α
The relative velocity or the diffusion velocity ( w α )   of the α t h constituent given by
w α = v α v
As the diffusion velocities are not independent, therefore (Bear and Bachmat [14])
α = 1 N ρ α w α = 0
Partial differentiation of Equations (1) and (2) with respect to X α and x α provide the material deformation gradient and the spatial deformation gradient as
F α = x α X α
H α = X α x α
It is worth to note that the relation between Equations (14) and (15) can be obtained through the ‘chain rule’ as follows
x i X j X j x k = X i x j x j X k = δ i k
where   δ represents the Kronecker delta. It is to note that the Jacobian determinant J = | x i X j | should be non-singular. Additionally, the material displacement gradient ( u i X j ) and the spatial displacement gradient ( u i x j ) can be obtained by partial differentiation of the displacement vector ( u i ) corresponding coordinate systems as follows (see also Bowen [27])
u i X j = x i X j δ i j
u i x j = δ i j X i X j
We assumed that at any time t = t ,   ψ be any quantity (e.g., scalar, vector, or tensor), and its instantaneous velocity can be defined as the material derivative as
D ψ i j ( x , t ) D t = ψ i j ( x , t ) t + ψ i j ( x , t ) x k x k t
In Equation (19), the right-hand side’s first part describes the local rate of change, while the second one represents the convective rate of change. Among others, Coleman and Noll [32] presented a functional form of ψ , associated theories and constitutive laws.
In the following section, we present balance equations and constitutive laws in this paper are presented in Appendix A.

2.2. Mass Balance Equation for Fluids

The mass balance equations for fluids in a multiphase porous media can be written as
t ( n ρ ) + · [ J F ] = Q ρ
where   is the divergence operator. n   is the porosity of the porous medium. ρ is the mixture density.   J F = J F A + J F D is the fluid flux while J F A   and J F D   are the advective flux and the diffusive flux, respectively (see Appendix A). Q ρ is the source or sink term.

2.3. Mass Balance Equation for the Fractional Mass Transport

To account thermodynamically admissible fluids flow and heat transfer in multiphase porous media, Hassanizadeh [24] presented the fractional mass transport as follows
n ρ w α f t + n ρ v · w α f · ( J α ) = Q w α f
where   J α is the dispersive mass flux. w α f the mole fraction for each component in the fluid. v is the mean velocity of any mixture. Q w α f is the source term in the fluid for the fractional mass transport.

2.4. Energy Balance Equation

For the solid phase and the fluid mixture, the heat transport equations can be obtained as, respectively,
( 1 n ) ρ s c p s T t + · [ J S T ] = Q T s
n ρ c p f T t + ρ c p f n v · T + · [ J F T ] = Q T f
where   J S T and J F T are the heat flux for the solid and the fluid, respectively, Q T s and Q T f are the heat source for the solid and the fluid. c p s and c p f are the specific heat capacity for the solid and fluid, respectively (see Appendix A).
We present detailed fluids properties and constitutive relations associated with the balance equations in Appendix A and Appendix B, respectively.

3. Coupled Finite Element Solutions

In matrix notation, the coupled equation is written as
M e X ˙ + ( A e + K e ) X = F e
where M e ,   A e and K e are the mass, advection, and Laplace matrices of elements, respectively, while F demonstrates the right-hand side vector. Elements of matrices in Equation (24) depend on solutions vectors X = { ψ , w α f   } T r , where T r represents the transpose and ψ = { p , T } T r . It is worth to note that for any time step, we solve p and T monolithically while following the staggered approach, we solve w k f for each component. Additionally, in Appendix C, we present the expanded form of Equation (24) and its elements. Using the summation convention, we also obtained global matrices.
Using the finite difference method in Equation (24), we find the time-stepping equation as follows (see also Wood [33])
[ C + θ Δ t K ] n t + θ X n t + 1 [ C ( 1 θ ) Δ t K ] n t + θ X n t Δ t F n t + θ = 0
where n t and ( n t + 1 ) represent time steps at time t n and t n + 1 , respectively. θ is the time integration parameter. Islam and Gnanendran [34] presented a detailed derivation of Equation (25) and the effect of θ in numerical results. In this paper, we assume θ 0.5 to obtain numerical stability. We used the SpBICSGSTAB solver (see Chen et al. [35]) and the PICARD solver (see Paniconi and Putti [36]), respectively, to solve linear and non-linear equations.
In the next section, to validate the performance of model formulations as discussed, we compare Ketzin CO2 storage measured data as reported by Kempka and Kühn [20] with the finite element model results presented herein and with results predicted using the TOUGH simulator (Pruess [21]).

4. Benchmark Validation

The Ketzin pilot project was the first European onshore pilot project to store CO2 in a deep saline aquifer (Kempka and Kühn [20]; Kempka et al. [37]). Among others, Martens et al. [38] and Chen et al. [39] presented descriptions of the geology at the Ketzin pilot project site. The project site is located in the Ketzin-Roskow anticline 25 km west of Berlin, and targets injection into the Stuttgart formation at a depth close to 630 m to 710 m. Above and below the Stuttgart formation, there are the Weser Formation and the Grabfeld formation, respectively, which are considered as no-flow boundaries (Kempka and Kühn [20]). Injection operations at this site ran from June 2008 through August 2013. Over that time, a total of 67,271 metric tons of CO2 were injected (see also Kempka et al. [37]). Three wells were drilled at a depth between 630 to 710 m, an injection well (Ktzi 201), and two observation wells (Ktzi 202 and Ktzi 200).
Observation wells Ktzi 202, and Ktzi 200 are offset from the injection well by 50 m and 112 m, respectively. Additionally, from Liebscher et al. [40], we obtained the injection rate (see also Figure 1) and injection well, and the observation well measured bottom hole pressures and temperatures.
From Kempka and Kühn [20], we obtain the initial pressure and the temperature of the domain as 6.2 MPa and 34 °C (=307.15 K); the storage reservoir was initially fully water saturated (i.e., initial CO2 gas phase saturation was zero). Along the injection boundary line, the CO2 gas phase saturation was assumed to be 1, while the temperature was 34 °C. Salinity in the aqueous phase is neglected in the model; total porosity of the solid medium (rock) is 0.25; the isotropic permeability is 4.6 × 10−14 m2; the thermal expansion is 1 × 10−5; the specific heat capacity is 1.2 × 103 J/kg/K and the thermal conductivity is 2.5 W/m/K. Pressure response to CO2 injection was predicted for a simulated period of approximately 27 months. For finite element geometry, we consider the well radius is 0.1 m (Chen et al. [39]). For simplicity, we also assumed the axisymmetric condition. We extended the width of the storage reservoir to 1500 m to minimize the boundary effect in fluid pressure calculation. In Figure 2, we compare 27 months of measured data with predicted data, which has a good agreement, testifying to the accuracy of the model prediction.
In the next section, we used the model in this paper to predict CO2 leakage through the faulted reservoir. The considered geometry is a simplified, hypothetical section. However, to explain the possible fluids flow through the faulted reservoir, we considered a condition approximately representative of the Ketzin field project (based on previously described literature). This example is useful to demonstrate the simulation of coupled fluids flow in a fault to consider the potential fluid migration and pressure response that may result from such a case.
For the Ketzin site, the carbon dioxide (CO2) is transported to the sequestration project by road transportation. Then, before starting the injection, CO2 is stored in the transitional storage tank, where approximate temperature and pressure are −18 °C and 2.1 MPa, respectively (see Ketzin project and Martens et al. [38]; Chen et al. [39]). This type of CO2 is also known as the cryogenic CO2 (see also Oldenburg [41]).
Initial reservoir conditions and model assumptions in a hypothetical model follow the same assumptions as presented in the Benchmark Validation section above. The rapid change of fluid conditions from −18 °C and 2.1 MPa in tank storage to 34 °C and 6.2 MPa in the geologic storage reservoir would result in a nearly three-fold volumetric increase (Liebscher et al. [40]). Additionally, the initial state of CO2 in the Ketzin field project is close to the critical point of CO2 (31 °C and 7.38 MPa). Hence, the possibility of CO2 phase transition in the injection well and the reservoir is also high. Thereby, for an identical field case with a fault, it is essential to couple fluids flow, thermal flow, the equation of state and thermo-physics in the model formulation to capture fluids leakage through the faulted reservoir. Additionally, we discuss these physics for an arbitrary case in the next section.

5. Fault Leakage Results

We assume that the study domain is located 1000 m below the ground surface to ensure CO2 is in the supercritical phase (Rutqvist and Tsang [10]). We consider a simplified geometry to reduce the computational time (Figure 3). The model consisted of a stacked system with two aquifers, two caprocks, and an underlying CO2 storage reservoir. We also assume a pre-existing fault in the geometry with minimal fault displacement (Figure 3). In this paper, we assume no fault propagation or change to fault properties during simulations. To minimize the computational time, we also consider a constant CO2 injection rate and monitor the porous media for 300 days. Additionally, we investigate three different injection rates effect, the thermal effect and fault aperture effect in three observation points.
We consider the hydrostatic condition with an average pressure gradient of 10,251.145 Pa/m (Ebigbo et al. [42]). Thereby, we also obtain the initial state pressure at the top and the bottom of the study domain. We also introduce such a condition at the left and the right lateral boundary (the Dirichlet boundary condition), which are constant during the simulation. We also accept that the storage reservoir is initially water saturated, and the CO2 saturation is zero. However, we assume the CO2 saturation in the injection point is 1. We assume that the Neumann boundary condition at the top and the bottom of the domain by restricting the flow along the boundary. We consider both the isothermal and the non-isothermal conditions. At the CO2 injection point (Figure 3), the temperature is 307.15 K. We introduce a low-temperature difference between the study domain boundaries to restrict thermal shock. We assume that caprock layers are impermeable (Ebigbo et al. [42]). We consider that the density, the thermal conductivity and the specific heat capacity of the porous medium are 2650 kg/m3, 3.5 W/(m K) and 750 J/(kg K), respectively (see also Class et al. [43]). We also assume that the porosity, permeability, and diffusion coefficient of the porous media are 0.15, 2 × 10−14 m2 and 6 × 10−7 m2/s respectivly (Ebigbo et al. [42]; Class et al. [43]). We consider that the porous domain is homogeneous and isotropic. Additionally, we obtain fluids properties from Appendix A and Appendix B; Figure A1 and Figure A2.
For the fault, we consider 1D line elements with constant fracture aperture ( b ) and corresponding intrinsic fault permeability is k f = b 2 12 (Bear and Berkowitz [44]). We also study the sensitivity of fault aperture considering three types of aperture. For all cases, we assumed fault porosity of 100%.
For finite element mesh generation, we used 103 line elements for the fault and 12,560 triangular elements for other sub-domains (see Figure 3). As the injection rate is slow enough to ensure the applicability of the small deformation theory, we considered the first-order finite elements (Zienkiewicz et al. [45]) with a small time increment ( Δ t ) to ensure the quasistatic condition (Bear and Bachmat [14]).
In the next section, we discuss the injection rate effect, thermal effects, and fault aperture effect.

5.1. Effect of Injection Rate

A safe CO2 injection rate is essential for any storage reservoir to avoid and or minimize fluids leakage. In a faulted storage reservoir, an excessively high CO2 injection rate may lead to in situ stress changes that can cause fault activation or dilation with associated effective permeability change, seal fracturing/fracture propagation, and potential increase in leakage rate. While the dynamics of coupled fault mechanics and leakage is not explicitly considered in this manuscript, the relationship between injection rate and leakage through a pre-existing fault of constant aperture is simulated to explore the range of potential leakage. Uncertainty in the sensitivity of the model to injection rate is explored using three different CO2 injection rates: 0.13773 m3/s, 0.27546 m3/s and 0.41319 m3/s. Simulation results are considered at several observation points (Figure 3) where we predict and compare distributions of pressure, CO2 mole fraction, density, viscosity, specific heat capacity and thermal conductivity. For the sensitivity analyses in this section, we consider fracture aperture is 100 µm; the temperature gradient is 0.02 °C per meter; the longitudinal dispersivity and transverse dispersivity are 5.0 m and 0.5 m, respectively.
In Figure 4a–c, we observe that close to the injection point pressure evolvement is high, and pressure decreases to the far-field observation points. We also notice in observation point 1, after 18 days, pressure starts to fall for all injection rates, while in observation points 2 and 3, with the injection rate increase, pressure decrease rate reduced. It is noteworthy that during the CO2 injection, the mean effective stress decreases [10]. In the early stage of CO2 injection, if the buildup pressure exceeds the allowable stress limit of the porous medium, a fracture may propagate along the compromised faulted reservoir [46]. Moreover, the ratio of the minimum to the maximum injection rate, we assume 3, while we observe that the maximum pressure build-up ratios are 1.128, 1.05 and 1.04, for observation points 1, 2 and 3, respectively. Furthermore, for all injection rates, along the observation points 1 and 2, incremental pressure approach to the steady-state condition, while in observation 3, after 300 days such a condition is not yet achieved, which also justifies the CO2 distribution as we illustrate in Figure 5.
Figure 5 shows the effect of varying CO2 injection rates on the distribution of CO2 saturation (fluid phase mole fraction). We observe that with injection rates 0.13773 m3/s, 0.27546 m3/s and 0.41319 m3/s, in the observation point 1, CO2 reaches full saturation after 231 days, 99 days and 60 days, respectively. CO2 saturation in the top formation observation point (Figure 3 and Figure 5), ranges from 0.013 and 0.944 for injection rates 0.13773 m3/s and 0.41319 m3/s, after 300 days.
To avoid the repetition of similar illustrations, in Figure 6, we present the spatial distribution of pressure, CO2 saturation, and fluids properties for the study model domain (Figure 3) after 300 days considering the largest CO2 injection rate (0.41319 m3/s).

5.2. Thermal Effect

The temperature of CO2 between the storage reservoir and the transitional storage tank is different. Therefore, it is essential to study the thermal effect of safe reservoir design. Oldenburg [41] also reported that during CO2 injection, due to the thermal effect, residual fluids in the storage reservoir may freeze and developed thermal stress could generate fracture in the storage reservoir. To explore the potential impact of this phenomenon, we compare injection under both the isothermal and the non-isothermal conditions (considering three thermal gradients: 0.02, 0.025, and 0.030 °C per meter).
Mathias et al. [47] also demonstrated that during injection, CO2 displaces reservoir fluids causing a gradual pressure increase; as CO2 expands from the injection point, the reservoir pressure begins to drop. We also observe a similar phenomenon in simulations herein. Figure 7 illustrates the thermal effect on the pressure distribution close to the injection point and the far-field observation points. We find that with the increase of temperature in the near field, there is a decrease in pressure. However, over time, a decrease in temperature and the corresponding increase in pressure is observed in the far-field. We also notice that the temperature gradient effect is significant in the CO2 mole fraction distribution, and with the increase of the thermal gradient, the CO2 mole fraction also increases (see Figure 8).
We also present the effect temperature distribution in Figure 9 and Figure 10 for the highest thermal gradient Δ T = 0.03 °C per meter, and the isothermal condition. We find that the temperature starts to drop for the non-isothermal state near the injection point. Comparing Figure 9 and Figure 10, we find that the distribution of fluids properties (e.g., density, viscosity) is influenced by the thermal gradient. Additionally, the maximum developed pressure and density in the isothermal condition are higher than in the non-isothermal condition.

5.3. Effect of Fracture Aperture Size

During and after CO2 injection, the fracture aperture is playing an important role in the fluid flow pathways. Therefore, in this section, we present the fracture aperture size effect in the near-field and the far-field through three observation points. We consider three different fracture apertures (100 µm, 200 µm, and 300 µm) and injection rates of 0.13773 m3/s, 0.27546 m3/s and 0.41319 m3/s. In the earlier section, we presented the thermal effect, considering the isothermal and the non-isothermal conditions. Therefore, herein we only consider the thermal gradient of 0.02 °C per meter.
In Figure 11, Figure 12 and Figure 13 for all observation points we present the pressure distribution effect on the injection rate and the fracture aperture size. We observe that for observation point 1 (near field to the CO2 injection) for relatively smaller aperture, the pressure is relatively high compared to the increase of fault aperture. In this paper, as we consider the fracture aperture is constant throughout the simulation; therefore, with the increase of aperture, the path to fluid flow experienced less resistance than a small conduit aperture. As a result, in the near field, the developed pressure for any specific thermal gradient is inversely proportional to the fault aperture. In contrast, for an identical thermal gradient and the near field condition, the developed pressure is proportional to the CO2 injection rate. However, in the relatively far field (e.g., Observation points 2 and 3), due to the energy dissipation during flow path in the space, the observed response for the pressure is opposite to the near field condition. Therefore, close to the injection point, the presence of fracture requires additional attention to avoid the instabilities due to the developed pressure in the near field.
In Figure 14, Figure 15, Figure 16 and Figure 17, we present CO2 distribution at three observation points. The CO2 leakage through the larger aperture is high compared to the small aperture. Therefore, we observe that in all observation points, the CO2 presence is inversely proportional to the aperture size. In contrast, CO2 distribution increase with the increase of the injection rate and the thermal gradient.

6. Conclusions

In this paper, we present a coupled finite element solution combining the mass balance equations of fluids phase, the energy balance equations, and the fractional mass transport equation. For the solution of pressure and temperature variables, we consider the monolithic coupling, while for the mass transport, we used the staggered approach [48]. Additionally, to obtain the thermo-physical fluids properties, we also use the Peng-Robinson Equation of State (EoS). For the validation of solution algorithms, we compare coupled finite element (FE) results with the Ketzin pilot project measured responses for 27 months and with the TOUGH predicted results, which show a good agreement and testify the performance of FE results.
We also present fluids flow through a faulted reservoir, demonstrating the effect of the CO2 injection rate, the thermal effect, and the fault aperture effect. We observe that during injection, the buildup of pressure and its distribution, the thermal gradients also change the thermo-physical fluids properties. Moreover, close to the injection point’s low temperature compared to the domain also reduce the temperature along the fluids flow path, which also justifies the importance of the thermal effect in the CO2 sequestration project. Herein, to avoid thermal shocks, we consider a low temperature gradient. However, if the temperature drop is too high, it may adversely impact the CO2 injection.
We observe that during CO2 injection, for a specific thermal gradient, the developed pressure and CO2 concentration is proportional to the injection rate. Again, for the non-isothermal condition, the developed pressure in the near field is high compared to the far field condition. In contrast, in the far field, the pressure is proportional to the thermal gradient, and the magnitude is lowest for the isothermal condition. We also observed that CO2 concentration for all observation points is proportional to the thermal gradient. For a specific injection rate, the density, the viscosity, and the thermal conductivity decrease with the thermal gradient increase. However, the specific heat capacity increase with the increase of the thermal gradient. Additionally, in the near field, for a particular thermal gradient and an injection rate, the developed pressure is inversely proportional to the fracture aperture. In contrast, for the far field, we noticed the opposite behavior. For all observation points, CO2 concentration is inversely proportional to the fracture aperture.
For the faulted reservoir, the safe injection rate is crucial to avoid the instabilities of the storage reservoir. In such a condition, developed pressure and CO2 saturation are both crucial. The temperature of CO2 in the transitional storage tank is relatively low compared to the storage reservoir. Therefore, it is essential to assess the effect of injection rate, temperature, and the aperture of fault in the near field and the far field condition to limit the CO2 leakage, which we addressed through a hypothetical faulted reservoir.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

This paper is dedicated to M. Massoudi.

Conflicts of Interest

The author declares no conflict of interest.

Appendix A. Constitutive Relations

Appendix A.1. For Fluids

In Equation (21), J F comprised of the advective flux ( J F A ) and the diffusive flux ( J F D ) . From the Darcy’s law, we obtain the expression of J F A while we find J F D from the Fick’s law (see also Bear and Bachmat [14]). Details are as follows
J F A = ρ n v = ρ k μ α ( p ρ g )
J F D = w α f ρ n ( v α v )
where v is the mean velocity of any mixture. k is the intrinsic permeability. μ α is the fluid phase dynamic viscosity. g is the gravitational acceleration. w α = v α v is known as the relative velocity or the diffusion velocity.
Applying the chain rule for fluids mixture pressure ( p ) , temperature ( T ) , and components mole fraction ( w α f ) , we find ρ t in Equation (1) as
ρ t = ρ p p t + ρ T T t + ρ w α f w α f t = ρ β p α f p t + ρ β T α f T t + ρ β w α f w α f t
where β p α f , β T α f and β w α f are the coefficient of isothermal compressibility, the isobaric thermal expansion and the solutal expansion, respectively.

Appendix A.2. For Fractional Mass Transport

In Equation (21), J α is defined as follows (Hassanizadeh [24])
J α = ρ D w α f
where D   i j     is the dispersion tensor and written as (Bear and Bachmat [14])
D   i j   = τ D k δ i j + α t | v | δ i j + ( α l α t ) v i v j | v |
where τ = n 3 is the tortuosity (see also Ghanbarian et al. [49]; Millington and Quirk [50]) and n is the porosity. Additionally, α l and α t are the longitudinal dispersivity and transverse dispersivity, respectively. D k is the diffusion coefficient of a mixture, which is a function of the binary diffusion coefficient, pressure and temperature.

Appendix A.3. For Energy

We assumed that the local thermal equilibrium between the solid phase and the fluids phase. Using the Fourier Law, combining the heat flux we obtained
J T = J S T + J F T = λ e f f T
where λ e f f is the coefficient of effective thermal conductivity and obtained as
λ e f f = λ f n + λ s ( 1 n )
where λ f and λ s   are the thermal conductivity of the fluid and the solid, respectively.
We also find, the effective heat capacity, ( ρ c p ) e f f as follows
( ρ c p ) e f f = ( 1 n ) ρ s c p s + n ρ c p f
The simplified form of the coupled energy balance equation for the porous medium can be written as
( ρ c p ) e f f T t + c p f ρ n v . T [ λ e f f T ] = ( n β T T p t + n β T T v . p ) n v . p + Q T
where Q T = Q T s + Q T f represents the heat conduction flux. β T defines the thermal expansivity. It is noteworthy that for the non-isothermal condition, in Equation (A9), we may ignore the right-hand side’s first part (also known as the Joule-Thomson Cooling) and the second part (also considered as the Viscous Heat Dissipation).
We present details of fluid properties for carbon dioxide (CO2) and water (H2O) in Appendix B.

Appendix B. Equations for Fluid Properties

Appendix B.1. Equation of State (EoS)

In the porous medium, fluids properties vary associated with pressure and temperature, which can be described by Equation of State (EoS). Following the ideal gas law, we obtained the mixture density as follows
ρ α = M α p z α ( p , T ) R T
where α represents CO2 and water. ρ α   is the density of a componential fluid. z α ( p , T ) is the compressibility factor of componential fluid. For the isothermal condition, z α is only the function of pressure. In the storage reservoir simulation, it is frequent to use z α ( p , T ) to demonstrate the gaseous mixture behavior. z α ( p , T ) is considered as the estimate of deviation from the ideal gas. For the ideal gas, z α ( p , T ) = 1 M α is the molar mass of componential fluid. R is the universal gas constant which is related to the specific heat at the constant volume and pressure. Additionally, T is the temperature. In Equation (A10), we obtain p from the Peng-Robinson Equation of State (see also Peng and Robinson [19]).
We present details of fluids constants for carbon dioxide (CO2) and water (H2O) (see also Peng and Robinson [19]).
Table A1. Constants for CO2 and H2O.
Table A1. Constants for CO2 and H2O.
Fluids ρ c T c p c ω M
Kg/m3KPa-kg/kmol
CO2467.6304.137,377,3000.2249144.01
H2O322647.09622,064,0000.34418.015
We assumed that the mixture’s molar volume is the summation of its componential partial volume (see also Sengers et al. [51]). Therefore, the compressibility factor (see Equation (A11)) of the mixture can be obtained as follows
z m i x = α w α f z α

Appendix B.2. Viscosity

Appendix B.2.1. Viscosity of Carbon Dioxide (CO2)

Fenghour et al. [52] presented a componential equation to demonstrate the viscosity of CO2 as follows
μ c ( ρ , T ) = μ 0 α ( T ) + Δ μ ( ρ , T ) + μ c α ( ρ , T )
where   μ 0 α ( T ) the viscosity corresponding to the zero-density limit. Δ μ ( ρ , T ) contribution of viscosity at a specific density and specific temperature and μ c k ( ρ , T ) is the viscosity at a critical point. Fenghour et al. [52] also presented a details formulations of μ 0 α ( T ) , Δ μ ( ρ , T ) and μ c α ( ρ , T ) (see Fenghour et al. [52]). It is worth mentioning that the effect of μ c α ( ρ , T ) is negligible and we assume μ c α ( ρ , T ) = 0 .

Appendix B.2.2. Viscosity of Water (H2O)

Huber et al. [53] presented the viscosity of water as
μ ¯ α = μ ¯ 0 ( T ¯ ) . μ ¯ 1 ( T , ¯ ρ ¯ ) . μ ¯ 2 ( T , ¯ ρ ¯ )
where T ¯ = T T , ρ ¯ = ρ ρ ,   p ¯ = p p and μ ¯ α = μ α μ α . Additionally,   T , ρ ,   p and μ α are reference constant and Huber et al. [53] also presented their values. In Equation (A13), μ ¯ 0 ( T ¯ ) is the viscosity in the zero density limit, μ ¯ 1 ( T , ¯ ρ ¯ ) represents the residual contribution of viscosity and μ ¯ 2 ( T , ¯ ρ ¯ ) demonstrate the critical component of the viscosity (see also Huber et al. [53]).

Appendix B.3. Thermal Conductivity

Appendix B.3.1. Thermal Conductivity of Carbon Dioxide (CO2)

Vesovic et al. [54] presented the thermal conductivity of CO2 as follows
λ α ( ρ , T ) = λ 0 α ( T ) + Δ λ ( ρ , T ) + λ c α ( ρ , T )
where λ 0 α ( T )   is the thermal conductivity corresponding to the zero-density limit. Δ λ ( ρ , T ) represents the thermal conductivity contribution at a specific density and specific temperature and λ c α ( ρ , T ) is the thermal conductivity at a critical point (see also Vesovic et al. [54]). We neglect λ c α ( ρ , T ) in Equation (A14) because of its negligible effect.

Appendix B.3.2. Thermal Conductivity Water (H2O)

Huber et al. [55] demonstrated the thermal conductivity of H2O as follows
λ ¯ α = λ ¯ 0 ( T ¯ ) . λ ¯ 1 ( T , ¯ ρ ¯ ) . λ ¯ 2 ( T , ¯ ρ ¯ )
where λ ¯ 0 ( T ¯ ) is the thermal conductivity in the zero density limit, λ ¯ 1 ( T , ¯ ρ ¯ ) represents the residual contribution of thermal conductivity and λ ¯ 2 ( T , ¯ ρ ¯ ) demonstrate the critical component of the thermal conductivity (see also Huber et al. [55]).

Appendix B.4. Free Helmholtz Energy

From Wagner and Pruß [18], the Helmholtz free energy for thermodynamics properties is given by
f ( ρ , T ) R T = ϕ ( δ , τ ) = ϕ 0 ( δ , τ ) + ϕ r ( δ , τ )
where ϕ 0 ( δ , τ ) and ϕ r ( δ , τ ) represent the ideal gas component and the residual component, respectively. Additionally, δ = ρ ρ c , τ = T c T , while Wagner and Pruß [18] presented reference constants ρ c ,   T c and R . Wagner and Pruß [18] also demonstrated the formulation of ϕ 0 ( δ , τ ) and ϕ r ( δ , τ ) as follows
ϕ 0 = l n   δ + n 1 0 + n 2 0 τ + n 3 0 l n   τ + i = 4 8 n i 0 l n ( 1 e γ i 0 τ )
ϕ r = i = 1 k 1 n i δ d i τ t i + i = k 1 + 1 k 2 n i δ d i τ t i e δ c i + i = k 2 + 1 k 3 n i δ d i τ t i e α i ( δ ε i ) 2 β ( τ γ i ) 2 + i = k 3 + 1 k 4 n i Δ b i δ e C i ( δ 1 ) 2 D i ( τ 1 ) 2 Δ = ( ( 1 τ ) + A i ( ( δ 1 ) 2 ) 1 2 β i ) 2 + B i ( ( δ 1 ) 2 ) a i
Additionally, by differentiation of Equation (A17) and rearranging, we obtain thermodynamic functions (see Table A2).
Table A2. Relationship of thermodynamics functions to the Helmholtz free energy and their derivatives.
Table A2. Relationship of thermodynamics functions to the Helmholtz free energy and their derivatives.
Thermodynamic FunctionRelation
Pressure p ( δ , τ ) ρ R T = 1 + δ ϕ δ r
Entropy s ( δ , τ ) R = τ ( ϕ τ 0 + ϕ τ r ) ϕ 0 ϕ r
Enthalpy h ( δ , τ ) R T = 1 + τ ( ϕ τ 0 + ϕ τ r ) + δ ϕ δ r
Isochoric heat capacity c v ( δ , τ ) R = τ 2 ( ϕ τ τ 0 + ϕ τ τ r )
Isobaric heat capacity c p ( δ , τ ) R = τ 2 ( ϕ τ τ 0 + ϕ τ τ r ) + ( 1 + δ ϕ δ r δ τ ϕ δ τ r ) 2 1 + 2 δ ϕ δ r + δ 2 ϕ δ δ r
Note: ϕ δ r = ϕ r δ , ϕ τ r = ϕ r τ , ϕ τ 0 = ϕ 0 τ , ϕ δ δ r = ϕ r δ δ ,   ϕ δ τ r = ϕ r δ τ ,   ϕ τ τ r = ϕ r τ τ ,   ϕ τ τ 0 = ϕ 0 τ τ .

Appendix B.4.1. Free Helmholtz Energy of Carbon Dioxide (CO2)

For CO2, the limiting values of k are k 1 = 7 , k 2 = 34 , k 3 = 39 and k 4 = 42 (see Span and Wagner [56], p. 1544). Additionally, from Span and Wagner [56], we obtain the exponents and coefficient of Equations (A17) and (A18) for CO2.

Appendix B.4.2. Free Helmholtz Energy of Water (H2O)

For H2O, the limiting values of k are k 1 = 7 , k 2 = 51 , k 3 = 54 and k 4 = 56 (see also Wagner and Pruß [18], p. 429). Additionally, from Wagner and Pruß [18] we obtain the exponents and coefficient of Equations (A17) and (A18) for H2O.
Fluids properties equation for mixture and components are summarized in Table A3.
Table A3. Properties of mixture and components.
Table A3. Properties of mixture and components.
Mixture EquationComponent’s EquationReference
ρ = α = 1 N ρ α w α f ρ α = M α p z α ( p , T ) R T Equation (A10)
μ = α = 1 N μ α ( ρ , T ) w α f μ c ( ρ , T ) = μ 0 α ( T ) + Δ μ ( ρ , T ) + μ c α ( ρ , T )
μ ¯ w = μ ¯ 0 ( T ¯ ) . μ ¯ 1 ( T , ¯ ρ ¯ ) . μ ¯ 2 ( T , ¯ ρ ¯ )
See Appendix B, Equations
(A12) and (A13)
λ = α = 1 N λ α ( ρ , T ) w α f λ c ( ρ , T ) = λ 0 α ( T ) + Δ λ ( ρ , T ) + λ c α ( ρ , T )
λ ¯ α = λ ¯ 0 ( T ¯ ) . λ ¯ 1 ( T , ¯ ρ ¯ ) . λ ¯ 2 ( T , ¯ ρ ¯ )
See Appendix B, Equations
(A14) and (A15)
β p = α = 1 N v α v β p α f β p α f = 1 ρ α ( ρ α p ) T , w α f = 1 v v p | T , w α f Equation (A3)
β T = α = 1 N v α v β T α f β T α f = 1 ρ α ( ρ α p ) p , w α f = 1 v v p | p , w α f Equation (A3)
β w = α = 1 N v α v β w α f β w α f = 1 ρ α ( ρ f w α f ) p , T = 1 v v p | p , T Equation (A3)
Enthalpy, Entropy and Heat Capacity,See Table A3
Note: c = Carbon dioxide, w = Water.
Figure A1. Properties of CO2 at 307.15 K (ag; see also CoolProp Package: Bell et al. [57]).
Figure A1. Properties of CO2 at 307.15 K (ag; see also CoolProp Package: Bell et al. [57]).
Geotechnics 02 00043 g0a1aGeotechnics 02 00043 g0a1b
Figure A2. Properties of H2O at 307.15 K (ag; see also CoolProp Package: Bell, Wronski, Quoilin and Lemort [57]).
Figure A2. Properties of H2O at 307.15 K (ag; see also CoolProp Package: Bell, Wronski, Quoilin and Lemort [57]).
Geotechnics 02 00043 g0a2aGeotechnics 02 00043 g0a2b

Appendix C. Coupled Finite Element Formulations

We assumed that (i) the finite element domain ( Ω ) comprised of sub-domains ( Ω ) , (ii) elements number in each domain is n   1   to   n e , where n e denotes the maximum number of elements.
Ω = n = 1 n e Ω n .
Ω = n = 1 n e Ω n .
We also considered similar assumptions for the surface ( Γ )   and the surface section ( Γ ) to account the boundary conditions.
We present the expanded form of Equation (24) as follows
M p p e p ˙ e + M p t e T ˙ e + ( A p p e + K p p e ) p e + ( A p t e + K p t e ) T e = F p e
M p p e p ˙ e + M p t e T ˙ e + ( A p p e + K p p e ) p e + M t p e p ˙ e + M t t e T ˙ e + ( A t p e + K t p e ) p e + ( A t t e + K t t e ) T e = F t e
M w α f e w ˙ α f + ( A w α f e + K w α f e ) w α f = F w α f e
We obtained the global matrices as follows
M i j = n = 1 n e M i j e ; K i j = n = 1 n e K i j e ; A i j = n = 1 n e A i j e ; F i j = n = 1 n e F i j e
It is worth to note that in the global matrices, “ij” does not represent tensorial notations, but matrix components of solutions vectors X = { ψ , w α f   } T r , where   ψ = { p , T } T r (see also Section 3: Coupled finite element solutions). Details of components are presented as follows
M p p = Ω N T r ( n 1 ρ α ρ α p ) N d Ω
M p t = Ω N T r ( n 1 ρ α ρ α T ) N d Ω
M t p = Ω N T r ( n β T T ) N d Ω
M t t = Ω N T r ( ρ c p ) e f f N d Ω
M w k f = Ω N T r ( n ρ α ) N d Ω
K p p = Ω ( N ) T r ( k μ α ) N d Ω
K p t = Ω ( N ) T r ( 0 ) N d Ω
K t p = Ω ( N ) T r ( 0 ) N d Ω
K t t = Ω ( N ) T r ( λ e f f ) N d Ω
K w k f = Ω ( N ) T r ( n ρ α D ) N d Ω
A p p = Ω N T r ( n v 1 ρ α ρ α p ) N d Ω
A p t = Ω N T r ( n v 1 ρ α ρ α T ) N d Ω
A t p = Ω N T r ( n v ( β T T 1 ) ) N d Ω
A t t = Ω ( N ) T r ( n v ρ α c p f ) N d Ω
A w k f = Ω ( N ) T r ( n v ρ α ) N d Ω
F p = Ω N T r ( n v ( β T T 1 ) ) N d Ω
F t = Ω ( N ) T r ( n v ρ α c p f ) N d Ω
F w k f = Ω ( N ) T r ( n v ρ α ) N d Ω

References

  1. Bachu, S. Sequestration of CO2 in geological media: Criteria and approach for site selection in response to climate change. Energy Convers. Manag. 2000, 41, 953–970. [Google Scholar] [CrossRef]
  2. Volk, T. CO2 Rising: The World’s Greatest Environmental Challenge; MIT Press: Cambridge, MA, USA, 2008. [Google Scholar]
  3. EPA. Determination of the Mechanical Integrity of Injection Wells: United States Environmental Protection Agency, Region 5: Underground Injection Control (UIC) Branch, Regional Guidance #5; US Environmental Prptection Agency: Washington, DC, USA, 2008. [Google Scholar]
  4. Metz, B.; Davidson, O.; de Coninck, H.; Loos, M.; Meyer, L. IPCC Special Report on Carbon Dioxide Capture and Storage; Cambridge University Press: New York, NY, USA, 2005. [Google Scholar]
  5. EPA. Inventory of U.S. Greenhouse Gas Emissions and Sinks: 1990–2015, EPA 430-P-17-001; US Environmental Protection Agency: Washington, DC, USA, 2017. [Google Scholar]
  6. Eggleton, T. A Short Introduction to Climate Change; Cambridge University Press: New York, NY, USA, 2013. [Google Scholar]
  7. Bruckner, T.; Bashmakov, I.A.; Mulugetta, Y.; Chum, H.; de la Vega Navarro, A.; Edmonds, J.; Faaij, A.; Fungtammasan, B.; Garg, A.; Hertwich, E.; et al. Energy Systems. In Climate Change 2014: Mitigation of Climate Change. Contribution of Working Group III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: New York, NY, USA, 2014. [Google Scholar]
  8. Nordbotten, J.M.; Celia, M.A. Geological Storage of CO2 Modeling Approaches for Large-Scale Simulation; John Wiley and Sons: Hoboken, NJ, USA, 2012. [Google Scholar]
  9. NETL. Carbon Storage Atlas, 5th ed.; National Energy Technology Laboratory (NETL): Pittsburgh, PA, USA, 2015. [Google Scholar]
  10. Rutqvist, J.; Tsang, C.-F. A study of caprock hydromechanical changes associated with CO2-injection into a brine formation. Environ. Geol. 2002, 42, 296–305. [Google Scholar] [CrossRef]
  11. Benson, S.M.; Cole, D.R. CO2 Sequestration in Deep Sedimentary Formations. Elements 2008, 4, 325–331. [Google Scholar] [CrossRef]
  12. Yielding, G.; Øverland, J.A.; Byberg, G. Characterization of Fault Zones for Reservoir Modeling: An Example from the Gullfaks Field, Northern North Sea. AAPG Bull. 1999, 83, 925–951. [Google Scholar]
  13. Helmig, R. Multiphase Flow and Transport Processes in the Subsurface: A Contribution to the Modeling of Hydrosystems; Springer: Berlin/Heidelberg, Germany, 1997. [Google Scholar]
  14. Bear, J.; Bachmat, Y. Introduction to Modeling of Transport Phenomena in Porous Media; Kluwer Academic Publishers: London, UK, 1990. [Google Scholar]
  15. Vasylkivska, V.; Dilmore, R.; Lackey, G.; Zhang, Y.; King, S.; Bacon, D.; Chen, B.; Mansoor, K.; Harp, D. NRAP-open-IAM: A flexible open-source integrated-assessment-model for geologic carbon storage risk assessment and management. Environ. Model. Softw. 2021, 143, 105114. [Google Scholar] [CrossRef]
  16. Rutqvist, J. Status of the TOUGH-FLAC simulator and recent applications related to coupled fluid flow and crustal deformations. Comput. Geosci. 2011, 37, 739–750. [Google Scholar] [CrossRef] [Green Version]
  17. Beck, M.; Rinaldi, A.P.; Flemisch, B.; Class, H. Accuracy of fully coupled and sequential approaches for modeling hydro- and geomechanical processes. Comput. Geosci. 2020, 24, 1707–1723. [Google Scholar] [CrossRef]
  18. Wagner, W.; Pruß, A. The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use. J. Phys. Chem. Ref. Data 2002, 31, 387–535. [Google Scholar] [CrossRef] [Green Version]
  19. Peng, D.Y.; Robinson, D.B. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15, 59–64. [Google Scholar] [CrossRef]
  20. Kempka, T.; Kühn, M. Numerical simulations of CO2 arrival times and reservoir pressure coincide with observations from the Ketzin pilot site, Germany. Environ. Earth Sci. 2013, 70, 3675–3685. [Google Scholar] [CrossRef] [Green Version]
  21. Pruess, K. ECO2N: A TOUGH2 Fluid Property Module for Mixtures of Water, NaCl, and CO2, LBNL-57952; Lawrence Berkeley National Laboratory: Berkeley, CA, USA, 2005. [Google Scholar]
  22. Witherspoon, P.A.; Wang, J.S.Y.; Iwai, K.; Gale, J.E. Validity of Cubic Law for fluid flow in a deformable rock fracture. Water Resour. Res. 1980, 16, 1016–1024. [Google Scholar] [CrossRef]
  23. Wilson, R.K.; Aifantis, E.C. On the theory of consolidation with double porosity. Int. J. Eng. Sci. 1982, 20, 1009–1035. [Google Scholar] [CrossRef]
  24. Hassanizadeh, S.M. Modeling species transport by concentrated brine in aggregated porous media. Transp. Porous Media 1988, 3, 299–318. [Google Scholar] [CrossRef]
  25. Truesdell, C.; Noll, W. The Non-Linear Field Theories of Mechanics, 1st ed.; Springer: New York, NY, USA, 1965. [Google Scholar]
  26. Atkin, R.J.; Craine, R.E. Continuum theories of mixtures: Basic theory and historical development. Q. J. Mech. Appl. Math. 1976, 29, 209–244. [Google Scholar] [CrossRef]
  27. Bowen, R.M. Porous Media Model Formulations by the Theory of Mixtures. In Fundamentals of Transport Phenomena in Porous Media; Bear, J., Corapcioglu, M.Y., Eds.; Springer: Berlin/Heidelberg, Germany, 1984. [Google Scholar] [CrossRef]
  28. Islam, M.N.; Bunger, A.P.; Huerta, N.; Dilmore, R. Bentonite Extrusion into Near-Borehole Fracture. Geosciences 2019, 9, 495. [Google Scholar] [CrossRef] [Green Version]
  29. Islam, M.N.; Gnanendran, C.T. Elastic-Viscoplastic Model for Clays: Development, Validation, and Application. J. Eng. Mech. 2017, 143, 04017121. [Google Scholar] [CrossRef] [Green Version]
  30. Islam, M.N.; Gnanendran, C.T.; Massoudi, M. Finite Element Simulations of an Elasto-Viscoplastic Model for Clay. Geosciences 2019, 9, 145. [Google Scholar] [CrossRef] [Green Version]
  31. Gurtin, M.E. On the thermodynamics of elastic materials. J. Math. Anal. Appl. 1967, 18, 38–44. [Google Scholar] [CrossRef] [Green Version]
  32. Coleman, B.D.; Noll, W. The thermodynamics of elastic materials with heat conduction and viscosity. Arch. Ration. Mech. Anal. 1963, 13, 167–178. [Google Scholar] [CrossRef]
  33. Wood, W.L. Practical Time-Stepping Schemes; Oxford University Press: Oxford, UK, 1990. [Google Scholar]
  34. Islam, M.N.; Gnanendran, C.T. Non-associated Flow Rule based Elasto-Viscoplastic Model for Clay. Geosciences 2020, 10, 227. [Google Scholar] [CrossRef]
  35. Chen, Z.; Huan, G.; Ma, Y. Computational Methods for Multiphase Flows in Porous Media; Society for Industrial and Applied Mathematics: Philadelphia, Pennsylvania, 2006. [Google Scholar]
  36. Paniconi, C.; Putti, M. A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems. Water Resour. Res. 1994, 30, 3357–3374. [Google Scholar] [CrossRef]
  37. Kempka, T.; Kühn, M.; Class, H.; Frykman, P.; Kopp, A.; Nielsen, C.M.; Probst, P. Modelling of CO2 arrival time at Ketzin—Part I. Int. J. Greenh. Gas Control. 2010, 4, 1007–1015. [Google Scholar] [CrossRef]
  38. Martens, S.; Kempka, T.; Liebscher, A.; Lüth, S.; Möller, F.; Myrttinen, A.; Norden, B.; Schmidt-Hattenberger, C.; Zimmer, M.; Kühn, M.; et al. Europe’s longest-operating on-shore CO2 storage site at Ketzin, Germany: A progress report after three years of injection. Environ. Earth Sci. 2012, 67, 323–334. [Google Scholar] [CrossRef] [Green Version]
  39. Chen, F.; Wiese, B.; Zhou, Q.; Kowalsky, M.B.; Norden, B.; Kempka, T.; Birkholzer, J.T. Numerical modeling of the pumping tests at the Ketzin pilot site for CO2 injection: Model calibration and heterogeneity effects. Int. J. Greenh. Gas Control. 2014, 22, 200–212. [Google Scholar] [CrossRef]
  40. Liebscher, A.; Möller, F.; Bannach, A.; Köhler, S.; Wiebach, J.; Schmidt-Hattenberger, C.; Weiner, M.; Pretschner, C.; Ebert, K.; Zemke, J. Injection operation and operational pressure–temperature monitoring at the CO2 storage pilot site Ketzin, Germany—Design, results, recommendations. Int. J. Greenh. Gas Control. 2013, 15, 163–173. [Google Scholar] [CrossRef] [Green Version]
  41. Oldenburg, C.M. Joule-Thomson cooling due to CO2 injection into natural gas reservoirs. Energy Convers. Manag. 2007, 48, 1808–1815. [Google Scholar] [CrossRef] [Green Version]
  42. Ebigbo, A.; Class, H.; Helmig, R. CO2 leakage through an abandoned well: Problem-oriented benchmarks. Comput. Geosci. 2007, 11, 103–115. [Google Scholar] [CrossRef]
  43. Class, H.; Ebigbo, A.; Helmig, R.; Dahle, H.K.; Nordbotten, J.M.; Celia, M.A.; Audigane, P.; Darcis, M.; Ennis-King, J.; Fan, Y.; et al. A benchmark study on problems related to CO2 storage in geologic formations. Comput. Geosci. 2009, 13, 409. [Google Scholar] [CrossRef]
  44. Bear, J.; Berkowitz, B. Groundwater flow and pollution in fractured rock aquifers. In Developments in Hydraulic Engineering; Novak, P., Ed.; Elsevier: New York, NY, USA, 1987; Volume 4, pp. 175–238. [Google Scholar]
  45. Zienkiewicz, O.C.; Taylor, R.L.; Zhu, J.Z. The Finite Element Method: Its Basis and Fundamentals, 6th ed.; Elsevier: New York, NY, USA, 2005. [Google Scholar]
  46. de Borst, R. Computational Methods for Fractute in Porous Media: Isogeometric and Extended Finite Element Methods; Elsevier: New York, NY, USA, 2018. [Google Scholar]
  47. Mathias, S.; McElwaine, J.N.; Gluyas, J.G. Heat transport and pressure buildup during carbon dioxide injection into depleted gas reservoirs. J. Fluid Mech. 2014, 756, 89–109. [Google Scholar] [CrossRef] [Green Version]
  48. Islam, M.N.; Huerta, N.; Dilmore, R. Effect of Computational Schemes on Coupled Flow and Geo-Mechanical Modeling of CO2 Leakage through a Compromised Well. Computation 2020, 8, 98. [Google Scholar] [CrossRef]
  49. Ghanbarian, B.; Hunt, A.G.; Ewing, R.P.; Sahimi, M. Tortuosity in Porous Media: A Critical Review. Soil Sci. Soc. Am. J. 2013, 77, 1461–1477. [Google Scholar] [CrossRef]
  50. Millington, R.J.; Quirk, J.P. Permeability of porous solids. Trans. Faraday Soc. 1961, 57, 1200–1207. [Google Scholar] [CrossRef]
  51. Sengers, J.V.; Kayser, R.F.; Peters, C.J.; White, H.J., Jr. Equation of State for Fluids and Fluids Mixtures, Part I: Experimental Thermodynamics; Elsevier: New York, NY, USA, 2000. [Google Scholar]
  52. Fenghour, A.; Wakeham, W.A.; Vesovic, V. The Viscosity of Carbon Dioxide. J. Phys. Chem. Ref. Data 1998, 27, 31–44. [Google Scholar] [CrossRef]
  53. Huber, M.L.; Perkins, R.A.; Laesecke, A.; Friend, D.G.; Sengers, J.V.; Assael, M.J.; Metaxa, I.N.; Vogel, E.; Mareš, R.; Miyagawa, K. New International Formulation for the Viscosity of H2O. J. Phys. Chem. Ref. Data 2009, 38, 101–125. [Google Scholar] [CrossRef] [Green Version]
  54. Vesovic, V.; Wakeham, W.A.; Olchowy, G.A.; Sengers, J.V.; Watson, J.T.R.; Millat, J. The Transport Properties of Carbon Dioxide. J. Phys. Chem. Ref. Data 1990, 19, 763–808. [Google Scholar] [CrossRef]
  55. Huber, M.L.; Perkins, R.A.; Friend, D.G.; Sengers, J.V.; Assael, M.J.; Metaxa, I.N.; Miyagawa, K.; Hellmann, R.; Vogel, E. New International Formulation for the Thermal Conductivity of H2O. J. Phys. Chem. Ref. Data 2012, 41, 033102. [Google Scholar] [CrossRef] [Green Version]
  56. Span, R.; Wagner, W. A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa. J. Phys. Chem. Ref. Data 1996, 25, 1509–1596. [Google Scholar] [CrossRef] [Green Version]
  57. Bell, I.H.; Wronski, J.; Quoilin, S.; Lemort, V. Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and the Open-Source Thermophysical Property Library CoolProp. Ind. Eng. Chem. Res. 2014, 53, 2498–2508. [Google Scholar] [CrossRef]
Figure 1. Injection rate and cumulative CO2 injection (Liebscher et al. [40], Kempka and Kühn [20], data also available in https://github.com/ufz/ogs5-benchmarks).
Figure 1. Injection rate and cumulative CO2 injection (Liebscher et al. [40], Kempka and Kühn [20], data also available in https://github.com/ufz/ogs5-benchmarks).
Geotechnics 02 00043 g001
Figure 2. Comparison of measured and predicted pressure response at Ktzin-201 (Kempka and Kühn [20], data also available in https://github.com/ufz/ogs5-benchmarks).
Figure 2. Comparison of measured and predicted pressure response at Ktzin-201 (Kempka and Kühn [20], data also available in https://github.com/ufz/ogs5-benchmarks).
Geotechnics 02 00043 g002
Figure 3. Finite element geometry of the faulted reservoir.
Figure 3. Finite element geometry of the faulted reservoir.
Geotechnics 02 00043 g003
Figure 4. Injection rate effect on the pressure distribution in observation points (ac).
Figure 4. Injection rate effect on the pressure distribution in observation points (ac).
Geotechnics 02 00043 g004
Figure 5. Injection rate effect on CO2 mole fraction distribution in observation points (ac).
Figure 5. Injection rate effect on CO2 mole fraction distribution in observation points (ac).
Geotechnics 02 00043 g005
Figure 6. Spatial distribution after 300 days of injection of pressure, density, CO2, viscosity, specific heat capacity, thermal conductivity and temperature for thermal gradient = 0.02, injection rate = 0.41319 m3/s, α l = 5.0 m and α t = 0.5   m and fracture aperture = 100 µm.
Figure 6. Spatial distribution after 300 days of injection of pressure, density, CO2, viscosity, specific heat capacity, thermal conductivity and temperature for thermal gradient = 0.02, injection rate = 0.41319 m3/s, α l = 5.0 m and α t = 0.5   m and fracture aperture = 100 µm.
Geotechnics 02 00043 g006
Figure 7. Thermal effect on the distribution of pressure in different observation points (ac).
Figure 7. Thermal effect on the distribution of pressure in different observation points (ac).
Geotechnics 02 00043 g007
Figure 8. Thermal effect on the distribution of CO2 mole fraction in different observation points (ac).
Figure 8. Thermal effect on the distribution of CO2 mole fraction in different observation points (ac).
Geotechnics 02 00043 g008
Figure 9. Distribution of pressure, density, CO2, viscosity, specific heat capacity, thermal conductivity, and temperature after 300 days of injection for the case with thermal gradient = 0.03, injection rate = 0.41319 m3/s, α l = 5 and α t = 0.5 and fracture aperture = 100 µm.
Figure 9. Distribution of pressure, density, CO2, viscosity, specific heat capacity, thermal conductivity, and temperature after 300 days of injection for the case with thermal gradient = 0.03, injection rate = 0.41319 m3/s, α l = 5 and α t = 0.5 and fracture aperture = 100 µm.
Geotechnics 02 00043 g009aGeotechnics 02 00043 g009b
Figure 10. Distribution of pressure, density, CO2, viscosity, specific heat capacity, thermal conductivity, and temperature after 300 days of injection for the case with the isothermal condition, injection rate = 0.41319 m3/s, α l = 5 and α t = 0.5 and fracture aperture = 100 µm.
Figure 10. Distribution of pressure, density, CO2, viscosity, specific heat capacity, thermal conductivity, and temperature after 300 days of injection for the case with the isothermal condition, injection rate = 0.41319 m3/s, α l = 5 and α t = 0.5 and fracture aperture = 100 µm.
Geotechnics 02 00043 g010
Figure 11. The pressure distribution (ac) in observation point 1.
Figure 11. The pressure distribution (ac) in observation point 1.
Geotechnics 02 00043 g011
Figure 12. The pressure distribution (ac) in observation point 2.
Figure 12. The pressure distribution (ac) in observation point 2.
Geotechnics 02 00043 g012
Figure 13. The pressure distribution (ac) in observation point 3.
Figure 13. The pressure distribution (ac) in observation point 3.
Geotechnics 02 00043 g013
Figure 14. The CO2 saturation at observation point 1 (ac).
Figure 14. The CO2 saturation at observation point 1 (ac).
Geotechnics 02 00043 g014
Figure 15. The CO2 saturation at observation point 2 (ac).
Figure 15. The CO2 saturation at observation point 2 (ac).
Geotechnics 02 00043 g015
Figure 16. The CO2 saturation at observation point 3 (ac).
Figure 16. The CO2 saturation at observation point 3 (ac).
Geotechnics 02 00043 g016
Figure 17. Distribution of pressure, density, CO2, viscosity, specific heat capacity, thermal conductivity, and temperature after 300 days of injection for the case with thermal gradient = 0.03, injection rate = 0.41319 m3/s, α l = 5 and α t = 0.5 and fracture aperture = 300 µm.
Figure 17. Distribution of pressure, density, CO2, viscosity, specific heat capacity, thermal conductivity, and temperature after 300 days of injection for the case with thermal gradient = 0.03, injection rate = 0.41319 m3/s, α l = 5 and α t = 0.5 and fracture aperture = 300 µm.
Geotechnics 02 00043 g017
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Islam, M.N. Finite Element Simulations of Fluids Leakage through the Faulted Reservoir. Geotechnics 2022, 2, 908-934. https://doi.org/10.3390/geotechnics2040043

AMA Style

Islam MN. Finite Element Simulations of Fluids Leakage through the Faulted Reservoir. Geotechnics. 2022; 2(4):908-934. https://doi.org/10.3390/geotechnics2040043

Chicago/Turabian Style

Islam, Mohammad Nurul. 2022. "Finite Element Simulations of Fluids Leakage through the Faulted Reservoir" Geotechnics 2, no. 4: 908-934. https://doi.org/10.3390/geotechnics2040043

APA Style

Islam, M. N. (2022). Finite Element Simulations of Fluids Leakage through the Faulted Reservoir. Geotechnics, 2(4), 908-934. https://doi.org/10.3390/geotechnics2040043

Article Metrics

Back to TopTop