Verification benchmark for a single-phase flow hydro - mechanical model comparison between COMSOL Multiphysics and DuMuX

Numerical modelling of hydromechanical processes in geological environments has become an invaluable tool in understanding and predicting system behaviour. However, due to the different algorithms and numerical schemes implemented in the different models, model reliability may vary considerably. Modelling of single- and multi-phase flow in porous media has been widely employed in various engineering applications such as geological disposal of nuclear waste, geological storage of carbon dioxide, high-temperature geothermal systems, or hydraulic fracturing for shale gas exploitation. Coupled hydro-mechanical (H-M) processes play a key role in the prediction of the behaviour of geological reservoirs during their development and testing operations. In this paper, we present a benchmark test on a single-phase flow problem in a hydrogeological reservoir with 5 horizontal layers of different properties. The aim is to compare two hydromechanical (H-M) models that use a vertex-centred finite-volume discretization and a finite element discretization. The first model is constructed with the free-open source simulator DuMuX, and the second with the commercial software COMSOL Multiphysics. The verification study suggests general confidence in the model reliability, but also highlights and discusses several areas of discrepancies between two models.


Introduction
Hydrogeological problems are generally difficult to study directly because the subsurface structures are not accessible to direct studies and characterisation. Thus, modelling is a widely utilized tool for understanding the influence of complexity in geometry, heterogeneity in the respective properties (mechanical, hydraulic, etc.) as well as the effect of process coupling, such as thermo-hydromechanical-chemical (THMC) processes. Various numerical simulators such as TOUGH-FLAC (Rutqvist, 2013), CODE-BRIGHT (UPC, 2002) and DuMu X (Flemisch et al. 2011) etc. were developed by different groups around the world, contributing to the enrichment of the field of hydrogeological modelling. However, due to different approaches in solving the respective differential equations, programming languages, and the design algorithms of the different models, differences between modelling results can be expected for the same problem. Several verification benchmarking initiatives for simulators of flow and transport in (fractured) porous media are reported in the literature, e.g.., (Ebigbo et  This paper presents the results of a benchmarking study on the hydro-mechanical (H-M) effects induced by a single-phase fluid injection into a fully saturated multi-layered geological system. In this sense, the same problem is implemented in two numerical simulators, i.e., COMSOL Multiphysics and DuMu X and the results are compared based on the state variables overpressure, vertical displacement and change in effective stress. DuMu X is an academic code which is available free of charge and public-domain. The hydro-mechanical model has been originally developed by Darcis (2013) and (Beck 2019 In contrast to the fully coupled DuMu X model, COMSOL model is sequentially coupled. First, at one time-step, the flow solution is calculated and then the effective pore pressure is used for the calculation of stresses and mechanical deformations.

Model setup
The benchmark problem, proposed initially by (Vilarrasa et al. 2020) considers flow and mechanical processes in a porous geological material perturbed by the injection of water into the multi-layered system. The mechanical stress-strain relationships assume linear elasticity, and the problem is isothermal. The reservoir and the other layers in the domain, i.e. cap-rock, under-burden, the upper aquifer and the basal aquifer ( The comparison of the two simulators is further provided as breakthrough curves (plots of parameter of interest per time) at different observation points, and as plots of parameter of interest over horizontal planes at different depth. ) are assumed to be homogenous and isotropic with respect to their hydraulic and mechanical properties. Water is injected into the initially fully saturated reservoir located between depths of 1500 m and 1530 m, with an injection rate of 30 kg/s. The injected water temperature is 60° C i.e. the same temperature as that in the reservoir and the entire domain. The reservoir and surrounding geological structures are assumed to be horizontal. The simulated domain is 1000m below the ground and the initial hydrostatic pressure on the top of the domain is 10 MPa, the densities of reservoir and surrounding rocks are set to a constant 2400 kg/m 3 , which means the lithostatic pressure imposed by the rocks above the domain is 24 MPa. The reservoir is bound by two impermeable layers, i.e. caprock and underburden. Since in reservoir application studies the assessment of caprock integrity is of higher importance than potential H-M induced changes in the reservoir formation, this model includes two more layers above the caprock and below the underburden. The initial water pressure is hydrostatic and maintains a constant pressure at the outer boundary, placed at 10 km distance from the injection well. The initial stress field coincides with the vertical effective stress distribution i.e. 24 MPa/km and the horizontal effective stresses equals 0.5 times the vertical effective stress (Rutqvist, 2013). A constant stress equal to the overburden is applied at the upper boundary. The other boundaries have no displacement perpendicular to the boundary.
The parameters used in the model and detailed geometry are presented in Table 1 and The comparison of the two simulators is further provided as breakthrough curves (plots of parameter of interest per time) at different observation points, and as plots of parameter of interest over horizontal planes at different depth. . The hydraulic and mechanical models are coupled, the coupling is done with through the pore pressure and the mechanical response to stress and pressure changes leads to changes in porosity and intrinsic permeability, for DuMu X and COMSOL, the correlations are the same, shown below: in which is the volumetric strain, which is equal to the sum of the axial strains, is porosity, K is the intrinsic permeability, the subscript eff stand for effective and indicates the mechanics-calibrated value, n is the empirical constant set to 15 in this model, e.g., the same as in Beck et al. (2016) and Figueiredo et al. (2020). With the effective porosity and intrinsic permeability, the mass conservation equation can be written: where is the compressibility of water, is the source and sink term, is the Darcy velocity, and and ρ are the dynamic viscosity and density of water, respectively, t is time, ∇ is the gradient of hydraulic head, and g, the gravitational acceleration.
The correlation between strain and pore water pressure is presented below: where is the Biot effective stress in porous media, is the Biot coefficient, ranged 0 (completely dry) to 1 (fully water saturated), is the water pressure, is Kronecker delta, means if i = j, = 1, else, = 0, here it donates that the water pressure only works on normal stress and axial strain ( = ), for shear stress and shear strain ( ≠ ) , water pressure is zero, and are the strain and elastic modulus, respectively, i and j donate the direction of axis, and are the displacement along the direction i and j.
A detailed description of the mathematical model and the numerical implementation in the DuMu X are given in The comparison of the two simulators is further provided as breakthrough curves (plots of parameter of interest per time) at different observation points, and as plots of parameter of interest over horizontal planes at different depth.

Overpressure
The comparison between the two sets of model results, essentially temporal and spatial distribution of water overpressure (pressure difference between initial status and injection) is the essential part of the benchmark study. The fully saturated system reached equilibrium before actual simulation commenced, the correlation between initial strain , and stress , is shown below: where the , is the stress from rock gravity and is steady pore water pressure, during the injection: where is the water overpressure, calculated by subtracting the initial water pressure at a point location, before injection start, from the water pressure at a calculation time t, ∆ is the additional strain due to the injection. In DuMu X , the initial equilibrium condition is set to be the reference, i.e. , = 0, since the system is fully saturated by water ( = 1). As shown in Eq. (7) and (9), is the only variable affecting the additional strain and displacements (elastic moduli are assumed to be constant during simulation), while in COMOSL, the initial strain is not zero i.e. , ≠ 0, indicating the effects from matrix gravity and initial pore water pressure cannot be ignored, the simulated results are determined by not only the overpressure but the initial stress and pore pressure. These two kinds of initial condition are the main difference in algorithm between the two simulators.
We selected four observation points (OPs) for the comparisons of the overpressure changes. Two of the observation points are located in the centre of the reservoir, while the other two are positioned 10 m above the caprock-reservoir interface. Detailed positions of OPs are shown in Table 2 and The comparison of the two simulators is further provided as breakthrough curves (plots of parameter of interest per time) at different observation points, and as plots of parameter of interest over horizontal planes at different depth.
. Figure 2 presents the comparison of overpressures at 100 m (OPs 1,3) and 1000 m (OPs 2,4) away from the injection well, respectively. The overpressures show very good agreement during the whole simulation.

Displacement
Another important variable considered for comparison is displacement obtained by Eq. (7). In DuMu X , the simulated displacement is fully caused by the overpressure ( , = 0 in Eq. (9)), while in COMSOL, stresses and overpressure are both the keys factor  affecting the displacement ( , ≠ 0 in Eq. (9)), for the comparison, the displacement presented below from COMSOL is the difference between the simulated results and initial displacement. The displacement is plotted over the interfaces of reservoir-caprock (RC) and caprockupper aquifer (CUA) at times 30 days, 365 days and 1000 days after the injection started. Figure 3 shows the comparisons between modelled vertical displacements. It can be observed that the profiles as well as the values of these two models are very similar.
Like the comparison of overpressure, the vertical displacements show very good agreement. The displacement is maximum near the injection well and declines exponentially with the increasing distance from the well. The maximum simulated displacements amount to 0.036 m and 0.015 m at the CUA interface and RC interface respectively. Discrepancies are found in nearwellbore regions, which dissipate towards the model boundary. Discrepancies between displacements are larger at early than at late times. Furthermore, the discrepancies are more pronounced at the RC interface compared to the CUA interface.

Stress
Temporal change in stress is assessed to investigate confidence of the two models.  Figure 4 presents the horizontal stress change over time for two OP locations. In OP 1, DuMu X and COMSOL agree well. However, in OP 2, the stress change modelled by COMSOL is lower by 0.2 MPa over the entire simulated period except at the beginning. The reason for the discrepancy is that, in DuMu X , the water overpressure is the only variable affecting the simulation results aforementioned, the strains and displacements are calculated firstly by the simulator and then the stresses are calculated later by the simulated results i.e. strains, however, in COMSOL, the reference is the condition that the stress is zero, i.e. , = 0, in this case, the influence from the initial stress and pore pressure need to be taken into consideration during the whole simulation, the strains are obtained by stresses and pore pressure, The opposite computation order leads to the large discrepancy in OP 2.

Conclusion
This paper provides a comparison study between two numerical simulators to model hydro-mechanical processes induced by fluid injection in a layered georeservoir. The simulators, DuMu X and COMSOL Multiphysics, are using different discretization schemes and are implemented in different frameworks.
A good agreement is obtained from the comparisons of water overpressure, displacement and horizontal stress change. This agreement increases the confidence in the reliability of these two hydro-mechanical simulators.
Benchmarking studies are useful to verify the reliability of models simulating multi-physics coupled processes. The more complex physical processes are involved, the less likely to find analytical solutions it becomes. At the same time, conducting experiments for achieving a model validation, is limited by increased technical complexity and cost. Generalisations based on a single benchmarking study using only two numerical simulators are limited. Different geometrical settings and boundary conditions, lead to differences between the simulated results. These differences can be considerably large. Therefore, further benchmarking studies are essential, and can be aimed to include higher geometrical complexity, or physical complexity, e.g. non-isothermal, multi-phase, etc. Numerical studies such as this one contribute to broaden the overall understanding of HM models and enhance the confidence in their ability to predict.