Keywords: Bottom-hole temperatures, kriging, Alsace, Soultz-sous-Forêts, hydrothermal convection, numerical modelling

Geothermal anomalies in sedimentary basins are strongly controlled by fluid circulation within permeable zones. Exploration of new targets requires the understanding of how and why hydrothermal circulation patterns establish in a particular zone. This study presents a new compilation of newly corrected bottom-hole temperature data in the French part of the Upper Rhine Graben, where the Soultz-sous-Forêts temperature anomaly is locked beneath a local horst structure. After a geostatistically constrained interpolation procedure, maps and crosssections are extracted from the 3D thermal block, together with the associated standard deviations. Thermal anomalies are preferentially associated with the thickest zones of the permeable fractured Buntsandstein (sandstones) formation, in apparent contradiction with previous models where two major fault zones were suggested to control fluid flow. The underlying fractured granitic basement hosts fluid circulation patterns which are apparently controlled at large-scale by the inclined basement-sediments interface. Based on these observations, numerical models of hydrothermal convection including an inclined basementsediments interface, a local horst structure, and realistic petrophysical properties have been carried out. The depth-decrease of permeability, the inclination of the interface and the fixed heat flow condition at the base of the model, explain why only a few upwellings can be triggered. Thermal anomalies and a measured temperature profile can be reproduced when fault permeability equals 10 -14 m². Interestingly, structure of convective patterns also exhibits a second and hotter upwelling, in the Rittershoffen area, 8 km east of the Soultz-sous-Forêts upwelling zone, where another geothermal exploration project is now underway. The understanding of hydrothermal convection with realistic fluid and rock properties clearly appears as a predictive tool for geothermal exploration strategies.

A C C E P T E D

M A N U S C R I P T

Introduction

Exploration of geothermal energy systems consists in the identification of hot reservoirs where fluids circulate naturally, and where best economical conditions are met. Apart from the required socioeconomic studies, characterization of a geothermal reservoir can be achieved with geological studies, geochemical analyses and geophysical field surveys (see the recent review of various methods by [START_REF] Bruhn | Exploration methods[END_REF]. In addition, the understanding of fluid dynamics in anomalously hot fractured area could also help in the identification of a potential reservoir [START_REF] Mckenna | Numerical modelling of transient Basin and Range extensional geothermal systems[END_REF][START_REF] Oliver | Numerical models of extensional deformation, heat transfer, and fluid flow across basement-cover interfaces during basin-related mineralization[END_REF][START_REF] Kuhn | Numerical investigation of the effect of heterogeneous permeability distributions on free convection in the hydrothermal system at Mount Isa, Australia[END_REF]Ingebritsen et al., 2010).

Although a number of exploration methods can be used to delineate favourable geological settings in terms of fluid circulation and thus related permeability (e.g. [START_REF] Muñoz | Exploring the Gross Schönebeck (Germany) geothermal site using a statistical joint interpretation of magnetotelluric and seismic tomography models[END_REF][START_REF] Bujakowski | Integrated seismic and magnetotelluric exploration of the Skierniewice, Poland, geothermal test site[END_REF][START_REF] Jousset | Hengill geothermal volcanic complex (Iceland) characterized by integrated geophysical observations[END_REF], characterization of anomalously hot zones at a given depth is very challenging. While permeability contrasts may locally reach several orders of magnitude [START_REF] Clauser | Permeability of crystalline rocks[END_REF][START_REF] Manning | Permeability of the continental crust: implications of geothermal data and metamorphic systems[END_REF], the required thermal anomalies may differ from the surroundings by only a few tens of °C. Actually, temperature in the first 5 kilometres of the subsurface is not easily predictable since equilibrium temperatures are only measured in basement rocks, through mining boreholes whose depth rarely exceeds 1500 m. In sedimentary basins, deeper temperature measurements are available through the numerous petroleum boreholes, but these data are strongly affected by transient disturbances. The estimate of equilibrium temperatures from these measurements is not straightforward [START_REF] Goutorbe | Comparison of several BHT correction methods: a case study of an Australian data set[END_REF] since it requires in particular scrupulous examination of each temperature data.

Previous studies demonstrated that in the subsurface, lateral temperature variations of several tens of °C can be found within a distance of a few kilometres [START_REF] Pribnow | Thermal tracking of Upper Crustal Fluid Flow in the Rhine Graben[END_REF][START_REF] Rühaak | Detecting thermal anomalies within the Molasse Basin, southern Germany[END_REF][START_REF] Bonté | Subsurface temperature maps in French sedimentary basins : new data compilation and interpolation[END_REF][START_REF] Garibaldi | Thermal anomalies and geological structures in the Provence basin : implications for hydrothermal circulations at depth[END_REF][START_REF] Guillou-Frottier | Heat flow and deep temperatures in the Southeast Basin of France: implications for local rheological contrasts[END_REF]. Moreover, downward extrapolation of subsurface temperature measurements through the use of an averaged regional temperature gradient may result in significant errors and could thus induce wrong decisions in geothermal exploration strategies.

A C C E P T E D M A N U S C R I P T [START_REF] Genter | Typology of potential Hot Fractured Rock resources in Europe[END_REF] analysed a map of presumed temperatures at 5 km depth in Europe and demonstrated that a number of potentially interesting zones were not confirmed by published heat flow data. In addition, interpolation of temperature data is seldom performed with adapted geostatistical analysis before elaborating temperature maps. Previous atlases of subsurface temperature in Europe [START_REF] Haenel | Atlas of subsurface temperatures in the European Community[END_REF][START_REF] Hurter | Atlas of geothermal resources in Europe[END_REF] contain numerous maps at different depths, but in most cases, essential information is not available (origin of data, method of correction, interpolation procedure, etc). Consequently, it is often delicate to work with existing temperature maps of the subsurface, since they contain temperature data which are not homogeneously corrected, and because resulting uncertainties at a few kilometres depths may easily reach 20-30 °C (e.g. [START_REF] Förster | Analysis of borehole temperature data in the Northeast German Basin: continuous logs versus bottom-hole temperatures[END_REF][START_REF] Goutorbe | Comparison of several BHT correction methods: a case study of an Australian data set[END_REF].

Nevertheless, the exploration of deep geothermal resources is still guided by the identification of hot and sufficiently permeable zones, at depths generally lower than 5 kilometres [START_REF] Fritz | On the way to the exploitation of deep geothermal resources in naturally fractured environments[END_REF]. Contrarily to natural permeability, temperature cannot be artificially increased, and the need of hot temperatures may require deeper drilling operations, involving additional costs. The Soultz-sous-Forêts geothermal project, located in the French Part of the Upper Rhine Graben (FPURG) began in the late 80's, and was subjected to such problematic choices (http://english.geothermie-soultz.fr/). At the beginning of the project, tens of shallow petroleum wells in the sedimentary basin had shown that temperature above 110 °C could be found at 1000 m depth, and at that time, temperatures of 150-170 °C were expected at depths lower than 2000 m [START_REF] Gerard | Proceedings of the first EEC/US workshop on geothermal Hot Dry Rocks technology[END_REF].

Actually, it appears today that temperature of 160 °C is only reached at 3600 m depth, because temperature gradient drastically decreases below 1000 m depth within a permeable sandstone formation, before being extremely low (~ 5 °C.km -1 ) within the permeable granitic basement. Below 3500 m depth, where the granitic basement is poorly permeable, a normal temperature gradient of 30 °C.km -1 is recovered. A more detailed description of the 5 kmdepth temperature profile measured in the GPK2 borehole can be found in [START_REF] Genter | Contribution of the exploration of deep crystalline fractured reservoir of Soultz to the knowledge of Enhanced Geothermal Systems (EGS)[END_REF].

A C C E P T E D M A N U S C R I P T

Such decrease in temperature gradient is typical of a convective signature created by a hot upwelling (or a cold downwelling) fluid (e.g. [START_REF] Benderitter | Calcul de l'effet thermique d'une remontée hydrothermale dans un socle fracture. Application à l'anomalie géothermique de Soultz-sous-Forêts (Nord Alsace)[END_REF]. Actually, deviations of subsurface temperatures from a regular increase can be explained by a number of heat transfer processes, including heat refraction effects due to thermal conductivity contrasts [START_REF] Guillou-Frottier | High heat flow in the Trans-Hudson Orogen, central Canadian Shield[END_REF], thermal convection within thick porous layers [START_REF] Simms | Thermal convection in faulted extensional sedimentary basins: theoretical results from finite-element modelling[END_REF] or anomalous heat flow coming from high heat producing rocks [START_REF] Mclaren | Long-term consequences of the redistribution of heat producing elements within the continental crust: Australian examples[END_REF]. In sedimentary layers, anisotropy of thermal conductivity may also lead to anomalous heat transfer processes resulting in thermal anomalies [START_REF] Vasseur | Convective and conductive heat transfer in sedimentary basins[END_REF].

Presence of evaporites or salty layers in sedimentary basins may also create heat refraction processes leading to large temperature variations (case of the south-east basin in France, [START_REF] Guillou-Frottier | Heat flow and deep temperatures in the Southeast Basin of France: implications for local rheological contrasts[END_REF]. Around and within faulted zones, fluid circulation may also create lateral temperature differences greater than 20 °C [START_REF] Clauser | Analysis of convective and conductive heat transfer in a sedimentary basin, demonstrated for the Rhinegraben[END_REF]Lopez andSmith, 1995, Fleming et al., 1998;[START_REF] Bächler | Impact of graben-parallel faults on hydrothermal convection -Rhine Graben case study[END_REF][START_REF] Wisian | Numerical modelling of Basin and Range geothermal systems[END_REF][START_REF] Yang | Factors controlling free thermal convection in faults in sedimentary basins: implications for the formation of zinc-lead mineral deposits[END_REF][START_REF] Appold | Numerical modelling of the origin of calcite mineralization in the Refugio-Carneros fault, Santa Barbara basin[END_REF][START_REF] Harcouët-Menou | Hydrothermal convection in and around mineralized fault zones: insights from twoand three-dimensional numerical modelling applied to the Ashanti belt[END_REF][START_REF] Magri | Deep geothermal groundwater flow in the Seferihisar-Balçova area, Turkey: results from transient numerical simulations of coupled fluid flow and heat transport processes[END_REF][START_REF] Garibaldi | Thermal anomalies and geological structures in the Provence basin : implications for hydrothermal circulations at depth[END_REF]. In other words, the various aspects of heat transfer processes in the shallow crust can lead to significant temperature variations.

The objective of this paper is to present a geostatistically-constrained 3D thermal block located in the French part of the Upper Rhine sedimentary graben, and to discuss the 3D geometry of the Soultz-sous-Forêts thermal anomaly in terms of fluid circulation. Starting with an updated and revisited database of temperature data coming from petroleum boreholes of the sedimentary area, a correction procedure is presented and applied.

Geostatistical analysis of the database provides error estimates for the entire thermal block.

Extracted maps and cross-sections can thus be illustrated with standard deviations related to the interpolation procedure. With a three-dimensional view of the present-day thermal regime of the area, some possible fluid circulation patterns at the basin scale are inferred. In order to better understand the suggested structure of hydrothermal convection in the FPURG, numerical simulations of heat and fluid transfer in the area, constrained by geological data and petrophysical properties, have been performed.

A C C E P T E D M A N U S C R I P T

Geology and geodynamics

Geological setting

The Upper Rhine Graben (URG) corresponds to one branch of the European Cenozoic Rift System (ECRIS), which has formed in the foreland of the Alps and Pyrenees [START_REF] Ziegler | Crustal evolution of Western and Central Europe[END_REF], as shown in Figure 1a. The URG started to form during Late Eocene in response to NNE trending Alpine compression. In the area of Soultz-sous-Forêts, the faulted granitic basement is overlain by 1400 m of Cenozoic, Mesozoic and Permian sediments. The rifting reached a climax in the Early Oligocene under WNW-ESE orientated extension, reactivating some Variscan structures. Moho topography shows an important rise up to 24 km depth in the southern URG [START_REF] Dèzes | Evolution of the European Cenozoic Rift System: interaction of the Alpine and Pyrenean orogens with their foreland lithosphere[END_REF]. The URG's deep crustal structure, which is likely to be related to mantle uplift, was investigated by deep seismic surveys [START_REF] Brun | Crustal scale structure of the southern Rhinegraben from ECORS-DEKORP seismic reflection data[END_REF].

The present-day stress field is a compressive event oriented NW-SE. Locally, based on borehole data, it has been shown that the present stress field characterizes a normal regime in the sediments and in the granite down to 3500 m depth [START_REF] Valley | The relation between natural fracturing and stress heterogeneities in deepseated crystalline rocks at Soultz-sous-Forêts (France)[END_REF]. Deeper, the magnitudes of the maximal horizontal stress and of the vertical stress (overburden) are very close, indicating a possible strike-slip regime below 3500 m [START_REF] Valley | The relation between natural fracturing and stress heterogeneities in deepseated crystalline rocks at Soultz-sous-Forêts (France)[END_REF].

In transverse cross-section around the Soultz-sous-Forêts area, the URG shows a strong dissymmetrical shape with a steeply sloping basement-sediments interface on the eastern side by comparison with a low angle one on the western side. Figure 1a illustrates thickness variations of the sedimentary cover within the URG (after [START_REF] Rotstein | The Upper Rhine Graben (URG) revisited: Miocene transtension and transpression account for the observed first-order structures[END_REF], indicating that the sediments-basement interface never represents a flat boundary.

Local geology at Soultz-sous-Forêts

The geothermal Soultz-sous-Forêts site is located on a local horst structure which is bounded by large-scale faults showing significant apparent vertical off-set (Figure 1b). Recent reprocessing of seismic reflection profiles highlighted the complex fault patterns in the horst

Physical and geological data

The physical parameters and geological data used in this study (thermal conductivity, permeability, heat production, formation thickness) were estimated from geophysical investigations made on different geological units (sediments, granite), or from laboratory measurements on rock and core samples.

In the Soultz-sous-Forêts area, a typical sedimentary sequence between the surface and 600 m corresponds to Tertiary soft sediments made of dolomitic marls, clayey marls, conglomerates marls with evaporites, and marls with intercalations of oily sand and sandstones. Those Tertiary layers are reputed to exhibit low permeability. Below 600 m depth, Tertiary formations directly overlay Jurassic limestone and Liassic marls. For this section (0-800 m), thermal conductivity ranges between 1.3 and 1.4 W.m -1 .K -1 (Flores- , 1992). Between 800 and 1000 m depth, Upper and Middle Triassic layers made of evaporites (Keuper) and dolomites (Muschelkalk) show thermal conductivity of 2.1 W.m -1 .K -1 .

Marquez

Deeper, the Permo-Triassic sandstones from Buntsandstein correspond to a thick porous unit of 400 m depth. They are naturally fractured with small-scale fractures associated with localized large-scale faults detected during drilling operations which are either sealed by mineral deposition or locally permeable [START_REF] Genter | Fracture analysis and reservoir characterization of the granitic basement in the HDR Soultz project (France)[END_REF]. Their thermal conductivity is 2.5 [START_REF] Haffen | Petrophysical characteristics of sandstones dating from the Buntsandstein in the Upper Rhine Graben: case of the borehole EPS1 (Soultz-sous-Forêts, France)[END_REF] and permeability of fracture zone is supposed to be around 5.10 -15 m².

A C C E P T E D M A N U S C R I P T W.m -1 .K -1
Within the nearby Landau geothermal field, permeability of Muschelkalk-Buntsandstein formations, fault and granite matrix were assumed to be of 10 -16 , 5 10 -13 and 10 -18 m 2 , respectively [START_REF] Bächler | Impact of graben-parallel faults on hydrothermal convection -Rhine Graben case study[END_REF]. Permian sandstone deposition show strong local variation in terms of thickness (from 10 to 100 m). However, these sandstones showing a lithology close to that of Buntsandstein formation, there is no reason to affect them a different permeability.

A number of studies on the Soultz-sous-Forêts geology show that the top basement (1400-2100 m) is much more fractured than the overlying Buntsandstein [START_REF] Genter | Contribution of the exploration of deep crystalline fractured reservoir of Soultz to the knowledge of Enhanced Geothermal Systems (EGS)[END_REF]. The deep basement is highly fractured and hydrothermally altered. All the natural fractures are systematically sealed with hydrothermal deposits, and only a minor part of them is naturally permeable [START_REF] Genter | Analysis of macroscopic fractures in granite in the HDR geothermal well EPS-1, Soultz-sous-Forêts, France[END_REF]. Linear fracture density based on core analysis in the Buntsandstein (1000-1400 m) and in the top granite (1400-2200 m) is 1.0 and 3.4 fracture.m -1 respectively [START_REF] Genter | Comparative analysis of direct (core) and indirect (borehole imaging tools) collection of fracture data in the Hot Dry Rock Soultz reservoir (France)[END_REF]. Moreover, detailed fracture analysis demonstrated that the section between 1400 and 2100 m is much more fractured than the deepest part [START_REF] Genter | Comparative analysis of direct (core) and indirect (borehole imaging tools) collection of fracture data in the Hot Dry Rock Soultz reservoir (France)[END_REF]. Based on Soultz granite sample measurements in laboratory under extrapolated in-situ pressure conditions, [START_REF] Hettkamp | Hydraulic properties in the Rhine Graben basement material[END_REF] have shown that matrix permeability and fracture permeability range from 10 -18 to 10 -20 m², and 10 -15 to 10 -17 , respectively. Permeability values derived from pressure oscillation tests in the Soultz wells having a 700 m open-hole section, range between 10 -14 and 10 -15 m² [START_REF] Hettkamp | Hydraulic properties in the Rhine Graben basement material[END_REF].

The porphyritic monzogranite has a high heat production rate decreasing from 6 µW.m -3 at 1.5 km to 2.7 µW.m -3 at 4.5 km [START_REF] Pribnow | Heat-and fluid-flow at the Soultz Hot Dry Rock system in the Rhine Graben[END_REF] compared to the value of 0.8-1.2 µW.m -3 for the sediments [START_REF] Clauser | Analysis of convective and conductive heat transfer in a sedimentary basin, demonstrated for the Rhinegraben[END_REF].

The heat flow anomaly

A C C E P T E D M A N U S C R I P T

During the last 30 Myrs, thermal and mechanical events have been affecting deformation of the European lithosphere [START_REF] Cloetingh | Lithosphere tectonics and thermo-mechanical properties: an integrated approach for Enhanced Geothermal Systems exploration in Europe[END_REF]. In France, undulations of the Moho can be attributed to response to compressive Alpine and Pyrenean events [START_REF] Lefort | Topography of the Moho undulations in France from gravity data: their age and origin[END_REF]. Thermal anomalies, such as those measured in the French Massif

Central and in the Rhine graben areas, may be related to shallow Moho depths. In the Rhine graben area, mantle upwelling and thinning of the lithosphere could have occurred since high heat flow values (greater than 110 mW.m -²; [START_REF] Lucazeau | Heat flow density data from France and surrounding margins[END_REF]) are presently recorded, some 30 Myrs after the beginning of the graben formation. Actually, any thermal event that occurred 30 Myrs ago at 30 km depth, diffuses toward the Earth's surface with a diffusive time constant ~ d²/ (where d is crustal thickness in m, and  is thermal diffusivity, equal to 10 -6 m².s -1 ). Hence, one part of the present-day high heat flow values may indeed correspond to thermal events having occurred 30 Myrs ago. However, the scale of lateral variations of surface heat flow values (a few kilometres) is much smaller than the crustalscale of thermal events leading to continental rifting (several tens of kilometres).

Consequently, it is probable that anomalously high temperature gradients in the subsurface are better explained by local geological features (e.g. [START_REF] Mareschal | Heat flow and deep thermal structure near southern edge of the Canadian Shield[END_REF]. This does not exclude the presence of an anomalously high mantle heat flow, as it was suggested for the French Massif Central [START_REF] Lucazeau | Interpretation of heat flow data in the French Massif Central[END_REF].

The high surface heat flow measured in the Rhine graben has been interpreted by [START_REF] Royer | Steady-state geothermal model of the crust and the problem of the boundary conditions: application to a rift system, the southern Rhinegraben[END_REF] as the result of steady-state heat refraction effects between distinct crustal units together with asymmetrical lithospheric heating. According to [START_REF] Clauser | Analysis of convective and conductive heat transfer in a sedimentary basin, demonstrated for the Rhinegraben[END_REF], it is a basin-scale fluid circulation pattern which can locally increase the background basal (below sediments) heat flow of 80 mW.m -² by more than 50%. Their convective models assume that fluid circulation predominantly occurs within two permeable fault zones (the Rhine and the Schwarzwald fault zones). Models by [START_REF] Kohl | Steps towards a comprehensive thermo-hydraulic analysis of the HDR test site Soultz-sous-Forêts[END_REF] explained the Soultz-sous-Forêts temperature profile with fluid circulation focused within a restricted permeable area (including the Kutzenhausen and the Soultz faults) in order to get the smallest size of a possible convective cell. Other numerical models of hydrothermal

A C C E P T E D M A N U S C R I P T
convection within the entire Rhine graben illustrate different convective patterns, where several convective cells develop within the permeable layers [START_REF] Carlier | Convective heat transfer at the Soultz-sous-Forêts geothermal site: implications for oil potential[END_REF].

However, according to recent permeability measurements [START_REF] Haffen | Petrophysical characteristics of sandstones dating from the Buntsandstein in the Upper Rhine Graben: case of the borehole EPS1 (Soultz-sous-Forêts, France)[END_REF], the sandstone formation overlying the fractured granite would be less permeable than initially thought by [START_REF] Carlier | Convective heat transfer at the Soultz-sous-Forêts geothermal site: implications for oil potential[END_REF] see also previous section. [START_REF] Clauser | Decoupled thermal and mantle helium anomalies : implications for the transport regime in continental rift zones[END_REF] developed additional models where, as suggested by [START_REF] Kohl | Steps towards a comprehensive thermo-hydraulic analysis of the HDR test site Soultz-sous-Forêts[END_REF], the role of fault permeability is of major importance. While topography of rift shoulders drives downward flow, thermal convection (buoyancy-driven flow) within permeable fault zones appears as a natural groundwater flow system. According to these studies, the Soultz-sous-Forêts thermal anomaly is clearly controlled by free convection within fault zones (see also a similar study for the Landau geothermal site by [START_REF] Bächler | Impact of graben-parallel faults on hydrothermal convection -Rhine Graben case study[END_REF], or two other ones on the Basin and Range geothermal system, by [START_REF] Mckenna | Numerical modelling of transient Basin and Range extensional geothermal systems[END_REF] and by [START_REF] Wisian | Numerical modelling of Basin and Range geothermal systems[END_REF]. However, the relative importance between fault zone permeability and pervasive (Buntsandstein formation and top granitic basement) permeability has not been addressed.

Numerical models presented here will tackle this issue since it is likely that fault zone permeability behaves as a time-dependent parameter [START_REF] Ingebritsen | Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism[END_REF].

The two distinct models of fluid circulation patterns at the graben scale (i.e. one scenario with several convective cells, and thus several hot upwellings, and one scenario with locally fault-controlled free convection, and thus a single upwelling zone) strongly depend on the chosen permeability structure as well as on thermal conditions beneath the graben. Before tackling this discussion, temperature data in the FPURG have been revisited

and gathered in a 3D thermal block, allowing visualization of anomalously hot and cold areas in the FPURG.

Data compilation and building of the 3D thermal block

The objective here is to build a 3D thermal model with an estimation of interpolation error, using geostatistical techniques. Compilation of all available temperature data and procedure for correcting them from transient disturbances are described below.

A C C E P T E D M A N U S C R I P T

Building of the database

Previous temperature maps of subsurface temperatures in the FPURG (Figure 2) were based on a temperature database where BHT (bottom-hole temperature) measurements were corrected with a purely statistical method [START_REF] Haenel | Atlas of subsurface temperatures in the European Community[END_REF][START_REF] Le Masne | Aquifères profonds d'Alsace : constitution d'une base de données à usage géothermique[END_REF].

More precisely, this statistical correction (i) was regionally defined, (ii) did not depend on the measurement history and (iii) was calibrated with available measurements from drill stem tests (DST) [START_REF] Gable | Acquisition et rassemblement de données géothermiques disponibles en France[END_REF]. For example, a correction of +20 °C was applied to all BHT measurements located at 1000m depth. Later, transient perturbations were more correctly accounted for [START_REF] Pribnow | Thermal tracking of Upper Crustal Fluid Flow in the Rhine Graben[END_REF], but the interpolation method was not detailed and error on the interpolation was not described. The recent country-scale temperature maps of [START_REF] Bonté | Subsurface temperature maps in French sedimentary basins : new data compilation and interpolation[END_REF] did not include more than 23 temperature data in the FPURG. Today, with an updated temperature database, and with an adapted correction method, we have significantly improved database quality. During this new compilation, we retrieved additional temperature measurements and temporal data, located on logging headers. In particular, we systematically collected, when available, the time delay between the end of mud circulation and the temperature measurement, in order to make analytical corrections possible [START_REF] Goutorbe | Comparison of several BHT correction methods: a case study of an Australian data set[END_REF].

We gathered four different compilations: a compilation of French BHT data, made by [START_REF] Gable | Acquisition et rassemblement de données géothermiques disponibles en France[END_REF], a later unpublished compilation (F. Lucazeau, pers. comm.) temperature profiles measured at Soultz-sous-Forêts in geothermal boreholes were not inserted in this database, but they will be used in the second part of this study to constrain the fluid circulation models.

As for the [START_REF] Bonté | Subsurface temperature maps in French sedimentary basins : new data compilation and interpolation[END_REF] study, the analytical correction called ICS method (for "instantaneous cylinder source", see [START_REF] Goutorbe | Comparison of several BHT correction methods: a case study of an Australian data set[END_REF] was applied to BHT data when temporal information was available on logging headers. This was the case for 111 BHT data.

The remaining set of BHT data (212 over 323) has been corrected with the AAPG statistical method (American Association of Petroleum Geologists, [START_REF] Bodner | Temperature variations in South Texas subsurface[END_REF][START_REF] Deming | Application of bottom-hole temperature corrections in geothermal studies[END_REF]. As detailed in [START_REF] Bonté | Subsurface temperature maps in French sedimentary basins : new data compilation and interpolation[END_REF], differences between ICS and AAPG correction methods do not exceed 5°C for depths lower than 2000m (i.e. for 95% of our dataset in this study, Figure 3b), and may reach 12 °C at 4000m depth.

Spatial distribution of temperature data shows 2 main regional clusters: one in the Soultz-sous-Forêts area (Pechelbronn oil field), and another one south of Colmar (Potash mines of Alsace). Since gathering of data may affect the interpolation results, this feature will be accounted for in the interpolation procedure. Likewise, the depth distribution shown in Figure 3a will affect interpolation uncertainty more importantly below 1500 m.

Methodology for building the 3D temperature model

As temperature increases with depth, the 3D interpolation of temperature data corresponds to a non-stationary problem for which two main options are possible: (1) use a non-stationary geostatistical model, where the drift (spatial trend) is included in the model and is locally interpolated (universal Kriging or Kriging in the "intrinsic random function of order k" framework theory, [START_REF] Chilès | Geostatistics: Modelling Spatial Uncertainty[END_REF], or (2) use a deterministic model for the drift, like a linear increase of temperature with depth, and work with the residual, the residual being the difference between the temperature and a global drift fitted to data (see Chilès and Delfiner, 1999, p. 164 and p. 233). The first approach is theoretically preferable because (i) the interpolation error includes the error on the drift, and (ii) the drift is fitted locally. The second approach assumes that the drift is fitted globally and is known everywhere with no

A C C E P T E D M A N U S C R I P T
error. However, in practice, the first approach is not straightforward and is difficult to tune: a generalised covariance must be fitted by automatic black box procedures based on trial and errors and cross validation tests, with little or no control of the user. At the interpolation step, determination of the local drift may lead to numerical instabilities. On the other hand the second approach is more intuitive, can be checked visually (fitting of the drift, variogram of the residuals) and can be tuned more easily. After a preliminary test of both options, we decided to choose the second approach.

More precisely, the temperature T(x,y,z) is expressed as :

) , , ( ) , , ( ) , , ( z y x R z y x D z y x T   (1)
where D(x,y,z) is the drift and R(x,y,z) a residual that is supposed to be stationary. In this paper, the drift was modeled as a function of the vertical coordinate (depth z) only, the possible horizontal components D(x,y) of the drift being in practice included in the residual.

Since there is no horizontal drift (regular increase or decrease) in temperature data, this is a reasonable hypothesis. Expression (1) was thus simplified as:

) , , ( ) ( ) , , ( z y x R z D z y x T   (2)
The procedure that was applied to obtain the 3D temperature model is the following:

1) From the initial data T(x,y,z), definition and fitting of the vertical drift. Let D*(z) be the fitted drift 2) Computation of the residuals as R(x,y,z) = T(x,y,z) -D*(z)

3) Building and fitting the variogram of the residuals R(x,y,z) 4) Interpolation of the residuals by kriging. Let R*(x,y,z) be the kriged result 5) Computation of the estimated temperature by T*(x,y,z) = D*(z) + R*(x,y,z)

The fitting geotherm

To estimate the vertical drift D*(z), different increasing trends of temperature with depth were tested (linear trend and polynomial trends fitted by least squares -Figure 3c). As observed in numerous studies, the averaged surface temperature (10 °C in our case) should not be imposed when fitting BHT data since it creates largest values of residual temperatures: a best fitting linear trend is obtained when surface temperature is about 10-15 °C larger than the real surface temperature (see similar cases in [START_REF] Husson | Thermal regime of the NW shelf of the Gulf of Mexico. Part A: thermal and pressure fields[END_REF]or in Bonté et al., 2010). This observation, which is not correlated with low thermal conductivity of shallow formations, might be due to overestimates of BHT correction at shallow levels. Both cases have however been tested, together with a non-linear trend (polynomial trend of degree 4). It appears that deep temperature data are better fitted with the polynomial trend (Figure 3d), since lower standard deviations for residual values are obtained. Note however that correlation coefficients (indicated on Figure 3c) are close together for the second linear trend and the polynomial one. This polynomial trend is therefore chosen as the representative fitting geotherm for the entire FPURG, and deviations from this trend correspond to residual temperatures. Since there is no theoretical criterion to define the best fitting geotherm at the scale of the graben, our choice could be discussed. It is however believed that the inferred anomalies might be underestimated if a locally averaged temperature gradient would be considered.

Declustering

Spatial distribution of temperature data is not homogeneous (Figure 3a), and the building of a representative variogram may be biased by clusters of data in a small area. To avoid the over-representation of some groups of data, a declustering method was applied to the entire set of data (e.g. [START_REF] Bourgault | Spatial declustering weights[END_REF][START_REF] Chilès | Geostatistics: Modelling Spatial Uncertainty[END_REF]. Two sizes of declustering area have been tested, 5 and 10 km. By comparing the average of standard deviations for residual values, the declustering size of 10 km appeared as the best choice. The result of the declustering is a weight that is assigned to each temperature data, and is a function of the inverse of local data density. This weight is then used for the variogram computation [START_REF] Rivoirard | Weighted Variograms[END_REF]. values (R(x i ), R(x i +h)) at distance of h is given by:

Geostatistical description of data: building of the variogram

    ) ( ) ( ) ( 2 1 ) ( 1 2 h N h x R x R h h N i i i       (3)
In expression (3), R(x i ) stands for R(x,y,z) at location i, and h is a 3D distance with 3 components.

The experimental variogram of the residuals was computed along the wells (vertical direction), and along the horizontal direction. Figure 4 shows the variogram obtained in both directions. For the same distance, the value of the variogram is very different along the vertical direction by comparison with the horizontal direction. This depicts a strong anisotropy of the residuals, the variability being higher along the vertical axis. It must be noted that no clear anisotropy was observed in the horizontal plane, this allowing the residuals to be considered as isotropic in this plane.

In order to use the features of the obtained experimental variograms in the interpolation procedure, a variogram model fitting the experimental variogram is required.

Several models were compared, each one including a nugget effect (corresponding to the sum of the local variability and of the variance of temperatures measurement error), and one to three stationary structures defined by their type, their sill and their range. By definition, the range of a structure is the distance beyond which there is no more spatial correlation between measurements, and the sill is the total contribution of the structure along the the variogram axis.

Table 1 describes a series of models that have been tested to fit the experimental variograms, for the three different increasing trends of temperature with depth. For each trend, several models of variogram (lines in Table 1) were tested. Among the different and 1600 m. We did not interpret these ranges physically, though it is possible that the largest horizontal range corresponds to the typical size of thermal anomalies. This point will be discussed later.

Kriging of the residuals and getting the 3D-thermal block

Prior to Kriging (see Chiles and Delfiner, 1999), duplicate data were removed in order to avoid numerical instabilities. This was done using the option "look for duplicates" in the Isatis ® geostatistical tool. When two data were closer than 10 m, only one was retained. After this pre-processing, the final set of data comprised 480 residual values separated by more than 10 m.

A kriging procedure was then applied on the selected 480 residual temperatures, which were calculated with the reference geotherm described in section 3.3 (the polynomial trend). The variogram model of Figure 4 (model C3 in Table 1) was used in the interpolation procedure. The output consists in a 3D thermal block with grid elements of 500 m x 500 m x 50 m, across which maps and cross-sections can be extracted, together with their respective standard deviations.

Figure 5 shows temperature maps of the FPURG at different depth levels, from 500 m to 3000 m below the surface. On each map are represented data points located at ± 500m from the depth level (black plus signs). The largest standard deviations of the interpolations (red colour in smallest maps) are the same for all maps (17 °C), and appear predominant at depths of 2500 m and 3000 m. The smallest interpolation errors (blue colour, 5-6 °C down to 2500 m) are predominant at depths of 1000 m and 1500 m. On temperature maps, the anomalously hot and cold areas (yellow to red, and dark blue zones, respectively) define thermal anomalies with respect to the polynomial trend. In other words, for a given depth, the

A C C E P T E D M A N U S C R I P T
background colour corresponds to the temperature value obtained with the polynomial trend while lower and greater temperature values represent thermal anomalies. It appears that these anomalies are not necessarily continuous with depth, and may become stronger or weaker at depth. Their amplitudes and their spatial evolution are described below.

Thermal anomalies

Definition and description of the anomalies

Figure 6 illustrate temperatures at a depth of 800 m, as it was done with another interpolation procedure by [START_REF] Pribnow | Thermal tracking of Upper Crustal Fluid Flow in the Rhine Graben[END_REF] for the entire URG (Figure 2b).

Thermal anomalies, as defined with this map, have been contoured and labelled from A to H.

The major hot (positive) anomalies are Soultz-sous-Forêts (A), south of Strasbourg (C) and Altkirch (H), and the major cold (negative) anomalies are north of Strasbourg (B) and Colmar-Mulhouse-Sélestat (E). Minor anomalies are also visible: west of Sélestat (D), south-west of Breisach (F), and Wittenheim (G). These anomalies are also visible in vertical cross-sections:

seven west-east cross-sections, numbered from 1 to 7 in Figure 6, and one south-north cross-section (from A to H in Figure 6) are shown in Figure 7. On each vertical cross-section, the red and black lines correspond to the upper limit of the Buntsandstein formation and to the base of Trias (interface between the granitic basement and sediments), respectively.

They were extracted from the maps of [START_REF] Munck | Synthèse géothermique du Fossé Rhénan SupérieurService Géologique Régional Alsace[END_REF].

The major positive thermal anomaly in the FPURG is the "Soultz-sous-Forêts

Thermal anomalies and geology

Contrary to thermal maps at constant depths, thermal cross-sections (Figure 7) help in the identification of possible convective processes (fluid circulation driven by buoyancy).

Whereas large-scale thermal undulations are not necessarily caused by fluid circulation, small-scale uplift of isotherms may reflect the presence of hot upwelling fluids. Section numbered 1 in Figure 7a crosses the Soultz-sous-Forêts anomaly, which looks like a hot bubble rising westwards along the inclined basement. Interestingly, one can note that the hottest zone of this upwelling is located within the permeable fractured Buntsandstein sandstone formation (depth range between the red and black lines), where thickness is important.

A C C E P T E D M A N U S C R I P T

The six other west-east cross-sections shown in Figure 7a illustrate distinct features.

The large-scale cold anomaly at shallow depths in section 2 (anomaly B) could correspond to downwelling fluids in response to the adjacent upwelling zone around Soultz-sous-Forêts.

However, at depth, isotherms around the thickest parts of the Buntsandstein formation (in between red and black lines) are not significantly disturbed. Because the anomaly is predominant in the upper 1500 m, it should better correspond to a recharge zone from cold meteoric fluids, preventing the possible development of a deeper positive anomaly.

On cross-section numbered 3, another large-scale positive anomaly of however less amplitude (anomaly C), is visible between 2000 and 3000 m, and it appears that thickness of the Buntsandstein formation is comparable to that of section 1. On the opposite, sections 4 to 7 do not show thick Buntsandstein formations, and no significant positive anomalies are visible. The small-scale shallow (~ 500 m depth) positive anomaly on section 5 corresponds to anomaly F. Despite few data are present, its width and its sharp geometry could correspond to the thermal signature of recent fluid flow within a fault. Indeed, models of fluid circulations have illustrated such sharp increase of temperature when hot fluids flow within a faulted zone (e.g. [START_REF] Yang | Factors controlling free thermal convection in faults in sedimentary basins: implications for the formation of zinc-lead mineral deposits[END_REF][START_REF] Wisian | Numerical modelling of Basin and Range geothermal systems[END_REF][START_REF] Mclellan | Modelling fluid convection stability in continental faulted rifts with applications to the Taupo Volcanic Zone, New Zealand[END_REF].

The north-south cross-section in Figure 7b shows that at the scale of the graben, and below 1500 m depth, positive anomalies develop more easily where the permeable fractured Buntsandstein formation is sufficiently thick. As it can be seen in Figure 7b, below 1500m, the large negative anomaly south of Strasbourg is of much smaller amplitude than those of positive anomalies to the north. At shallower depths, the anomalies can be due to downwelling fluids coming from the surrounding topographic highs, thus defining a recharge zone in the graben, or to fluid circulation in permeable faulted zones.

A focus on the Soultz-sous-Forêts anomaly

In order to investigate more thoroughly the possible fluid circulation process within the thick parts of the permeable fractured Buntsandstein formation, additional cross-sections through the Soultz-sous-Forêts anomaly (anomalie (A) in Figure 6) have been extracted from the 3D thermal block. Figure 8 In accordance with previous observations, sections A-B and B-C show that the hot zone is confined within the thickest parts of Buntsandstein sandstones. The two white arrows point to sharp faults which limit the extent of the permeable formation. It must be emphasized that the spreading trend of the anomaly along the inclined roof of the Buntsandstein formation towards higher levels (shallower depths), which was described with the west-east cross-section, is also visible in this A-B-C cross-section. The role of the slope (inclined basement-sediments interface) may be of major importance since in the D-F cross-section, where slope angle is very low, spreading of the anomaly is not visible and a quasi-symmetric blob-like pattern is observed.

One may also note that Soultz-sous-Forêts is not located above the hottest zone of the anomaly (see cross-section A-B), but as shown in the 3D view (bottom part of Figure 8), the 140 ± 9 °C isotherm is reached at shallowest depths around Soultz-sous-Forêts. If temperatures above 160 °C were to be explored, then one should preferentially move eastwards to Rittershoffen and drill down to 1800 m (see cross-section A-B-C).

In order to test these possible different fluid circulation patterns within the graben, numerical models of hydrothermal convection within permeable zones, at different scales, are described in the next section.

Fluid circulation models

Objectives of the simulations

A C C E P T E D M A N U S C R I P T

Thanks to the updated database and to the resulting 3D thermal block, we have performed numerical simulations where the primary objective is to reproduce the measured temperature field. First, the temperature profile measured in the geothermal GPK2 borehole will be used to discriminate between different models, even if it is only one-dimensional information. As it will be shown, several models can reproduce this temperature profile. Consequently, construction of numerical simulations will also be driven by the reproduction of the abovedescribed two-dimensional features of the Soultz-sous-Forêts anomaly (Figure 8). Apart from simulations, numerical experiments have also been performed to understand specific effects on the structure of hydrothermal convection. In particular, the role of the inclined basementsediments interface will be tested, as well as the role of the horst structure beneath Soultzsous-Forêts.

In order to simulate a realistic phenomenon, one must account for a representative permeability distribution, and for realistic fluid and rock properties. Numerical models presented below are two-dimensional, which is justified when the geometrical structure of the Rhine graben is considered. However, because the basement-sediments interface from north to south of the graben show strong topographic variations (see Figure 1a and black lines in cross-sections of Figure 7), different possibilities will be discussed in the frame of a threedimensional circulation pattern.

Governing equations

Thermal convection within permeable and porous media can be numerically studied by coupling heat transfer and fluid flow equations with appropriate fluid and rock physical properties. Physically, fluid flow in fractured or porous media can be described by the Darcy law. Here, the sedimentary rocks are considered as more or less permeable layers, within which fluids can circulate if permeability is sufficiently large, whereas the granitic basement is a fractured permeable medium, where fluids can also circulate. Heat transfer equation writes:

  T u C T t T C L L p             (4) A C C E P T E D M A N U S C R I P T
where T is temperature, t is time and u the fluid filtration (or Darcy) velocity. Properties  C p and  correspond to density, specific heat and thermal conductivity of the saturated porous medium, respectively;  L and C L are fluid density and fluid specific heat. Properties of the saturated porous medium obey to a simple mixing rule between the fluid (L) and the solid (s):

    s p L p p C C C      ) 1 (    (5) s L      ) 1 ( ) (    (6)
where  is porosity. Fluid is supposed incompressible and mass conservation equation writes:

0    u (7)
For porous and permeable formations or for a fractured medium, fluid motion is driven by pressure gradient and buoyancy, and fluid filtration velocity obeys Darcy's law:

  g p k u L       (8)
where k is the permeability of the medium,  is the fluid dynamic viscosity, p is pressure and g is acceleration of gravity. Hence, coupling between heat transfer and fluid motion occurs through the velocity field, which appears in equations ( 4), ( 7) and (8). In addition, temperature (the other coupling variable) appears in equations ( 4) and (8) since fluid properties ( L and ) are temperature-dependent, as it is detailed later.

Geometries and boundary conditions

Several geometrical configurations have been tested and two series of experiments have been performed. The first series is dedicated to simple rectangular box simulations (Figure 9a) where (i) the effect of a depth-dependent permeability and (ii) the effect of an inclined basement-sediments interface, on convective patterns and distribution of isotherms are investigated. The second series of experiments includes a horst structure within a basinscale model (Figure 9b), where depth-dependent permeability is accounted for. Figure 9 shows dimensions of the models and the chosen boundary conditions. Thermal boundary conditions correspond to a fixed temperature (10 °C) at the surface and a fixed heat flow at the base of the model (90 mW.m -2 at 6 km depth). For one model, a fixed temperature condition of 200°C at 5 km depth was imposed. Lateral boundaries are thermally insulating. A fixed atmospheric pressure is imposed at the surface while a zero flux condition is imposed on all other boundaries (

0   n u
, where n is a unit vector normal to the boundary).

Petrophysical and fluid properties

Fluid density and viscosity depend on temperature. The following fitting polynomial trend for density has been chosen (www.engineeringtoolbox.com):

  2 * 0025 . 0 * 1905 . 0 4 . 1002 T T T L     (9)
where T is temperature in °C and  L the fluid density in kg.m -3 . For the fluid dynamic viscosity, the law used by [START_REF] Smith | On the thermal effects of groundwater flow -1. Regional flow systems[END_REF] or by [START_REF] Rabinowicz | Two-and three-dimensional modeling of hydrothermal convection in the sedimented Middle Valley segment, Juan de Fuca Ridge[END_REF] As emphasized by many authors, rock permeability may vary by several orders of magnitudes in the first 5 kilometres [START_REF] Ingebritsen | Geological implications of a permeability-depth curve for the continental crust[END_REF][START_REF] Saar | Depth dependence of permeability in the Oregon Cascades inferred from hydrogeologic, thermal, seismic, and magmatic modelling constraints[END_REF][START_REF] Stober | Hydraulic properties of the crystalline basement[END_REF] and should also be considered as a time-dependent property [START_REF] Ingebritsen | Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism[END_REF]. In this study, we only considered steady-state thermal signatures due to a heterogeneous crustal structure and composition, where hydrothermal circulation occurs. As suggested by permeability measurements and estimates (see section 2.3), a constant low permeability (10 -18 to 10 -17 m²) will be assigned to the sedimentary cover of the upper kilometres,. Permeability of the two main faults defining the horst structure (the Soultz and the Kutzenhausen faults) will be varied from 10 -16 to 3 10 -14 m². The underlying

A C C E P T E D M A N U S C R I P T

Buntsandstein sandstone formation will be considered as a permeable layer (from 5 10 -15 to 10 -14 m²) but with no depth-dependence (its thickness being lower than 500-600 m). For the granitic basement down to 3700 m, the following depth-dependence of permeability is chosen:

         z k z k G G exp ) ( 0 (11)
as in [START_REF] Garibaldi | Thermal anomalies and geological structures in the Provence basin : implications for hydrothermal circulations at depth[END_REF]. The depth-scale  may vary according to the intensity of fracture density and to the lithology. Because the Soultz-sous-Forêts granite shows a pervasive fracturing content down to ~ 3700 m [START_REF] Pribnow | Thermal tracking of Upper Crustal Fluid Flow in the Rhine Graben[END_REF][START_REF] Genter | Contribution of the exploration of deep crystalline fractured reservoir of Soultz to the knowledge of Enhanced Geothermal Systems (EGS)[END_REF], the depth-scale should be chosen so that the depth-dependent decrease is not too strong.

Determination of k G 0 and  are detailed below. Below 3700 m depth, the granitic basement is much less permeable and a constant value of 10 -18 m² is assumed. All other petrophysical properties (thermal conductivity, heat production rate, heat capacity, porosity) are taken from literature and previous hypotheses [START_REF] Flores-Marquez | Transferts de chaleur et de masse en milieu sédimentaire et fracturé. Modélisation numérique de la convection naturelle autour du site de Soultz (Graben du Rhin)[END_REF][START_REF] Carlier | Convective heat transfer at the Soultz-sous-Forêts geothermal site: implications for oil potential[END_REF][START_REF] Clauser | Decoupled thermal and mantle helium anomalies : implications for the transport regime in continental rift zones[END_REF], and are indicated in Table 2.

Numerical method

Coupling of equations 4, 7 and 8 has been performed with the Comsol Multiphysics™ software (finite element method), with which temperature (or space) dependent physical properties can be accounted for (equations 9, 10 and 11). Benchmark experiments (including the reproduction of numerical simulations by [START_REF] Gerdes | Convective fluid flow through heterogeneous country rocks during contact metamorphism[END_REF][START_REF] Rabinowicz | Two-and three-dimensional modeling of hydrothermal convection in the sedimented Middle Valley segment, Juan de Fuca Ridge[END_REF] have already been performed and published [START_REF] Eldursi | Numerical investigation of transient hydrothermal processes around intrusions: heat-transfer and fluid-circulation controlled mineralization patterns[END_REF][START_REF] Garibaldi | Thermal anomalies and geological structures in the Provence basin : implications for hydrothermal circulations at depth[END_REF].

For the largest models (at the graben-scale with a horst structure), mesh consisted in 29,600 triangles, whose size was not larger than 20 m in highly permeable zones (faults and top of granitic basement in particular). A total of 120,000 degree of freedom was reached.

During computation, the initial thermal field (a steady-state conductive regime) evolves towards a steady-state convective regime in a few Myrs; the initial pressure field increases

A C C E P T E D M A N U S C R I P T

linearly with depth. The final steady-state convective regime is then considered as a possible solution.

After each experiment, the first validation test consisted in checking the temperature profile near the Soultz-sous-Forêts area. A number of petrophysical parameters could have been varied to get closer to the measured temperature profile. However, as it is well-known, permeability and its depth variation represent the main unknown properties. Consequently, physical parameters of Table 2 were imposed for each model, while k G 0 and the associated depth-scale  have been varied (i) in order to keep a permeability value of ~ 10 -14 m² at 1.4 km depth in the Soultz-sous-Forêts area, (ii) with the objective of reproducing the typical thermal signatures illustrated in the 3D thermal block, and (iii) in order to retrieve the measured temperature profile at Soultz-sous-Forêts. Achievement of these three criteria allows selecting the result as a possible solution.

Results

Simple models reproducing the Soultz-sous-Forêts temperature profile

A number of simplified models can easily reproduce the typical temperature profile of Soultzsous-Forêts. Some examples have already been published (e.g. [START_REF] Clauser | Analysis of convective and conductive heat transfer in a sedimentary basin, demonstrated for the Rhinegraben[END_REF][START_REF] Carlier | Convective heat transfer at the Soultz-sous-Forêts geothermal site: implications for oil potential[END_REF][START_REF] Kohl | Steps towards a comprehensive thermo-hydraulic analysis of the HDR test site Soultz-sous-Forêts[END_REF][START_REF] Clauser | Decoupled thermal and mantle helium anomalies : implications for the transport regime in continental rift zones[END_REF][START_REF] Guillou-Frottier | La convection hydrothermale et les ressources associées[END_REF], the last one being computed with several horizontal homogeneous layers of varying permeability.

However, in these studies, the obtained convective patterns do not show a single and large upwelling as it seems to be the case according to the 3D thermal block (Figure 8). It must be noted that even if uncertainties of Figure 5 are accounted for (± 9°C at 2000 m depth, ± 12°C at 3000 m depth), this single and large upwelling cannot disappear from thermal crosssections. Since convective patterns are largely governed by permeability values, but also by the geometry of the permeable layer, simple tests with a depth-dependent permeability and with an inclined basement-sediments interface are presented below.

The role of a permeability decrease

A C C E P T E D M A N U S C R I P T

Figure 10a illustrates several cases where permeability structure is varied from a constant value of 10 -14 m² (top experiment) to a strong decrease from 10 -14 m² to 3 10 -16 m² (bottom case). Values 0 G k and  have been adjusted to get ~ 10 -14 m² at the basement-sediment interface, and to reproduce the linear decrease shown on the left of each convective pattern.

The upper homogeneous unit is here considered as impervious (10 -19 m²). For a constant permeability, the steady-state convective regime shows 4 fully-developed upwelling plumes and 2 half plumes at box edges. A slight decrease of permeability with depth enlarges convective wavelength since one half-plume has disappeared. For more important decrease of permeability with depth, one gets 4 upwellings (third case) and 3.5 upwellings for the last case.

One may note that for each case, temperature at mid-depth of the convective layer oscillates around 140 °C, which is understandable since thermal boundary conditions are identical for each experiment. However, convective patterns are strongly affected since velocity field is mainly controlled by permeability. For example, the averaged Darcy velocity at mid-depth equals 7 cm.yr -1 in the top case, while it is decreased by a factor 3 in the bottom experiment. Actually, when permeability decreases with depth, location of maximum values of velocity is shifted towards the interface. Thus, the cold downwelling fluid can circulate more easily while the hot fluid at depth is surrounded by less permeable rocks. This vertical asymmetry in the flow conditions of the three bottom experiments results, from top to bottom, in a decreasing lateral width of upwellings (and an increasing lateral width of downwellings), involving more focused upwellings and more diffuse downwellings. This phenomenon has important implications when a 3D system is considered, since upwellings would probably be less numerous and with a plume-like geometry, while geometry of downwellings should occupy much larger volumes.

The role of an inclined topography at the top of the basement

A C C E P T E D M A N U S C R I P T

Figure 10b illustrates several cases where the slope of the basement-sediment interface has been increased from 0 to 15 %, knowing that the W-E cross-section in Figure 8 indicates a slope of 11%. Here, 0 G k equals 10 -13 m² and  equals 800 m, so that a permeability of 4 10 -15

m² is obtained at the basement-sediments interface, at 2500 m depth. In all cases, interface crosses the 2500 m depth at the centre of the model. With no slope (top experiment in Figure 10b), two full upwelling plumes and two half upwelling plumes (at box edges) are observed.

With a slope of 5%, the two half plumes have disappeared and one fully-developed plume is obtained, the other one being less developed than in the case of a flat interface. In the case of a larger slope (two bottom experiments), only one fully-developed upwelling remains while the less developed plume disappears. It must be noted that permeability values at the basement-sediments interface differs from one case to the other, involving higher permeability on the left side at shallow depths for a 15% slope than for a 10% slope. An interesting result is illustrated by the asymmetry of the obtained convective pattern due to the inclined interface. In particular, isotherms around the isolated fully-developed plume are distorted towards shallow levels (left side in Figure 10b), where permeability is higher (see thick white isotherms in the last case).

In other words, an inclined interface tends to favour the development of a single and asymmetric upwelling plume. From these 2D experiments, one could infer that a threedimensional topography (thus including inclined interface in several directions) would probably reinforce this effect. In other words, it may not be surprising to get only one strong anomalously hot zone in the FPURG.

Another observation deals with the predominance of a cold downwelling beneath the thinnest sedimentary cover. This can be simply explained by the fact that the hot fluid flowing upwards along the interface is more and more efficiently cooled from above, since it encounters colder and colder rocks. It follows that density contrast at these shallow depths is larger than density contrasts of fluids moving downwards along the interface. Consequently,

A C C E P T E D M A N U S C R I P T

cold downwelling is promoted on the left side of the box (thin sedimentary cover), while hot upwelling is favoured beneath thicker parts of the sedimentary cover.

Additional experiments similar to those of Figure 10b have been carried out with a slope of 11%, and with different 0 G k values, ranging from 5 10 -14 to 4 10 -13 m². Same results

were observed: a single narrow hot upwelling was locked around the centre of the box and a large cold downwelling was promoted on the left side.

5.6.4 Comparison with the west-east cross section extracted from the 3D thermal block

Figure 11 shows a comparison between modelled temperatures obtained with a slope of 11%

and temperature data on a W-E cross-section near Soultz-sous-Forêts. The four models are computed with distinct permeability distributions (shown in Figure 11b). As emphasized by thick white lines, distortion of the 140 °C isotherm observed in Figure 11c can be reproduced with the first experiment (top of Figure 11a), while the 160 °C isotherm is partly reproduced with the most permeable granite (bottom case). Clearly, permeability distribution in the FPURG should differ from that of the model, but it appears that accounting for both an inclined interface and for a depth-dependent permeability can explain the shapes of isotherms inferred from W-E thermal cross-sections.

5.6.5 Models at graben-scale with a horst structure

The models illustrated in Figure 12 and 13 account for the peculiar topography of the basement-sediment interface through an entire east-west cross section, as well as for the horst structure. In Figure 12, permeability of the two main faults which define the horst is low (10 -17 m²), in order to separate the distinct sources (the Buntsandstein fractured formation and the major faults) for temperature fluctuations. In Figure 13, this value is varied from 10 -15 to 3 10 -14 m². Several permeability distributions have been tested, and to avoid numerical divergence at shallow depths, a low permeability (10 -17 m²) was imposed at external edges of the box, where basement-sediment interface is shallower than 1600 m. When the

A C C E P T E D M A N U S C R I P T

temperature profile measured in the GPK2 well at Soultz-sous-Forêts was not reproduced, results from the simulations were not retained.

Two examples are shown in Figure 12a and12b, where maximum surface heat flow reaches 120 mW.m -² above the horst structure. In Figures 12a,a fixed temperature of 200 °C is imposed at 5 km depth while a fixed heat flow of 90 mW.m -² is imposed at the base of the model in Figure 12b. Actually, for such basin-scale models, there is no physical reason to impose a fixed temperature at a depth of a few kilometres. We however show the differences between these two conditions. As it can be seen, upwellings are more easily developed when a fixed temperature is imposed. In some cases, a negative temperature gradient can be obtained at 3 km depth beneath the horst structure because thermal constraints are too strong. On the opposite, a fixed basal heat flow (which is a more reasonable thermal boundary condition) allows temperature variations at depth (Figure 12b), preventing a regular juxtaposition of upwellings. It follows that less upwellings are permitted when a fixed heat flow is imposed.

Figure 12c shows the permeability distribution used in Figure 12b, together with the velocity field (white arrows) and several isotherms (white lines). Figure 12d illustrates a vertical profile of permeability at the GPK2 borehole. Figure 12e shows the temperature profiles obtained with experiments shown in Figure 12a (blue) and 12b (green). One note that the used thermal boundary condition does not play a significant role in the shape of the modelled geotherm, while it is of major importance on convective pattern, as explained above. The measured temperature profile (red) and the modelled ones are comparable but not exactly superimposed. In the upper 1000 m, the model apparently underestimates temperature values. This can be due to the simplification of a single lithology in the sedimentary pile overlying the Buntsandstein formation, or to the too strong hypothesis of impervious major faults. The role of fault zone permeability on temperature profile is studied in Figure 13.

Maximum Darcy velocity reaches 5 cm/yr in the Buntsandstein formation while fluid circulates within the fractured granite at an averaged Darcy velocity of ~10-15 cm.yr -1 . These

A C C E P T E D M A N U S C R I P T

values are in agreement with those found by previous studies (Le [START_REF] Carlier | Convective heat transfer at the Soultz-sous-Forêts geothermal site: implications for oil potential[END_REF][START_REF] Clauser | Analysis of convective and conductive heat transfer in a sedimentary basin, demonstrated for the Rhinegraben[END_REF][START_REF] Kohl | Steps towards a comprehensive thermo-hydraulic analysis of the HDR test site Soultz-sous-Forêts[END_REF], and fit the hydraulic tests in geothermal boreholes [START_REF] Jung | Hydraulic fracturing and hydraulic testing in the granitic section of borehole GPK1, Soultz-sous-Forêts[END_REF]. One can also see on Figure 12c that temperatures of 150 °C are located near 2500 m depth beneath the thickest sediments, while they are reached much deeper beneath Soultz-sous-Forêts, at 3500 m depth.

In Figure 13a, permeability of major faults (k F ) is increased from top to bottom. Figure 13b shows temperature profiles at GPK2 location. In the first case (10 -15 m²), convective pattern in the granitic basement is similar to previous cases (Figure 12) where faults were impervious. For a fault permeability of 10 -14 m², the suction effect of the faults becomes significant and attracts the hot upwelling beneath the horst structure, thus reheating slightly the overlying layers (see the green temperature profile in Figure 13b). For a larger value (3.

10 -14 m²), a smaller convective cell establishes between the two faults, inducing a downwelling beneath the horst structure, and thus a too cold temperature profile (blue curve in Figure 13b).

Discussion

The obtained results on 2D numerical models of hydrothermal convection within the permeable fractured granite and its overlying permeable sandstones, show that the 1D measured temperature profile at Soultz-sous-Forêts (Figure 12e and 13b) can be correctly reproduced, and that similar 2D thermal features are illustrated both in Figure 7 and 8 (crosssections from the 3D thermal block, constrained by corrected BHT data) and in Figure 12 and 13 (2D numerical models). Temperature profile and 2D thermal cross-sections have thus been reproduced by a 2D numerical model accounting for realistic physical properties and geological information. Consequently, at this stage, a 3D numerical model does not appear necessary since temperature data have been reproduced. It is interesting to note that a fault permeability value of 10 -14 m² yields a better reproduction of the temperature profile.

As explained in section 3, the interpolation procedure may not be the most accurate one, since a universal kriging method would be based on a local fitting of the drift. However,

A C C E P T E D M A N U S C R I P T

choosing an average temperature increasing trend provide us with thermal anomalies defined within the entire FPURG. This approach appears more useful than another one where anomalies would not refer to the same vertical trend.

During data treatment, a number of temperature data have been simply removed (inconsistent data, too many similar data at the same locations, etc). The interpolation procedure could have used these data by assigning them a less important weight, but we considered that a database with 420 corrected temperature data was largely sufficient to get rid of possible biased data.

The obtained 3D thermal block must be provided with the used variogram, that is to say that maps and cross-sections have to be presented with the error on interpolation (Figure 5). It must be pointed out that in well-covered areas, standard deviation does not exceed 5°C

(Figure 5). In the model variogram, the largest horizontal range reaches 20 km, a distance which can be associated with the size of thermal anomalies (see Figures 5 and8). It must be noted that similar values were found in the Paris basin and in the south-east basin of France [START_REF] Bonté | Subsurface temperature maps in French sedimentary basins : new data compilation and interpolation[END_REF][START_REF] Garibaldi | Thermal anomalies and geological structures in the Provence basin : implications for hydrothermal circulations at depth[END_REF]. Vertically, the size of thermal anomalies would be around 1200-1600m, indicating that a vertical superimposition of cold and hot anomalies is not excluded (see examples around faulted zones of the south-east basin, [START_REF] Garibaldi | Thermal anomalies and geological structures in the Provence basin : implications for hydrothermal circulations at depth[END_REF]. However, in the FPURG, anomalies are more developed where the permeable fractured sandstone formation is sufficiently thick, as it is particularly illustrated in Figure 7b (from distance 0 to ~100 km), and not necessarily correlated with faulted zones, as previously suggested by [START_REF] Kohl | Steps towards a comprehensive thermo-hydraulic analysis of the HDR test site Soultz-sous-Forêts[END_REF] or by [START_REF] Clauser | Decoupled thermal and mantle helium anomalies : implications for the transport regime in continental rift zones[END_REF]. The fact that main thermal features can be reproduced with low fault permeability (e.g. Figure 12b) indicates that convective patterns are predominantly controlled by the thickness of the permeable zone, including both the Buntsandstein formation and the permeable top granitic basement.

Note however that details in the temperature profile are better fitted when the two major faults are considered permeable (Figure 13).

One major key point of our study lies in the fact that basin geometry, depth-dependent permeability, and realistic basal thermal boundary condition explain why the hot thermal

A C C E P T E D M A N U S C R I P T
anomaly is located near Soultz-sous-Forêts, where a horst structure is present. First, the inclined basement-sediments interface, the depth-dependent permeability and the fixed heat flow condition help to decrease the number of possible upwellings. Second, when fault permeability is low, the upwelling is attracted by regions where the sedimentary cover is locally the thickest. Since the horst structure involves some topography of the interface, an upwelling is locked beneath the thickest sedimentary unit, about 1 km north-west of GPK2. It follows that Soultz-sous-Forêts is not necessarily the hottest zone in the area, as suggested by 2D cross-sections in Figure 8 and illustrated by models of Figure 12. This point is in agreement with recent 2D magnetotellurics data where a conductive structure (the Kutzenhausen fault) was also visible west of the horst structure [START_REF] Geiermann | 2-D Magnetotellurics at the geothermal site at Soultz-sous-Forêts: Resistivity distribution to about 3000 m depth[END_REF]. However, when fault permeability is accounted for, the upwelling tends to be shifted beneath the Soultz-sous-Forêts fault. Since permeability is probably a time-dependent parameter controlled by mineral dissolution and precipitation processes, one can expect the hot upwelling to be located between the thickest sedimentary cover and the Soultz fault.

Another interesting feature appears on the models of Figure 12b and12c. A second upwelling develops about 8 km east of the Soultz-sous-Forêts area (Rittershoffen area).

Although it is less-developed in terms of lateral extent, it shows higher temperatures at relatively shallow depths: a temperature of 160 °C is reached at 2300 m depth, whereas such temperatures are found below 3700 m in the Soultz-sous-Forêts area. This second upwelling is less visible when the 140 °C isotherm is examined (see the 3D view in Figure 8). It must be emphasized that this second thermal anomaly, which was already suspected several decades ago [START_REF] Laurent | Etudes préliminaires des possibilités géothermiques de la Basse Alsace, entre Schirrhein et Scheibenhard (Bas-Rhin)[END_REF] was not observed in previous numerical models. More precisely, the three upwellings that were found by [START_REF] Carlier | Convective heat transfer at the Soultz-sous-Forêts geothermal site: implications for oil potential[END_REF] within a 10 km lateral distance had decreasing amplitudes when going eastward. In their models, permeability was not depth-dependent and the basal thermal boundary condition was a varying heat flow condition. Today, the area of Rittershoffen corresponds to a future new geothermal target (http://www.lalsace.fr/actualite/2011/06/07/roquette-alimentee-en-2014par-la-geothermie-profonde) and our models, which account for realistic hypotheses and

A C C E P T E D M A N U S C R I P T
conditions, tend to better describe how hydrothermal convection settles in the FPURG.

However, as for any modelling study, our suggested interpretations and models might be non-unique, especially if 3D numerical modelling was undertaken. It must be emphasized that numerical modelling of hydrothermal convection within a basin-scale model with representative features of the FPURG resulted in the development of an upwelling around the Soultz-sous-Forêts area. Besides reproducing results obtained in Figure 8, numerical models also help to understand how fluid circulation establishes, and why a second hotter upwelling can be generated at a few kilometres eastwards.

In the German part of the URG, the Landau geothermal anomaly is also located on the western side of a similarly inclined asymmetrical basin structure (see thickness variations in Figure 1a). The horizontal distance between Soultz-sous-Forêts and Landau (~ 34 km) may actually correspond to the critical horizontal wavelength of 3D hydrothermal convection patterns in such peculiar geometries. Three-dimensional models including an inclined basement-sediment interface should bring new insights on spatial distribution of hydrothermal plumes within permeable zones.

Several hypotheses have not been discussed in our study. The large value of the basal heat flow, which partly explains why hydrothermal circulation can occur, has been chosen so that surface heat flow value of 120 mW.m -² can be recovered above the Soultzsous-Forêts area. The origin of this high basal heat flow remains unknown, but as discussed earlier, one part of this high basal heat flow could correspond to the mantle uplift signature, whose imprint should have diffused towards the surface 30 Ma after the beginning of rifting.

The relatively high heat production of the top granitic basement may also indicate that middle crust could be particularly radiogenic, thus providing an excess heat from natural radioactivity. According to magneto-telluric experiments by [START_REF] Spichak | Deep temperature extrapolation in the Soultz-sous-Forêts geothermal area using magnetotelluric data[END_REF], a deeper convective cell would be present beneath 5 km depth, but in that case, basement should be more permeable than between 3500 and 5000 m. Since local stress field seems to change at

A C C E P T E D M A N U S C R I P T

seems difficult to get high permeability at deeper levels since this parameter generally decreases with depth (Manning and Ingebristen, 1999).

Our models have been performed in the frame of a steady-state convective regime.

Likewise, permeability is assumed to be time-independent, whereas it is well known that fractures and faults can seal and re-open with time, through multiple tectonic constraints (e.g.

Ingebristen and [START_REF] Ingebritsen | Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism[END_REF]. Present-day temperature profile is in agreement with fluid inclusion studies on secondary quartz veins [START_REF] Dubois | Détermination des conditions de précipitation des carbonates dans une zone d'altération du granite de Soultz (soubassement du fossé Rhénan, France) : l'enregistrement des inclusions fluides[END_REF], which correspond to present-day sealing of permeable zones [START_REF] Genter | Over 10 years of geological investigations within the European Soultz HDR project, France[END_REF]. In other words, the hydrothermal systemin the depth range of 1000-2200 m investigated by fluid inclusion studieswould be in a sealing phase. It is clear that further modelling studies should now account for time-dependent physical properties.

Conclusion

This study was based on a new compilation (and correction) of temperature data, allowing to build a 3D thermal block, where error on interpolation is computed. 3D imaging of thermal anomalies illustrate strongly uplifted isotherms, even when uncertainty is accounted for; these anomalously hot zones can only be explained by fluid circulation processes, which are mainly favoured in areas where the permeable Buntsandstein formation is particularly thick.

Accounting for realistic petrophysical and fluid properties, a series of 2D numerical models helped to understand how convective cells established beneath an inclined basementsediments interface. Physics of hydrothermal convection within a medium of varying permeability is quite complex, but it also obeys to simple physical rules as soon as buoyancy and temperature contrasts are not evenly distributed (depth-dependent). Besides

  produced a sharp definition of the basement-sediments interface[START_REF] Place | Decoupling of deformation in the Upper Rhine graben sediments. Seismic reflection and diffraction on 3-component vertical seismic profiling (Soultz-sous-Forêts area)[END_REF]. Several boreholes penetrated Cenozoic, Mesozoic and Permian sediments which have been deposited on the Paleozoic basement. The Tertiary filling sediments are marine and lacustrine limestones, marls and evaporites, including the petroleum layers of Pechelbronn, unconformably overlying the Jurassic limestones, the German Triassic layers (Keuper, Muschelkalk, Buntsandstein) and the Permian sandstones. The Buntsandstein formation represents about 400 m thick of fluvial sandstones. The unconformity between Cenozoic and Mesozoic is outlined by the absence of Cretaceous and late Jurassic sediments. The Paleozoic basement includes porphyritic monzo-granite and two-mica granite.

  , a compilation made by the local geological survey, and the recent compilation of digitalized logging headers of all petroleum boreholes. Among the 324 boreholes present in the FPURG, 251 possess one or more temperature data. After removal of inconsistent or duplicate data, the final database contained 509 temperature data, among which 34 correspond to equilibrium temperatures, 152 are close to equilibrium (such as temperature data from drill-stem tests (DST)), and 323 BHT data, that have to be corrected for transient disturbances, such as circulation of cold mud a few hours before measurement. Statistics on absolute values and vertical repartition of the corrected temperature data is shown in Figure 3. The deep A C C E P T E D M A N U S C R I P T

  A C C E P T E D M A N U S C R I P TTo account for spatial correlation of the residuals of temperature data after removing the vertical drift, an experimental variogram of the residuals was computed. The variogram is defined by a correlation, represented by the variogram value (h) versus the distance h between two values. Calculation of the variogram for N(h) pairs of measured (here residuals)

  A C C E P T E D M A N U S C R I P T models of variograms, model C3 (comprising a nugget effect and two Gaussian structures) was the one which had the best fit to the experimental variogram. This model (Figure 4) has a nugget effect of 25 °C². The sills of the two Gaussian structures are respectively 60 °C² and 200 °C². The horizontal ranges are 1,800 and 20,000 m, and the vertical ranges are 1200 m

  anomaly", which actually extends over a horizontal distance greater than 20 km below 1500 m. Maximum absolute temperature of this anomaly reaches 82 °C at 500 m depth, 109 °C at 800 m depth (120 °C according to[START_REF] Pribnow | Thermal tracking of Upper Crustal Fluid Flow in the Rhine Graben[END_REF], 127 °C at 1000 m depth and 152 °C at 1500 m depth. Being centred below Soultz-sous-Forêts city down to 1000 m depth, the anomaly extends eastwards towards Oberroedern below 1500 m and then northwards (towards Scheibenhard) from 2000 m to 3000 m. This asymmetric 3D character of the anomaly is presented and discussed later.A C C E P T E D M A N U S C R I P TThe positive anomaly south of Strasbourg (C) is weak (a few °C) at 500-1000 m depth but slightly increases in amplitude (up to ~10°C) down to 2500 m depth, without significant changes in geometry and location. The Altkirch anomaly (H) is rather strong (~15°C) at 500 m depth but disappears below 1500 m. It must be noted that this anomaly is located at the edge of the model and because it is constrained by only 2 data points, the existence of the anomaly might be questioned.The amplitude of the negative anomaly north of Strasbourg (B) is almost as strong as the one at Soultz-sous-Forêts (A) (see north-south cross section in Figure7). Lateral temperature variations of 70-80 °C are reached at depths of 1200-1500 m. Another common feature between (A) and (B) anomalies is their spatial evolution with depth. While the hot anomaly extends eastwards and then northwards (thus defining an anticlockwise rotation), here, the cold anomaly shows a clockwise rotation with increasing depth (elongated in a NE-SW direction at 500 m depth, it evolves regularly towards a NW-SE direction at 2500 m depth).The negative anomaly of Colmar-Mulhouse-Sélestat (E) stays quite large above 1500 m depth (~-10 to -20°C), but decreases in size and amplitude below. The other minor anomalies (D, F and G) are linked to a small number of data and can thus be questioned.

  presents a temperature map at 2000 m depth (top left), four crosssections (two of them being already presented in Figure7) and a 3D view of the 140 ± 9 °C isotherm. Cross-sections A-B, B-C and D-F show new pictures of the hot zone, together with topography of the sediment-basement interface. As it can be seen in maps of standard deviations (Figure5), uncertainties around the Soultz-sous-Forêts anomaly never exceed 12°C. Standard deviation is around 8-9°C down to 2000 m depth and slightly increases to 12°C at 3000 m depth.

  where T is in °C and dynamic viscosity is in Pa.s. Both laws reproduce experimental data measured on pure water between 4 and 250 °C.

Figure 1

 1 Figure captions

Figure 2 :

 2 Figure 2: Previous subsurface temperature maps with no indication on interpolation

Figure 3 :

 3 Figure 3: Statistics on uncorrected temperature data (a,b) and depth distribution of the entire

Figure 4 :

 4 Figure 4: Plane and vertical variogram for the temperature residuals (see text).

Figure 5 :Figure 6 :

 56 Figure 5: Extractions from the 3D thermal block: temperature maps at different depths with

Figure 7 :

 7 Figure 7: a) west-east cross-sections (see locations in Figure 6). Green vertical lines indicate

Figure 8 :

 8 Figure 8: Maps, cross-sections and 3D view of thermal field in the Soultz-sous-Forêts area.Note the asymmetrical shape of the hot zone when the basement interface is inclined,

Figure 9 :

 9 Figure 9: Geometrical and boundary conditions studied by numerical modelling of

Figure 10 :

 10 Figure 10: a) Effect of a depth-decrease of permeability on convective patterns (isotherms

Figure 12 :

 12 Figure 12: Models at basin-scale with parameters listed in Table2, with low permeability for

Figure 13 :

 13 Figure 13: Models at basin-scale with a varying permeability (k F ) for the two main faults. Left: Thermal and velocity fields. The associated temperature profiles at GPK2 location are shown

Figure 2 A

 2 Figure 2

Figure 3 A

 3 Figure 3



  Figure 10

Acknowledgements

One part of this work corresponds to an improved version of Clément Carré Master's thesis.

This study was funded by ADEME (French Environment and Energy Management Agency) and BRGM (French Geological Survey) within the framework of the CLASTIQ-2 research project. We would like to thank F. Lucazeau for reviewing an earlier draft of this manuscript, and Simon Lopez for his advice regarding geostatistics and data treatment. We thank P.

Jousset and an anonymous reviewer for their valuable comments.