Numerical Simulation of Multiphase Hydromechanical Processes Induced by CO 2 Injection into Deep Saline Aquifers

Numerical Simulation of Multiphase Hydromechanical Processes Induced by CO2 Injection into Deep Saline Aquifers — In this paper, the conceptual modeling and the numerical simulation of two-phase flow during CO2 injection into deep saline aquifers is presented. The work focuses on isothermal short-term processes in the vicinity of the injection well. Governing differential equations are based on balance laws for mass and momentum, and completed by constitutive relations for the fluid and solid phases as well as their mutual interactions. Constraint conditions for the partial saturations and the pressure fractions of CO2 and brine are defined. To characterize the stress state in the solid matrix, the effective stress principle is applied. The coupled problem is solved using the inhouse scientific code OpenGeoSys (an open source finite element code) and verified with benchmarks. Deep Saline Aquifers for Geological Storage of CO2 and Energy Stockage géologique du CO2 et de l'énergie en aquifères salins profonds IFP Energies nouvelles International Conference Rencontres Scientifiques d’IFP Energies nouvelles 106 Oil & Gas Science and Technology – Rev. IFP Energies nouvelles, Vol. 66 (2011), No. 1


INTRODUCTION
Among other mitigation concepts Carbon dioxide Capture and Storage (CCS) is worldwide under discussion as an emerging transition technology to reduce anthropogenic greenhouse gas emissions into the atmosphere.In 2005 the Intergovernmental Panel on Climate Change (IPCC) published a special report [1] addressing the state of the art, perspectives and various knowledge gaps related to the longterm storage (sequestration) of CO 2 in the underground.
Three types of geological formations are particularly considered for the safe storage of CO 2 : (nearly) depleted hydrocarbon reservoirs, deep saline aquifers and unminable coal seams.In cases of hydrocarbon reservoirs and deep aquifers, carbon dioxide is injected in a dense form into porous rock formations filling out the pore space and partially displacing in situ residing fluids.According to various studies, deep saline aquifers provide the most substantial carbon dioxide storage capacity [1][2][3][4][5][6][7], and are often located near possible CO 2 sources such as coal-fired power plants.
The migration of the carbon dioxide in the subsurface, and its interaction with the formation fluids as well as with the porous reservoir media is characterized by various complex transport, reaction and deformation phenomena.Several trapping mechanisms, which prevent the migration of the buoyant CO 2 back to the surface, are based on these phenomena.Particularly during injection, advection dominated multiphase flow can be observed, which is mainly driven by buoyant, viscous as well as capillary effects due to different physical properties of the injected carbon dioxide and the formation fluids.These transport processes are restricted by stratigraphic trapping below low-permeable caprocks, structural trapping in folded or fractured rocks, and hydrodynamic trapping of the carbon dioxide as an immobile fraction (residual CO 2 saturation).Dissolution and mixing effects result in the solubility trapping of carbon dioxide in the subsurface, and, finally, geochemical reactions with the reservoir rock (carbonate precipitation) provide the foundation for the mineral trapping of CO 2 .Usually, the time scales associated with different trapping mechanisms are believed to be quite different.Within the context of safety aspects, long-term carbonate formation is the most favorable trapping mechanism, but it can be accompanied by unfavorable changes of the porosity and permeability of the reservoir media if it occurs during injection.
Within the context of mechanical loading, the injection of carbon dioxide into the subsurface results in high pressure in the vicinity of the injection well.Due to the injection pressure, the stress distribution in this reservoir region can be changed significantly.To consider high pressure induced medium deformation is very important for the integrity of potential geological storage sites, which is affected by several geomechanical factors (cf.[8,9] and others).
The modeling and simulation of the injection and the spreading of carbon dioxide in the underground is essential for the proper understanding of the physical and chemical processes at different length and time scales, to ascertain migration and trapping of CO 2 in the porous formations, and in assessing the capacity as well as the safety (possible leakage) of the reservoir.Sophisticated mathematical models, numerical algorithms and software tools have to be developed taking into account all relevant physicochemical phenomena during migration and storage of CO 2 in the subsurface, such as flow and transport of multiple phases including dissolution and mixing effects, hydrodynamic instabilities (e.g.viscous fingering), rock deformation as well as fracturing, heat transport and phase changes.In most cases it is not necessary to be able to describe all these processes for the whole simulation time period or for the entire model domain.The formulation of a model concept has to consider the spatial and temporal scale of the problem and the dominating processes in each case.Then, the model complexity can be adapted to the relevant processes for the given problem while other processes having less contribution to a certain feature or event can be neglected for the sake of computational efficiency.
While transport and deformation processes in porous media have been studied for several decades considering various applications (e.g.groundwater flow and consolidation problem), the modeling of carbon dioxide migration in geological formations is a fairly new subject of investigations in different research areas (e.g.hydrology, geotechnology, computational mechanics and mathematics).As about ten years ago only a few publications have been dedicated to this topic (cf.[10][11][12]), recently an increasing number of numerical (see e.g.[13][14][15][16][17][18][19][20]), semi-analytical (cf.[21][22][23]) and analytical (cf.[7,[24][25][26]) studies have been published.Within the context of numerical simulations, more complex problems can be treated (e.g.coupling of different processes, consideration of heterogeneities and various geological conditions), although the development of efficient and stable algorithms is challenging.Among hydromechanical site specific simulations of CO 2 injection, hypothetical cases are reported (e.g., [8], saline aquifer) as well as real field studies (e.g., [27], depleted hydrocarbon reservoir).
Currently numerical studies of carbon dioxide storage are mostly based on simulators developed for the use in the oil, gas and geothermal energy production.They represent convenient starting points for specialized model and code adaptations targeted at modeling the geological storage of CO 2 (cf.[28][29][30][31][32][33]). For the state of the art in corresponding model and software developments see also [1,12] and literature cited there.
Commercial as well as scientific codes are available for multiphase flow processes, geomechanical deformations including fracture and damage evolution, various nonisothermal effects and chemical reactions.As most of the numerical codes are focused on a subset of these coupled processes, at present, capabilities for a comprehensive treatment of the different phenomena are still restricted.To compare several numerical simulators with respect to their capabilities, efficiency and accuracy code comparison studies have been conducted [34,35].Within this context, benchmark problems and real site studies have been defined addressing various aspects of CO 2 storage at different reservoir conditions.
Modeling of CO 2 storage on a reservoir scale for feasibility studies and risk analyses is very demanding with respect to the computational costs due to the complex geometry that needs to be described and due to the diversity of interacting hydraulic, thermal, mechanical, and geochemical processes.According to the evolving processes at various time and length scales we propose a successive model development with the final goal of a comprehensive coupled simulation of all relevant physico-chemical effects, but starting with the most relevant effects at CO 2 injection near the injection well.Since the model complexity is adapted to these dominating processes, other processes having less contribution to the event under consideration are neglected for computational efficiency (e.g.dissolution and mineral trapping).Within this context, in this paper, the investigation is restricted to hydromechanical effects of the CO 2 injection into deep saline aquifers at isothermal conditions.
In the present study, we utilize numerical methods to analyze the stress changes caused by the interaction with the two-phase fluids during the injection.To this purpose, we first set a conceptual injection model, which represents the two-phase flow process of CO 2 and brine, and also the deformation process in the near field in deep saline aquifers.To enable a stable and efficient simulation, the coupling terms in the two-phase flow and the deformation processes are handled by a mixture approach of monolithic and staggered schemes.Moreover a pressure-pressure scheme is applied to the two-phase flow model for convenience in stress calculation in porous media [36].With the present assumption, the changes in flow and deformation fields of the model are simulated by using the standard Galerkin finite element method, which is realized in an object oriented scientific tool, Open-GeoSys (an open source finite element code) developed by the authors.Finally, we present some simulation results and discussion.
In the following, vectors and higher order tensors will be denoted by boldface characters in direct notation.Their scalar product is characterized by a single dot while double dots indicate the summation product (double inner product) a • • b = a i j b ji .A superposed dot indicates the material rate of a vector or a tensor, a superscript T the transposed tensor.

Preliminary Remarks
It is widely accepted to consider Biot's studies on the macroscopic theory of saturated wet soils (cf.[37]) as crucial for the development of models and algorithms for the numerical simulation of flow processes in porous media.Based on these studies, the Theory of Mixtures as one of the basic approaches to model the complex behavior of porous media has been developed (concerning basic assumptions see e.g.[38,39]).As the Theory of Mixtures does not incorporate any information about the microscopic structure of the material (1) , it has been combined with the Concept of Volume Fractions by e.g.[40][41][42][43].Within the context of this enhanced Theory of Mixtures (also known as Theory of Porous Media), all kinematical and physical quantities can be considered at the macroscale as local statistical averages of their values at the underlying microscale.
Concerning a detailed overview about the history of the modeling of the behavior of multiphase multicomponent porous media, the reader is referred to e.g.[44].Comprehensive studies about the theoretical foundation and numerical algorithms for the simulation of coupled problems of multiphase continua are given in e.g.[42,44,45] and the quotations therein.Recent developments of enhanced numerical approaches for the porous media modeling are reported particularly in the classical applications of geomechanics.They are frequently focused on coupled multiphysics issues for partially saturated media (cf.[46] and others).In e.g.[47,48] modified Galerkin type variational formulations as well as stabilized approaches are presented for the simulation of consolidation problems.Numerical strategies for the treatment of thermo-hydromechanical processes are introduced in e.g.[49][50][51][52].
Following the assumptions about porous media, the subsurface formation designated as carbon dioxide reservoir is considered as a mixture of a solid skeleton and a pore fluid content, which can either be: -a single liquid or gas (single-phase flow in porous media), -an immiscible fluid mixture of gas and liquids (multiphase flow in porous media) or, -a miscible fluid mixture of different reacting constituents allowing phase transitions due to evaporation, condensation, precipitation (multiphase multicomponent flow in porous media).The individual constituents ϕ α of a porous material represent the phases of the overall aggregate or components within a phase.Below, α = s marks one immiscible solid phase (no sorption processes are considered), and α = γ denote several immiscible pore fluid phases.
Within the framework of the Concept of Volume Fractions, scalar variables like volume fractions and saturations are defined to describe the microstructure of a porous medium in a macroscopic manner neglecting the real topology and distribution of the pores.These variables serve as (1) Within the context of the Theory of Mixtures the ideal mixture of all constituents of a multiphase medium is postulated.Consequently, the realistic modeling of the mutual interactions of the constituents is difficult.
measures of local fractions of the individual constituents.
The volume fractions n α represent the ratio of the partial volume dv α of a given constituent ϕ α of a multiphase body with respect to the overall volume dv of a Representative Elementary Volume (REV) of the control domain Ω under consideration.Consequently, based on the definitions of the overall volume of the control domain and the corresponding partial volumes of the individual constituents the volume fractions provide some information about the local volume distribution of the individual constituents.
One of the most characteristic media properties of a porous material is the porosity, the local amount of fluid volume fractions.
Since, in general, the overall medium is completely filled with matter, from Equation (2) follows the saturation condition regarding the overall aggregate.
If multiphase flow occurs, it is more convenient for various applications to use the (partial) fluid saturations S γ instead of the volume fractions.These local functions are given by obviously fulfilling the saturation condition regarding the pore content.
Usually, constraint conditions addressing real physical effects are formulated to simplify complex mathematical and numerical models.Within the context of porous media, it is reasonable in most applications to assume the (material) incompressibility of constituents as a substantial constraint condition.The issue of (in)compressibility of a material is closely connected to the possible temporal evolution of its mass density.
Within the framework of the Concept of Volume Fractions, two different formulations of mass density related to the constituents of a porous medium are introduced.The so-called material (effective, realistic) density ρ αR is defined as the ratio of the mass fraction dm α of the given individual constituent ϕ α with respect to its partial volume fraction.
In contrast, the so-called partial (global, bulk) density is given by the ratio of the mass fraction of the constituent under consideration with respect to the volume fraction of the overall aggregate.
Based on the definition of the volume fractions Equation ( 3), the material and the partial densities are correlated to each other.
If the volume fractions change with time under external loading, from Equation (11) follows that for an intrinsically incompressible individual constituent (constant material mass density) compressibility referred to the overall aggregate is observed.
Obviously, the mass density of the porous medium (homogenized overall aggregate) is defined as the sum of the partial densities of its constituents.
The conceptual idea behind the formulations and relations presented above consists in the assumption that the mass fractions of all constituents of the multiphase medium are simultaneously present and statistically uniformly distributed over the entire control domain.Within this context, the material body under consideration is theoretically substituted by an aggregate completely and continuously filled by superimposed (overlapping) homogenized partial continua.In other words, all constituents of a porous medium are characterized as smeared substitute continua with reduced mass densities.Consequently, the motion and physics of the individual constituents as well as the overall aggregate can be specified by well-accepted phenomenological methods of continuum mechanics.
Describing the transport and deformation of the constituents of porous media within the framework of continuum mechanics it is assumed that the geometry of the control domain under consideration is characterized at each time by the solid skeleton, whereas the fluid pore content is able to flow across the boundary of the surface.This assumption serves as conceptual nucleus for the simulation of complex, coupled physical processes in multiphase porous media, particularly if a deformable solid skeleton is observed.Within this context, it proves to be reasonable not to model the absolute motion state of the pore content, but its motion relative to the motion of the solid phase, considering the porous medium as a local thermodynamic open system with the solid skeleton as volume under observation.
The macroscopic characterization of the physical processes considering the real microstructural situation in a statistically averaged manner is completely adequate for most engineering, geotechnological and biomechanical problems under consideration (cf.[53] and others).

Basic Assumptions
In the following, isothermal transport and deformation processes during carbon dioxide injection into deep saline aquifers will be studied using two-phase flow models in deformable porous media, which are based on the concepts of the Theory of Porous Media.Within this context, brine, as native pore fluid content and CO 2 represent the fluid constituents.At this stage, for the solid skeleton it is assumed that no fracturing and damage occurs.Consequently, it is justified to postulate that small deformations can be observed, whereas the stress response is analyzed based on linear elastic material models.
It is generally accepted that, in order to maximize the use of the pore space, CO 2 should be injected into saline aquifers as supercritical fluid, because under these conditions it behaves like a gas but having the density of a liquid (however, being less dense and less viscous than brine).Within this context, carbon dioxide and brine are assumed to be immiscible.Although controversially discussed, most of the authors assume that the dissolution of CO 2 into brine is a very slow process (cf.[1,2,5,6,54] and others).Therefore, as we are at this stage interested in short-term migration processes of CO 2 in the vicinity of the injection well, dissolution and other phase transition effects are neglected.
To displace the brine to its residual saturation, CO 2 is injected at high pressure.After injection CO 2 will form a plume around the injection well and tends to migrate upward due to the density difference of the fluids involved in the two-phase flow processes.This upward flow continues until the CO 2 plume reaches the caprock assumed to be impermeable.
The high injection pressure of carbon dioxide could result in unfavorable deformation of the porous rock formation possibly causing fractures or the opening of failures along fault planes.Consequently, the geomechanical modeling of deformation processes in CO 2 is necessary for storage site assessment focusing on maximum acceptable formation pressure.

Governing Equations
The governing field equations for the modeling of transport and deformation processes of multiphase flow in deformable porous media are formulated based on the specific (local) individual balance relations of the constituents, particularly mass and momentum balances.Within this context, the application of the fundamental balance relations on the analysis of multiphase materials is based on Truesdell's metaphysical principles (2) .The crucial idea behind these principles is the assumption that the balance relations of the constituents as well as the balance relations of the overall aggregate of a porous medium can be formulated in accordance to the corresponding classical relations of singlephase continuum mechanics.Additionally, to account for the interaction mechanisms between the constituents, socalled production terms are introduced for the individual balance relations of the constituents (cf.[45]).Following, superposition is used to define the balance relations of the overall aggregate based on the individual balance relations of the constituents.
Neglecting mass exchange between the phases (no dissolution and sorption processes), the local mass balance for the individual constituent ϕ α of the porous medium is given by with the velocity u α of the constituent under consideration, and the usual divergence operator ∇ • ().From the velocitydisplacement relation for the solid skeleton follows with the solid displacement vector u s .The derivative with the usual gradient operator ∇() denotes the material time derivative of an arbitrary variable a with respect to the motion of a material point of the constituent ϕ α .It consists of a local (diffusive) part and a convective part associated with the velocity of the constituent.As mentioned above, the transport processes of the fluid constituents of a porous medium are considered as their relative motion with respect to the motion of the deformable (2) Truesdell [55] formulated the following principles within the context of the Theory of Mixtures: 1.All properties of the mixture must be mathematical consequences of properties of the constituents.2. So as to describe the motion of a constituent, we may in imagination isolate it from the rest of the mixture, provided we allow properly for the actions of the other constituents upon it.3. The motion of the mixture is governed by the same equations as is a single body.These statements were later being called metaphysical principles.By now, this is a term, which is frequently used by many authors in literature about the Theory of Porous Media.
solid skeleton.Consequently, the relations between the material time derivatives (here, of an arbitrary variable a) with respect to the solid skeleton, and with respect to the individual fluid constituent ϕ γ is of crucial interest in terms of a unified numerical characterization of the different processes.
Here, u γs = u γ − us is the so-called seepage velocity describing the fluid motions with respect to the deforming skeleton material.
According to the generalized formulation Equation ( 14), considering Equations ( 5) and (11), the local solid phase mass balance is given by Following the same procedure, additionally considering Equations ( 7) and ( 17), the mass balance relations for the fluid constituents ϕ γ with can be defined with respect to the solid phase motion.
Assuming material incompressibility of the solid phase, i.e. d s ρ sR /dt = 0, and applying the solid phase mass balance Equation (18), Equation ( 20) can be represented in a more detailed description. Here is usually known as filter velocity of the motion of the pore fluid constituent ϕ γ .In geotechnical problems, the internal fluid friction forces can be neglected in comparison to the interaction terms between fluid and skeleton motions.Thus, the total Cauchy stress tensor σ referring to the local loading state of the overall aggregate is given by the sum of all partial stresses of the constituents.Consequently, the stress tensor is defined according to the well-known effective stress concept [56]: with the solid effective stress tensor σ s E and the identity tensor I .Therein, the pore pressure p is given in analogy to Dalton's law by with the fluid pressure fractions p γ .For the sign convention, it should be noted that, according to the usual convention in continuum mechanics, compressive coefficients of the stress tensors σ s E and σ are assumed to be negative.Thus, the overall pore pressure as well as the pressure fractions related to the fluid constituents is assumed to be positive.
It can be shown that the decomposition of the stress tensor Equation ( 23) results thermodynamically consistently (satisfying the balance relations of thermodynamics) from the assumptions of the Theory of Porous Media (see [41,44,45]).As mentioned above, among others, the material (intrinsic) incompressibility of the solid phase is assumed, which is justified for the case that the relative matrix compressibility can be neglected compared to the compressibility of the overall aggregate due to changes of the pore space content of the porous medium.In the literature, different definitions of the overall stress tensor characterizing the porous medium mechanics are presented.Particularly in rock mechanics, Biot's poroelastic approach with the effective stress component σ and the scaling factors α γ is frequently used (originally, Biot proposed an according stress decomposition in case of saturated single-phase flow in porous media with the well-known single Biot's coefficient α).Different effective stress concepts within the poroelastic framework of saturated and unsaturated soils are discussed in [57].As we used the Theory of Porous Media as theoretical foundation of the numerical analyses presented here, in the following, the representation Equation ( 23) of the overall stress tensor is preferred.
Completing the set of necessary balance relations according to the problem under consideration, the stress field in the saline aquifers is governed by the overall local momentum balance relation where ρg is the volume force with the gravity vector g.

Constitutive Relations
If the balance relations characterize fundamental physical and thermodynamical properties of the matter independently of specific material properties, in real applications the response of a physical body on similar interactions with the external environment differs for various materials.Thus, socalled constitutive relations have to be defined to characterize the specific material behavior.In terms of the mathematical modeling of physico-chemical processes this observation is equivalent to the formulation of closed systems of equations, which should consist of balance as well as constitutive relations.
Within the context of the multiphase problem under consideration, constitutive equations are required for selected production terms of the specific balance relations of the individual constituents, as mentioned above, for pore fluid properties like pressure and saturation, and for the partial effective stress tensor of the solid skeleton.
Interacting during their motion (transport and deformation), the fluid constituents exchange linear momentum fractions, which are considered defining appropriate production (coupling) terms for the equations of momentum balance for the fluid constituents ϕ γ .These terms are introduced to characterize effects on the microscale within the context of macroscopic, phenomenological approaches.It can be shown that, based on appropriate constitutive relations for the momentum production terms, the fluid momentum balances are represented by modified Darcy's law for multiphase flow as a constitutive flux expression (cf.[45] and quotations cited therein), with k γ rel denoting the relative permeability, μ γ the dynamic viscosity, and k the intrinsic permeability of the solid skeleton.
We assume all constituents to behave materially incompressible and define the capillary pressure in the usual way as the difference between the nonwetting and wetting fluid pressure fractions.The capillary pressure-saturation functions as well as the relations between relative permeability and saturation are substantial constitutive equations required for multiphase flow.Within this context, usually algebraic expressions are fitted to the corresponding experimentally observed curves.Among the widely-used of these algebraic expressions are the Brooks-Corey [58] and van Genuchten [59] relations.If both are realized within the scientific software code developed by the authors, the numerical results presented in this paper are based on Brooks-Corey's approach.
The Brooks-Corey equations relating the saturation to the capillary pressure are where p D is usually known as entry pressure, λ is a poresize distribution index, and S eff is a normalized wetting fluid saturation defined following [60]: where S l res is the wetting phase residual or irreducible saturation, and S CO 2 res is the nonwetting phase residual saturation.The constitutive parameters p D , λ, S l res and S CO 2 res are identified by fitting Equation (30) to experimental data.Within this context, the entry pressure is to be understood as the minimum pressure that the nonwetting fluid must have to enter the largest pores.The relations between the relative permeability and the saturation are given by Additionally to the presented constitutive relations characterizing media properties, appropriate equations of state for pressure dependent fluid properties (density, viscosity) are considered.
As mentioned above, the host rock of the reservoirs under consideration is assumed to be linear elastic.Following the geometrically linear approach (small deformations) to elasticity, the overall strain tensor ε s is defined as: In the isotropic case, the solid effective stress is governed by the generalized Hookean law where μ s and λ s are the Lamé constants of the porous material.The Lamé constants can be expressed by the Young's modulus E, the Poisson's ratio ν and the shear modulus G (so-called engineering constants) which can be obtained experimentally.
Typically, real geological formations are characterized by heterogeneous, frequently layered soil and rock structures, which requires the application of anisotropic constitutive models for skeleton deformation.Details of wellknown anisotropic elastic models, like transverse isotropy or orthotropy are given in a huge number of textbooks, monographs and scientific papers.

NUMERICAL SCHEME
The numerical treatment of the coupled problem of twophase flow in deformable porous media is based on the governing field equations together with discretization methods in the space and time domains.Some general representations of corresponding numerical approaches to solve the problem under consideration can be found in [42,47,61,62].The method of weighted residuals is applied to derive the weak formulations of all the governing equations given above.Within the framework of a standard Galerkin procedure, the local overall momentum balance equation Equation ( 26) is multiplied by an arbitrary test function ū = ū(x) (with ū ∈ (H 1 0 (Ω)) 3 , ū = 0 on the Dirichlet boundary Γ u of the overall medium), and integrated over the current domain Ω bounded by the solid skeleton.
Therein t is the external load vector acting on the Neumann boundary Γ t of the overall medium.If the test function is interpreted as a virtual displacement vector function, Equation (38) represents the principle of virtual work.Following the same procedure, the mass balance relations of the pore fluids Equation ( 21) are multiplied by arbitrary test functions qγ ∈ L 2 (Ω) with qγ = 0 on the Dirichlet boundary Γ p of the overall medium, and integrated over the current domain bounded by the solid skeleton φγ qγ dΓ q (39) with φγ characterizing the efflux of fluid mass of the constituent ϕ γ through the Neumann boundary Γ φ .
Concerning the coupled problem under consideration, the filter velocities can be eliminated by use of the Darcy's law Equation (27).Thus, w γs loses the status of an independent field variable, and pressure variables are considered instead.Additionally, considering the pressure-saturation condition and constitutive relations between the material mass densities and pressure fractions, solid skeleton displacements and fluid pressure fractions remain as independent variables in Equation (39).
For two-phase flow problems, the choice of primary variables is of crucial importance for the stability of the employed numerical scheme [36].Within the context of the example presented below, capillary pressure and CO 2 pressure (pressure-pressure scheme) have been selected as primary variables for the two-phase flow modeling.As an alternative, a pressure-saturation scheme has been realized using the pressure of the wetting phase (here: brine) and the saturation of the nonwetting phase (here: CO 2 ).Naturally, the solid phase displacements serve as the primary variable modeling the deformation processes.
As usual, in the finite element space the continuous functions of the selected primary variables are interpolated based on their nodal values and appropriately defined shape functions.After discretizing the weak forms of the balance relations a system of nonlinear algebraic equations can be obtained.
with the mass matrix M and the Laplacian K .Equation ( 40) has three sets of unknowns, i.e. p c , p CO 2 , u at element nodes.
For the required time discretization we use the generalized first order difference scheme with the time step Δ t = t n+1 − t n (42) and the weighting factor developed for the solution of ordinary differential equations For a compact representation, the finite element system Equation ( 40) is given symbolically by Considering Equation (45) like a system of differential equations in terms of Equation (44), the discretization scheme Equation ( 41) can be applied.After some straightforward algebra we get In most applications, a fully implicit Euler scheme is used for the time discretization (i.e., α = 1).The implicit Euler scheme is unconditionally stable and can be applied to complex problems such as of multiphase flow in porous media.However, numerical diffusion is introduced into the system by time discretization.This problem is reduced by an adaptive time-stepping scheme.
The nonlinear coupled boundary value problem is solved iteratively using the Picard linearization within the context of the finite element method [52].Within this context, usually all unknowns can be solved at the same time with a so called monolithic scheme [42].However, solving the whole system of equations monolithically may lead to memory problems when employing a fine mesh due to four or more degrees of freedom per element node of the coupled problem under consideration.We adopt a mixture of monolithic and staggered schemes to avoid this bottleneck without losing the accuracy of the solutions.The coupling procedure is as follows: 1. Solve ( p c , p CO 2 ) monolithically.2. Solve the equation of u. 3. Go to step 1. if the converged solution of p c , p CO 2 , u is not achieved.
At each linearization step, the systems of linear algebraic equations are solved using an iterative BiCGSTAB solver in conjunction with ILU preconditioning.
The presented numerical scheme is realized in an object oriented scientific software tool, OpenGeoSys, developed by the authors, and has been verified by several classic benchmarks such as Buckley-Leverett [63], Liakopoulos [64], McWhorter-Sunada [65], Kueper et al. [62,66] (for the OpenGeoSys applications cf.[67]), and the DECOVALEX project [68].dioxide is in supercritical state.As mentioned above, this is the most advantageous state for CO 2 to be injected into deep saline aquifers.Depending on the real geological conditions (i.e., pressure and geothermal gradient), these conditions are found at depths starting from about 600 . . . 1 000 m [69].Considering that carbon dioxide is injected into a deep confined aquifer in supercritical state, for standard operating conditions it remains supercritical and in thermal equilibrium with the residing pore fluid (cf.[7,69] and quotations cited therein).

NUMERICAL EXAMPLE: MODEL RESERVOIR
Our example represents a hypothetical (synthetic) problem of the CO 2 injection into the subsurface based on selected data (geometry, reservoir conditions) of a realistic formation.We assume that the deep saline aquifer is located at a depth of 770 m from the ground surface and has a thickness of 6 m as illustrated in Figure 1.It is assumed that the saline aquifer is fully saturated with brine before injection.The reservoir temperature is 34 • C.
In the present study, the near field of the injection well is the domain of interest.To this purpose, we assume that the problem is axisymmetrical in both geometry and physics.Taking the injection well's radius as 0.2 m and cutting the 3D domain at the radius of 200 m, we generate a finite mesh for axisymmetrical analysis as depicted in Figure 2. as special case of Equation ( 21), and the stress decomposition Equation ( 23).The Brooks-Corey's model is employed to characterize the hydraulic properties of supercritical CO 2 and brine in the porous medium.Material parameters for the two-phase flow and the deformation processes are taken from literature [34] and the equations of state for CO 2 and water from [67].Within this context, for the hypothetical situation under consideration (one homogeneous reservoir layer, short distance near field model) isotropic linear elastic material behavior is assumed, and the constitutive relation Equation ( 35) is applied for the characterization of deformation processes.All values of the material parameters are given in Tables 1 and 2. Initially, the stress tensor coefficients in the deep saline aquifer are assumed to be caused by the gravity force only, the distribution of which is calculated by solving Equation (26) with the volume force term ρg (lithostatic conditions).This initial distribution analysis is used as the initial stress condition for modeling the CO 2 injection.
For the carbon dioxide influx into the reservoir, a constant Neumann boundary condition of 0.4475 • 10 −5 m/s (0.03374 L/s accordingly) is applied at the injection well boundaries for the hydraulic equations Equation ( 21) to represent the injection rate of CO 2 .If the injection borehole does not have any water flow, the maximum water saturation and residual CO 2 saturation are assigned in the terms of Dirichlet boundary conditions at the end of the domain in the radial direction.On both of the top and bottom surfaces, the displacement in vertical direction is fixed.In the horizontal direction, displacement is fixed only on the outer boundary.The inner boundary is assumed to be free to move.
In order to investigate the impact of solid deformation on the hydraulic field, the simulations of the two-phase flow process in the deep saline aquifer were performed considering and neglecting the coupled deformation process, respectively.To demonstrate this impact, in Figure 3 the change of CO 2 saturation during a period of 1 000 hours at two specified observation points is plotted, which are at a distance of 20 m and 50 m from the injection well, respectively.
Figure 3 portrays that the propagation of CO 2 is slightly enhanced by deformation.On the other hand, the stress field is significantly altered by the CO 2 injection pressure in the near field.We can clearly see stress changes at the two observation points as depicted in Figure 4.The tangential stress decreases at the beginning of the injection due to the extension at the well surface, and then increases due to Spatial distribution of the CO 2 saturation 1 and 1 000 hours respectively after injection started.
the propagation of the injection pressure.Since we assume that the initial stress is induced by the gravitational force exclusively, the initial stress distribution in the analyzed domain has a vertical gradient.Due to the injection pressure of CO 2 , the stress in the tangential direction increases significantly.Figure 5 shows the distribution of the tangential stress near the injection well at the time of one hour and one thousand hours after the beginning of the injection.
Corresponding to the stress field demonstrated in Figure 5, the distribution of carbon dioxide in the vicinity of the injection well at the same time is provided in Figure 6.Comparing Figure 5 and Figure 6, we see that the distribution Spatial distribution of tangential stresses 1 and 1 000 hours respectively after injection started.
patterns of tangential stress and CO 2 are quite similar.This implies the coupling effect between hydraulic and mechanical processes.

CONCLUSION
In this paper, we presented the conceptual model and the numerical algorithm for the simulation of isothermal twophase flow in deformable porous media.The discussed numerical scheme has been implemented in the open source scientific finite element code OpenGeoSys.The study is part of an ongoing development of the theoretical and numerical framework for the solution of Thermo-Hydro-Mechanical-Chemical (THMC) coupled problems related to CO 2 storage in geological formations.The capabilities of the method and the software tools were shown on the example of the simulation of short-term near field migration processes during injection of supercritical carbon dioxide into a deep saline aquifer.With this fairly simple axisymmetrical example of the coupling of two-phase flow and deformation processes, it was demonstrated that the stress change in the vicinity of the CO 2 injection well is distinct from the initial stress state.Consequently, such coupling must be considered as an important issue for the safety assessment of carbon dioxide injection and storage.In the presented test example, a constant CO 2 pressure prescribed against the injection well is assumed from the beginning of the injection.This leads to an instant increase of stress at the well, which may intrigue mechanical damage of the surrounding rock such as fracture development.
On the other hand, the stress distribution in the near field approaches to a kind of 'homogenous' state, which makes a safer environment for the near field aquifer.The present stress results prompt that a gradually increasing injection rate can be a safe process of injection.More details of near field stress state taking into account material nonlinearities and fractures in porous media would be future research topics.The numerical scheme presented above proved to be stable, and the numerical results showed a good agreement with given benchmark results.
To enlarge the range of applications regarding time scales and reservoir characteristics, hydraulic and mechanical anisotropy, fracture and damage analyses, non-isothermal processes, diffusion (mixing of fluid components), supercritical fluid behavior and phase changes due to various mechanisms will be considered in future work.These aspects are particularly important for deep reservoirs (e.g.depleted gas fields) where the thermodynamical state variables vary significantly during CO 2 injection and migration.The presented generalized concept is also suitable for other geotechnical applications such as geothermal reservoir analyses and safety assessment of nuclear waste repositories.

1
Figure 1 Axisymmetric model of CO 2 injection into deep saline aquifer analyzing near field liquid spreading.

Figure 2 Triangular
Figure 2Triangular mesh for finite element analyses Figure 3Temporal evolution of the CO 2 saturation in two selected points.

45 Figure 5
Figure 4Temporal evolution of tangential stresses in two selected points.