Using highly accurate land gravity and 3D geologic modeling to discriminate potential geothermal areas: Application to the Upper Rhine Graben, France

New land gravity data results acquired in northern Alsace were presented. Compared to the available old Bouguer anomaly, we recovered an accurate Bouguer anomaly field showing data uncertainties <65 μgal. A qualitative data analysis using pseudotomographies reveals several negative anomalies suggesting a decrease of the bulk density at the depth of geothermal interest. We have performed a quantitative study on the basis of the existing 3D geologic model derived from a reinterpretation of the vintage seismics. The theoretical gravity response indicates a great mismatch with the observed Bouguer anomaly. The stripping approach was applied, and the stripped Bouguer anomaly indicates that the density values of the Jurassic, but especially for the Triassic, the Buntsandstein, and the upper part of the basement, were overestimated even using the density values measured in the deep geothermal borehole. This suggests that the borehole density values do not reflect the density variations occurring at larger scale. To reduce the Bouguer anomaly during stripping, a negative density contrast should be affected to the Buntsandstein layer overlaying the basement, suggesting that the part located between the Buntsandstein and the upper part of the basement presents a low-density value compared to the reference density, which is not necessarily expected and is not observed in the densities measured in the borehole. Interestingly, a correlation is found between the gravity analyses and the thermal gradient boreholes in the northern part of the study area. For two boreholes, the gravity interpretation suggests a huge density decrease in the Buntsandstein, which may arise from a combination of high-density fracturing and the important quantity of geothermal fluid significantly affecting the bulk density. Analysis of the thermal borehole data suggests that these two boreholes indicate higher geothermal potential compared with the other boreholes.

ABSTRACT New land gravity data results acquired in northern Alsace were presented. Compared to the available old Bouguer anomaly, we recovered an accurate Bouguer anomaly field showing data uncertainties <65 μgal. A qualitative data analysis using pseudotomographies reveals several negative anomalies suggesting a decrease of the bulk density at the depth of geothermal interest. We have performed a quantitative study on the basis of the existing 3D geologic model derived from a reinterpretation of the vintage seismics. The theoretical gravity response indicates a great mismatch with the observed Bouguer anomaly. The stripping approach was applied, and the stripped Bouguer anomaly indicates that the density values of the Jurassic, but especially for the Triassic, the Buntsandstein, and the upper part of the basement, were overestimated even using the density values measured in the deep geothermal borehole. This suggests that the borehole density values do not reflect the density variations occurring at larger scale. To reduce the Bouguer anomaly during stripping, a negative density contrast should be affected to the Buntsandstein layer overlaying the basement, suggesting that the part located between the Buntsandstein and the upper part of the basement presents a low-density value compared to the reference density, which is not necessarily expected and is not observed in the densities measured in the borehole. Interestingly, a correlation is found between the gravity analyses and the thermal gradient boreholes in the northern part of the study area. For two boreholes, the gravity interpretation suggests a huge density decrease in the Buntsandstein, which may arise from a combination of high-density fracturing and the important quantity of geothermal fluid significantly affecting the bulk density. Analysis of the thermal borehole data suggests that these two boreholes indicate higher geothermal potential compared with the other boreholes.

INTRODUCTION
For any subsurface georesource project (oil and gas, mining, geothermal, etc.), choosing the right drilling locations is crucial, especially for those with limited financial support, which is often the case for (deep) geothermal energy projects. An adequate exploration program is critical to project success by accurately delineating highpotential geothermal areas. Classically, the geophysical methods used in deep geothermal exploration are a combination of surface seismic and nonseismic methods (e.g., magnetotellurics, gravity, and magnetics). Each method has strengths, weaknesses, and costs; hence, the choice of geophysical exploration method depends not only on the scientific and technical challenges but also on the available budget.
Recent geothermal projects developed in the Upper Rhine Graben (URG) tend to target fault zones at the sediment-basement interface (e.g., Vidal et al., 2015), where there is sufficient geothermal water Manuscript received by the Editor 8 February 2019; revised manuscript received 22 November 2019; published ahead of production 11 December 2019; published online 12 February 2020. 1 for industrial use. Seismic methods provide an excellent image of the sediments; nonseismic methods document the nature and structure of basement rocks and are also sensitive to porosity changes. Permeable fractured zones in the basement and in the deepest overlying sediments, at Soultz-sous-Forêts, often show density contrasts up to 20% compared to unfractured basement, due to the intense hydrothermal alteration, high fracture density, and natural permeability (Genter et al., 2000). Therefore, their geophysical characterization should benefit from nonseismic methods such as gravity.
In the URG, gravity surveys have been used to identify density heterogeneities in the basement (e.g., Edel et al., 2018), in the crystalline basement of Switzerland (Klingelé and Schwendener, 1984), and more recently to delineate Permo-Carboniferous grabens in the crystalline basement of the Swiss Molasse Basin (e.g., Abdelfettah et al., 2014). The fracture porosity at the Sankt-Gallen geothermal project was estimated using gravity and 3D seismic interpretation (Altwegg et al., 2015). However, gravity anomalies have often been interpreted as lithologic effects (e.g., Campos-Enriquez et al., 1992;Edel and Weber, 1995;Rotstein et al., 2006) located mostly in the overlying sedimentary units. This later interpretation is supported by field observations to the west in the Vosges mountains and to the east in the Black Forest, as well as in several deep geothermal and oil wells within the graben around the towns of Soultz-sous-Forêts and Pechelbronn.
Gravity modeling was also combined with thermal modeling to better understand the distribution and the origin of the thermal anomalies in the URG (e.g., Baillieux et al., 2013;Freymark et al., 2016). A major conclusion was that the thermal anomalies are not necessarily correlated with the sediment thickness, but they originate mainly from radiogenic structures at much greater depths, that is, from the lower crust-upper mantle boundary (Freymark et al., 2016). Areas showing a high thermal anomaly from geophysical exploration are not necessarily economically valuable: They must also contain sufficient amounts of hot water. These features have a strong effect on the bulk density by hydrothermal alteration as well as by the presence of large water volumes and natural fracturing and faulting. Consequently, gravity studies could provide a useful discrimination tool of economically interesting geothermal zones.
Several authors present interpretations of gravity data in the URG in terms of geophysical and geologic features. The high heterogeneity of density values of the basement is well-documented in Campos-Enriquez et al. (1992) and Edel and Weber (1995), who propose a classification of basement rocks following their density values from low to high. Rotstein et al. (2006) add the infill sedimentary covers to the 2D geologic model; that is, 3D effects especially in the basement were not modeled. Edel and Schulmann (2009) extend this approach to 3D, assuming that the infill sedimentary layers were heterogeneous.
The major limitation of the aforementioned gravity studies in the URG (also Baillieux et al., 2014) is the large data uncertainty. The regional database is assembled by merging data from several surveys (e.g., Rotstein et al., 2006;Baillieux et al., 2014) acquired since 1947. The data uncertainty of these surveys could reach 1 mGal at some points. Earlier studies point out the difficulty of merging different databases, for example, to assess the actual location of the gravity measurements (e.g., Rotstein et al., 2006).
We present here new high-accuracy gravity data acquired in 2013 and 2016 with average uncertainties at approximately 0.02 mGal. Qualitative and quantitative interpretation based on an existing 3D geologic model will be shown and discussed, with an emphasis on misfit (understood as the discrepancy between the observations and the model responses) interpretation. We then present the results of the infill sediment stripping and how they allow us to estimate the gravity effect of the porosity arising from fractures and faults. Our work leads to a quantitative approach allowing discrimination of the potential geothermal areas combining gravity and shallow borehole temperature data. In this paper, we will show what we can achieve using accurate gravity data combined with geologic modeling in the geothermal environment. We give an overview of an accurate data processing, 3D forward modeling using finite-element meshing, multiangle horizontal and vertical gradient, pseudotomography added value, and 3D stripping. The benefit of this approach will be demonstrated and discussed by the results obtained from new observed data in northern Alsace.

GEOLOGIC CONTEXT
The URG has been the object of hydrocarbon (e.g., Sittler, 1972;Heritier, 1994) and more recently of geothermal exploration and production (e.g., Genter et al., 2010;Dezayes et al., 2011). The URG is located between northeast France and southwest Germany, extending approximately 350 km from north to south. It is limited to the east by the Black Forest and to the west by the Vosges mountains. The URG underwent an east-west extension from the Late Eocene until the Oligocene (e.g., Dèzes et al., 2004;Edel et al., 2006). It was during this extension phase that the Pechelbronn oil-rich layers were deposited (Villemin et al., 1986). The present-day stress regime is dominated by northwest-southeast compression, favoring the strike-slip behavior of the preexisting faults (e.g., Edel et al., 2006;Valley and Evans, 2007).
Our study area is located in the French portion of the URG, bounded to the south by the city of Haguenau, to the north and east by the French-German border, and to the west by the Vosges mountains ( Figure 1). Geologic structures in the sedimentary cover present the north-south and north-northeastsouth-southwest strike, along the extension of the graben. Several main faults within the study area, for example, the Soultz, Rittershoffen, and Kutzenhausen faults, and basement fractures hydraulically connect these main faults with north-south extension. Genter et al. (1995Genter et al. ( , 1997) present a detailed fracture analysis of the granitic basement of Soultz-sous-Forêts and its surrounding area. The chronostratigraphy of the study area is relatively well-known from several seismic lines acquired for oil and gas exploration and from deep oil wells and recent geothermal wells around Soultz-sous-Forêts and Rittershoffen (Genter et al., 2010;Baujard et al., 2017;Vidal et al., 2017). The geometries of the main structures were defined from 2D seismic interpretation and by reinterpretation of earlier vintage seismic surveys. These seismic lines, though, targeted the Tertiary sediments, and their penetration depth and recording times were insufficient to image the deeper structures. Consequently, the interfaces related to the Jurassic and Triassic formations were not accurately imaged.
Shallow subsurface structures (<5 km) could be subdivided into two main units: sediments and basement. These sediments mainly include (1) the Cenozoic formations (Quaternary, Pliocene, Oligocene, Eocene, and Paleocene) (2) the Mesozoic formations mainly deposited in the Jurassic known as Lias and Dogger, and in the Triassic, known as Keuper, Muschelkalk, and Buntsandstein; and (3) Paleozoic formations, namely, the Permian. The Variscan basement is granitic and dated at approximately 333 Ma, that is, Paleozoic (Cocherie et al., 2004). More detailed geologic information about the URG can be found in Brunnacker and Boenigk (1983) and Genter et al. (2010).

Gravity measurement
In 2016, we acquired 500 new gravity data points as part of the EGS Alsace project. These data cover an area of approximately 20 × 20 km that includes the Soultz-sous-Forêts and Rittershoffen geothermal areas and extends north to Wissembourg (see the red dots in Figure 1). To control the data quality, we established 18 secondary bases (Figure 1), over the whole area, linked to the primary base located close to the GPK-1 Soultz geothermal borehole, where an absolute gravity measurement was also acquired (Hinderer et al., 2015). We used a Scintrex CG5#1317 gravimeter to perform the measurements. As common practice when using a spring gravimeter, the data acquisition followed a loop pattern. Before closing the loop at the secondary base, we remeasured one point on the loop. We then computed the standard deviation for all of the points of the loop and removed the instrumental drift in the same run. The loops were also interconnected ensuring continuity between points and enhancing data quality control. The number of repeat measurements is 97 (78 intraloop and 19 interloop), which represent 19% of all points. A second set of 300 gravity points (the green dots in Figure 1) were acquired in 2013 by a contractor completing the 234 earlier measurements in the southern part of the study area (the blue dots in Figure 1). Overall, 1033 gravity measurements were used in this study (Figure 1).
Positioning of gravity measurement point using a differential global navigation satellite system Elevation is the predominant and critical parameter in high-resolution land gravity measurements; we measured accurate geographic positions at each point. We used a differential global navigation satellite system (DGNSS) to acquire positions with the remote station close to the Rittershoffen geothermal plant. A control point was measured every day before starting the loops to ensure proper measurements. High-accuracy elevations (AE5 cm maximum) were obtained for all measurements leading to gravity uncertainties AE15 μGal.

Geologic modeling
As part of the EGS Alsace and ANR Cantare projects, a geologic model approximately over our survey has been built ( Figure 1) from interpretation of reprocessed legacy seismic data . These surveys targeted potentially oil-rich Tertiary layers; hence, we have little information about deeper seismic velocities (i.e., from the Tertiary), so this geologic model had an average depth uncertainty of 30 m (maximum 103 m) at the Jurassic-Triassic Figure 1. Geographic map of the study area in northern Alsace including the Soultz-sous-Forêts and Rittershoffen geothermal areas. Three gravity surveys were used in this study: the most recent data set acquired in 2016, a data set acquired in 2013, and the earlier gravity database comprising surveys acquired since 1947. Labels 1-7 represent the thermal gradient boreholes where 1 shows the location of borehole F3, 2: F4.1, 3: F4.2, 4: F4.3, 5: F5, 6: F6.2, and 7: F6.3, respectively. interface, assessed by comparing its time-to-depth conversion with the depth measured at the boreholes. The first-order faults as well as major sedimentary layers (Schists, Tertiary unconformity, top Trias, top Muschelkalk, and top Buntsandstein) have been mapped in 3D. The footprint of the geologic model is shown in Figure 1, and its 3D perspectives are shown in Figure 2. The modeling was achieved using Petrel software. The 3D geologic model was built by interpolating existing 2D seismic sections, constrained by borehole information located within the study area. A description of the outcrop geology can be found in Vidal and Genter (2018).

GRAVITY DATA PROCESSING AND MODELING
A high-quality gravity anomaly data set was obtained and used in the data analysis. See Appendix A for more details on how it was produced.
From our data processing workflow, the data uncertainty for each measurement is the sum of (1) the uncertainty during data acquisition and instrumental drift, <10 μGal from PyGrav code (Hector and Hinderer, 2016), (2) the uncertainty from DGNSS (vertical) location, <15 μGal (see above), and (3) the uncertainty from topographic correction, 40 μGal at most in high-topography areas mostly outside our area of interest.
Consequently, the total data uncertainty is AE65 μGal, which allows us to obtain accurate Bouguer anomalies. In the past, this uncertainty was AE700 μgal (a nonpublic report).
After completing the processing workflow, we obtain the Bouguer anomaly map using a density value of 2300 kg:m −3 (Figure 3a). Over the same study area, the previous Bouguer anomaly map is also available but with high uncertainties falling within AE1.5 mGal. The difference between our new Bouguer anomaly and the earlier map is shown in Figure 3b. Note that we show the comparison only in the western part of the study area where the 2016 survey was conducted (see the red dots in Figure 1). We observe that within this area, differences reach a value of 1.5 mGal in the basin and 3.5 mGal in the Vosges mountains. High mismatch values are mostly located in areas of large topography gradients. Because we used a high-resolution digital elevation model (DEM) (2 m), we consider that the topography effect is sufficiently removed by our topography correction. Our new Bouguer map also shows improvements in the western edge of the URG (Figure 3b). We observe differences, mainly for amplitude values, at other locations within the study area, but we will not discuss them here.
For gravity modeling, we used the Grav3Dfem code developed by Abdelfettah et al. (2014). This code uses a finite-element meshing with a homogeneous density value for each finite volume. We adapted the algorithm developed by Pohanka (1988). The main idea is to reduce by one order the volumetric integrals defining the gravity effect for a body located at a distance R to the measurement point. This reduction, from volume to surface integrals, is obtained through the Gauss theorem. In addition, this surface integral can be expressed as a line integral along a closed curve, representing the boundary of the exterior surface. After resolving several mathematical singularities avoiding values dividing by zero, we can compute the gravity attraction value g by ϕðU k;l ðrÞ;V k;l ðrÞ;W k;l ðrÞ;z k ðrÞ;εÞ; (1) where r is the distance between the observation point and the unit element k, G is the universal gravity constant, δ is either the absolute density or the density contrast value in SI units, n k is the normal vector to the surface k delimited by l edges, and LðkÞ is the total number of the edges for a surface k. The geometric functions U, V, and W are built according to the x-, y-, and z-directions taking into account the geometry of the element in the three directions. To avoid (numerical) singularities, we introduce an infinitesimal number ε, which is basically ≤10 −6 representing only 1 μGal on the total gravity value. The accuracy of the computation can be quantified by the parameter Δg ε LðkÞ; (2) where gðrÞ is the theoretical value. The computed gravity anomaly is then bounded by the value AEΔg ε . The complete modeling algorithm is presented by Pohanka (1988). For the finite-element meshing, we used the version 3.0.2 of gmsh software (Geuzaine and Remacle, 2009).

Pseudotomographies
In this section, we qualitatively interpret the results. Our interpretation will focus on the shape and the location of the anomalies rather than on their origin, depth, vertical, and horizontal extension. We discuss these topics in the quantitative interpretation in the section below.
The qualitative interpretation follows the approach developed and used by Abdelfettah et al. (2014Abdelfettah et al. ( , 2016. We applied this methodology to the observed Bouguer anomaly map ( Figure 4a) to build pseudotomographies using different wavelength filters (Figures 4 and 5).
The term "pseudotomography" was introduced by Abdelfettah et al. (2014), and it means one residual anomaly extracted from several residuals analyzed and presented in the same view-panel. The idea behind pseudotomography is the successive interpretation of (b) Difference between the new and old Bouguer anomalies shown by the isovalues superimposed onto the topography of the study area (color scale). High difference values correlate spatially with the topography gradient. The unit of gravity contour is 0.25 mGal. We used Generic Mapping Tools-GMT for drawing figures (e.g., Wessel et al., 2013). residuals computed using different filters according to the target depth, for example, using wavelengths from band-pass filters of 10-20, 10-30, 10-40, and 10-60 km to remove the shallower effects and focus on the deeper part. This approach gives more qualitative information than using only one residual anomaly through the interpretation process (e.g., by subtracting polynomial terms to the Bouguer anomaly and so on). To make a connection to the anomaly source depth, we need a quantitative study using a 3D geologic model, and we compare the theoretical and observed pseudotomographies, both being analyzed with the same filters.
Two kinds of negative anomalies are observed in the pseudotomographies: shallow and intermediate-to-deep anomalies. The shallower one, located at the western side of the study area extends from the towns of Drachenbron to Wissembourg (Figure 4b). The absolute amplitude of this short-wavelength, negative anomaly reaches 2 mGal. This elongated north-northeast-south-southwest feature is hosted on the west by the western border fault of URG. This anomaly remains visible through the other residuals (Figures 4c-4f and 5c-5f) as expected, but at longer wavelengths (i.e., much deeper) it is combined with another negative anomaly located east of Drachenbron and south of Wissembourg reaching the west side of Soultz-sous-Forêts (Figures 4c-4e and 5c-5e), with an amplitude of −5 mGal. Nevertheless, we observe that its amplitude decreases again (Figures 4f and 5f), which means that at its corresponding depth, most likely in the basement, this deep negative anomaly faded away. This behavior could be explained by the presence of a deep positive anomaly below this negative anomaly, so that their combination generates a negative anomaly with decreasing amplitude from −5 to −2 mGal. This deeper positive anomaly suggests that high-density material could be present at a much greater depth (approximately >5 km). The discussed features are visible in Figures 4 and 5.
For the intermediate and deeper parts, two other negative anomalies are observed: one located at Rittershoffen and the second one at Seebach (Figure 4c-4f). The former anomaly shows an amplitude increase from approximately 1 ( Figure 4c) to 2 mGal ( Figure 4d) with relatively the same amplitude reaching a maximum of approximately 2.5 mGal at the deeper part ( Figure 4f). This means that it is likely located at an intermediate depth. The latter anomaly seems different from the former anomaly. Its amplitude increases through the progression of the pseudotomographies reaching a maximum of approximately 6 mGal at much longer wavelengths (from Figure 4c-4f). Its horizontal boundaries are well-defined, suggesting a clear (sub-) vertical interface with a large density contrast. Two other small negative anomalies, located at the north end of the study area, outside the outline of the geologic model, can be observed. One may be the continuity of the Seebach negative anomaly discussed above. The second anomaly is located at the southern east, southeast of Leutenheim (Figure 4c-4f). Consequently, these two anomalies will not be discussed further.
Several positive anomalies are also highlighted. The Soultz and Soufflenheim horsts are clearly identified. At longer wavelengths, the most important positive anomaly is observed at the eastern part of the study area, from Leutenheim to the south to Mothern to the east, passing through Seltz. This positive anomaly coincides with a high-seismic-velocity zone revealed by passive seismic tomography (Lehujeur et al., 2018).

Multiangle horizontal and vertical gradient analysis
Horizontal and vertical gravity gradients are common products used by potential-field interpreters (e.g., Grauch and Cordell, 1987;Marson and Klingelé, 1993). The objective of using a horizontal gradient is to emphasize vertical and subvertical density contrasts, for instance to highlight the abrupt lateral variations caused by faults, basin edges, salt domes, and dikes. A useful feature of the vertical gradient is that it relocates the gravity anomaly tops right above their generating sources for a local structure, that is for structures generating short wavelengths.
We applied multiangle horizontal derivatives (MAHDs) to the complete Bouguer anomaly map shown in Figure 4a. In the case in which the geology follows a preferred orientation, which is the case of almost all geologic situations, we can perform a geologicstrike analysis. The obtained horizontal derivatives computed using angles of 45°, 90°, 180°, and 270°are shown in Figure 6. The angles are taken counterclockwise from the north. Multiangle analysis shows features that appear and disappear according to the analysis angle (e.g., compare Figure 6a with 6b and 6c). A complete interpretation is then possible combining different angle analysis.
Several horizontal contrasts are highlighted by the MAHD results ( Figure 6) and show mainly north-south contrast orientation, which is an agreement with the extension of the known geologic units (e.g., Georg Project Team, 2013). The MAHD results also coincide with several fault locations. We also note that the targeted faults of GRT-1 and GRT-2 (the thick continuous blue line in Figure 6) are offset with respect to the lateral contrast from MAHD (Figure 6b and 6d). This offset can be explained by the known western dip of the fault; therefore, the gravity effect of this lateral contrast will be shifted laterally to the depth of the top of basement.
In Figure 6c, where the Bouguer anomaly is analyzed from the south (i.e., 180°), we identify two bodies: MAG and HdPZ. These bodies also produce magnetic anomalies as shown in Baillieux et al. (2014). The MAG body shows low to intermediate densities whereas the HdPZ body shows high-density material. Their effects are visible in the pseudotomography shown in Figure 4. These bodies are only revealed when the data are analyzed from the south. Note that these magnetic bodies were interpreted as basement heterogeneities.
In Figure 7, the horizontal derivatives are shown, using the 180°angle obtained from the existing data ( Figure 7a) and from merging our more recent data with the existing database (Figure 7b). We identify the east-west structures in the right side of the study area (see the black arrows in Figure 7) and northwest-southeast orientation on the east side (the black arrow in Figure 7a). These shadow effects disappear when we use only the new Bouguer anomaly map. The HdPZ body is not delineated clearly as shown in Figure 6c, when using only the accurate Bouguer anomaly. From here, we can infer that the accurate data acquisition and processing, and the resulting Bouguer anomaly map, helped avoid possible over-and misinterpretations.
The first vertical derivative (FVD) has also been computed from the complete Bouguer anomaly ( Figure 8). The aim of applying the FVD is to highlight and relocate the anomalies according to their origins; it is sensitive to deep and shallow anomalies, but it is more sensitive to shallower effects. In Figure 8, the negative anomaly observed around Rittershoffen (see also Figures 4 and 5) correlates with the observed faults. The fault located at the east of Rittershoffen and midway to Niederroedern (Nieder in Figure 8) also shows a north-south sinuous shape confirmed by the FVD analysis. The same behavior is observed for the fault located halfway between Soultz and Rittershoffen (i.e., Soultz fault).
In addition, the maximum negative values of the FVD (<−3 mGal∕km) seem to correlate with the thickness of Plio-Quaternary sediments. Two thick zones (>50 and <100 m) are identified from the surface (Georg Project Team, 2013): (1) between Hatten and Rittershoffen, which corresponds to the negative gravity trend going from Rittershoffen eastward ( Figure 8) and (2) at Seebach and northward, which corresponds to the negative anomaly shown from Seebach and northward. These gravity effects are also visible in the pseudotomography shown in Figure 4b and interpreted as a shallow gravity effect.
The maximum negative values of the FVD located at the western side of the study area, west of Soultz crossing Drachenbronn north- Figure 5. Pseudotomography results obtained from the observed data. (a-d) Pseudotomographies computed using band-pass filters of wavelengths 10-20, 10-30, 10-40, and 10-60 km, respectively. (e and f) have been obtained using wavelengths of 20-30 and 20-60 km, respectively. The black points show the location of the gravity stations, and the red polygon shows the horizontal extension of the 3D geologic model. ward to Wissembourg, show areas with important lateral contrasts. These lateral and vertical contrasts (Figures 4 and 8) are greater because in many (sub) zones in this area, that is at northwest of the study area, Plio-Quaternary sediments are absent (Georg Project Team, 2013). The succession of the absence and presence of this youngest sedimentary unit creates at some locations a large density contrast at short wavelengths.

Quantitative study
One of the main aims of this study is to quantify the gravity response of the geologic model through postprocessing data analysis and to compare the results with the observed data. The mismatch between the observations and the theoretical responses can be interpreted in terms of gravity anomalies. This allows us to (1) control the geologic model and then provide areas where it accurately satisfies the data and (2) quantify the gravity effect for each geologic unit and quantify their individual contribution to the total Bouguer anomaly.

Theoretical pseudotomographies
In this section, we quantitatively interpret the results. Our approach is to compute the 3D gravity response of the geologic model (presented above) and compare it to the observed Bouguer anomaly map and to the pseudotomography results (Figures 4 and 5). The 3D geologic model gravity effect was computed using the method discussed in Abdelfettah et al. (2014) incorporating the density values presented in Table 1. During the 3D forward modeling, the computation points are coincident to as those of the database, including their elevations (Farr et al., 2007). The computed Bouguer anomaly and the theoretical pseudotomographies are presented in Figures 9 and 10, respectively.
The computed Bouguer anomaly map (Figure 9a) is mainly dominated by two effects: (1) the deepening of the basement eastward and, consequently, (2) the thickening of the sediments from west to east. These two effects produce high values on the west side and lower values on the east side. Their cumulative gravity effect keeps the same positive (higher) and negative (lower) trends; hence, the theoretical response shows a dominant east-west trend. Although the dynamics of the computed ( Figure 9a) and observed (Figure 4a) Bouguer anomalies are similar (approximately 36 mGal), the comparison between them reveals that the gravity features arising from structures cannot completely explain (alone) the observed anomalies and the features of the pseudotomography. Because we used the best model available, structures from seismics, and density values from logging, we expected a theoretical gravity response close to the observed anomaly. Even if the density values are not correct for the whole area, their variations should not greatly change the overall gravity effects. Our results show that this is not the case (Figures 9a and 10a): There are substantial differences between the computed and observed anomalies. This suggests that the structures alone cannot explain this mismatch and other anomalous sources must be taken into account, for instance, density variations within the same geologic unit, originating from lithologic heterogeneities and/or (primary and secondary) porosity. To analyze this mismatch at different depth levels, we compute pseudotomographies using the same wavelengths as those used for the observed data (Figures 9 and 10). At short wavelengths, the observed negative anomaly (Figure 4b and 4c) elongated from northwest of to north of the study area is not reproduced by the model (compare with Figure 9b and 9c). In these observed anomalies, the gravity values are dominated by the density contrast generated by the succession of the presence or absence of the Plio-Quaternary sediments, and because this geologic unit was not included in the model, the mismatch in this area is expected. In addition, this area is located at the edge of the geologic model, so it is not well-constrained and the boundary conditions of the gravity modeling may affect the calculated gravity response there. The pseudotomography mismatch continues at greater depths (Figures 4d-4f, 9d-9f, and 10c). A small but relatively deep basin is also reported in this zone (Edel and Fluck, 1989).
The negative anomaly around Rittershoffen is locally well-reproduced in the shallower part but with a lower amplitude (Figures 4b, 4c and 9b, The thick continuous blue line shows the surface projection of the targeted fault of GRT-1 and GRT-2 shown at the bottom of Buntsandstein, that is, the interface between sediments and bedrock (Baujard et al., 2017). The black lines show the location of the major faults recovered from 2D seismic reinterpretation (from Baillieux et al., 2014) at the top of the basement. The blue continuous lines show the same faults but located at an 800 m depth. Some contrasts are only visible on one panel (i.e., observed from a specific angle), and the combination of other different views reveals the remaining contrast, confirming and refining the other contrasts. The MAG and HdPZ bodies are the magnetic arc granitoids and the high-density phyllite zone, respectively, as shown in Baillieux et al. (2014) recovered from several data compilation. The gravity footprint of the HdPZ is clearly positive as shown in Figure 4d-4f. These two bodies are (clearly) identified only when analyzing the data from 180°(i.e., from the south).  Figure 10 to Figure 5). At Soultzsous-Forêts, the modeled response suggests a positive anomaly from the surface down to the basement. This pattern is observed until about mid-depth ( Figure 4c) and it becomes a negative anomaly deeper (Figures 4d-4f and 5a-5d), which is not reproduced by modeling. The negative anomaly located immediately north of Seebach is the one negative anomaly reproduced from the shallow down to the deep parts; nevertheless, at deeper parts, its amplitude was underestimated (Figures 9e, 9f, 4e, 4f, 10b-10f, and 5b-5f). The negative anomalies between the west and the center of the study area (i.e., between Drachenbron and Seebach) and from south to north (i.e., between Rittershoffen and Seebach northward) are not reproduced by the modeled response (compare Figures 4a-4d and 9a-9d).
A complementary approach helping the quantitative interpretation is the analysis of the misfit resulting from the observed and computed Bouguer anomalies. We can consider this misfit as "an anomaly," and we analyze it according to the depth using successive wavelengths. After computing the misfit, we built the corresponding pseudotomographies as was already done for the modeled response and for the real data. The obtained results are shown in Figures 11 and 12 where several features especially around Soultz northward to Wissembourg and around Rittershoffen and northern of Trimbach are revealed.
The total misfit shown in Figure 11a is computed by subtracting the theoretical modeled gravity response from the real Bouguer anomaly.
The differences could then be interpreted as gravity anomalies. The negative (positive) anomalies mean that we removed more (less) of the gravity effect than was needed to satisfy the observations. The misfit (Figure 11a) shows that the eastern part of the study area is underestimated whereas the western part, around Benheim, Seltz, and Mothern, is overestimated. The peudotomographies built from the misfit reveal that the effect of the unmodeled Plio-Quaternary is actually small (<2.5 mGal) where its maximum effect is located northwest of the study area, between Drachenbroen and Wissembourg (Figure 11b). This effect could also be affected by shallow (<100−200 m) density variations, possibly in the Tertiary. This misfit also reveals that at intermediate depth, probably <1 km, the model does not satisfy the observations around Soultz, Surbourg, and Soufflenheim where the misfit reaches −2 mGal (Figures 11c and 12a). Much deeper (<1.5−2 km), the misfit shows large discrepancies north of Soultz and northward to Wissembourg, where it reaches −4 mGal (Figures 11d and 12b). In this area and at this depth, we probably image the upper part of the basement at a 1.4 km depth (e.g., Georg Project Team, 2013). Continuing to image downward, the misfit reveals strong mismatch values at Soultz-sous-Forêts reaching −5 mGal (Figures 11e and 12c). This interpretation is an agreement with that obtained from 3D gravity inversion (Schill et al., 2010), in which an important low-density area was recovered under Soultz-sous-Forêts. At the same depth level, significant negative anomalies are observed around Rittershoffen and north of Trimbach reaching −2.5 mGal. Downstairs, the misfit highlighted three local areas where the model did not satisfy the data: The first is located at Soultz and northward until approximately 5 km south of Wissembourg city, the second is located between Rittershoffen and Hatten (i.e., immediately east of Rittershoffen), and the third one is located at Trimbach and northward (Figure 12d). Their misfits reach a negative value of −5 mGal. Much deeper, a large negative anomaly of 5-6 mGal is obtained, which indicates the density overestimation at larger depth, probably >3 km (i.e., completely in the basement) (Figures 11f, 12e, and 12f).
The comparison of the misfit and its pseudotomographies with the faults located at the top of the basement (Baillieux et al., Figure 7. Horizontal derivatives computed from 180°using (a) existing data (provided by BRGM, G. Martelet, personal communication, 2016) and (b) when merging new and old data. These two derivatives can be directly compared to that obtained from new data shown in Figure 6c. We can mainly observe that the shadow anomalies (black arrows) disappear when we use only the more accurate Bouguer anomaly. Land gravity and 3D geologic modeling 2014) reveals a good correlation between the fault-density and the negative misfits (Figures 11  and 12). In the areas where the number of faults is higher, the gravity misfit is strongly negative, for instance, at Soultz and northward ( Figure 12c) and at Rittershoffen where the faults' intersection was observed (Figure 12d). Because the rock bulk density is strongly affected by (secondary) permeability arising from fractures and faults, which are not imaged by seismics, and because the misfit correlates with the faulted zones, these negative misfits could indicate areas of higher porosity suggesting higher permeability. Note that in the geologic model, only vertical displacements of the faults were included but not the contribution of the porosity and permeability effects. Oppositely, and at this depth level, other faults do not match with the higher negative misfit although the density of faults is higher, for instance, the faults located immediately western of Rittershoffen, at the same latitude ( Figure 12e). This behavior could mean that not all faults have higher permeability and then not all faults are geothermally interesting.

Stripping
To quantitatively understand the origin of the observed and computed anomalies, better understand the misfit, but also quantify the gravity effect of each geologic unit and understand the corresponding footprint, we conducted a stripping approach. We computed the theoretical gravity effect from the geologic model and subtracted it from the observed Bouguer anomaly. This stripping can be done for the whole geologic model, but it can also be done for a specific geologic unit (e.g., only sediments), or several geologic units (e.g., Jurassic + Trias) as needed. Figure 13 shows the stripping results for individual geologic units, mainly the upper and lower Rupelian, Jurassic, Trias, Buntsandstein, and basement. Table 1. Density values assigned to the geologic model used in the forward modeling. The names in brackets indicate the simplified name used in the body of the paper. Stripping highlights several features. The more we remove the gravity effect, starting from the surface downward, the more the resulting gravity anomaly becomes negative (Figure 13d, 13f,  13h, 13j, and 13l). This means that the greater the depth, the more we overestimate the density values; that is, we remove more gravity effects than needed. Consequently, the density values used in the forward modeling should be decreased at specific locations, for instance, at Soultz, Surbourg, Rittershoffen, and around Seebach. The resulting stripped Bouguer anomalies shown in Figure 13b, 13d, 13f, 13h, 13j, and 13l are directly compared with an initial Bouguer anomaly (Figure 4a) map. The scale is deliberately set constant from −25 to 25 mGal to highlight the evolution and the changes in the stripped (Bouguer) anomalies.
We confirm also the hypothesis that the dominant gravity effect comes from the upper Rupelian and the basement (Figure 13k). Their stacked gravity effect removed from the Bouguer anomaly resulted in a large negative anomaly located at the center and to the west of the study area (Figure 13l).
From the stripping results, we also understand that the range of the observed Bouguer anomalies, approximately 38 mGal, could be reproduced by summing only the gravity effects of the upper Rupelian and of the basement. The misfit between the observed anomaly and this stacked effect is approximately 37 mGal, which means that although the value range is reproduced, the recovered anomalies are of opposing sign (compare Figures 13k and 4a).
Another point highlighted by the stripping concerns the negative anomaly located west of Soultz and Surbourg elongated from northwest to north (see Figures 4 and 5). The gravity effects of the lower Rupelian alone are negative with an amplitude of approximately −2 mGal, whereas the gravity effects of Jurassic, Trias, and Buntsandstein show positive values with an amplitude of 1.0-1.5, approximately 1.5, and 2.5-3.0 mGal, respectively (Figure 13c-13e and 13g); consequently, their cumulative gravity effect will be positive. The gravity effect of the upper Rupelian is negative (approximately −1 mGal, Figure 14a), which suggests that this negative elongated anomaly can indeed be explained by the effect of the upper and lower Rupelian, but its density value should be lower than what we originally used (i.e., 2100 and 2250 kg:m −3 , respectively, according to a reference density of 2300 kg:m −3 ). This elongated negative anomaly is then explained partly by the geometric variation because the recovered amplitude remains lower than the observed magnitude. We need an additional approximately 2-3 mGal to reach the observed amplitude. The lost gravity effect could be due to Plio-Quaternary effects, or much shallower effects, which are not modeled here (compare Figure 13a to Figures 4b and 5b).
The gravity effect for each unit, except for the upper and lower Rupelian, shows positive values under Soultz and Rittershoffen, which does not agree with the observed values recovered by the pseudotomography anomalies (compare Figures 4 and 13). This behavior is also observed in stripped Bouguer anomalies (Figure 13). Although the used Bouguer reference density is 2300 kg:m −3 , the stripped effects remain higher because the resulting new Bouguer anomalies remain negative (Figure 13b, 13d, and 13f-13l), which reinforces our hypothesis that the density values are generally overestimated with respect to the observed values.
The upper Rupelian stripping results show a stripped Bouguer anomaly with a small-scale dynamic range (approximately −15-0 mGal, Figure 14b); that is, the resulting Bouguer amplitude is lower than the original Bouguer dynamic range (Figure 13b). This suggests that the gravity effect of the upper Rupelian unit is accounted for, except for the area west of Surbourg where the Bouguer values remain strongly negative (approximately −18 mGal, Figure 14b). This last behavior could also be generated by heterogeneities that may be found in the upper part of the basement or Figure 10. Pseudotomography results obtained from the computed data. (a-d) Pseudotomography computed using band-pass filters of wavelengths 10-20, 10-30, 10-40, and 10-60 km, respectively. (e and f) have been obtained using wavelengths of 20-30 and 20-60 km, respectively. The black dots show the location of the gravity stations, and the red polygon shows the horizontal extension of the 3D geologic model. at the Buntsandstein as revealed by the former multiangle horizontal derivative analysis (Figure 6c).
The stripping of the lower Rupelian gravity effect reduces the range of the Bouguer anomaly values compared with those obtained from the upper Rupelian stripping (approximately −15-0 mGal, Figure 13b), which is expected. Nevertheless, the values under Rittershoffen and Seebach remain lower than the minimum negative values (approximately −19 mGal) indicating that the origin of this negative effect should be deeper than the lower Rupelian unit.
North of Seebach, the upper Rupelian shows negative values (approximately −6 mGal, Figure 14a) as well as the Jurassic (Figure 13b) and Trias (Figure 13c), whereas the lower Rupelian shows positive values of approximately 1 mGal (Figure 13a). Their cumulative gravity effect yields negative values similar to the observations but with amplitudes lower by approximately 1-2 mGal. Moreover, in this area, the horizontal and vertical derivatives show lateral contrasts, which are also observed in the upper and lower Rupelian, Jurassic, and Buntsandstein gravity effects.
The above analysis gives us a quantitative idea about the amplitude and origin of the anomalies, as well as each individual unit's gravity contribution to the total Bouguer anomaly. However, sequential and cumulative stripping also helps us to quantitatively understand the evolution of the anomalies. The results obtained from the sequential stripping are shown in Figure 14.
Analyzing the cumulative stripping results reinforces some observations done during sequential stripping but also reveals complementary information. The first concerns the total stripped Bouguer results (Figure 14l), where negative values are obtained ranging from −15 mGal for minimum negative values to −26 mGal for maximum negative values for the whole study area. We also observe that the shallower effects of the upper and lower Rupelian as well as the Jurassic are relatively wellrecovered and then removed and that the amplitude of the stripped Bouguer anomaly is clearly reduced by approximately 10 mGal between Niederroedern and Soultz northward approximately 5 km north of Seebach, except around Surbourg (Figure 14a and 14b) where a maximum negative anomaly of approximately −20 mGal remains (Figure 14c). However, after removing the Trias effect cumulated from the upper and lower Rupelian and Jurassic effects, the discrepancy of the stripped Bouguer anomaly (Figure 14h) increases to an approximately −20 mGal maximum amplitude. The affected area is greater than the observed Bouguer anomaly (Figure 4a) and more widespread, from 2-3 km east of Rittershoffen to approximately 3-4 km north of Soultz until approximately 6-7 km south of Surbourg, westward at the end of the study area. This feature started when removing the cumulative effect up to Jurassic but with small mismatching areas (Figure 14f).
After adding the Buntsandstein effect to those computed and stripped previously, the discrepancy zone does not increase in the maximum amplitude compared to the previous stripped anomaly, but it does increase in size. The amplitude of the maximum negative values remains approximately −20 mGal, and the area is located west of the line joining Leutenheim to Niederroedern and Trimbach northward until approximately 5-6 km east of Seebach (Figure 14e) to the western border of the study area. East of this threshold line, the stripped values remain at approximately −10 mGal and increase smoothly to be null at the eastern end of the study area.
The stripping of the whole model effectthat is we added the stripping of the basement to the previous results (Figure 14l) shows, however, very negative values reaching −30 mGal under Soultz, Surbourg, and their surroundings until the western edge of the study area (Figure 14l). Under Rittershoffen, Seebach, Drachenbronn, Figure 11. Discrepancy between the observed and computed anomalies. (a) Misfit obtained by subtracting the theoretical response to the observed Bouguer anomaly and (b-f) pseudotomography computed using high-pass filters of wavelength cutoffs of 10, 20, 30, 40, and 60 km, respectively, applied on the resulted misfit (a). The black dots show the location of the gravity stations, and the red polygon shows the horizontal extension of the 3D geologic model. The white lines indicate the location of the faults at the top of the basement identified by 2D seismic lines. and between them, the stripped anomaly nevertheless shows negative values of approximately −20 mGal, which remain unexplained by the model. Eastward, the stripped values increase smoothly to reach positive values at Mothern, Seltz, and the surrounding areas ( Figure 14l).
These observations mean that the effect of the sedimentary layers (i.e., from the surface to the bottom of the Lower Rupelian) is wellrecovered by the modeling and their densities (contrasts) are relatively well-estimated, whereas for the Jurassic, Trias, and especially for the Buntsandstein and basement, the anomalies are not recovered and their densities are clearly overestimated.

Fracture-related porosity effect on the gravity anomaly
The objective in this section is to assess the fracture-related porosity effect on gravity values. This is necessary to demonstrate the capability of gravity to assess and delineate the faulted zones using gravity data and a 3D geologic model. In Soultz-sous-Forêts and at a 2000 m depth, which corresponds to the upper reservoir, porosity of 3% and 20% was observed for matrix and fracture-related porosity (Genter, 1990), respectively. From that, we will explore fracture-related porosity values of 1%, 5%, and 10% to a maximum value of 15% and compute their effect on the density values. We compute the bulk density following the porosity-density relationship: Rieke and Chilingarian, 1974), where ϕ is the porosity, δ M is the matrix density, δ B is the bulk density, and δ F is the fluid density. The area is considered to be completely saturated, and the density of the geothermal fluid is set to 1.060 kg:m −3 (Baujard and Bruel, 2006). Although, in a recent study, the density of a geothermal brine was estimated to 977 kg:m −3 at 167°C (Baujard et al., 2017), we deliberately choose the higher density value of 1060 kg:m −3 , which gives a smaller fracture-related porosity effect. The density values affected by the fracturerelated porosity obtained from the above porositydensity relationship are reported in Table 2.
We extracted a 2D cross section from the 3D geologic model shown in Figure 2. We chose its location crossing the GRT1 and GRT2 wells of the Rittershoffen EGS project (see line AB in Figure 4), where a deep negative gravity anomaly was observed and the model did not give satisfactory results (Figures 4, 5, 11, and 12). The cross section (Figure 15e) includes the same geologic units as those used in 3D modeling: upper and lower Rupelian, Jurassic, Trias, Buntsandstein, and basement. The lateral extension of the faulted area (Figure 15d) is determined approximately according to the negative gravity anomaly as observed in Figures 4, 5, 11, and 12. Because the stripping approach shows that the model does not explain the observed values, from Jurassic until the top of the basement, we follow this hypothesis and spread the vertical extension of the faulted zone from the Jurassic down to the top of the basement (Figure 15d). We extend this fractured zone only 200 m in the basement. In Figure 15c, we show the location of the faulted zone in the model after merging it in the initial geologic model. We use the density values reported in Table 1 for the upper and lower Rupelian, that is without fracture-related porosity effect, and those of Table 2 for the fractured units. We compute the gravity effect for the original geologic model, that is, without a faulting area as shown in Figure 15e (reported as "Computed Bouguer," see the gray line in the same figure) and the fracture-related porosity effect on the gravity anomaly according to 5%, 10%, and 15% porosity (Figure 15a). We focus our analysis on the center of the profile where the fractured zone is embedded, mainly between the x-coordinates of 1010 and 1018 km.
From the residual anomalies (Figure 15b), we see that the "observed anomaly" is much more negative than the computed anomaly. We need an additional approximately 2 mGal to reduce the mismatch Figure 12. Discrepancy between the observed and computed anomalies. (a-f) Pseudotomography computed using band-pass filters of wavelengths 10-20, 10-30, 10-40, 10-50, 10-55, and 10-60 km, respectively. The black dots show the location of the gravity stations, and the red polygon shows the horizontal extension of the 3D geologic model. The white lines indicate the location of the faults at the top of the basement identified by the 2D seismic lines. between them. The residual anomaly computed using 5% of the fracture-related porosity decreases this misfit by a half (the red line in Figure 15b). The gravity anomaly arising from 15% of the fracture-related porosity is much higher, reaching a mean value of 1 mGal below the observed values. The best fit, which explains the observed anomaly, is the residual anomaly computed using 10% of fracture-related porosity (compare the gray, red, and black lines in Figure 15b).
These results suggest that (1) the observed negative anomaly at the Rittershoffen geothermal area not explained by the modeling could be caused by low density values compared to those obtained from borehole data. The decrease in density could be caused by approximately 10% of fracture-related porositythis could explain why 280 m 3 :h −1 of geothermal brine are produced from GRT-2 (Baujard et al., 2017), and it could confirm the reason for the project's success, Figure 13. Stripping results, where panel (a) is the gravity effect of the Lower Rupelian and (b) shows the corresponding stripped "Bouguer" anomaly, that is the Bouguer anomaly shown in Figure 4a where the effect of the lower Rupelian is removed. Panel (c) is the gravity effect of the Jurassic, and (d) shows the corresponding stripped Bouguer anomaly. Panel (e) shows the unity gravity effect of the Triassic, and (f) shows the resulting stripped Bouguer anomaly. Panel (g) shows the gravity effect of the Buntsandstein, and (h) shows the resulting stripped Bouguer anomaly. Panel (i) shows the gravity effect of the basement, and (j) shows the corresponding stripped Bouguer anomaly. Panel (k) shows the two stacked gravity effects of the upper Rupelian and the basement, and (l) shows the resulting stripped Bouguer anomaly. Labels: SU, Surbourg; RI, Rittershoffen; SF, Soufflenheim; LE, Leutenheim; BE, Beinheim; SZ, Seltz; NI, Niederroedern; MO, Mothern; TR, Trimbach; SE, Seebach; SO, Soultz-sous-Forêts; DB, Drachenbronn; and WI, Wissembourg. and (2) the density of the fracturing of the crystalline basement at Soultz, which reaches 20% (Genter, 1990), is much higher than the fracturation density at Rittershoffen, which was assessed by our calculations to be only 10%. This value of 10% could be higher if the density of the geothermal fluid decreases to reach the estimated value of 977 kg:m −3 by Baujard et al. (2017). Note that these density decreases could also be generated by a petrophysical change within the basement affecting the gravity value by the same amount.

Gravity analysis vs. thermal gradient
The interpretations made in the above sections stipulated that the gravity effect of the upper and lower Rupelian are well-recovered by the model, unlike the Jurassic, Trias, and especially the Buntsandstein and basement effects, which are clearly visible in the cumulative stripping visualization shown in Figure 16. We focus now on an analysis based on the thermal gradient boreholes (see F3, F4.1, F4.2, F4.3, F5, F6.2, and F6.3 in Figures 1 and 14) in the northern part of the study area. The stripping of the gravity effect of the upper and lower Rupelian was correctly done (according to the qualitative dashed area in Figure 16), whereas the stripping of the effect of the Jurassic, but mainly the Triassic and the Buntsandstein, was not correct. This suggests that for these three geologic units, the combined effect of their thickness and their density value was not correct and, consequently, the resulting (stripped) anomalies do not reflect the observed data. Because the geometry of the structures used in the modeling are well-constrained, our interrogation is then around the density values. From that observation, we change the density contrasts according to the reference density to improve the stripping results.
The first test used borehole density values (Figure 17a) measured in the GRT-1 geothermal borehole in 2012 at Rittershoffen, except for the upper Rupelian and the upper part of the lower Rupelian where average values were taken because no borehole data were available. The stripping curves (Figure 17a) show that they can be logically subdivided into three parts: (1) part B where the stripping effect seems overestimated because the stripped values reach zero only after removing the upper and lower Rupelian Formations; (2) part A where the obtained stripped values are smaller than those obtained before, which means that the used density contrast (positive here) should be reversed and then a negative contrast should be used to take into account the gravity effect of the Jurassic, the Trias, and the Buntsandstein units; and (3) part C where the obtained stripped values are relatively close to zero, and so the effect of the basement is well-assessed.
Two other density tests are shown in Figure 17b and 17c, which aim to test the reference density. We used the same density values except for the upper Rupelian for which it was increased by 100 kg:m −3 in Figure 17c to keep the density contrast at 200 kg:m −3 . We can observe that for both tests, the stripped values for the upper and lower Rupelian (only the upper Rupelian in Figure 17b) were reduced compared to the Bouguer anomaly, which is visible in the surface values in Figure 17. Part A for both tests (Figure 17b and 17c) shows that the density contrast should be inverted and then negative density contrasts should be considered. Part C shows that the stripped values obtained  Figure 14 at the location of the thermal gradient boreholes named F3, F4.1, F4.2, F4.3, F5, F6.2, and F6.3 shown on the same figure (see the section below for these boreholes). The gray area shows a qualitative zone that could follow during the stripping approach to reduce the Bouguer anomaly because the aim of the stripping step is to reduce the Bouguer anomaly after each step. (b) The depth and thickness of the geologic units used in the 3D geologic model. We can mainly infer that the gravity response, or the stripped gravity values, is not necessarily correlated to the thickness of the geologic units, for instance, under wells F3 and F4.1.
in Figure 17c are well-corrected, whereas the stripped values obtained in Figure 17b remain unsatisfactory. This means that the used density contrast in part C is correct in sign (i.e., it should be positive) but not in amplitude because it should be higher.
From these three tests and others (not shown), we vary the densities and propose the results shown in Figure 17d. The objective is to reduce the real Bouguer anomaly at each iteration by changing the density values. A single iteration represents one stripping step, for instance, removing the upper Rupelian Formation or removing the upper and lower Rupelian with Jurassic and Triassic are two different iterations.
The stripped curves recovered in Figure 17d show that from the surface to Buntsandstein, the density contrast should be negative rather than positive. Three (relatively) important negative density contrasts were identified: (1) For the upper Rupelian with its −150 kg:m −3 , this negative density contrast is well-observed in all of the tests but also in the sampling measurements and in previous work (e.g., Baillieux et al., 2014). This negative density contrast can be considered as realistic and its gravity effect well-recovered.
(2) The second observed negative density contrast of −100 kg:m −3 is recovered at Jurassic formation. A negative density contrast is also observed in the borehole data between Jurassic and Trias, which is the reverse of what is observed here using stripping (Figure 17a and 17d). This negative density contrast reveals that the positive gravity effect of the Trias, especially for Buntsandstein observed at the center of the study area located between Soultz-Surbourg-Rittershoffen-Trimbach (Figure 13e and 13g), should be negative. Its cumulated gravity effect could be approximately þ6 to þ7 mGal (see Figure 13e and 13g). (3) The third negative density contrast of −100 kg:m −3 is located in the Buntsandstein Formation. This is required in all the tests (see the arrows in Figure 17), which stipulates negative density contrast to reduce de Bouguer anomaly; otherwise, the resulting stripped anomaly becomes even more negative (Figures 14 and 16).
The cumulative stripping approach suggested that the density contrast, mainly from Jurassic down to the bottom of basement, should be changed to reduce the observed Bouguer anomaly (Figure 17d). We now take these density values (Table 3) and apply them to the cross section shown in Figure 15e, achieving sensitivity analysis. The computed response (see the dashed gray line in Figure 18) shows a better fit than an obtained residual anomaly using the mean densities, reported in Table 1 (the gray continuous line of Figure 18). We note that the mismatch between the observed and computed anomalies was reduced when using the densities proposed by the stripping, that is, Table 3. Moreover, the positive anomalies observed in the computed anomaly when using the mean densities (those of Table 1), for instance, at an x-distance of 1012.8 and 1017.5 km, are not yet visible. The computed response shows a smooth curve like that of the observed anomaly, whereas the amplitude remains underestimated. A modeling test with 10% of the fracture-related porosity shows a reduced mismatch with the observed anomaly, although the location of the fractured zone should be moved Figure 17. Stripped gravity anomalies according to the depth/thickness extracted at the same location of the seven gradient boreholes . These curves were obtained for several density values: (a) using borehole values obtained in GRT-1 geothermal borehole (e.g., Baujard et al., 2017), (b and c) using average density values according to different density reference, and (d) improved borehole density values (kg:m −3 ) aimed to reduce the stripped values to reach at the end minimum values around zero. Note that the densities are given with respect to a reference density and can be changed. to the east by approximately 2500 m to obtain the best misfit (compare the green and black lines of Figure 18). We also note that the maximum modeled amplitude is higher than the observed maximum amplitude by approximately 0.2 mGal, whereas the maximum modeled amplitude using 10% of the fracture-related porosity without taking into account the Jurassic is lower than the maximum observed amplitude by only approximately 0.05 mGal (Figure 18). This suggests that the fractured area within the Jurassic is located only in its lower part. We can state that the density values recovered by the stripping at the northern part under the (geo-) thermal boreholes better explains the observations despite the poorly located negative anomaly, for instance, that crossing GRT-1 and 2 is not explained only by the geometry and density values; high fracturing could be the reason.
We return to the relationship between the (negative) gravity anomalies and the thermal gradient boreholes. The temperature profiles in the deep geothermal wells drilled in the URG, and in particular in Northern Alsace, show a linear temperature curve in the sedimentary cover down to the top of the Muschelkalk Formation, that is, Trias. Below this depth, the gradient decreases sharply, indicating the transition between a conductive gradient in the upper part of the formation and a convective reservoir below. The top of the Muschelkalk is defined as the caprock of the convection loop of the geothermal brine, which produces a large geothermal gradient of up to 10°C/100 m in this region (Vallier et al., 2018). From this observation, a new exploration method emerged, which consists of drilling shallow wells in the sedimentary part of the graben to measure a temperature profile calculating the gradient in the sedimentary cover and extrapolating the temperature at the caprock and in the geothermal reservoir. This method of exploration was found to be relevant because of its valuable contribution to the estimation of the temperatures at the geothermal target.
A comparative analysis of all of the investigated zones shows strong heterogeneity in the measured thermal gradient values from one zone to another. The gradients range from 6.3°C/100 m to 7.6°C/100 m. The analysis of the thermal values obtained from the gradient boreholes shows that the wells F3 and F5 have a larger gradient than the others . The gravity values show, for these same wells, residual negative values even after stripping (see the dashed circle in Figure 17): This feature was obtained in the examples shown in Figure 17 but also in other examples not shown here. This means that independent of the density values used for modeling, these two wells show lower density in the lower part of the sediments, mainly in the Triassic and Buntsandstein Formations. The decrease in the density values could be related to petrological heterogeneity as revealed by magnetic data analysis (Edel and Fluck, 1989) or by the presence of a highly fractured area with porosity large enough to significantly decrease the bulk density, as shown above in the fracture-related porosity analysis.
From the density, the areas showing high geothermal potential are depicted in Figure 12d and 12e. This potential can be compared to that obtained at the Rittershoffen geothermal project where the area presents a thick sedimentary column (which is the case of the highly negative anomaly shown in Figure 12e at Rittershoffen and east of Seebach southward until Trimbach), but it can also be compared to that obtained at Soultz-sous-Forêts where the area presents a thin sedimentary column (which is the case for the high negative anomaly observed in Figure 12d around Soultz and northward until well F3). These three locally strong negative anomalies ( Figure 12d) compared with those shown in Figure 11f could concern a bigger zone than the one limited by the isovalue of −4 mGal. This is probably true for the zones under the F3 and F5 wells, but it can also reach the F6.2 and F6.3 wells. These last two wells did not show a density deficit like those of F3 and F5. Probably under F6 (2 and 3) the density effect is rather explained by the geometric variation (see the theoretical gravity response around these wells in Figures 9 and 10). To accurately discriminate the best potential between the wells F6.2 and F6.3, additional high-density gravity measurements should be done at this specific area, which could help to better discriminate between these two wells.

CONCLUSION
We presented the results obtained from a new gravity data acquired in northern Alsace. We obtained a Bouguer anomaly using a new data processing approach, which shows small data uncertainties (<0.06 mGal) compared to the older Bouguer anomaly. The western edge of the URG is also accurately delineated compared to the existing one. The mismatch between the new Bouguer anomaly and the older one is located mostly in the areas of significant topography gradient, which indicates our data processing improvement, mainly the topography correction.
We performed a qualitative analysis using mainly a pseudotomography approach, which shows several areas of negative anomaly stipulating the decrease of the bulk density at the depth  of geothermal interest. The comparison between the obtained results of the multiangle horizontal and vertical analysis with the geologic known faults, located mainly at the top of the basement, leads interestingly to a good correlation, but the horizontal resolution is not enough to be able to draw the location of the faults directly from the gravity results. in addition, and as the normal faults showing western (e.g., the Rittershoffen, Soultz, and Kutzenhausen faults) and eastern (e.g., the Hermerswiller fault and other faults located east of the Rittershoffen fault) dips in the same area, this layout somewhat complicates the interpretation of the horizontal derivatives in term of faulting characterization because the gravity response is basically horizontally shifted according to the surface projection of the top of the fault.
We also achieved a quantitative study on the basis on the existing 3D geologic model derived from vintage seismics that shows a high discrepancy with the observed residual anomalies when using the density values measured in boreholes. The discrepancy between the observed and computed Bouguer anomalies reveals areas with lowdensity values. This density decrease could be explained either by the variation of petrography within the basement and/or highly fractured zones associated with geothermal fluid affecting the bulk density values around the known geothermal sites of Soultz-sous-Forêts and Rittershoffen.
We presented also the fracture-related porosity modeling effect on the density values and then on the gravity anomaly. The resulting Bouguer anomaly was significantly affected, and this effect is clearly measurable. The Bouguer anomaly obtained with a 10% fracture-related porosity better explains the observed negative Bouguer anomaly, which decreases the computed anomaly by a maximum value of approximately 2 mGal. Undoubtedly, the fracturerelated porosity between 8% and 10% should be taken into account to explain the negative observed anomalies, for instance, around Rittershoffen, Soultz and northward, Trimbach and northward, and others.
The interpretation of the stripped Bouguer anomaly showed that the density values of the Jurassic but especially for the Triassic and Buntsandstein were overestimated even using the density values measured in the GRT-1 borehole. This means that the borehole density values do not reflect the horizontal density variations, which occur on larger scales. To reduce the Bouguer anomaly, the stripping reveals that a negative density contrast must be taken into account in the Buntsandstein part, which overlays the basement. This means that the part located between the Buntsandstein and the upper part of the basement presents a low density compared with the reference density, which is not necessarily expected.
We also observed this feature when analyzing the stripped values according to thickness and depth. The negative density contrast of −100 kg:m −3 , at least, should be considered to take into account the effect of the Buntsandstein. We revealed through the stripping approach that the lower part of the sediment, mainly the Triassic and Buntsandstein, presents a negative density contrast, which is not observed by borehole measurements. Interestingly, we found a good correlation between the gravity analyses and the thermal gradient boreholes done in the same study area. Under boreholes F3 and F5, gravity interpretation suggests huge density decreases in the Buntsandstein (maybe also in the upper part of the basement), which may arise from high-density fracturing and important geothermal water affecting significantly the bulk density. Analysis of the thermal borehole data also suggests that these two boreholes show a high geothermal potential compared with the other boreholes. This hypothesis will be confirmed only after achieving a deep geothermal well, which may be done in the next couple of years.

ACKNOWLEDGMENTS
This work was conducted in the framework of the EGS Alsace research project, which was co-funded by ADEME (French Agency for Environment and Energy), Electricité de Strasbourg, and the University of Strasbourg. This study is also financially supported by the COGEOS project and LABEX "G-EAU-THERMIE PRO-FONDE." The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under grant ANR-15-CE06-0014-04 called Cantare Alsace. Thanks go to P. Baillieux for the fruitful discussion and for providing us the fault locations. Thanks also go to G. Marquis for helping us to improve the English language.

DATA AND MATERIALS AVAILABILITY
Data associated with this research are available and can be obtained by contacting the corresponding author.

APPENDIX A GRAVITY DATA PROCESSING
Before data processing, the preprocessing step was done using PyGrav (Hector and Hinderer, 2016). The objective is selecting the acceptable values from the measured data. The selection was done for each point, where only accurate values were chosen and kept. This step is very important because it is the quality control step. The instrumental drift was also removed during this step, and a standard deviation was obtained for each point.
Several data processing steps are required before data analysis and interpretation (e.g., Parasnis, 1986). The most critical step is certainly the topography correction, where the gravity effect of the surrounding topography should be assessed and removed from the observed values (e.g., Hinze et al., 2005). This step was done using PGraviFor3D software (Y. Abdelfettah, personal communication, 2017). This software provides several benefits and allows for accurate computation of the topography and the Bouguer plateau in the same run.
To accurately compute the topography effect affecting the observed data, PGraviFor3D uses two DEMs with different resolutions: one having high resolution around the measured point (<10 m) and the second one having a lower (but not bad!) resolution elsewhere (between 30 and 100 m). The gravity effect of the topography until a distance of 1.5°(approximately 167 km) from the gravity station was assessed. Basically, the whole area is subdivided by the software in four zones: (1) inner, (2) near, (3) intermediate, and (4) remote zones. The predominant effect is that generated first by the inner zone and secondarily by the near zone, but the gravity effects of the intermediate and remote zones are not negligible.
In practice, in the inner zone, a 2 m high-resolution DEM was used (CIGAL, 2011). In this zone, the gravity effect is computed by converting the DEM to a finite element mesh (e.g., Abdelfettah et al., 2014). The radius of this zone is 200 m around the gravity station. This kind of meshing allows for accurately approaching the real topography and avoiding under or overestimation using, for instance, the prism assumption. In the near zone, the free DEM of G54 approximately 30 m provided by NASA (SRTM) was used (Farr et al., 2007). The exact prism formula (e.g., Parasnis, 1986) was used in this zone to compute its gravity effect. This zone was extended more than 20 km around the gravity point. In the intermediate zone, a new DEM of cell size of 500 m is generated by the code where the mean elevation is affected to each prism and the gravity effect is also computed using the exact prism formulas. This zone is generally extended until 80 km. The remote zone was extended until 167 km away from the station. The gravity effect was computed by generating a new DEM of 1 km cell size. A geometric parameter exists in the code taking into account the elliptical shape of the earth, allowing the computation of the spherical Bouguer anomaly.