UNCERTAINTY ANALYSIS OF GRAVITY DISTURBANCE FOR GEOTHERMAL EXPLORATION IN THE GENEVA BASIN (SWITZERLAND)

Gravity data from the International Gravimetric Bureau and the Gravimetric Atlas of Switzerland have been used to evaluate their application and limitations as a subsurface investigation tool to constrain key geological structures in support of the geothermal exploration in the Geneva Basin (GB). In this context, the application of an effective processing workflow able to produce a topography-free gravity disturbance and quantify its uncertainty is a crucial first step to delineate gravitational effects correlated to geologic sources of geothermal interest. This study focuses on an innovative processing workflow applied to publicly available gravity data in order to produce stochastic realizations of topography free gravity disturbance residuals anomalies. The resulting realizations are compared to a standard interpolation technique and correlate to the available density data and geologic knowledge, demonstrating the potential of the gravity method for geothermal exploration in sedimentary basins.


Introduction
The deployment of renewable energy sources for both power and heat production is accelerating in Switzerland, a trend that will continue thanks to the 2050 Swiss Energy Strategy (SFOE, 2018) that aims at gradually phasing out nuclear power by reducing the energy consumption and increasing heat and electric power generation from renewable energy sources. Geothermal energy will be, therefore, an important resource to supply heat and power for industrial, agricultural, and domestic use.
Increased energy demand, together with the political vision of reducing the use of fossil fuels for heat production in the Canton of Geneva, triggered the development of medium to long-term activities under the umbrella of the GEothermies program (Moscariello, 2016(Moscariello, , 2019. This program aims to identify potential geological targets at shallow/medium (0.5km to 3 km) to large depth (>3 km) to combine heat, power production and potentially mineral extraction.
The Geneva Basin (GB) has been intensively explored for hydrocarbon exploration since the 1960s and for geothermal exploitation in the 1990s, but only non-economically viable production of geo-resources has been recorded. The Thônex-01 well (Jenny et al., 1995), drilled for geothermal heat production, was not commercially productive despite the favourable bottom hole temperature of 88°C at 2530 m TVD. However, the geothermal well GEo-01, drilled in 2018, proved to be successful with a discharge of 50 l/s of geothermal water at 34°C from the upper Mesozoic units (i.e., Lower Cretaceous and Upper Jurassic) and an 8 bars wellhead pressure (SIG, 2018;Guglielmetti et al., 2020;Moscariello et al., 2020).
The geothermal conditions in the Geneva area have been reconstructed by thermal modeling (Chelle-Michou et al., 2017) and geochemical data  revealing a geothermal gradient in the area ranging from 25°C to 30°C. Such a gradient predicts a temperature up to 150°C at the top of the basement in the southern part of the Geneva area. Technically and economically extractable geothermal resources are found at many depths in the Geneva area. At only few tens of meters in the Quaternary deposits, the potential for low enthalpy ground-source heat pump installations is already broadly exploited. From a depth of 0.5 km to 3 km the porous Cenozoic Molasse and fractured Mesozoic sequence offers possibilities for both heat extraction and storage (GeoMol Team, 2015) . At 4 km to 5 km the Triassic carbonates, the Permo-Carboniferous (PC) sediments, and the crystalline basement have potential for cogeneration of power, heat, and metal extraction.
The identification and characterization of these geological structures before drilling are crucial to target potential geothermal reservoirs. The geophysical method that has demonstrated the best results in many geothermal contexts is 2D seismic. This is true both in sedimentary basins such as the Bavarian region (Lüschen et al., 2014), Eastern Switzerland (Heuberger et al., 2016), and also in volcanic regions such as Larderello (Casini et al., 2010). Seismic 2D was mostly acquired in the Geneva Basin (GB) in the 1980s to the 1990s for hydrocarbon exploration (Clerc et al., 2015;Moscariello, 2019) but shows some limitations in delineating shallow Quaternary deposits and deep geologic structures such as the Permo-Carboniferous (PC) grabens which can be associated with lateral density contrast with respect to the surrounding geologic formations. Gravity can therefore contribute to reducing such limitations, this being a standard geophysical subsurface exploration technique in different geological settings (Blakely, 1995;Reynolds, 2011;Telford et al., 1990), commonly applied in the early stages of geothermal exploration programs to identify regions of potential interest (Allis et al., 2000;Uwiduhaye et al., 2018;Guglielmetti et al., 2013) and for monitoring production operations (Portier et al., 2018).
In the framework of a continued desire to improve the understanding of the subsurface in the GB and reduce the uncertainty related to its petrophysical properties, the aim of this study is to assess the value of existing gravity data. This is achieved through an innovative processing workflow to obtain reliable stochastic realizations of gravity disturbance within the geothermal horizon targets, which allow quantifying gravity disturbance uncertainty and to better correlate gravitational effects to geologic sources of geothermal interest.

Geological setting
The GB is the westernmost part of the North Alpine Foreland Basin that extends from the Savoy region in France to Linz in the Austrian area. The GB covers about 2000 km 2 from the town of Nyon to the NE, down to Vuache Mountain to the SW and it is limited by the Jura Haute-Chaine to the NW and by the subalpine nappes to the SE (Fig. 1), and (Kempf and Pfiffner, 2004;Kuhlemann and Kempf, 2002;Mazurek et al., 2006). The basin originated as a lithospheric flexure response to the topographic load of the Alpine orogeny (Pfiffner et al., 2002) which crustal structure in the study area can be described as a gently dipping surface towards the Alpine orogen with increasing thickness from an average of 20km to 30 km in the Molasse Basin to 40km to 50 km in the high topographic region of the Alpine orogen and to 50km to 60 km south of the Pennidic front.
In the GB four major lithostratigraphic units are recorded at depth (Jenny et al., 1995;Rusillon, 2018;Brentini, 2018;Moscariello, 2019). From bottom to top, these are 1) the crystalline basement including PC troughs at the bottom and its sedimentary cover composed respectively of 2) Mesozoic carbonate units and, at the top, 3) Cenozoic and 4) Quaternary sediments. These units can be approximated to a so-called layer-cake model with sub-parallel formations gently plunging towards SE with an average dip of 3 degrees (Fig. 1).
The crystalline basement is only exposed in the Alps, to the South of the GB, and has been drilled in the study area only by a limited number of wells Rusillon (2018). In the GB and, more generally, across the entire Molasse Basin, the basement is often affected by SW-NE oriented tectonic depressions originated as pull-apart basins during Carboniferous (McCann et al., 2006). These pull-apart structures were further developed during Permian in a wrench faulting regime with syntectonic sedimentation filling them with several thousand meters of deposits (P. Diebold, 1994;Ziegler, 1990). In the Geneva Basin, Permo-Carboniferous (P-C) grabens have been mapped using reflection seismic data (Moscariello et al., 2014); this reveals a set of structures located below the Bornes Plateau, at the northern side of the Saleve ridge and the Jura foothills. The Mesozoic sequence consists, at its base, of Triassic units formed by siliciclastic (Buntsandstein), dolomitic, evaporitic (Mushelkalk), and alternations of evaporitic and dolomitic to marno-dolomitic sequence (Rhetian) that can reach up to 500 m in thickness. The Jurassic is mostly composed of competent, often massive, carbonate deposits, intercalated with marls (Dogger), and towards the upper part (Malm) this interval locally shows enhanced porosity values thanks to the presence of biothermal reef facies (base Malm) and to fault corridors that cut through sectors of the GB. The Lower Cretaceous, represents the top of the Mesozoic sequence in the GB subsurface and is composed of fine-grained limestones makes it an important geothermal target in the GB. In fact, it can be intensely fractured as shown by the drilling results of the GEo-01 well, where more than 70% of the total flow rate is discharged from the Cretaceous level. The Cenozoic sequence is composed of the Lower Freshwater Molasse (LFM), mostly composed of alternating fine-grained sandstone, marls, clays, and subordinately lacustrine carbonates. Borehole records reveal that the LFM reaches a thickness of 1300 m in the southern part of the Geneva region, where the Thônex-01 well is located. This unit is mostly impermeable and therefore considered to form the cover of the geothermal system. However, utilization for heat storage installation in the sand-rich intervals of the LFM unit of the areas of sandstone levels is under investigation . The top of the sequence is mostly composed of heterogeneous, up to 120 m thick glacial, fluvio-glacial, and glacio-lacustrine sequences of Quaternary age (Moscariello, 1996;Fiore, 2007). These Quaternary units host the main freshwater resources of the Canton of Geneva and are of great interest for shallow, low-enthalpy geothermal installations.

Gravity dataset of the study area
The Swiss Molasse Plateau is part of a large negative anomaly that characterizes the entire northern sector of the Alpine orogen and is associated with the flexural response to Alpine loading. In the Western Alps, the gravity disturbance is characterized by a positive-negative trend like many other mountain ranges (Karner and Watts, 1983). The trend is NW-SE trend, decreasing gently from the Bresse Graben towards the Pennidic Nappes and then increases abruptly to form the Ivrea-body (IB) anomaly (Fig. 2). Over the Swiss Molasse Plateau and the Swiss-French border, land and airborne gravity data have been collected in the past Group et al. (1989); Kahle et al. (1976); Verdun et al. (2003); Massona et al. (1999) and the gravity disturbance shows a NW-SE regional trend, controlled by crustal thickening (Klingelé, 2006). Several gravity studies have been conducted over the Swiss Plateau in conjunction with geothermal exploration. Particularly in Eastern and Central Switzerland gravity studies have been used to delineate geologic structures in the top 4 km to 5 km as geothermal targets. Integration of gravity data and 3D geologic modeling based on reflection seismic data have been proposed to identify deep PC structures in the St. Gallen area (Altwegg, 2015) and Neuchâtel area (Mauri et al., 2017). Furthermore, pseudo-tomography using sequential Butterworth filtering was applied in Northern Switzerland in a very similar geologic context as in Geneva to match the wavelength content of the gravity signal with geologic features located at different depths (Abdelfettah et al., 2014). Integration of 2D seismic wells data has been proposed to calibrate the inversion of gravity data in the GB (Carrier et al., 2020). Gravity data were collected in the Geneva area and surrounding regions for hydrocarbon-resource exploration and research studies. Gravity surveys carried out across the whole Geneva canton reveal the potential of gravity methods in delineating shallow Quaternary features as well as deeper structures such as the transgressive contact of the Cenozoic Molasse on the Mesozoic units (Adamer and Montandon, 2000;Poldini et al., 1963). For this study, a total of 3558 public gravity data from the gravimetric Atlas of Switzerland (Olivier et al., 2002) and from the International Gravimetric Bureau (Wilmes et al., 2009) have been used (Fig 2). The minimum distance between stations is 100 m.

Gravity disturbance vs. gravity anomaly
In geophysics, gravity measurements are used to learn about the density variations of the Earth's interior (Li and Götze, 2001). Since the gravity vector is the sum of the gravitational and centrifugal accelerations felt by a body, geophysicists are usually only interested in the gravitational component of the observed gravity because that is what reflects the Earth's internal density distribution. The first step toward isolating the gravitational component is to remove the effects of Earth tides and instrumental drift among others. If these effects are properly removed, the resulting observations are considered to be solely the sum of the centrifugal acceleration due to the Earth's rotation and the gravitational acceleration produced by the internal density distribution of the whole Earth. The isolation of this particular gravitational component, and its subsequent use for estimating density distributions related to geological structures in the subsurface, are the main goals of the application of the gravity method (Blakely, 1995). The gravity anomaly, which is defined as the difference between the Earth's gravity on the geoid and the normal gravity on the reference ellipsoid, depends only on longitude and latitude and is not a function of the height. Moreover, since it is not defined at the same point, gravity anomaly contains centrifugal accelerations and cannot be considered a harmonic function (Barthelmes, 2009). In contrast, the gravity disturbance is defined as the difference between normal gravity and the Earth's gravity at the same point. As a result, the centrifugal accelerations can be neglected and the disturbance can be considered as a harmonic function (Li and Götze, 2001;Hofmann-Wellenhof and Moritz, 2006;Barthelmes, 2009). For these reasons, the gravity disturbance is the logical choice for representing the gravitational effects produced by the heterogeneous density distribution of the Earth.  (Brentini, 2018;Clerc and Moscariello, 2020); c) Cross-section cutting through the GB (modified from (Moscariello, 2019)) indicating the main geothermal targets. The first step of our approach is to compute the ellipsoidal-produced normal gravity. Nowadays, it is possible to calculate the exact theoretical gravity analytically at any latitude and height. (Li and Götze, 2001) provide an efficient formulation for the closed formula which can calculate normal gravity at any latitude (Θ) and height (h): where R is the sphere radius, G is the gravitational constant, M is the mass of the sphere, and ω is the angular velocity.
The closed formula has been coded in Boule (Uieda and Soler, 2020) a Python library developed for calculating gravity fields for the reference WGS84 ellipsoid as defined by the values given in (Hofmann-Wellenhof and Moritz, 2006). By subtracting normal gravity from observed gravity, we obtain the gravity disturbance at each point station.

Topographic correction
Before a gravimetric survey can be interpreted for anomalous signals, the effect of the topographic masses on the gravity measurements must be calculated and reduced (Hirt et al., 2019). This correction is denoted as a topographic mass reduction (Jacoby and Smilde, 2009) or gravimetric terrain correction (Li and Sideris, 1994). Instead of using the traditional approaches such as planarization that neglect the effect of topographic masses beyond some fixed integration radius (Pasteka et al., 2017), we propose an approach that directly models Earth's topographic masses and their gravity effects using a high-resolution digital elevation model. The gravitational field generated by a layer of prisms in Cartesian coordinates, representing the Earth topography in the study area, has been computed through the analytical solutions given by (Nagy et al., 2000) and (Nagy et al., 2002). In particular, we used Harmonica , a Python library for processing and modeling gravity and magnetic data, that makes use of the modified arctangent function proposed by (Fukushima, 2020) which could be used to apply the topographic correction for small regions where the curvature of the Earth could be neglected. The spatial resolution of terrain correction computations is important to better resolve and detect small-scale or near-surface mass-density anomalies, e.g., in the context of geophysical exploration (Hirt et al., 2019). The elevation data used in this work comes from the Shuttle Radar Topography Mission (SRTM) Version 3.0 for which the voids are filled with non-commercial Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM 002). This elevation data set has a vertical absolute height error of less than 16 m and a circular absolute geolocation error of less than 20 m. For computational reasons, the SRTM elevation model has been upsampled to a cell size of 100 m x 100 m. By subtracting the topographic effect from the gravity disturbance, we obtain the topography-free gravity disturbance.

Trend estimation
In the previous section, we describe the approach to remove the normal gravity and the effect of topography from the measured gravity. The final result is a value of topography-free gravity disturbance (GDTOPOFREE) at each point station. Generally, the next step consists of removing the regional trend from the GDTOPOFREE and interpolate the residuals to produce a grid that can be imaged or contoured. Removal of a trend is an important step when dealing with geophysical data since most interpolation method can struggle with long-wavelengths trends in the data (Uieda, 2018). In this work, we fit a 2D polynomial trend estimated through weighted least-squares regression, to the GDTOPOFREE data.

Uncertainty quantification
In the literature, we can find several gridding techniques such as the minimum curvature method (Briggs, 1974), kriging (Marcotte and Chouteau, 1993), or even better the equivalent layer gridding technique (Dampney, 1969;Cooper, 2000) that has the advantage to take into account the 3D nature of the observations, not just their horizontal positions and the fact that potential fields are harmonic functions (Rubén Soler and Uieda, 2020). However, a major limitation of the gridding techniques is that they only preserve the mean value observed in the data. To overcome this limitation we propose to apply geostatistical simulation methods that have the advantage of preserving the variance observed in the data, instead of just the mean value, as in interpolation. Their stochastic approach allows the calculation of many equally probable solutions (realizations), which can be post-processed to quantify and assess uncertainty. In particular, we propose to use a sequential Gaussian simulation (Deutsch et al., 1992;Pyrcz and Deutsch, 2014) approach to simulate 50 equally probable realizations grids of the gravity residuals. As described in the literature, the general sequential simulation process could be resumed in the following step: 1. Perform a normal-score transformation of the raw data (residuals in our case); 2. Randomly select the next node to simulate in the grid; 3. Estimate the local conditional probability distribution function (pdf) from the observations or already simulate values by simple kriging under Gaussian hypothesis; 4. Sample a random value from the pdf and assign it to the grid cell; 5. Repeat steps 1 to 4 until all the grid cells are simulated, taking into account, for every new cell grid, the values previously simulated.
By changing the seed at each loop, we obtain a new realization. From these realizations, we can then compute first-order statistics that allow quantifying spatial uncertainties.
The stochastic simulations are then compared with an interpolated map of residuals obtained through the python's Harmonica package . In particular, we used an equivalent layer technique in which the point sources are located beneath the observed potential-field measurement points by default Cooper (2000). Coefficients associated with each point source are estimated through linear least-squares with damping (Tikhonov 0 th ) order) regularization. The Green's function for point mass effects used is the inverse Cartesian distance between the grid coordinates and the point source: wherex andx, are the coordinate vectors of the observation point and the source, respectively.  Figure 3 compares the results of a traditional gravity processing approach (Fig. 3b to 3e) that use approximation formulas for computing the gravity corrections, with the proposed approach (Fig. 3f to 3i) that allow computing the normal gravity analytically and that use a forward approach for modeling the effect of topographic masses over the survey area. The far-right column of Figure 3 shows the difference between the two approaches, for each processing step.

Gravity disturbance
The normal gravity, computed using the Somigliana's formula and the free-air correction (Fig. 3b) slightly underestimates the normal gravity computed using the analytical form (Fig. 3f), by almost 1 mGal in regions where the topography is more pronounced, as shown in (Fig. 3j). The minimum, maximum, mean, and standard deviation values for the normal gravity using both approaches are resumed in Table 1.
By subtracting the normal gravity from the observed gravity, we obtain the so-called free-air anomaly for the traditional approach (Fig. 3c) and the gravity disturbance for the proposed approach (Fig. 3g). Both approaches show similar results; the gravity disturbance is positive where the topography is greater and negative where the topography is lower than the WGS84 ellipsoid as defined by the values given in (Hofmann-Wellenhof and Moritz, 2006). As for the normal gravity, the higher the topography, the greater the difference reaching almost 1 mGal over the Jura mountains.

Topographic correction
The results for the gravitational effect of topographic masses using the prisms layer forward modeling are presented in Figure 3h. Figure 3d shows the results obtained with a planar approximation. The minimum, maximum, mean, and standard deviation values for the topographic correction for both approaches are resumed in Table 1. It is interesting to note that the mean value of the two methods differs by about 5 mGal with a maximum value of more than 200 mGal when the terrain corrections rely on planar approaches versus a maximum value of about 156 mGal when the effect of topographic masses is generated by a layer of prisms with higher resolution. The greater differences are generally located where the topography is more pronounced, which is normal since the elevation models used for the two approaches are of different resolutions and by using a coarser grid (as in the planar approximation approach), the higher values tend to be smoothed.
Finally, by subtracting the topographic (terrain) correction from the gravity disturbance (free-air anomaly) we obtain the topography-free gravity disturbance (GDTOPOFREE) and the complete Bouguer anomaly (CBA), as shown in figure 3i and 3e respectively. Figure 3m shows the differences between CBA and GDTOPOFREE, which can reach values of +/-5 mGal.

Trend removal
The last step of the processing workflow consists to remove the regional trend from the GDTOPOFREE to obtain the residuals. A 2D polynomial of fourth-degree has been fitted to the GDTOPOFREE data and the trend is then estimated through weighted least-squares regression using Verde (Uieda, 2018). The trend and residuals are shown in Figures 4a  and 4b respectively. Figure 4c shows the distributions of the residuals.

Imaging the residuals
To image the residuals over the study area, we need to interpolate them over a regular grid. As described above, we propose to compare a well-established gridding method, known as equivalent layer gridding Dampney (1969); Cooper (2000) with a sequential Gaussian simulation technique, a stochastic approach that gives a better representation of the natural variability of the property to image and provides a means for quantifying uncertainty.  The geostatistical approach relies on the quantification of the spatial continuity of the variables we estimate or simulate. The first step for computing Sequential Gaussian simulation is to compute the experimental variogram for a variety of directions (Pyrcz and Deutsch, 2014) and identify the direction for the spatial continuity. As shown in Figure 5b, the azimuth 135 is the major direction of the spatial continuity. By fitting a variogram model we obtain the kriging parameters (nugget effect, variogram structure, and the range) that are needed to compute the sequential Gaussian stochastic realizations. Figures 6a and 6b show the results of two among 50 equiprobable residuals realizations obtained with the SGS approach. Figure 6c shows the conditional variance of the 50 residual realizations. The high values (reddish colors) indicate a high variability on the value simulated among the 50 realizations. The conditional variance increases with the distance from data locations (i.e., the conditional variance is greater when the cell to simulate the residual value is located far away from the gravity stations and is 0 at data points). Figure 7 compares the e-type of 50 realizations (7a) and a gridding result using the equivalent layer technique (7b). The e-type is the local expectation and, for an unlimited number of realizations, it corresponds to a kriging estimation. The equivalent layer gridding is obtained after cross-validate the input parameters. The best result is obtained using a relative depth at which the point sources are placed beneath the observation points of 3000 m and a damping regularization factor of 1.0.
Both the realizations, the e-type, and the equivalent layer gridding shows a set of positive and negative anomalies. The main positive anomalies are associated with the regions where the Mesozoic units are exposed as in the Jura mountains, Figure 5: xperimental variogram computed on residuals for the spatial continuity major direction.
and Saleve Ridge, and in minor parts the Vuache. Three main negative anomalies are observed in the north-west (1), south-west (2), and south-east (3) parts of the study area.
1. The NW anomaly is located in the Jura Mountains, where the Cretaceous and Jurassic limestones are the only exposed lithologies. This anomaly became narrower in the S-SW. The source of this anomaly is unclear. Small outcrops of Triassic limestones and anhydrites are exposed in the surrounding region (Figure 1a), and likely extends at 1-2km in depth below the Cretaceous and Jurassic limestones as indicated by seismic data in the area (Gorin et al., 1993). Triassic lithologies can have lower density compared to the surrounding Cretaceous and Jurassic limestones and are often affected by decoupling along with the main detachment, with subsequent repetitions as observed in the Charmont-1 well, located about 10km to the West from Anomaly A (outside of this map). The contribution of P-C sediments, drilled at rather shallow depth (1793m-2291m MD) can also be considered a possible source of this anomaly, at least partially.
2. This anomaly is located at the Rumilly Basin, where the main formations cropping out are the Quaternary and Tertiary sediments and their total thickness is estimated to be around 500 m. The Musiege-1 well is located at the south-eastern boundary of this anomaly in the proximity of the Vuache Fault, where the well encountered 150 m of Molasse sediments above the Lower Cretaceous limestones Moscariello et al. (2014).
3. The third negative anomaly can be associated with the thick sequence of Molasse sediments filling the Bornes Plateau where the Saleve-2 well reached the base of the Cenozoic sediments at 1750 m and the Faucigny-1 well at 2915 m.
It is not surprising that the aforementioned anomalies show high variability, due to the heteroscedastic nature of geological parameters that commonly show a higher variability for high values and a lower variability for low values (Linsel et al., 2020). This is illustrated in single stochastic realizations ( Fig. 6a and b) as well as in the residuals profile plot shown in figure 8. This figure shows the profile for 2 randomly chosen realizations, the e-type, and the equivalent layer technique along the A-B line. Broadly, the 4 profiles show the same trend; however, small-scale variability is absent in the equivalent layer compared to the SGS realizations. Even the e-type, which for an unlimited number of realizations is equal to a kriging interpolation, shows some small-scale heterogeneity. This figure also shows the profile for the conditional variance. It worth noting that when the variance is higher (0m to 7km and 35km to 40 km), the residual values are higher (both positives and negatives). It is also important to remark that when the conditional variance is close to 0 the Equivalent layer (blue) and the e-type (dark red) profile shows a very similar trend with values close to 0. In fact, when the residuals are close to 0, there is not much variability between realizations and the results of a standard interpolator and a stochastic one are close. Figure 6: Results of the stochastic simulation: a) and b) two randomly realizations; c) Conditional variance map.

Discussion
The properly learn about the Earth's interior, we should compare the observed gravity to that of ellipsoidal-produced normal gravity values at each observation station (Li and Götze, 2001). In exploration geophysics, the normal gravity is commonly approximated using the Somigliana formula (Somigliana, 1929) and the height correction (also called freeair). The difference between the observed gravity and this approximation is the so-called free-air anomaly. Traditionally, geophysicists use the elevation as the vertical position of the gravity station and the topographic model. The elevation is used in all the corrections including the height correction and the complete Bouguer correction. However, as pointed out by (Gumert, 1985), the free-air varies significantly with horizontal position and can affect the reduction of observed gravity data. Land gravity measurements made at varying elevations in an area of rugged topography, processed using the standard accepted free-air factor, can produce highly erroneous maps. Our approach for processing gravity data demonstrates that the free-air correction step is outdated. Thanks to the increased computational performance, nowadays it is possible to use open source software such as Boule , which uses a closed-form expression to compute normal gravity anywhere outside of the ellipsoid. This means that free-air approximation is completely unnecessary. The same is valid for the topographic (or terrain) correction. Even if the Bouguer correction is fast and practical, you need a prism model to correct the flaws of the Bouguer correction. Instead, with the proposed approach we are able to directly compute the effect of topography through the prism layer forward modeling approach, and we do not need to do a Bouguer correction at all.  Imaging gravimetric residuals using a stochastic simulation instead of an interpolation technique has two main advantages for uncertainty quantification. In fact, since stochastic simulation preserves both the mean and the variance observed in the data, they can reproduce the extreme value encountered in the variable distribution. And since the higher variability (uncertainty) is encountered in the extreme values of the data distributions (Linsel et al., 2020), stochastic techniques are best fitted to quantify this variability through variance (or standard deviation) maps. Instead, classical interpolation techniques, preserve only the mean value observed in the data and tends to produce images that are smoothed in which the interpolated value is close to the mean value, and in which the extreme values are underrepresented.

Conclusions
A new approach allowing the processing of gravimetric data in a two-step workflow using analytical approaches instead of approximations is presented. Firstly, the normal gravity is computed analytically and subtracted to the observed gravity in order to obtain the gravity disturbance. Then, a forward approach is used to compute the gravitational effect of topographic masses that is removed to the gravity disturbance. By using Boule and Harmonica Python packages, we demonstrate that gravity processing could be accurately done using open source tools. Secondly, we present a stochastic approach to image the gravimetric residuals. The main advantage of this approach is that each stochastic realization reproduces the distribution of the data, and not only the mean value. This is a key aspect when uncertainties need to be taken into account in the context of decision-making processes supporting geo-energy exploration activities in sedimentary basins.