The OGS-Eclipse code for simulation of coupled multiphase flow and geomechanical processes in the subsurface

This paper presents a numerical simulation tool for the analysis of coupled processes related to subsurface operations. The tool combines the open-source scientific code OpenGeoSys with the reservoir simulator Eclipse enabling the coupling of thermal, hydraulic, mechanical and geochemical processes. While the coupling of multiphase flow with heat and reactive geochemical component transport has been already implemented, OpenGeoSys-Eclipse is now extended for the coupling of multiphase flow and deformation. By this, OpenGeoSys-Eclipse is capable of addressing the impact of pore pressure changes on rock stability and deformation as well as the feedback effects of geomechanical processes on multiphase flow via pore volume coupling and porosity and permeability update. The coupling is verified by several test cases of gas storage scenarios and compared with reference simulations of OpenGeoSys. The results are in good agreement regarding the general effects of geomechanical feedback on pore pressure as well as porosity and permeability changes. Differences in the results are only observed for the pore volume coupling arising from the different implementation of rock compressibility models in the two simulators. The simulations are furthermore used to investigate the relevance of addressing geomechanical feedback in numerical scenario simulations for the assessment of subsurface operations. The results show clearly, that, depending on the given storage site conditions and rock types, the feedback of deformation on pore pressure can be significant and should therefore be accounted for in the assessment.


Introduction
Long-time monitoring of the earth's climate has shown an increasing impact of the anthropogenic CO 2 emissions on global warming [13]. CO 2 is one of the main greenhouse gases enhancing atmospheric radiation and warming of the ground surface. Due to the large industrialization and population growth, CO 2 emission has increased significantly in the last decades enhancing global temperature rise and therefore climate change [25]. This has emerged public and politic concern to climate change worldwide. Therefore, international climate protection agreements has been adopted to limit global temperature increases by reducing industrial carbon emissions and passing to a sustainable energy supply from renewables [24].
One method to reduce CO 2 release from industries to the atmosphere is the permanent storage or sequestration of CO 2 in the subsurface. CO 2 is hereby captured from fossil-fuel power plants as one of the main emission sources for CO 2 and stored in deep geological porous formations. Impermeable sealing cap rock formation above the storage formation inhibit the CO 2 from migrating upward into the shallow subsurface. Potential global storage capacity is estimated to be of about 8000 to 15,000 Gt CO 2 [23]. With actual 31 Gt CO 2 global emissions per year, the amount of CO 2 released to the atmosphere could be significantly decreased by this method [25]. Further efforts in reducing CO 2 emissions are done in the field of energy production from renewable sources, which has been significantly enhanced in the last decades [14]. One disadvantage of these energy sources is their natural fluctuation, e.g., wind and solar power, resulting in a fluctuating electricity delivery to the grid [53]. To overcome this gap in renewable energy production, new storage technologies have been considered for the temporary underground storage of surplus energy [22]. These include compressed air energy storage and power-to-X technologies, where electricity from renewables is used to produce hydrogen by hydrolysis, which can be stored underground or in a second step, can be converted to methane and then stored [34].
Besides new technologies for subsurface gas and energy storage, also existing subsurface application like oil and gas production, groundwater withdrawal and retention, nuclear waste disposal or the reinjection of waste and thermal water has increased due to the rising demand of energy and water consumption worldwide [71]. This increase in new and old subsurface engineering and geotechnical applications have evoked growing awareness of the environmental and drinking water protection, as any subsurface activity is accompanied by environmental impacts on the subsurface including pressurization and stresses, heating of the ground, migration of fluids or changes in the geochemical and biological system [2]. To assess these potential thermal, hydraulic, geochemical and geomechanical processes and impacts related to subsurface applications, a profound understanding of the complex subsurface system is essential.
Especially the storage of mass or energy in the subsurface can induce large-scale hydraulic and geomechanical impacts on the subsurface environment [72]. The injection hereby causes an overpressure at the injection well, which propagates far into the surrounding formations, both laterally and vertically as shown, e.g., for the case of long-term CO 2 storage [7,47,54]. Similar effects were also investigated for the topics of natural gas and hydrogen storage [53,64] or compressed air energy storage [67]. The pressure increase also affects the earth's stress field. Depending on the site-specific geomechanical conditions and the injection-induced overpressure, changes in the subsurface stress field may lead to large-scale deformation processes. This can result in ground uplift, hydraulic fracturing or fault reactivation within the cap rock risking the long-term integrity of the storage site [36,38,50]. Therefore, the injection rate needs to be limit to a maximum overpressure allowed to be applied in the target storage formation to avoid any rock damage [50]. The injection-induced deformation processes in turn affect injection performance as they cause an increase or decrease of the pore volume of the porous rock leading to changes in the fluid pressure within the pores as well as in the relevant fluid flow parameters permeability and porosity [21]. This strong interaction between fluid flow and deformation reinforces the importance of a thorough process assessment and determination of geomechanical effects during subsurface storage operations.
Coupled numerical process simulations help to consider the link between fluid flow and deformation. They allow to predict the individual processes and their feedback to the system and with that to adequately assess the subsurface operation. Two main coupling components have to be considered in the process coupling: volume coupling, accounting for pore volume changes due to stress changes, and flow properties coupling, considering strain and stressrelated changes in permeability and porosity [56]. Several numerical methods exist capable of solving the complex differential equation system that describes the process coupling. The fully coupled method solves for all processes simultaneously through one system of equations. It offers a high internal consistency and accurate solution, but requires long computation times and large development efforts [57]. The iterative coupling method is used when each process is solved separately. Information between the flow simulator and the geomechanics module is transferred back and forth until convergence of fluid flow and stress unknowns is reached. The iterative coupling is less accurate and need some more iterations, but results in a faster run time [63]. It is therefore a good compromise between accuracy and feasibility and more flexible than fully coupled as it enables the coupling of various reservoir simulators and geomechanics codes [45,62]. The explicit coupling method is the loosest coupling method. Results of the fluid flow are sent to the geomechanical process to calculate deformation, but there is no feedback to the fluid flow.
Various simulation codes have been developed for the numerical modelling of coupled processes. These include STARS [8], FEMH [10], CODE-BRIGHT [40,65], DuMuX [16], DYNAFLOW [46], or OpenGeoSys [28,30], to name a few. These codes can be used for all types of coupling, especially fully coupled simulations, as all processes are solved within one code. Though these codes provide high numerical accuracy, the development effort is high. Another more efficient option is the coupling of two existing codes, where each code solves for one process and the codes are then coupled using the iterative or explicit coupling. This has been done, e.g., for TOUGH-FLAC [49,52], TOUGH2-Code-Aster [48], Eclipse-Visage [41,42], NUFT-GEODYN-L [38], or Sierra [35]. Most of these codes were originally flow simulators, which has been coupled to e.g. a mechanical module. This is of huge benefit instead of developing a complete new code as the coupling is more flexible and can be adjusted more easily to the individual requirements.
However, many codes, whether single or coupled codes, are subject to individual limits regarding process spectrum, available coupling methods or simulation efficiency. Therefore, a new coupling approach is presented here, which combines two well-developed and well-proven simulation codes in order to cover the entire complexity of subsurface processes by one simulation code and to benefit furthermore from the individual strengths of both codes, e.g. robustness, flexibility, run time speed and numerical stability. For this, the scientific open-source code Open-GeoSys is coupled to the commercial Eclipse software suite [55]. OpenGeoSys provides the simulation of coupled thermal-hydraulic-mechanical-chemical processes in the subsurface with flexible process coupling and has been applied to a variety of scientific problems, e.g. energy storage [9,33], reactive transport [58,59], geochemical modelling [12] or heat storage [39]. Eclipse is well known for its stable and robust simulation capability for multiphase and multicomponent flow [11,26]. By coupling OpenGeoSys with Eclipse, the process capabilities and flexibility of OpenGeoSys can be combined with the efficiency, the high numerical stability and the computation speed of Eclipse [19,44]. Up to now, the OpenGeoSys-Eclipse simulator can handle single and multiphase flow, component transport, heat transport and geochemical reactions [44] and has been successfully applied to coupled THCproblems, e.g. CO 2 storage or cyclic hydrogen injection and extraction [19,33,37,43]. Coupled hydromechanical processes, especially the geomechanical feedback on flow are not yet considered.
Therefore, OpenGeoSys-Eclipse is now extended for coupled flow and geomechanical processes including geomechanical feedback. Fluid flow is hereby solved in Eclipse and deformation is solved in OpenGeoSys. The new implementation includes the whole complexity of the process feedback between the two processes including pore volume update and changes in the fluid flow parameters porosity and permeability and can be run in an iterative or explicit manner to enable the individually optimal method for the given problem.

OpenGeoSys
The open-source scientific software OpenGeoSys is a finite element code [27,28,68]. It uses an objectoriented and process-oriented approach that allows the solution of partial differential equations for different physical problems in the subsurface using a generic object structure [27]. OpenGeoSys has been already used for multiple underground storage applications, e.g. gas storage or nuclear waste disposal [4,9,18,29] and could be successfully compared and verified to other simulation codes [3,11,17,21,69]. Coupled processes can be solved sequentially (iterative and explicit coupling) or monolithic (fully coupled). Fluid flow can be solved in a pressure-pressure or in a pressure-saturation formulation, using pressure of the non-wetting phase p nw and the capillary pressure p c or pressure of the wetting phase p w and saturation of the non-wetting phase S nw as primary variables, respectively. For the simulation of geomechanical process, several constitutive models are implemented, e.g. a poroelastic model using Biot formulation solving for solid displacement, stresses and strains. OpenGeoSys therefore provides a wide range of geomechanical process applications.

Eclipse
Eclipse is a conventional reservoir simulator suite mainly used by the oil and gas industry. It offers high robustness as well as computational speed. The simulator suite involves two flow simulators E100 and E300 [55] for the simulation of black oil and compositional system, respectively. Both Eclipse simulators use the finite volume method. The implemented numerical methods provide a fully implicit formulation or an implicit pressure, explicit saturation (IMPES) formulation. Besides traditional oil and gas reservoir simulations, Eclipse has been also applied to further energy-related problems as unconventional resources and underground storage of CO 2 , H 2 , and compressed air energy [26,43,54,67]. The advantages of Eclipse regarding simulation speed and numerical stability have been demonstrated by several authors [1,11]. To predict geomechanical processes related to subsurface operations, Eclipse can be coupled to the finite element code Visage. Visage has been a stand-alone simulation code [41], which has been incorporated into the Schlumberger Software Package as Reservoir Geomechanics Tool in 2010. More details and application examples of this Eclipse-Visage code will be presented in Section 5. In this study, only Eclipse is used for the simulation of multiphase flow within the coupled OpenGeoSys-Eclipse code.

Mathematical model
The coupling of fluid flow and deformation is described by a set of basic equations solved for the variables fluid pressure p and solid displacement u. The equation formulations are based on the continuum approach representing the porous medium on an averaged, macroscopic scale with continuously distributed constituents [20].

Geomechanical process
The solid phase mechanics can be described by the following linear momentum equation based on the concept of effective stress [5,61].
This concept states that the effective stress σ , which acts on the subsurface rock is composed of the total stress σ and the mean fluid pressurep. Using the conventional stress definition with positive stress for tensile stress, the effective stress is defined by Biot's constant b defines the amount of change in bulk volume related to changes in fluid pressure [6].
where K b is the drained bulk modulus and K s is the bulk modulus of solid. For K s = ∞, the solid matrix is incompressible. To account for multiple phases (n phases), mean fluid pressurep and mean densityρ are calculated bȳ where p α is pressure and S α is saturation of fluid phase α.
For the stress-strain relationship in terms of the effective stress, the following constitutive law is applied [70].
where the material tensor D describes the linear elastic material constitutive. For an isotropic material, it is defined as [70]: With the solid dependent material constants defined as λ = E·ν (1+ν)·(1−2·ν) and G = E 2(1+ν) . The strain-displacement relationship for linear elasticity can be written as [31]: Here, displacement vector u is the primary variable to be solved. Strain ε is then derived from u.

Fluid flow process
Following [32], the balance equation for one phase in a multiphase deformable porous medium based on Biot's consolidation concept [5] is governed by The compressibility of the system is expressed by modulus K α for the fluid phase and K s for the solid. Geomechanical feedback on fluid flow is considered by the volumetric strain rate ∂ε ∂t . It causes a change in the pore volume, which accordingly influences the pressure of the pore filling fluids. The ratio of mechanical feedback to fluid flow is expressed by Biot's coefficient. Effects of non-isothermal processes are neglected. All phases are related to each other by capillary pressure-saturation-relationships. The individual flux of one phase can be written as an extended Darcy's law [32] where k r,α expresses the relative permeability of one phase derived from permeability-saturation relationships.

Hydromechanical parameters
Deformation-induced changes of permeability k and porosity φ can be approximated using empirical relationships between stress and porosity [51].
whereσ is calculated from the three effective principle stresses Bases on the porosity change, permeability k can be updated as follows: The presented formulations of deformation and fluid flow form a complex and strongly coupled differential equation system, which can be solved by numerical approximation methods. Depending on the problem, different coupling methods can be used as presented above.

General coupling scheme
The general process coupling structure between Open-GeoSys and Eclipse is based on the process-oriented approach of OpenGeoSys. A detailed description can be found in [44]. All processes including flow, transport, deformation and geochemical reactions (THMC processes) are solved consecutively in a coupled or uncoupled manner. For the coupled OpenGeoSys-Eclipse simulator, Eclipse is integrated into the OpenGeoSys process structure within the flow and mass transport processes and can be used as an alternative simulator for these processes [44]. Flow processes in Eclipse can be either single-phase or multiphase (E100) or multiphase multicomponent (E300). Eclipse is coupled to OpenGeoSys using an operator splitting approach. Hereby, OpenGeoSys controls the simulation configuration and schedule and defines the time step lengths. For each time step, the Eclipse executable is called within the OpenGeoSys simulation to perform the flow process. The results of Eclipse are then passed through the interface to OpenGeoSys, where the remaining processes are calculated for the current time step. For each individual process, a feedback to the flow process in Eclipse can be incorporated through the THMC feedback unit. Before running the next time step in Eclipse, the data entries of Eclipse are updated accounting for the chosen process feedbacks.
As mentioned above, OpenGeoSys and Eclipse use different numerical schemes. OpenGeoSys uses the finite element method, where the numerical grid consists of either 1D (lines), 2D (triangles, quadrilaterals), or 3D (tetrahedron, hexahedron) elements to represent a certain geometry. Eclipse on the other hand uses the finite volume method. It provides high flexible grids consisting of regular or irregular (collapsed or distorting) hexahedrons [55] and also allows for non-neighbor connections. For the coupled OpenGeoSys-Eclipse simulation, an identical mesh for both simulators is required. Therefore, a mesh converter has been developed which converts even complex Eclipse grids to the OpenGeoSys grid structure [19,66]. Although both simulators use the same grid for a coupled OpenGeoSys-Eclipse simulation, they store the numerical solutions at different element positions. In OpenGeoSys, values are stored at the element nodes, whereas in Eclipse, values are stored at the element center. Therefore, the transferred data at the element centers of Eclipse is interpolated to the element nodes of OpenGeoSys using the inverse volume weighting approach [44]. Also phase velocities from the element faces of the Eclipse grid are transferred to the Gauss points of the respective OpenGeoSys grid elements. Correspondingly, the results of OpenGeoSys are interpolated back to the Eclipse grid before the next Eclipse run.

Geomechanical feedback to the flow equation in OpenGeoSys
First, the method of geomechanical feedback on fluid flow in OpenGeoSys is presented. The feedback is done through the volumetric strain rate dε dt , which is added as source term to the fluid flow equation (see Eq. 9). Following [32], the strain rate is defined as where ε is the volumetric strain and v s the solid velocity. The strain rate is therefore the change in solid velocity, meaning the change of the rate of displacement u of the porous medium with time.
The equation system in OpenGeoSys is solved for the primary variable displacement u using the Galerkin finite element method [31]. This method of weighted residuals converts the differential equations (Eqs. 1 and 9) to weak formulations, which are then approximated by trial solution using the values at the element nodes and interpolation functions (shape functions). For a precise solution of the deformation process, displacement u is approximated using quadratic shape functions.
whereũ i is the approximated value of displacement u at node i, N u i is the quadratic shape function of node i, and NN is the number of quadratic nodes [31]. The discretized weak form of strain change is then calculated by where N j is the linear shape function of element j . By inserting Eq. 16 into Eq. 17, the calculation of the change of the strain rate ∂ε for one element j can be written as Applying the method of weighted residuals, the numerical approximation of Eq. 18 is then (19) with w k as the weighting function, j k the Jacobian matrix (mapping global to local coordinates), and NG the number of Gauss points of the finite element. The interpolation functions are merged into a large coupling matrix. It is therefore appropriable for monolithic and sequential coupling.
The strain rate at each node for the current time step or iteration is calculated by Eq. 20 and added as source term to the fluid flow equation system of the next time step or iteration.

Geomechanical feedback to the flow equation in OpenGeoSys-Eclipse
The effect of strain changes on fluid flow, which is expressed by bS α ∂ε ∂t in Eq. 9, is not considered in the fluid flow equation of a pure Eclipse simulation. Therefore, to account for geomechanical feedback on fluid flow in the OpenGeoSys-Eclipse interface, the volumetric strain rate calculated by OpenGeoSys is translated into a fluid pressure change in Eclipse. This is done by a strain equivalent pressure correction, which is derived for the mass balance equation for one phase α in a deformable porous medium (see Eq. 9).
By summarizing compressibilities to a storage term S T and neglecting mass flux and source terms, Eq. 21 reduces to For one time step t, Eq. 22 can be rewritten as Solving Eq. 23 for the new time k + 1, it becomes where the right term represents the strain equivalent pressure correction. As the compressibility is the denominator of the pressure correction term, a compressible system is essential for the numerical stability of this coupling approach.
As Eclipse uses the finite volume method, the strain rate ε (Eq. 24) calculated in the geomechanical process of OpenGeoSys is sum up for each element instead of interpolating the values to the finite element nodes. Thus, Eq. 19 reduces to The strain changes are calculated using Eq. 25 and inserted into Eq. 24 to calculate the pressure correction. The corrected pressures are then added to the current pressure values in Eclipse. For multiple phases, the pressure correction is applied to all phase pressures factorized by their phase saturation. The updated pressure values are then transferred back to Eclipse to perform the next time step or iteration.
The pressure correction procedure is implemented within the interface of OpenGeoSys-Eclipse. It can be switched on or off depending on the coupling method applied.

Porosity and permeability update in OpenGeoSys-Eclipse
Besides the strain-related pressure correction, mechanical feedback on flow also includes the update of hydromechanical parameters, because mechanical load leads to deformation of the porous medium and therefore causes changes in porosity and permeability. To incorporate this effect in the hydromechanical coupling of OpenGeoSys-Eclipse, new porosity and permeability models have been implemented in OpenGeoSys accordingly to the presented porosity-mean effective stress relationship of Eq. 11 and permeability-porosity relationship of Eq. 13. These empiric formulations can be used either for a coupled OpenGeoSys-Eclipse simulation or for a single OpenGeoSys simulation. Figure 1 illustrates the implemented coupling scheme of coupled hydromechanical simulations in OpenGeoSys-Eclipse. Fluid flow and deformation are hereby solved accordingly to the general coupling scheme presented in [44]. First, fluid flow is performed in Eclipse and the updated pressures, saturations, and densities are transferred to OpenGeoSys. OpenGeoSys then calculates the stress field taking into account changes of the mean fluid pressure. Based on the stress and strain changes, pressures, porosity, and permeability are corrected and transferred back to Eclipse. For the next time step or iteration in Eclipse, the corrected pressure values are loaded by a restart functionality and the newly calculated porosities and permeabilities are imported through a grid parameter file.

Numerical procedure
The implemented coupling routine enables multiple coupling methods with increasing coupling accuracy: The explicit coupling, where geomechanical feedback on flow is neglected, so that no feedback parameters are transferred back to Eclipse; the non-iterative coupling, which incorporates the correction of fluid pressure, porosity and permeability after each time step and the iterative coupling, where iterations are performed within each time step and the feedback parameters are updated with every iteration until a defined convergence criterion is reached. The last one is the most accurate method, because

Model verification
The hydromechanical coupling method implemented in OpenGeoSys-Eclipse is verified by several test cases in 2D and 3D using the explicit and iterative coupling methods (Table 1). All OpenGeoSys-Eclipse results are compared with those of a pure OpenGeoSys simulation.
As OpenGeoSys has been extensively tested before, this comparison satisfies the validation process [21,29,30]. Both, single-phase and multiphase flow are tested to verify the correct data transfer and averaging of the pressure and saturation data. In addition to the method verification, the different coupling methods are discussed with regard to their impact on the simulation results. For an adequate comparison, the same mesh dimensions and elements are used in both simulators for all test cases. As Eclipse can only handle hexahedron element types, the 2D test case grids consist of three-dimensional hexahedrons, but can be treated as quasi-2D models as they comprise only one cell in the third direction (Y-direction). The 3D model contains multiple cells in all directions.

Test case 1: Brine injection into a saline aquifer (2D)
The first test case represents a brine injection into a deep non-faulted saline aquifer. The model has an extent of 200 m and a thickness of 6 m using a discretization of 5 m and 1.2 m, respectively (Fig. 2). Petrophysical parameters and fluid properties are given in Table 2. Impermeable formations are assumed above and below the aquifer, so that the injection-induced pressure propagation is restricted to the horizontal direction. The upper, left and lower model boundaries are set to zero displacement conditions for the X-and Z-directions, respectively. A constant hydrostatic pressure gradient is defined at the right model boundary. Initially, hydrostatic pressure gradient is given as listed in Table 3. Initial stresses are neglected as the test cases focuses on injection-induced effective stress changes. For In a first step, test case 1 is simulated using the explicit coupling without rock compressibility to compare the simulators for a simple incompressible fluid flow system. As deformation-induced pressure changes are compensated though compressibility (see Eq. 24), test case 1 is simulated in a second step using both the explicit and the implicit coupling method including compressibility. The simulation results of fluid pressure, effective stress, and strain are evaluated and compared for both simulators OpenGeoSys and OpenGeoSys-Eclipse.
The results of pressure and effective stress in X-direction of OpenGeoSys-Eclipse are depicted in Figs. 3 and 4. For the explicit coupling without rock compressibility, water injection leads to a pressure increase within the aquifer and water displacement towards the open model boundary on the right side (Fig. 3a). A steady-state pressure distribution is reached immediately after injection starts with a maximum pressure of 33.5 MPa at the injection location and a maximum effective stress change of 16 MPa correspondingly to the injection-induced pressure increase (Fig. 4a). By setting bulk modulus K b to 2.22 MPa, pressure increase after 10 days only amounts to 18.5 MPa due to pressure compensation by rock compressibility (Fig. 3b). With time, pressure increases and propagates further into the model domain. After 40 days, pressure is increased to 19.6 MPa at the injection well (Fig. 3c). This is an increase of 2.3 MPa, i.e., 14% of the pressure build-up found in the simulation without compressibility. Effective stress in Xdirection reacts correspondingly, with the highest changes of 2.3 MPa at the injection well after 40 days (Fig. 4b  and c).  The results of OpenGeoSys and OpenGeoSys-Eclipse are compared in Fig. 5. They agree well for all three simulations, whereby the best comparison is obtained for the explicit coupled simulation without compressibility (Fig. 5a). For the explicit coupled simulation with compressibility, the simulators calculate the same decreasing pressure gradient along the model domain. At 5-m distance from the injection well, a small difference of 0.085 MPa between the values indicates a slightly stronger pressure compensation by compressibility in OpenGeoSys (Fig. 5b). This is a difference between the results of 3.7% relative to the absolute pressure build-up. At 30-m distance from the injection well, a similar difference between OpenGeoSys and OpenGeoSys-Eclipse of 0.078 MPa is observed, i.e.,  9.3% relative to the absolute pressure increase at this point. Similar differences are observed for the iterative coupling with compressibility, as the deformation-induced pressure changes are compensated by compressibility, which systematically leads to deviating results between the two codes.
The results of effective stress in X-direction are identical to those of the pressure with the best agreement for the explicit coupling and similar differences for the explicit and iterative coupling with compressibility ( Fig. 5 a and b). The comparison between explicit and iterative coupling including compressibility indicates a generally low impact of the mechanical feedback on flow for the given model setup. At 5-m and 30-m distance from the injection well, the differences in pressures between the explicit and iterative coupling are less than 1% relative to the general pressure build-up at these points in both simulators.

Test case 2: Radial brine injection into a saline aquifer (3D)
To test the implemented feedback for the three-dimensional case, a 3D model is prepared representing a radial injection scenario into a deep saline aquifer. It has a lateral extent of 10 by 10 m with increasing cell sizes from 0.2 m at the injection well to 1 m at the model boundaries and a thickness of 1m with a cell size of 0.5 m. The model dimensions and boundary conditions are illustrated in Fig. 6. The fluid and rock parameters are taken from Table 2. Petrophysical parameters are isotropic and homogeneously distributed within the model domain. Constant hydrostatic boundary conditions are assumed around the model domain simulating infinite open boundaries. For the geomechanical process, the lower and lateral model boundaries are set to zero displacement conditions. The upper boundary is free to move. The initial hydrostatic pressure gradient and total stress gradients in vertical and horizontal direction are defined accordingly to Table 3 assuming compressive stress regime. Water is injected in the model center for 10 days with an injection rate of 30 m 3 per day.
Test case 2 is simulated in OpenGeoSys and OpenGeoSys-Eclipse using both the explicit and iterative coupling methods. The iterative simulations are performed with varying rock materials E and ν (Table 2, column 2) Fig. 6 Model setup of test case 2 showing model geometry, boundary and initial conditions (blue: pressure gradient p; brown: vertical stress gradient σ zz ) and the injection location. The dashed blue circles represent the schematic pressure build-up generated by the injection representing two different rock types, with m1 (E = 10 MPa, ν = 0.3) representing a low deformable rock and m2 (E = 1 MPa, ν = 0) a highly deformable rock. The latter one increases rock deformation and with that geomechanical feedback to fluid flow. All fluid and rock parameters required for this simulation are reported in Table 2, column 2. The explicit coupled simulation is only performed with rock material m1, as the geomechanical feedback is neglected and different rock parameters therefore do not affect fluid pressure.
The results of both simulators after 10 days of simulation are shown in Fig. 7. As water is injected, pressure builds up at the injection location and propagates radially into the aquifer. For the explicit coupling, a maximum pressure of 26.25 MPa is reached at the injection well after 10 days (Fig. 7a). The injection-induced increase in pressure results in a general reduction of the effective stresses as can be seen for the vertical effective stress in Z-direction (Fig. 7b). It reaches a minimum of − 7.86 MPa at the injection location after 10 days. For the iterative coupled Fig. 7 Results of OpenGeoSys (symbols) and OpenGeoSys-Eclipse (lines) for test case 2 showing the distribution of a pressure and b effective stress in Z-direction along a line through the well along an X-transection at 1500.5-m depth after 10 days of simulation simulations, pressure and stress changes are smaller as the geomechanical feedback leads to an increase in pore volume for the given test case. As a result, fluid pressure increase within the pores is less, which in turn results in less decrease of the effective stress. This is most significant for the high deformable medium m2. Here, the maximum pressure at the injection well is reduced by 1.1 to 25.15 MPa and effective stress is decreased to − 9.25 MPa. In contrast, the pressure reduction by the iterative coupling using material m1 amounts to 0.07 MPa compared with the explicit coupled simulation.
The results of OpenGeoSys and OpenGeoSys-Eclipse show the same behavior in pressure increase and effective stress decrease with time ( Fig. 7 a and b). Small differences between the results arise out of the different pressure compensation by compressibility within the two simulators as already seen for test case 1. For the explicit coupling, the difference at 1-m distance from the injection location amounts to 0.05 MPa, which is 0.62% of the overall pressure build-up of 8.05 MPa. For the iterative coupling, the difference is 0.08 MPa for material m1 and 0.34 MPa for material m2 corresponding to 1.04 and 4.77%, respectively. This shows an increase in error between the simulator results with increasing coupling strength. For model setups with low impact of compressibility on the overall pressure build-up (the explicit coupling and iterative coupling m1), this effect is small, whereas for setups with strong geomechanical feedback and high influence of the compressibility, e.g., for the iterative coupling m2, the effect of compressibility and therefore the error between OpenGeoSys and OpenGeoSys-Eclipse gets more visible.

Test case 3: Radial CO 2 injection into a saline aquifer (3D)
Test case 3 is used to test the geomechanical feedback in a coupled multiphase flow and deformation simulation. Hereby, mean pressurep used for the calculation of the deformation process is a mixture of the phase pressures and saturations of all phases. Therefore, this test case checks for correctness of pressure averaging and data transfer between OpenGeoSys and Eclipse. The same conceptual model is used as in test case 2 (Fig. 6). The injected phase is a gas phase simulating CO 2 injection into a saline aquifer. Phase and material parameters are given in Table 2, column 3. All phase properties, i.e., densities and viscosities of brine and CO 2 , are assumed to be constant and dissolution is not considered. These assumptions are not realistic for site specific applications but are made to be able to clearly discern the coupling effects between fluid pressure and rock deformation, as these are within the focus of this manuscript. Boundary and initial conditions for the fluid flow and deformation processes are adopted from Fig. 6 and Table 3. Initially, the aquifer is filled with water. CO 2 is injected with a rate of 1.3 m 3 per day for 10 days. The model is simulated with OpenGeoSys and OpenGeoSys-Eclipse using again the explicit and iterative coupling methods. To increase deformation and therefore geomechanical feedback on fluid flow, a high deformable rock is assumed.
As CO 2 is injected, it spreads radially into the aquifer and leads to a concentric pressure increase around the injection location ( Fig. 8 a and b). Due to the increase in pressure, effective stress in Z-direction is concurrently reduced with Fig. 8 Comparison of the simulation results of OpenGeoSys (symbols) and OpenGeoSys-Eclipse (lines) for test case 3. a Pressure, b gas saturation, and c effective stress in Z-direction along a line through the well along the X-direction at 1500.5-m depth after 10 days of injection. The result of the gas saturation of the iterative coupling case in OpenGeoSys-Eclipse is additionally marked by red crosses as there is only a small, invisible difference to the explicit coupling case a maximum decrease at the injection location (Fig. 8c). Towards the outer model boundaries, pressure build-up tends towards zero due to the open boundary conditions and so effective stress changes are zero. After 10 days, the CO 2 saturation at the injection well reaches a maximum of 0.42 in OpenGeoSys-Eclipse and 0.51 in OpenGeoSys.
Comparing the iterative and explicit coupling methods, pressure build-up and stress reduction is slightly reduced due to the pore volume coupling. This effect is strongest at the injection well and decreases towards the far field correspondingly to the magnitude of pressure and stress changes. In OpenGeoSys, the difference between iterative and explicit coupling is 0.8 MPa at the injection well and less than 0.1 MPa in the far field of the aquifer. Because changes in pressures are small between iterative and explicit coupling, no visible differences in the CO 2 saturation distribution can be observed. A stronger change in pressure would consequently impact on the saturation distribution as capillary pressure is affected by these changes.
Comparing OpenGeoSys and OpenGeoSys-Eclipse, a similar pressure build-up and stress reduction is observed for the explicit coupling (Fig. 8a, c). Small differences between the results occur due to a steeper pressure and saturation gradient in OpenGeoSys indicating less numerical dispersion compared with OpenGeoSys-Eclipse. This results in a higher maximum CO 2 saturation at the injection well in OpenGeoSys. Regarding the iterative coupling, the effect of geomechanical feedback on flow is stronger in OpenGeoSys as pressure increase is less than in OpenGeoSys-Eclipse (Fig. 8a) demonstrating the stronger impact of compressibility in OpenGeoSys. At 1-m distance from the injection well, the difference in pressure between both simulators is 0.01 MPa for the explicit coupling corresponding to 3.18% of the absolute pressure build-up, and for the iterative coupling, the difference amounts to 0.048 MPa, which is 11.72%. Also, the decrease in effective stress changes causes by the iterative coupling is less for the OpenGeoSys simulation ( Fig. 8c) according to higher pressure compensation.
However, although the impact of the geomechanical feedback is generally lower in OpenGeoSys-Eclipse due to the different compressibility models, both simulators calculate the same dampening effects by the iterative pore volume coupling.

Porosity and permeability update
The implementation of stress-dependent porosity and permeability changes within OpenGeoSys-Eclipse is verified in test case 4 using the model setup depicted in Fig. 6. Porosity and permeability values are given in Table 2 (column 4), parameters for permeability and porosity update in Eqs. 11 and 13 are taken from the study of [51] representing properties of a typical sandstone. The test case is performed in both simulators with and without porosity and permeability update, respectively. The implemented update equations will be verified as well as the data transfer of porosity and permeability between OpenGeoSys and Eclipse. Compressibility is neglected for all simulations to avoid overlapping of the hydromechanical coupling effects.
The results of pressure, mean effective stress, vertical displacement and deformation-induced changes in porosity and permeability after 10 days of injection are shown in Fig. 9. Pressure is increased by the injection process and reaches a maximum pressure at the injection well of 20.13 MPa in the OpenGeoSys-Eclipse simulation and 20.22 MPa in the OpenGeoSys simulation (Fig. 9a, simulations without update). Effective stresses are reduced correspondingly to the pressure increase as can be seen for the mean effective stress in Fig. 9b, which is decreased to − 7.39 MPa in the OpenGeoSys simulation and − 7.48 MPa in the OpenGeoSys-Eclipse simulation. As the model top is set to unconstrained conditions for the geomechanical process (Fig. 6), uplift can be observed being highest at the injection well. Vertical displacement at the well amounts to 0.0043 m in OpenGeoSys and 0.0041 m in OpenGeoSys-Eclipse (Fig. 9c).
Calculations of porosity and permeability change show an increase proportionally to the change in effective stress ( Fig. 9 d and e). This is most significant near the injection well, where the strongest changes in pressure and effective stresses appear. In OpenGeoSys, porosity increases by 0.0068 from 0.206 to 0.213, i.e., 3.3%, and permeability by 1.1·10 −13 m 2 from 1·10 −13 to 2.10·10 −13 m 2 , i.e., 110%. Due to the slightly smaller pressure increase in OpenGeoSys-Eclipse, also porosity and permeability changes are smaller. Porosity is increased by 0.0062 corresponding to 3.0% and permeability by 9.48·10 −14 m 2 , i.e., 94.8%. At 3-m distance from the injection well, the impact of water injection and therefore pressure increase is less. As a result, the differences of the results of pressure changes between both simulators are smaller and with that also differences in flow properties. Porosity is increased by 0.0013 in OpenGeoSys and by 0.0012 in OpenGeoSys-Eclipse, i.e., 0.64% and 0.58%, respectively. Permeability is increased by 1.15·10 −13 m 2 in OpenGeoSys and 1.14·10 −13 m 2 in OpenGeoSys-Eclipse corresponding to 15% and 14%. Towards the outer model boundaries, changes in pressure and effective stresses are less, therefore also the changes in porosity and permeability are closed to zero.
Regarding the effect of porosity and permeability update on fluid flow, a generally smaller pressure increase is observed for the simulations including the update (Fig. 9a). Compared with the simulation without update, pressure increases up to 19.02 MPa in OpenGeoSys-Eclipse and Fig. 9 a Pressure, b mean effective stress, c vertical displacement, d porosity changes, and e permeability changes along a line through the well along the X-direction at 1500.5-m depth after 10 days comparing the results of the simulations with and without the porosity and permeability update in OpenGeoSys (symbols) and OpenGeoSys-Eclipse (lines) to 19.06 MPa in OpenGeoSys, which is a reduction of 38.76% and 39.27%, respectively. As a result of the smaller pressure increase, also the decrease in mean effective stress is less (Fig. 9b). It can be seen from Fig. 9 that both simulators OpenGeoSys and OpenGeoSys-Eclipse produce similar results, which verifies the implemented update method. Only small differences occur near the injection well arising from slightly different pressures at the well as well as the interpolation method and data transfer between OpenGeoSys and Eclipse, as already detected for different test cases by [44].

Computation speed
As mentioned above, OpenGeoSys is coupled to Eclipse to benefit from its fast and efficient computation time of flow processes. To evaluate the time savings by this coupling, simulation times of the presented test cases are compared in Table 4. For small simple single-phase models (case 1), no simulation speed-up is archived by the OpenGeoSys-Eclipse simulator. The data transfer through the Eclipse interface is time consuming compared with a single OpenGeoSys simulation without any data transfer. Increasing the dimensions and phases in the model, OpenGeoSys-Eclipse shows significant faster run times. Test case 2 is speed up by a factor of about 1.5 and test case 3 by a factor of 5.6 for the iterative coupling. Regarding potential large-scale and long-term scenario simulations of gas storage, this time speed-up archived by the OpenGeoSys-Eclipse coupling will significantly reduce simulation times.

Discussion and conclusion
Numerical simulation codes are a useful tool to assess the complex process system related to subsurface operations such as gas storage. The simulation and prediction of geomechanical processes is of particular importance for the operation safety as pressure disturbances associated with the injection or extraction of mass can significantly change the subsurface stress regime and causes rock deformation. Adequate simulation tools are required capable of taking into account each individual process as well as the process coupling to investigate the impact of process feedback to the system. Therefore, the coupled simulator OpenGeoSys-Eclipse, developed to simulate coupled thermal, hydraulic, geomechanical and geochemical processes in the subsurface, has been extended for the hydromechanical process coupling, whereby Eclipse solves for fluid flow and Open-GeoSys for deformation. The newly developed process feedback incorporates two feedback components: The pore volume coupling, which is translated into a pressure correction term in Eclipse and the stress-dependent porosity and permeability update based on changes of the mean effective stress. Two different coupling methods, the explicit and the iterative coupling method, are implemented to simulate the feedback. The feedback component as well as the type of process coupling can be chosen for the individual OpenGeoSys-Eclipse simulation demonstrating the very flexible and effective handling of this simulator.
The hydromechanical coupling in OpenGeoSys-Eclipse is verified using several test cases in 2D and 3D with increasing complexity. This is important as the transfer parameters between OpenGeoSys and Eclipse vary with the number of fluid phases in the model, the feedback components and the model dimensions. Differences between the simulation results arise for the volume coupling, which are mainly caused by different handling of compressibility in the individual simulators. Thus, pressure build-up by injection is compensated differently leading to varying pressure signals in the simulation results, whereby Eclipse shows a generally lower pressure compensation by compressibility than OpenGeoSys for all test cases. As the impact of compressibility depends on the model setup, differences between OpenGeoSys and OpenGeoSys-Eclipse vary between 1 and 10% relative to the overall pressure build-up for the given test cases. The difference between the simulation results caused by compressibility systematically increases with increasing impact of geomechanical feedback on flow as seen, e.g., for test case 2. However, the geomechanical coupling causes the same dampening effect on the pressure build-up in both simulators.
Besides the verification of coupled hydromechanical simulations in OpenGeoSys-Eclipse, the test cases are used to evaluate the differences between explicit and iterative coupling and their impact on the simulation results. By this, the advantages and disadvantages of each method can be evaluated, especially with regard to simulation accuracy and efficiency.
The simulations including pore volume coupling (test cases 1 to 3) show generally less differences between the explicit and the iterative coupling for the given model setups. This is due to the fact, that the model parameterizations are related to a deep underground gas storage featuring low compressibility and a less deformable rock. The induced changes in effective stresses and strains are small and with that also the geomechanical feedback on flow. By increasing the deformability of the rock (e.g., test case 2, material m2), a stronger geomechanical feedback to the flow is observed leading to a significantly smaller pressure build-up. Similar results are found when increasing the compressibility of the porous medium (not shown here). The parameter sets used in the test cases were chosen to enable a clear model verification and comparison between the different coupling approaches. Both the use of an incompressible CO 2 phase and the use of independent values for bulk modulus, Youngs modulus, and Poisson ratio, which are related to each other through E = 3K(1 − 2ν) [15], would not be valid for a site-specific study. In test case 2 and for material 2, this would lead to a higher compressibility and thus a somewhat lower pressure increase, while for the other results shown here changes would be small. It can be seen, that the impact of deformation on flow is strongly controlled by the individual reservoir rock and reservoir conditions. Therefore, iterative coupling should be preferred for high compressible and high deformable media to archive adequate accuracy in the simulation results. For low compressible and low deformable systems, the explicit coupling produces sufficient results within a faster run time. For test case 2 with material m1, the simulation time is 0.05 h for the explicit coupling and 0.32 h for the iterative coupling, although the difference between the simulation results is less than 0.1 MPa. However, it is probable that the accuracy of the explicit coupling may decrease for bigger time step sizes as the change in pore volume per time step increases as shown by [60].
In contrast to the simulation with pore volume coupling, the porosity and permeability update has a stronger effect on the pressure build-up using the same model setup and parameterization (test case 4). An increase of permeability by maximum factor of about 2 is observed within the model area as the injection-induced changes in effective stresses cause an elastic expansion of the rock. This leads to a less pressure increase compared with the simulation without update. Similar changes in permeability of a factor of 1.3 to 1.7 were already obtained by [51] using the TOUGH-FLAC3D simulator.
In general, geomechanical feedback should be considered when assessing subsurface operations. Depending on the given storage site conditions and rock types, the feedback of deformation on the pressure can be significant. This is an important fact regarding the risk analysis of underground storage processes, e.g., the potential of ground uplift and hydraulic fracturing. The presented test cases show, that the feedback has an overall positive effect to the storage operation. For both coupling components, the pore volume and the stress-dependent porosity and permeability, pressure build-up in the reservoir tends to be smaller. This would allow for a higher injection flow rate without risking an increase of the injection pressure and with that hydraulic fracturing. Rutqvist et al. [51], for example, presented a reduction of 3 MPa in pressure in the long term of a realistic storage scenario, when considering geomechanical feedback. The reduction of the pressure build-up will positively act on geomechanical processes as changes of the stress field are reduced and therefore also deformation and vertical uplift (e.g., Fig. 9 b and c). To investigate these effects on a larger scale, more realistic storage scenario simulations are needed. For this purpose, the OpenGeoSys-Eclipse simulator offers an appropriate tool, as it provides flexibility regarding the process coupling and a fast run time speed, especially for multiphase flow simulations.