Get 20M+ Full-Text Papers For Less Than $1.50/day. Subscribe now for You or Your Team.

Learn More →

Locating Geothermal Resources: Insights from 3D Stress and Flow Models at the Upper Rhine Graben Scale

Locating Geothermal Resources: Insights from 3D Stress and Flow Models at the Upper Rhine Graben... Hindawi Geofluids Volume 2019, Article ID 8494539, 24 pages https://doi.org/10.1155/2019/8494539 Research Article Locating Geothermal Resources: Insights from 3D Stress and Flow Models at the Upper Rhine Graben Scale Antoine Armandine Les Landes , Théophile Guillon, Mariane Peter-Borie, Arnold Blaisonneau, Xavier Rachez, and Sylvie Gentier BRGM, Orléans, France Correspondence should be addressed to Antoine Armandine Les Landes; a.armandineleslandes@brgm.fr Received 20 December 2018; Revised 1 March 2019; Accepted 12 March 2019; Published 12 May 2019 Guest Editor: Victor Vilarrasa Copyright © 2019 Antoine Armandine Les Landes et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. To be exploited, geothermal resources require heat, fluid, and permeability. These favourable geothermal conditions are strongly linked to the specific geodynamic context and the main physical transport processes, notably stresses and fluid circulations, which impact heat-driving processes. The physical conditions favouring the setup of geothermal resources can be searched for in predictive models, thus giving estimates on the so-called “favourable areas.” Numerical models could allow an integrated evaluation of the physical processes with adapted time and space scales and considering 3D effects. Supported by geological, geophysical, and geochemical exploration methods, they constitute a useful tool to shed light on the dynamic context of the geothermal resource setup and may provide answers to the challenging task of geothermal exploration. The Upper Rhine Graben (URG) is a data-rich geothermal system where deep fluid circulations occurring in the regional fault network are the probable origin of local thermal anomalies. Here, we present a current overview of our team’sefforts to integrate the impacts of the key physics as well as key factors controlling the geothermal anomalies in a fault-controlled geological setting in 3D physically consistent models at the regional scale. The study relies on the building of the first 3D numerical flow (using the discrete-continuum method) and mechanical models (using the distinct element method) at the URG scale. First, the key role of the regional fault network is taken into account using a discrete numerical approach. The geometry building is focused on the conceptualization of the 3D fault zone network based on structural interpretation and generic geological concepts and is consistent with the geological knowledge. This DFN (discrete fracture network) model is declined in two separate models (3D flow and stress) at the URG scale. Then, based on the main characteristics of the geothermal anomalies and the link with the physics considered, criteria are identified that enable the elaboration of indicators to use the results of the simulation and identify geothermally favourable areas. Then, considering the strong link between the stress, fluid flow, and geothermal resources, a cross-analysis of the results is realized to delineate favourable areas for geothermal resources. The results are compared with the existing thermal data at the URG scale and compared with knowledge gained through numerous studies. The good agreement between the delineated favourable areas and the locations of local thermal anomalies (especially the main one close to Soultz-sous-Forêts) demonstrates the key role of the regional fault network as well as stress and fluid flow on the setup of geothermal resources. Moreover, the very encouraging results underline the potential of the first 3D flow and 3D stress models at the URG scale to locate geothermal resources and offer new research opportunities. 1. Introduction it is environmentally friendly, has low emissions (greenhouse gases), and uses a low number of raw materials such as rare Worldwide, deep geothermal energy offers an enormous elements and strategic metals. In addition, unlike many other potential for future electricity and heat productions. Geother- renewable energies depending on climate conditions (wind, mal energy is sustainable, clean, and renewable, as the tapped solar, hydraulic regimes, and the water levels of rivers), the heat from an active reservoir is continuously restored [1]. geothermal energy from the Earth’s crust is nonintermittent. Geothermal energy has a number of positive characteristics: However, only a fraction of the deep geothermal energy of the 2 Geofluids the widespread deployment of deep geothermal energy to Earth’s crust is currently exploited. In practice, the exploita- tion is mostly limited to “brown fields,” i.e., areas where a res- other areas, for instance, the Limagne Graben. ervoir is already sited or where there are targeted structures Several researchers have suggested that the stresses asso- and parameters characterized that are already known [2]. ciated with faults significantly affect the flow properties and The wider development of deep geothermal resource exploi- that there is a strong correlation between the groundwater tation now requires focusing on little to unknown areas, flow and geothermal anomalies. A regional heat-flow study called “green fields.” In these fields, geothermal resource requires an analysis of the basin-scale flow distribution (due exploration (including exploratory drillings) involves a high to the strong link between the regional groundwater flow degree of uncertainty and financial risk and requires reliable and anomalies of the geothermal heat) and should account data to provide convincing arguments in the decision- for the active faulting caused by stresses. Stresses can play making process [3]. Consequently, the improvement of an important role in the setup of thermal anomalies by exploration methods is required to lower the risk and to over- favouring preferential flow paths through the rock mass [6]. come the barriers to the wider development of deep geother- For example, shear stresses on fractures have been observed mal projects. to favour the fluid circulation perpendicular to the shear To be exploited, geothermal resources require heat, fluid, direction [22, 23]. Then, stresses can favour the appearance and permeability. These favourable geothermal conditions of thermal anomalies by enabling vertical upflows of a fluid are strongly linked to the specific geodynamic context and previously heated during deep circulation [6, 24]. Thus, key the main physical transport processes, notably the stresses factors to understanding thermal anomalies are faults and permeability heterogeneities [25], which are influenced by and fluid circulations that impact heat-driving processes [4–6]. The Upper Rhine Graben (URG) is a relatively well- the stresses (a key physic) that control both the fluid circula- known example of a favourable faulted geothermal system tion (a key physic) [6, 25] and the location of exploitable hot showing thermal anomalies at the regional scale. It has been fluid upwelling. That is the reason why geothermal studies of largely investigated, providing various data such as a tectonic an extended area, such as exploration for geothermal history of the Rhine Graben [7, 8], in situ observations from resources, must consider these key factors and must be based deep boreholes, and numerical models and geochemical on the study of these physics. studies [9–15]. Moreover, a European project (the GeORG In addition to other characterization methods (geologi- project) has resulted in an important structural database cal, geophysical, and geochemical), exploration could benefit (http://www.geopotenziale.org/home?lang=3) that provides from numerical models and then a predictive analysis to locate geothermal resources. Numerical models could allow a map of seismic faults at the URG scale. The formation of geothermal resources in the URG results from the interac- an integrated evaluation of the physical processes in the tions between different coupled processes comprising URG with adapted time and space scales, considering 3D mechanical constraints, groundwater flows, and heat trans- effects. They constitute a useful tool to improve the under- fer. Large-scale structures of the graben and its associated standing of the dynamic context and may provide answers to the challenging task of geothermal exploration. Indeed, faults may act as a primary structural control on the upwell- ing geothermal fluids [16]. The rifting context and the asso- the physical conditions favouring the setup of the geother- ciated thinned crust and elevated mantle favour higher mal resources can be searched for in predictive models, temperatures at shallower depths [4, 17]. A series of thermal hence giving estimates of the so-called “favourable areas.” anomalies with temperatures above 140 C at a 2 km depth Some major phenomena were highlighted, and some local anomalies were explained using 1D hydrothermal models characterizes the URG. These anomalies are mainly located in the western part of the URG [18]. Recalling that thermal and 2D models based on continuum approaches [10, 12, anomalies are mostly affected by the regional groundwater 13] but without any hope for an industrial application to circulation [19], the high potential of the URG for thermal exploration. These models are indeed restricted to local anomalies is widely interpreted as a signature of the scales and do not consider major features such as preferen- tial pathways within fault zones or 3D effects. Stress and regional-scale fluid migrations hosted by the numerous fault zones cutting through the rock mass (multiscale fracture- numerical flow models have also been used to study geo- controlled systems, [20]). Arguing that the fault zones are thermal systems ([6, 26] and references therein), but either critically stressed and mostly under shearing conditions, the scale is not appropriate or the approach is not holistic. Baillieux et al. [21] estimated that the regional fluid migra- Thus, to the authors’ knowledge, there is currently no numerical approach available to contribute to geothermal tions could be responsible for up to 75-85% of the thermal anomalies in the URG. The radiogenic heat production due resource exploration. to the crystalline composition of the basement may explain The present work aims at demonstrating the potential of the remaining 15-25%. The geochemical data tend to confirm the combined use of 3D flow and stress models at the this assumption: Sanjuan et al. [14] conclude that the thermal regional scale (URG scale) to locate/understand the location of hidden thermal anomalies. Considering the large amount anomaly in the Soultz-sous-Forêts surroundings results from the regional circulation through a complex, but still poorly of data and knowledge, we use the URG as an application defined, system of deep faults (probably along NE-SW fault case. The results of the GeORG project constitute a relevant zones). Thus, the deep circulation is probably the origin of basis to build a structural model, which is essential to devel- the thermal anomalies and governs their geographical distri- oping the subsequent numerical physical models and then such a numerical approach. This enables us to compare butions. Their understanding is therefore fundamental for Geofluids 3 results with observed physical evidence. The main objective (b) N170 sinistral strike-slip faults reactivated from of this paper is to investigate the potential of an approach the late Miocene up to now under a NW-SE based on the use of numerical models (stress and flow compression of the graben as a strike-slip system models) to locate hidden thermal anomalies. The paper is (ii) The F2 set, ~NNW-SSE striking and E-dipping, is organized as follows: first, we give a brief overview of the also commonly interpreted as the result of the Oligo- study area. Then, we present the 3D geometry for the models cene deformation (N0/20 striking, E-dipping). according to the faulted geological setting. The major faults However, it could also be considered as different sets of the URG are considered using a DFN model (detailed in gathering the Hercynian N20 sinistral strike-slip Section 3) consistent with the geological knowledge of the and the N135 dextral strike-slip reactivated during area (see Section 2). This DFN (discrete fracture network) the late Eocene and to N170 sinistral strike-slip model is declined in 3D flow and stress models at the URG faults scale, which are, respectively, detailed in Sections 4 and 5. The results of the physical models are interpreted as func- (iii) The F3 and F4 sets, respectively, NE-SW and tions of their thermal signatures through the definition of NW-SE striking, are consistent with the Hercy- geothermal criteria related to the physics considered and nian orientations observed regionally the setting up of the geothermal resources. These criteria (iv) A fifth fault set (F5) is introduced to incorporate the enable the elaboration of indicators to use the results of the “background noise” simulation and then identify geothermally favourable areas in relation with each physic. Then, considering the strong Considering the stress regime, the sinistral strike-slip link between the stress, fluid flow, and geothermal resources, regime that started at the onset of the Miocene has been a cross-analysis of the results is realized to delineate the persisting until today, and the present-day stress state is favourable areas for geothermal resources. The results are assumed to be driven by the African-Eurasian plate conver- compared with existing thermal data at the URG scale and gence [30]. with knowledge gained through numerous studies (see Sec- As a consequence of its polyphaser tectonic history, the tion 6). Before concluding, we discuss the potential of the URG is a complex environment setting in which the rifting present work. process played a major role by providing the space for sedi- mentation, creating a series of basins, and providing a net- work of faults along the active rift margins that enabled 2. Geological Setting of the Study Area: fluid circulation and conditioned water exchanges between the URG basins [8, 14]. The transnational geological model of the URG based on boreholes and new processed seismic data The URG is an NNE-trending crustal-scale rift of the constitutes an important piece of the puzzle for the numerical European Cenozoic Rifting System (ECRIS) [8]. Its length models by providing essential data. The map of the faults at reaches approximately 300 km from Frankfurt (Germany) the top of the basement can be used to develop numerical to Basel (Switzerland), and its width is 30-40 km (Figure 1). models, and the in situ temperature data at the URG scale The URG development started from the late Eocene, mainly constitute a relevant means to assess the fault potential by the reactivation of the late Variscan, Permo-Carbonifer- (Figure 1). ous, and Mesozoic fracture systems [8]. Schumacher [8], Edel et al. [7], and Rotstein et al. [27], among others, detail the tectonic history and the resulting structure. Only some key 3. DFN Model at the URG Scale aspects are recalled here. 3.1. Why a Discrete Approach. The very first step of our work Since the Visean age during which the granitic intrusion is to build a model geometry consistent with the geology and took place, a repeatedly changing stress field has occurred the physical processes involved. The large heterogeneities in the URG, leading to the reactivation of a complex set of introduced in the physical system by the presence of the fault crustal discontinuities and the creation of a complex fault zones rule out any chance for an analytical solution, thus network in the sedimentary blanket [8]. Based on the tectonic making numerical tools necessary. As argued by Witter and history and on knowledge acquired from previous works Melosh [31], the choice of the numerical tools is a challenge [28, 29], four main fault sets corresponding to the main sta- in itself since it will impact the results. The need to account tistical models have been individualized: for real 3D effects raises another challenge, which is the avail- able panel of suitable tools: the numerical tools for solving (i) The F1 set, mainly ~N-S striking and W-dipping, is large-scale models in reasonable computation times while widely observed in the Rhine Graben and commonly offering the possibility to take into account geological struc- interpreted as the result of the Oligocene deforma- tures are few. In the case of faulted settings, the presence of tion (N0/20 striking, W-dipping). However, this natural discontinuities that can result in heterogeneous stress set could also be related to fields and localized fluid flow pathways [32, 33] has pro- (a) the lower Carboniferous to Permian N20 sinis- moted the development of fracture network models for the tral faulting reactivated during the late Eocene numerical simulation of fractured rocks [34]. Here, the key [8] or role of the regional fault network is taken into account due 4 Geofluids Worms Mannheim Mannheim United Kingdom Germany Belgium Heidelberg Speyer Heidelberg France Switzerland Landau Speyer Italy Bruchsal Wissembourg Spain Karlsruhe Landau Soultz Haguenau Bruchsal Baden-Baden Wissembourg Karlsruhe Germany Strasbourg France Soultz Obernai Offenburg Lahr Haguenau Selestat Baden-Baden Colmar Freiburg Strasbourg Mulhouse Obernai Offenburg Lurrach Switzerland Lahr km km Legend Temperature at 2 km depth Cities Value High: 167 Faults_GeORG_basement top Low: 61 Country outlines Figure 1: Location of the URG and of the case study in the URG. The map shows the faults at the top of the basement and the temperature map at 2 km depth delivering the main thermal anomalies (GeORG project). expertise on geological settings similar to the one studied (in to a discrete numerical approach (DFN) for fault zones at the regional scale that remains a challenge alone. The model to our case, rifting systems), an understanding of the tectonic represent the regional fault network is generated with the history, and the integration of generic structural tools (such 3DEC™ software [35]. as Riedel systems, see [41]). Furthermore, the numerical constraints and specificities of the codes used have to be 3.2. Geometry Building. The case study is a 130 × 150 km taken into account. area centred on the main thermal anomalies in the URG Following the above-defined concept and using the geo- (Figure 1). The lateral east-west boundaries are chosen large logical knowledge of the URG, the discontinuities of the enough to incorporate a large portion of the rift shoulders to DFN model are the geological faults observed at different limit the boundary effects in the subsequent computations (at scales, from plurikilometric to kilometric scales. The largest discontinuities that cross the model are on the regional scale least in this direction). The geometry building is focused on the conceptualization of the 3D fault zone network into a (primary faults of hundreds of kilometres known as major DFN at the studied scale; more specifically, the primary faults of the Variscan or tertiary rifting history). As disconti- source of the data is the seismic fault map of the transna- nuities split the model blocks, the resulting block sizes tional database [36], where any valuable data (e.g., seismic decrease, allowing the integration of more local-scale discon- reflection results, borehole data, outcrop traces, and scanline tinuities (faults of several kilometres from seismic profiles or window sampling) should be integrated into the structural and lineament maps). The main faults integrated in the interpretation to provide estimates on fault zone orienta- DFN model are the faults delimiting the URG and those tions, dip angles, and spatial extents [37–40]. Additional inherited from the main Variscan NE-SW trending sinistral generic geological knowledge is required to complete the shearing (Lubine and Baden-Baden faults) ([7, 8, 42], among missing information (e.g., the dips of some deep fault zones) numerous others), as shown in Figure 2(a). The secondary and to guide the conceptualization into a DFN (e.g., the fault fault network within the basement of the URG has been built zone priority rule dictating “who stops on who”). This from data of the GeORG project (within the basement and includes the following nonexhaustive list of items: structural sedimentary layers, see [36]) or from data found in the Geofluids 5 010 20 km Secondary fault network within the basement on both sides of the URG Primary fault network: Secondary fault network within the basement of the URG DFN Major fault inherited of the  GeORG fault network Variscan orogeny European digital terrain model DFN Faults delimiting the URG Lubine(LB) / Baden‐Baden fault (Variscan orogeny) Secondary fault network Fault network (a) (b) (c) Figure 2: Illustration of the DFN building method: (a) main geological structures (example of integration of structures highlighted by [8]), (b) secondary fault network from the GeORG database within the URG [36], and (c) secondary fault network in the basement on both sides of the URG from DTM analysis [90]. literature [7, 27] (Figure 2(b)). The secondary fault network The same 3D geometry will be used for the 3D flow and in the basement on both sides of the URG is from geological stress models. maps and digital terrain model (DTM) analysis (Figure 2(c)). Note that the fault network is denser in the URG than on 4. 3D Flow Model at the URG Scale the borders, as the study is focused on the URG. The dips of the faults are chosen to be equal to the fault dip at the top of Regionally moving groundwater is the basic and common the basement (estimated using the fault offset). It should be cause of a wide variety of natural processes and phenomena noted that additional geometric simplifications have been that makes the gravity-driven groundwater flow a key pro- applied to fulfil the 3DEC™ topological requirements: cess [19, 44]. The key reason for the geological scale effects of the groundwater flow is that the flow systems are ubiqui- (i) Planar surfaces are considered in 3DEC™, so the tous and active over large spectra of time and space. Among initial curved fault zones were flattened the numerous processes associated with regional groundwa- ter flow systems, the transport of heat is the most visible (ii) 3DEC™ does not accept discontinuity tips within and best understood [45, 46]. Many geothermal systems blocks, so the fault zones were extended until are associated with high densities of faults and fractures they stopped on a boundary or a previously created [25, 47–50]. Geological discontinuities such as faults and fault zone stratigraphic layers play a key role in controlling the ground- The fault zone extension requires the definition of a water flow within faulted systems [51]. priority rule to assess which fault zone stops on which one. The numerical simulation of groundwater flows has The priority rule was defined according to the tectonic become a common hydrogeological tool to investigate a vari- history first and then using a Riedel system for the contem- ety of geological settings. Significant advances in regional porary fault zones ([43] evidenced a Riedel system for the flow system analysis have been made through the application Variscan fault zones). of 3D groundwater flow models [52]. Considering that Following these main steps based on the use of the regional groundwater flows enable the approximation of the geological structures and structural expertise, the final geom- anomalies of the geothermal heat, we used a flow model to etry of the DFN is composed of 351 fault zones in the region locate the preferential discharge areas. The flow model is of interest (generated with the 3DEC™ software, Figure 3). built with the conviction that the patterns of basin-scale 6 Geofluids 130 km Mannheim Speyer Landau Karlsruhe Soultz-sous-Forêts 150 km Strasbourg 10 km Z: up Y: north X: east Figure 3: Final DFN considered in the models. groundwater flow are mainly influenced by (1) the relief of Sedimentary deposits are modelled by 3D elements. We the water table and (2) permeability heterogeneities of the used the results provided by a transnational database [36] basin’s rock framework (more specifically, permeability con- to build a geometric model of the basin. The thickness of trasts associated with faults). With that in mind, the geome- the sedimentary layers ranges between 2 km and 6 km (in try of the 3D fault network and its hydraulic connections the eastern part), with an average value of 4 km. As a result, with 3D elements (sedimentary units) need to be treated with a 3D block domain (composed of tetrahedral elements) is special care. The modelling of a regional groundwater system built from digital elevation models (DEMs) of the basement requires large sets of data, stored in various forms and scales top surface and the surface topography of the domain area. (maps, tables, spreadsheets, etc.). Thus, a GIS (geographical Note that the 3D block domain is discretized into more than information system) software was used to improve the pre- 10,000 elements (Figure 4(e)). and postprocessing of the data (ArcGIS). In the following subsections, we describe the main steps to build the flow 4.2. Hydraulic Properties. The 3D fault network (Figure 4(d)) model at the URG scale and the simulations realized to sim- and 3D block domain (sedimentary cover, Figure 4(e)) have ulate the flow paths related to geothermal resources. The 3D to be physically connected (to allow flow transfer). To do flow model was built using the 3FLO software (Itasca®), that, nodes are added at the intersection between the pipes which is designed to solve groundwater flow and transport and faces of the 3D elements. These two types of components equations in porous or fractured media (see the description (1D-3D elements) were joined in accordance with the tec- of the code in Appendix A: 3FLO Numerical Code). tonic history. Following this, we obtain a 3D model of the URG. The 4.1. 3D Geometry: DFN and 3D Volumetric Elements. The rock matrix of the basement (between faults) is considered starting point to build the 3D DFN flow model is the 3D impermeable, and the fluid flow takes place within faults. regional fault network of the URG developed in Section 3. As a first approximation, the hydraulic properties of the fault -7 2 We used information provided by the 3DEC™ model (rele- sets are uniform, with a transmissivity of 10 m /s. The sed- vant data for each fault plane: strike, dip, coordinates of the imentary cover (3D continuum elements) is considered per- mass centre, and intersections with other planes) to repro- meable and hydraulically connected to faults (pipes) in duce the same 3D network in 3FLO. In 3FLO (Appendix A: accordance with the tectonic history. The sedimentary cover 3FLO Numerical Code), the flow within faults is supported of the URG is composed of different lithology units. As a first by 1D elements (called pipes). These 1D elements are used approximation, we considered two main different units to build a 2D grid, using two sets of pipes that are each regu- where the different units are combined according to their larly spaced. Both sets are oriented perpendicularly, with ori- permeability values. The values of the physical parameters entations of 0 and 90 (with respect to the plane dip), used in the model such as the permeability were estimated respectively, forming a grid, which is generated on each flow from geophysical investigations and hydraulic tests made on different geological units [53, 54] and, for instance, have plane. According to the 3DEC structural model, the whole model domain of the URG for the flow model is a parallele- been used in the hydrothermal model of [13]. The deepest piped measuring 130 km in width (west to east), 150 km in one corresponds to Keuper, Muschelkalk, and Buntsandstein -14 2 length (north to south), and approximately 10 km in height units with a permeability of 4.10 m and a porosity of 0.17 (depending on the surface topography). (see the blue layer in Figure 4(e)). The second one (above, see Based on an analysis of the fault sets (F1 to F5) in line the gray layer in Figure 4(e)) corresponds to Jurassic, Oligo- with the knowledge of the regional tectonic history (Section cene, Miocene, Pliocene, and Quaternary units, considered -18 2 3.2 and Figure 4(d)), each fault has been assigned to a fault less permeable with a permeability of 1.10 m and a poros- set as a function of its orientation. ity of 0.15 [13]. Note that the selection of 3D elements in Geofluids 7 Geological data & 3D numerical flow model Particle tracking conceptualization Geometry Boundary conditions Steady‐state groundwater flow 3DEC DFN (a) DEM topography Input of particles (d) Initial positions of particles (g) Hydraulically Tectonic history favourable Regional fault network Resampling areas (b) 3D DFN Regionalization Discharge areas of Simulation Flow paths deep regional ground water loops Connection based on Fixed head boundary (h) tectonic history conditions Steady‐state (e) Flow path analysis xN groundwater flow 3D volumetric (c) DEMs (f) elements Elaboration of hydraulic criteria related to geothermal anomaly characteristics Sedimentary units Elaboration of indicators (i) Time (ii) Distance (iii) Depth (iv) Flow direction (v) Concentration of relevant particles Step 3 Step 4 Step 1 Step 2 Figure 4: Main steps to build the flow model and locate preferential discharge areas at the URG scale. order to apply the respective hydraulic properties is done groundwater flow system is simulated. Under steady-state conditions, the mathematical formulation of the model is using the DEMs provided by the transnational database [36], such as the position of the base and top of main geolog- highly simplified. The hydraulic conductivity is the only ical units. required parameter, contrary to other parameters such as the storability and specific yield related to transient responses, which are not necessary. 4.3. Initial and Flow Boundary Conditions. A hydraulic no- To study flow paths, a particle-tracking methodology that flow boundary condition is specified for the lateral (north, accounts for advective and dispersive transports is intro- south, east, and west) and bottom (z = −10 km) model duced (Figures 4(g) and 4(h)). Particles (over a thousand) boundaries. For the top of the model, constant hydraulic are positioned close to the border faults of the URG at the heads are imposed (fixed). The relief of the water table at top of the model that can be possibly identified as recharge the regional scale can be either topography-controlled or areas [9, 13] (see Figure 4(g)). A series of simulations with recharge-controlled [55]. Here, we assume a topography- different values for the random number generator is carried controlled water table system, which is justified by the fact out to capture the different possible flow pathways. For that the system is in a humid region with a subdued topogra- the flow system modelled, several thousand flow paths are phy and low hydraulic conductivity [56]. The current topo- constructed. The simulated time is chosen large enough to graphic elevation has been regionalized (averaged by zone), allow particles to move out of the system or to be stuck in and the corresponding hydraulic head values have been a zone. At the end of each simulation, an ensemble of flow applied as boundary conditions. As an initial condition, the paths (from recharge to discharge zones) is constructed hydraulic heads correspond to the minimum valley centre’s (Figure 4(h)). For each particle, the following characteristics elevation (100 m). are extracted at each time step: position (x, y, and z), absolute time (t), properties of the element containing the particle 4.4. Numerical Simulations to Assess Favourable Areas such as its type (1D or 3D for fault or sedimentary unit, 4.4.1. Particle-Tracking Considerations. Once the flow respectively), its permeability, the name of the domain (fault model is built, the steady-state dynamic equilibrium of the sets or lithology units), etc. 8 Geofluids distance). The aim is to highlight the particles that cover 4.4.2. Elaboration of Criteria and Indicators. To use the results of the particle-tracking method (flow paths) for geo- greater distances. A weighting factor of 2 is attributed to each thermal resource characterization, we need to choose rele- particle covering a distance of more than three-quarter of the maximal distance, which represents around 1% of the vant hydraulic criteria that constitute a good signature of the geothermal anomalies. Then, based on these main char- total particles. A weighting factor equal to 1 is attributed acteristics of geothermal anomalies (hydraulic criteria), indi- to each particle covering a distance of more than one-half cators are elaborated and extracted from each particle and less than three-quarter of the maximal distance, which trajectory (Figure 4-Step 3) to use the results of the simula- corresponds to 14.5% of total particles. By the same reason- ing, two threshold values are calculated related to the max- tion that enable the identification of the geothermally favour- able areas. The following assumptions are used: imal travel time; then, a weighting factor of 2 is attributed to each particle with a travel time of more than three-quarter (1) Geothermal anomalies result from a state of imbal- of the maximal travel time, and a weighting factor equal ance over a long period where the only efficient to 1 is attributed to each particle having a travel time of mechanism (ubiquitous and active over large spectra more than half and less than three-quarter of the maximal of time and space) is the regional groundwater flow. travel time. For the indicator related to the depth reached Then, the flow path length (long enough to be by the particle along its trajectory, we considered that a par- regional) and the time required to cover this distance ticle, which has reached a sufficient depth of more than are used as the two main indicators to highlight 2.5 km (related to the average geothermal gradient), is inter- regional groundwater loops esting and deserves a weighting factor. Then, a weighting factor of 1 is attributed to each particle that reached a depth (2) Geothermal anomalies reflect positive thermal anom- of more than 2.5 km along its trajectory. To finish, the sum alies resulting from upward flow in discharge areas of weighting factors is calculated for each particle. The idea due to ascending warm waters. Therefore, the particle is that particles with the higher scoring are likely to produce must come from a deeper part and has to reach a suf- thermal anomalies and correspond to the more relevant ficient depth along its trajectory to carry warm(er) particles and inversely, the lower scoring to the less relevant water until the upflow area. The depth reached by particles. To eventually constitute a preferential target area the particle along its trajectory, the final position (or (thermal anomaly), a final requirement has to be fulfilled: a position slightly upstream from it) of the particle asufficient density of relevant particles must be localized along its track, and the flow direction are thus also in a “restricted” spatial area. used as main indicators 4.5. Favourable Area Estimates: The Flow Model. Figure 5(a) As mentioned just above, the indicators that are used presents the results of the flow path analysis at the URG scale are (1) the total distance travelled by the particle from where particles are plotted at a 2 km depth (equal to that of the recharge area (starting point) to the discharge area the available temperature map, see Section 6.2) as a function (ending point) during time t, (2) the time spent by the of the weighting factors. Particles with a longer distance and particle in the system to travel between the recharge and time (in accordance with regional groundwater loops) and discharge areas, and (3) the depth reached by the particle reaching a depth over 2 km along the trajectory (in accor- throughout the trajectory. Note that the relevant informa- dance with the temperature) are selected. The target areas tion associated with each track is extracted to be analysed are then delineated by defining the contours encompassing in a GIS software. these relevant particles (Figure 5(a), green lines). Five differ- 4.4.3. Flow Path Analysis. The main objective of the ground- ent preferential areas are highlighted from the results of the water flow path analysis is to cross-correlate the properties of 3D flow model, which are distributed as follows: two areas each particle trajectory (distance, time, depth, localization, to the north of the URG (surrounding Speyer and between and flow direction) to locate the preferential discharge areas Landau and Mannheim), one area clearly located to the (with specific features of geothermal anomalies) of regional south of Soultz-sous-Forêts, and two areas located near groundwater loops. and south of Strasbourg (Figure 5(a)). Figure 5(b) presents After an individual flow path analysis, a spatial and qual- the trajectories associated with the most relevant particles itative analysis is realized to localize the preferential dis- providing the preferential areas. The northern areas are charge areas of regional groundwater loops. A preferential mainly related to particles from the eastern shoulder of the discharge area corresponds to a zone with a high concentra- URG and along a SE-NW direction. The second area to tion of relevant particles (related to the deep regional flow the south of Soultz-sous-Forêts is mainly related to particles circulation), which is likely to present thermal anomalies from the western edge of the URG and moving in a NW-SE (Figure 4-Step 4). To select relevant particles, a preliminary direction. Note that few particles come from the eastern edge data analysis is performed. The statistical distribution is and in a SE-NW direction. The southernmost areas (sur- observed to define threshold values that are used to assign rounding Strasbourg) are related to particles from the east- weighting factors to particles. For example, to select particles ern and western edges. that correspond to regional groundwater loops (related to a The results of the 3D flow model provide preferential distance indicator), two threshold values are determined discharge areas related to geothermal resources. Since the (equal to one-half and three-quarter of the maximal ultimate aim is to support expensive exploration campaigns, Geofluids 9 Mannheim Speyer Landau Soultz-sous-Forêts Strasbourg Classification of particles Trajectories associated with most relevant particles (see from above) Main flow directions related to relevant particles More Less relevant relevant Favourable area km Fault network (a) (b) Figure 5: (a) Classification of particles resulting from the flow path analysis based on selection criteria and the respective delineated preferential target areas obtained (favourable areas are delineated by the green circles). Note that the particles are plotted at a 2 km depth upstream from the ending point of their trajectory. (b). Particle trajectories of the most relevant particles (red crosses) and the associated main flow directions. any information improving the confidence in the area loca- ([29, 58, 59], among numerous others) and the segmentation tion is valuable. When no in situ data is available, cross- of the medium (the more fault zones, the more deviation). validating the interpretation with results from other models To build a 3D stress model of the URG, we use the code is a way to increase the confidence level. In this study, we pro- 3DEC™ [35]. pose to build a 3D stress model to extrapolate preferential areas from the mechanical perspective and to finally compare 5.1. Mechanical Properties. A full description of the code and its constitutive equations is provided in Appendix B: them to the flow model estimates. 3DEC Numerical Code, and only a short outline is given here. Since we focus on the presence of the fault zones, 5. 3D Stress Model at the URG Scale the matrix behaviour is kept simple and assumed to be elas- tic (Table 1). The fault zones, however, follow an elastoplas- tic law (Coulomb slip), with a possible dilation (Table 2). In We propose to investigate the heterogeneities in the stress distribution that may favour the setup of thermal anoma- addition, we assume that any large deformation occurring lies. The stress distribution in rock masses is significantly within the medium is localized into the joints and not into affected by structural and material heterogeneities (including the rock matrix (see Appendix B: 3DEC Numerical Code mechanical behaviours/parameters and fault zones) and for more details). Due to the absence of data at the considered scale, the anthropic activity. Yale [57] gives a review of the in situ evi- dence of fault impact on the local deviation of stress from the parameters are extrapolated from estimates obtained by fit- far-field trend. The conclusions are that the main factors ting numerical models on in situ experiments [60, 61]. We locally affecting the stress are the proximity to the fault zones extrapolated the parameters qualitatively considering fault 10 Geofluids Table 1: Mechanical parameters for the rock matrices. 5800000 Notation Value Value Parameter (unit) (sediment) (basement) −3 –1 Density ρ (kg·m ) 2400 2680 0.2 mm·y Bulk modulus K (GPa) 16.7 42.8 Shear G (GPa) 7.7 22.1 modulus Table 2: Mechanical parameters for the fault zones. Parameter Notation (unit) Value −1 Normal stiffness k (GPa·m)10 −1 Shear stiffness k (GPa·m)5 Cohesion c (GPa) 0 –1 Internal friction angle ϕ ()15 0.2 mm·y Dilation angle ψ ()5 –1 0.56 mm·y Critical shear displacement u (m) 1 200000 300000 400000 500000 600000 Easting (UTM, m) zones with noncohesive infills (weaker elastic modulus, zero Figure 6: Far-field displacements and embedding model (grey) for cohesion) and large wavelength profiles (smaller dilation setting the faulted model (red frame) boundary conditions. angle, larger critical shear displacement) that are altered by Modified after Buchmann and Connolly [62]. UTM: Universal numerous tectonic stages (smaller friction). The parameters Transverse Mercator. are on the orders of magnitude of the parameters from models built at similar scales in the URG [62, 63]. are applied to the model boundaries, while the gravity is still 5.2. Boundary Conditions. The displacement boundary con- active within the model: ditions are set in accordance with the present-day tectonic context of the URG and are derived from the model of Buchmann and Connolly [62]. The faulted geometry pre- z = − ρgdz , v d sented in Section 3 is embedded in a continuous material 0 having the dimensions of Buchmann and Connolly’s model S z = S z = S z , H max d H min d V d (Figure 6). Following a procedure used by Peters [64], this 1 − v embedding enables simulating the far-field nature of the boundary conditions. As a first approximation, the embed- where S (Pa), S (Pa), and S (Pa) are the vertical, V H max H min ding material is assumed to have the same properties as maximum horizontal, and minimum horizontal principal −2 the crystalline basement. In addition to these lateral bound- stresses (respectively); z (m) is the depth; g (m·s )isthe −3 ary conditions, the normal displacement is fixed at the vertical component of the gravity acceleration; ρ (kg·m )is model base (roller-type boundary condition), and the top the rock mass density, and ν (−) is Poisson’s ratio. is left free. 5.4. Simulations to Assess Favourable Areas. Once the 5.3. Initial State. In addition to the present-day forces, the assembly of blocks and joints is balanced (see Appendix stress field is affected by past tectonic events that leave B: 3DEC Numerical Code for more details), the model residual stresses through the rock mass. The prestressing gives an estimation of the mechanical equilibrium through method is commonly used to simulate the impact of past the system. In particular, displacement and stress heteroge- events [62, 65–67]. Assuming that the previous tectonic epi- neities are obtained and can be inspected to highlight areas sode occurred a sufficiently long time ago to allow the relax- where the mechanical state can favour the presence of ther- ation of the medium, the initial state can be approximated mal anomalies. As previously, this requires the definition of using a gravitational loading. Gravity is activated within the an appropriate indicator. The definition of a mechanical blocks, and a uniaxial stress condition (equation (1)) is pre- indicator is difficult since stresses can act on flows in sev- scribed on the lateral boundaries of the model to account eral ways. However, since stresses can act on the ground- for the influence of the surrounding rock mass. Note that water heads directly (impact on pressure gradients) or on due to the high density of the fault zone network, the initial the pore/fracture network (impact on medium effective condition makes the stress distribution very heterogeneous permeability), most mechanical indicators can be catego- through the rock mass. Once the initial equilibrium is rized following the volumetric and deviatoric decomposi- reached, the conditions described in the previous section tion of the stress tensor: Northing (UTM, m) Geofluids 11 Geological data & 3D numerical geomechanical model conceptualization Geometry Initial stress state 3DEC DFN Elaboration of stress criteria Gravity balance related to geothermal anomaly characteristics Elaboration of indicators - Mean stress Tectonic history Regional fault network 3D DFN Boundary conditions Regional tectonic trend Result analysis using Geomechanically aforedefined indicators favourable areas Connection based on tectonic history 3D continuum Results of the stress model DEMs Sediment pile and basement Figure 7: Main steps of the geomechanical simulations (from data to delineation of preferential areas). (i) The volumetric part of the stress tensor (e.g., mean Low values of the indicators are thus used to identify favour- stress) is directly linked with the groundwater heads: able areas. flow will happen from the more compressed areas to The different steps to delineate the favourable areas with the less compressed zones [68]. “Volumetric” indica- the geomechanical model are summarized in Figure 7. tors could also relate to modifications of the pore net- 5.5. Favourable Area Estimates: The Stress Model. The equi- work (e.g., fracturing), but since we consider hard librium of the mechanical system is computed through rocks under compressive states, the pore network 3DEC™. We obtain a stress map that can be used for delin- modification will mostly happen under the deviatoric eating the preferential target areas. Prior to postprocessing influence the stress, the depth at which we search for favourable areas (ii) The deviatoric part of the stress tensor does not must be set. Rather than a fixed depth, we decided to inter- affect the pressure heads but can affect the effective pret the results close to the sediment/basement interface. permeability of the medium. In the rock matrix, Since the DFN was built according to the fault zone traces the stress deviator can be responsible for shear- observed at the top of the basement, the level of uncertainty induced fractures that greatly contribute to the in the geometry is more limited close to the interface. The permeability increase of the medium. For fault mean stress distribution is presented in Figure 8(a) (the zones zones, deviatoric stresses result in shear displace- located in the horst are ignored). ments perpendicular to which upflow is favoured Five different preferential areas are highlighted from the [24]. In addition, the shear stress maintains the results of the 3D stress model. They are distributed as follows: fault zones in an active state, which is argued to two areas in the northwestern part of the URG (to the west reduce the fault sealing and hence favour the per- of Speyer), one area located around Soultz-sous-Forêts, and meability [69] two areas located in the southern part (one to the south- west of Strasbourg and one to the east) and to the south In the first study, we selected a single scalar indicator of Strasbourg (Figure 8). to be chosen according to either the volumetric or the deviatoric part of the stress tensor. Denoting σ (Pa) the 6. From Regional-Scale Models towards the stress tensor (see Appendix B: 3DEC Numerical Code for more details), we chose the mean stress σ = Tr σ /3 as Prediction of Favourable Areas for the mechanical indicator: by releasing the efforts acting Geothermal Resources on the rock mass, the less compressed areas (i.e., low σ ) 6.1. Cross-Analysis of Favourable Areas. The final step of the facilitate fluid discharges, including vertical upflows of fluids that were possibly heated up to the rock mass below. work is to estimate areas where the physical conditions ++ 12 Geofluids σ compared to the existing thermal map of the URG at a 2 km depth. From north to south, three main local thermal Mannheim anomalies can be distinguished from the temperature map: one in the north, between Landau and Speyer; one in the cen- Speyer tre, around Soultz-sous-Forêts; and one in the south, to the west of Strasbourg (Figure 9(b)). Then, from ~east to west, Landau the temperature map shows that the highest temperatures are preferentially located in the western half of the rift. The major local thermal anomaly is located around Soultz-sous- Soultz-sous-Forêts Forêts in the centre of the URG, where the highest tempera- tures are located to the east of Soultz-sous-Forêts and then in the centre of the basin. The second most important local thermal anomaly is located between Landau and Speyer, Strasbourg where the temperature decreases from south (Landau) to north (Speyer), and the third local thermal anomaly is located to the west of Strasbourg (Figure 9(b)). From a global perspective, the cross-analysis of the results of the stress and flow models delineates three distinct Fault network km favourable areas: one in the north (west of Speyer), one in the centre (south of Soultz-sous-Forêts), and one in the south Figure 8: Mean stress distribution in the model and the delineated (south of Strasbourg). This distinct distribution (from north preferential target areas obtained. to south) in three main areas is in good agreement with the distribution of the main local thermal anomalies known in the URG (Figure 9). favour the setup of geothermal resources. Favourable areas were delineated for both key physics (flow and stress models) based on related selection criteria defined in accordance with 6.3. Discussion of the Cross-Analysis Results the geothermal characteristics and the physics considered. The most favourable areas are assumed to be at the intersec- 6.3.1. Accuracy of the Areas Delineated through the Cross- tion of flow-favourable and stress-favourable areas. They are Analysis. From a global perspective, the results obtained from estimated qualitatively by overlapping the different favour- the cross-analysis are very encouraging, with the delineation able areas obtained with each model (flow and stress models) of three favourable areas distributed from north to south (Figure 9(a)). along the rift, which is in good agreement with the distribu- The cross-analysis highlights three main favourable tion of the main local thermal anomalies known in the areas: one is in the northern part of the URG, close to Speyer; URG. However, the accuracy of the results against the ther- a second one is to the south of Soultz-sous-Forêts, and a third mal map varies from one delineated area to another. From one is to the south of Strasbourg (composed of three different north to south: subareas) (Figure 9(a)). The northern delineated area (west of Speyer) results from the overlapping of the two main (i) the area between Landau and Speyer partially high- favourable areas of the stress (blue circle orientated in the lights the local thermal anomaly. The cross-analysis east-west direction) and flow (green circle orientated in the results delineate an area to the west of Speyer, NS direction) models. Close to Soultz-sous-Forêts, the two approximately 20 kilometres away from the in situ distinct areas (blue circle orientated in a NE-SW direction area and green circle orientated in a NW-SE) result in one main (ii) the area close to Soultz-sous-Forêts is the best favourable area. Note that these two main favourable areas correlated one by the model. Still, the final estima- (west of Speyer and south of Soultz) are located in the west- tion is to the south of Soultz, while the actual area ern part of the rift (Figure 9(a)). For the four areas (two is most visible to its east respective areas for each model corresponding to two green and two blue circles, respectively) located to the south of (iii) the model results for the area close to Strasbourg are Strasbourg, the cross-analysis results in three subareas shifted eastwards compared to the thermal map. (including two small ones), where the largest one is orien- However, as opposed to the two other areas, the in tated north-south in the western part of the URG situ data in this area is much more diffuse. The (Figure 9(a)). To summarize, the cross-analysis highlights model results follow the natural tendency of the two main favourable areas that are located in the western part thermal anomalies, and in this regard, we conclude of the URG, to the south of Soultz-sous-Forêts, and to the that this area is not the least precise one west of Speyer. A third area can be distinguished to the south of Strasbourg, where three subareas are delineated. The uneven precision of the model results steams from the cross-analysis of the stress and flow results. While the 6.2. Delineated Favourable Areas versus the Thermal Map. two hydraulically favourable areas tend to correlate with the actual thermal anomaly located between Landau and Speyer The main favourable areas from the cross-analysis are Compression Geofluids 13 Mannheim Mannheim Speyer Speyer Landau Landau Soultz-sous-Forêts Soultz-sous-Forêts Strasbourg Strasbourg (Temperature map from GeORG) Temperature at 2 km depth Stress Flow Value High: 167 Cross-analysis Low: 61 km Fault network (a) (b) Figure 9: (a) The results of the cross-analysis based on the overlapping of the respective favourable areas issued from the stress and flow models. (b) Map of temperature at a 2 km depth obtained from the GeORG project. (Figure 9), the mechanical results are less precise and move reasonably well estimated from a regional perspective. From the final area estimation northwards. In the end, the northern the 3D flow model, this favourable area results from particles delineated area partially highlights the local thermal anom- originating mainly from the western edge of the URG aly. Conversely, the mechanical results appear to better high- (located to the northwest of Soultz-sous-Forêts, see the large blue circle Figure 10(a)) and in a main NW-SE flow direction light the thermal anomaly located on the western edge of the rift than the flow model’s ones. According to the 3D flow (Figure 5(b)). The 3D cross-sectional view (Figure 10(b)) model results, these hydraulically favourable areas are related shows that particles move from the west (downward flow) to particles originating from the eastern border (see yellow along the URG east-dipping boundary fault to the east circle, Figure 10(b)) and are characterized by regional flow (upward flow), mainly along west-dipping faults and associ- ated subvertical faults that are connected to these west- pathways that mainly occur through NW-SE faults (F4 fault set) in a global direction SE to NW (Figure 5(b)). These few dipping faults (Figure 10(b)). These results suggest deep particles go through the centre of the URG (where the sedi- regional groundwater loops that occur through east-dipping mentary basin reaches more than 4 km deep) and then (downward flow) and west-dipping (upward flow) faults that migrate through the regional fault network up to the favour- intersect at great depth and along the main NW-SE flow direction and where the main recharge area is located to the able areas. These deep regional flow paths appear consistent with the geochemical evidence [14]. Indeed, the geochemical northwest of Soultz-sous-Forêts. These preferential flow signatures of geothermal brines collected from the granite pathways along the west-dipping fault zones characterize basement (Landau, Rittershoffen, Soultz, and Insheim at the horst structures and are consistent with some previous works graben’s NW borders) suggest that they all reacted with deep [9, 70]. Moreover, according to the results (particle flow paths), we can observe deep flow paths (few particles) from sedimentary rock at temperatures close to 225 Cat depths ≥ 4 km [14]. This signature involves a migration of geother- the SE (small blue circle in Figure 10(a)) and in a global mal brines from the centre of the URG to the graben’sNW SE-NW flow direction emerging close to the Soultz and borders through a complex system of deep faults. Rittershoffen areas. These few particles go through the centre For the centre area, the delineated area (close to Soultz) of the URG (where the depth of the sedimentary basin reaches more than 4 km) and then migrate through the that corresponds to the major thermal anomaly known is 14 Geofluids a. Areas of origin for main favorable areas km Fault network South of Soultz Favourable areas Between Mannheim and Speyer (hydraulic) Around Strasbourg b. Soultz Temperature at 2 km depth Trajectories associated with relevant particles Main flow direction Value Fault network High: 167 Start/middle/end Low: 61 Figure 10: (a) Plan view of trajectories associated with relevant particles and areas of origin of particles associated with favourable areas. (b) 3D cross-view (south-side view) of trajectories in the favourable area located close to Soultz-sous-Forêts according to the blue rectangle drawn in (a). regional fault network up to the favourable area (which could analysis because it presents the key characteristics from a also be consistent with the results of [14]). The favourable hydraulic point of view and also from a mechanical point area close to Soultz is the best located one by the cross- of view. Indeed, the results of the mechanical model show a Geofluids 15 –4 network affects the redistribution of stresses and regional groundwater fluid circulation. Additionally, the stresses sig- Δu (m) nificantly affect the flow properties, and there is a strong cor- relation between the regional groundwater flow and geothermal anomalies. The groundwater flow paths are Opening mainly influenced by the regional fault network (structural assemblages and heterogeneity of permeability created by faults) as well as by the relief of the water table. Thus, the key factors are the faults and permeability heterogeneities [25], which are influenced by the stresses (a key physic) that control both the fluid circulation (a key physic) [6, 25] and the location of hot fluid upwelling. The respective favourable areas delineated by the flow and stress models demonstrate the strong link between the deep regional fluid circulation Closing within the regional fault network, the stresses, and the distri- bution of the geothermal resources. These previous observa- tions lead us to ask the following question: is it essential to –4 –10 consider the thermal physical processes to highlight the ther- Stress Opening mal anomalies in the context of geothermal exploration where the ultimate goal is to limit the area to explore and Figure 11: Normal displacement variations of fault zones and reduce the associated costs? Building a hydrothermal model corresponding favourable areas (in purple). Blue circles are a at the regional scale could facilitate the interpretation of the reminder from the mean stress interpretation. results by providing the distribution of the temperature and consequently local thermal anomalies. This type of model could serve to avoid approximations that can appear through less compressed region (Figure 8(b)) where the stress state the definition of criteria and then have an impact on the can favour upflow. Thus, the major thermal anomaly delineation of the favourable areas. Thermal processes such known in the URG is the best located one by the numer- as the heat transfer between faults and the surrounding ical approach due to the combination of numerous key sit- matrix and at the intersection of the base of the thermal blan- uations: structurally, hydraulically, and mechanically. ket and the fault (in which the upward transfer of heat Structurally, this region corresponds to a horst with a occurs) result in redistributions of heat. These redistributions west-dipping fault (enabling the intersection at depth with of heat could provide slightly shifted favourable areas in com- east-dipping faults). Hydraulically, this area corresponds to parison with those currently delineated in this study. Thus, a preferential discharge area of deep regional fluid circula- this type of numerical model could be relevant and, more spe- tion. Geomechanically, the mechanical stress in this region cifically, useful to move towards a better overall understand- favours fluid circulation and then hot fluid upwelling in ing of the URG geothermal systems and should provide some the rock matrix and fault. constraints on the relative influence of each physics. In any To finish, for the southern area, the local thermal anom- case, building a hydrothermal model at the regional scale aly is particularly well delineated by the geomechanical constitutes a real challenge and then requires further analysis. model (less compressed area and faults exhibiting opening The promising results also demonstrate the key role of the tendencies, see Figure 8). Second, preferential pathways regional fault network, which is at the basis of the two numer- (from the western border fault) that occur along the west- ical physical models. Then, a better delineation of the favour- dipping fault zone also characterize the southern area (to able areas requires an improvement of each numerical model the west of Strasbourg), and we can observe particles from and thereby an improvement of the DFN geometry. the eastern border fault that emerge in this favourable area (see red circles, Figure 10(a)). All of this information seems (2) Cross-Analysis versus Coupling. One of the main assump- consistent with the knowledge provided by the numerous tions of the study presented in this paper is that the estima- studies realized within the URG (geophysics, geochemistry, tion of favourable areas can be obtained by cross-analysing etc.). The results demonstrate the potential of this type of results from separate single-physic models. Several authors numerical approach to provide constraints for the challenge use either one-way coupled or fully coupled models as an of geothermal exploration and to move towards a better over- attempt to more accurately address the physical features all understanding of the URG geothermal systems. involved in large-scale phenomena (>kilometric characteris- tic length) [71–75]. 6.3.2. Relevancy of the Model Cross-Analysis In our case, a solution for achieving one-way coupling (1) Choice of the Physics Described by the Models. The pro- would be to use mechanical results for parameterizing the posed approach is based on the choice of the single physics groundwater flow model. Indeed, mechanical results predict described by the models. Here, we propose an approach huge heterogeneities in stress distribution, including in its deviatoric part (which can be related to second-order based on these main underlying hypotheses: the regional fault 16 Geofluids carried out by hierarchizing them based on their impact on fracturation). As a result, the rock mass is expected to exhibit heterogeneous permeabilities (increase of permeability with the results. Most notably, the boundary conditions and the fracturation). More complex effects could also occur in the use of nonheterogeneous parameters for the groundwater flow and mechanical models, respectively, are expected to fault zones such as permeability decreases in fault zones that are less shear active [69] or increases in permeability in the be mostly responsible for the uneven accuracy of the areas direction perpendicular to the shear displacement [22, 24]. they delineate (see Section 6.3.1.). Full coupling for the considered geometry (large scale Nevertheless, the reasonable agreement of the model and DFN) seems, however, to be unrealistic with the tools results with the thermal observations leads us to believe that the physical model accuracy, both in terms of state equa- at hand today. The fully coupled models found in the litera- ture consider much more reduced scales [74, 75]. Conversely, tions and parameters, is a secondary aspect in the approach. models exhibiting scales comparable to ours (~100 km char- Rather, the geometry description seems to be of more acteristic length) incorporate either only a single physics importance, as is usually the case with models based on a and/or unfaulted media [26, 62, 65, 76]. CPU- or GPU- dense DFN. based parallel simulations show promising results for han- dling fully coupled models but need further progression to 6.4.2. Geometry Description. The DFN geometry (regional reach the scale and the medium segmentation considered in fault network) requires a more detailed description. The geo- this study [71–73]. dynamic evolution of the rift system and existing conceptual At this stage of development of the study, we anticipate models of the rifting zones are important sources of informa- tion to improve the DFN construction [77–79]. The interpre- the use of coupled models to introduce more uncertainties than improvements in the overall approach. Given the tation of the two deep seismic lines illustrates the structural absence of information on coupled phenomena at large spa- asymmetry of the Rhine Graben, which is characterized by tial and time scales, the coupling laws extrapolated from thickness variations of the tertiary deposits and the presence smaller scales might be erroneous, and the estimation of the of boundary faults in the north along the eastern flank and in the south along the western flank. The asymmetric evolution coupling parameters would be very difficult (even qualita- tively). In contrast, focusing on the prevailing processes of the rift from north to south that affects the deep structural related to the setup of geothermal resources for each model features has to be considered. Similarly, horizontal bedding is expected to better constrain the results of the approach planes and a simple conjugate fault pattern characterize the (i.e., positions of the favourable areas). In our opinion, the Cenozoic fill in the north, while in the south, a more complex structural pattern exists, with the presence of tilted fault cross-correlation of relatively well-trusted models brings more objectivity in the results than estimation stemming blocks from west to east, conjugate normal faults, and a from a coupled but more uncertain model. The case study smooth anticlinal structure. Then, the steeply dipping nor- mal fault that controls the western margin of the Cenozoic presented in this paper proves that the cross-analysis tend to estimate reasonably well the thermal map obtained from graben and merges at the midcrustal levels provides con- straints on the fault system and its evolution with the depth. in situ data and demonstrates how the overlapping results from uncoupled models can be used to refine the final esti- Building the geometry by integrating this type of constraint mates of the favourable areas. could directly affect the connectivity of the medium as well as the stress redistribution due to the model segmentation 6.4. Perspectives of Improvement for the Proposed Approach. and consequently the groundwater flow paths and the geo- graphical distribution of the discharge areas. Even if there is acceptable agreement between the global distribution of the preferential areas delineated and the major thermal anomalies, we can notice some differences 6.4.3. Result Analyses. The final estimations of the proposed approach depend highly on the way in which the model with the temperature map at a 2 km depth, which underlines the need for improvements that will constitute the next results are interpreted. First, considering the two models sep- steps of this work. Based on the results and discussion, we arately, the relation between the model results and the char- conclude that improvements can be made at three main acteristics of the geothermal resource rely on the definition of levels: the geometry description, model precision, and result the indicators for each model. Recalling that a particle- tracking method is used for the groundwater flow, the (cross-)analyses. From this perspective, the following para- graphs detail the opportunities for improvement and their retained indicators are the residence time, the travelled dis- relative importance. tance, the depth reached by the particles, the flow direction, and the concentration of relevant particles. The first three 6.4.1. Physical Model Formulations. Both the groundwater indicators are linked with the temperature (fluid heating up the rock mass), the fourth one indicates resource upwelling, flow and mechanical models rely on several simplifying assumptions and could be refined in several ways (equations, and the final one relates to the discharge areas. Because of parameters, boundary conditions, and initialstates). However, the rather straightforward effect of flows on the thermal building almost exact models is an unachievable task given anomalies, we are rather confident in the chosen indicator. the large spatial and time scales involved in light of the However, the relative interpretation of the particles could absence of parameter measurements and evolution of tec- be improved through a rigorous statistical treatment. In this tonic states and related mechanical boundary conditions paper, and as a first approximation, we chose to set the [8]. Targeting a limited number of improvements could be weighting factors qualitatively to focus on the overall Geofluids 17 systems might have to be set and run to progressively refine approach. Regarding the mechanical model, the impact on the thermal anomalies is only secondhand, through affecting the estimation precision. If the proposed approach is used the flow and/or thermal processes. In this paper, we retained as a support for an exploration campaign whose aim is to site the mean stress σ , arguing that less compressed areas favour exploration wells, the results of a first region-scaled model fluid migration. The mean stress seemed to be a reliable indi- might be too diffuse or uncertain by, e.g., providing several distinct areas in the same region. If more precise results are cator for identifying the true thermal anomalies. We remind here that, in accordance with the whole approach, the expected (e.g., if the range of the final estimates is too large mechanical criterion can only be interpreted qualitatively, or if there are concerns about the locations of the favourable i.e., so as to highlight less compressed areas close to more areas), steps one to four could be repeated on progressively compressed ones. In such cases, fluid migration is assumed smaller-scaled geometries until the desired size is achieved or more confidence is obtained for the final estimates of the to move towards less compressed areas, then move upwards due to stress reduction when getting closer to the Earth’s sur- favourable areas. In this case, the results from larger-scaled face. In our opinion, a quantitative analysis is not realistic at models can provide boundary conditions for the smaller- this stage given the lack of information on the rock mass scaled ones (see, e.g., [76]). properties at the considered scale (law of behaviour and asso- ciated parameters, behaviour heterogeneities, second-order 7. Conclusion fracturation…). σ cannot be guaranteed to be the best one. Other indicators were tested such as the preferentially This paper is a current overview of our team’sefforts to inte- opened fault zones. Taking advantage of the discrete nature grate the impacts of the key physics as well as key factors con- of the model, the opening tendency can be estimated through trolling the geothermal anomalies in a fault-controlled the fault zone normal displacement u (m). More specifically, geological setting in 3D physically consistent models at the and to focus on the impact of the present-day stress state, the geological system scale. The study relies on the building of init displacement increment can be plotted: Δu the first 3D numerical flow (using the discrete-continuum = u − u , n n n init where u (m) is the normal displacement due to the initial method) and mechanical uncoupled models (using the dis- equilibrium. Figure 11 presents the contours of the fault zone tinct element method) at the URG scale (faulted geological setting). Studying the interactions of these main processes averages Δu = 1/s ∬ Δu dS (S (m ) being the area of the n n fault zone) in order to interpret the tendencies at the fault while coupling them with the effects of the regional structural elements (fault zones) helps us to understand their relative zone scale. NW-SE faults exhibit opening tendencies in con- trast to EW and NE-SW faults that exhibit a tendency to impacts on the distribution of geothermal resources. Here, close. Estimated favourable areas tend to be more widespread we applied this approach in the URG context (Figure 12), where available data (such as seismic lines) allow the building (light blue circles in Figure 11). Interestingly, the results in the surroundings of Soultz- of these types of models and a direct assessment of the results sous-Forêts correlate well with the σ (using, for instance, the temperature map). Based on the fact interpretation in this region, tending to reinforce the idea that this region com- that there exist strong links between the stress, fluid flow, and bines numerous key situations. Comparison of different geothermal resources, we followed these main steps. First, the key role of the regional fault network is taken mechanical results highlights the difficulty in choosing a mechanical criterion, which in the end might be a cross- into account using a discrete numerical approach. The geom- correlation between rock matrix and fault zone criteria. On etry building is focused on the conceptualization of the 3D a broader perspective, these results question how the fault zone network based on structural interpretation and mechanical results should be interpreted: in a straightfor- generic geological concepts and is consistent with the geolog- ical knowledge of the URG. This major step constitutes a ward manner (as is done in this paper) or to parameterize the flow model? Thorough investigations must be carried challenge alone and marks the beginning of new opportuni- out in the future to shed light on these interrogations. ties to understand the key effect of the regional fault network Second, the cross-analysis of the model results itself could on the geothermal resources. Indeed, this question remains be improved. Here, as a first approximation, we proposed to an open question and still constitutes a challenge for geother- set the favourable areas by choosing the overlapped regions mal engineers, who attempt to explain hydrothermal pro- that were qualitatively delineated. In the studies found in cesses in this faulted geological context. the literature (see, e.g., [80, 81]), the authors propose to Then, we decline the DFN model of the regional fault net- quantitatively establish favourability maps by interpolating work to build the first 3D flow and stress models of the URG in situ data and using weighting factors according to the using adapted numerical tools that explicitly consider the key measurement confidence. In our case, because the physical role of the faults. The respective results issued from each models are built under several assumptions, it seems slightly numerical model are indirectly interpreted through the defi- too early to propose quantitative favourability maps. The nition of criteria (using the link between the physics consid- more precise data at our disposal are the structural data. ered and the setting up of geothermal resources) that enable A first step could be to weight the results according to the elaboration of indicators to identify geothermally favour- how close (and thus less extrapolated) they are to our con- able areas related to the considered physics. For the flow fident structures. model, the hydraulic criteria are related to the preferential Finally, depending on the ratio between the obtained and discharge areas of the deep regional fluid circulation, and desired scales for the preferential target areas, additional the indicators used are related to the properties of the particle 18 Geofluids Available geological data Geological maps Conceptualization 3D numerical models Tectonic Results Interpretation Geological profiles history 3D flow Favourable model areas Elaboration of indicators Seismic Elaboration of criteria Favourable related to geothermal Cross-analysis areas Faults anomalies & physics 3D stress Borehole data model Elaboration of indicators Favourable Layers areas Digital elevation Results Interpretation models Figure 12: Methodology employed to delineate preferential target areas for geothermal resources using numerical models. Equilibrium of blocks F u ext Block displacements Joint forces affect block affect joint external forces Joint constitutive equations u F ext Figure 13: DEM solution scheme at each time step: the solver loops until u and F are balanced (inspired by Itasca, 2016). ext trajectories (distance, depth, etc.). For the stress model, the considered (stress and fluid flow) in the geothermal resource stress criteria are related to the mechanical state that can setup, and the method to interpret the results (criteria and favour the presence of thermal anomalies, and then, the associated indicators) are relevant to capturing the thermal related indicator chosen is the stress tensor (mean stress). processes and then can help to locate thermal anomalies. In As a final step, we proposed to compare the hydraulically other words, the geographical distribution of the geothermal and geomechanically favourable areas to locate thermal resources is strongly correlated with the preferential dis- anomalies. A cross-analysis of the results of the 3D flow charge areas of the deep regional fluid circulation occurring and stress models within the URG provides three preferential within the regional fault network and related to the less com- target areas for geothermal projects. These preferential areas pressed areas that favour these upwellings of deep hot fluids. are distributed along the rift as follows from north to south: Moreover, a detailed analysis of each favourable area to the west of Speyer, to the south of Soultz-sous-Forêts, demonstrates the ability of this type of approach to move and to the south of Strasbourg. The flow and stress models towards a better overall understanding of the URG geother- result on a common, rather condensed, overlap for two of mal systems. For example, the best favourable area delin- the three areas (Speyer and Soultz-sous-Forêts) which are eated by the proposed approach corresponds to the major related to major local thermal anomalies currently identified thermal anomaly known within the URG. A more detailed in the western part of the rift. The overlapping is more diffuse investigation of the delineated favourable area close to for the third area where the in situ data is more scattered as Soultz-sous-Forêts enables the highlighting of the key situa- well. Given that no parameter fitting was done, the delineated tions (structurally, geomechanically, and hydraulically) that favourable areas and local known thermal anomalies tend to make it a special location for geothermal resources. Struc- be in reasonable agreement, which underlines the potential of turally, the west-dipping fault resulting in the intersection such a cross-numerical approach. These promising results at depth with east-dipping border faults forms a key struc- suggest that the numerical approach (discrete approach to tural feature. This V-shaped structure favours loops of deep explicitly consider the effect of the fault network), the physics fluid circulation. Hydraulically, this area corresponds to a Geofluids 19 transport equations can be used to assess the favour- preferential discharge area of deep regional fluid circulation. In addition, the mechanical stress in this region favours fluid able geothermal areas circulation and then hot fluid upwelling in the rock matrix. (iii) flow interactions between fractures and porous In conclusion, the Soultz-sous-Forêts area is probably the media by allowing the superimposition of 1D ele- major thermal anomaly because this area merges structur- ments generated on disk-shaped fracture planes (to ally, geomechanically, and hydraulically favourable condi- simulate fractures) on 3D elements (to discretize tions. This evidence is consistent with previous studies and the rock matrix) demonstrates the primary role of the regional fault network. Then, the constraints provided by the respective models (iv) transport simulated by the particle-tracking method demonstrate what we can learn from these first models The flow equation solved in the 1D classical elements is and how we can improve the proposed approach. This study provides suggestions for future research on the dynamic context of the geothermal resource setup. The ∂ h ∂h C = A S , A 1 c S S very encouraging results underline the potential of such a i i ∂x ∂t cross-numerical approach to help locate geothermal 3 −1 resources. This type of approach could ensure the localiza- where C (m s ) is the pipe conductivity; h (m) the hydrau- tion of four-fifths of all thermal anomalies. Indeed, the lic head; x (m) the abscissa along an element; A (m ) the sec- −1 hydrothermal convection along faults (activated due to the tion of the element i; S (m ) the specific storage coefficient tectonic context) may explain 75-85% of all temperature of the element i; and t (s) the time. A pipe is formed by two anomalies [21], as demonstrated in the URG context. This infinite parallel ribbons, open enough so that a flow regime method proposes practical help in locating and estimating like the flow between two parallel plates can take place. To the extension of hidden thermal anomalies, with the pros- satisfy the laminar flow regime which dominates in rock pects of reducing (1) the required number of exploration masses, the pipe conductivity is correlated to a cubic law. wells through the optimization of the exploration area and The pipe section is then proportional to the cubic root of (2) the cost of each exploration well by minimizing the dril- its conductivity. ling depth needed to reach the target temperatures. More- For a flow plane with a regularly spaced grid of channels, over, the opportunity for a better match between resources 2 the transmissivity of the fracture T (m /s) is related to the and needs offered by the 3D regional models will facilitate pipe conductivity and d (m), the spaced grid size: the deployment of deep geothermal energy in the territory and its effective integration in local and regional distribution 3 −1 C m ⋅ s networks while minimizing the environmental impacts. To Tm /s = A 2 dm conclude, the first 3D numerical flow (using the discrete- continuum method) and mechanical models (using the dis- The diffusivity equation is solved in 3D elements: tinct element method) at the URG scale offer new research opportunities. The ultimate goal of our team’s work is to pro- ∂h vide innovative exploration technologies/methodologies to div K∇h = S + q, A 3 limit the area to explore and reduce the risks associated with ∂t geothermal exploration. −1 where K (m·s ) is the permeability of the porous media; h −1 (m) the hydraulic head; and q a source term (s ). Appendix Solute transport is solved by using the random walk method. The flow is one-dimensional everywhere in pipes A. 3FLO Numerical Code and tubes, except at intersections. The dispersion is therefore only longitudinal and is “completed” by the full mixing The 3FLO code is based on the finite element method. 3FLO occurring at intersections. A pipe porosity n (-) can be spec- can solve different types of problems, such as [82] ified, as well as a pipe dispersivity (m). In a pipe, the particle displacement is rectilinear and uniform: (i) flow in fracture networks, represented by a 3D net- work of pipes. A network of channels (called pipes) uΔt is generated on each fracture. The connection of x = x + , A 4 the channels from one plane to another through the fracture plane intersections (called tubes) results −1 where x (m), u (m·s ), Δt (s), and n are the particle position in a 3D network of 1D elements along the pipe, the fluid speed, the time step, and the pipe (ii) flow (saturated or unsaturated) in porous media by porosity, respectively. The 3FLO code enables the simulation using mixed-hybrid 3D finite elements. These ele- of the particle displacement in channels and in 3D elements. Indeed, in 3D elements, once the velocity is known, the ments are more precise than classical (Galerkin) 3D elements and allow a proper computation of face movement of a particle can be separated into a longitudinal fluxes. This helps minimize solute transport compu- displacement along a flow line (corresponding to advection tation errors, as will be detailed later on. Solute and longitudinal dispersion) and two orthogonal transverse 20 Geofluids of the assembly of blocks and joints. To be solved, the state displacements (simulating transverse dispersion). Moreover, direct exchanges are possible between 1D and 3D elements; equation must be completed with the constitutive equations, particles can originate from a pipe and enter a 3D element the boundary conditions, and the initial state. and vice versa. B.2. Constitutive Equations. The mechanical constitutive equations, or laws of behaviour, give the relation between dis- B. 3DEC Numerical Code placements and stresses. For the DEM, the constitutive equa- B.1. Overview of 3DEC Software. The mechanical software tions must be provided for both the blocks and the joints. In this study, because we work in a highly segmented medium (3DEC™ for the three-dimensional Distinct Element Code, see Itasca 2016) relies on the distinct element method, where where the fault network is expected to have the greatest the contact forces between solids are explicitly incorporated. impact on the stress redistribution, the emphasis is set on the fault zone characterization. Initially proposed to study slopes in jointed rock masses [83], the DEM has been developed ever since, especially in recent B.2.1. Rock Matrix (Blocks). A simple behaviour is considered years. Today, its scope of application has extended to rock for the rock matrix, and blocks are assumed to be linearly mass scales [84–86] down to microscopic scales [87–89]. In elastic, homogeneous, and isotropic materials: addition to supporting discontinuous displacement and stress fields, the particularity of the DEM is to consider the σ − σ = C ε, B 2 full mechanical behaviours for the contacts. When the con- tacts depict fault zones, the law of behaviours embraces the where C (Pa) is the elastic stiffness tensor, σ (Pa) and σ (Pa) behaviour of the infills and of the hanging and foot walls 0 the current and initial stress tensors, respectively, and ε (−) (which belong to the rock mass), and the law of behaviours the strain tensor. A convenient way of formulating the equa- can be chosen accordingly. tion splits the strain tensor into volumetric ε (−) and devia- The DEM thus offers the possibility to V toric ε (−) parts: (i) describe discontinuous displacement and stress σ − σ = Kε δ +2Gε , B 3 fields 0 V d (ii) explicitly account for mechanically active disconti- where K (Pa) and G (Pa) are the bulk and shear moduli, nuities and their impact on blocks respectively. The volumetric strain is the trace of the strain tensor, and the deviatoric part is given by ε = ε − ε /3δ. (iii) use complex mechanical laws for the joints d V The primary unknown in the mechanics is the displacement The solution phase of the DEM is similar to that of any u (m), and its relation with the strain is yielded by the com- continuous numerical method: at each time step, the solver patibility condition: runs until mechanical equilibrium is achieved within the sys- tem. For the deformable blocks, the problem unknowns are grad u+ grad u ε = B 4 the displacements u (m) at the mesh grid points. u is com- puted by solving the balance in B.2.2. Fault Zones (Joints). A more sophisticated behaviour is considered for the fault zones. Fault zones are mechanical div σ dV + F + mg = mu, B 1 ext weaknesses in the rock mass and accommodate more defor- mation than the surrounding matrix. Because of the larger deformations involved, irreversible processes can more easily where Ω (m ) and m (kg) are the volume and mass for the take place within the fault zones and should be considered in considered block, σ (Pa) the stress tensor, F (N) the sum ext −2 the model. Recalling that joints are flat surfaces in 3DEC™, of external forces other than gravity forces, g (m·s ) the the joint behaviour can be described along the normal and gravity vector, and u (m) the displacement, and a top dot parallel directions to address opening/closing and shearing denotes a time derivative. phenomena, respectively. With the Coulomb slip law, the With the DEM, the external forces F account for the ext shear behaviour is characterized by a bilinear relation: an interactions with the contiguous blocks. These interaction elastic ramp and a plateau as soon as the onset of plasticity forces are obtained through the joint constitutive equations is reached: which, given a prescribed displacement, return the resulting forces. That is, in addition to the constitutive equations σ = yield k u , u < u , f u that must be given for the continuous methods, the s s s s τ = B 5 DEM also requires joint constitutive equations to solve equa- yield yield τ , u ≤ u , s s s tion B.1. The complete solution scheme at each time step is then achieved by balancing the displacements resulting from −1 where k (Pa·m ) is the joint shear stiffness, τ (Pa) the shear s s the block deformations with the forces resulting from the yield stress, τ (Pa) the plastic yield stress, u (m) the shear dis- block interactions (Figure 13). s yield yield At the whole-system scale (here, the URG model), the placement, and u = τ /k (m) the plastic yield displace- s s state equation solved in 3DEC™ is the mechanical balance ment. The plastic yield stress depends on the fault zone Geofluids 21 internal parameters and on the normal stress σ (Pa) acting Since the mechanical balance of the blocks explicitly on it: depends on F (see equation (B.1)), the impact of the joints ext on the blocks is directly expressed through this vector. Con- yield versely, the impact of the blocks on the joints is obtained τ = c + σ tan ϕ, B 6 s 0 n once the block mechanical balance is achieved: the balance modifies the displacement field u through the model, which where c (Pa) is the cohesion and ϕ ( ) the angle of internal yield will affect the forces acting on the joints through equations friction. When the onset of plasticity is reached u ≤ u , (B.5) and (B.7). yield the exceeding displacement u − u is permanent, hence the “plastic” denomination. When the joint lies in the plastic Data Availability phase, another mechanical feature is triggered: dilation. Dila- tion is the opening of the fault zones under shearing and is The data used to build the models and support the find- related to the fault zones actually being irregular surfaces ings of this study are included within the article or can rather than flat planes. The normal behaviour of the joint is be found within references cited (articles, database issued thus elastoplastic as well and reads from the GeORG project results: http://www.geopotenziale. org/home?lang=2). el dil σ = k u + u , B 7 n n n n Conflicts of Interest −1 el where k (Pa·m ) is the joint normal stiffness, u (m) the The authors declare that they have no conflicts of interest. n n dil elastic part of the normal displacement, and u (m) the irre- versible part of the normal displacement due to dilation. The Acknowledgments dilation relation introduces two additional material parame- This work was conducted in the framework of the IMAGE ters, the dilation angle ψ ( ) and the critical shear displace- project (Integrated Methods for Advanced Geothermal ment u (m) above which dilation vanishes: Exploration, grant agreement No. 608553). This work was completed under the REFLET and TEMPERER projects yield 0, u < u , (French National Research Agency grant agreement Nos. dil yield yield c ANR_10-IEED-0802-02 and ANR-10-IEED-0801-02, respec- u = B 8 u − u tan ψ, u ≤ u < u , n s s s s s tively) which are supported by the French Government through the Investments for the Future programmes. This 0, u ≤ u work also received internal funding from the French Geologi- cal Survey. Both dilation parameters are conditioned by the geome- try of the fault zone irregularities: ψ relates to the shape, References and u depends on the maximum height of the irregularities. The joint formulation remains valid for large displace- [1] B. Goldstein, G. Hiriart, R. Bertani et al., “Geothermal energy,” ments. The block equations however are restricted to in IPCC Special Report on Renewable Energy Sources and infinitesimal strains (see equation (B.4)), implying that Climate Change Mitigation, O. Edenhofer, R. Pichs-Madruga, any large strain at the model scale would be localized into Y. Sokona, K. Seyboth, P. Matschoss, S. Kadner, T. Zwickel, the joints. P. Eickemeier, G. Hansen, S. Schlomer, and C. von Stechow, Eds., Cambridge University Press, Cambridge, 2011. B.2.3. Block and Joint Interactions. For each contact joint, the [2] J.-D. VanMees, J. Hopman, G. Pall Hersir et al., IMAGE normal and shear stresses add up to the joint contact force F Beyond a Standardized Workflow for Geothermal Exploration, (N), which depicts the interaction force between the two blocks sharing the joint: [3] K. Young, T. Reber, and K. Witherbee, “Hydrothermal explo- ration best practices and geothermal knowledge exchange on F = σ nA + τ sA , B 9 c n c s c OpenEI,” in Proc. 37th Work. Geotherm. Reserv. Eng, pp. 1455–1469, Stanford, California, 2012. where A (m ) is the contact area and n (−) and s (−) the unit c [4] L. Rybach, “The geothermal conditions in the Rhine graben - a normal vectors in the normal and shear directions, respec- summary,” Bulletin Für Angewandte Geologie, vol. 12, pp. 29– tively, along the joint. For each block, the external force vec- 32, 2007. tor F is expressed as the sum of the contact forces due to [5] D. L. Siler, J. E. Faulds, B. Mayhew, and D. D. McNamara, ext the n joints surrounding the block: “Analysis of the favorability for geothermal fluid flow in 3D: Astor Pass geothermal prospect, Great Basin, northwestern nc Nevada, USA,” Geothermics, vol. 60, pp. 1–12, 2016. F = 〠 F , B 10 ext c [6] M. W. Swyer, T. T. Cladouhos, C. Forson, J. L. Czajkowski, i=1 N. C. Davatzes, and G. M. Schmalzle, “Permeability potential modeling of geothermal prospects combining regional crustal where F (N) is the contact force of the i-th joint surrounding strain rates with geomechanical simulation of fault slip and the block. volcanic center deformation: a case study for Washington 22 Geofluids [21] P. Baillieux, E. Schill, J.-B. Edel, and G. Mauri, “Localiza- State geothermal play fairways,” in 50th US Rock Mechanics/- Geomechanics Symposium 2016, Houston, Texas, 2016. tion of temperature anomalies in the Upper Rhine Graben: insights from geophysics and neotectonic activity,” Interna- [7] J.-B. Edel, K. Schulmann, and Y. Rotstein, “The Variscan tional Geology Review, vol. 55, no. 14, pp. 1744–1762, tectonic inheritance of the Upper Rhine Graben: evidence of reactivations in the Lias, Late Eocene–Oligocene up to the recent,” International Journal of Earth Sciences, vol. 96, no. 2, [22] H. Auradou, G. Drazer, A. Boschan, J. P. Hulin, and J. Koplik, pp. 305–325, 2007. “Flow channeling in a single fracture induced by shear displacement,” Geothermics, vol. 35, no. 5-6, pp. 576–588, [8] M. E. Schumacher, “Upper Rhine Graben: role of preexisting structures during rift evolution,” Tectonics, vol. 21, no. 1, pp. 6-1–617, 2002. [23] S. Gentier, E. Lamontagne, G. Archambault, and J. Riss, “Anisotropy of flow in a fracture undergoing shear and [9] P. Baillieux, E. Schill, Y. Abdelfettah, and C. Dezayes, “Possible its relationship to the direction of shearing and injection natural fluid pathways from gravity pseudo-tomography in the pressure,” International Journal of Rock Mechanics and geothermal fields of Northern Alsace (Upper Rhine Graben),” Mining Sciences, vol. 34, no. 3-4, pp. 94.e1–94.e12, Geothermal Energy, vol. 2, no. 1, p. 16, 2014. [10] C. Clauser and H. Villinger, “Analysis of conductive and [24] J. V. Rowland and R. H. Sibson, “Structural controls on convective heat transfer in a sedimentary basin, demonstrated hydrothermal flow in a segmented rift system, Taupo for the Rhein graben,” Geophysical Journal International, Volcanic Zone, New Zealand,” Geofluids, vol. 4, no. 4, 283 vol. 100, no. 3, pp. 393–414, 1990. pages, 2004. [11] S. Gentier, X. Rachez, T. Dung et al., 3D Flow Modelling [25] D. L. Siler, N. H. Hinz, J. E. Faulds, and J. Queen, “3D analysis of the Medium-Term Circulation Test Performed in the of geothermal fluid flow favorability: Brady’s, Nevada, USA,” Deep Geothermal Site of, World Geothermal Congress, in The 41st Workshop on Geothermal Reservoir Engineering, Stanford University, Standford, California, 2016. [12] L. Guillou-Frottier, C. Carrė, B. Bourgine, V. Bouchot, and [26] D. E. Dempsey, S. F. Simmons, R. A. Archer, and J. V. A. Genter, “Structure of hydrothermal convection in the Rowland, “Delineation of catchment zones of geothermal Upper Rhine Graben as inferred from corrected tempera- systems in large-scale rifted settings,” Journal of Geophysical ture data and basin-scale numerical models,” Journal of Research: Solid Earth, vol. 117, no. B10, 2012. Volcanology and Geothermal Research, vol. 256, pp. 29–49, [27] Y. Rotstein, J.-B. Edel, G. Gabriel, D. Boulanger, M. Schaming, and M. Munschy, “Insight into the structure of the Upper [13] C. Le Carlier, J.-J. Royer, and E. L. Flores, “Convetive Rhine Graben and its basement from a new compilation of heat transfer at the Soultz-sous-Forets Geothermal Site: Bouguer Gravity,” Tectonophysics, vol. 425, no. 1-4, pp. 55– implications for oil potential,” First Break, vol. 12, 70, 2006. no. 1285, 1994. [28] S. Gentier, X. Rachez, M. Peter-Borie, and M. Loubaud, “A [14] B. Sanjuan, R. Millot, C. Innocent, C. Dezayes, J. Scheiber, and flow model of the deep geothermal reservoir of Soultz-sous- M. Brach, “Major geochemical characteristics of geothermal Forêts (France),” in 47 US Rock Mechanics/Geomechanics brines from the Upper Rhine Graben granitic basement with Symposium, San Francisco, 2013. constraints on temperature and circulation,” Chemical Geol- ogy, vol. 428, pp. 27–47, 2016. [29] B. C. Valley, “The relation between natural fracturing and stress heterogeneities in deep-seated crystalline rocks at [15] I. Spottke, E. Zechner, and P. Huggenberger, “The southeast- Soultz-sous-Forêts (France),” in ETH Reprozentrale Höngger- ern border of the Upper Rhine Graben: a 3D geological model berg, HIL C45, Zürich (2007), Swiss Federal Institute of Tech- and its importance for tectonics and groundwater flow,” Inter- national nology Zurich, 2007. Journal of Earth Sciences, vol. 94, no. 4, pp. 580–593, [30] T. Plenefisch and K.-P. Bonjer, “The stress field in the Rhine Graben area inferred from earthquake focal mechanisms and [16] D. King and E. Metcalfe, “Rift zones as a case study for advancing geothermal occurence models,” PROCEEDINGS: estimation of frictional parameters,” Tectonophysics, vol. 275, no. 1-3, pp. 71–97, 1997. Thirty-Eighth Workshop on Geothermal Reservoir Engineer- ing, Standford University, p. 11, 2013. [31] J. B. Witter and G. Melosh, “The value and limitations of 3D models for geothermal exploration,” in 43rd [17] J. Freymark, J. Sippel, M. Scheck-wenderoth et al., The Thermal Field of the Upper Rhine Graben–Temperature Pre- Workshop on Geothermal Reservoir Engineering, no. article SGP-TR-213, 2018Standford Univeristy, Standford, California, dictions Based on a 3D Model, Eur. Geotherm. Congr. 2016, 2016. 2018. [18] P. Baillieux, Multidisciplinary Approach to Understand the [32] D. D. Pollard and P. Segall, “Theoretical displacements and stresses near fractures in rock: with applications to faults, Localization of Geothermal Anomalies in the Upper Rhine Gra- ben from Regional to Local Scale, Fac. Sci. Hydrogeol. Geotherm, joints, veins, dikes, and solution surfaces,” in Fracture Mechanics of Rock, pp. 277–349, Academic Press, San Diego, [19] J. Tóth, “Geothermal phenomena in the context of gravity- driven basinal flow of groundwater,” Central European Geol- [33] C.-F. Tsang and I. Neretnieks, “Flow channeling in heteroge- ogy, vol. 58, no. 1-2, pp. 1–27, 2015. neous fractured rocks,” Reviews of Geophysics, vol. 36, no. 2, pp. 275–298, 1998. [20] J. Vidal and A. Genter, “Overview of naturally permeable fractured reservoirs in the central and southern Upper Rhine [34] A. Herbert, Modelling Approaches for Discrete Fracture Graben: insights from geothermal wells,” Geothermics, Network Flow Analysis, Coupled Thermo-Hydro-Mechanical vol. 74, pp. 57–73, 2018. Processes of Fractured Media, Elsevier, Amsterdam, 1996. Geofluids 23 [52] Y. Zhou and W. Li, “A review of regional groundwater flow [35] Itasca Consulting Group Inc, 3DEC–Three-Dimensional Distinct Element Code, Ver. 5.2, Itasca Consulting Group modeling,” Geoscience Frontiers, vol. 2, no. 2, pp. 205–214, Inc., Minneapolis, Minnesota, 2017. 2011. [36] GeORG project, 2012, [WWW Document]. [53] F. Rummel, “Physical properties of the rock in the granitic sec- tion of borehole GPK1 borehole, Soultz-sous-Forêts,” in Geo- [37] G. B. Baecher, “Statistical analysis of rock mass fracturing,” thermal Energy in Europe: The Soultz Hot Dry Rock Project, Journal of the International Association for Mathematical pp. 199–216, Gordon & Breach Science Publishers Ltd, 1992. Geology, vol. 15, no. 2, pp. 329–348, 1983. [54] I. Stober and K. Bucher, “Hydraulic and hydrochemical [38] P. H. S. W. Kulatilake and T. H. Wu, “Estimation of mean trace properties of deep sedimentary reservoirs of the Upper Rhine length of discontinuities,” Rock Mechanics and Rock Engineer- Graben, Europe,” Geofluids, vol. 15, no. 3, 482 pages, 2015. ing, vol. 17, no. 4, pp. 215–232, 1984. [55] T. Gleeson, L. Marklund, L. Smith, and A. H. Manning, [39] S. D. Priest and J. A. Hudson, “Estimation of discontinuity “Classifying the water table at regional to continental scales,” spacing and trace length using scanline surveys,” International Geophysical Research Letters, vol. 38, no. 5, 2011. Journal of Rock Mechanics and Mining Science and Geomecha- nics Abstracts, vol. 18, no. 3, pp. 183–197, 1981. [56] H. M. Haitjema and S. Mitchell-Bruker, “Are water tables a subdued replica of the topography?,” Ground Water, vol. 43, [40] L. Zhang and H. H. Einstein, “Estimating the intensity of rock no. 6, pp. 781–786, 2005. discontinuities,” International Journal of Rock Mechanics and Mining Sciences, vol. 37, no. 5, pp. 819–837, 2000. [57] D. P. Yale, “Fault and stress magnitude controls on variations in the orientation of in situ stress,” Geological Society, London, [41] W. Riedel, “Zur Mechanik Geologischer Brucherscheinun- Special Publications, vol. 209, no. 1, pp. 55–64, 2003. gen,” Zur Mechanik Geologischer Brucherscheinungen B, pp. 354–368, 1929. [58] L. Dorbath, K. Evans, N. Cuenot, B. Valley, J. Charléty, and M. Frogneux, “The stress field at Soultz-sous-Forêts from focal [42] J. H. Behrmann, O. Hermann, M. Horstmann, D. C. Tanner, mechanisms of induced seismic events: cases of the wells and G. Bertrand, “Anatomy and kinematics of oblique conti- GPK2 and GPK3,” Comptes Rendus Geoscience, vol. 342, nental rifting revealed: a three-dimensional case study of the no. 7-8, pp. 600–606, 2010. southeast Upper Rhine Graben (Germany),” American Associ- ation of Petroleum Geologists Bulletin, vol. 87, no. 7, pp. 1105– [59] T. Hettkamp, J. Baumgärtner, R. Baria et al., Electricity Produc- 1121, 2003. tion from Hot Rocks, in: Proceedings, Twenty-Ninth Workshop on Geothermal Reservoir Engineering, Stanford University, [43] N. Harthill, No Title. Tecton. Basis Geotherm. Reg. Oberrhein- Stanford, California, 2018. graben, 20 Jahre Tiefe Geotherm. Deutschland. 7. Geotherm. Fachtagung, Waren, 2002. [60] Y. Guglielmi, D. Elsworth, F. Cappa et al., “In situ observations on the coupling between hydraulic diffusivity and displace- [44] P. A. Domenico and V. V. Palciauskas, “Theoretical analysis of ments during fault reactivation in shales,” Journal of Geophys- forced convective heat transfer in regional ground-water flow,” ical Research: Solid Earth, vol. 120, no. 11, pp. 7729–7748, Geological Society of America Bulletin, vol. 84, no. 12, p. 3803, [61] X. Rachez and S. Gentier, 3D-Hydromechanical Behavior of a [45] M. P. Anderson, “Heat as a ground water tracer,” Ground Stimulated Fractured Rock Mass, vol. 8, World Geotherm. Water, vol. 43, no. 6, pp. 951–968, 2005. Congr. 2010, 2010. [46] M. O. Saar, “Review: geothermal heat as a tracer of large-scale [62] T. J. Buchmann and P. T. Connolly, “Contemporary kinemat- groundwater flow and as a means to determine permeability ics of the Upper Rhine Graben: a 3D finite element approach,” fields,” Hydrogeology Journal, vol. 19, no. 1, pp. 31–52, 2011. Global and Planetary Change, vol. 58, no. 1-4, pp. 287–309, [47] J. Faulds, M. Coolbaugh, V. Bouchot, I. Moeck, O. Kerem, and O. Cedex, Characterizing Structural Controls of Geothermal [63] A. Hosni, Modélisation par la Méthode des Éléments Dis- Reservoirs in the Great Basin, USA, and Western Turkey: tincts du champ des contraintes à l’échelle du Fossé Rhénan Developing Successful Exploration Strategies in Extended Ter- et de Soultz-sous-Forêts, Institu National Polytechnique de ranes, pp. 25–29, 2010. Lorraine, 1997. [48] J. E. Faulds, N. H. Hinz, M. F. Coolbaugh et al., “Assessment of [64] G. Peters, Active Tectonics in the Upper Rhine Graben, p. 298, favorable structural settings of geothermal systems in the Great Basin, Western USA,” Geothermal Resources Council Transactions, vol. 35, pp. 777–783, 2011. [65] Y. Gunzburger and V. Magnenet, “Stress inversion and basement-cover stress transmission across weak layers in the [49] N. H. Hinz, J. E. Faulds, and D. L. Siler, “Developing systematic Paris basin, France,” Tectonophysics, vol. 617, pp. 44–57, 2014. workflow from field work to quantitative 3D modeling for successful exploration of structurally controlled geothermal [66] T. Hergert, O. Heidbach, K. Reiter, S. B. Giger, and systems,” Geothermal Resources Council Transactions, P. Marschall, “Stress field sensitivity analysis in a sedimentary vol. 37, pp. 275–279, 2013. sequence of the Alpine foreland, northern Switzerland,” Solid Earth, vol. 6, no. 2, pp. 533–552, 2015. [50] I. C. Wallis, D. McNamara, J. V. Rowland, and C. Massiot, “The nature of fracture permeability in the basement grey- [67] K. Reiter and O. Heidbach, “3-D geomechanical-numerical wacke at Kawerau Geothermal Field, New Zealand,” in Pro- model of the contemporary crustal stress state in the Alberta ceedings Thirty-seventh Workshop on Geothermal Reservoir Basin (Canada),” Solid Earth, vol. 5, no. 2, pp. 1123–1149, Engineering, Stanford University, California, 2012. 2014. [51] R. A. Freeze and P. A. Witherspoon, “Theoretical analysis of [68] P. Connolly and J. Cosgrove, “Prediction of static and dynamic regional groundwater flow: 1. Analytical and numerical solu- fluid pathways within and around dilational jogs,” Geological tions to the mathematical model,” Water Resources Research, Society, London, Special Publications, vol. 155, no. 1, pp. 105– vol. 2, no. 4, pp. 641–656, 1966. 121, 1999. 24 Geofluids [84] V. Bonilla-Sierra, L. Scholtès, F. V. Donzé, and M. K. Elmouttie, [69] D. Curewitz and J. A. Karson, “Structural settings of hydro- thermal outflow: fracture permeability maintained by fault “Rock slope stability analysis using photogrammetric data and propagation and interaction,” Journal of Volcanology and Geo- DFN–DEM modelling,” Acta Geotechnica, vol. 10, no. 4, thermal Research, vol. 79, no. 3-4, pp. 149–168, 1997. pp. 497–511, 2015. [70] M. Cathelineau and M.-C. Boiron, “Downward penetration [85] K. Farahmand, I. Vazaios, M. S. Diederichs, and N. Vlachopoulos, “Investigating the scale-dependency of the and mixing of sedimentary brines and dilute hot waters at 5km depth in the granite basement at Soultz-sous-Forêts geometrical and mechanical properties of a moderately jointed (Rhine graben, France),” Comptes Rendus Geoscience, vol. 342, rock using a synthetic rock mass (SRM) approach,” Computers no. 7-8, pp. 560–565, 2010. and Geotechnics, vol. 95, pp. 162–179, 2018. [71] M. Cacace and A. B. Jacquey, “Flexible parallel implicit model- [86] J. H. Wu, W. K. Lin, and H. T. Hu, “Assessing the impacts of a ling of coupled thermal-hydraulic-mechanical processes in large slope failure using 3DEC: the Chiu-fen-erh-shan residual fractured rocks,” Solid Earth, vol. 8, no. 5, pp. 921–941, 2017. slope,” Computers and Geotechnics, vol. 88, pp. 32–45, 2017. [72] Z. Q. Huang, P. H. Winterfeld, Y. Xiong, Y. S. Wu, and [87] E. Ghazvinian, M. S. Diederichs, and R. Quey, “3D random J. Yao, “Parallel simulation of fully-coupled thermal-hydro- Voronoi grain-based models for simulation of brittle rock mechanical processes in CO2 leakage through fluid-driven damage and fabric-guided micro-fracturing,” Journal of fracture zones,” International Journal of Greenhouse Gas Rock Mechanics and Geotechnical Engineering, vol. 6, no. 6, Control, vol. 34, pp. 39–51, 2015. pp. 506–521, 2014. [73] A. Lisjak, O. K. Mahabadi, L. He et al., “Acceleration of a [88] N. P. Kruyt and L. Rothenburg, “A micromechanical study of 2D/3D finite-discrete element code for geomechanical simula- dilatancy of granular materials,” Journal of the Mechanics tions using general purpose GPU computing,” Computers and and Physics of Solids, vol. 95, pp. 411–427, 2016. Geotechnics, vol. 100, pp. 84–96, 2018. [89] W. Li, Y. Han, T. Wang, and J. Ma, “DEM micromechanical [74] S. Salimzadeh, A. Paluszny, H. M. Nick, and R. W. Zimmerman, modeling and laboratory experiment on creep behavior of salt “A three-dimensional coupled thermo-hydro-mechanical rock,” Journal of Natural Gas Science and Engineering, vol. 46, model for deformable fractured geothermal systems,” Geother- pp. 38–46, 2017. mics, vol. 71, pp. 212–224, 2018. [90] T. G. Farr and M. Kobrick, “Shuttle Radar Topography Mis- [75] B. Vallier, V. Magnenet, J. Schmittbuhl, and C. Fond, “THM sion produces a wealth of data,” Eos, Transactions American modeling of hydrothermal circulation at Rittershoffen geo- Geophysical Union, vol. 81, no. 48, pp. 583–583, 2000. thermal site, France,” Geothermal Energy, vol. 6, no. 1, 2018. [76] M. O. Ziegler, O. Heidbach, J. Reinecker, A. M. Przybycin, and M. Scheck-Wenderoth, “A multi-stage 3-D stress field model- ling approach exemplified in the Bavarian Molasse Basin,” Solid Earth, vol. 7, no. 5, pp. 1365–1382, 2016. [77] J. P. Brun, M.-A. Gutscher, and dekorp-ecors teams, “Deep crustal structure of the Rhine Graben from dekorp-ecors seis- mic reflection data: a summary,” Tectonophysics, vol. 208, no. 1-3, pp. 139–147, 1992. [78] M. Schwarz and A. Henk, “Evolution and structure of the Upper Rhine Graben: insights from three-dimensional thermomechanical modelling,” International Journal of Earth Sciences, vol. 94, no. 4, pp. 732–750, 2005. [79] F. Wenzel, J. P. Brun, and ECORS‐DEKORP team, “Crustal scale structure of the Rhine Graben from ECORS‐DEKORP deep seismic reflection data,” in SEG Technical Program Expanded Abstracts 1991, pp. 163–166, Houston, Texas, 1991. [80] J. R. Harris, E. Grunsky, P. Behnia, and D. Corrigan, “Data- and knowledge-driven mineral prospectivity maps for Canada’s North,” Ore Geology Reviews, vol. 71, pp. 788–803, 2015. [81] J. Iovenitti, F. H. Ibser, M. Clyne, J. Sainsbury, and O. Callahan, “The Basin and Range Dixie Valley Geothermal Wellfield, Nevada, USA—a test bed for developing an Enhanced Geothermal System exploration favorability meth- odology,” Geothermics, vol. 63, pp. 195–209, 2016. [82] B. Paris, Investigation of the Correlation between Early-Time Hydraulic Response and Tracer Breakthrough Times in Fractured Media, Äspö Hard Rock Laboratory, 2002. [83] P. A. Cundall, “A computer model for simulating progres- sive large scale movements in blocky rock systems,” in Proc. Symp. Int. Soc. rock Mech. (Nancy, Fr. 1971) 1, Nancy, France, 1971. International Journal of The Scientific Advances in Advances in Geophysics Chemistry Scientica World Journal Public Health Hindawi Hindawi Hindawi Hindawi Publishing Corporation Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 http://www www.hindawi.com .hindawi.com V Volume 2018 olume 2013 www.hindawi.com Volume 2018 Journal of Environmental and Public Health Advances in Meteorology Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Submit your manuscripts at www.hindawi.com Applied & Environmental Journal of Soil Science Geological Research Hindawi Volume 2018 Hindawi www.hindawi.com www.hindawi.com Volume 2018 International Journal of International Journal of Agronomy Ecology International Journal of Advances in International Journal of Forestry Research Microbiology Agriculture Hindawi Hindawi Hindawi Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 International Journal of Journal of Journal of International Journal of Biodiversity Archaea Analytical Chemistry Chemistry Marine Biology Hindawi Hindawi Hindawi Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geofluids Hindawi Publishing Corporation

Locating Geothermal Resources: Insights from 3D Stress and Flow Models at the Upper Rhine Graben Scale

Loading next page...
 
/lp/hindawi-publishing-corporation/locating-geothermal-resources-insights-from-3d-stress-and-flow-models-QYvrs1aaRK

References (95)

Publisher
Hindawi Publishing Corporation
Copyright
Copyright © 2019 Antoine Armandine Les Landes et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
ISSN
1468-8115
eISSN
1468-8123
DOI
10.1155/2019/8494539
Publisher site
See Article on Publisher Site

Abstract

Hindawi Geofluids Volume 2019, Article ID 8494539, 24 pages https://doi.org/10.1155/2019/8494539 Research Article Locating Geothermal Resources: Insights from 3D Stress and Flow Models at the Upper Rhine Graben Scale Antoine Armandine Les Landes , Théophile Guillon, Mariane Peter-Borie, Arnold Blaisonneau, Xavier Rachez, and Sylvie Gentier BRGM, Orléans, France Correspondence should be addressed to Antoine Armandine Les Landes; a.armandineleslandes@brgm.fr Received 20 December 2018; Revised 1 March 2019; Accepted 12 March 2019; Published 12 May 2019 Guest Editor: Victor Vilarrasa Copyright © 2019 Antoine Armandine Les Landes et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. To be exploited, geothermal resources require heat, fluid, and permeability. These favourable geothermal conditions are strongly linked to the specific geodynamic context and the main physical transport processes, notably stresses and fluid circulations, which impact heat-driving processes. The physical conditions favouring the setup of geothermal resources can be searched for in predictive models, thus giving estimates on the so-called “favourable areas.” Numerical models could allow an integrated evaluation of the physical processes with adapted time and space scales and considering 3D effects. Supported by geological, geophysical, and geochemical exploration methods, they constitute a useful tool to shed light on the dynamic context of the geothermal resource setup and may provide answers to the challenging task of geothermal exploration. The Upper Rhine Graben (URG) is a data-rich geothermal system where deep fluid circulations occurring in the regional fault network are the probable origin of local thermal anomalies. Here, we present a current overview of our team’sefforts to integrate the impacts of the key physics as well as key factors controlling the geothermal anomalies in a fault-controlled geological setting in 3D physically consistent models at the regional scale. The study relies on the building of the first 3D numerical flow (using the discrete-continuum method) and mechanical models (using the distinct element method) at the URG scale. First, the key role of the regional fault network is taken into account using a discrete numerical approach. The geometry building is focused on the conceptualization of the 3D fault zone network based on structural interpretation and generic geological concepts and is consistent with the geological knowledge. This DFN (discrete fracture network) model is declined in two separate models (3D flow and stress) at the URG scale. Then, based on the main characteristics of the geothermal anomalies and the link with the physics considered, criteria are identified that enable the elaboration of indicators to use the results of the simulation and identify geothermally favourable areas. Then, considering the strong link between the stress, fluid flow, and geothermal resources, a cross-analysis of the results is realized to delineate favourable areas for geothermal resources. The results are compared with the existing thermal data at the URG scale and compared with knowledge gained through numerous studies. The good agreement between the delineated favourable areas and the locations of local thermal anomalies (especially the main one close to Soultz-sous-Forêts) demonstrates the key role of the regional fault network as well as stress and fluid flow on the setup of geothermal resources. Moreover, the very encouraging results underline the potential of the first 3D flow and 3D stress models at the URG scale to locate geothermal resources and offer new research opportunities. 1. Introduction it is environmentally friendly, has low emissions (greenhouse gases), and uses a low number of raw materials such as rare Worldwide, deep geothermal energy offers an enormous elements and strategic metals. In addition, unlike many other potential for future electricity and heat productions. Geother- renewable energies depending on climate conditions (wind, mal energy is sustainable, clean, and renewable, as the tapped solar, hydraulic regimes, and the water levels of rivers), the heat from an active reservoir is continuously restored [1]. geothermal energy from the Earth’s crust is nonintermittent. Geothermal energy has a number of positive characteristics: However, only a fraction of the deep geothermal energy of the 2 Geofluids the widespread deployment of deep geothermal energy to Earth’s crust is currently exploited. In practice, the exploita- tion is mostly limited to “brown fields,” i.e., areas where a res- other areas, for instance, the Limagne Graben. ervoir is already sited or where there are targeted structures Several researchers have suggested that the stresses asso- and parameters characterized that are already known [2]. ciated with faults significantly affect the flow properties and The wider development of deep geothermal resource exploi- that there is a strong correlation between the groundwater tation now requires focusing on little to unknown areas, flow and geothermal anomalies. A regional heat-flow study called “green fields.” In these fields, geothermal resource requires an analysis of the basin-scale flow distribution (due exploration (including exploratory drillings) involves a high to the strong link between the regional groundwater flow degree of uncertainty and financial risk and requires reliable and anomalies of the geothermal heat) and should account data to provide convincing arguments in the decision- for the active faulting caused by stresses. Stresses can play making process [3]. Consequently, the improvement of an important role in the setup of thermal anomalies by exploration methods is required to lower the risk and to over- favouring preferential flow paths through the rock mass [6]. come the barriers to the wider development of deep geother- For example, shear stresses on fractures have been observed mal projects. to favour the fluid circulation perpendicular to the shear To be exploited, geothermal resources require heat, fluid, direction [22, 23]. Then, stresses can favour the appearance and permeability. These favourable geothermal conditions of thermal anomalies by enabling vertical upflows of a fluid are strongly linked to the specific geodynamic context and previously heated during deep circulation [6, 24]. Thus, key the main physical transport processes, notably the stresses factors to understanding thermal anomalies are faults and permeability heterogeneities [25], which are influenced by and fluid circulations that impact heat-driving processes [4–6]. The Upper Rhine Graben (URG) is a relatively well- the stresses (a key physic) that control both the fluid circula- known example of a favourable faulted geothermal system tion (a key physic) [6, 25] and the location of exploitable hot showing thermal anomalies at the regional scale. It has been fluid upwelling. That is the reason why geothermal studies of largely investigated, providing various data such as a tectonic an extended area, such as exploration for geothermal history of the Rhine Graben [7, 8], in situ observations from resources, must consider these key factors and must be based deep boreholes, and numerical models and geochemical on the study of these physics. studies [9–15]. Moreover, a European project (the GeORG In addition to other characterization methods (geologi- project) has resulted in an important structural database cal, geophysical, and geochemical), exploration could benefit (http://www.geopotenziale.org/home?lang=3) that provides from numerical models and then a predictive analysis to locate geothermal resources. Numerical models could allow a map of seismic faults at the URG scale. The formation of geothermal resources in the URG results from the interac- an integrated evaluation of the physical processes in the tions between different coupled processes comprising URG with adapted time and space scales, considering 3D mechanical constraints, groundwater flows, and heat trans- effects. They constitute a useful tool to improve the under- fer. Large-scale structures of the graben and its associated standing of the dynamic context and may provide answers to the challenging task of geothermal exploration. Indeed, faults may act as a primary structural control on the upwell- ing geothermal fluids [16]. The rifting context and the asso- the physical conditions favouring the setup of the geother- ciated thinned crust and elevated mantle favour higher mal resources can be searched for in predictive models, temperatures at shallower depths [4, 17]. A series of thermal hence giving estimates of the so-called “favourable areas.” anomalies with temperatures above 140 C at a 2 km depth Some major phenomena were highlighted, and some local anomalies were explained using 1D hydrothermal models characterizes the URG. These anomalies are mainly located in the western part of the URG [18]. Recalling that thermal and 2D models based on continuum approaches [10, 12, anomalies are mostly affected by the regional groundwater 13] but without any hope for an industrial application to circulation [19], the high potential of the URG for thermal exploration. These models are indeed restricted to local anomalies is widely interpreted as a signature of the scales and do not consider major features such as preferen- tial pathways within fault zones or 3D effects. Stress and regional-scale fluid migrations hosted by the numerous fault zones cutting through the rock mass (multiscale fracture- numerical flow models have also been used to study geo- controlled systems, [20]). Arguing that the fault zones are thermal systems ([6, 26] and references therein), but either critically stressed and mostly under shearing conditions, the scale is not appropriate or the approach is not holistic. Baillieux et al. [21] estimated that the regional fluid migra- Thus, to the authors’ knowledge, there is currently no numerical approach available to contribute to geothermal tions could be responsible for up to 75-85% of the thermal anomalies in the URG. The radiogenic heat production due resource exploration. to the crystalline composition of the basement may explain The present work aims at demonstrating the potential of the remaining 15-25%. The geochemical data tend to confirm the combined use of 3D flow and stress models at the this assumption: Sanjuan et al. [14] conclude that the thermal regional scale (URG scale) to locate/understand the location of hidden thermal anomalies. Considering the large amount anomaly in the Soultz-sous-Forêts surroundings results from the regional circulation through a complex, but still poorly of data and knowledge, we use the URG as an application defined, system of deep faults (probably along NE-SW fault case. The results of the GeORG project constitute a relevant zones). Thus, the deep circulation is probably the origin of basis to build a structural model, which is essential to devel- the thermal anomalies and governs their geographical distri- oping the subsequent numerical physical models and then such a numerical approach. This enables us to compare butions. Their understanding is therefore fundamental for Geofluids 3 results with observed physical evidence. The main objective (b) N170 sinistral strike-slip faults reactivated from of this paper is to investigate the potential of an approach the late Miocene up to now under a NW-SE based on the use of numerical models (stress and flow compression of the graben as a strike-slip system models) to locate hidden thermal anomalies. The paper is (ii) The F2 set, ~NNW-SSE striking and E-dipping, is organized as follows: first, we give a brief overview of the also commonly interpreted as the result of the Oligo- study area. Then, we present the 3D geometry for the models cene deformation (N0/20 striking, E-dipping). according to the faulted geological setting. The major faults However, it could also be considered as different sets of the URG are considered using a DFN model (detailed in gathering the Hercynian N20 sinistral strike-slip Section 3) consistent with the geological knowledge of the and the N135 dextral strike-slip reactivated during area (see Section 2). This DFN (discrete fracture network) the late Eocene and to N170 sinistral strike-slip model is declined in 3D flow and stress models at the URG faults scale, which are, respectively, detailed in Sections 4 and 5. The results of the physical models are interpreted as func- (iii) The F3 and F4 sets, respectively, NE-SW and tions of their thermal signatures through the definition of NW-SE striking, are consistent with the Hercy- geothermal criteria related to the physics considered and nian orientations observed regionally the setting up of the geothermal resources. These criteria (iv) A fifth fault set (F5) is introduced to incorporate the enable the elaboration of indicators to use the results of the “background noise” simulation and then identify geothermally favourable areas in relation with each physic. Then, considering the strong Considering the stress regime, the sinistral strike-slip link between the stress, fluid flow, and geothermal resources, regime that started at the onset of the Miocene has been a cross-analysis of the results is realized to delineate the persisting until today, and the present-day stress state is favourable areas for geothermal resources. The results are assumed to be driven by the African-Eurasian plate conver- compared with existing thermal data at the URG scale and gence [30]. with knowledge gained through numerous studies (see Sec- As a consequence of its polyphaser tectonic history, the tion 6). Before concluding, we discuss the potential of the URG is a complex environment setting in which the rifting present work. process played a major role by providing the space for sedi- mentation, creating a series of basins, and providing a net- work of faults along the active rift margins that enabled 2. Geological Setting of the Study Area: fluid circulation and conditioned water exchanges between the URG basins [8, 14]. The transnational geological model of the URG based on boreholes and new processed seismic data The URG is an NNE-trending crustal-scale rift of the constitutes an important piece of the puzzle for the numerical European Cenozoic Rifting System (ECRIS) [8]. Its length models by providing essential data. The map of the faults at reaches approximately 300 km from Frankfurt (Germany) the top of the basement can be used to develop numerical to Basel (Switzerland), and its width is 30-40 km (Figure 1). models, and the in situ temperature data at the URG scale The URG development started from the late Eocene, mainly constitute a relevant means to assess the fault potential by the reactivation of the late Variscan, Permo-Carbonifer- (Figure 1). ous, and Mesozoic fracture systems [8]. Schumacher [8], Edel et al. [7], and Rotstein et al. [27], among others, detail the tectonic history and the resulting structure. Only some key 3. DFN Model at the URG Scale aspects are recalled here. 3.1. Why a Discrete Approach. The very first step of our work Since the Visean age during which the granitic intrusion is to build a model geometry consistent with the geology and took place, a repeatedly changing stress field has occurred the physical processes involved. The large heterogeneities in the URG, leading to the reactivation of a complex set of introduced in the physical system by the presence of the fault crustal discontinuities and the creation of a complex fault zones rule out any chance for an analytical solution, thus network in the sedimentary blanket [8]. Based on the tectonic making numerical tools necessary. As argued by Witter and history and on knowledge acquired from previous works Melosh [31], the choice of the numerical tools is a challenge [28, 29], four main fault sets corresponding to the main sta- in itself since it will impact the results. The need to account tistical models have been individualized: for real 3D effects raises another challenge, which is the avail- able panel of suitable tools: the numerical tools for solving (i) The F1 set, mainly ~N-S striking and W-dipping, is large-scale models in reasonable computation times while widely observed in the Rhine Graben and commonly offering the possibility to take into account geological struc- interpreted as the result of the Oligocene deforma- tures are few. In the case of faulted settings, the presence of tion (N0/20 striking, W-dipping). However, this natural discontinuities that can result in heterogeneous stress set could also be related to fields and localized fluid flow pathways [32, 33] has pro- (a) the lower Carboniferous to Permian N20 sinis- moted the development of fracture network models for the tral faulting reactivated during the late Eocene numerical simulation of fractured rocks [34]. Here, the key [8] or role of the regional fault network is taken into account due 4 Geofluids Worms Mannheim Mannheim United Kingdom Germany Belgium Heidelberg Speyer Heidelberg France Switzerland Landau Speyer Italy Bruchsal Wissembourg Spain Karlsruhe Landau Soultz Haguenau Bruchsal Baden-Baden Wissembourg Karlsruhe Germany Strasbourg France Soultz Obernai Offenburg Lahr Haguenau Selestat Baden-Baden Colmar Freiburg Strasbourg Mulhouse Obernai Offenburg Lurrach Switzerland Lahr km km Legend Temperature at 2 km depth Cities Value High: 167 Faults_GeORG_basement top Low: 61 Country outlines Figure 1: Location of the URG and of the case study in the URG. The map shows the faults at the top of the basement and the temperature map at 2 km depth delivering the main thermal anomalies (GeORG project). expertise on geological settings similar to the one studied (in to a discrete numerical approach (DFN) for fault zones at the regional scale that remains a challenge alone. The model to our case, rifting systems), an understanding of the tectonic represent the regional fault network is generated with the history, and the integration of generic structural tools (such 3DEC™ software [35]. as Riedel systems, see [41]). Furthermore, the numerical constraints and specificities of the codes used have to be 3.2. Geometry Building. The case study is a 130 × 150 km taken into account. area centred on the main thermal anomalies in the URG Following the above-defined concept and using the geo- (Figure 1). The lateral east-west boundaries are chosen large logical knowledge of the URG, the discontinuities of the enough to incorporate a large portion of the rift shoulders to DFN model are the geological faults observed at different limit the boundary effects in the subsequent computations (at scales, from plurikilometric to kilometric scales. The largest discontinuities that cross the model are on the regional scale least in this direction). The geometry building is focused on the conceptualization of the 3D fault zone network into a (primary faults of hundreds of kilometres known as major DFN at the studied scale; more specifically, the primary faults of the Variscan or tertiary rifting history). As disconti- source of the data is the seismic fault map of the transna- nuities split the model blocks, the resulting block sizes tional database [36], where any valuable data (e.g., seismic decrease, allowing the integration of more local-scale discon- reflection results, borehole data, outcrop traces, and scanline tinuities (faults of several kilometres from seismic profiles or window sampling) should be integrated into the structural and lineament maps). The main faults integrated in the interpretation to provide estimates on fault zone orienta- DFN model are the faults delimiting the URG and those tions, dip angles, and spatial extents [37–40]. Additional inherited from the main Variscan NE-SW trending sinistral generic geological knowledge is required to complete the shearing (Lubine and Baden-Baden faults) ([7, 8, 42], among missing information (e.g., the dips of some deep fault zones) numerous others), as shown in Figure 2(a). The secondary and to guide the conceptualization into a DFN (e.g., the fault fault network within the basement of the URG has been built zone priority rule dictating “who stops on who”). This from data of the GeORG project (within the basement and includes the following nonexhaustive list of items: structural sedimentary layers, see [36]) or from data found in the Geofluids 5 010 20 km Secondary fault network within the basement on both sides of the URG Primary fault network: Secondary fault network within the basement of the URG DFN Major fault inherited of the  GeORG fault network Variscan orogeny European digital terrain model DFN Faults delimiting the URG Lubine(LB) / Baden‐Baden fault (Variscan orogeny) Secondary fault network Fault network (a) (b) (c) Figure 2: Illustration of the DFN building method: (a) main geological structures (example of integration of structures highlighted by [8]), (b) secondary fault network from the GeORG database within the URG [36], and (c) secondary fault network in the basement on both sides of the URG from DTM analysis [90]. literature [7, 27] (Figure 2(b)). The secondary fault network The same 3D geometry will be used for the 3D flow and in the basement on both sides of the URG is from geological stress models. maps and digital terrain model (DTM) analysis (Figure 2(c)). Note that the fault network is denser in the URG than on 4. 3D Flow Model at the URG Scale the borders, as the study is focused on the URG. The dips of the faults are chosen to be equal to the fault dip at the top of Regionally moving groundwater is the basic and common the basement (estimated using the fault offset). It should be cause of a wide variety of natural processes and phenomena noted that additional geometric simplifications have been that makes the gravity-driven groundwater flow a key pro- applied to fulfil the 3DEC™ topological requirements: cess [19, 44]. The key reason for the geological scale effects of the groundwater flow is that the flow systems are ubiqui- (i) Planar surfaces are considered in 3DEC™, so the tous and active over large spectra of time and space. Among initial curved fault zones were flattened the numerous processes associated with regional groundwa- ter flow systems, the transport of heat is the most visible (ii) 3DEC™ does not accept discontinuity tips within and best understood [45, 46]. Many geothermal systems blocks, so the fault zones were extended until are associated with high densities of faults and fractures they stopped on a boundary or a previously created [25, 47–50]. Geological discontinuities such as faults and fault zone stratigraphic layers play a key role in controlling the ground- The fault zone extension requires the definition of a water flow within faulted systems [51]. priority rule to assess which fault zone stops on which one. The numerical simulation of groundwater flows has The priority rule was defined according to the tectonic become a common hydrogeological tool to investigate a vari- history first and then using a Riedel system for the contem- ety of geological settings. Significant advances in regional porary fault zones ([43] evidenced a Riedel system for the flow system analysis have been made through the application Variscan fault zones). of 3D groundwater flow models [52]. Considering that Following these main steps based on the use of the regional groundwater flows enable the approximation of the geological structures and structural expertise, the final geom- anomalies of the geothermal heat, we used a flow model to etry of the DFN is composed of 351 fault zones in the region locate the preferential discharge areas. The flow model is of interest (generated with the 3DEC™ software, Figure 3). built with the conviction that the patterns of basin-scale 6 Geofluids 130 km Mannheim Speyer Landau Karlsruhe Soultz-sous-Forêts 150 km Strasbourg 10 km Z: up Y: north X: east Figure 3: Final DFN considered in the models. groundwater flow are mainly influenced by (1) the relief of Sedimentary deposits are modelled by 3D elements. We the water table and (2) permeability heterogeneities of the used the results provided by a transnational database [36] basin’s rock framework (more specifically, permeability con- to build a geometric model of the basin. The thickness of trasts associated with faults). With that in mind, the geome- the sedimentary layers ranges between 2 km and 6 km (in try of the 3D fault network and its hydraulic connections the eastern part), with an average value of 4 km. As a result, with 3D elements (sedimentary units) need to be treated with a 3D block domain (composed of tetrahedral elements) is special care. The modelling of a regional groundwater system built from digital elevation models (DEMs) of the basement requires large sets of data, stored in various forms and scales top surface and the surface topography of the domain area. (maps, tables, spreadsheets, etc.). Thus, a GIS (geographical Note that the 3D block domain is discretized into more than information system) software was used to improve the pre- 10,000 elements (Figure 4(e)). and postprocessing of the data (ArcGIS). In the following subsections, we describe the main steps to build the flow 4.2. Hydraulic Properties. The 3D fault network (Figure 4(d)) model at the URG scale and the simulations realized to sim- and 3D block domain (sedimentary cover, Figure 4(e)) have ulate the flow paths related to geothermal resources. The 3D to be physically connected (to allow flow transfer). To do flow model was built using the 3FLO software (Itasca®), that, nodes are added at the intersection between the pipes which is designed to solve groundwater flow and transport and faces of the 3D elements. These two types of components equations in porous or fractured media (see the description (1D-3D elements) were joined in accordance with the tec- of the code in Appendix A: 3FLO Numerical Code). tonic history. Following this, we obtain a 3D model of the URG. The 4.1. 3D Geometry: DFN and 3D Volumetric Elements. The rock matrix of the basement (between faults) is considered starting point to build the 3D DFN flow model is the 3D impermeable, and the fluid flow takes place within faults. regional fault network of the URG developed in Section 3. As a first approximation, the hydraulic properties of the fault -7 2 We used information provided by the 3DEC™ model (rele- sets are uniform, with a transmissivity of 10 m /s. The sed- vant data for each fault plane: strike, dip, coordinates of the imentary cover (3D continuum elements) is considered per- mass centre, and intersections with other planes) to repro- meable and hydraulically connected to faults (pipes) in duce the same 3D network in 3FLO. In 3FLO (Appendix A: accordance with the tectonic history. The sedimentary cover 3FLO Numerical Code), the flow within faults is supported of the URG is composed of different lithology units. As a first by 1D elements (called pipes). These 1D elements are used approximation, we considered two main different units to build a 2D grid, using two sets of pipes that are each regu- where the different units are combined according to their larly spaced. Both sets are oriented perpendicularly, with ori- permeability values. The values of the physical parameters entations of 0 and 90 (with respect to the plane dip), used in the model such as the permeability were estimated respectively, forming a grid, which is generated on each flow from geophysical investigations and hydraulic tests made on different geological units [53, 54] and, for instance, have plane. According to the 3DEC structural model, the whole model domain of the URG for the flow model is a parallele- been used in the hydrothermal model of [13]. The deepest piped measuring 130 km in width (west to east), 150 km in one corresponds to Keuper, Muschelkalk, and Buntsandstein -14 2 length (north to south), and approximately 10 km in height units with a permeability of 4.10 m and a porosity of 0.17 (depending on the surface topography). (see the blue layer in Figure 4(e)). The second one (above, see Based on an analysis of the fault sets (F1 to F5) in line the gray layer in Figure 4(e)) corresponds to Jurassic, Oligo- with the knowledge of the regional tectonic history (Section cene, Miocene, Pliocene, and Quaternary units, considered -18 2 3.2 and Figure 4(d)), each fault has been assigned to a fault less permeable with a permeability of 1.10 m and a poros- set as a function of its orientation. ity of 0.15 [13]. Note that the selection of 3D elements in Geofluids 7 Geological data & 3D numerical flow model Particle tracking conceptualization Geometry Boundary conditions Steady‐state groundwater flow 3DEC DFN (a) DEM topography Input of particles (d) Initial positions of particles (g) Hydraulically Tectonic history favourable Regional fault network Resampling areas (b) 3D DFN Regionalization Discharge areas of Simulation Flow paths deep regional ground water loops Connection based on Fixed head boundary (h) tectonic history conditions Steady‐state (e) Flow path analysis xN groundwater flow 3D volumetric (c) DEMs (f) elements Elaboration of hydraulic criteria related to geothermal anomaly characteristics Sedimentary units Elaboration of indicators (i) Time (ii) Distance (iii) Depth (iv) Flow direction (v) Concentration of relevant particles Step 3 Step 4 Step 1 Step 2 Figure 4: Main steps to build the flow model and locate preferential discharge areas at the URG scale. order to apply the respective hydraulic properties is done groundwater flow system is simulated. Under steady-state conditions, the mathematical formulation of the model is using the DEMs provided by the transnational database [36], such as the position of the base and top of main geolog- highly simplified. The hydraulic conductivity is the only ical units. required parameter, contrary to other parameters such as the storability and specific yield related to transient responses, which are not necessary. 4.3. Initial and Flow Boundary Conditions. A hydraulic no- To study flow paths, a particle-tracking methodology that flow boundary condition is specified for the lateral (north, accounts for advective and dispersive transports is intro- south, east, and west) and bottom (z = −10 km) model duced (Figures 4(g) and 4(h)). Particles (over a thousand) boundaries. For the top of the model, constant hydraulic are positioned close to the border faults of the URG at the heads are imposed (fixed). The relief of the water table at top of the model that can be possibly identified as recharge the regional scale can be either topography-controlled or areas [9, 13] (see Figure 4(g)). A series of simulations with recharge-controlled [55]. Here, we assume a topography- different values for the random number generator is carried controlled water table system, which is justified by the fact out to capture the different possible flow pathways. For that the system is in a humid region with a subdued topogra- the flow system modelled, several thousand flow paths are phy and low hydraulic conductivity [56]. The current topo- constructed. The simulated time is chosen large enough to graphic elevation has been regionalized (averaged by zone), allow particles to move out of the system or to be stuck in and the corresponding hydraulic head values have been a zone. At the end of each simulation, an ensemble of flow applied as boundary conditions. As an initial condition, the paths (from recharge to discharge zones) is constructed hydraulic heads correspond to the minimum valley centre’s (Figure 4(h)). For each particle, the following characteristics elevation (100 m). are extracted at each time step: position (x, y, and z), absolute time (t), properties of the element containing the particle 4.4. Numerical Simulations to Assess Favourable Areas such as its type (1D or 3D for fault or sedimentary unit, 4.4.1. Particle-Tracking Considerations. Once the flow respectively), its permeability, the name of the domain (fault model is built, the steady-state dynamic equilibrium of the sets or lithology units), etc. 8 Geofluids distance). The aim is to highlight the particles that cover 4.4.2. Elaboration of Criteria and Indicators. To use the results of the particle-tracking method (flow paths) for geo- greater distances. A weighting factor of 2 is attributed to each thermal resource characterization, we need to choose rele- particle covering a distance of more than three-quarter of the maximal distance, which represents around 1% of the vant hydraulic criteria that constitute a good signature of the geothermal anomalies. Then, based on these main char- total particles. A weighting factor equal to 1 is attributed acteristics of geothermal anomalies (hydraulic criteria), indi- to each particle covering a distance of more than one-half cators are elaborated and extracted from each particle and less than three-quarter of the maximal distance, which trajectory (Figure 4-Step 3) to use the results of the simula- corresponds to 14.5% of total particles. By the same reason- ing, two threshold values are calculated related to the max- tion that enable the identification of the geothermally favour- able areas. The following assumptions are used: imal travel time; then, a weighting factor of 2 is attributed to each particle with a travel time of more than three-quarter (1) Geothermal anomalies result from a state of imbal- of the maximal travel time, and a weighting factor equal ance over a long period where the only efficient to 1 is attributed to each particle having a travel time of mechanism (ubiquitous and active over large spectra more than half and less than three-quarter of the maximal of time and space) is the regional groundwater flow. travel time. For the indicator related to the depth reached Then, the flow path length (long enough to be by the particle along its trajectory, we considered that a par- regional) and the time required to cover this distance ticle, which has reached a sufficient depth of more than are used as the two main indicators to highlight 2.5 km (related to the average geothermal gradient), is inter- regional groundwater loops esting and deserves a weighting factor. Then, a weighting factor of 1 is attributed to each particle that reached a depth (2) Geothermal anomalies reflect positive thermal anom- of more than 2.5 km along its trajectory. To finish, the sum alies resulting from upward flow in discharge areas of weighting factors is calculated for each particle. The idea due to ascending warm waters. Therefore, the particle is that particles with the higher scoring are likely to produce must come from a deeper part and has to reach a suf- thermal anomalies and correspond to the more relevant ficient depth along its trajectory to carry warm(er) particles and inversely, the lower scoring to the less relevant water until the upflow area. The depth reached by particles. To eventually constitute a preferential target area the particle along its trajectory, the final position (or (thermal anomaly), a final requirement has to be fulfilled: a position slightly upstream from it) of the particle asufficient density of relevant particles must be localized along its track, and the flow direction are thus also in a “restricted” spatial area. used as main indicators 4.5. Favourable Area Estimates: The Flow Model. Figure 5(a) As mentioned just above, the indicators that are used presents the results of the flow path analysis at the URG scale are (1) the total distance travelled by the particle from where particles are plotted at a 2 km depth (equal to that of the recharge area (starting point) to the discharge area the available temperature map, see Section 6.2) as a function (ending point) during time t, (2) the time spent by the of the weighting factors. Particles with a longer distance and particle in the system to travel between the recharge and time (in accordance with regional groundwater loops) and discharge areas, and (3) the depth reached by the particle reaching a depth over 2 km along the trajectory (in accor- throughout the trajectory. Note that the relevant informa- dance with the temperature) are selected. The target areas tion associated with each track is extracted to be analysed are then delineated by defining the contours encompassing in a GIS software. these relevant particles (Figure 5(a), green lines). Five differ- 4.4.3. Flow Path Analysis. The main objective of the ground- ent preferential areas are highlighted from the results of the water flow path analysis is to cross-correlate the properties of 3D flow model, which are distributed as follows: two areas each particle trajectory (distance, time, depth, localization, to the north of the URG (surrounding Speyer and between and flow direction) to locate the preferential discharge areas Landau and Mannheim), one area clearly located to the (with specific features of geothermal anomalies) of regional south of Soultz-sous-Forêts, and two areas located near groundwater loops. and south of Strasbourg (Figure 5(a)). Figure 5(b) presents After an individual flow path analysis, a spatial and qual- the trajectories associated with the most relevant particles itative analysis is realized to localize the preferential dis- providing the preferential areas. The northern areas are charge areas of regional groundwater loops. A preferential mainly related to particles from the eastern shoulder of the discharge area corresponds to a zone with a high concentra- URG and along a SE-NW direction. The second area to tion of relevant particles (related to the deep regional flow the south of Soultz-sous-Forêts is mainly related to particles circulation), which is likely to present thermal anomalies from the western edge of the URG and moving in a NW-SE (Figure 4-Step 4). To select relevant particles, a preliminary direction. Note that few particles come from the eastern edge data analysis is performed. The statistical distribution is and in a SE-NW direction. The southernmost areas (sur- observed to define threshold values that are used to assign rounding Strasbourg) are related to particles from the east- weighting factors to particles. For example, to select particles ern and western edges. that correspond to regional groundwater loops (related to a The results of the 3D flow model provide preferential distance indicator), two threshold values are determined discharge areas related to geothermal resources. Since the (equal to one-half and three-quarter of the maximal ultimate aim is to support expensive exploration campaigns, Geofluids 9 Mannheim Speyer Landau Soultz-sous-Forêts Strasbourg Classification of particles Trajectories associated with most relevant particles (see from above) Main flow directions related to relevant particles More Less relevant relevant Favourable area km Fault network (a) (b) Figure 5: (a) Classification of particles resulting from the flow path analysis based on selection criteria and the respective delineated preferential target areas obtained (favourable areas are delineated by the green circles). Note that the particles are plotted at a 2 km depth upstream from the ending point of their trajectory. (b). Particle trajectories of the most relevant particles (red crosses) and the associated main flow directions. any information improving the confidence in the area loca- ([29, 58, 59], among numerous others) and the segmentation tion is valuable. When no in situ data is available, cross- of the medium (the more fault zones, the more deviation). validating the interpretation with results from other models To build a 3D stress model of the URG, we use the code is a way to increase the confidence level. In this study, we pro- 3DEC™ [35]. pose to build a 3D stress model to extrapolate preferential areas from the mechanical perspective and to finally compare 5.1. Mechanical Properties. A full description of the code and its constitutive equations is provided in Appendix B: them to the flow model estimates. 3DEC Numerical Code, and only a short outline is given here. Since we focus on the presence of the fault zones, 5. 3D Stress Model at the URG Scale the matrix behaviour is kept simple and assumed to be elas- tic (Table 1). The fault zones, however, follow an elastoplas- tic law (Coulomb slip), with a possible dilation (Table 2). In We propose to investigate the heterogeneities in the stress distribution that may favour the setup of thermal anoma- addition, we assume that any large deformation occurring lies. The stress distribution in rock masses is significantly within the medium is localized into the joints and not into affected by structural and material heterogeneities (including the rock matrix (see Appendix B: 3DEC Numerical Code mechanical behaviours/parameters and fault zones) and for more details). Due to the absence of data at the considered scale, the anthropic activity. Yale [57] gives a review of the in situ evi- dence of fault impact on the local deviation of stress from the parameters are extrapolated from estimates obtained by fit- far-field trend. The conclusions are that the main factors ting numerical models on in situ experiments [60, 61]. We locally affecting the stress are the proximity to the fault zones extrapolated the parameters qualitatively considering fault 10 Geofluids Table 1: Mechanical parameters for the rock matrices. 5800000 Notation Value Value Parameter (unit) (sediment) (basement) −3 –1 Density ρ (kg·m ) 2400 2680 0.2 mm·y Bulk modulus K (GPa) 16.7 42.8 Shear G (GPa) 7.7 22.1 modulus Table 2: Mechanical parameters for the fault zones. Parameter Notation (unit) Value −1 Normal stiffness k (GPa·m)10 −1 Shear stiffness k (GPa·m)5 Cohesion c (GPa) 0 –1 Internal friction angle ϕ ()15 0.2 mm·y Dilation angle ψ ()5 –1 0.56 mm·y Critical shear displacement u (m) 1 200000 300000 400000 500000 600000 Easting (UTM, m) zones with noncohesive infills (weaker elastic modulus, zero Figure 6: Far-field displacements and embedding model (grey) for cohesion) and large wavelength profiles (smaller dilation setting the faulted model (red frame) boundary conditions. angle, larger critical shear displacement) that are altered by Modified after Buchmann and Connolly [62]. UTM: Universal numerous tectonic stages (smaller friction). The parameters Transverse Mercator. are on the orders of magnitude of the parameters from models built at similar scales in the URG [62, 63]. are applied to the model boundaries, while the gravity is still 5.2. Boundary Conditions. The displacement boundary con- active within the model: ditions are set in accordance with the present-day tectonic context of the URG and are derived from the model of Buchmann and Connolly [62]. The faulted geometry pre- z = − ρgdz , v d sented in Section 3 is embedded in a continuous material 0 having the dimensions of Buchmann and Connolly’s model S z = S z = S z , H max d H min d V d (Figure 6). Following a procedure used by Peters [64], this 1 − v embedding enables simulating the far-field nature of the boundary conditions. As a first approximation, the embed- where S (Pa), S (Pa), and S (Pa) are the vertical, V H max H min ding material is assumed to have the same properties as maximum horizontal, and minimum horizontal principal −2 the crystalline basement. In addition to these lateral bound- stresses (respectively); z (m) is the depth; g (m·s )isthe −3 ary conditions, the normal displacement is fixed at the vertical component of the gravity acceleration; ρ (kg·m )is model base (roller-type boundary condition), and the top the rock mass density, and ν (−) is Poisson’s ratio. is left free. 5.4. Simulations to Assess Favourable Areas. Once the 5.3. Initial State. In addition to the present-day forces, the assembly of blocks and joints is balanced (see Appendix stress field is affected by past tectonic events that leave B: 3DEC Numerical Code for more details), the model residual stresses through the rock mass. The prestressing gives an estimation of the mechanical equilibrium through method is commonly used to simulate the impact of past the system. In particular, displacement and stress heteroge- events [62, 65–67]. Assuming that the previous tectonic epi- neities are obtained and can be inspected to highlight areas sode occurred a sufficiently long time ago to allow the relax- where the mechanical state can favour the presence of ther- ation of the medium, the initial state can be approximated mal anomalies. As previously, this requires the definition of using a gravitational loading. Gravity is activated within the an appropriate indicator. The definition of a mechanical blocks, and a uniaxial stress condition (equation (1)) is pre- indicator is difficult since stresses can act on flows in sev- scribed on the lateral boundaries of the model to account eral ways. However, since stresses can act on the ground- for the influence of the surrounding rock mass. Note that water heads directly (impact on pressure gradients) or on due to the high density of the fault zone network, the initial the pore/fracture network (impact on medium effective condition makes the stress distribution very heterogeneous permeability), most mechanical indicators can be catego- through the rock mass. Once the initial equilibrium is rized following the volumetric and deviatoric decomposi- reached, the conditions described in the previous section tion of the stress tensor: Northing (UTM, m) Geofluids 11 Geological data & 3D numerical geomechanical model conceptualization Geometry Initial stress state 3DEC DFN Elaboration of stress criteria Gravity balance related to geothermal anomaly characteristics Elaboration of indicators - Mean stress Tectonic history Regional fault network 3D DFN Boundary conditions Regional tectonic trend Result analysis using Geomechanically aforedefined indicators favourable areas Connection based on tectonic history 3D continuum Results of the stress model DEMs Sediment pile and basement Figure 7: Main steps of the geomechanical simulations (from data to delineation of preferential areas). (i) The volumetric part of the stress tensor (e.g., mean Low values of the indicators are thus used to identify favour- stress) is directly linked with the groundwater heads: able areas. flow will happen from the more compressed areas to The different steps to delineate the favourable areas with the less compressed zones [68]. “Volumetric” indica- the geomechanical model are summarized in Figure 7. tors could also relate to modifications of the pore net- 5.5. Favourable Area Estimates: The Stress Model. The equi- work (e.g., fracturing), but since we consider hard librium of the mechanical system is computed through rocks under compressive states, the pore network 3DEC™. We obtain a stress map that can be used for delin- modification will mostly happen under the deviatoric eating the preferential target areas. Prior to postprocessing influence the stress, the depth at which we search for favourable areas (ii) The deviatoric part of the stress tensor does not must be set. Rather than a fixed depth, we decided to inter- affect the pressure heads but can affect the effective pret the results close to the sediment/basement interface. permeability of the medium. In the rock matrix, Since the DFN was built according to the fault zone traces the stress deviator can be responsible for shear- observed at the top of the basement, the level of uncertainty induced fractures that greatly contribute to the in the geometry is more limited close to the interface. The permeability increase of the medium. For fault mean stress distribution is presented in Figure 8(a) (the zones zones, deviatoric stresses result in shear displace- located in the horst are ignored). ments perpendicular to which upflow is favoured Five different preferential areas are highlighted from the [24]. In addition, the shear stress maintains the results of the 3D stress model. They are distributed as follows: fault zones in an active state, which is argued to two areas in the northwestern part of the URG (to the west reduce the fault sealing and hence favour the per- of Speyer), one area located around Soultz-sous-Forêts, and meability [69] two areas located in the southern part (one to the south- west of Strasbourg and one to the east) and to the south In the first study, we selected a single scalar indicator of Strasbourg (Figure 8). to be chosen according to either the volumetric or the deviatoric part of the stress tensor. Denoting σ (Pa) the 6. From Regional-Scale Models towards the stress tensor (see Appendix B: 3DEC Numerical Code for more details), we chose the mean stress σ = Tr σ /3 as Prediction of Favourable Areas for the mechanical indicator: by releasing the efforts acting Geothermal Resources on the rock mass, the less compressed areas (i.e., low σ ) 6.1. Cross-Analysis of Favourable Areas. The final step of the facilitate fluid discharges, including vertical upflows of fluids that were possibly heated up to the rock mass below. work is to estimate areas where the physical conditions ++ 12 Geofluids σ compared to the existing thermal map of the URG at a 2 km depth. From north to south, three main local thermal Mannheim anomalies can be distinguished from the temperature map: one in the north, between Landau and Speyer; one in the cen- Speyer tre, around Soultz-sous-Forêts; and one in the south, to the west of Strasbourg (Figure 9(b)). Then, from ~east to west, Landau the temperature map shows that the highest temperatures are preferentially located in the western half of the rift. The major local thermal anomaly is located around Soultz-sous- Soultz-sous-Forêts Forêts in the centre of the URG, where the highest tempera- tures are located to the east of Soultz-sous-Forêts and then in the centre of the basin. The second most important local thermal anomaly is located between Landau and Speyer, Strasbourg where the temperature decreases from south (Landau) to north (Speyer), and the third local thermal anomaly is located to the west of Strasbourg (Figure 9(b)). From a global perspective, the cross-analysis of the results of the stress and flow models delineates three distinct Fault network km favourable areas: one in the north (west of Speyer), one in the centre (south of Soultz-sous-Forêts), and one in the south Figure 8: Mean stress distribution in the model and the delineated (south of Strasbourg). This distinct distribution (from north preferential target areas obtained. to south) in three main areas is in good agreement with the distribution of the main local thermal anomalies known in the URG (Figure 9). favour the setup of geothermal resources. Favourable areas were delineated for both key physics (flow and stress models) based on related selection criteria defined in accordance with 6.3. Discussion of the Cross-Analysis Results the geothermal characteristics and the physics considered. The most favourable areas are assumed to be at the intersec- 6.3.1. Accuracy of the Areas Delineated through the Cross- tion of flow-favourable and stress-favourable areas. They are Analysis. From a global perspective, the results obtained from estimated qualitatively by overlapping the different favour- the cross-analysis are very encouraging, with the delineation able areas obtained with each model (flow and stress models) of three favourable areas distributed from north to south (Figure 9(a)). along the rift, which is in good agreement with the distribu- The cross-analysis highlights three main favourable tion of the main local thermal anomalies known in the areas: one is in the northern part of the URG, close to Speyer; URG. However, the accuracy of the results against the ther- a second one is to the south of Soultz-sous-Forêts, and a third mal map varies from one delineated area to another. From one is to the south of Strasbourg (composed of three different north to south: subareas) (Figure 9(a)). The northern delineated area (west of Speyer) results from the overlapping of the two main (i) the area between Landau and Speyer partially high- favourable areas of the stress (blue circle orientated in the lights the local thermal anomaly. The cross-analysis east-west direction) and flow (green circle orientated in the results delineate an area to the west of Speyer, NS direction) models. Close to Soultz-sous-Forêts, the two approximately 20 kilometres away from the in situ distinct areas (blue circle orientated in a NE-SW direction area and green circle orientated in a NW-SE) result in one main (ii) the area close to Soultz-sous-Forêts is the best favourable area. Note that these two main favourable areas correlated one by the model. Still, the final estima- (west of Speyer and south of Soultz) are located in the west- tion is to the south of Soultz, while the actual area ern part of the rift (Figure 9(a)). For the four areas (two is most visible to its east respective areas for each model corresponding to two green and two blue circles, respectively) located to the south of (iii) the model results for the area close to Strasbourg are Strasbourg, the cross-analysis results in three subareas shifted eastwards compared to the thermal map. (including two small ones), where the largest one is orien- However, as opposed to the two other areas, the in tated north-south in the western part of the URG situ data in this area is much more diffuse. The (Figure 9(a)). To summarize, the cross-analysis highlights model results follow the natural tendency of the two main favourable areas that are located in the western part thermal anomalies, and in this regard, we conclude of the URG, to the south of Soultz-sous-Forêts, and to the that this area is not the least precise one west of Speyer. A third area can be distinguished to the south of Strasbourg, where three subareas are delineated. The uneven precision of the model results steams from the cross-analysis of the stress and flow results. While the 6.2. Delineated Favourable Areas versus the Thermal Map. two hydraulically favourable areas tend to correlate with the actual thermal anomaly located between Landau and Speyer The main favourable areas from the cross-analysis are Compression Geofluids 13 Mannheim Mannheim Speyer Speyer Landau Landau Soultz-sous-Forêts Soultz-sous-Forêts Strasbourg Strasbourg (Temperature map from GeORG) Temperature at 2 km depth Stress Flow Value High: 167 Cross-analysis Low: 61 km Fault network (a) (b) Figure 9: (a) The results of the cross-analysis based on the overlapping of the respective favourable areas issued from the stress and flow models. (b) Map of temperature at a 2 km depth obtained from the GeORG project. (Figure 9), the mechanical results are less precise and move reasonably well estimated from a regional perspective. From the final area estimation northwards. In the end, the northern the 3D flow model, this favourable area results from particles delineated area partially highlights the local thermal anom- originating mainly from the western edge of the URG aly. Conversely, the mechanical results appear to better high- (located to the northwest of Soultz-sous-Forêts, see the large blue circle Figure 10(a)) and in a main NW-SE flow direction light the thermal anomaly located on the western edge of the rift than the flow model’s ones. According to the 3D flow (Figure 5(b)). The 3D cross-sectional view (Figure 10(b)) model results, these hydraulically favourable areas are related shows that particles move from the west (downward flow) to particles originating from the eastern border (see yellow along the URG east-dipping boundary fault to the east circle, Figure 10(b)) and are characterized by regional flow (upward flow), mainly along west-dipping faults and associ- ated subvertical faults that are connected to these west- pathways that mainly occur through NW-SE faults (F4 fault set) in a global direction SE to NW (Figure 5(b)). These few dipping faults (Figure 10(b)). These results suggest deep particles go through the centre of the URG (where the sedi- regional groundwater loops that occur through east-dipping mentary basin reaches more than 4 km deep) and then (downward flow) and west-dipping (upward flow) faults that migrate through the regional fault network up to the favour- intersect at great depth and along the main NW-SE flow direction and where the main recharge area is located to the able areas. These deep regional flow paths appear consistent with the geochemical evidence [14]. Indeed, the geochemical northwest of Soultz-sous-Forêts. These preferential flow signatures of geothermal brines collected from the granite pathways along the west-dipping fault zones characterize basement (Landau, Rittershoffen, Soultz, and Insheim at the horst structures and are consistent with some previous works graben’s NW borders) suggest that they all reacted with deep [9, 70]. Moreover, according to the results (particle flow paths), we can observe deep flow paths (few particles) from sedimentary rock at temperatures close to 225 Cat depths ≥ 4 km [14]. This signature involves a migration of geother- the SE (small blue circle in Figure 10(a)) and in a global mal brines from the centre of the URG to the graben’sNW SE-NW flow direction emerging close to the Soultz and borders through a complex system of deep faults. Rittershoffen areas. These few particles go through the centre For the centre area, the delineated area (close to Soultz) of the URG (where the depth of the sedimentary basin reaches more than 4 km) and then migrate through the that corresponds to the major thermal anomaly known is 14 Geofluids a. Areas of origin for main favorable areas km Fault network South of Soultz Favourable areas Between Mannheim and Speyer (hydraulic) Around Strasbourg b. Soultz Temperature at 2 km depth Trajectories associated with relevant particles Main flow direction Value Fault network High: 167 Start/middle/end Low: 61 Figure 10: (a) Plan view of trajectories associated with relevant particles and areas of origin of particles associated with favourable areas. (b) 3D cross-view (south-side view) of trajectories in the favourable area located close to Soultz-sous-Forêts according to the blue rectangle drawn in (a). regional fault network up to the favourable area (which could analysis because it presents the key characteristics from a also be consistent with the results of [14]). The favourable hydraulic point of view and also from a mechanical point area close to Soultz is the best located one by the cross- of view. Indeed, the results of the mechanical model show a Geofluids 15 –4 network affects the redistribution of stresses and regional groundwater fluid circulation. Additionally, the stresses sig- Δu (m) nificantly affect the flow properties, and there is a strong cor- relation between the regional groundwater flow and geothermal anomalies. The groundwater flow paths are Opening mainly influenced by the regional fault network (structural assemblages and heterogeneity of permeability created by faults) as well as by the relief of the water table. Thus, the key factors are the faults and permeability heterogeneities [25], which are influenced by the stresses (a key physic) that control both the fluid circulation (a key physic) [6, 25] and the location of hot fluid upwelling. The respective favourable areas delineated by the flow and stress models demonstrate the strong link between the deep regional fluid circulation Closing within the regional fault network, the stresses, and the distri- bution of the geothermal resources. These previous observa- tions lead us to ask the following question: is it essential to –4 –10 consider the thermal physical processes to highlight the ther- Stress Opening mal anomalies in the context of geothermal exploration where the ultimate goal is to limit the area to explore and Figure 11: Normal displacement variations of fault zones and reduce the associated costs? Building a hydrothermal model corresponding favourable areas (in purple). Blue circles are a at the regional scale could facilitate the interpretation of the reminder from the mean stress interpretation. results by providing the distribution of the temperature and consequently local thermal anomalies. This type of model could serve to avoid approximations that can appear through less compressed region (Figure 8(b)) where the stress state the definition of criteria and then have an impact on the can favour upflow. Thus, the major thermal anomaly delineation of the favourable areas. Thermal processes such known in the URG is the best located one by the numer- as the heat transfer between faults and the surrounding ical approach due to the combination of numerous key sit- matrix and at the intersection of the base of the thermal blan- uations: structurally, hydraulically, and mechanically. ket and the fault (in which the upward transfer of heat Structurally, this region corresponds to a horst with a occurs) result in redistributions of heat. These redistributions west-dipping fault (enabling the intersection at depth with of heat could provide slightly shifted favourable areas in com- east-dipping faults). Hydraulically, this area corresponds to parison with those currently delineated in this study. Thus, a preferential discharge area of deep regional fluid circula- this type of numerical model could be relevant and, more spe- tion. Geomechanically, the mechanical stress in this region cifically, useful to move towards a better overall understand- favours fluid circulation and then hot fluid upwelling in ing of the URG geothermal systems and should provide some the rock matrix and fault. constraints on the relative influence of each physics. In any To finish, for the southern area, the local thermal anom- case, building a hydrothermal model at the regional scale aly is particularly well delineated by the geomechanical constitutes a real challenge and then requires further analysis. model (less compressed area and faults exhibiting opening The promising results also demonstrate the key role of the tendencies, see Figure 8). Second, preferential pathways regional fault network, which is at the basis of the two numer- (from the western border fault) that occur along the west- ical physical models. Then, a better delineation of the favour- dipping fault zone also characterize the southern area (to able areas requires an improvement of each numerical model the west of Strasbourg), and we can observe particles from and thereby an improvement of the DFN geometry. the eastern border fault that emerge in this favourable area (see red circles, Figure 10(a)). All of this information seems (2) Cross-Analysis versus Coupling. One of the main assump- consistent with the knowledge provided by the numerous tions of the study presented in this paper is that the estima- studies realized within the URG (geophysics, geochemistry, tion of favourable areas can be obtained by cross-analysing etc.). The results demonstrate the potential of this type of results from separate single-physic models. Several authors numerical approach to provide constraints for the challenge use either one-way coupled or fully coupled models as an of geothermal exploration and to move towards a better over- attempt to more accurately address the physical features all understanding of the URG geothermal systems. involved in large-scale phenomena (>kilometric characteris- tic length) [71–75]. 6.3.2. Relevancy of the Model Cross-Analysis In our case, a solution for achieving one-way coupling (1) Choice of the Physics Described by the Models. The pro- would be to use mechanical results for parameterizing the posed approach is based on the choice of the single physics groundwater flow model. Indeed, mechanical results predict described by the models. Here, we propose an approach huge heterogeneities in stress distribution, including in its deviatoric part (which can be related to second-order based on these main underlying hypotheses: the regional fault 16 Geofluids carried out by hierarchizing them based on their impact on fracturation). As a result, the rock mass is expected to exhibit heterogeneous permeabilities (increase of permeability with the results. Most notably, the boundary conditions and the fracturation). More complex effects could also occur in the use of nonheterogeneous parameters for the groundwater flow and mechanical models, respectively, are expected to fault zones such as permeability decreases in fault zones that are less shear active [69] or increases in permeability in the be mostly responsible for the uneven accuracy of the areas direction perpendicular to the shear displacement [22, 24]. they delineate (see Section 6.3.1.). Full coupling for the considered geometry (large scale Nevertheless, the reasonable agreement of the model and DFN) seems, however, to be unrealistic with the tools results with the thermal observations leads us to believe that the physical model accuracy, both in terms of state equa- at hand today. The fully coupled models found in the litera- ture consider much more reduced scales [74, 75]. Conversely, tions and parameters, is a secondary aspect in the approach. models exhibiting scales comparable to ours (~100 km char- Rather, the geometry description seems to be of more acteristic length) incorporate either only a single physics importance, as is usually the case with models based on a and/or unfaulted media [26, 62, 65, 76]. CPU- or GPU- dense DFN. based parallel simulations show promising results for han- dling fully coupled models but need further progression to 6.4.2. Geometry Description. The DFN geometry (regional reach the scale and the medium segmentation considered in fault network) requires a more detailed description. The geo- this study [71–73]. dynamic evolution of the rift system and existing conceptual At this stage of development of the study, we anticipate models of the rifting zones are important sources of informa- tion to improve the DFN construction [77–79]. The interpre- the use of coupled models to introduce more uncertainties than improvements in the overall approach. Given the tation of the two deep seismic lines illustrates the structural absence of information on coupled phenomena at large spa- asymmetry of the Rhine Graben, which is characterized by tial and time scales, the coupling laws extrapolated from thickness variations of the tertiary deposits and the presence smaller scales might be erroneous, and the estimation of the of boundary faults in the north along the eastern flank and in the south along the western flank. The asymmetric evolution coupling parameters would be very difficult (even qualita- tively). In contrast, focusing on the prevailing processes of the rift from north to south that affects the deep structural related to the setup of geothermal resources for each model features has to be considered. Similarly, horizontal bedding is expected to better constrain the results of the approach planes and a simple conjugate fault pattern characterize the (i.e., positions of the favourable areas). In our opinion, the Cenozoic fill in the north, while in the south, a more complex structural pattern exists, with the presence of tilted fault cross-correlation of relatively well-trusted models brings more objectivity in the results than estimation stemming blocks from west to east, conjugate normal faults, and a from a coupled but more uncertain model. The case study smooth anticlinal structure. Then, the steeply dipping nor- mal fault that controls the western margin of the Cenozoic presented in this paper proves that the cross-analysis tend to estimate reasonably well the thermal map obtained from graben and merges at the midcrustal levels provides con- straints on the fault system and its evolution with the depth. in situ data and demonstrates how the overlapping results from uncoupled models can be used to refine the final esti- Building the geometry by integrating this type of constraint mates of the favourable areas. could directly affect the connectivity of the medium as well as the stress redistribution due to the model segmentation 6.4. Perspectives of Improvement for the Proposed Approach. and consequently the groundwater flow paths and the geo- graphical distribution of the discharge areas. Even if there is acceptable agreement between the global distribution of the preferential areas delineated and the major thermal anomalies, we can notice some differences 6.4.3. Result Analyses. The final estimations of the proposed approach depend highly on the way in which the model with the temperature map at a 2 km depth, which underlines the need for improvements that will constitute the next results are interpreted. First, considering the two models sep- steps of this work. Based on the results and discussion, we arately, the relation between the model results and the char- conclude that improvements can be made at three main acteristics of the geothermal resource rely on the definition of levels: the geometry description, model precision, and result the indicators for each model. Recalling that a particle- tracking method is used for the groundwater flow, the (cross-)analyses. From this perspective, the following para- graphs detail the opportunities for improvement and their retained indicators are the residence time, the travelled dis- relative importance. tance, the depth reached by the particles, the flow direction, and the concentration of relevant particles. The first three 6.4.1. Physical Model Formulations. Both the groundwater indicators are linked with the temperature (fluid heating up the rock mass), the fourth one indicates resource upwelling, flow and mechanical models rely on several simplifying assumptions and could be refined in several ways (equations, and the final one relates to the discharge areas. Because of parameters, boundary conditions, and initialstates). However, the rather straightforward effect of flows on the thermal building almost exact models is an unachievable task given anomalies, we are rather confident in the chosen indicator. the large spatial and time scales involved in light of the However, the relative interpretation of the particles could absence of parameter measurements and evolution of tec- be improved through a rigorous statistical treatment. In this tonic states and related mechanical boundary conditions paper, and as a first approximation, we chose to set the [8]. Targeting a limited number of improvements could be weighting factors qualitatively to focus on the overall Geofluids 17 systems might have to be set and run to progressively refine approach. Regarding the mechanical model, the impact on the thermal anomalies is only secondhand, through affecting the estimation precision. If the proposed approach is used the flow and/or thermal processes. In this paper, we retained as a support for an exploration campaign whose aim is to site the mean stress σ , arguing that less compressed areas favour exploration wells, the results of a first region-scaled model fluid migration. The mean stress seemed to be a reliable indi- might be too diffuse or uncertain by, e.g., providing several distinct areas in the same region. If more precise results are cator for identifying the true thermal anomalies. We remind here that, in accordance with the whole approach, the expected (e.g., if the range of the final estimates is too large mechanical criterion can only be interpreted qualitatively, or if there are concerns about the locations of the favourable i.e., so as to highlight less compressed areas close to more areas), steps one to four could be repeated on progressively compressed ones. In such cases, fluid migration is assumed smaller-scaled geometries until the desired size is achieved or more confidence is obtained for the final estimates of the to move towards less compressed areas, then move upwards due to stress reduction when getting closer to the Earth’s sur- favourable areas. In this case, the results from larger-scaled face. In our opinion, a quantitative analysis is not realistic at models can provide boundary conditions for the smaller- this stage given the lack of information on the rock mass scaled ones (see, e.g., [76]). properties at the considered scale (law of behaviour and asso- ciated parameters, behaviour heterogeneities, second-order 7. Conclusion fracturation…). σ cannot be guaranteed to be the best one. Other indicators were tested such as the preferentially This paper is a current overview of our team’sefforts to inte- opened fault zones. Taking advantage of the discrete nature grate the impacts of the key physics as well as key factors con- of the model, the opening tendency can be estimated through trolling the geothermal anomalies in a fault-controlled the fault zone normal displacement u (m). More specifically, geological setting in 3D physically consistent models at the and to focus on the impact of the present-day stress state, the geological system scale. The study relies on the building of init displacement increment can be plotted: Δu the first 3D numerical flow (using the discrete-continuum = u − u , n n n init where u (m) is the normal displacement due to the initial method) and mechanical uncoupled models (using the dis- equilibrium. Figure 11 presents the contours of the fault zone tinct element method) at the URG scale (faulted geological setting). Studying the interactions of these main processes averages Δu = 1/s ∬ Δu dS (S (m ) being the area of the n n fault zone) in order to interpret the tendencies at the fault while coupling them with the effects of the regional structural elements (fault zones) helps us to understand their relative zone scale. NW-SE faults exhibit opening tendencies in con- trast to EW and NE-SW faults that exhibit a tendency to impacts on the distribution of geothermal resources. Here, close. Estimated favourable areas tend to be more widespread we applied this approach in the URG context (Figure 12), where available data (such as seismic lines) allow the building (light blue circles in Figure 11). Interestingly, the results in the surroundings of Soultz- of these types of models and a direct assessment of the results sous-Forêts correlate well with the σ (using, for instance, the temperature map). Based on the fact interpretation in this region, tending to reinforce the idea that this region com- that there exist strong links between the stress, fluid flow, and bines numerous key situations. Comparison of different geothermal resources, we followed these main steps. First, the key role of the regional fault network is taken mechanical results highlights the difficulty in choosing a mechanical criterion, which in the end might be a cross- into account using a discrete numerical approach. The geom- correlation between rock matrix and fault zone criteria. On etry building is focused on the conceptualization of the 3D a broader perspective, these results question how the fault zone network based on structural interpretation and mechanical results should be interpreted: in a straightfor- generic geological concepts and is consistent with the geolog- ical knowledge of the URG. This major step constitutes a ward manner (as is done in this paper) or to parameterize the flow model? Thorough investigations must be carried challenge alone and marks the beginning of new opportuni- out in the future to shed light on these interrogations. ties to understand the key effect of the regional fault network Second, the cross-analysis of the model results itself could on the geothermal resources. Indeed, this question remains be improved. Here, as a first approximation, we proposed to an open question and still constitutes a challenge for geother- set the favourable areas by choosing the overlapped regions mal engineers, who attempt to explain hydrothermal pro- that were qualitatively delineated. In the studies found in cesses in this faulted geological context. the literature (see, e.g., [80, 81]), the authors propose to Then, we decline the DFN model of the regional fault net- quantitatively establish favourability maps by interpolating work to build the first 3D flow and stress models of the URG in situ data and using weighting factors according to the using adapted numerical tools that explicitly consider the key measurement confidence. In our case, because the physical role of the faults. The respective results issued from each models are built under several assumptions, it seems slightly numerical model are indirectly interpreted through the defi- too early to propose quantitative favourability maps. The nition of criteria (using the link between the physics consid- more precise data at our disposal are the structural data. ered and the setting up of geothermal resources) that enable A first step could be to weight the results according to the elaboration of indicators to identify geothermally favour- how close (and thus less extrapolated) they are to our con- able areas related to the considered physics. For the flow fident structures. model, the hydraulic criteria are related to the preferential Finally, depending on the ratio between the obtained and discharge areas of the deep regional fluid circulation, and desired scales for the preferential target areas, additional the indicators used are related to the properties of the particle 18 Geofluids Available geological data Geological maps Conceptualization 3D numerical models Tectonic Results Interpretation Geological profiles history 3D flow Favourable model areas Elaboration of indicators Seismic Elaboration of criteria Favourable related to geothermal Cross-analysis areas Faults anomalies & physics 3D stress Borehole data model Elaboration of indicators Favourable Layers areas Digital elevation Results Interpretation models Figure 12: Methodology employed to delineate preferential target areas for geothermal resources using numerical models. Equilibrium of blocks F u ext Block displacements Joint forces affect block affect joint external forces Joint constitutive equations u F ext Figure 13: DEM solution scheme at each time step: the solver loops until u and F are balanced (inspired by Itasca, 2016). ext trajectories (distance, depth, etc.). For the stress model, the considered (stress and fluid flow) in the geothermal resource stress criteria are related to the mechanical state that can setup, and the method to interpret the results (criteria and favour the presence of thermal anomalies, and then, the associated indicators) are relevant to capturing the thermal related indicator chosen is the stress tensor (mean stress). processes and then can help to locate thermal anomalies. In As a final step, we proposed to compare the hydraulically other words, the geographical distribution of the geothermal and geomechanically favourable areas to locate thermal resources is strongly correlated with the preferential dis- anomalies. A cross-analysis of the results of the 3D flow charge areas of the deep regional fluid circulation occurring and stress models within the URG provides three preferential within the regional fault network and related to the less com- target areas for geothermal projects. These preferential areas pressed areas that favour these upwellings of deep hot fluids. are distributed along the rift as follows from north to south: Moreover, a detailed analysis of each favourable area to the west of Speyer, to the south of Soultz-sous-Forêts, demonstrates the ability of this type of approach to move and to the south of Strasbourg. The flow and stress models towards a better overall understanding of the URG geother- result on a common, rather condensed, overlap for two of mal systems. For example, the best favourable area delin- the three areas (Speyer and Soultz-sous-Forêts) which are eated by the proposed approach corresponds to the major related to major local thermal anomalies currently identified thermal anomaly known within the URG. A more detailed in the western part of the rift. The overlapping is more diffuse investigation of the delineated favourable area close to for the third area where the in situ data is more scattered as Soultz-sous-Forêts enables the highlighting of the key situa- well. Given that no parameter fitting was done, the delineated tions (structurally, geomechanically, and hydraulically) that favourable areas and local known thermal anomalies tend to make it a special location for geothermal resources. Struc- be in reasonable agreement, which underlines the potential of turally, the west-dipping fault resulting in the intersection such a cross-numerical approach. These promising results at depth with east-dipping border faults forms a key struc- suggest that the numerical approach (discrete approach to tural feature. This V-shaped structure favours loops of deep explicitly consider the effect of the fault network), the physics fluid circulation. Hydraulically, this area corresponds to a Geofluids 19 transport equations can be used to assess the favour- preferential discharge area of deep regional fluid circulation. In addition, the mechanical stress in this region favours fluid able geothermal areas circulation and then hot fluid upwelling in the rock matrix. (iii) flow interactions between fractures and porous In conclusion, the Soultz-sous-Forêts area is probably the media by allowing the superimposition of 1D ele- major thermal anomaly because this area merges structur- ments generated on disk-shaped fracture planes (to ally, geomechanically, and hydraulically favourable condi- simulate fractures) on 3D elements (to discretize tions. This evidence is consistent with previous studies and the rock matrix) demonstrates the primary role of the regional fault network. Then, the constraints provided by the respective models (iv) transport simulated by the particle-tracking method demonstrate what we can learn from these first models The flow equation solved in the 1D classical elements is and how we can improve the proposed approach. This study provides suggestions for future research on the dynamic context of the geothermal resource setup. The ∂ h ∂h C = A S , A 1 c S S very encouraging results underline the potential of such a i i ∂x ∂t cross-numerical approach to help locate geothermal 3 −1 resources. This type of approach could ensure the localiza- where C (m s ) is the pipe conductivity; h (m) the hydrau- tion of four-fifths of all thermal anomalies. Indeed, the lic head; x (m) the abscissa along an element; A (m ) the sec- −1 hydrothermal convection along faults (activated due to the tion of the element i; S (m ) the specific storage coefficient tectonic context) may explain 75-85% of all temperature of the element i; and t (s) the time. A pipe is formed by two anomalies [21], as demonstrated in the URG context. This infinite parallel ribbons, open enough so that a flow regime method proposes practical help in locating and estimating like the flow between two parallel plates can take place. To the extension of hidden thermal anomalies, with the pros- satisfy the laminar flow regime which dominates in rock pects of reducing (1) the required number of exploration masses, the pipe conductivity is correlated to a cubic law. wells through the optimization of the exploration area and The pipe section is then proportional to the cubic root of (2) the cost of each exploration well by minimizing the dril- its conductivity. ling depth needed to reach the target temperatures. More- For a flow plane with a regularly spaced grid of channels, over, the opportunity for a better match between resources 2 the transmissivity of the fracture T (m /s) is related to the and needs offered by the 3D regional models will facilitate pipe conductivity and d (m), the spaced grid size: the deployment of deep geothermal energy in the territory and its effective integration in local and regional distribution 3 −1 C m ⋅ s networks while minimizing the environmental impacts. To Tm /s = A 2 dm conclude, the first 3D numerical flow (using the discrete- continuum method) and mechanical models (using the dis- The diffusivity equation is solved in 3D elements: tinct element method) at the URG scale offer new research opportunities. The ultimate goal of our team’s work is to pro- ∂h vide innovative exploration technologies/methodologies to div K∇h = S + q, A 3 limit the area to explore and reduce the risks associated with ∂t geothermal exploration. −1 where K (m·s ) is the permeability of the porous media; h −1 (m) the hydraulic head; and q a source term (s ). Appendix Solute transport is solved by using the random walk method. The flow is one-dimensional everywhere in pipes A. 3FLO Numerical Code and tubes, except at intersections. The dispersion is therefore only longitudinal and is “completed” by the full mixing The 3FLO code is based on the finite element method. 3FLO occurring at intersections. A pipe porosity n (-) can be spec- can solve different types of problems, such as [82] ified, as well as a pipe dispersivity (m). In a pipe, the particle displacement is rectilinear and uniform: (i) flow in fracture networks, represented by a 3D net- work of pipes. A network of channels (called pipes) uΔt is generated on each fracture. The connection of x = x + , A 4 the channels from one plane to another through the fracture plane intersections (called tubes) results −1 where x (m), u (m·s ), Δt (s), and n are the particle position in a 3D network of 1D elements along the pipe, the fluid speed, the time step, and the pipe (ii) flow (saturated or unsaturated) in porous media by porosity, respectively. The 3FLO code enables the simulation using mixed-hybrid 3D finite elements. These ele- of the particle displacement in channels and in 3D elements. Indeed, in 3D elements, once the velocity is known, the ments are more precise than classical (Galerkin) 3D elements and allow a proper computation of face movement of a particle can be separated into a longitudinal fluxes. This helps minimize solute transport compu- displacement along a flow line (corresponding to advection tation errors, as will be detailed later on. Solute and longitudinal dispersion) and two orthogonal transverse 20 Geofluids of the assembly of blocks and joints. To be solved, the state displacements (simulating transverse dispersion). Moreover, direct exchanges are possible between 1D and 3D elements; equation must be completed with the constitutive equations, particles can originate from a pipe and enter a 3D element the boundary conditions, and the initial state. and vice versa. B.2. Constitutive Equations. The mechanical constitutive equations, or laws of behaviour, give the relation between dis- B. 3DEC Numerical Code placements and stresses. For the DEM, the constitutive equa- B.1. Overview of 3DEC Software. The mechanical software tions must be provided for both the blocks and the joints. In this study, because we work in a highly segmented medium (3DEC™ for the three-dimensional Distinct Element Code, see Itasca 2016) relies on the distinct element method, where where the fault network is expected to have the greatest the contact forces between solids are explicitly incorporated. impact on the stress redistribution, the emphasis is set on the fault zone characterization. Initially proposed to study slopes in jointed rock masses [83], the DEM has been developed ever since, especially in recent B.2.1. Rock Matrix (Blocks). A simple behaviour is considered years. Today, its scope of application has extended to rock for the rock matrix, and blocks are assumed to be linearly mass scales [84–86] down to microscopic scales [87–89]. In elastic, homogeneous, and isotropic materials: addition to supporting discontinuous displacement and stress fields, the particularity of the DEM is to consider the σ − σ = C ε, B 2 full mechanical behaviours for the contacts. When the con- tacts depict fault zones, the law of behaviours embraces the where C (Pa) is the elastic stiffness tensor, σ (Pa) and σ (Pa) behaviour of the infills and of the hanging and foot walls 0 the current and initial stress tensors, respectively, and ε (−) (which belong to the rock mass), and the law of behaviours the strain tensor. A convenient way of formulating the equa- can be chosen accordingly. tion splits the strain tensor into volumetric ε (−) and devia- The DEM thus offers the possibility to V toric ε (−) parts: (i) describe discontinuous displacement and stress σ − σ = Kε δ +2Gε , B 3 fields 0 V d (ii) explicitly account for mechanically active disconti- where K (Pa) and G (Pa) are the bulk and shear moduli, nuities and their impact on blocks respectively. The volumetric strain is the trace of the strain tensor, and the deviatoric part is given by ε = ε − ε /3δ. (iii) use complex mechanical laws for the joints d V The primary unknown in the mechanics is the displacement The solution phase of the DEM is similar to that of any u (m), and its relation with the strain is yielded by the com- continuous numerical method: at each time step, the solver patibility condition: runs until mechanical equilibrium is achieved within the sys- tem. For the deformable blocks, the problem unknowns are grad u+ grad u ε = B 4 the displacements u (m) at the mesh grid points. u is com- puted by solving the balance in B.2.2. Fault Zones (Joints). A more sophisticated behaviour is considered for the fault zones. Fault zones are mechanical div σ dV + F + mg = mu, B 1 ext weaknesses in the rock mass and accommodate more defor- mation than the surrounding matrix. Because of the larger deformations involved, irreversible processes can more easily where Ω (m ) and m (kg) are the volume and mass for the take place within the fault zones and should be considered in considered block, σ (Pa) the stress tensor, F (N) the sum ext −2 the model. Recalling that joints are flat surfaces in 3DEC™, of external forces other than gravity forces, g (m·s ) the the joint behaviour can be described along the normal and gravity vector, and u (m) the displacement, and a top dot parallel directions to address opening/closing and shearing denotes a time derivative. phenomena, respectively. With the Coulomb slip law, the With the DEM, the external forces F account for the ext shear behaviour is characterized by a bilinear relation: an interactions with the contiguous blocks. These interaction elastic ramp and a plateau as soon as the onset of plasticity forces are obtained through the joint constitutive equations is reached: which, given a prescribed displacement, return the resulting forces. That is, in addition to the constitutive equations σ = yield k u , u < u , f u that must be given for the continuous methods, the s s s s τ = B 5 DEM also requires joint constitutive equations to solve equa- yield yield τ , u ≤ u , s s s tion B.1. The complete solution scheme at each time step is then achieved by balancing the displacements resulting from −1 where k (Pa·m ) is the joint shear stiffness, τ (Pa) the shear s s the block deformations with the forces resulting from the yield stress, τ (Pa) the plastic yield stress, u (m) the shear dis- block interactions (Figure 13). s yield yield At the whole-system scale (here, the URG model), the placement, and u = τ /k (m) the plastic yield displace- s s state equation solved in 3DEC™ is the mechanical balance ment. The plastic yield stress depends on the fault zone Geofluids 21 internal parameters and on the normal stress σ (Pa) acting Since the mechanical balance of the blocks explicitly on it: depends on F (see equation (B.1)), the impact of the joints ext on the blocks is directly expressed through this vector. Con- yield versely, the impact of the blocks on the joints is obtained τ = c + σ tan ϕ, B 6 s 0 n once the block mechanical balance is achieved: the balance modifies the displacement field u through the model, which where c (Pa) is the cohesion and ϕ ( ) the angle of internal yield will affect the forces acting on the joints through equations friction. When the onset of plasticity is reached u ≤ u , (B.5) and (B.7). yield the exceeding displacement u − u is permanent, hence the “plastic” denomination. When the joint lies in the plastic Data Availability phase, another mechanical feature is triggered: dilation. Dila- tion is the opening of the fault zones under shearing and is The data used to build the models and support the find- related to the fault zones actually being irregular surfaces ings of this study are included within the article or can rather than flat planes. The normal behaviour of the joint is be found within references cited (articles, database issued thus elastoplastic as well and reads from the GeORG project results: http://www.geopotenziale. org/home?lang=2). el dil σ = k u + u , B 7 n n n n Conflicts of Interest −1 el where k (Pa·m ) is the joint normal stiffness, u (m) the The authors declare that they have no conflicts of interest. n n dil elastic part of the normal displacement, and u (m) the irre- versible part of the normal displacement due to dilation. The Acknowledgments dilation relation introduces two additional material parame- This work was conducted in the framework of the IMAGE ters, the dilation angle ψ ( ) and the critical shear displace- project (Integrated Methods for Advanced Geothermal ment u (m) above which dilation vanishes: Exploration, grant agreement No. 608553). This work was completed under the REFLET and TEMPERER projects yield 0, u < u , (French National Research Agency grant agreement Nos. dil yield yield c ANR_10-IEED-0802-02 and ANR-10-IEED-0801-02, respec- u = B 8 u − u tan ψ, u ≤ u < u , n s s s s s tively) which are supported by the French Government through the Investments for the Future programmes. This 0, u ≤ u work also received internal funding from the French Geologi- cal Survey. Both dilation parameters are conditioned by the geome- try of the fault zone irregularities: ψ relates to the shape, References and u depends on the maximum height of the irregularities. The joint formulation remains valid for large displace- [1] B. Goldstein, G. Hiriart, R. Bertani et al., “Geothermal energy,” ments. The block equations however are restricted to in IPCC Special Report on Renewable Energy Sources and infinitesimal strains (see equation (B.4)), implying that Climate Change Mitigation, O. Edenhofer, R. Pichs-Madruga, any large strain at the model scale would be localized into Y. Sokona, K. Seyboth, P. Matschoss, S. Kadner, T. Zwickel, the joints. P. Eickemeier, G. Hansen, S. Schlomer, and C. von Stechow, Eds., Cambridge University Press, Cambridge, 2011. B.2.3. Block and Joint Interactions. For each contact joint, the [2] J.-D. VanMees, J. Hopman, G. Pall Hersir et al., IMAGE normal and shear stresses add up to the joint contact force F Beyond a Standardized Workflow for Geothermal Exploration, (N), which depicts the interaction force between the two blocks sharing the joint: [3] K. Young, T. Reber, and K. Witherbee, “Hydrothermal explo- ration best practices and geothermal knowledge exchange on F = σ nA + τ sA , B 9 c n c s c OpenEI,” in Proc. 37th Work. Geotherm. Reserv. Eng, pp. 1455–1469, Stanford, California, 2012. where A (m ) is the contact area and n (−) and s (−) the unit c [4] L. Rybach, “The geothermal conditions in the Rhine graben - a normal vectors in the normal and shear directions, respec- summary,” Bulletin Für Angewandte Geologie, vol. 12, pp. 29– tively, along the joint. For each block, the external force vec- 32, 2007. tor F is expressed as the sum of the contact forces due to [5] D. L. Siler, J. E. Faulds, B. Mayhew, and D. D. McNamara, ext the n joints surrounding the block: “Analysis of the favorability for geothermal fluid flow in 3D: Astor Pass geothermal prospect, Great Basin, northwestern nc Nevada, USA,” Geothermics, vol. 60, pp. 1–12, 2016. F = 〠 F , B 10 ext c [6] M. W. Swyer, T. T. Cladouhos, C. Forson, J. L. Czajkowski, i=1 N. C. Davatzes, and G. M. Schmalzle, “Permeability potential modeling of geothermal prospects combining regional crustal where F (N) is the contact force of the i-th joint surrounding strain rates with geomechanical simulation of fault slip and the block. volcanic center deformation: a case study for Washington 22 Geofluids [21] P. Baillieux, E. Schill, J.-B. Edel, and G. Mauri, “Localiza- State geothermal play fairways,” in 50th US Rock Mechanics/- Geomechanics Symposium 2016, Houston, Texas, 2016. tion of temperature anomalies in the Upper Rhine Graben: insights from geophysics and neotectonic activity,” Interna- [7] J.-B. Edel, K. Schulmann, and Y. Rotstein, “The Variscan tional Geology Review, vol. 55, no. 14, pp. 1744–1762, tectonic inheritance of the Upper Rhine Graben: evidence of reactivations in the Lias, Late Eocene–Oligocene up to the recent,” International Journal of Earth Sciences, vol. 96, no. 2, [22] H. Auradou, G. Drazer, A. Boschan, J. P. Hulin, and J. Koplik, pp. 305–325, 2007. “Flow channeling in a single fracture induced by shear displacement,” Geothermics, vol. 35, no. 5-6, pp. 576–588, [8] M. E. Schumacher, “Upper Rhine Graben: role of preexisting structures during rift evolution,” Tectonics, vol. 21, no. 1, pp. 6-1–617, 2002. [23] S. Gentier, E. Lamontagne, G. Archambault, and J. Riss, “Anisotropy of flow in a fracture undergoing shear and [9] P. Baillieux, E. Schill, Y. Abdelfettah, and C. Dezayes, “Possible its relationship to the direction of shearing and injection natural fluid pathways from gravity pseudo-tomography in the pressure,” International Journal of Rock Mechanics and geothermal fields of Northern Alsace (Upper Rhine Graben),” Mining Sciences, vol. 34, no. 3-4, pp. 94.e1–94.e12, Geothermal Energy, vol. 2, no. 1, p. 16, 2014. [10] C. Clauser and H. Villinger, “Analysis of conductive and [24] J. V. Rowland and R. H. Sibson, “Structural controls on convective heat transfer in a sedimentary basin, demonstrated hydrothermal flow in a segmented rift system, Taupo for the Rhein graben,” Geophysical Journal International, Volcanic Zone, New Zealand,” Geofluids, vol. 4, no. 4, 283 vol. 100, no. 3, pp. 393–414, 1990. pages, 2004. [11] S. Gentier, X. Rachez, T. Dung et al., 3D Flow Modelling [25] D. L. Siler, N. H. Hinz, J. E. Faulds, and J. Queen, “3D analysis of the Medium-Term Circulation Test Performed in the of geothermal fluid flow favorability: Brady’s, Nevada, USA,” Deep Geothermal Site of, World Geothermal Congress, in The 41st Workshop on Geothermal Reservoir Engineering, Stanford University, Standford, California, 2016. [12] L. Guillou-Frottier, C. Carrė, B. Bourgine, V. Bouchot, and [26] D. E. Dempsey, S. F. Simmons, R. A. Archer, and J. V. A. Genter, “Structure of hydrothermal convection in the Rowland, “Delineation of catchment zones of geothermal Upper Rhine Graben as inferred from corrected tempera- systems in large-scale rifted settings,” Journal of Geophysical ture data and basin-scale numerical models,” Journal of Research: Solid Earth, vol. 117, no. B10, 2012. Volcanology and Geothermal Research, vol. 256, pp. 29–49, [27] Y. Rotstein, J.-B. Edel, G. Gabriel, D. Boulanger, M. Schaming, and M. Munschy, “Insight into the structure of the Upper [13] C. Le Carlier, J.-J. Royer, and E. L. Flores, “Convetive Rhine Graben and its basement from a new compilation of heat transfer at the Soultz-sous-Forets Geothermal Site: Bouguer Gravity,” Tectonophysics, vol. 425, no. 1-4, pp. 55– implications for oil potential,” First Break, vol. 12, 70, 2006. no. 1285, 1994. [28] S. Gentier, X. Rachez, M. Peter-Borie, and M. Loubaud, “A [14] B. Sanjuan, R. Millot, C. Innocent, C. Dezayes, J. Scheiber, and flow model of the deep geothermal reservoir of Soultz-sous- M. Brach, “Major geochemical characteristics of geothermal Forêts (France),” in 47 US Rock Mechanics/Geomechanics brines from the Upper Rhine Graben granitic basement with Symposium, San Francisco, 2013. constraints on temperature and circulation,” Chemical Geol- ogy, vol. 428, pp. 27–47, 2016. [29] B. C. Valley, “The relation between natural fracturing and stress heterogeneities in deep-seated crystalline rocks at [15] I. Spottke, E. Zechner, and P. Huggenberger, “The southeast- Soultz-sous-Forêts (France),” in ETH Reprozentrale Höngger- ern border of the Upper Rhine Graben: a 3D geological model berg, HIL C45, Zürich (2007), Swiss Federal Institute of Tech- and its importance for tectonics and groundwater flow,” Inter- national nology Zurich, 2007. Journal of Earth Sciences, vol. 94, no. 4, pp. 580–593, [30] T. Plenefisch and K.-P. Bonjer, “The stress field in the Rhine Graben area inferred from earthquake focal mechanisms and [16] D. King and E. Metcalfe, “Rift zones as a case study for advancing geothermal occurence models,” PROCEEDINGS: estimation of frictional parameters,” Tectonophysics, vol. 275, no. 1-3, pp. 71–97, 1997. Thirty-Eighth Workshop on Geothermal Reservoir Engineer- ing, Standford University, p. 11, 2013. [31] J. B. Witter and G. Melosh, “The value and limitations of 3D models for geothermal exploration,” in 43rd [17] J. Freymark, J. Sippel, M. Scheck-wenderoth et al., The Thermal Field of the Upper Rhine Graben–Temperature Pre- Workshop on Geothermal Reservoir Engineering, no. article SGP-TR-213, 2018Standford Univeristy, Standford, California, dictions Based on a 3D Model, Eur. Geotherm. Congr. 2016, 2016. 2018. [18] P. Baillieux, Multidisciplinary Approach to Understand the [32] D. D. Pollard and P. Segall, “Theoretical displacements and stresses near fractures in rock: with applications to faults, Localization of Geothermal Anomalies in the Upper Rhine Gra- ben from Regional to Local Scale, Fac. Sci. Hydrogeol. Geotherm, joints, veins, dikes, and solution surfaces,” in Fracture Mechanics of Rock, pp. 277–349, Academic Press, San Diego, [19] J. Tóth, “Geothermal phenomena in the context of gravity- driven basinal flow of groundwater,” Central European Geol- [33] C.-F. Tsang and I. Neretnieks, “Flow channeling in heteroge- ogy, vol. 58, no. 1-2, pp. 1–27, 2015. neous fractured rocks,” Reviews of Geophysics, vol. 36, no. 2, pp. 275–298, 1998. [20] J. Vidal and A. Genter, “Overview of naturally permeable fractured reservoirs in the central and southern Upper Rhine [34] A. Herbert, Modelling Approaches for Discrete Fracture Graben: insights from geothermal wells,” Geothermics, Network Flow Analysis, Coupled Thermo-Hydro-Mechanical vol. 74, pp. 57–73, 2018. Processes of Fractured Media, Elsevier, Amsterdam, 1996. Geofluids 23 [52] Y. Zhou and W. Li, “A review of regional groundwater flow [35] Itasca Consulting Group Inc, 3DEC–Three-Dimensional Distinct Element Code, Ver. 5.2, Itasca Consulting Group modeling,” Geoscience Frontiers, vol. 2, no. 2, pp. 205–214, Inc., Minneapolis, Minnesota, 2017. 2011. [36] GeORG project, 2012, [WWW Document]. [53] F. Rummel, “Physical properties of the rock in the granitic sec- tion of borehole GPK1 borehole, Soultz-sous-Forêts,” in Geo- [37] G. B. Baecher, “Statistical analysis of rock mass fracturing,” thermal Energy in Europe: The Soultz Hot Dry Rock Project, Journal of the International Association for Mathematical pp. 199–216, Gordon & Breach Science Publishers Ltd, 1992. Geology, vol. 15, no. 2, pp. 329–348, 1983. [54] I. Stober and K. Bucher, “Hydraulic and hydrochemical [38] P. H. S. W. Kulatilake and T. H. Wu, “Estimation of mean trace properties of deep sedimentary reservoirs of the Upper Rhine length of discontinuities,” Rock Mechanics and Rock Engineer- Graben, Europe,” Geofluids, vol. 15, no. 3, 482 pages, 2015. ing, vol. 17, no. 4, pp. 215–232, 1984. [55] T. Gleeson, L. Marklund, L. Smith, and A. H. Manning, [39] S. D. Priest and J. A. Hudson, “Estimation of discontinuity “Classifying the water table at regional to continental scales,” spacing and trace length using scanline surveys,” International Geophysical Research Letters, vol. 38, no. 5, 2011. Journal of Rock Mechanics and Mining Science and Geomecha- nics Abstracts, vol. 18, no. 3, pp. 183–197, 1981. [56] H. M. Haitjema and S. Mitchell-Bruker, “Are water tables a subdued replica of the topography?,” Ground Water, vol. 43, [40] L. Zhang and H. H. Einstein, “Estimating the intensity of rock no. 6, pp. 781–786, 2005. discontinuities,” International Journal of Rock Mechanics and Mining Sciences, vol. 37, no. 5, pp. 819–837, 2000. [57] D. P. Yale, “Fault and stress magnitude controls on variations in the orientation of in situ stress,” Geological Society, London, [41] W. Riedel, “Zur Mechanik Geologischer Brucherscheinun- Special Publications, vol. 209, no. 1, pp. 55–64, 2003. gen,” Zur Mechanik Geologischer Brucherscheinungen B, pp. 354–368, 1929. [58] L. Dorbath, K. Evans, N. Cuenot, B. Valley, J. Charléty, and M. Frogneux, “The stress field at Soultz-sous-Forêts from focal [42] J. H. Behrmann, O. Hermann, M. Horstmann, D. C. Tanner, mechanisms of induced seismic events: cases of the wells and G. Bertrand, “Anatomy and kinematics of oblique conti- GPK2 and GPK3,” Comptes Rendus Geoscience, vol. 342, nental rifting revealed: a three-dimensional case study of the no. 7-8, pp. 600–606, 2010. southeast Upper Rhine Graben (Germany),” American Associ- ation of Petroleum Geologists Bulletin, vol. 87, no. 7, pp. 1105– [59] T. Hettkamp, J. Baumgärtner, R. Baria et al., Electricity Produc- 1121, 2003. tion from Hot Rocks, in: Proceedings, Twenty-Ninth Workshop on Geothermal Reservoir Engineering, Stanford University, [43] N. Harthill, No Title. Tecton. Basis Geotherm. Reg. Oberrhein- Stanford, California, 2018. graben, 20 Jahre Tiefe Geotherm. Deutschland. 7. Geotherm. Fachtagung, Waren, 2002. [60] Y. Guglielmi, D. Elsworth, F. Cappa et al., “In situ observations on the coupling between hydraulic diffusivity and displace- [44] P. A. Domenico and V. V. Palciauskas, “Theoretical analysis of ments during fault reactivation in shales,” Journal of Geophys- forced convective heat transfer in regional ground-water flow,” ical Research: Solid Earth, vol. 120, no. 11, pp. 7729–7748, Geological Society of America Bulletin, vol. 84, no. 12, p. 3803, [61] X. Rachez and S. Gentier, 3D-Hydromechanical Behavior of a [45] M. P. Anderson, “Heat as a ground water tracer,” Ground Stimulated Fractured Rock Mass, vol. 8, World Geotherm. Water, vol. 43, no. 6, pp. 951–968, 2005. Congr. 2010, 2010. [46] M. O. Saar, “Review: geothermal heat as a tracer of large-scale [62] T. J. Buchmann and P. T. Connolly, “Contemporary kinemat- groundwater flow and as a means to determine permeability ics of the Upper Rhine Graben: a 3D finite element approach,” fields,” Hydrogeology Journal, vol. 19, no. 1, pp. 31–52, 2011. Global and Planetary Change, vol. 58, no. 1-4, pp. 287–309, [47] J. Faulds, M. Coolbaugh, V. Bouchot, I. Moeck, O. Kerem, and O. Cedex, Characterizing Structural Controls of Geothermal [63] A. Hosni, Modélisation par la Méthode des Éléments Dis- Reservoirs in the Great Basin, USA, and Western Turkey: tincts du champ des contraintes à l’échelle du Fossé Rhénan Developing Successful Exploration Strategies in Extended Ter- et de Soultz-sous-Forêts, Institu National Polytechnique de ranes, pp. 25–29, 2010. Lorraine, 1997. [48] J. E. Faulds, N. H. Hinz, M. F. Coolbaugh et al., “Assessment of [64] G. Peters, Active Tectonics in the Upper Rhine Graben, p. 298, favorable structural settings of geothermal systems in the Great Basin, Western USA,” Geothermal Resources Council Transactions, vol. 35, pp. 777–783, 2011. [65] Y. Gunzburger and V. Magnenet, “Stress inversion and basement-cover stress transmission across weak layers in the [49] N. H. Hinz, J. E. Faulds, and D. L. Siler, “Developing systematic Paris basin, France,” Tectonophysics, vol. 617, pp. 44–57, 2014. workflow from field work to quantitative 3D modeling for successful exploration of structurally controlled geothermal [66] T. Hergert, O. Heidbach, K. Reiter, S. B. Giger, and systems,” Geothermal Resources Council Transactions, P. Marschall, “Stress field sensitivity analysis in a sedimentary vol. 37, pp. 275–279, 2013. sequence of the Alpine foreland, northern Switzerland,” Solid Earth, vol. 6, no. 2, pp. 533–552, 2015. [50] I. C. Wallis, D. McNamara, J. V. Rowland, and C. Massiot, “The nature of fracture permeability in the basement grey- [67] K. Reiter and O. Heidbach, “3-D geomechanical-numerical wacke at Kawerau Geothermal Field, New Zealand,” in Pro- model of the contemporary crustal stress state in the Alberta ceedings Thirty-seventh Workshop on Geothermal Reservoir Basin (Canada),” Solid Earth, vol. 5, no. 2, pp. 1123–1149, Engineering, Stanford University, California, 2012. 2014. [51] R. A. Freeze and P. A. Witherspoon, “Theoretical analysis of [68] P. Connolly and J. Cosgrove, “Prediction of static and dynamic regional groundwater flow: 1. Analytical and numerical solu- fluid pathways within and around dilational jogs,” Geological tions to the mathematical model,” Water Resources Research, Society, London, Special Publications, vol. 155, no. 1, pp. 105– vol. 2, no. 4, pp. 641–656, 1966. 121, 1999. 24 Geofluids [84] V. Bonilla-Sierra, L. Scholtès, F. V. Donzé, and M. K. Elmouttie, [69] D. Curewitz and J. A. Karson, “Structural settings of hydro- thermal outflow: fracture permeability maintained by fault “Rock slope stability analysis using photogrammetric data and propagation and interaction,” Journal of Volcanology and Geo- DFN–DEM modelling,” Acta Geotechnica, vol. 10, no. 4, thermal Research, vol. 79, no. 3-4, pp. 149–168, 1997. pp. 497–511, 2015. [70] M. Cathelineau and M.-C. Boiron, “Downward penetration [85] K. Farahmand, I. Vazaios, M. S. Diederichs, and N. Vlachopoulos, “Investigating the scale-dependency of the and mixing of sedimentary brines and dilute hot waters at 5km depth in the granite basement at Soultz-sous-Forêts geometrical and mechanical properties of a moderately jointed (Rhine graben, France),” Comptes Rendus Geoscience, vol. 342, rock using a synthetic rock mass (SRM) approach,” Computers no. 7-8, pp. 560–565, 2010. and Geotechnics, vol. 95, pp. 162–179, 2018. [71] M. Cacace and A. B. Jacquey, “Flexible parallel implicit model- [86] J. H. Wu, W. K. Lin, and H. T. Hu, “Assessing the impacts of a ling of coupled thermal-hydraulic-mechanical processes in large slope failure using 3DEC: the Chiu-fen-erh-shan residual fractured rocks,” Solid Earth, vol. 8, no. 5, pp. 921–941, 2017. slope,” Computers and Geotechnics, vol. 88, pp. 32–45, 2017. [72] Z. Q. Huang, P. H. Winterfeld, Y. Xiong, Y. S. Wu, and [87] E. Ghazvinian, M. S. Diederichs, and R. Quey, “3D random J. Yao, “Parallel simulation of fully-coupled thermal-hydro- Voronoi grain-based models for simulation of brittle rock mechanical processes in CO2 leakage through fluid-driven damage and fabric-guided micro-fracturing,” Journal of fracture zones,” International Journal of Greenhouse Gas Rock Mechanics and Geotechnical Engineering, vol. 6, no. 6, Control, vol. 34, pp. 39–51, 2015. pp. 506–521, 2014. [73] A. Lisjak, O. K. Mahabadi, L. He et al., “Acceleration of a [88] N. P. Kruyt and L. Rothenburg, “A micromechanical study of 2D/3D finite-discrete element code for geomechanical simula- dilatancy of granular materials,” Journal of the Mechanics tions using general purpose GPU computing,” Computers and and Physics of Solids, vol. 95, pp. 411–427, 2016. Geotechnics, vol. 100, pp. 84–96, 2018. [89] W. Li, Y. Han, T. Wang, and J. Ma, “DEM micromechanical [74] S. Salimzadeh, A. Paluszny, H. M. Nick, and R. W. Zimmerman, modeling and laboratory experiment on creep behavior of salt “A three-dimensional coupled thermo-hydro-mechanical rock,” Journal of Natural Gas Science and Engineering, vol. 46, model for deformable fractured geothermal systems,” Geother- pp. 38–46, 2017. mics, vol. 71, pp. 212–224, 2018. [90] T. G. Farr and M. Kobrick, “Shuttle Radar Topography Mis- [75] B. Vallier, V. Magnenet, J. Schmittbuhl, and C. Fond, “THM sion produces a wealth of data,” Eos, Transactions American modeling of hydrothermal circulation at Rittershoffen geo- Geophysical Union, vol. 81, no. 48, pp. 583–583, 2000. thermal site, France,” Geothermal Energy, vol. 6, no. 1, 2018. [76] M. O. Ziegler, O. Heidbach, J. Reinecker, A. M. Przybycin, and M. Scheck-Wenderoth, “A multi-stage 3-D stress field model- ling approach exemplified in the Bavarian Molasse Basin,” Solid Earth, vol. 7, no. 5, pp. 1365–1382, 2016. [77] J. P. Brun, M.-A. Gutscher, and dekorp-ecors teams, “Deep crustal structure of the Rhine Graben from dekorp-ecors seis- mic reflection data: a summary,” Tectonophysics, vol. 208, no. 1-3, pp. 139–147, 1992. [78] M. Schwarz and A. Henk, “Evolution and structure of the Upper Rhine Graben: insights from three-dimensional thermomechanical modelling,” International Journal of Earth Sciences, vol. 94, no. 4, pp. 732–750, 2005. [79] F. Wenzel, J. P. Brun, and ECORS‐DEKORP team, “Crustal scale structure of the Rhine Graben from ECORS‐DEKORP deep seismic reflection data,” in SEG Technical Program Expanded Abstracts 1991, pp. 163–166, Houston, Texas, 1991. [80] J. R. Harris, E. Grunsky, P. Behnia, and D. Corrigan, “Data- and knowledge-driven mineral prospectivity maps for Canada’s North,” Ore Geology Reviews, vol. 71, pp. 788–803, 2015. [81] J. Iovenitti, F. H. Ibser, M. Clyne, J. Sainsbury, and O. Callahan, “The Basin and Range Dixie Valley Geothermal Wellfield, Nevada, USA—a test bed for developing an Enhanced Geothermal System exploration favorability meth- odology,” Geothermics, vol. 63, pp. 195–209, 2016. [82] B. Paris, Investigation of the Correlation between Early-Time Hydraulic Response and Tracer Breakthrough Times in Fractured Media, Äspö Hard Rock Laboratory, 2002. [83] P. A. Cundall, “A computer model for simulating progres- sive large scale movements in blocky rock systems,” in Proc. Symp. Int. Soc. rock Mech. (Nancy, Fr. 1971) 1, Nancy, France, 1971. International Journal of The Scientific Advances in Advances in Geophysics Chemistry Scientica World Journal Public Health Hindawi Hindawi Hindawi Hindawi Publishing Corporation Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 http://www www.hindawi.com .hindawi.com V Volume 2018 olume 2013 www.hindawi.com Volume 2018 Journal of Environmental and Public Health Advances in Meteorology Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 Submit your manuscripts at www.hindawi.com Applied & Environmental Journal of Soil Science Geological Research Hindawi Volume 2018 Hindawi www.hindawi.com www.hindawi.com Volume 2018 International Journal of International Journal of Agronomy Ecology International Journal of Advances in International Journal of Forestry Research Microbiology Agriculture Hindawi Hindawi Hindawi Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 International Journal of Journal of Journal of International Journal of Biodiversity Archaea Analytical Chemistry Chemistry Marine Biology Hindawi Hindawi Hindawi Hindawi Hindawi www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018 www.hindawi.com Volume 2018

Journal

GeofluidsHindawi Publishing Corporation

Published: May 12, 2019

There are no references for this article.