CO2 storage in deep saline aquifers: evaluation of geomechanical risks using integrated modeling workflow

CO2 injection into a saline aquifer crossed by a tectonic fault is studied with coupled fluid mechanics - geomechanics modeling. The simulation approach is based on coupling of the MUFITS reservoir simulator and the FLAC3D mechanical simulator via an in-house API (i.e., an algorithm for data transfer between simulators). MUFITS simulates the non-isothermal multiphase flow of CO2 and brine in rock formation accounting for phase transitions and thermal effects. The modeling workflow is sequential, so that hydrodynamical simulations are carried out at a certain time interval, after which pressure, temperature, and density distributions are passed to FLAC3D, which calculates the equilibrium mechanical state. Computed deformations and stresses are utilized to update the porosity and permeability fields for the subsequent hydrodynamic modeling. In particular, we focus on the tectonic fault and its behavior during CO2 injection. We distinguish the damage zone and core inside the fault and derive the closure relations for their permeability alteration analytically. The coupled approach developed here is applied to simulate CO2 injection into synthetic and realistic reservoir models. For the former one, we study the effect of formation depth and presence of the tectonic stresses at the initial mechanical state, while for the latter, we consider different injection modes (bottomhole pressure). In each numerical experiment, we describe the evolution of the fault permeability due to the slip along its plane and the development of plastic deformations leading to the loss of reservoir integrity and CO2 leakage. Sensitivity analysis of the coupled model to realistic values of input parameters to assess the fault stability is carried out.


Introduction
Carbon dioxide (CO 2 ) is a greenhouse gas contained in the atmosphere.The disbalance of CO 2 in the atmosphere forms predominantly due to burning of fossil fuels as a result of human activities.Continuous growth in CO 2 concentration in the atmosphere leads to the greenhouse effect.Consequently, the mean temperature grows contributing to an increase in the frequency and severity of catastrophic climate events (Pörtner et al., 2022).Emission of carbon dioxide can be reduced by using low-carbon fuels and improving energy efficiency.However, application of these techniques is insufficient to reduce the concentration of CO 2 significantly and prevent mean temperature growth (IEA, 2020).Resolving this issue requires the development and application of CO 2 sequestration technologies (Kazemifar, 2022).They include carbon dioxide capture at emission sources, its liquefaction, and transportation to fields, where CO 2 is injected into deep saline aquifers or depleted oil/gas formations in the supercritical state (Bickle, 2009).During the injection stage, CO 2 displaces pore fluid (e.g., natural gas or water) and forms a plume propagating due to the pressure difference between the injector and far-field pore pressure.In addition, the buoyancy force affects the plume dynamics since carbon dioxide is usually lighter as compared to pore fluids (Bryant et al., 2008); therefore, the target reservoir for secure CO 2 storage has to be covered by a low-permeable caprock.Flow of CO 2 in the reservoir rock is accompanied by various phase transitions and chemical reactions including dissolution of carbon dioxide in brine as well as water vapor in CO 2 (Huppert and Neufeld, 2014), precipitation of insoluble salts due to interactions of carbonated water with minerals composing the rock (de Coninck and Benson, 2014;De Silva et al., 2015;Cai et al., 2021).
Development of an efficient technology of CO 2 storage in underground formations requires solving a number of different problems, namely, (i) evaluation of the formation capacity or the maximum CO 2 volume that can be injected into the geological formation (Bachu et al., 2007;Bradshaw et al., 2007;Zhou et al., 2008;De Silva and Ranjith, 2012); (ii) estimation of the reservoir injectivity, which is the maximum injection rate of carbon dioxide based on operating parameters, rock permeability, properties of CO 2 and pore fluid (Rutqvist et al., 2007;Stauffer et al., 2009;Burton et al., 2009;Mathias et al., 2013); (iii) assessment of geomechanical effects leading to undesired events (Lucier and Zoback, 2008;Newmark et al., 2010;Zoback and Gorelick, 2012;Rutqvist, 2012;White and Foxall, 2016;Gholami et al., 2021).The current work is devoted to the latter problem, namely, evaluation of geomechanical risks of CO 2 injection into underground formations.
Four types of undesirable phenomena related to geomechanical effects are identified (Hawkes et al., 2005;Pawar et al., 2015;Pan et al., 2016): (i) loss of CO 2 storage integrity and leakage of carbon dioxide out of the target aquifer to the upper layers, (ii) activation of tectonic faults and fractures intersecting the storage reservoir near the overpressured zones, (iii) development and uncontrolled growth of hydraulic fractures due to overpressured and cooled zone around the injector; (iv) deformation of the Earth surface located above the storage area.These events can contribute to the seismic activity (earthquakes) in the neighbor regions, pollution of fresh water aquifers due to leakage of CO 2 or pore brine out of the target layer as well as damage to surface infrastructure.When the reservoir is intersected by a tectonic fault, and the CO 2 plume at high pore pressure reaches it, the effective normal stress at the fault plane declines leading to fault activation (Streit and Hillis, 2004;Konstantinovskaya et al., 2020).The slip promotes not only seismic activity and earthquakes but also the enhancement of the fault zone permeability (Rinaldi et al., 2014(Rinaldi et al., , 2015;;Guglielmi et al., 2017Guglielmi et al., , 2021)), which can lead to opening of the major crack in the fault core and natural fractures in the damage zone.Consequently, the tectonic fault forms a conduit along which CO 2 can flow out of the target reservoir.Summarizing the outlined geomechanical risks, we would like to highlight that it is crucial to model accurately CO 2 storage in underground reservoir during the planning stage.Mathematical modeling can be carried out to estimate the favorable injection regimes to prevent undesirable geomechanical phenomena running in the target reservoir and surrounding layers.
The problem of CO 2 injection into deep geological formation includes coupled thermal-hydraulic-mechanical processes.Two types of mathematical algorithms are developed to solve the set of governing equations: fully coupled and sequentially coupled (Dean et al., 2006;Kim, 2010;Kim et al., 2011;Ferronato et al., 2010;Rutqvist, 2011).The former approach involves simultaneous solution of the thermal-hydraulic and mechanical sets of governing equations, which is computationally intensive.In the latter method, the governing equations of the sub-problems are solved sequentially, and the required data is passed in between solvers via the coupling algorithm.The sequentially coupled approach is less computationally expensive as compared to fully coupled one and it is more flexible in terms of calculation algorithm as it allows one to (i) choose the degree of coupling in between thermal-hydraulic and mechanical models, namely, every iteration at the current time step (iterative coupling), a single coupling procedure at each time step (non-iterative coupling); (ii) the coupling frequency (time intervals, when coupling is applied); (iii) to utilize different spatial domains and time-stepping algorithms in each sub-problem; (iv) to utilize the existing open-source codes and commercial simulators as hydrodynamical and/or mechanical solvers.If the numerical solution converges, both approaches provide identical results of simulations as discussed by Kim et al. (2011).
Besides FLAC3D, the TOUGH family codes were coupled with other geomechanical simulators and packages, in particular, the open-source library PyLith (Aagaard et al., 2013) as demonstrated by Blanco-Martín et al. (2022).Various TOUGH-based geomechanical models are also summarized in a review paper (Rutqvist, 2017) and references therein.Rutqvist (2011) provided an overview of reservoir and mechanical simulators applied for modeling of coupled thermal-hydraulic-mechanical processes in geological formations.We would like to mention that certain software packages, namely, CMG (module GEM, CMG (2018)) and CODE-BRIGHT (Olivella et al., 1994(Olivella et al., , 1996) ) allow simulating CO 2 storage accounting for a non-isothermal multiphase flow of carbon dioxide and pore fluids accompanied by geomechanical effects within the frame of the single simulator (Vilarrasa et al., 2010;Jahandideh et al., 2021;Zheng et al., 2021).
In the current study, we develop a coupled hydrogeomechanical model of CO 2 storage in a saline aquifer intersected by a single tectonic fault.The model is based on the MUFITS academic reservoir simulator (Afanasyev, 2020), the FLAC3D commercial mechanical simulator (Itasca, 1997), and our in-house API, i.e., an algorithm for data transfer between simulators.Furthermore, we propose an analytical model describing permeability evolution in the tectonic fault zone of rock formation and embed it into the coupled modeling workflow.The closure relations include several parameters, which can be evaluated in lab experiments and via solving the inverse problem.The latter one is tuning of input parameters to approximate field observations by the results of simulations using the coupled hydro-geomechanical model.Our main goal is to study the fault behavior during CO 2 injection at different bottomhole pressure values, formation configurations, stress conditions, and reservoir properties to evaluate the associated geomechanical risks including fault activation and CO 2 leakage out of the target aquifer.
We organize the paper in the following way.Section 2 describes the coupled hydro-geomechanical model based on the MUFITS reservoir simulator, the FLAC3D mechanical simulator, and the algorithm executing the data transfer between simulators.In Section 3, we derive the analytical relations governing the permeability evolution in the damage zone and the core of the tectonic fault.Section 4 presents the results of CO 2 storage simulations based on synthetic and realistic reservoir models supplemented by the discussion with the focus on geomechanical risk assessment.Finally, we summarize the findings of the paper in Section 5.

Modeling approach
We develop a coupled hydro-geomechanical model to evaluate the geomechanical risks associated with CO 2 sequestration in a deep saline aquifer.The model is based on the MUFITS freely distributed reservoir simulator (Afanasyev (2020), MUFITS -Reservoir Simulation Software) and the FLAC3D commercial mechanical simulator (Itasca, 1997).

Hydrodynamical model
We employ the MUFITS reservoir simulator (Afanasyev et al., 2016;Afanasyev and Vedeneeva, 2021) for modeling the Darcy flow caused by the injection of CO 2 .We set the software to account for the dissolution of CO 2 in brine and the presence of the water vapor in the gas phase.The brine salinity is characterized by the NaCl concentration.Furthermore, we account for the temperature changes caused by the Joule-Thomson effect in the supercritical CO 2 , the convective heat transfer and the heat conduction.Thus in our study, the non-isothermal flow of CO 2 -H 2 O-NaCl fluid is governed by the following equations (Aziz and Settari, 1979;Pruess and Garcia, 2002;Pruess et al., 2004;Fanchi, 2005): (5) where  is the porosity,  is the density,  () is the th component mass fraction in th phase,  is the phase saturation,  is the Darcy velocity, k = diag(  ,   ,   ) is the permeability tensor,   is the relative permeability of the th phase,  is the fluid viscosity, g is the gravity acceleration,  and ℎ are the specific energy and enthalpy,  is the heat conductivity of saturated porous medium,  is the temperature,   is the gas-water capillary pressure, and   is the specific heat capacity of rock.The subscripts ,  and  refer to the parameters of the gas, liquid (brine) and rock phases, while the subscript () denotes the parameters of the th component, where  = 1, 2 and 3 correspond to CO 2 , H 2 O and NaCl, respectively.Eq. ( 1) is the mass conservation equation for the th component, Eq. ( 2) is Darcy's law, and Eq. ( 3) is the energy conservation equation.Eqs. ( 1)-( 3) are supplemented by the closing relations for the relative permeability and phase pressures in Eq. ( 5) and saturations in Eq. ( 6).The changes in temperature governed by Eqs. ( 3), ( 4) can be significant in the domains characterized by rapid pressure variations.The thermal effects can be important due to the difference between the reservoir temperature and that of the injected CO 2 , and they can affect stresses and strains induced by the injection.Therefore, we track the temperature distribution by simulating a non-isothermal flow.
Eq. ( 4) shows schematically the relations used for modeling the fluid phase equilibria.We assume that NaCl is present only in brine, i.e., its concentration in gas is zero ( (3) = 0).According to Eq. ( 4), the parameters of the fluid including the phase densities and viscosities are parameterized as functions of the pressure, temperature, and the brine salinity  𝑙(3) .Generally, we follow the methodology of Spycher and Pruess (2005) and Pruess and Spycher (2007) for predicting the fluid properties and the mutual solubilities of CO 2 and H 2 O. Here, we avoid further complications of the hydrodynamic model that can also account for the halite precipitation near the injection well (Pruess et al., 2004).The precipitation is not observed in the simulation scenarios considered in this study.Therefore, we present a simplified version of the governing equations to keep the presentation short.
The relative permeability and capillary pressure curves in Eq. ( 5) are given by where   is the saturation of the liquid phase,   is the saturation of the gas phase,   is the irreducible saturation of brine,   is the critical gas saturation,   ,   are the exponents,  0 is the strength coefficient, correspondingly.We use the Corey (1954) curve for the relative permeability of gas, while the relative permeability of brine and capillary pressure are calculated by the model proposed by Van Genuchten (1980).
In the present study, we do not consider chemical processes occurring in the storage reservoir and their impact on rock filtration-storage properties.It is known that chemical reactions between CO 2 dissolved in formation water and reactive minerals composing the rock can contribute to the dissolution of primary minerals and precipitation of secondary minerals resulting in the alteration of porosity and permeability.We focus on modeling of geological storage of CO 2 in sandstone formations.For this rock type, chemical reactions proceed very slowly, typically over hundreds of years, so that the characteristic time scale of chemical processes is much larger as compared to that of CO 2 injection (decades) (Xu et al., 2004;Ranganathan et al., 2011;Liu et al., 2011).Moreover, the impact of dissolution and precipitation of minerals on rock porosity and permeability manifests itself in the vicinity of the injector, and the corresponding effect on CO 2 injection is relatively small (Shiraki and Dunn, 2000;Tarkowski et al., 2015;De Silva et al., 2015;Zhang et al., 2019).

Mechanical model
We utilize the FLAC3D simulator for computing the mechanical equilibrium state of the fluid-saturated formation in terms of stresses and deformations corresponding to predefined pore pressure, temperature, and density distributions.While the MUFITS hydrodynamics simulator solves the dynamic problem, the FLAC3D mechanical simulator deals with a static one in the framework of the quasi-static approximation.Two constitutive models implemented in FLAC3D are utilized in the present study, namely, linear elastic isotropic model based on Hooke's law and elastoplastic one based on the Drucker-Prager yield condition (Itasca, 1997).
FLAC3D calculates the distributions of stresses and deformations via the numerical solution of the system of governing equations formulated with respect to mechanical (stresses) and kinematic (strain increment, velocity) variables as follows: where  is the saturated rock density,  is the velocity distribution,  is the stress tensor,  is the infinitesimal strain tensor,  is the displacement distribution, ε is the infinitesimal strain rate tensor,  ′ =  +  is the effective stress tensor,  is the pore pressure,  is the Biot coefficient,  denotes material functions, Δ is the time increment.The strain increment, Δ, can be represented by the sum of elastic, plastic, and thermal parts.Equation (8) describes momentum conservation law, Eq. ( 9) are Cauchy relations, while Eq. ( 10) are constitutive relations.FLAC3D approximates the deformable solid medium by elementary tetrahedrals, and their behavior is governed by Eqs. ( 8)-( 10) in accordance with the applied forces and boundary conditions.The system of equations is solved for the specified geometry and material properties at prescribed boundary and initial conditions.Below we describe several features of numerical solution to governing equations as implemented into FLAC3D: • The finite difference technique is applied.The first order time and spatial derivatives are represented by the finite differences based on the assumption that the variables alter linearly over a spatial segment and throughout a time interval; • The continuous medium is replaced by an equivalent discrete one, in which all forces (external and internal) are evaluated at the nodes of the three-dimensional grid used to approximate the deformable medium; • The dynamic solution approach is used.The inertial terms in the equation of motion are utilized as an indicator for the asymptotic approximation of the system mechanical equilibrium state.
Within the described framework, the motion equation for the continuous medium is transformed into the discrete form of Newton law formulated at the grid nodes.The system of ordinary differential equations is solved numerically using an explicit finite-difference time advance scheme.The definition of the strain rate tensor through the velocities at nodes includes the spatial derivatives.For each time step, the calculation procedure is as follows: 1. calculation of the updated deformations based on the velocities approximated at grid nodes; 2. computation of the stresses using the deformations, stresses at the previous time moment, and constitutive relations; 3. update the velocities and displacements based on the motion equation.
The sequence 1-3 is repeated at each internal FLAC step, at which the maximum unbalanced force is evaluated at the grid nodes.When the force becomes less as compared to the tolerance value, the mechanical system is assumed to be in the equilibrium state.When the unbalanced force reaches a constant non-zero value, it means that the entire system or its part is in the steady-state plastic flow.Calculations can be interrupted at any FLAC step to analyze the solution behavior.
For the convenience, we utilize the thermoelastic model, in which stresses and deformations are linked with the temperature  total distribution only.Here, the temperature equals to the sum of the actual temperature  actual and pseudo-temperature  pseudo .The latter one is introduced to describe the effect of the pore pressure  , and it can be determined from the analogy of poroelastic and thermoelastic problems: where  is the thermal-expansion coefficient, and  is the bulk modulus.Hence, stresses and deformations corresponding to the pressure  and temperature  actual fields can be calculated by the thermoelastic model where the formation is heated up to the temperature  total =  actual +  pseudo .In the current work, the Biot coefficient  is set to 1.

Coupling algorithm
The governing equations implemented in the MUFITS reservoir simulator and the FLAC3D mechanical simulator are solved sequentially, and the coupling parameters are transferred between simulators at certain time instants.In the current study, we consider identical spatial meshes for hydrodynamical and mechanical simulations and implement the in-house algorithm for the data transfer between the simulators using the approach proposed by Rutqvist et al. (2002).
Coupling procedure is summarized in Fig. 1.During time interval  ∈ [ −1 ,   ], MUFITS carries out the simulation of the multiphase Darcy non-isothermal flow at the current rock porosity (,  −1 ) and permeability (,  −1 ) distributions providing the pore pressure  (,   ), temperature  (,   ), and saturated rock density (,   ) fields.The later parameter is calculated as follows: where   is the liquid phase density (water with/without dissolved CO 2 ),   is the gas phase density (CO 2 with/without dissolved water vapor),   is the liquid phase saturation, and   is the rock density.Then, the hydrodynamical simulation is paused, and the calculated fields  (,   ),  (,   ), (,   ) are passed to the FLAC3D mechanical simulator.We should note that pore pressure and temperature are approximated in the centers of the mesh cells in MUFITS.However, in FLAC3D, temperature  total is approximated in the mesh nodes so that an interpolation procedure is applied.FLAC3D computes the components of the stress tensor (,   ) and infinitesimal strain tensor (,   ), corresponding to the mechanical equilibrium at the pore pressure  (,   ), temperature  (,   ), and density (,   ) fields based on the embedded geological model of the formation.
Coupling algorithm updates the reservoir porosity field (,   ) using the total volumetric strain distribution   (,   ) = tr (,   ) (both plastic and elastic), i.e., the trace of the strain tensor, calculated by FLAC3D and according to the relation (Chin et al., 2000): where  and  0 are the current and initial porosity values.
The permeability field (,   ) is found using the power-law relation (van Golf-Racht, 1982): where  and  0 denote the current and initial permeability values.The power-law exponent varies typically between 3 and 8 (Yehya et al., 2018), and we take  = 5 for the calculations.
In the current study, we pay particular attention to the fault zone and permeability alteration there due to the slip along the fault plane and plastic deformations.Section 3 outlines the mathematical model governing the permeability alteration in fault zone, and the procedure of its implementation into the hydrodynamical model is presented in Appendix B.
As a result, we estimate the changes in the rock porosity and permeability corresponding to the pore pressure, temperature variations, and the tectonic fault state contributing to both elastic and plastic deformations.Next, updated distributions (,   ) and (,   ) are passed back to the MUFITS simulator (interpolation is not required since the values of the required mechanical parameters used for porosity and permeability estimation are defined in the centers of the mesh cells), where the hydrodynamic simulation continues for the next time interval  ∈ [  ,  +1 ] on which  and  are fixed and equal (,   ), (,   ), respectively.

Alteration of permeability of the rock with a tectonic fault due to variations in pore pressure
Injection of CO 2 into an underground formation is accompanied by a local change in pore pressure and temperature leading to the alteration of the stress state, which can contribute to the activation of a tectonic fault located at a certain distance from the well.One can note that pressure perturbation propagates faster than the CO 2 plume meaning that the activation of faults and fractures can occur not only in the vicinity of the injection well.Moreover, fault slip and rock deformations in the fault zone improve its permeability.Therefore, a coupled hydro-geomechanical model of CO 2 injection has to describe the permeability modification due to changes in pore pressure and temperature.In this section, we present mathematical sub-models implemented into the coupled model (Section 2), which allow us to describe the corresponding geomechanical effects and risks: (i) activation of the tectonic fault and slip along its plane due to variation in a stress state; (ii) alteration of the fault zone permeability due to opening of the major crack on asperities in the core and natural fractures in the damage zone; (iii) disintegration of CO 2 storage zone and carbon dioxide leakage to the upper collectors along the fault zone being a highly permeable conduit.

Effect of inelastic deformations along the principal fault slip line on permeability of damage rock zone
Fault zones determine many processes running in the Earth crust and affect its mechanical properties significantly.Tectonic faults are most active structural formations through which an energy exchange in between tectonic blocks is carried out (Rice et al., 2005;Kanamori and Rivera, 2006), they also play a significant role in underground fluid movement (Townend and Zoback, 2000).Fault-block structure of Earth crust and sedimentation layer is one of key factors determining the stress state of underground formation.

Fault structure
Fault structure usually includes relatively narrow zone of large deformations surrounded by the transitional zone of fractured rock usually named as damage zone (Wibberley et al., 2008).The damage zone is surrounded by a host rock (see Fig. 2).The damage zone width is determined by a set of parameters including the thickness of rock layer being deformed, fault length and the density of fractures.Quantitative evaluation of damage zone is usually carried out by determining the density of rock fractures as a function of distance to the fault slip line.According to existing studies, for rocks with low porosity, there is an exponential decrease in density of fractures with an increase in the distance to the fault slip line (Anders and Wiltschko, 1994;Vermilye and Scholz, 1998;Wilson et al., 2003;Mitchell and Faulkner, 2009).
For fault zones in high-porosity rocks, the distribution of fracture density in the vicinity of the fault core is not clear: exponential law is confirmed in some cases (Anders and Wiltschko, 1994), while in other cases no correlation can be established (Shipton and Cowie, 2001).
Correlations in between the width of damage zone and key fault parameters (fault slip, length, fault displacement, throw) were discussed by geologists for several decades (Wibberley et al., 2008;Faulkner et al., 2010).Existing correlations differ significantly, and the reason is that authors define the width of damage zone differently, for example, it can be (i) a length measured from one the sides of the fault core; (ii) a single measurement or mean of the measurements; (iii) maximum width of the zone confined by the damage zone envelope (Shipton et al., 2006).
The structure of rock damage zone and fault core are closely related with the permeability.Caine et al. (1996) identified four types of the fault zone permeability structure based on the analysis of studies (Chester and Logan, 1987;Forster and Evans, 1991;Moore and Vrolijk, 1992;Newman and Mitra, 1994).Depending on dimensions and structure elements of the fault core and damage zones it was suggested to consider faults with localized conduits, distributed conduits, combined conduit-barriers and localized barriers as shown in Fig. 3 (Caine et al., 1996).

Localization of inelastic deformations
Tectonic faults are zones of localized irreversible deformations so that their initiation can be considered as a result of bifurcation of deformation process developed due to rheological instability of Earth crust (Rudnicki and Rice, 1975;Garagash and Nikolaevskii, 1989).Development of instability is associated with fracturing and dilatancy of rock under applied stress as well as with the effect of pressure on inelastic deformations.
Laboratory experiments on rock samples showed that their deformation is accompanied by the development of existing microfractures and pores as well as initiation of new fractures leading to alteration of effective mechanical properties of rock (Nikolaevskii et al., 1978;Paterson and Wong, 2005).This process depends both on applied stress and interaction of fracture walls.Its distinctive feature is dilatancy, which is irreversible increase in rock volume due to increase in size of pores and fracture aperture (Fig. 4).The most intensive change in rock structure occurs in the vicinity of peak stress before initiation of microscopic fracture-like defects.

Internal instability in the framework of non-associated plastic deformation law
Inelastic deformation of rock is carried out due to slip along existing fractures and initiation of new defects.The process can be described using exiting laws of plastic deformation, e.g., in the Prandtl-Rice form, accounting for the dilatancy and interaction of fracture walls.
Generalization to Prandtl-Rice constitutive relations to a medium with internal friction and dilatancy is first formulated by Nikolaevskii (1971) in the form of dependence of deformation increment on stress increment, while its alternative form (stress increment in terms of deformation increment) is given by Rudnicki and Rice (1975).We provide these relations in Appendix A for convenience.These relations are used by Rudnicki and Rice (1975) to study localization of inelastic deformations.It was found that under the Drucker-Prager yield criterion (Drucker and Prager, 1952) the formation of ordered structure of fractures occurs at the critical value of plastic hardening modulus   expressed as In Eqs. ( 14) and ( 15),  is the shear stress intensity,  is the rock cohesion,  is the mean stress,  and  are shear modulus and Poisson's ratio, respectively, Λ is dilatancy coefficient;  is internal friction coefficient.Formed fractures are ordered along the axis of mean principal stress  2 at the angle  with respect to direction of principal confining stress  3 : where  1 >  2 >  3 are principal components of the deviatoric stress.

Evaluation of horizontal tectonic stress and dilatancy coefficient
We evaluate the dilatancy coefficient Λ and horizontal stress leading to formation of a fault inclined by angle  under certain parameters of rock strength, namely, cohesion  and internal friction coefficient .
The Drucker-Prager yield criterion is formulated for fluid-saturated porous medium as follows: where   is the pore pressure.We consider the elastic rock layer located at the depth ℎ under vertical stress  33 and horizontal stresses due to lateral rock repulsion   11 and   22 as well as tectonic stresses   11 ,   22 .Total horizontal stresses due to lateral repulsion in elastic fluid-saturated layer are expressed as follows (Eaton's solution, Eaton (1969)): Expression ( 19) was obtained in the absence of thermalinduced stresses in the rock formation.Temperature of the rock increases with an increase in the depth according to geothermal gradient, which depends on rock composition, thermal conductivity and density of heat flux.Usually geothermal gradient takes values in the range between 0.5 • C up to 20 • C with the average value of 3 • C for 100 m depth.
Constitutive relation for a heated elastic layer has the following form (Timoshenko and Goodier, 1970) where   is the temperature and  is the thermal expansion coefficient.
Following the derivation of expressions ( 19) we consider elastic half-space with temperature distribution along the depth   =   ( 3 ).In this case, the deformations are expressed as follows so that equilibrium equations  , = 0 are satisfied.Now stress components according to (20) are expressed as follows: Eaton's solution ( 19) and the Drucker-Prager yield criterion ( 17) can be generalized using expressions ( 22) to describe the stress state of a heated fluid-saturated medium as follows: where The obtained stress components (23) allow calculating the stress intensity  and mean stress  under the applied tectonic stress   11 and   22 (see expressions in Appendix A below Eq. ( 51) and Eq. ( 18)).Parameters  and  are substituted into limiting condition (24) and assuming that   22 =   11 we find the critical horizontal stress   11() , at which the horizontal layer turns into inelastic state.Next, assuming that the localization of deformations is formed at the stress   11() , consider expression ( 16).Substituting all the known parameters and the angle of fault inclination , we find the corresponding dilatancy coefficient Λ.

Evaluation of rock permeability in the damage zone
The calculated dilatancy coefficient allows one to determine the increase in permeability of rock damage zone in the vicinity of the tectonic fault due to plastic deformations along its plane.Consider the expression for permeability of the fracture network (Basniev et al., 1993): where  is the mean fracture aperture,   is the fracture porosity, which is the ratio of the volume of fractures to the total rock volume.These parameters are related with each other by the following expression: where  is the fracture intensity determined experimentally as the ratio of total fracture length to the rock cross-section area.
According to the definition of dilatancy coefficient Λ the following expression holds where Γ  is the intensity of plastic shear deformations (the second invariant of the shear plastic deformation).
We assume that the ratio of volume of fractures opened due to rock dilatancy to the geometric volume of the rock is equal to an increase in the rock volume due to inelastic deformations   .Therefore, according to Eq. ( 27), we can formulate the following expression for   : (28) Substituting ( 28) into (25) and using Eq. ( 26) we find Expression (29) allows one to determine an increase in the permeability of near fault damage zone in the presence of plastic deformations.Note that all the parameters required to use the obtained expression are determined via standard laboratory measurements of rock mechanical properties.
The permeability alteration model described above is embedded into the framework of coupled hydro-geomechanical model (Section 2).Plastic deformations due to change in stress state of the rock during CO 2 injection are calculated directly using the FLAC3D mechanical simulator with the help of the relation: where   are elastic deformations,   is the total volumetric deformation calculated using FLAC3D, Δ ′ is the change in mean effective stress between the current and initial rock state, and  is the bulk modulus.The fracture intensity  varies between 10 and 50 m −1 (Golf-Rakht, 1986), and in the simulations shown below we assume  = 30 m −1 , which corresponds to sandstone.

Effect of shear slip on permeability of fractures closed on natural asperities
Modeling of a rock as elastoplastic medium with internal friction and dilatancy allows determining deconsolidation of the fault zone due to shear deformations as described above (see Eq. ( 29)).This approach can be used to describe the alteration of the fault dynamic influence zone (damage zone) containing microfractures, while it does not allow calculating the aperture of the main crack in the fault core.
Major fracture opening due to shear deformations is a result of applied shear stress being larger as compared to the friction force and resistant force of interaction of natural asperities as shown in Fig. 5.  Barton and Choubey (1977).
In the framework of Barton-Bandis model (Barton and Choubey, 1977), fracture aperture due to slip   is expressed as follows: where   is the slip displacement, and   is the dynamic angle of dilatancy.Instead of the tangent of   , we utilize the dilatancy coefficient Λ = tg(  ).
Laboratory experiments on rock samples show that the dilatancy angle depends on the normal stress, limiting shear stress, friction coefficient at the walls and residual friction angle.Zhigulskii and Tikhotskii (2020) show that 3D geomechanical modeling allows evaluating mechanical and hydraulic fracture aperture based on the rock stress state parameters and shear deformations.
Consider the plane problem of linear fracture in an elastic orthotropic medium.The fracture of the length 2 is oriented along the principal axes of anisotropy and is loaded at infinity with horizontal  1 , vertical  3 and tangential  stresses (see Fig. 6a).Fracture walls are confined with the stress  3 and interact according to dry friction law.Friction coefficient in static state   is maximum.If the shear stress increases up to its critical value   =    3 + , where  is the cohesion, then the fracture walls start to move and friction coefficient drops to the value   .As a result, the shear stress decreases to the value   =    3 and under the applied stress Δ (see Fig. 6b) the fracture walls are displaced by   =  1 ( 1 , 0).We define the shear deformation following the study (Garagash and Osiptsov, 2021).Equilibrium conditions are formulated as follows The condition (35) can be satisfied by using the constitutive relations for transversal isotropic body at plane strain state  22 = 0: 11 =  11  11 +  13  33 ,  33 =  13  11 +  33  33 , (36) The stiffness tensor components for isotropic medium have the following form: Solving Eq. ( 39) can be solved using Fourier transformation (Nowacki, 1975), and the following expression for displacement along the fracture axis at | 1 | ≤  is obtained (Garagash and Osiptsov, 2021): where Finally substituting displacement  1 ( 1 , 0) into expression (31) we find the fracture aperture   : In the mechanical model, we utilize  *  value corresponding to the averaged displacement profile  1 ( 1 , 0) along the fracture  1 ∈ [−, ]: By using the fracture aperture determined by Eq. ( 42) we calculate the permeability of the main fracture in the fault core closed on natural asperities as follows: where   is the hydraulic aperture of the main fracture determined according to Barton and Choubey (1977) as follows: where   and   are in microns.In Eq. ( 44), JRC is the joint roughness coefficient determined in the laboratory experiments.We take JRC = 4 in our study.As it is reported in (Barton and Choubey, 1977), JRC varies in the range between 1 and 20 The model of fracture permeability closed on natural asperities is embedded into the coupling procedure described in Section 2 for the description of the activation of the fault core as follows.If the deformations in the mesh cells containing fracture core are pure elastic, and the shear stress   on the fault plane reaches the critical value: where the shear stress   and normal effective stress  ′  on the fault plane with inclination angle  provided plane strain conditions (see Fig. 7) are given by then the fracture aperture  *  is calculated using Eq. ( 42).Stress difference Δ is determined based on the stress state evaluated by FLAC3D: We set   = 0.1,   = 0.05,  = 0 in Eq. ( 45), (46) during simulations (see typical values of the friction angle and cohesion in the fault zone in (Ikari et al., 2011;Carpenter et al., 2015;Treffeisen, 2021)).For estimation of the prefactor before Δ in Eq. ( 42), we assume that: (i) the fracture length 2 equals the height of the mesh cell, which is close to 10 m, so that  ∼ 5 m; (ii) elastic modulus of the transversal isotropic body is given by Bessmertnykh and Dontsov (2018):  11 = 20 GPa,  12 = 6.8 GPa,  13 = 7.6 GPa,  33 = 13 GPa,  44 = 3 GPa; (iii) the dilatancy coefficient is set to Λ = 0.5.As a result, we obtain where Δ is in Pa.Note that the dilatancy coefficient Λ varies in between 0 and 1 (Alejano and Alonso, 2005).If plastic deformations are observed, then the mechanical fracture aperture   is calculated using dilatancy and shear deformations as follows: where  is the shear-strain increment (the square root of the second invariant of the strain increment deviator), and  is the mesh cell length in direction perpendicular to the fault.

Verification of the coupled hydro-geomechanical model
In the current section, we verify the implementation of the in-house algorithm performing coupling between the MUFITS reservoir simulator and the FLAC3D mechanical simulator via data transfer.Both simulators have been verified previously.Results of the benchmark tests of MUFITS are given in papers (Afanasyev et al., 2016;De Lucia et al., 2016;Afanasyev, 2017) and are available at the web-site of the simulator (Afanasyev, 2020).At the same time, various tests of the FLAC3D commercial simulator are provided in the manual.
We consider the transient fluid flow to a vertical well fully penetrating the infinite-acting aquifer with thickness ℎ.We introduce the cylindrical coordinate system (, , ) with an origin located at the bottom of the perforation interval.Fig. 8 shows the schematic representation of the model.Initially, the pore fluid pressure  0 is uniform, and the reservoir is under isotropic stress  0  .The vertical well produces water at a constant rate .The elastic porous medium is taken as homogeneous and isotropic with drained bulk modulus , shear modulus , Biot coefficient , Biot modulus , porosity , permeability .Pore fluid (water) is described by the viscosity  and bulk modulus   .We consider an incompressible solid constituent  = 1 so that  =   ∕.The vertical stress is assumed constant during the water production   (, ) =  0  , and horizontal strains,   ,   , are negligibly small as compared to   .
The solution to the formulated problem in terms of the spatial-temporal parameters (pore fluid pressure , stress tensor components   ,   ,   , and vertical displacement   ) can be derived analytically and is formulated as follows: ) , ) , where We solve the formulated problem numerically using two approaches: (i) developed coupled hydro-geomechanical model based on MUFITS and FLAC3D and (ii) FLAC3D (FLAC3D can simulate single phase flow of slightly compressible liquid in addition to the mechanical calculations).Subsequently, model (ii) is applied to verify the implementation of porosity and permeability alteration in the hydrodynamical model according to Eq. ( 12), ( 13).
Table 1 outlines the values of model parameters used in the simulations.Note that it is necessary to put the adjusted fluid modulus into the hydrodynamic model in order to preserve the real diffusivity (in the expression, we take into account the Biot coefficient value  = 1).We model the water production within 60 days, and the outer boundary of the reservoir located in the numerical models at a distance of 1 km from the producer is not reached by the pressure wave during the simulation so that the condition of the infinite-acting reservoir is satisfied.
Fig. 9 compares two numerical solutions computed by MUFITS+FLAC3D and FLAC3D with the analytical solution given by Eq. ( 49).We look at the distributions of pressure, vertical displacement, radial and vertical stresses over the coordinate interval  ≤ 500 m at different time moments.Since the vertical stress is constant over time, we depict the numerical solution for this parameter at the end of the simulation.One can conclude that the outcome of the proposed coupled hydro-geomechanical model based on MUFITS and FLAC3D matches acceptably with the analytical solution of the problem, and we make similar inference regarding the numerical model built on FLAC3D.
In the numerical experiment shown in Fig. 9, we check the data transfer from MUFITS to FLAC3D.However, we should also verify the reverse data flow between simulators (i.e., from FLAC3D to MUFITS).For that purpose, we amend the numerical models by adding the alteration of porosity and permeability in the hydrodynamical part after each mechanical calculation.We embed the following relations for porosity and permeability:  = 1 − (1 −  0 ) −50⋅  ,  =  0 (∕ 0 ) 8 , where  0 = 0.1,  = 1 mD, so that we modify artificially Eq. ( 12) to increase the impact of volumetric deformations on porosity, while the value of the power-law exponent in Eq. ( 13) is within the typical range.Moreover, in the current numerical experiment, it is not required to adjust the fluid modulus.As a result, for comparison purposes, we can demonstrate the offset of the numerical solution from an analytical one as described by Eq. ( 49), which corresponds to the storage coefficient  = ∕  .
Fig. 10 shows the obtained results.We demonstrate the distributions of pressure, vertical displacement, radial stress, and permeability over the coordinate interval  ≤ 100 m at different time instants (permeability evolution is shown within the area  ∈ [0, 500] m).The numerical solutions deviate from the analytical ones after a couple of days of water production, the discrepancy is observed at  = 5 days.We would like to stress that the analytical curves in Fig. 10 are not the solution to the problem under consideration, and they correspond to the water flow in an incompressible porous medium with constant porosity  0 = 0.1 and permeability  0 = 1 mD.Moreover, one can observe a satisfactory match between the results of simulations obtained via MUFITS+FLAC3D and FLAC3D.In other words, Fig. 10 confirms the correct implementation of the data transfer from FLAC3D to MUFITS performed via the porosity and permeability modification in the hydrodynamical model implemented in MUFITS based on deformations computed by FLAC3D.

Coupled hydro-geomechanical modeling of CO 2 sequestration in an aquifer on the example of the synthetic formation case
The current section presents the simulation results of CO 2 injection into the target aquifer intersected by a tectonic fault.We construct a synthetic model of two-dimensional multilayered formation.During the interpretation of the calculations, we focus on the development of undesired mechanical processes, namely, slip along the fault plane resulting in the crack opening on the asperities, the plastic deformations in the fault zone, target aquifer, and caprock, as well as carbon dioxide leakage along the fault zone towards the overlying collector.Multilayered reservoir structure with plane layers is typical of aquifers located at a flat terrain, in contrast to depleted gas reservoirs with a pronounced structural trap formed by a caprock.Moreover, we suppose that in the coupled hydrogeomechanical model it is sufficient to consider only the fault closest to the injection well to evaluate geomechanical risks linked with the fault stability.The risks associated with more distant faults are assumed to be less likely to occur.

Model description
Fig. 11 shows the schematic representation of the utilized synthetic model of a two-dimensional multilayered reservoir with the following geometrical parameters: the lateral size of the reservoir is 1500 m, while its height equals 600 m.Two cases of the reservoir depth are considered: the upper bound locates at a depth of (i) 600 m and (ii) 1700 m.Case 1 corresponds to the formation for CO 2 storage located at the minimal depth where the injected carbon dioxide remains in the supercritical state (Ringrose et al., 2013).In turn, Case 2 represents deeper saline aquifers suitable for CO 2 storage, which can be located at depths down to 1.5-2.5 km (Konstantinovskaya et al., 2023).
Carbon dioxide is injected into the target aquifer of thickness 100 m surrounded by the low permeable caprock and basement layers, each of 150 m height.These three layers are intersected by the tectonic fault.The fault zone thickness is 20 m, the inclination angle relative to the vertical direction is 15 • .The distance between the left edge of the formation and the fault at a depth corresponding to the middle of the storage aquifer is 550 m.The fault consists of the damage zone and core, the thickness of latter one is 1 m.The upper and basal aquifers are placed above and below the caprock and basement layers.The injector has a vertical completion coinciding with the left border of the reservoir.Before CO 2 injection, the formation is assumed to be water-saturated, the pressure distribution is hydrostatic (in the general case, it is computed from the phases distribution and gravitational-capillary equilibrium), and the temperature field corresponds to the geothermal gradient of 25 • C/km.FLAC3D computes the in-situ mechanical state of the reservoir using the prescribed pressure and temperature at the initial time instant and the geological model (distributions of density and elasticity modulus).The calculated initial displacements are set to zero, so that the subsequent deformations appearing due to CO 2 injection are measured from the in-situ state.The vertical well injects carbon dioxide at constant bottomhole pressure: (i) 150 bar and (ii) 300 bar.We choose the bottomhole pressure values relying on the minimal principal stresses  min observed at a depth of perforations in the first (i) and second (ii) cases as follows: bottomhole pressure should be lower than  min (e.g., by 10 %), to prevent the initiation of hydraulic fractures.In the framework of the problem formulation under consideration (see boundary conditions and permeability values below) pore pressure in the storage aquifer on the left hand of the fault increases relatively fast.Consequently, a constant flow rate control at the injection well leads to reaching the pressure limit after a short time period.We neglect this initial time span and utilize constant bottomhole pressure condition.The top and left edges of the reservoir are impermeable (solid black lines in Fig. 11).At the bottom and right boundaries, we specify a constant pore pressure equal to the initial one (dashed black lines in Fig. 11).At the left and right edges of the formation, we fix zero displacements along x-axis   = 0, while the displacements along x and z-axis are prohibited at the bottom boundary   =   = 0.In addition, we fix the displacement along z-axis at the right border   = 0.At the top boundary, we apply constant loading   corresponding to the lithostatic pressure created by a layer of thickness 600 m (Case 1) and 1700 m (Case 2) with a density 2400 kg/m 3 .
We set the specific heat capacity of rock   = 0.81 kJ / (kg ⋅ K) and heat conductivity of saturated porous medium  = 3 W / (m ⋅ K) for the entire domain.Pore water salinity is neglected.In our simulations, we utilize the same filtration-storage properties of the reservoir in Cases 1 and 2. It is done intentionally to simplify the comparison of the simulation results.Permeability value impacts significantly on the dynamics of CO 2 plume migration.In turn, in Case 2 the pressure distribution does not change significantly in the simulated zone when the permeability value of the storage reservoir decreases.There are examples of reservoirs located at the depth of 2 km and characterized by high values of porosity and permeability similar to that considered in the current section (Tawiah et al., 2020).
Note that the chosen permeability values of the damage zone and fault core are consistent with the experimental measurements and estimations provided in papers (Faulkner et al., 2003;Wibberley and Shimamoto, 2003;Scibek, 2020).The authors of these studies conclude that the fault core permeability is typically lower than that of the damage zone.As a result, the fault permeability along its plane is larger than the permeability in the normal direction.Fig. 12 shows the utilized relative permeabilities for the gas (black line) and liquid (red line) phases, as well as the capillary pressure curve (blue line).These parameters depend on the liquid phase saturation   according to Eqs. ( 7   Mechanical and flow properties of the synthetic reservoir model depicted in Fig. 11.The values of cohesion, angle of internal friction, and dilatancy corresponding to the upper aquifer, basement, and basal aquifer are absent, since these layers are considered as elastic.We set the values of the elastic modulus and strength parameters in the fault zone in such a way as to reproduce its complex structure, and one can find the description of this procedure in the main text. The spatial meshes in the MUFITS reservoir simulator and the FLAC3D mechanical simulator are identical.The dimensions of mesh cells out of the fault zone are Δ = 20 m, Δ = 10 m (Fig. 13) with slightly variation towards the fault, where the columns of cells become parallel to the fault plane.We apply the grid refinement in the fault zone, which allows reproducing its complex structure including the core surrounded by the damage zone.Seven columns of cells are used to approximated the fault zone, and their thickness decreases towards the fault center.The central column denotes the fault core and contains the main crack, which is initially healed.The remaining 6 columns (3 to the left and right of the fault core) describe the damage zone.Elastic modules (bulk modulus  and shear modulus ), cohesion , and angle of internal friction  decrease towards the fault core in accordance with the studies (Gudmundsson, 2004;Faulkner et al., 2006;Treffeisen, 2021).We assume that the values of these parameters at the fault core comprise 70% of those corresponding to the host rock at the considered depth (see typical values of the contrast in the mechanical properties between the host rock and fault core in Holdsworth (2004); Collettini et al. (2009);Treffeisen (2021)).The trend for variation of the dilatancy coefficient Λ is similar, but it increases towards the fault center so that its values at the fault core are 30% larger than that of the host rock.CO 2 injection is simulated for 30 years.During the first 5 years, we compute the mechanical equilibrium state every 6 months using the FLAC3D simulator, and the reservoir porosity and permeability are updated according to the current stress state and deformations.After that period, the FLAC3D simulator is called once a year for 25 years.

Modeling results
Case 1: target aquifer at 950m depth We start with discussion of the results of pure hydrodynamical modeling of CO 2 injection using the MUFITS simulator (see Fig. 14).After 30 years of injection, the CO 2 plume reaches the fault zone and crosses it (see Fig. 14a).On the right side of the fault, CO 2 flows along the upper part of the target aquifer and passes 400 m.Carbon dioxide also flows along the tectonic fault where the CO 2 plume uplifts at a distance of 50 m.From Fig. 14b, one can notice that the storage aquifer leftward to the fault zone contains the pressure plume.Increase in pore pressure (the difference in pore pressure values at the final and initial time instants) is homogeneous in the target aquifer due to the high contrast in permeability values corresponding to the CO 2 storage domain and surrounding layers, where the pressure disturbance gradually decreases.On the left hand of the fault, the pore pressure increase is observed in the caprock and basement layers only.Consequently, according to the hydrodynamic simulation, there is no leakage of CO 2 from the target layer to the upper aquifer.
In the framework of coupled simulations we consider two cases of the initial mechanical state:   =   (no tectonic stresses) and   ∕  ∼ 2.7 at a depth of the target aquifer.
x'We begin with the case of no tectonic stresses, and Fig. 15 presents the results of simulations.Carbon dioxide flows along the fault, and the CO 2 plume reaches the upper aquifer as it can be noted in Fig. 15a.The opening of the main crack facilitates this behavior.During CO 2 injection, the slip along the fault plane in an elastic mode occurs, and the fracture opens at natural asperities.After 30 years of CO 2 sequestration, the main crack opens along its entire length.However, the fracture aperture at the depth interval corresponding to the caprock and storage aquifer is smaller as compared to that at the basement layer.Plastic deformations do not develop in the reservoir so that the permeability increase in the fault zone is observed only in the direction parallel to the fault plane due to opening of the main fracture.The appearance of the conduit leads to CO 2 flow along the fault to the upper aquifer and the carbon dioxide leakage out of the disposal zone.Moreover, rightward to the fault zone, CO 2 flow occurs along the shorter distance compared to that obtained using pure hydrodynamical modeling (Fig. 14a).Fig. 15d demonstrates the mass of the injected CO 2 in the coupled model (red line) and hydrodynamical model (blue line), and the former one is higher.The pore pressure increment shown in Fig. 15b is similar to that obtained using hydrodynamical model (Fig. 14b).However, the pressure increase on the right side of the fault in the caprock and basement layers is more pronounced in the coupled model.Volumetric strain is maximal in the target aquifer leftward to the fault and declines towards the upper and basal aquifers in the caprock and basement layers (see Fig. 15c).
Next, we discuss the results of the coupled hydrogeomechanical modeling of CO 2 injection into the target aquifer with pronounced tectonic stresses at the initial state as shown in Fig. 16.The ratio   ∕  ∼ 2.7 is set to observe the development of plastic deformations in the fault zone.After 30 years of carbon dioxide injection, the main crack opens along the entire length, and the fracture aperture decreases towards the basement layer.Intense plastic deformations develop in the fault zone at the depth of the target aquifer so that the natural fractures open in the damage zone of the tectonic fault.Thereby, the permeabilities of the fault along and perpendicular to its plane are increased.As a result, we observe a significant CO 2 leakage into the upper aquifer and CO 2 flow in the target aquifer on the right side of the fault (see Fig. 16a).From Fig. 16d, it is clear that the injected mass of CO 2 in the coupled hydro-geomechanical model is much higher as compared to the hydrodynamic simulation.The distributions of the pore pressure increment and volumetric strain (Figs.16b and c) are similar to that obtained in the case of no tectonic stresses (Fig. 15) except for the domain of the upper aquifer where both parameters grow due to the leakage of carbon dioxide.Additionally, note the significant volumetric strain and pore pressure in the fault zone at a depth of the caprock layer.
Case 2: target aquifer at 2000m depth Now we discuss the results of the modeling of CO 2 injection into the target aquifer of the formation with the upper boundary located at the depth of 1700 m.We begin with the hydrodynamic simulation (see Fig. 17).At the end of the simulation period (30 years), the CO 2 plume reaches the fault zone and crosses it (see Fig. 17a).During the flow rightward to the fault, CO 2 passes 900 m reaching approximately the right boundary of the reservoir.We also observe CO 2 flow along the fault zone, and carbon dioxide rises at a distance of 150 m.The shape of the pore pressure plume demonstrated in Fig. 17b is similar to that obtained in the previous case of the reservoir depth of 950 m (Fig. 14b).The difference is that the pressure increase with respect to the initial distribution is higher due to larger bottomhole pressure value (300 bar versus 150 bar).Thus, the hydrodynamic computation demonstrates that the CO 2 plume almost reaches the upper aquifer.
Let us consider the results of the coupled simulations in the absence of tectonic stress at the initial mechanical state   =   (see Fig. 18).Plastic deformations do not develop in the reservoir, and the main crack opens along the entire length in the elastic mode.The formed conduit contributes to the leakage of carbon dioxide into the upper aquifer.The CO 2 plume crosses the fault zone and spreads along the target layer on the right side of the fault over a shorter distance as compared to that obtained using pure hydrodynamical calculations (Fig. 18a).Fig. 18d compares the mass of injected CO 2 in the coupled and hydrodynamical simulations, and the former one is higher due to the major crack opening on asperities.The shape of the pressure plume shown in Fig. 18b is similar to that obtained in the previous configuration of the formation with no tectonic stresses (Fig. 15b).The target aquifer exhibits large values of the volumetric strain (Fig. 18c).We also observe the reduced volumetric deformations in the vicinity of the injector.We attribute this to the temperature effect since the injected carbon dioxide is colder as compared to the reservoir water inside the target aquifer.In the present case, CO 2 cools the reservoir in the vicinity of injector by 40 degrees, while in the previous case (Fig. 15c), we do not observe noticeable influence of the non-isothermal flow on the mechanical equilibrium state of the formation due to small difference in between CO 2 and pore fluid temperatures.
Fig. 19 shows the results of coupled hydro-geomechanical modeling of CO 2 injection into the formation in the presence of pronounced tectonic stresses.The ratio of principal stresses at the target layer   ∕  ∼ 1.8 is set to facilitate the development of plastic deformations in the fault zone.Note the plastic deformations are formed at a lower value of the ratio   ∕  as compared to the previous case of target aquifer located at the depth of 950 m.After 30 years of  carbon dioxide injection, the main fracture opens along the entire length of the fault zone.We obtained smaller values of the crack aperture at the depth interval corresponding to the bottom segment of the storage aquifer.Moreover, we observe the development of the plastic deformations at the fault zone at the caprock layer and target aquifer.The opened major crack in the core and natural fractures in the damage zone of the fault contribute to the CO 2 leakage into the upper aquifer and CO 2 flow towards the right border of the reservoir (see Fig. 19a).The latter effect results in the larger volume of CO 2 plume rightward to the fault as compared to that obtained using the hydrodynamic model.The mass of injected CO 2 is substantially larger in the case of the coupled model (red curve in Fig. 19d) as compared to that obtained using the hydrodynamic model.The pore pressure in the plume differs from that obtained in the case of no tectonic stresses by a tangible pressure increase in the upper aquifer due to the carbon dioxide leakage (Fig. 19b).Volumetric deformation is maximal in the fault zone at the depth of the caprock layer (Fig. 19c).Still they are large in the target aquifer on the left side of the fault.The reduced values of the volumetric strain are attributed to the thermal effects.

Coupled hydro-geomechanical modeling of CO 2 sequestration in a real aquifer
In the current section, we show the results the modeling of CO 2 injection into an aquifer located in Volga-Ural oil province, Russia, using the proposed coupled hydrogeomechanical approach.We consider a slice of the reservoir sector that contains a tectonic fault.Similar to Section 4.2, the model is two-dimensional.For the model construction we use the field data collected from well logging, well test, seismic survey, and laboratory experiments.

Model description
Fig. 20 shows the schematic representation of the formation under consideration.Here, we illustrate its structure, geometrical characteristics, and the tectonic fault placement relative to the injector.
The reservoir has the length and height of 1500 m and 1350 m, respectively.The upper boundary of the formation is located at the depth of 1350 m.Carbon dioxide is injected into the upper and lower aquifers through a vertical well located at the left boundary of the reservoir; perforations are distributed along the interval  ∈ [−2300, −2100] m except for the 10 m thick caprock layer.From Fig. 20 one can observe a layer of anhydrite and a massive layer of salt above the upper aquifer.A tectonic fault starts at the salt layer and finishes at the basement crossing the anhydrite layer as well as aquifers and caprock.The fault is almost vertical with an inclination angle of 3 • with respect to the vertical direction.The distance between the injector and the fault equals 360 m at  = −2300 m corresponding to the middle of the lower aquifer.We embed the same fault structure into the coupled model as considered in Section 4.2, namely, the fault has thickness of 20 m, while its core thickness is set to 1 m.
Carbon dioxide is injected into the upper and lower aquifers at fixed bottomhole pressure.We consider two cases of injection with the following bottomhole pressures: (i) 350 bar (base case) and (ii) 600 bar (upper limit).The boundary conditions in the hydrodynamic and mechanical models are similar to the synthetic reservoir model.The difference is in the constant loading applied at the upper border of the formation: in the current model,   corresponds to the lithostatic pressure created by a layer of thickness 1350 m with a density 2400 kg/m 3 .
We describe the mechanical and flow properties of the formation in Table 3.The relative permeabilities and capillary pressure curve are the same as provided in Fig. 12.
Principal components of the tectonic strain tensor at the initial state are    = −10 −4 ,    = −3 ⋅ 10 −4 at a depth of the lower aquifer.Using the relations we compute the principal components of the tectonic stress tensor.In these relations, Young's modulus  and Poisson ratio correspond to the lower aquifer.Since we are interested in the tectonic stresses observed at the plane of the examined reservoir slice, we utilize the following expressions: where  is the angle between the principal x-direction and the slice plane (see Fig. 21).In the current case, the angle  equals 15 • .
In the mechanical reservoir model, we describe the rheology of all layers by the elastoplastic model with the Drucker-Prager yield condition.The spatial meshes in MU-FITS and FLAC3D are identical.The cell dimensions are Δ = 20 m and Δ = 10 m in the domain outside of the fault zone, where we utilize the same grid refinement in the fault zone as in the synthetic reservoir model (Section 4.2) to represent its complex structure.Carbon dioxide is injected for 30 years.

Modeling results
We start with the results of the base case.Here, carbon dioxide is injected into the upper and lower aquifers at a constant bottomhole pressure of 350 bar.Fig. 22 shows the gas saturation, pore pressure increment, and volumetric strain distributions at the end of the simulation period.
After 30 years of CO 2 injection, the carbon dioxide plume does not reach the fault zone (Fig. 22a).Its maximum lateral size and height are 250 m and 450 m, respectively.Large pore pressure increase is observed inside the domain occupied by the CO 2 plume (Fig. 22b).The pore pressure perturbations extend across the entire thickness of the reservoir and along the distance of 1 km from the injector in the lateral direction.Thus, the pore pressure plume dimensions

Table 3
Mechanical and flow properties of the realistic formation shown in Fig. 20); all layers are governed by the elastoplastic constitutive model; the distribution of the mechanical properties inside the fault zone is similar to that described in Section 4.2.are much larger as compared to that of CO 2 plume.We observe the large values of the volumetric strain in both aquifers and in the lower part of the salt deposit (Fig. 22c).Similar to the synthetic reservoir model (Section 4.2), the reduced values of the volumetric strain near the perforations are attributed to the cooling effect.Plastic deformations are not developed in the reservoir.The main fracture does not open so that the condition Eq. ( 45) is not satisfied along the entire crack.
In the current case, permeability can increase in the formation due to the volumetric deformations only according to Eq. ( 13).Since deformations are relatively small, the changes in porosity and permeability are also small.For example, permeability increase is less than 1%.Comparing dynamics of the mass of injected carbon dioxide in the case of coupled modeling and hydrodynamic simulation, we obtain a negligible difference between them.
Next, we move on to the results of the case in which the bottomhole pressure is fixed to 600 bar, which exceeds the minimum principal stress.Fig. 23 illustrates the distributions of the gas saturation, pore pressure increment, and volumetric strain after 30 years of CO 2 injection.
From Fig. 23a it is clear that the CO 2 plume reaches the fault zone, crosses it in both aquifers, and distributes partially along the damage zone and core of the fault.Rightward to the fault, CO 2 flows in the bottom part of the upper aquifer, and throughout the entire thickness of the lower aquifer, where the maximum lateral size of the CO 2 plume is about 550 m.The shape of the pore pressure plume shown in Fig. 23b is similar to the previous case, in which bottomhole pressure equals 350 bar (Fig. 22b): pore pressure diffuses over the entire reservoir thickness and along the domain of 1 km length in horizontal direction.Moreover, the volumetric strain field is qualitatively similar to that obtained in the previous case (Fig. 22c), while the deformations are larger by an about an order of magnitude due to the higher pore pressure (Fig. 23c).Plastic deformations do not develop in the fault zone, but we observe them in the vicinity of the perforation interval located in the upper aquifer and in the bottom part of the salt deposit near the left boundary of the formation.The former indicate the possible appearance of hydraulic fractures near the injector since the bottomhole pressure exceeds the minimal stress.The main crack does not open, and despite the increased value of the bottomhole pressure (600 bar versus 350 bar), we still obtain that the slip along the fault plane does not occur.The maximum permeability increase observed in the aquifers is about 2.5%, while this value in the damage zone of the fault reaches 4%.Due to the improvement of permeability in the target layers, the mass of the injected CO 2 is higher in the coupled simulation as compared to that obtained using the hydrodynamic simulation (red line compared to the blue line in Fig. 24).However, the difference in the injected mass at the end of the simulation period is insignificant.

Sensitivity analysis of the coupled hydro-geomechanical model
In the current section, we perform the sensitivity analysis of the coupled hydro-geomechanical model by varying mechanical and flow properties of the reservoir.We analyze the fault stability as well as the development of plastic deformations leading to the loss of integrity of the storage domain.The simulations are carried out for realistic parameters determining CO 2 storage in the aquifers of at the bottomhole pressure of 600 bar.
We assume that the mechanical properties of both aquifers and caprock layer between them can vary in certain ranges.In Section 4.3.2,we put into the FLAC3D simulator their average values (Table 3).However, one can suggest that the minimal values of the elastic modulus and strength parameters can facilitate the undesired mechanical effects related to the tectonic fault.We carry out the coupled simulations using the decreased values of the mechanical parameters outlined in Table 4.In Fig. 24, we compare the dynamics of the injected carbon dioxide mass during the second half of the simulation period obtained by two coupled models: with the average (red line) and reduced (green line) values of the mechanical properties (Table 4).We find that in the case of modified properties the injected mass is larger as compared to than obtained in the base case.The larger mass of the injected CO 2 is associated predominantly with the larger volumetric strains (see Fig. 25b) leading to the improvement of porosity and permeability (see Eqs. ( 12), ( 13)).In contrast to the base model, we do not observe the development of plastic deformations near the perforation interval in the upper aquifer, while they are localized to the lower left part of the salt layer only.The main fracture opens on asperities along the depth intervals corresponding to the lower aquifer and the top part of the upper aquifer.The CO 2 plume shape in the modified In Section 4.2, we describe how the mechanical properties (bulk modulus, shear modulus, cohesion, angle of internal friction, and dilatancy) are set in the fault zone.All parameters except for dilatancy decrease by 30% towards the fault core as compared to the corresponding parameters of the host rock at the considered depth, while the dilatancy grows by the same quantity.In the current numerical experiment for the sensitivity analysis, we carry out the coupled simulations of CO 2 injection assuming the alteration of the mechanical properties in the fault zone by 80% towards its center.Analyzing the results of simulations we conclude that there are plastic deformations developed in the fault zone along its entire length in addition to the similar domains near the perforation interval in the upper aquifer and in the bottom left part of the salt layer as in the base case.The main fracture opens along the entire length; however, its aperture is smaller as compared to the case of the modified coupled model with the reduced values of the mechanical properties.We find that the improvement of permeability in the fault zone due to the mechanical effects is negligible leading to the same shapes of the CO 2 plume and close values of the injected mass of carbon dioxide in the modified and base cases.The comparison of the solutions in terms of the injected CO 2 mass is shown in Fig. 24 (dashed blue line versus red line).
Finally, we conduct a coupled modeling of CO 2 injection into storage aquifers considering the fault zone with the uniform permeability set to 10 mD so that in the current experiment, we do not distinguish the damage and core zones in terms of the flow properties and present the fault zone as a conduit.The carbon dioxide injected mass obtained using the modified model exceeds that obtained in the base case (see 24, orange line compared to the red line).In the modified model, plastic deformations are developed in the same zones as in the base case, namely, near the perforation interval in the upper aquifer and in the bottom part of the salt layer close to the left border of the formation.Similar to the base case, the main fracture does not open.The increased permeability in the fault zone contributes to the flow of a larger volume of CO 2 in parallel and perpendicular directions to the main crack resulting in the slightly greater size of the carbon dioxide plume rightward to the fault as compared to the base case.

Summary and conclusions
In this paper, we developed a coupled hydro-geomechanical model for the simulation of CO 2 injection (storage) into a saline aquifer intersected by a tectonic fault.The simulation approach is based on the coupling of MUFITS reservoir simulator and FLAC3D mechanical simulator via the developed in-house algorithm performing the two-way coupling of the simulators: pressure, temperature, and density distributions are calculated in MUFITS and transferred to FLAC3D, while porosity and permeability fields calculated on the basis of deformations and stresses calculated in FLAC3D are passed to MUFITS and used in calculations during the next time step.Dynamics of rock filtration properties are calculated according to the novel mathematical model proposed in the current study, which describe the dependence of porosity and permeability inside the the host rock, damage zone and fault core.During the modeling of the CO 2 injection, MUFITS simulates a transient non-isothermal multiphase flow in a rock formation, while FLAC3D is applied to solve a quasi-static problem and computes the mechanical equilibrium of the reservoir at the fixed distribution of parameters (pore pressure, fluid saturations and temperature).
We verified the proposed coupled model by solving the problem of transient flow of slightly compressible fluid to a vertical well fully penetrating the infinite-acting reservoir.Results of numerical simulations are compared with an analytical solution in terms of distributions of pressure, stress tensor components, and vertical displacement preserving the true diffusivity in the hydrodynamic simulations.We also accounted for the alteration of porosity and permeability in the hydrodynamic model depending on the volumetric strain evaluated by the mechanical model.In this numerical experiment, we compared the results of coupled MUFITS+FLAC3D simulations against that carried out in FLAC3D only in terms of the parameters listed above, and a good match is obtained.
The developed coupled model is used to simulate the CO 2 storage in the framework of two-dimensional synthetic and realistic reservoir models.
In the former series of simulations, we examined two locations of the synthetic storage aquifer at a depth of 950 m and 2 km.For each configuration, we demonstrated the results of modeling at missing or pronounced tectonic stresses applied to the reservoir.We observed that in the absence of tectonic plastic deformations do not develop in the reservoir, and the major fracture in the fault core opens on asperities in the elastic mode contributing to an increase in the fault zone permeability along its plane and CO 2 leakage out of the target aquifer.When the tectonic stresses are pronounced, we found that the plastic deformations are developed in the fault zone in addition to the major crack opening.It results in the permeability increase in the directions parallel and perpendicular to the fault, CO 2 leakage into the upper aquifer, and considerably larger mass of the injected carbon dioxide as compared to that obtained in the case with no tectonic stresses.
The second set of simulations is carried out using the set of parameters corresponding to a real aquifer.We considered a vertical slice of the formation sector and analyzed the CO 2 injection with constant bottomhole pressure of 350 bar (base value) and 600 bar (upper limit).We determined the absence of undesirable mechanical effects in the base case.For an increased bottomhole pressure, we demonstrated that the fault remains stable while plastic deformations are developed in the vicinity of the perforation interval indicating the possible initiation of hydraulic fractures.We performed the sensitivity analysis of the coupled model to the input parameters describing the fault behavior by varying the mechanical properties of the storage layers and fault zone as well as the fault permeability.It is shown that the reduced values of the mechanical properties in the target layers contribute to an increase in the volumetric deformations (leading to an increase in porosity and permeability) and a partial opening of the main fracture.Strengthened contrast in the mechanical parameters of the host rock and fault core yields insignificant plastic deformations in the fault zone and the main fracture opening with a negligibly small aperture.Finally, an increased permeability in the fault zone results in a noticeable increase in the injected CO 2 mass.In the current section, we describe how parameters ,   , and   are combined in the hydrodynamical model.
We begin with the cells related to the damage zone (Fig. 27a).Each cell includes two interpenetrating isotropic continua, namely, host rock and network of natural fractures.We apply the equation describing the total permeability of the layered formation in the direction to the stratification: where   and ℎ  are the permeability and the thickness of each layer in the layered reservoir.As a result, we estimate the permeability of the cells located in the damage zone of the fault as follows: where   is plastic volumetric strain.Next, we move to the cells related to the fault core (Fig. 27b).Each of them is intersected by a main crack.In the derivations, we do not account for the fracture inclination assuming that the tectonic fault is approximately vertical and parallel to z-axis.Leftward and rightward to the fracture, the reservoir structure is similar to the damage zone, which is the host rock containing the network of the natural fractures.
The main crack opening impacts on the cell permeability in the vertical direction only.Based on Eq. ( 54), we derive the expressions for the total permeability along x and z-axis as follows: where   is the geometrical aperture of the main crack,  is the cell size along x-axis.Fig. 28 summarizes the implementation of models describing permeabilities of the damage zone and fault core into the hydrodynamic model.Note that the rock formation beyond the fault zone is considered undamaged, and its permeability is equal to that of the host rock.

Figure 1 :
Figure 1: Schematic representation of the coupling procedure between MUFITS and FLAC3D simulators applying for hydrogeomechanical simulations.

Figure 2 :
Figure 2: Conceptual representation of the fault structure, after Shipton and Cowie (2003); damage zone and fault core are shown by DZ and FC, respectively.

Figure 4 :
Figure 4: Diagram measured during triaxial stress loading of rock sample in laboratory conditions, after Jaeger et al. (2007).

Figure 6 :
Figure 6: Plane rock domain containing fracture (a) with the walls interacting according to dry friction law (b).

Figure 7 :
Figure 7: Stress components   and   at the fault plane inclined by angle  to the vertical direction under plane stress conditions.

Figure 8 :
Figure 8: Schematic picture of a vertical production well fully penetrating the infinite-acting aquifer.

Figure 9 :
Figure 9: Results of the analytical and numerical modeling of transient fluid flow to a vertical well fully penetrating the infiniteacting poroelastic reservoir in terms of pressure (a), vertical displacement (b), radial stress (c) and vertical stress (d) distributions; The analytical solution is shown by solid lines, while the numerical solutions are given by markers, MUFITS+FLAC3D by crosses and FLAC3D by circles; the solutions correspond to the following time instants:  = {1, 5, 15, 30, 60} day(s); zoomed domain in the vicinity of the well  ∈ [0, 100] m is shown in plots (a) -(c).

Figure 10 :
Figure 10: Results of the numerical modeling of transient axisymmetric fluid flow to a vertical well fully penetrating the infinite-acting reservoir accounting for the porosity and permeability alterations based on the mechanical calculations: pressure (a), vertical displacement (b), radial stress (c) and permeability (d) distributions; The numerical solutions are given by markers, MUFITS+FLAC3D by crosses and FLAC3D by circles; analytical solution corresponding to the pore fluid flow in the incompressible porous medium with the constant porosity and permeability (dashed lines); we demonstrate the solutions along the spatial domain  ∈ [0, 100] m ( ∈ [0, 500] m in plot d) at the following time instants:  = {1, 5, 15, 30, 60} day(s).

Figure 11 :
Figure 11: Synthetic reservoir model; solid and dashed black lines denote impermeable and constant pressure boundaries, respectively; the red arrow marks the perforation interval through which CO 2 is injected into the target aquifer; by red and purple colors we show core and damage zones in the fault domain.

Figure 12 :
Figure 12: Relative permeabilities and capillary pressure curves embedded into the hydrodynamic reservoir model.

Figure 13 :
Figure 13: Reservoir spatial discretization; the mesh structure in the fault zone is zoomed; colors highlight different layers and fault zone (see Table2).
Figure 13: Reservoir spatial discretization; the mesh structure in the fault zone is zoomed; colors highlight different layers and fault zone (see Table2).

Figure 14 :
Figure 14: Results of hydrodynamical modeling of CO 2 injection into the target aquifer with the upper boundary located at a depth of 600 m; plots (a) and (b) show the gas saturation and pore pressure increment (as compared to the initial hydrostatic distribution) fields at the end of the simulation period (30 years).

Figure 15 :
Figure 15: Results of coupled hydro-geomechanical modeling of CO 2 injection into the target aquifer with the upper boundary located at a depth of 600 m in the absence of tectonic stresses (  =   ); plots (a) -(c) show the gas saturation, pore pressure increment (as compared to the initial hydrostatic distribution), and volumetric strain fields at the end of the simulation period (30 years), respectively; plot (d) depicts the dynamics of injected mass of carbon dioxide in pure hydrodynamical model (red curve) and coupled model (blue curve).

Figure 16 :
Figure 16: Results of coupled hydro-geomechanical modeling of CO 2 injection into the target aquifer with the upper boundary located at a depth of 600 m at the pronounced tectonic stresses (  ∕  ∼ 2.7 in the target aquifer); plots (a) -(c) show the gas saturation, pore pressure increment (compared to the initial hydrostatic distribution), and volumetric strain fields at the end of the simulation period (30 years), respectively; plot (d) depicts the dynamics of injected mass of carbon dioxide in pure hydrodynamical model (red curve) and coupled model (blue curve).

Figure 17 :
Figure 17: Results of hydrodynamical modeling of CO 2 injection into the target aquifer with the upper boundary located at the depth of 1700 m; plots (a) and (b) show the gas saturation and pore pressure increment (compared to the initial hydrostatic distribution) fields at the end of the simulation period (30 years), respectively.

Figure 18 :
Figure 18: Results of coupled hydro-geomechanical modeling of CO 2 injection into the target aquifer with the upper boundary located at a depth of 1700 m with no tectonic stresses (  =   ); plots (a) -(c) show the gas saturation, pore pressure increment (compared to the initial hydrostatic distribution), and volumetric strain fields at the end of the simulation period (30 years), respectively; plot (d) depicts the dynamics of injected mass of carbon dioxide in pure hydrodynamical model (red curve) and coupled model (blue curve).

Figure 19 :
Figure 19: The figure presents the results of coupled hydro-geomechanical modeling of CO 2 injection into the target aquifer with the upper boundary located at a depth of 1700 m accounting for the tectonic stresses (  ∕  ∼ 1.8 in the target aquifer).Panels (a) -(c) show the gas saturation, pore pressure increment (compared to the initial hydrostatic distribution), and volumetric strain fields at the end of the simulation period (30 years).Panel (d) depicts the dependence of the injected mass of carbon dioxide on time.

Figure 20 :
Figure 20: The slice of the sector of the real formation; by solid and dashed black lines we show impermeable and constant pressure boundaries, respectively; the yellow bars mark the perforation intervals through which CO 2 is injected into the upper and lower aquifers.

Figure 21 :
Figure 21: Components of the tectonic stress tensor in the coordinate axes  ′  ′ (   ) and in the principal axes  (Σ   ).

Figure 22 :
Figure 22: Results of coupled hydro-geomechanical modeling of CO 2 injection into the upper and lower aquifers at a constant bottomhole pressure of 350 bar; plots (a) -(c) show the gas saturation, pore pressure increment (compared to the initial hydrostatic distribution), and volumetric strain fields at the end of the simulation period (30 years), respectively.

Figure 23 :
Figure 23: Results of coupled hydro-geomechanical modeling of CO 2 injection into the upper and lower aquifers at a constant bottomhole pressure of 600 bar; plots (a) -(c) show the gas saturation, pore pressure increment (compared to the initial hydrostatic distribution), and volumetric strain fields at the end of the simulation period (30 years), respectively.

Figure 24 :
Figure24: Dynamics of the CO 2 mass injected at a constant bottomhole pressure 600 bar computed via the coupled hydrogeomechanical (red line) and hydrodynamic (blue line) models during the second half of the injection period (15-30 years); results of simulations using modified coupled models (see details in Section 4.3.3)are also shown: with reduced values of the mechanical characteristics in the aquifers and caprock (green line), with strengthened contrast in the mechanical properties between the host rock and the fault core in the fault zone (dashed blue line), with increased permeability in the fault zone (green line).

Figure 25 :
Figure 25: Results of coupled hydro-geomechanical modeling of CO 2 injection into the upper and lower aquifers at a fixed bottomhole pressure of 600 bar and decreased (as compared to base case described in Section 4.3.1)values of the mechanical properties of the storage aquifers and caprock as described in Table4; plots (a) and (b) show the gas saturation and volumetric strain fields at the end of the simulation period (30 years), respectively.

Figure 26 :
Figure 26: The flowchart summarizes the procedure for calculating the host rock , fracture network   , and major crack   permeabilities.

Figure 27 :
Figure 27: The figure shows the representations of the cells belonging to the damage zone of the tectonic fault (panel a) and to its core (panel b).

Figure 28 :
Figure 28: The flowchart summarizes the approach for implementing the permeabilities outlined in Fig. 26 into the hydrodynamical model.

Table 1
Parameters of the vertical well model considered in the numerical experiments.

Table 4
Decreased values of mechanical properties of the aquifers and caprock in the realistic formation model shown in Fig. 20.