THM modeling of gravity anomalies related to deep hydrothermal circulation at Soultz-sous-Forêts (France)

Gravity measurements in the Upper Rhine Graben evidence spatial variations at the regional scale and close to the geothermal sites. They are classically interpreted as linked to the geology. We aim to bring new insights on another potential origin of these gravity changes. Our approach is to quantify gravity anomalies related to the deep hydrothermal circulation. A thermo-hydro-mechanical model is developed at the reservoir scale for the Soultz-sous-Forêts site (Soultz), France. A finite element method is used in 2D and 3D. The size of the representative elementary volume is assumed to be 100 m with no regional fault in the reservoir. Surface gravity profiles and maps associated to the large-scale hydrothermal circulation are computed from the variations of the effective density through the model. Synthetic spatial gravity variations in 2D and 3D models are shown to have an amplitude of 0.02 mGal. They are shown to be mostly linked to the convective system. Their wavelength is about 7.5 km, consistent with the width of the hydrothermal convection cells. The anomaly maximum is located at the top of the maximum surface heat flux. However, gravity anomaly observations show much higher amplitude and heterogeneity. Spatial gravity variations linked to the hydrothermal circulation are shown to be smaller than the observed gravity spatial variations, but still measurable with very sensitive instruments (absolute gravimeters).

However, other contributions such as the hydrothermal circulation have been investigated to explain the origin of the spatial gravity changes in the URG. These recent studies include hybrid gravimetric surveys (Portier et al. 2018) or reinterpretations of the current potential method database (Baillieux et al. 2014;Abdelfettah et al. 2018).  Rotstein et al. (2006) Baillieux et al. (2013) highlights the potential relationship between the gravity lows and positive temperature anomalies due to the hydrothermal circulation. Butterworth filtering of the Bouguer anomalies at different spatial wavelength has been used in order to cancel the regional trend of the gravity signature ). The filtering shows fault parallel negative anomalies linked to regional faults identified as fluid pathways by geochemical tracers and MT measurements.
The exact origin of the spatial gravimetric variations measured in the URG is still unclear. The aim of the current study is to quantify the contributions from hydrothermal circulation in the observed gravity anomalies. To address the issue, we propose here a thermo-hydro-mechanical (THM) modeling of one deep geothermal reservoir in the URG (Magnenet et al. 2014;Vallier et al. 2018Vallier et al. , 2019. We focus on the spatial gravity variations observed at the deep geothermal site called Soultz-sous-Forêts (Soultz) located in the Northern Alsace (France). The Enhanced Geothermal System (EGS) project is associated to an available and important database allowing a precise characterization of the geothermal reservoir (Bresee 1992;Kappelmeyer et al. 1991;Gérard and Kappelmeyer 1987;Olasolo et al. 2016;Schaming et al. 2016;GeORG 2013).
For the first part of the study, a comparison has been conducted between the 2D and 3D simulations of the hydrothermal circulation for the Soultz reservoir. The 3D model is following the same parametrization as in the 2D approach. The goal of the comparison is to show the consistency between the two approaches. For the second part, we quantify the gravity changes from density variations due to the poro-elastic effects or to the hydrothermal circulation. Based on a laterally homogeneous and simplified geology, the simulated gravimetric signal from 2D and 3D simplified reservoir models are compared with the gravity observations in the URG and we finally discuss the exact contribution of the hydrothermal circulation sensitivity on the gravity signal.  Vallier et al. Geotherm Energy (2020) 8:13 Overview of the gravimetric data Gravity anomalies have been measured in the URG since 1930 and a large network of gravity stations is currently settled in the region (Rotstein et al. 2006). The existent database includes: 12,250 gravity stations from GGA Leibniz Institute in German national gravity archive; over 12,400 stations from BRGM database; 11,000 stations from MDPA (Potash Mines of Alsace); 100 stations on salt domes in the south of the URG (Baillieux et al. 2014;Rotstein et al. 2006). Spatial variations of the gravity have been observed as illustrated in Fig. 1 for the whole Graben and in Fig. 2 for a close-up around the Soultz-Rittershoffen area. In Fig. 2, the Bouguer anomaly map is compared to the spatial distribution of the isotherms and the regional fault system in the URG. The observed Bouguer anomaly is negative with lateral variations between − 39 and − 23 mGal in the area of Soultz and Rittershoffen. The highest Bouguer local anomalies are located near the two geothermal sites, Soultz. At larger scale, spatial variations of gravity are also observed (see Fig. 1). The amplitude of the measured gravity increases globally with the distance from the center of the Graben with variations from − 90 to 20 mGal.
The gravimetric studies have classically interpreted the spatial variations of gravity in the URG as mainly due to the lithology and the geological structure of the region. In Rotstein et al. (2006), the gravity measurements are compared with the URG geomorphology. The comparison shows an accordance between the gravity spatial distribution and the shape of the continental rift. Baillieux et al. (2014) and Rotstein et al. (2006) both interpret some specific gravity lineament to be linked with main regional faults in the URG. Edel et al. (1982) and Rotstein et al. (2006) relate the gravity signal to the variation of composition in the granitic basement. For example, the positive anomaly near Karlsruhe is associated to Paleozoic schists included in the basement having a higher density than the other rocks (Rotstein et al. 2006).
All studies assumed that the contribution of lithology is the prominent one to explain the spatial variations of the gravity in the URG. However, recent studies highlight a potential relationship between the negative gravity and positive temperature anomalies due to the hydrothermal circulation (Baillieux et al. 2013). The largest temperature anomaly at Soultz corresponds to the largest Bouguer anomaly. Time-lapse relative microgravimetry and absolute gravimetry measurements have been performed since 2014 in order to study the mass redistribution in the Soultz and Rittershoffen reservoirs during their industrial exploitations (Hinderer et al. 2015;Hector and Hinderer 2016;Portier et al. 2018). This hybrid gravimetry monitoring reveals a temporal and spatial correlation between the fluid exploitation and the gravity change of the order of − 31.0 ± 9.0 µGal near the injection area. Our modeling approach aims to quantify what are the contributions differing from the lithology to the spatial variations of the gravity anomalies in the URG.

Simplified large-scale representation of the reservoir
To address this issue, gravity anomalies are simulated here using a simplified geological model without taking into account regional faults. We assume that the deep hydrothermal circulation in the upper crust is mainly influenced by the closely connected small-scale fracture network and not by the large-scale sparse faults (Vallier et al. 2019). The validation of this main hypothesis in our modeling has been discussed and supported in Vallier et al. (2018). The recent work of Freymark et al. (2019) highlights a weaker role than expected of the major border faults on the hydrothermal circulation at the Central URG scale.
The simplified representation of the geothermal reservoir is shown for the Soultz case in Fig. 3b and compared to a more extended geological model in Fig. 3a. The underground is divided into four homogenized units: the upper and lower sediments; upper and lower granites. The depth of the transition between the upper and lower sediments is estimated at 0.1 km after the back-analysis of the observed temperature and stressdepth profiles (Vallier et al. 2019). The depths of the sediments-granite interface and the transition between the two granites are, respectively, at 1.4 km and 3.9 km in depth.

Governing equations of the THM model
We described the thermo-hydro-mechanical (THM) coupling assuming that the four units are homogenized as a porous medium, fully saturated with a single-phase brine in the thermo-elastic regime. The equations governing the THM coupling are developed from the reference book of Coussy (2004) and detailed in Appendix 1. Here, several assumptions are made: (i) the Cauchy stress tensor σ is composed of an effective stress σ ′ and a hydraulic stress σ p 1 ( 1 being the unit tensor); (ii) the thermodynamic flows are linearly related to thermodynamic forces (linearized strain ǫ , gradient of pore pressure a b Fig. 3 a 2D conceptual geology at Soultz [modified after Dezayes et al. (2005a, Fig. 28, 2005b and after Aichholzer et al. (2016, Fig. 3)]. Trajectories of GPK1, EPS1, GPK2, GPK3 and GPK4 wells are shown with colored lines. b Temperature-depth profiles after drilling operations superimposed to the simplified reservoir model with four main geological units (Cuenot et al. 2008). Here, e 1 , e 2 , e 3 and e 4 correspond to the unit thicknesses, e 1 is inverted from a back-analysis (Vallier et al. 2019) using: e 2 = 1.4 km − e 1 ; e 1 + e 2 + e 3 = 3.9 km ; e 1 + e 2 + e 3 + e 4 = 5.4 km ∇p w , gradient of temperature ∇T ). Most of the homogenized properties in Hooke's law, Darcy's law, and Fourier's law depend on porosity, fluid pressure and temperature; (iii) the brine rheology is extrapolated from experimental results for artificial brines at different salinities (Zaytsev and Aseyev 1992;Kestin et al. 1981;Rowe and Chou 1970). The fluid is assumed to be a pure NaCl solution with a mean specific mass content of 100 g L −1 (Magnenet et al. 2014;Genter et al. 2010).

Numerical method to simulate the hydrothermal circulation in 2D and 3D approaches
To solve the constitutive equations, the open-source finite element software Code_ Aster is used (EDF 2014). Some particular developments have been supplemented in the solver concerning the radioactive sources, the brine rheology and the search of stationary solutions (Magnenet et al. 2014). An Euler implicit scheme is used to vanish the nodal residues and the Newton-Raphson method is used to solve the system of nodal unknowns (Vallier et al. 2019). The hydrothermal circulation at Soultz has been simulated after inversion of the temperature and stress-depth profiles in a 2D approach by using the PEST software (Doherty 2005). We checked that the parametrization of the 2D model is similar to the 3D extension of the model. The size of the finite elements is 100 m × 100 m , i.e., the size of the supposed representative elementary volume. The 2D model has been considered as a vertical cross-section of the reservoir (see Fig. 4a). Its height and width are, respectively, 5.35 km and 30 km. The 3D model is an extension of the 2D representation along the out-of-plane horizontal direction (see  The bottom is set at 5.35 km in depth, but the two horizontal dimensions are both set to 10 km for the 3D model and 30 km in the 2D approach.
The THM boundary conditions are also illustrated in Fig. 4 for the 2D and 3D approaches: • The temperatures are set on the upper and lower boundaries. The lateral boundaries are taken as adiabatic. • The fluid pressure is set at the value of the atmospheric pressure on the upper boundary. The other boundaries are assumed to be impermeable. • No normal displacement is considered on the lower and lateral boundaries and the upper boundary is stress free.
Uniform temperature and fluid pressure fields are initially assumed for the whole geothermal reservoir. Three consecutive steps in the calculation are considered to facilitate the convergence of the calculation (Magnenet et al. 2014): • during 1000 years, the boundary conditions and gravity are progressively applied.
• during 100,000 years, the system freely evolves along constant boundary conditions. • in one last increment, the system is set at a stationary state by canceling the nonstationary terms from the THM system of equations.

Forward simulation of the gravity anomaly
By solving the THM problem in 2D and 3D, hydrothermal circulation has been simulated at the reservoir scale. The circulation induces a spatial change of the total homogenized specific mass which is equivalent to a spatial variation of the density of the total homogenized medium. The density of the total homogenized medium depends on the porosity and is linked to the fluid and solid grains densities via a classical mixing law: where "f " and "s", respectively, refer to the fluid and solid grains, φ is the porosity and ρ , the specific mass. The gravity signal is deduced from the distribution of density variations in the reservoir. The gravity effect is computed from the point mass anomaly formulation (Torge 1989). We denote by x the horizontal direction in the 2D model, z the out-of-plane horizontal direction in 3D and y the vertical one. For each node of the elements along the x-axis on the surface, the gravity effect of every Gauss interpolation point of the finite elements is calculated as follows: with N, the total number of interpolation points in the spatial domain; ω i , the weight of the i th interpolation point; dρ i , the difference between the density of the total homogenized medium extracted at the i th point and the reference value ρ ref ; r i , the distance between the point on the surface and the i th point and G , the gravitational constant.
Here, the reference density of the total homogenized medium ρ ref is taken as the initial which is shown to be similar to the mean value in the spatial domain ρ m . The difference of gravity effect obtained by considering ρ 0 or ρ m is discussed in details in Appendix 2.

Comparison between the 2D and 3D models of hydrothermal circulation
The simulated circulation has been obtained after inversion of the observed temperature and stress-depth profiles at Soultz and Rittershoffen site (Vallier et al. 2018(Vallier et al. , 2019. Both studies are limited to a 2D representation of the deep geothermal reservoir as a vertical cross-section. Here, we focus on the case of the Soultz geothermal reservoir. The parametrization for the 2D model inverted from the observed profiles (temperature, stress) is introduced in a 3D approach to check the representativeness of the 2D results. The simulated hydrothermal circulation from the 3D model is then compared to the results of the 2D approach.  Figure 5a, b illustrates the temperature spatial distributions from the 2D and 3D models, respectively. The spatial temperature distribution from the 3D approach consists of a vertical cross-section along the x-axis in the middle of the reservoir. The temperature maps from 2D and 3D approaches show strong similarities. The hydrothermal circulation occurs in the upper granitic reservoir as well as in most of the sediments. The hydraulic cap-rock is at about 100 m in depth for both approaches. As in the 2D model, four convection cells are developing in the 10 km of the reservoir along the x-axis as well as the y-axis. Their heights and widths are, respectively, about 3.8 km and 3.0 km as in the 2D approach. Figure 5c pictures the comparison of the temperature-depth profiles between the 2D and 3D models. These temperature-depth profiles have been extracted at the location of the surface heat flux, i.e., the middle of the two vertical cross-sections. The temperature evolution with depth in the 3D model is very well reproduced compared to the one from the 2D approach in particular until 3.0 km deep. The temperatures from the 3D model are slightly lower than those of the 2D approach in the lower granitic basement but the discrepancy is never more than 5 °C.
To sum up, the 3D model reproduces a similar hydrothermal circulation compared to the 2D model when using the same parametrization. The same number of convection cells exists in the reservoir with the same dimensions for both approaches. The evolution of the temperature with the depth is nearly identical at the location of the maximum surface heat flux between the 2D and 3D models. Thus, the similarities between 2D and 3D models validate that the 2D approach is adequate to describe the hydrothermal circulation for a simplified large-scale representation.

Profiles of the synthetic gravimetric signal
The gravimetric signal associated to the large-scale hydrothermal circulation is calculated from the two-dimensional model considering the Soultz reservoir. The study does not aim at reproducing the average value of the observed Bouguer anomaly in the Soultz area. Indeed, our reservoir model is limited at the first 5 km in depth. Much deeper structures in the crystalline crust and the lithospheric mantle until 80 km deep are known to highly contribute to the gravity anomalies observed at the surface in the Upper Rhine Graben (Freymark et al. 2017;Pavlis et al. 2012). Moreover, our idealized reservoir representation does not include any lateral structural variation for the Soultz reservoir. As developed in the overview of the gravimetric data, the lithology spatial changes are currently interpreted as the main contribution of the gravity regional spatial distribution (Rotstein et al. 2006;Edel et al. 1982). The current study aims at estimating the potential spatial variations linked to the large-scale hydrothermal circulation and at characterizing their amplitude and wavelength. To address the issue, we consider a 2D model of 30 km width and 5.35 km height and study the calculated gravity signal in the central part, between 5 and 25 km to discard the boundary effects. The signal is then compared to the associated hydrothermal characteristics of the Soultz reservoir. Figure 6a features the simulated relative variation of the total homogenized specific mass, i.e., the density of the total homogenized medium associated to the hydrothermal circulation at the Soultz reservoir. The spatial variations of the density of the total homogenized medium are included between − 0.125 and 0.125%. A density contrast is located at the geological transition between the sediments and the granites. However, the spatial distribution of the density changes is independent of the lateral geological changes. Indeed, lateral density variations are obtained, whereas no lateral lithology change is included in the reservoir model. Figure 6b shows the associated profile of the simulated gravimetric signal at the surface (along the x-axis). Spatial variations in the simulated gravity signal exhibit oscillations around the mean value of − 2.27 mGal. The amplitude of the oscillations is about 0.02 mGal (i.e., 20.0 µGal ). The wavelength is 7.5 km, equal to the convection cell horizontal extension. The largest gravity anomaly is located at the middle of the profile, at 15 km. Figure 6c shows the vertical component of the computed surface heat flow along the x-axis associated to the hydrothermal circulation. The heat flow is about 150 mW m −2 ± 85 mW m −2 . These values are consistent with the estimated ones at Soultz and from the Upper Rhine Graben region (Bachler et al. 2003;Pribnow et al. 1999;Pribnow and Schellschmidt 2000;Clauser et al. 2002). Gravity anomalies are shown to be opposed in phase to heat flux anomalies.
To sum up, lateral gravity variations not linked to lateral lithology changes but to the large-scale hydrothermal circulation are quantified. They are very weak, of the order of 20 µGal and having an anti-phase relationship with the heat flux anomalies that are relatively much larger ( 85 mW m −2 ). a b c Fig. 6 a Map of the simulated relative variation of the total homogenized specific mass, i.e., the density of the homogenized medium associated to the hydrothermal circulation at the Soultz reservoir. b Associated synthetic gravimetric profile along the x-axis. c Profile of the computed vertical surface heat flow along the x-axis. The dotted line corresponds to the mean value of the heat flow

Map of the gravimetric signal
Observed gravity maps of the URG are available and are compared to the results of our simulations (Baillieux et al. 2014;Rotstein et al. 2006). The same methodology is applied for the three-dimensional approach. The depth of the Soultz reservoir is set at 5.35 km as in the two-dimensional approach. The two horizontal dimensions are both set at 10 km instead of 30 km due to the cost of CPU time. Based on Equation 2, but extended to a 3D context, synthetic maps of gravity signals are computed from the THM model. Gravity profiles are extracted from the 3D model to be compared to the ones from the 2D approach. Comparisons with heat flow map and profiles are also done. Figure 7a shows the map of the simulated gravimetric signal obtained from the same parametrization as in the 2D model. The variations of the gravimetric signal are between around − 2.28 mGal and − 2.26 mGal. The maximum of the gravity anomaly is located at the center of the domain. Figure 7b illustrates the spatial distribution of the calculated vertical component of the total heat flux extracted from the surface. The heat flow varies similarly to the 2D simulation, from 235 mW m −2 at the middle of the reservoir to less than 50 mW m −2 at the descending part of the convection cells (see Fig. 5b).
The comparison of the gravity signal between the 2D and 3D approaches is shown in Fig. 7c. The gravity profile from the 3D model is extracted at the middle of the y-axis along the x-axis. The two gravity signals show many consistencies between each other in terms of variations and average value but some discrepancies are noticeable. Both signals in 2D and 3D present spatial variations between − 2.28 and − 2.26 mGal. The wavelength is consistent but influenced by larger finite size effects since the model has

Contribution of the hydrothermal circulation compared to a fully thermal diffusive case
To quantify more precisely the sensitivity of the gravity signal to the hydrothermal circulation, we compare the results to a configuration where we cancel the convective heat transfer in the model. To obtain a fully thermal diffusive case, the homogenized permeability is set to 10 −20 m 2 for the sediments and the granitic basement. Figure 8 features the comparison of the gravimetric profiles between the fully thermal diffusive case (as orange line) and the one associated to the large-scale hydrothermal convection (as blue line) for the 2D Soultz model. For the fully thermal diffusive case, there is no variation of the gravity along the x-axis, contrary to the case where the large-scale hydrothermal convection occurs in the reservoir. The simulated gravity main value is negative and about − 2.275 mGal as estimated with the case associated to the large-scale convective system. This comparison confirms that the gravity spatial variations are related to the large-scale convective system. The sensitivity of the gravity signal to the hydrothermal circulation can be clearly quantified as oscillations around the mean value of about 0.02 mGal.

Comparison with observed gravity anomalies
Similarities are highlighted between the observed gravity signal and the simulated one: (i) the simulated gravity signal in 2D and 3D (see Figs. 6 and 7) is negative like the observations (see Figs. 1 and 2) with a maximum value closely located at the maximum of the thermal heat flow (Rotstein et al. 2006;Baillieux et al. 2014;Abdelfettah et al. 2012); (ii) both observed and calculated gravity signals shows spatial variations with some consistencies with the hydrothermal circulation (Baillieux et al. 2014). However, the range of values between the calculated and observed gravity anomalies oscillations is clearly different: [− 39 ; − 35] mGal for the observations (Baillieux et al. 2014) and [− 2.28 ; − 2.26] mGal for the calculated signal. We aim to further compare the amplitude and wavelength of the gravity signal between the simulated and the observed ones. To address the issue, the average value of the simulated gravity signal has been artificially shifted in order to match the expected range of values from the observations (Baillieux et al. 2014;Abdelfettah et al. 2012) assuming that deep contributions from the lower crust are dominant but not described in our model. The shifted profile is then compared to a profile oriented north-south extracted from the gravity map of Baillieux et al. (2014) (see Fig. 2). Figure 9 illustrates the comparison between the artificially shifted simulated profile and the observed one oriented north-south. The shift applied to the calculated signal is about − 34.7 mGal to obtain a match with observations. The observed amplitude is about 1.5 mGal, a much higher value than for the simulated one about 0.02 mGal. The observed gravity map (cf Fig. 2) and the calculated one (cf Fig. 7) show different spatial distributions. Both 2D and 3D approaches do not reproduce either the amplitude or the complete spatial distribution of the observed signal. The higher amplitude and spatial heterogeneity in the observations may be due to a contribution of the lithology changes at regional scale (Rotstein et al. 2006;Baillieux et al. 2014). The main origin of the gravity anomaly is most certainly not located in the first 5 km in depth. Geological structures in the deeper crust and the lithospheric mantle until 80 km in depth greatly influence the gravimetric signal in the Upper Rhine Graben (Freymark et al. 2017;Pavlis et al. 2012). The necessary shift of − 34.7 mGal probably is associated to the lack of this main deep contributions in the simulated gravity signal. Other contributions which are not included in our current model can also potentially explain the discrepancy between the observed and simulated profiles. Baillieux et al. (2014) investigated the potential signature of fluid pathways through major fault systems on the observed gravity. Further studies still need to be yet carried out to build a 3D THM model that includes the major faults to fully interpret the observed signal. Geochemical reactions between geothermal waters and the sediments have also to be considered since they have been shown to influence the Bouguer anomalies in hydrothermal systems (Goldstein and Carle 1986;Allis 1990).

Conclusion
The comparison between the 2D and 3D THM models of the deep geothermal reservoir at Soultz-sous-Forêts or Rittershoffen reveals similar large-scale hydrothermal circulation. The same number and dimension of convection cells have been obtained in both models as well as very close temperature-depth profiles at the ascending part of the convection cells. This study supports that the 2D approach is sufficient to simulate the largescale hydrothermal circulation. The second part of the study deals with a forward THM modeling to quantify the contribution from hydrothermal circulation in the observed gravity anomalies. The amplitude of these oscillations is about 0.02 mGal (i.e., 20.0 µGal ) around a average value of − 2.27 mGal. Their wavelength is about 7.5 km, consistent with the width of the convection cells. The oscillations in the gravity signal show an anti-phase relationship with the surface maximum heat flux. The 2D and 3D approaches are very consistent. The amplitude of the gravity signal is slightly lower and the wavelength higher in the 2D model, but it includes larger finite size effects.
A comparison with the fully thermal diffusive case clearly supports that the gravity oscillations of 20.0 µGal in amplitude are only due to the large-scale convective system in the reservoir. The simulated signal could not reproduce the important amplitude nor the heterogeneity of the gravity spatial distribution observed in the URG. This suggests that lateral geological changes or close by thermal anomalies or structures in the deeper crust and mantle do have an higher influence on the gravity signal than the hydrothermal circulation.
with K 0 the bulk modulus of the skeleton and c s the specific heat at constant stress. Note that ρ w and µ w as well as and c s are functions of temperature and/or pore pressure. For the Fourier's law, the thermal conductivity of the dry rock is described by the classical mixing law: with s (resp. air ) the thermal conductivity of solid grains (resp. air). We assume that the thermal conductivity of air is negligible. Consequently, the thermal conductivity of solid grains can be written as: Thermal conductivity of the dry material is assumed to depend linearly on temperature: with a dry and b dry empirical constant parameters obtained from experimental measurements. Finally, the homogenized thermal conductivity of the saturated porous media is expressed by using the same kind of mixing law as previously: The specific heat for the dry medium is defined using a similar experimental correlation as Eq. (16): As proposed for the homogenized thermal conductivity, we can define the specific heat capacity and the initial specific mass as: with c air the specific heat capacity of air.

Validation by the V-sheet analytical case
We reproduce the gravity effects calculated from robust analytical expressions. The first analytical case concerns the gravity anomaly due to a buried vertical finite sheet body (Hinze et al. 2013). Figure 10 illustrates the conceptual representation of the geometric parameters for modeling the gravity effect of the V-sheet. The analytical equation of the gravity effect noted g vs from the V-sheet source is as following (Hinze et al. 2013): where g vs is the vertical gravity anomaly in Gal ( m s −2 = 10 5 mGal ) calculated at each point along x-axis on the surface due to the density contrast dρ ( kg m −3 ) at distances x (m) and z (m). The other parameters t, z 2 1 and z 2 2 are illustrated in Fig. 10.
For the analytical cases considered here, the heat transfer is assumed to be only driven by diffusion. All THM couplings are canceled including the temperature and fluid pressure dependencies for the fluid and rock properties. The reservoir is considered as a 2D cross-section of 10 km in width and 5.35 km in height and the V-sheet is 200 m in width and 1000 m in height. The sheet is set to have a specific mass contrast of 300 kg m −3 with the rest of the reservoir. Figure 11 shows the comparison of the gravity effect profiles along the surface between the analytical expression and the simulation for the V-sheet case. First the dependence to element types is investigated. The simulated gravity anomaly is obtained from simulations with two different kinds of finite elements: triangle and quad. No significant effect of the element type is shown. The issue concerning the reference specific mass is addressed for the V-sheet analytical case. For the modeling with quad elements, two cases are considered: (i) the reference total homogenized specific mass ρ ref is taken as the initial value ρ 0 ; (ii) ρ ref corresponds to the mean value in the spatial domain ρ m . The two simulations reproduce the analytical solution, demonstrating that the assumption " ρ ref ≈ ρ 0 ≈ ρ m " is validated for low gravity effect in large-scale reservoir. Finally, all (21) g vs = 2Gtdρ ln z 2 2 + x 2 z 2 1 + x 2 ,

Validation by the Mogi's sphere analytical case
The Mogi's sphere has been first introduced in the works of Mogi (1958) to study the link between the ground deformations and volcanic eruption and can be applied in a geothermal context (Heimlich 2017;Portier et al. 2018). Figure 12 illustrates the conceptual model of the buried Mogi's sphere. The analytical equation of the gravity effect noted g m is (Hinze et al. 2013): where the same notations and the geometric parameters are shown in Fig. 12.
The hypotheses and reservoir geometry for the Mogi's sphere model are the same than in the V-sheet case. The radius of the sphere is 500 m and its center is at 2.65 km in depth. It uses a density contrast of 300 kg m −3 . Figure 13 illustrates the comparison between the analytical solution and the simulated gravity effect. The simulated gravity anomaly profile is well reproducing again the analytical solution. Our modeling (22) g m = 4π G dρR 3 3z 2 (1 + ( x z ) 2 )