On the use of gravity data in delineating geologic features of interest for geothermal exploration in the Geneva Basin (Switzerland): prospects and limitations

Gravity data retrieved from the Bureau Gravimétrique International and the Gravimetric Atlas of Switzerland have been used to evaluate their applicability as a subsurface investigation tool to assess key geological features in support of the geothermal exploration in the Geneva Basin (GB). In this context, the application of an effective processing workflow able to produce reliable residual gravity anomalies was implemented as a crucial first step to investigate whether and to what level gravity anomalies can be correlated to geologic sources of geothermal interest. This study focusses on the processing workflow applied to publicly available gravity data, including the quantification of the uncertainty. This was then also used for first-order 2D forward gravity modelling. The resulting residual anomalies demonstrate the potential use of gravity investigations for geothermal exploration in sedimentary basins, and also reveal areas of significant, irreparable misfit, which calls for the use of complementary data and 3D subsurface structural knowledge. The results of such investigations will be presented in subsequent studies.


Introduction
The deployment of renewable energy sources for both power and heat production is accelerating in Switzerland, a trend that will continue as a result of the 2050 Swiss Energy Strategy (Swiss Federal Office of Energy, 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 (500-3000 m) to large depth (> 3000 m) to combine heat, power production and, potentially, mineral extraction.
The Geneva Basin (GB) and its surrounding French region (Plateau des Bornes, Bellegarde, Fig. 1) has been intensively studied for hydrocarbon exploration since the 1960s, and for geothermal exploitation in the 1990s, but only economically non-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 bottomhole temperature of 88 °C at 2530 m (TVD), but 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) L. Guglielmetti , A. Moscariello and 8 bar wellhead pressure Moscariello et al., 2020;Services Industriels de Genève, 2018). The geographic positions of the wells in the study are reported in Fig. 1 and their stratigraphy is reported in Table 1.
The geothermal conditions in the study area have been reconstructed by thermal modelling (Chelle-Michou et al., 2017) and geochemical data  revealing a geothermal gradient in the area ranging from 25-30 °C/km. 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 different depths in the Geneva area (Moscariello, 2016;Fig. 2). At only a few tens of meters in the Quaternary deposits there is potential for shallow, low enthalpy ground-source heat-pump installations. Between a depth of 0.5 to 3 km the porous Cenozoic Molasse and fractured Mesozoic sequence (Allenbach et al., 2017) offers possibilities for both heat extraction  Table 1 Stratigraphy of the exploration wells located in the study area (modified from Rusillon, (2018) and Guglielmetti et al., (2020  and storage as observed at the GEo-01 well where 50 l/s of thermal water at 34 °C is discharged from the upper Mesozoic. At 4 to 5 km depth the fractured Triassic carbonates, porous Permo-Carboniferous (PC) sediments, and fractured crystalline basement have the potential for combined applications such as generation of power, heat and metal extraction from geothermal brines. The reliable identification and characterization of these geological structures is crucial to target productive geothermal reservoirs and thus increase the chance of success of drilling activities. The geophysical investigation method that has demonstrated the best results in many geothermal contexts is reflection seismics. 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).
2D seismic data was mostly acquired in the Geneva Basin (GB) in the 1980s to 1990s for Mesozoic hydrocarbon exploration targets located at 2000-3000 m below the ground surface (Clerc et al., 2015;Moscariello, 2019) and therefore shows some limitations in delineating both shallow Quaternary deposits and deeper (> 3500 m) geologic structures such as the Permo-Carboniferous (PC) grabens. In this context gravity surveys, a standard geophysical subsurface exploration technique in different geological settings (Blakely, 1995;Reynolds, 1996;Telford et al., 1990) can contribute to reduce such limitations. Gravimetry is commonly applied, in the early stages of geo-resources exploration programs, to identify regions of potential interest (Allis et al., 2000;Guglielmetti et al., 2013;Uwiduhaye, 2018) and monitoring production operations (Eysteinsson, 2000;Mariita, 2000;Narasimhan & Goyal, 1984;Portier et al., 2018).
In the framework of an ongoing geothermal exploration program to improve the understanding of the subsurface in the GB, the use of existing and newly acquired gravity data has therefore been proposed. To this end, this study aims at evaluating the possible added value of analysing existing gravity data which can be used as a subsurface investigation tool to assess the basin structure and identify gravity anomalies within the geothermal stratigraphic targets.
In this paper, we aim at (i) assessing the potential of processing the gravity data across the GB; (ii) quantifying the uncertainty associated with the available dataset; (iii) delineating gravity anomalies associated with the main geological features of geothermal interest in the study region; (iv) comparing the gravity signal to the gravity response of geologic models along three cross-sections to highlight the correlations and the misfits between the two signals; (v) discussing the limitations of using gravity data alone and point towards recommended use of gravimetry and complementary dataset for subsequent investigations.

Geological setting
The study area covers an area of about 2000 km 2 extending 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 is the westernmost part of the North Alpine Foreland Basin (NAFB) that extends from the Savoy region in France to Linz in Austria (Kempf & Pfiffner, 2004;Kuhlemann & Kempf, 2002;Mazurek et al., 2006). The NAFB originated as a lithospheric flexure response to the topographic load of the Alpine orogen (Pfiffner et al., 2002) whose crustal model in the study area can be described as a gently dipping surface towards the Alpine orogen with increasing thickness from an average of 20-30 km in the Molasse Basin to 40-50 km in the high topographic region of the Alpine orogen and to 50-60 km south of the Pennidic front.
In the study area, four major lithological units are observed at depth according to drilling records and seismic data (Brentini, 2018;Clerc & Moscariello, 2020;Jenny et al., 1995;Moscariello, 2019;Rusillon, 2018). From bottom to top these are (1) the crystalline basement including PC troughs at the bottom, (2) sedimentary cover composed respectively of Mesozoic carbonate units and, at the top, (3) the Cenozoic and (4) Quaternary sediments. These units can be approximated to a "layercake" model with subparallel formations gently plunging towards SE with an average dip of 3 degrees.
The crystalline basement is only exposed in the Alps, to the South of the study area but has never been drilled by exploration boreholes in the study area. The basement is often affected by SW-NE oriented tectonic depressions originated as pull-apart basins during the Carboniferous (McCann, et al., 2006) and filled with several thousand of meters of PC sediments (Diebold, 1990;Mann et al., 1983;Ziegler, 1990). In the Geneva Basin, PC grabens have been mapped using reflection seismic data (Brentini, 2018;Moscariello, 2014Moscariello, , 2019; this reveals a set of structures located below the Bornes Plateau, at the northern side of the Saleve ridge and at the Jura foothills (Fig. 2).
The Mesozoic sequence spans from the Triassic to the Lower Cretaceous and is bound at the base by the regional Hercynian unconformity. The Triassic lithologies, have limited exposures in the NW part of the study and, according to borehole data, are represented by the German stratigraphic units consisting of: (i) the siliciclastic Buntsandstein, (ii) dolomitic and evaporitic Muschelkalk, (iii) alternance of evaporitic, dolomitic, marno-dolomitic sequence, shale and sandstones sequences of Rhaetian age. The Jurassic sequence is mostly composed by carbonates deposition accumulated in different depositional environments including shallow to deep platform environments, reef and peri-reefal (lagoonal) settings, tidal (wave-dominated) systems (Brentini, 2018;Rusillon, 2018). The Lower Cretaceous is dominated by carbonate sedimentation locally marked by an increase in intercalation of siliciclastic material likely associated with the nearby development of a subaerial landscape during the Valanginian. This subaerial exposure resulted in the development of an erosive and deeply karstified surface (Siderolithic) which, combined with tectonic structures might play a major role in controlling geothermal fluid circulations at depth.
The Cenozoic sequence consists of clastic rocks represented by the autochthonous Lower Freshwater Molasse (LFM), mostly composed of alternating sandstones, marls, gypsum-rich clays and, at the base, thin lacustrine carbonates. Borehole records reveal that the LFM attains a thickness of 1300 m in the southern part of the Geneva Basin, where the Thônex-01 well is located. The LFM is separated from the Subalpine Molasse (SAM) by the Subalpine Molasse Frontal Thrust (SAFT) (Chaudhary et al., 2014) and has been crossed by the Mt. Boisy-1 well (Rusillon, 2018). At the front of the Prealps, the SAM has been involved in the displacement of the Ultra-Helvetic nappe, thrusting over the autochthonous Cenozoic-Mesozoic sequence as observed in the Faucigny-1 wells (Rusillon, 2018).
The top of the sequence is composed of Quaternary (Würmian), heterogeneous, up to 120 m thick (Fiore, 2007), glacial, fluvio-glacial and glacio-lacustrine sequences hosting the main freshwater resources of the Canton of Geneva.
The tectonic evolution of the study area is associated with the alpine compressional phase that caused the decoupling of the sedimentary succession from the basement by a detachment surface occurring on the Triassic evaporites (Affolter & Gratier, 2004;Arn et al., 2005;Guellec et al., 1990;Sommaruga, 1999). Additionally, the inherited basement relief and normal faults bounding Permo-Carboniferous troughs may have played a role in the nucleation of the Mesozoic north-westward thrusts observed in the SE sector of the Geneva Basin and Bornes Plateau (Gorin et al., 1993;Signer & Gorin, 1995).
Shortening affected the entire Mesozoic and Cenozoic sedimentary cover of the study area and was absorbed through fold and thrust reliefs of the Jura arc mountains during the late Miocene and Early Pliocene (Affolter & Gratier, 2004;Homberg et al., 2002;Meyer, 2000;Moscariello, 2019) and by the coeval formation of strikeslip faults. The most relevant surface evidence of such structures is the NW-SE Vuache fault ( Fig. 1) (Charollais et al., 2013;Moscariello, 2021), which crosscuts the entire basin and marks the western boundary of the study area ( Fig. 1). The central part of the GB may have undergone local uplift controlled by a sequence of SW-NE oriented intra-basin blind thrusts that created an anticlinal/ synclinal flexure at the top of the Mesozoic (Moscariello, 2019;Signer & Gorin, 1995). The Molasse sequence has been affected by horizontal compaction associated with the horizontal shortening which acted along the main Alpine thrusts (Burkhard & Sommaruga, 1998), as highlighted by local high velocities observed on seismic data (Kälin et al., 1992).
Towards the northeast of the basin, the structural configuration is dominated by E-W striking faults (Fig. 1). NW-SE and E-W strike-slip faults occur as a series of sub-vertical individual faults often affecting most of the Mesozoic sequence, down to the Triassic detachment surface, and their extension through the Cenozoic interval often appears as flower structures (Angelillo, 1983;Charollais et al., 2007;Clerc & Moscariello, 2020;Moscariello, 2019).

Gravity and density data in the study area
The Swiss Molasse Plateau is part of a large negative gravity anomaly that characterizes the entire Alpine orogen and is associated with crustal thickening and the flexural response to Alpine loading. Crossing the main anomaly in a perpendicular way from NW to SE, it decreases gently from the Bresse Graben towards the Pennidic Nappes and then increases abruptly due the Ivrea-body anomaly ( Fig. 3) (Kissling, 1984;Scarponi et al., 2020Scarponi et al., , 2021. The use of gravity data in delineating lithospheric structures across the Alps has been demonstrated by 3D gravity modelling (e.g., Spooner et al., 2019).
Over the Swiss Molasse Plateau and the Swiss-French border, land and airborne gravity data have been collected in previous studies (Bayer et al., 1989;Kahle et al., 1976;Massona et al., 1999;Verdun et al., 2003) and the BA shows a NW-SE oriented regional trend, controlled by crustal thickening (Klingele, 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-5 km as geothermal targets. Altwegg et al. (2015) integrated gravity data and 3D modelling to distinguish between Quaternary deposits and deep PC strata based on density variations. A similar modelling approach based on the combination of gravity data, 3D litho-constrained gravity forward modelling, and stripping processing aimed at delineating a deep PC structure in the Neuchâtel area, was applied by Mauri et al. (2017). Furthermore, "pseudo-tomography" consisting of the application sequential Butterworth bandpass filtering to the gravity signal (Abdelfettah et al., 2014), was applied in Northern Switzerland in a very similar geologic context as in Geneva. This study allowed matching the gravity anomaly wavelengths associated with the different geologic features of the region: shallow Cenozoic deposits were identified in the residual data applying 10 km cut-off wavelength and the top of crystalline basement was constrained by applying the 70 km wavelength limit.
Gravity data were collected in the Geneva area and surrounding regions for hydrocarbon-resources exploration and research studies. Gravity surveys carried out across the whole Geneva canton (Adamer & Montadon, 2000;Poldini et al., 1963), revealed the potential of gravity in delineating shallow, Quaternary features, and also deeper structures such as the morphological features of the transgressive contact of the Cenozoic Molasse on the Mesozoic units, defined by sequences of NE-SW oriented negative and positive anomalies crossing the study area. Over the study area, public gravity data from the gravimetric Atlas of Switzerland (Olivier et al., 2002) for the Swiss sector and from the International Gravimetric Bureau (BGI) for the surrounding French areas are available as well as density, porosity, and permeability values from literature (Table 2).

Methods
The processing and analysis of the existing gravity data across the GB was carried out using the software Geosoft Oasis Montaj (2018). The study aimed at applying a processing workflow to assess the uncertainty on the gravity data and produce a residual gravity anomaly, that can be possibly representative of the gravity effect of the main geologic features of geothermal interest and also highlight areas where significant improvements are required for a comprehensive understanding of structures. The specific processing steps undertaken were: (a) Combination of the absolute gravity data from the Gravimetric Atlas of Switzerland available from the Swiss Federal Office of Topography (Swisstopo https:// www. swiss topo. admin. ch/) and from the Bureau Gravimétrique International (BGI http:// bgi. obs-mip. fr), in order to produce a harmonised dataset. Quality control on observed gravity data was performed in order to preserve as much as possible the data available by (1) (Sheriff, 1984), free-air, Bouguer slab corrections to respectively produce FAA (Sheriff, 2002), and Bouguer anomalies BA, as well as the terrain correction to generate a complete Bouguer anomaly CBA. The 1967 International Gravity Formula (IGF67) for datum reference was used for free-air and Bouguer slab corrections. The standard reduction density of 2670 kg/m 3 , assumed to represent an average upper crustal density (Hinze, 2003), was used for the Bouguer and terrain corrections. The BA was computed according to the correction formula: used in Geosoft Oasis Montaj, where: BA: Bouguer anomaly (mGal); FAA: free-air anomaly (mGal); ρ: reduction density (kg/m 3 ); Z s : station elevation (meters); ρ w : water density, set to 1000 kg/m 3 . The terrain correction is calculated by the software using a combination of the methods described by Nagy (1966) and Kane (1962). The algorithm sums the effects of four tilted triangular sections, which  describe a surface between the gravity station and the elevation at each diagonal corner. In the inner zone, the terrain effect is calculated for each point using the flat-topped square prism approach of Nagy (1966). In the far zone, the terrain effect is derived based on the annular ring segment approximation to a square prism as described by Kane (1962). For more processing efficiency, the far zone calculation is be optimized by de-sampling the outer zone to a coarser averaged grid having cell size of 960 m. The calculation is carried within the local correction distance, which was set to 1 km at full resolution of 30 meters, up to the specified outer correction distance which was set to 167 km. (c) Gridding of the gravity data: spectral filtering (vertical derivatives and bandpass) requires gridded data, therefore setting the optimal cell size of the grid is a crucial step in the transition from scattered gravity data to a regular grid. Gravity stations in the study area show an irregular spacing, ranging from 50 m along the Geneva Lake to about one station every 2500 m at the boundaries of the study area, with a mean value of 250 m calculated using a nearest neighbour algorithm over the entire region. A sensitivity study was carried out on the observed gravity data (G_Obs) using different cell sizes configurations (i.e., 100, 250, 500 and 1000 m), gridding algorithms (minimum curvature and kriging algorithms). The misfit between the gridded and nongridded G_Obs values was calculated, and the associated uncertainty was quantified. (d) Uncertainty quantification on gravity data: gravity uncertainty can be considered as the combination of different components including data acquisition (geographic positioning, elevation), instrumental drift, and processing such as gridding methods. In our case, original data include only the coordinates, elevation, observed gravity, free-air and Bouguer anomalies, therefore it was not possible to analyse the uncertainties related to the acquisition, drift and Earth tides corrections. In order to identify the regions affected by the highest level of uncertainty we analysed the elevation data and carried out a sensitivity study about the gridding method. With respect to elevation data, we calculated the misfit between the CBA resulting from using elevation data from the original gravity dataset and the CBA using elevation from the ASTER GDEM sampled at each gravity station. Free-air, Bouguer anomalies as well as topographic correction were calculated using Geosoft Oasis Montaj to have a common processing platform and avoid possible sources of inconsistencies in the BA datasets from BGI and Swisstopo. Taking into account the elevation from the data and the elevation from the ASTER GDEM used for the terrain correction allowed to consider a broad source of error as the elevation from the acquisition should be more accurate than the ASTER data. For that reason, we used the elevation from the data to perform the free-air correction, and only used the ASTER for the terrain correction. To quantify the uncertainty due to the gridding method and cell size the misfit between irregularly spaced and gridded G_Obs resulting from two gridding algorithms (minimum curvature and kriging) and four different cell sizes (100, 250, 500 and 1000 m) was analysed via sensitivity study. (e) Residual anomaly calculation: the goal of the study was to identify gravity anomalies generated by geological structures located in the upper 5 km in the subsurface, therefore a trend removal was applied to produce a residual gravity anomaly. The regionalresidual separation is commonly accomplished using polynomial trend surfaces, wavelength filtering and upward continuation, and the minimum curvature techniques (Dobrin & Savit, 1988;Gupta & Ramani, 1980;Mickus et al., 1991). We accomplished this step by wavelength filtering after having carried out a sensitivity study where different trends have been computed (including 3rd degree polynomial, upward continuation) over a region extending at least 50 km from the study area and the resulting trends were compared to the trend observed on the CBA in the GB and the best fitting was eventually identified. Therefore, to separate long and short wavelength potential field anomalies, the cutoff wavelengths was be obtained from the calculated radially averaged power spectrum of the data using fast Fourier transform (FFT) (Spector & Grant, 1969;Bhattacharya, 1965). The power spectrum (Fig. 8) of the gravity data shows several linear segments of the curve decaying with increasing wavenumber and when plotted against the wavenumber allows to interpret major density variation zones from the slopes of the different linear segments observed. The slope of each straight-line curves fitted to the different spectral segments indicates the average depth of the different causative bodies (Hahn et al., 1976;Rama Rao et al., 2011;Reeves, 2005). Eventually, power spectral analysis was also used to identify and filter the noise portion of the signal, constrained by the shorter wavelengths unsubstantiated by exposed geology. This ensures preserving the gravity signal from the local geologic sources of geothermal relevance while removing the effect of regional structures and the noise component.
(f ) Calculation of the gravity response of 2D models via forward modelling along three cross-sections passing by the main deep wells in the study area: forward modelling is a useful tool to validate the interpretation of gravity data. Forward models were computed via the GM-SYS profile modelling tool embedded in Geosoft Oasis Montaj. The methods used to calculate the gravity model response in Geosoft Oasis Montaj is based on GM-SYS program (2007) and are based on the methods of Talwani et al., 1959, Talwani & Heirtzler, 1964, and make use of the algorithms described in Won & Bevis, 1987. GM-SYS Profile uses a two-dimensional, flat-earth model for the gravity calculations, where (1) each structural unit or block extends to ± infinity in the direction perpendicular to the profile, (2) the Earth is assumed to have topography but no curvature, (3) the model also extends to the far-field (practically: ± 30,000 km) along the profile to eliminate edge effects. Geologic data from the publicly available GEOMOL model were used as constraints for the geologic contacts, and the calculated gravity response was compared to the residual anomaly. GEOMOL data are available as 3D grids which for the scope of this study, were sampled along the 3 selected profiles. Homogeneous density was assigned to each formation, according to the average values reported in Table 2

Uncertainty quantification
Gravity data are affected by uncertainties associated with both acquisition and the processing. Elevation data are some of the most relevant parameters that control gravity observations; therefore, it is important to assess the uncertainty of such data and in order to identify the regions affected by the highest level of elevation uncertainty. In the Geneva area we compared the elevation data from the gravity dataset and the elevation retrieved from the ASTER GDEM that was deployed in this study for the terrain correction of the gravity data. ASTER GDEM has a pixel size of 0. To assess the effect of elevation uncertainty in the CBA results used in this study, a CBA using elevation data from the original dataset and a second CBA using the ASTER GDEM elevation data were computed and the difference between the two anomalies ("CBA original data" vs "CBA ASTER") was calculated. The results reveal an overall low misfit between the two anomalies (average difference of 0.025 mGal) and the larger misfit up to − 11.6 mGal is observed in mountainous regions (Jura in particular most probably affected by ASTER GDEM large uncertainties) at the boundaries of the study area and of only marginal interest for geothermal exploration in the study area. The large majority of stations (49.9%) fit within a ± 0.75 mGal uncertainty, 97.4% within a ± 3 mGal uncertainty.
The uncertainty associated with the gridding algorithm and cell-size highlights how the 250 m cell size minimum curvature algorithm results in an average misfit between irregularly spaced and gridded data of 0.003 mGal, with 75% of the stations within ± 0.75mGal, 98% within a ± 3mGal uncertainty. As expected, gridded data are affected by data loss as the gridding cell-size increases and the misfit between gridded and non-gridded increases accordingly. Eventually, the grid resulting from a minimum curvature interpolation using a 250 m cell size resulted in the optimal compromise. In order to assess the uncertainty due to gridding cell size, the selected 250 m cell size minimum curvature algorithm was analyzed revealing an average misfit between nongridded and gridded data of 0.003 mGal, with 73% of the stations within ± 0.75 mGal, 96% within a ± 3 mGal uncertainty.
According to these results the quantification of the overall uncertainty shows that 57.5% of the stations are within ± 0.75 mGal, 94% within ± 3 mGal and higher values of uncertainty can be identified at the mountainous regions and at the boundaries of the study (Figs. 4 and 5).

Gravity data anomalies
The elevation, observed gravity, free-air anomaly (FAA), BA, terrain correction, CBA, and regional trend maps are shown in Fig. 6a-g.
The FAA in Fig. 6c shows values between − 71 mGal and 98 mGal and the anomaly distribution is controlled by the elevation of the gravity stations, hence the local topography, and by the local geology. Positive values of the FAA are controlled by the main topographic heights of the Jura and Vuache Mountains, the Salève Ridge, le Signal des Voirons on the East and the Mt. Sur Cou peak L. Guglielmetti , A. Moscariello on the south-east, which also represent the main areas where the Upper Mesozoic carbonates are exposed. In contrast, the FAA lows correspond to the main basin structures filled with lower density Cenozoic sediments such as the Geneva area, the Rumilly basin, and the Arve Valley.
The BA (Fig. 6d) shows a NW-SE oriented linear trend, ranging from − 136 to − 37 mGal. Local anomalies with respect to the regional trend can be observed associated with the main peaks such as Crêt de la Neige and Crêt d'Eau peaks across the Jura Mountains, the Salève Ridge and the Signal des Voirons Mountain. This effect was removed by applying the terrain correction (Fig. 6e) to produce the complete Bouguer anomaly (CBA) in Fig. 6f. CBA ranges between − 113 and − 34 mGal and is characterised by a NW-SE linear trend (Fig. 6g), coherent with the regional trend observed for the Molasse Plateau and for the Western Alps at the regional scale.

Residual anomaly
In order to separate the regional trend produced by the regional geologic structures and gravity effect of more shallow structures including those of geothermal interest, a sensitivity study was conducted by computing a CBA over a region extending at least 50 km outside the study area. The resulting CBA data have been gridded using a minimum curvature algorithm with the cell size of 250 m to be consistent with the resolution of the grid covering the Geneva area. The SW part of the region is poorly constrained by only scattered and broadly spaced stations, but still provides constraints to define the regional trend. The CBA shows values ranging between − 179 and 58 mGal, with higher values in the NW, lower values in the central part and a sharp increase observed on the SE part of the area (Fig. 7). The power spectrum of the CBA reveals that the wavelength content over the entire region is within 315 km (Fig. 8) and by analyzing the inverse of cut-off wavenumber 0.01 km −1 (purple line in Fig. 8), it gives a cut-off wavelength of 80 km, therefore it can be estimated that the wavelengths between 80 and 315 km (the largest wavelength being controlled by the size of the investigation area) are associated with sources located between 21 and 30 km depth that is consistent with the crustal thickness value of 25-30 km (Bassin et al., 2012;Mooney et al., 1998;Reguzzoni & Sampietro, 2015). Wavelengths between 80 and 35 km are associated with sources located below 10 km depth whereas sources between 1.5 and 10 km generate anomalies with wavelengths larger than 3.5 km. Such geologic sources are coherent with the depth interval where the lower Mesozoic, the Permo-Carboniferous structures are located in the Geneva basin. Shallower sources within the waveband between 1 and 3.5 km can be generated by sources located between 0.2 and 1.5 km matching with the shallow sources in the GB corresponding to the Quaternary, Cenozoic and upper Mesozoic. Wavelengths shorter than 1 km which, in the GB do not show any correlation with known exposed geologic features (i.e. Quaternary glacial valleys) are therefore considered as noise. As a next step, 80 km cut-off low-pass Butterworth filtering was applied to the CBA and the results show that a good fit with both the trend of the regional and the Geneva area CBAs is observed (Figs. 9 and 10) which is coherent with the results from the power spectrum and previous studies on similar regions (Abdelfettah et al., 2014regions (Abdelfettah et al., , 2020. The residual gravity anomaly resulting from the separation of the regional trend ranges between − 9.7 and 9.7 mGal and shows strong correlation with the Fig. 6 a Elevation map; b observed gravity; c free-air anomaly; d Bouguer anomaly; e terrain correction; f complete Bouguer anomaly; g gravity regional trend L. Guglielmetti , A. Moscariello local geology in the GB (Fig. 11). Positive values are observed in regions where the Mesozoic units are exposed, in particular across the Jura and Vuache Mountains and along the Salève Ridge; four main negative anomalies also can be identified: • Anomaly A has an amplitude of − 7.5 mGal, shows a 20 km SW-NE oriented elongated shape and is located in the Geneva Basin where the residual anomaly is larger in the NW and SW sectors and smaller in the SE part. This distribution correlates with the thickness of the Cenozoic units as constrained by the GEo-01 well (408 m), the Humilly-1 well (439 m) and the Thonex-01 well (1331 m). In addition, the presence of a 3-5 km wide PC trough extending along the Salève front has been identified on 2D seismic data (Gorin et al., 1993;Moscariello et al., 2014) and was drilled at the base of the Humilly-2 well (Brentini, 2018;Rusillon, 2018); it can contribute the observed negative anomaly. Density data available reveal that the density of the PC sediments is highly variable, and the density contrast up to 0.15 kg/m 3 with respect to crystalline rocks suggests that such structures can be the source of the larger wavelengths composing this negative anomaly. • Anomaly B is located in the Jura Mountains, where the Cretaceous and Jurassic limestones are the only exposed lithologies. The anomaly has a rather circular shape of about 10 km in diameter and an amplitude of − 5.1 mGal. The source of this anomaly is unclear. Small outcrops of Triassic limestones and anhydrites are exposed in the surrounding region (Fig. 1a), and likely extend at 1-2 km 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 the main detachment, with subsequent repetitions as observed in the Charmont-1 well. The Charmont-1 well, where repetitions of Jurassic-Triassic units were drilled, is located about 10 km to the West from Anomaly B, potentially contributing to the observed negative anomaly. The contribution of PC sediments, drilled at rather shallow depth (1793-2291 m MD) can also be considered as a possible source of this anomaly, at least partially. • Anomaly C shows an amplitude of − 3.2 mGal and 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 500 m. The Musiege-1 well is located at the 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). • Anomaly D is located in the Bornes Plateau and extends for more than 20 km in the SW-NE direction and is about 10 km from NW to SE and has an amplitude of − 6.8 mGal. Three wells are located at the boundaries of this anomaly: the Salève-2 well, on the SW, drilled 1750 m of Tertiary Molasse sediments after 25 m of Quaternary deposits; the Faucigny-1 well, the deepest of the region, drilled 1576 m of Cenozoic Molasse, after an interval of 1339 m composed by Quaternary deposits and Prealpine carbonate units and also penetrated 359 m of Permo-Carboniferous sediments at bottomhole (Brentini, 2018;Rusillon, 2018); La Balme-1 well drilled 508 m of Cenozoic sediments before entering the Sub-Alpine carbonate series. The source of Anomaly D is most probably a coupled contribution of the Cenozoic sediments with a PC graben at greater depth, as constrained by seismic data and the Faucigy-1 well.

Comparison between residual anomaly and results from 2D forward modelling
The analysis of the residual anomaly focussed on the correlation of the gravity signal to the known geological features in terms of horizontal distribution, amplitude and wavelength. For this study, borehole stratigraphy and the data retrieved from the GEOMOL model, provide vertical lithologic constraints regarding thickness variations of the different units. These data have been used to compute 2D gravity forward models along 3 profiles and the results have been compared to the residual gravity anomaly (Fig. 12). Profile 1 in Fig. 12, runs NW-SE from the GEo-01 to the Faucigny-1 wells, intersecting the Thônex-01 well. The residual anomaly shows a full amplitude variation of 11.5 mGal and the gravity response of the model of 18.4 mGal. Both the residual anomaly and the response of the model show a decreasing trend from higher values in the NW and lower values in the SE in the Thonex-01 well area. This observation is coherent with the increase of thickness of the Molasse sediments along this axis as constrained by the stratigraphy of the GEo-01 well, where the Mesozoic carbonates are located at 408 m b.g.l. and the Thônex-01 well where the top of the carbonates is at 1331 m b.g.l. At Thônex-01 a relative gravity minimum having 3 km in wavelength is observed in both anomalies corresponding at Anomaly A. On the residual signal the amplitude of this anomaly is 2.8 mGal and on the gravity response it shows a value of 2.7 mGal; however they are offset by ca. 2-3 mGals as in the entire central part of the profile. On the SE of the Thônex-01 well, the GEOMOL data reconstruct a syncline at the level of the Upper Mesozoic, representing the NE extension in the subsurface of the Saleve ridge. This geologic feature is visible as a relative maximum in the gravity signal of both anomalies. However, the residual anomaly does not decrease further and stays at higher levels than the GEMOL response. In the proximity of the Faucigny-1 well the two gravity anomalies diverge with the residual anomaly showing a relative maximum of 2.8 mGal in amplitude and 4 km in wavelength, whereas the response of the model a large negative anomaly of 10.1 mGal and a wavelength of 10 km. The offset reaches 9 mGal. This can be explained by the fact that the GEOMOL model did not took into account the presence of the UHN units, which crop out in the vicinities of the well and have been drilled in the upper 1911 m before entering the Molasse down to 2452 m TVD. Being the SAM composed by a complex stack of carbonates and compacted marls and shales thus resulting denser than the LFM sediments (Kälin et al., 1992), these units produce an increase in the residual anomaly, which is not appearing in the modelled data. Overall, the GEOMOL and density-based 2D model does Gravity in Geneva Basin Profile 2 in Fig. 12, runs SW-NE from the Musiège well to the Mt. Boisy well, cutting the large 20 km wide Anomaly A. The residual anomaly has a full amplitude variation of 9.5 mGal and gravity response of the model of 22.4 mGal, with again significant offsets approaching one end of the profile. Along this profile, both the residual anomaly and the model results show a gravity maximum with amplitude of 2.8 mGal in the residual anomaly and 7.9 mGal in the results of the model, in correspondence of the Musiege well, where the Mesozoic carbonates along the Vuache mountain are exposed. Both anomalies show a decrease in the SW-NE direction reaching a minimum at the Thônex-01 well where the thickness of the Cenozoic sediment is the highest in the region. This trend is constrained by the increase in depth of the top Mesozoic from the Humilly-2 well (439 m TVD) towards the Thônex-1 well (1331 m TVD). In the residual signal a negative anomaly, corresponding to Anomaly A, with wavelength of 10.5 km and amplitude of 3.7 mGal is observed. This feature is not fully constrained in the gravity response of the model, with ca. 2-3 mGal offset along the central part of the profile. Additionally, at the Humilly-2 well relative maximum is observed in both anomalies, having amplitude of 0.5 mGal with wavelength of 5 km in the residual anomaly, and of amplitude of 4.6 mGal with wavelength of 3.5 km in the results of the model. This gravity feature can be explained by the presence of the antiform structure affecting the Upper Mesozoic. Between the Thônex-01 and Mt. Boisy wells the residual anomaly and the gravity response of the model diverge, up to as much as 12-13 mGal. The residual anomaly shows an increase whereas the gravity response of the model shows a minimum due to the large thickness of the Cenozoic sediments modelled in this area, and also constrained by the Mt. Boisy stratigraphy. One possible explanation of the increase trend observed in the residual anomaly can be identified in the lateral contrast at the level of the Cenozoic sediments between the LFM and the SAM, thrusted along the SAFT. The SAM faced a higher degree of compaction (Kälin et al., 1992) due to the Alpine tectonics, possibly reflecting on a higher density of the SAM with respect to the LFM. However, the lack of density data available on this region does not allow to push further the interpretations, and non-negligible misfit remains between processed data and the GEOMOL-based 2D model.
Profile 3 in Fig. 12, covers the southern part of the study region from the Musiège-1 to Faucigny-1 wells. Along this profile, the residual anomaly and the results of the modelling show a good similarity with respect to the distribution of the gravity anomalies, and three main negative anomalies can be identified, however there is room for improvement near the Musiège well and at anomaly D. The residual anomaly has a full amplitude range of 9.7 mGal and gravity response of the model of 16.6 mGal. The first negative anomaly, located on the West of the Musiége-1 well, is controlled by the Cenozoic sediments filling the Rumilly basin. In the residual signal the negative anomaly shows amplitude of 3.2 mGal compared to the 1.3 mGal observed in the gravity response of the model and, in both anomalies, the wavelength is 3.5 km. Between the Musiège-1 and Salève-2 wells, a second gravity anomaly is observed in both gravity signals. In the gravity response of the model, it has amplitude of 6.6 mGal and wavelength of 9 km and shows a rectangular shape controlled by the modelled steep geometry of the Cenozoic/top Mesozoic contact in the Vuache and Saleve areas. The same negative anomaly has amplitude of 3.2 mGal, wavelength of 7 km and displays a regular trend towards SW in the residual anomaly, which is more consistent with the expected geometry of the foredeep limited to the SE by the Salève mountain. The latter controls the large gravity maximum observed in both anomalies on the SW of the Salève-2 well constrained by amplitude of 3.2 mGal in the residual anomaly and 3.6 mGal in the response of the model and wavelengths of 3.6 and 2.9 km respectively. A large negative anomaly can be observed in both residual and modelled gravity signals between the Salève -2 and the Faucigny-1 wells and corresponds to Anomaly D, though with nearly 3 mGal offset over more than 5 km distance. This is 15 km wide and 5 mGal in amplitude for the residual anomaly and 8 mGal for the modelled response. The source of Anomaly D is believed to be a coupled contribution of the Cenozoic and the deep PC graben infill sediments as identified by Moscariello et al. (2014) and Moscariello (2019). A relative maximum of 1.4 mGal in amplitude and 2.3 km in wavelength is observed on the residual anomaly around the Faucigny-1 well, but the gravity response of the model shows a relative minimum 1.8 mGal in amplitude and 2.3 km in wavelength, confirming the opposite behaviour of the two anomalies already observed on Profile 2.
Overall, the three presented sections clearly demonstrate that the GEOMOL and local-density based geometric-physical models are not able to explain the observed residual anomalies sufficiently well. In some regions, the differences between the observed and calculated anomalies are several mGals over several kilometres, therefore interpreting the geothermally interesting structures with such simple 2D approaches could lead to misleading interpretations. Although gravimetric data are essential to provide a control on density structure over the entire zone of interest, it seems that complementary datasets are required to jointly constrain properties and also to provide more reliable, 3D structural input. Re-interpretation of seismic data, with particular focus on the Quaternary, Molasse and Upper Mesozoic units, with gravity data facilitating the interpolation between directly imaged lines, will clearly improve the fit of the models to the gravity anomalies derived from observations. Such studies will be addressed in a different paper.

Conclusions
The results of this study are based on the application of a processing workflow of public gravity data. This includes uncertainty assessment, Bouguer anomaly processing and 2D forward modelling. This latter has shown encouraging results in assessing geological structures of interest for geothermal development, but also highlighted serious limitations with respect to geologic models when considering heterogeneous density distribution across the different geologic units.
The uncertainty assessment of the gravity data confirmed that elevation and gridding methods must be analysed carefully before proceeding to the further steps of processing and modelling, to identify those regions where interpretations can be reliable developed and regions where data shows limitations.
The regional trend separation achieved with wavelength filtering technique allowed to produce a residual gravity anomaly which is coherent with the current knowledge of the main geologic structure of the Geneva Basin. Four main areas of negative anomalies have been identified. Anomalies A, C and D, located on the main sedimentary structures have a double source composed mainly by the Molasse sediments and by a deeper source located between the base of the Mesozoic and the PC units. Anomaly B, located in the Jura Mountains, has a rather shallow source, most probably some thick evaporitic level in the Triassic units which locally crop out in the surrounding area. However, a secondary contribution from a PC structure cannot be excluded. Major positive anomalies are associated with the Mesozoic units that crop out in the whole region and are deeply rooted at the level of the Triassic units, but low amplitude positive anomalies at shallower depths can be identified as indicating anticlinal structures associated with blind thrust structures rooted in the Middle Jurassic levels.
The comparison between the residual anomaly and gravity response of the GEOMOL model resulting from 2D forward modelling, highlight the general similarity between the two signals, but also, the limitations related to 2D modelling. Significant differences are observed, as expected, in particular in terms of amplitude of the anomalies and opposite trend between observed and computed gravity signals. These are mostly related to both the simple GEOMOL model utilized for this work and the first-pass gravity modelling presented in this study. These shortcomings are associated to an additional element of uncertainty that could be addressed by considering a more detailed understanding of the lateral changes in thickness of stratigraphic units and of their composition controlling their density distribution heterogeneities. Future studies should include complementary datasets, 3D structural input with possibly reprocessed seismic data to which gravimetry provides interpolation, and quantitative uncertainty assessment due to structural elements and physical properties of units.
In conclusion, the thorough workflow of gravity data processing as presented in this study is an effective tool that can be deployed at the early stage of geothermal exploration activities at a regional scale. When such data is available in sufficient spatial coverage, it has potential, in any sedimentary basin, to help delineating and assessing the most relevant subsurface geological features of geothermal interest. This is best performed in conjunction with other, complementary datasets and geophysical tools, and in the frame of 3D modelling including uncertainty analysis.