Modelling ground displacement and gravity changes with the MUFITS simulator

We present an extension of the MUFITS reservoir simulator for modelling the ground displacement and gravity changes associated with subsurface flows in geologic porous media. Two different methods are implemented for modelling the ground displacement. The first approach is simple and fast and is based on an analytical solution for the extension source in a semi-infinite elastic medium. Its application is limited to homogeneous reservoirs with a flat Earth surface. The second, more comprehensive method involves a one-way coupling of MUFITS with geomechanical code presented for the first time in this paper. We validate the accuracy of the development by considering a benchmark study of hydrothermal activity at Campi Flegrei (Italy). We investigate the limitations of the first approach by considering domains for the geomechanical problem that are larger than those for the fluid flow. Furthermore, we present the results of more complicated simulations in a heterogeneous subsurface when the assumptions of the first approach are violated. We supplement the study with the executable of the simulator for further use by the scientific community.

1 Introduction

Reservoir simulation remains an essential area for forecasting various parameters of subsurface exploration and natural flows. The capabilities of the reservoir simulators, i.e., the computer programs for modelling the flows in geologic porous media, are constantly improving. The modern simulators can account for additional physical phenomena and technological processes. For example, the modelling of Darcy flows is often coupled with rather sophisticated approaches for geomechanics, multiphase transport in wellbores, and other processes (Fig. 1). More common in petroleum reservoir simulation is the development and utilisation of universal software packages allowing for conveniently integrated workflows for coupled thermo-hydro-mechanical modelling and history matching of the models (Fanchi, 2006).

Figure 1 Sketch of typical processes in hydrothermal systems (left) and petroleum reservoirs (right). Magma degassing results in a plume of hot magmatic fluid. Near the surface, the fluid can mix with colder meteoric water. The observation points, where the ground displacement and gravity changes can be measured, are shown.

The numerical modelling of hydrothermal systems differs in several respects from that of petroleum reservoirs (Fig. 1). First, the flows occur under a much wider range of temperatures. This necessitates the application of sophisticated models accounting for phase transitions and reactive transport under both low and high temperatures. Also, the plastic deformations of rocks in the brittle-ductile transition zone of hydrothermal systems should be considered. Secondly, the transport in hydrothermal systems, especially the natural one, can be observed through a very limited number of parameters. Usually, only observations at the surface are available. This is incomparable with petroleum reservoirs, where many more observable parameters are available through extensive drilling and seismic studies. The parameters of heat and mass transfer in hydrothermal systems can be estimated by measuring the surface fluxes of magmatic gas and heat. Also, the parameters of the flows can be estimated by measuring the gravity changes and ground displacement. The periods of unrest of more intense magma degassing into overlying rocks can result in the fluid density increasing in extended regions leading to significant gravity changes at the surface, which can be observed by gravimeters mounted on the surface or aircraft (Fig. 1). Similarly, the redistribution of strains and stresses causes significant ground displacement that can reach tens of centimeters in particular cases (Chiodini et al., 2003; Hurwitz et al., 2007; Hutnak et al., 2009; Rinaldi et al., 2010). Such values of the ground uplift (or subsidence) can be observed by either surface measurements or satellite radar interferometry.

In the numerical modelling of hydrothermal systems, a common approach for estimating gravity changes and ground displacement assumes the application of different computer programs that are not directly compatible with each other. The fluid flow is simulated by using hydrodynamic code. Then, the gravity changes are calculated by post-processing the simulation results in separate code, i.e., by summing the contributions of every grid block to the strength of the gravitational field (Todesco, 2009). Such calculations, although very straightforward, can require a significant effort to develop the interface between the computer programs. A similar approach based on summing the contributions of the grid blocks to the magnitude of ground displacement exists for predicting the subsidence or uplift of Earth's surface (Rinaldi et al., 2011). It also requires considerable efforts to post-process the simulation results. For ground displacement, a more comprehensive approach assumes the coupling of hydrodynamic and geomechanical simulators, which again requires time-consuming efforts to program the interface between the simulators (Todesco et al., 2004; Rutqvist, 2011). Therefore, simulation software allowing for all the calculations described above in an integrated framework is in demand. The goal of this work is the development of such software by extending MUFITS. We supplement the simulator with built-in capabilities for calculating gravity changes and ground displacement and verify the development against a benchmark study.

2 The MUFITS extension

2.1 Review of the simulator and new developments

MUFITS has been developed over the past decade. Initially, it was designed for modelling the subsurface storage of CO₂ in saline aquifers and petroleum reservoirs at relatively low temperatures (Afanasyev, 2013, 2015). Then, it was extended for the modelling at much higher reservoir temperatures including supercritical parameters. The simulator was applied for modelling the natural convection at Campi Flegrei caldera down to a depth of 5 km (Afanasyev et al., 2015). A 3-D model of the subsurface near the Solfatara crater was created and calibrated against the measurements of the surface fluxes of magmatic gas and heat as well as the temperature profiles in a few boreholes. Also, MUFITS was applied to the investigation of the high-temperature flows in the kimberlite pipes during their cooling (Afanasyev et al., 2014) and the formation of lenses of hypersaline fluid above degassing magma bodies (Afanasyev et al., 2018). Until now, the functionality of the simulator did not allow for easy computations of gravity changes and ground displacement, which could help in a better understanding of the flows in the described applications and more reliable history matching of the corresponding reservoir models. The development of such a functionality would be a good step towards transforming MUFITS into a universal software package for modelling hydrothermal systems.

As presented in Sect. 2.3, MUFITS now allows for the built-in calculation of gravity changes. The developed extension of the software includes two methods, A and B, for calculating ground displacement (Fig. 2). Option A is built into the hydrodynamic simulator. It employs an analytical solution for the extension source in a semi-infinite elastic medium corresponding to Earth's subsurface. Every grid block is considered such a point source, and the ground displacement is calculated by summing the contributions of all grid blocks, as described in Sect. 2.4. Thus, method A does not require any interface between hydrodynamics and geomechanics, because it is built into the hydrodynamic code. Option B assumes that MUFITS is coupled with geomechanical code through external files. Now, MUFITS is supplemented with Matlab code for axisymmetric modelling the distributions of strains and stresses and a command-line tool (an external script) that converts MUFITS output files into a simpler binary format readable by Matlab. This tool allows Matlab code to read the external files produced by MUFITS; thus, additional efforts for programming the interface are not needed. Furthermore, this open-source tool can be useful for a quick MUFITS coupling with other software. The equations implemented in the geomechanical code are discussed in Sect. 3.

Figure 2 Diagram showing the software options for modelling gravity changes and ground displacement. The interface between the simulators is organised through external files.

2.2 Observation points

To extend MUFITS for the modelling of gravity changes and ground displacement, a new primitive (element) of the reservoir model, namely the “observation point” (OP), is introduced. Every OP is assigned a character name and is characterised by 3 coordinates in space. For example, an OP can correspond to surface or airborne measurements (Fig. 1). The simulator can be set up for automated calculation of the gravity changes and ground displacement in every OP at every moment of time. To ease the reporting of space distributions, OPs can be grouped into networks of equally spaced points along the x, y, and z axes.

The developed options for simulating gravity changes and ground displacement (by method A) are generally designed for 3-D simulations with domains of arbitrary complexity (Fig. 3a). However, a domain symmetry can often be implied in 2-D simulations to speed them up. The influence of the symmetries on parameters in an OP can automatically be accounted for in the following cases. First, axisymmetric simulations are allowed when the flow is calculated only in a fraction of the full circle (corresponding to the opening angle φ). In this case, the x and y coordinates of a point E belonging to the vertical axis of rotation must be specified (Fig. 3b). Secondly, translational symmetry of the domain is allowed when the 2-D flow is calculated only in a plane fraction (of thickness h) of the 3-D space (Fig. 3c). In this case, the x and y coordinates of the horizontal translational vector n_t must be specified. Thirdly, in addition to the previous option, reflection symmetry is allowed by specifying a point B belonging to the reflection plane and the normal to the plane n_b (Fig. 3d). For these symmetries, the simulator can automatically account for other parts of 3-D space that do not belong to the domain.

Figure 3 Possible symmetries that can be accounted for in the calculation of gravity changes and ground displacement by method A. The simulation domain (red) occupies a region of the subsurface reservoir (black).

We denote by ${\mathit{r}}_o=\mathit{\{}{x}_o,y_o,z_o\mathit{\}}$ and ${\mathit{r}}_i=\mathit{\{}{x}_i,y_i,z_i\mathit{\}}$ the position vectors of an OP and a grid block centre, respectively. In accordance with the finite volume method implemented in MUFITS, the gravity changes and ground displacement in an OP are calculated by summing the contributions of every grid block to the magnitudes of these fields. Herewith, the finite size of the blocks is neglected by considering them as point sources placed at r_i. The corresponding equations are discussed below.

2.3 Modelling gravity changes

According to Newton's law of universal gravitation, the gravity change in every OP is calculated using the following equation:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }\mathit{g}=\sum _{i=\mathrm{1}}^N\mathit{\gamma }{\displaystyle \frac{({\mathit{\rho }}_i-{\mathit{\rho }}_{i,\mathrm{0}})V_i}{{\left|{\mathit{r}}_i-{\mathit{r}}_o\right|}^{\mathrm{3}}}}({\mathit{r}}_i-{\mathit{r}}_o),\end{array}$$

where Δg is the gravity change, N is the number of grid blocks in the simulation, γ is the gravity constant, ρ is the bulk density of the saturated porous medium, ρ₀ is the reference bulk density, V is the volume of a grid block, and the subscript i denotes the parameters of the ith block. The change Δg is calculated against the reference state with distribution of density ρ₀. By default, the distribution of ρ at the initial moment of time is set as the reference distribution ρ₀. However, the user can override this setting by specifying moments of time at which the distribution of ρ is copied to ρ₀. Thus, one can choose (or redefine) the moment of time of the reference state.

2.4 Modelling ground displacement (method A)

Method A is restricted by the following assumptions:

• The mechanical properties of the saturated porous medium are homogeneous.

• The top boundary of the domain, z=z_top, corresponding to Earth's surface, is flat, horizontal, and free of stresses.

• The domain for modelling the displacement is a semi-infinite region z≥z_top filled with an elastic medium.

Given that these assumptions are satisfied, the displacements are calculated by employing an analytical solution for the centre of compression (or dilatation) placed in the interior of the semi-infinite solid (Mindlin, 1936; Mindlin and Cheng, 1950). Every grid block is considered such a centre of compression, and the relative change in the grid block volume is calculated as (Rinaldi et al., 2010)

$$\begin{array}{}\text{(2)}& {\mathit{\theta }}_i={\displaystyle \frac{\mathrm{\Delta }{V}_i}{V_i}}={\displaystyle \frac{P_i-P_{i,\mathrm{0}}}{H}}+{\mathit{\alpha }}_{\mathrm{s}}(T_i-T_{i,\mathrm{0}}),{\displaystyle \frac{\mathrm{1}}{H}}={\displaystyle \frac{\mathrm{1}}{K}}-{\displaystyle \frac{\mathrm{1}}{K_{\mathrm{s}}}},\end{array}$$

where P and T are the pressure and temperature, P₀ and T₀ are their values in the reference state, $\mathrm{\Delta }V=V-V_{\mathrm{0}}$ is the change in grid block volume against the reference state, 1∕H is Biot's constant, K is the isothermal drained bulk modulus, K_s is the bulk modulus of the solid phase, and α_s is the volumetric thermal expansion coefficient.

Using the analytical solution of Mindlin and Cheng (1950), the displacement in every OP is calculated as:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{u}_{\mathit{\xi }}& =\sum _{i=\mathrm{1}}^N({\mathit{\xi }}_i-{\mathit{\xi }}_o){\displaystyle \frac{{\mathit{\theta }}_i{V}_i(\mathrm{1}+\mathit{\nu })}{\mathrm{12}\mathit{\pi }(\mathrm{1}-\mathit{\nu })}}\left({\displaystyle \frac{\mathrm{1}}{R_{\mathrm{1}}^{\mathrm{3}}}}+{\displaystyle \frac{\mathrm{3}-\mathrm{4}\mathit{\nu }}{R_{\mathrm{2}}^{\mathrm{3}}}}\right.\\ & \left.-{\displaystyle \frac{\mathrm{6}{\stackrel{\mathrm{̃}}{z}}_i({\stackrel{\mathrm{̃}}{z}}_i+c_i)}{R_{\mathrm{2}}^{\mathrm{5}}}}\right),\mathit{\xi }=x,y,\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}{u}_z& =\sum _{i=\mathrm{1}}^N{\displaystyle \frac{{\mathit{\theta }}_i{V}_i(\mathrm{1}+\mathit{\nu })}{\mathrm{12}\mathit{\pi }(\mathrm{1}-\mathit{\nu })}}\left({\displaystyle \frac{{\stackrel{\mathrm{̃}}{z}}_i-c_i}{R_{\mathrm{1}}^{\mathrm{3}}}}-(\mathrm{3}-\mathrm{4}\mathit{\nu }){\displaystyle \frac{\stackrel{\mathrm{̃}}{z}+c}{R_{\mathrm{2}}^{\mathrm{3}}}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{2}{\stackrel{\mathrm{̃}}{z}}_i}{R_{\mathrm{2}}^{\mathrm{3}}}}-{\displaystyle \frac{\mathrm{6}{\stackrel{\mathrm{̃}}{z}}_i({\stackrel{\mathrm{̃}}{z}}_i+c_i{)}^{\mathrm{2}}}{R_{\mathrm{2}}^{\mathrm{5}}}}\right),\end{array}\end{array}$$

where $\mathit{u}=\mathit{\{}{u}_x,u_y,u_z\mathit{\}}$ is the displacement and the following notations are introduced:

$$\begin{array}{c}{\displaystyle }{\stackrel{\mathrm{̃}}{z}}_i=z_i-z_{\mathrm{top}},c_i=\mathrm{2}{z}_{\mathrm{top}}-z_i,\mathit{\nu }={\displaystyle \frac{\mathrm{3}K-\mathrm{2}\mathit{\mu }}{\mathrm{2}(\mathrm{3}K+\mathit{\mu })}},\\ {\displaystyle }{R}_{\mathrm{1}}^{\mathrm{2}}=(x_i-x_o{)}^{\mathrm{2}}+(y_i-y_o{)}^{\mathrm{2}}+({\stackrel{\mathrm{̃}}{z}}_i-z_o{)}^{\mathrm{2}},\\ {\displaystyle }{R}_{\mathrm{2}}^{\mathrm{2}}=(x_i-x_o{)}^{\mathrm{2}}+(y_i-y_o{)}^{\mathrm{2}}+(c_i-z_o{)}^{\mathrm{2}}.\end{array}$$

Here, ν and μ are Poisson's ratio and the shear modulus, respectively. Equations (3) and (4) are those presented by Rinaldi et al. (2010) with the exception of the sign of the second term within the brackets in Eq. (3).

The displacement u is calculated against the reference state characterised by the distributions of pressure P=P₀ and temperature T=T₀. By default, the distributions of P and T at the initial moment of time are set as the reference. This setting can be overridden by specifying the moments of time at which MUFITS copies the distributions of P and T to P₀ and T₀, respectively. Thus, one can choose (or redefine) the moment of time of the reference state.

3 Geomechanical code (method B)

Method A provides a fast option for calculating ground displacement. However, this mathematical model is restricted to the case of a homogeneous distribution of thermoporoelastic moduli and is subject to other assumptions described in Sect. 2.4. This limits the range of potential applications of the built-in ground deformation model. Several studies (Trasatti et al., 2005; Chiarabba and Moretti, 2006; Manconi et al., 2010) show that the mechanical properties of rocks within hydrothermal systems can exhibit significant heterogeneity, which method A cannot deal with.

To allow for the modelling of hydrothermal systems with heterogeneous mechanical properties, we have developed numerical code for calculating stresses and displacements by the finite volume method in axisymmetric domains. The solution of the elastic problem is uncoupled from the equations governing the fluid flow, i.e., deformations of the solid phase do not influence the fluid pressure or the porosity and permeability of the rock matrix. This is justified by the assumption of small deformations. In addition, uncoupling the ground displacement from the fluid flow makes it possible to implement method B in a post-processing module of the simulator, removing the necessity to re-run MUFITS simulations in order to change the mechanical properties of the rocks.

The constitutive equations for linear isotropic thermoporoelastic medium are expressed as (Coussy, 2004)

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& P_{\mathrm{tot}}-P_{\mathrm{tot},\mathrm{0}}=-K{\mathit{\epsilon }}_{kk}+b(P-P_{\mathrm{0}})+{\mathit{\alpha }}_{\mathrm{s}}K(T-T_{\mathrm{0}}),\\ & b=\mathrm{1}-{\displaystyle \frac{K}{K_{\mathrm{s}}}},\\ & {\mathit{\tau }}_{ij}-{\mathit{\tau }}_{ij,\mathrm{0}}=\mathrm{2}\mathit{\mu }({\mathit{\epsilon }}_{ij}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\mathit{\epsilon }}_{kk}),i,j=\mathrm{1},\mathrm{2},\mathrm{3}.\end{array}\end{array}$$

Here, P_tot and τ_ij are the hydrostatic and deviatoric components of the total stress tensor ${\mathit{\sigma }}_{ij}={\mathit{\tau }}_{ij}-P_{\mathrm{tot}}{\mathit{\delta }}_{ij}$, ε_ij is the infinitesimal strain tensor, b is the Biot–Willis coefficient, and the subscript 0 denotes the quantities in the reference state.

We seek a numerical solution of the equations of static equilibrium:

$$\begin{array}{}\text{(6)}& \mathbf{\nabla }\cdot \mathit{\sigma }=\mathrm{0}.\end{array}$$

There is no gravity contribution in Eq. (6) because we consider the deviations of stresses and displacements from the reference state. For axisymmetric problems, system (6) is solved in a cylindrical coordinate system $(r,\mathit{\phi },z)$, and the displacements in the direction of the angular coordinate, u_φ, are assumed to be zero. For 2-D problems formulated in Cartesian coordinates, system (6) is solved under conditions of plane strain.

The equations of elastic equilibrium (6) are integrated using the pseudo-transient (PT) method (Duretz et al., 2018). In accordance with this method, instead of solving the elliptic problem directly, which usually requires assembling matrices for a discretised system, the pseudo-time derivative is added to the right-hand side of Eq. (6). Then, the system is advanced in pseudo-time from some initial distribution of the parameters until it reaches a steady state.

The advantages of the PT method in comparison to other numerical techniques that require assembling matrices are the brevity of the corresponding code and the simplicity of its development and modification. However, since the PT method involves only explicit time integration, the number of iterations required for convergence scales quadratically with the number of grid blocks due to stability conditions. This is acceptable only for 1-D problems and can be too restrictive for 2-D and 3-D problems. To overcome this limitation and accelerate convergence, we use the second-order pseudo-transient method (Räss et al., 2019), which implies adding both first and second time derivatives to the right-hand side of Eq. (6). The second-order PT method is mathematically equivalent to the method of successive over-relaxation (Frankel, 1950).

Figure 4 Average number of pseudo-transient iterations n_iter required to satisfy different convergence criteria ϵ (a) and relative deviation of the radial component of displacements calculated with method B ($u_r^{\mathrm{B}}$) from those calculated using method A ($u_r^{\mathrm{A}}$) for various grid extension factors (b).

As shown in Fig. 4a, the second-order PT method scales linearly with problem size which is sufficient for the present study. In all benchmarks considered in the subsequent sections, the total computation time of geomechanical code doesn't exceed 10 % of the running time of the hydrodynamic code, if the same hardware is used both for MUFITS and the code implementing method B.

4 Benchmark study

4.1 Problem statement for the hydrodynamic modelling

To validate the developed modelling options, we consider an axisymmetric study of hydrothermal activity at Campi Flegrei, a caldera located in southern Italy. An accelerated ground deformation and heating is observed in this densely populated area of Naples causing interest to understanding and predicting the hydrothermal activity. The corresponding flows in the hydrothermal system and associated observable parameters at the surface were broadly investigated by Chiodini et al. (2003, 2016), Todesco et al. (2003, 2004), Todesco (2009), Rinaldi et al. (2010, 2011), Troiano et al. (2011), and Coco et al. (2016) among others. Thus, the parameters of the hydrothermal activity are well constrained, and they can be used for benchmarking. In this paper, we consider the problem statement that is most consistent with that described by Rinaldi et al. (2011).

We simulate the non-isothermal flow of a CO₂–H₂O binary mixture in the axisymmetric domain $r\in [\mathrm{0},\mathrm{10}]$ km, $z\in [\mathrm{0},\mathrm{1.5}]$ km, where r is the distance to the axis of rotation and z is the depth (Fig. 5). The domain near r=0 corresponds to the region below the Solfatara crater. At the initial moment of time, $t=-\mathrm{4000}$ years, the homogeneous porous medium is saturated with pure water under a hydrostatic distribution of P and a linear distribution of T corresponding to a geothermic gradient of 50 ^∘C km⁻¹. A fixed atmospheric pressure and temperature, P_atm=1 bar and T_atm=20 ^∘C, are imposed at the opened upper boundary $z=z_{\mathrm{top}}=\mathrm{0}$ corresponding to Earth's surface. The system can be recharged with pure H₂O through z=0. The constant T of 95 ^∘C, which is consistent with the geothermic gradient and T_atm, is maintained at the impermeable lower boundary z=1.5 km. The side boundary r=10 km is impermeable and adiabatic. The simulation results indicate that the conditions at r=10 km do not influence the flow near r=0.

Figure 5 The simulation domain and the boundary conditions for the hydrodynamic simulation. The distribution of the gas saturation (s_g) and the isotherms are shown at t=120 months.

The influx of fluid from a deep magmatic source is simulated with a point source placed at r=0 km and z=1.5 km (Fig. 5). A hot mixture of CO₂ and H₂O, which should be regarded as a proxy for magmatic fluid, is injected into the domain through this source. The mixture enthalpy corresponds to T=350 ^∘C at P=180 bar. The pressure near the fluid source increases up to ≈180 bar; thus, the temperature near the source is close to 350 ^∘C. Consequently, it is assumed that the domain includes only a region above the brittle-ductile transition, where T does not exceed 350 ^∘C (Hayba and Ingebritsen, 1997). Therefore, we consider only elastic deformations of the saturated rocks. First, we simulate 4000 years of degassing with a constant injection rate of 3400 t d⁻¹ and a CO₂∕H₂O molar ratio of 0.17 (i.e., the molar concentration of CO₂ is 0.1453). This interval is considered to be a period of dormancy over which a quasi-steady state is reached. After 4000 years, the plume of hot fluid forms near the axis of rotation r=0. It contains two distinct gas zones: one near the injection point and the other at shallow depths. Then, we simulate an event of unrest by temporally increasing the injection rate to 12 100 t d⁻¹ and the CO₂∕H₂O molar ratio to 0.4. The unrest begins at t=0 and lasts over 20 months. After the unrest, at t=20 months, the point source is returned to the quiet state by reducing the injection rate and the molar ratio to their initial values. Other parameters of the study are summarised in Table 1.

Table 1 Parameters of the benchmark study.

For modelling the non-isothermal flow associated with the formation of the plume of hot magmatic fluid, we use a standard system of governing equations. It includes the continuity equations for CO₂ and H₂O, the heat equation, and Darcy's law. The governing equations are described in Afanasyev et al. (2015). They are identical to those used by Chiodini et al. (2003), Todesco (2009), Rinaldi et al. (2011), and others. The only difference is the approach for predicting the vapor-liquid equilibria and the parameters of the fluid phases for the CO₂–H₂O mixture. Chiodini et al. (2003), Todesco (2009), and Rinaldi et al. (2011) applied the TOUGH2/EOS2 simulator (Pruess et al., 1999), which up to the critical point of H₂O employs Henry's law for the solubility of CO₂ in liquid and other empirical correlations. Unlike them, we apply a fully consistent thermodynamic model of the mixture based on a cubic equation of state (EoS). The EoS is used for predicting both the mutual solubilities of CO₂ and H₂O and the phase densities and enthalpies in the whole range of P and T occurring in the study. The details of the EoS and the discussion of its accuracy can be found in Afanasyev et al. (2015).

We use radial non-uniform grids with the r and z factors equal to 1.025 and 1.05, respectively. Thus, the grids become denser near the axis of rotation (r=0) and the surface (z=0) to better resolve the plume. We consider grids 1X (45×8), 2X (60×15), 4X (90×30), etc., where the number of grid blocks along the r and z axes, n_r×n_z, are given in the brackets, respectively. The thickness of the grid blocks is less than 5 m in the case of the most refined grid, 12X (160×90).

For modelling the gravity changes and ground displacement by method A, we specify a network of equally spaced OPs along the straight line z=0 (Fig. 5). The simulator automatically calculates Δg and u for every OP in the network at every moment of time in accordance with the problem symmetry shown in Fig. 2b. Only one OP placed at r=0 and z=0 is needed for reporting the temporal evolution of Δg and u at the axis of rotation (i.e., at Solfatara). All OPs are used for reporting the space distributions of Δg and u at a fixed moment of time. The distributions of P, T, and the bulk density ρ at the beginning of the unrest (at t=0) are chosen as the reference values P₀, T₀, and ρ₀ in Eqs. (1)–(3), (5), and (6).

4.2 Problem statement for the geomechanical modelling (method B)

We apply an extended grid to reduce the influence of the boundary conditions on the elastic equilibrium (Fig. 6). This makes possible the comparison and additional validation of methods A and B. For Eq. (6), we consider the domain $r\in [\mathrm{0},r_{max}]$, $z\in [\mathrm{0},z_{max}]$, where $r_{max}\ge \mathrm{10}$ km and $z_{max}\ge \mathrm{1.5}$ km. The traction-free condition, $\mathit{\sigma }\cdot \mathit{n}=\mathrm{0}$, is imposed at the upper boundary z=0, where n is the normal to the boundary. The “roller” boundary condition, $\mathit{u}\cdot \mathit{n}=\mathrm{0}$, is imposed at the lower and side boundaries, $r=r_{max}$ and $z=z_{max}$. When modelling the ground displacement, we account for the fluid pressure and temperature changes only in the region $r\in [\mathrm{0},\mathrm{10}]$ km, $z\in [\mathrm{0},\mathrm{1.5}]$ km, where the fluid flow is simulated. The changes in P and T are assumed to be zero outside this region.

Figure 6 The domain and the boundary conditions for the geomechanical simulation by method B. The mesh used for the fluid flow modelling is shown in black. The distribution of total displacement is shown at t=20 months.

If the domains for the hydrodynamic and geomechanical modelling coincide, i.e., $r_{max}=\mathrm{10}$ km and $z_{max}=\mathrm{1.5}$ km, then the displacements u calculated by methods A and B differ significantly. This is explained by the influence of the boundary conditions at $r=r_{max}$ and $z=z_{max}$ on u near r=0 because they are not consistent with the assumption of the semi-infinite domain z≥0. If the boundaries are set progressively farther away from the centre region, i.e., $r_{max}>\mathrm{10}$ km and $z_{max}>\mathrm{1.5}$ km, then their influence on u decreases. We investigated optimal values of r_max and z_max that allow for the simulations to be consistent with the assumptions in Sect. 2.4 without excess computational cost. The investigation shows that a three times larger domain along both the r and z axes produces a satisfactory agreement between methods A and B (Fig. 4b). Thus, $r_{max}=\mathrm{30}$ km and $z_{max}=\mathrm{4.5}$ km are used in the following presentation. Herewith, the mesh outside the domain for hydrodynamic modelling should not be very dense. A coarse grid such as that shown in Fig. 6 is sufficient, and further grid refinement doesn't change the displacements in the region of interest.

4.3 Simulation results

The vertical components of the gravity change and ground displacement, Δg_z and u_z, at r=0, z=0 against t are shown in Fig. 7. The positive (or negative) values of u_z correspond to the surface uplift (or subsidence). The positive Δg_z corresponds to the elevated gravity field as compared to the reference state. Over the unrest, the average pressure and temperature as well as the bulk density increase because the hydrothermal system is inflated at the higher injection rate. As a consequence, the surface rises by about u_z=10 cm and the gravity increases by Δg_z=20 µGal. Then, as the fluid source reverts to the quiet state, u_z and Δg_z decrease because the fluid evaporates and is released through the upper boundary z=0, leading to a reduction in pressure and bulk density. The quantity Δg_z reaches a minimum of −120 µGal at t=75 months and then increases back to 0 as the system returns to the quasi-steady state.

Figure 7 Vertical components of the ground displacement (u_z) and gravity change (Δg_z) at the centre OP (r=0, z=0) against time.

An agreement between the u values predicted by methods A and B and those reported by Rinaldi et al. (2011) validates the correctness of the numerical algorithms implemented for modelling the ground displacement (Figs. 7 and 8). In general, Δg_z changes with time similarly to that in Rinaldi et al. (2011). The intervals of Δg_z increase and decrease are identical. However, the absolute values of Δg_z are almost 1.5 times smaller than those reported by Rinaldi et al. (2011). This discrepancy can be caused by different EoS for the CO₂–H₂O mixture implemented in MUFITS and TOUGH2/EOS2. The different approaches for predicting the fluid properties can result in slightly different distributions of the bulk density (ρ) and, thus, the gravity change (Δg).

The results shown in Figs. 7 and 8 are simulated using the 4X grid (see Sect. 4.1). This choice is based on the mesh dependency study, the results of which are presented in Fig. 9. The utilisation of coarse grids, such as 1X and 2X, results in a slight underestimation of both u_r and u_z. The values of u_r and u_z become higher with increasing grid resolution. The convergence is achieved for the 8X grid, and further refinement is not needed. However, u for the 4X grid is satisfactorily close to that for 8X (and 12X) and can be used to reduce the computational cost. Thus, the simulations using method B were performed only for the grid resolutions up to 4X, as shown in Fig. 9c, d.

The results of Rinaldi et al. (2011) are in between our simulation results with the 2X and 4X grids (Fig. 9a, b). Presumably, this is caused by different methods for predicting the fluid properties employed in TOUGH2/EOS2 simulator, which up to the critical point of H₂O employs Henry's law for the solubility of CO₂ in liquid, and MUFITS, which employs a fully consistent thermodynamic model based on a cubic equation of state (EoS). Consequently, the distributions of P, T, and u obtained in this study can be slightly different from those in Rinaldi et al. (2011). However, the distributions are identical in qualitative sense confirming the correctness of the algorithms implemented in MUFITS.

Figure 8 Radial distribution of the vertical (a) and horizontal (b) ground displacement at sequential times. The coloured and dash-dotted lines correspond to methods A and B, respectively.

The distribution of u calculated by method A exhibits an oscillatory behaviour near the maximum of u_r (Fig. 9b). Such an artificial oscillation is caused by the space discretisation of the domain and the finite volume method. When Eqs. (2)–(4) are applied, the changes in P_i and T_i in a grid block, which are distributed continuously over the entire grid block, are pulled to its centre to produce the point source of extension. Therefore, if the OP is closer to such a nearby centre in the uppermost row of grid blocks (i.e. if R₁ and R₂ are smaller), then u is higher than if the OP is between such centres (i.e. if R₁ and R₂ are larger). This results in the oscillatory behaviour because, according to Eqs. (2)–(4), there is inverse relationship between u and R₁ (and R₂). The quantity Δg can exhibit similar behaviour.

Figure 9 Displacements u_z (a, c) and u_r (b, d) at t=20 months calculated using different methods and grids. The oscillatory behaviour of u shown in (b) is caused by the space discretisation.

4.4 Modelling ground displacement in a heterogeneous reservoir

In addition to validating the geomechanical code against the semi-analytical approach, we consider another benchmark study to demonstrate the capabilities of method B for modelling hydrothermal systems with spatial heterogeneities in elastic properties of rocks. The study is identical to that described in Sect. 4.2 except for the values of μ and K. Now, the distributions of μ and K are taken to qualitatively reflect the structure of the shallow subsurface at Campi Flegrei as it is presented by Manconi et al. (2010) based on the seismic tomography study by Chiarabba and Moretti (2006). Shear modulus μ shown in Fig. 10 varies in the range of 2–14 GPa, and Poisson's ratio ν lies in the range of 0.2–0.4. The general trend is that μ rises exponentially with depth. Closer to the axis of symmetry, i.e., inside the caldera, the material is softer. A thoroidally shaped inclusion composed of much stiffer material, which models the caldera rim, is located approximately 5.5 km from the axis. Assuming that the stiffer material is also less compressible, the Poisson's ratio ν is linearly interpolated between 0.2 and 0.4 to match the minimal and maximal values of μ, respectively. The drained bulk modulus K is then calculated according to the relation

$$K={\displaystyle \frac{\mathrm{2}\mathit{\mu }(\mathrm{1}+\mathit{\nu })}{\mathrm{3}(\mathrm{1}-\mathrm{2}\mathit{\nu })}}.$$

The distributions of u_r and u_z at z=0 calculated with the described distribution of mechanical properties do not differ much from those calculated under the assumption of a homogeneous medium (Fig. 10). The deviation from the results presented in Sect. 4.3 does not exceed 20 %. This can be explained by the synthetic nature of the benchmark, in that K is artificially constrained to enforce consistency with the range of ν without direct relation to the seismic data. Furthermore, the results reported by Rinaldi et al. (2011) and other authors were thoroughly calibrated against observations, so the effect of incorporating heterogeneities into the numerical model is not expected to be significant in this particular case.

Figure 10 Distributions of u_r and u_z at z=0 km and t=20 months (a) and shear modulus (b).

5 Summary

The developed extension of MUFITS allows for the convenient built-in calculations of ground displacement and gravity changes. These calculations are performed automatically by the simulator without the involvement of any external post-processing utilities. As a consequence, the reservoir models developed and simulated with MUFITS can now be history matched to the observations of gravity changes and ground displacement. This software development, although quite straightforward, makes MUFITS closer to a universal package for modelling flows in hydrothermal systems. The developed modelling options, applied here in a study of hydrothermal activity, can also be utilised in other applications, including oil and gas extraction (Fig. 1), by employing a different fluid property module of the simulator (Afanasyev, 2015).

The simulation results also demonstrate an acceptable accuracy of the semi-empirical approach (method A) for predicting ground displacement. Given that the necessary assumptions in Sect. 2.4 are satisfied, the displacements predicted by utilising analytical solution given by Eqs. (3) and (4) are indistinguishable from those predicted by the more comprehensive approach (method B). This makes method A very attractive because it is fast and can easily deal with 3-D flows. However, if the assumptions in Sect. 2.4 are violated, then the coupling with the geomechanical code remains the only reliable approach for the ground displacement modelling. For example, the thermo-hydromechanical coupling is the only approach to assess the mechanical integrity of reservoirs consisting of the layers with significantly different poroelastic properties (e.g., sandstone layers alternating with shales) or complicated by the presence of faults. To assess the reservoir integrity, one needs to know the distribution of stresses in the entire domain, which is readily available with method B. The influence of topography, e.g., a load of volcanic edifice, on the subsurface flows and the ground displacement can be estimated only by method B. Besides the space heterogeneity, the semi-empirical method, unlike method B, does not allow for generalization to inelastic behaviour of saturated rocks, e.g., at high temperatures (>350 ^∘C) that can be achieved in hydrothermal systems. The extension of MUFITS onto these cases remains a subject of further code development.