Abstract
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’s efforts 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
Worldwide, deep geothermal energy offers an enormous potential for future electricity and heat productions. Geothermal energy is sustainable, clean, and renewable, as the tapped heat from an active reservoir is continuously restored [1]. Geothermal energy has a number of positive characteristics: it is environmentally friendly, has low emissions (greenhouse gases), and uses a low number of raw materials such as rare elements and strategic metals. In addition, unlike many other renewable energies depending on climate conditions (wind, solar, hydraulic regimes, and the water levels of rivers), the geothermal energy from the Earth’s crust is nonintermittent. However, only a fraction of the deep geothermal energy of the Earth’s crust is currently exploited. In practice, the exploitation is mostly limited to “brown fields,” i.e., areas where a reservoir is already sited or where there are targeted structures and parameters characterized that are already known [2]. The wider development of deep geothermal resource exploitation now requires focusing on little to unknown areas, called “green fields.” In these fields, geothermal resource exploration (including exploratory drillings) involves a high degree of uncertainty and financial risk and requires reliable data to provide convincing arguments in the decision-making process [3]. Consequently, the improvement of exploration methods is required to lower the risk and to overcome the barriers to the wider development of deep geothermal projects.
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 the stresses and fluid circulations that impact heat-driving processes [4–6]. The Upper Rhine Graben (URG) is a relatively well-known example of a favourable faulted geothermal system showing thermal anomalies at the regional scale. It has been largely investigated, providing various data such as a tectonic history of the Rhine Graben [7, 8], in situ observations from deep boreholes, and numerical models and geochemical studies [9–15]. Moreover, a European project (the GeORG project) has resulted in an important structural database (http://www.geopotenziale.org/home?lang=3) that provides a map of seismic faults at the URG scale. The formation of geothermal resources in the URG results from the interactions between different coupled processes comprising mechanical constraints, groundwater flows, and heat transfer. Large-scale structures of the graben and its associated faults may act as a primary structural control on the upwelling geothermal fluids [16]. The rifting context and the associated thinned crust and elevated mantle favour higher temperatures at shallower depths [4, 17]. A series of thermal anomalies with temperatures above 140°C at a 2 km depth characterizes the URG. These anomalies are mainly located in the western part of the URG [18]. Recalling that thermal anomalies are mostly affected by the regional groundwater circulation [19], the high potential of the URG for thermal anomalies is widely interpreted as a signature of the regional-scale fluid migrations hosted by the numerous fault zones cutting through the rock mass (multiscale fracture-controlled systems, [20]). Arguing that the fault zones are critically stressed and mostly under shearing conditions, Baillieux et al. [21] estimated that the regional fluid migrations could be responsible for up to 75-85% of the thermal anomalies in the URG. The radiogenic heat production due to the crystalline composition of the basement may explain the remaining 15-25%. The geochemical data tend to confirm this assumption: Sanjuan et al. [14] conclude that the thermal anomaly in the Soultz-sous-Forêts surroundings results from the regional circulation through a complex, but still poorly defined, system of deep faults (probably along NE-SW fault zones). Thus, the deep circulation is probably the origin of the thermal anomalies and governs their geographical distributions. Their understanding is therefore fundamental for the widespread deployment of deep geothermal energy to other areas, for instance, the Limagne Graben.
Several researchers have suggested that the stresses associated with faults significantly affect the flow properties and that there is a strong correlation between the groundwater flow and geothermal anomalies. A regional heat-flow study requires an analysis of the basin-scale flow distribution (due to the strong link between the regional groundwater flow and anomalies of the geothermal heat) and should account for the active faulting caused by stresses. Stresses can play an important role in the setup of thermal anomalies by favouring preferential flow paths through the rock mass [6]. For example, shear stresses on fractures have been observed to favour the fluid circulation perpendicular to the shear direction [22, 23]. Then, stresses can favour the appearance of thermal anomalies by enabling vertical upflows of a fluid previously heated during deep circulation [6, 24]. Thus, key factors to understanding thermal anomalies are 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 exploitable hot fluid upwelling. That is the reason why geothermal studies of an extended area, such as exploration for geothermal resources, must consider these key factors and must be based on the study of these physics.
In addition to other characterization methods (geological, geophysical, and geochemical), exploration could benefit from numerical models and then a predictive analysis to locate geothermal resources. Numerical models could allow an integrated evaluation of the physical processes in the URG with adapted time and space scales, considering 3D effects. They constitute a useful tool to improve the understanding of the dynamic context and may provide answers to the challenging task of geothermal exploration. Indeed, the physical conditions favouring the setup of the geothermal resources can be searched for in predictive models, hence giving estimates of the so-called “favourable areas.” Some major phenomena were highlighted, and some local anomalies were explained using 1D hydrothermal models and 2D models based on continuum approaches [10, 12, 13] but without any hope for an industrial application to exploration. These models are indeed restricted to local scales and do not consider major features such as preferential pathways within fault zones or 3D effects. Stress and numerical flow models have also been used to study geothermal systems ([6, 26] and references therein), but either the scale is not appropriate or the approach is not holistic. Thus, to the authors’ knowledge, there is currently no numerical approach available to contribute to geothermal resource exploration.
The present work aims at demonstrating the potential of the combined use of 3D flow and stress models at the regional scale (URG scale) to locate/understand the location of hidden thermal anomalies. Considering the large amount of data and knowledge, we use the URG as an application case. The results of the GeORG project constitute a relevant basis to build a structural model, which is essential to developing the subsequent numerical physical models and then such a numerical approach. This enables us to compare results with observed physical evidence. The main objective of this paper is to investigate the potential of an approach based on the use of numerical models (stress and flow models) to locate hidden thermal anomalies. The paper is organized as follows: first, we give a brief overview of the study area. Then, we present the 3D geometry for the models according to the faulted geological setting. The major faults of the URG are considered using a DFN model (detailed in Section 3) consistent with the geological knowledge of the area (see Section 2). This DFN (discrete fracture network) model is declined in 3D flow and stress models at the URG scale, which are, respectively, detailed in Sections 4 and 5. The results of the physical models are interpreted as functions of their thermal signatures through the definition of geothermal criteria related to the physics considered and the setting up of the geothermal resources. These criteria enable the elaboration of indicators to use the results of the simulation and then identify geothermally favourable areas in relation with each physic. Then, considering the strong link between the stress, fluid flow, and geothermal resources, a cross-analysis of the results is realized to delineate the favourable areas for geothermal resources. The results are compared with existing thermal data at the URG scale and with knowledge gained through numerous studies (see Section 6). Before concluding, we discuss the potential of the present work.
2. Geological Setting of the Study Area: the URG
The URG is an NNE-trending crustal-scale rift of the European Cenozoic Rifting System (ECRIS) [8]. Its length reaches approximately 300 km from Frankfurt (Germany) to Basel (Switzerland), and its width is 30-40 km (Figure 1). The URG development started from the late Eocene, mainly by the reactivation of the late Variscan, Permo-Carboniferous, 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 aspects are recalled here.
Since the Visean age during which the granitic intrusion took place, a repeatedly changing stress field has occurred in the URG, leading to the reactivation of a complex set of crustal discontinuities and the creation of a complex fault network in the sedimentary blanket [8]. Based on the tectonic history and on knowledge acquired from previous works [28, 29], four main fault sets corresponding to the main statistical models have been individualized: (i)The F1 set, mainly ~N-S striking and W-dipping, is widely observed in the Rhine Graben and commonly interpreted as the result of the Oligocene deformation (N0/20° striking, W-dipping). However, this set could also be related to (a)the lower Carboniferous to Permian N20° sinistral faulting reactivated during the late Eocene [8] or(b)N170° sinistral strike-slip faults reactivated from the late Miocene up to now under a NW-SE compression of the graben as a strike-slip system(ii)The F2 set, ~NNW-SSE striking and E-dipping, is also commonly interpreted as the result of the Oligocene deformation (N0/20° striking, E-dipping). However, it could also be considered as different sets gathering the Hercynian N20° sinistral strike-slip and the N135° dextral strike-slip reactivated during the late Eocene and to N170° sinistral strike-slip faults(iii)The F3 and F4 sets, respectively, NE-SW and NW-SE striking, are consistent with the Hercynian orientations observed regionally(iv)A fifth fault set (F5) is introduced to incorporate the “background noise”
Considering the stress regime, the sinistral strike-slip regime that started at the onset of the Miocene has been persisting until today, and the present-day stress state is assumed to be driven by the African-Eurasian plate convergence [30].
As a consequence of its polyphaser tectonic history, the URG is a complex environment setting in which the rifting process played a major role by providing the space for sedimentation, creating a series of basins, and providing a network of faults along the active rift margins that enabled fluid circulation and conditioned water exchanges between basins [8, 14]. The transnational geological model of the URG based on boreholes and new processed seismic data constitutes an important piece of the puzzle for the numerical models by providing essential data. The map of the faults at the top of the basement can be used to develop numerical models, and the in situ temperature data at the URG scale constitute a relevant means to assess the fault potential (Figure 1).
3. DFN Model at the URG Scale
3.1. Why a Discrete Approach
The very first step of our work is to build a model geometry consistent with the geology and the physical processes involved. The large heterogeneities introduced in the physical system by the presence of the fault zones rule out any chance for an analytical solution, thus making numerical tools necessary. As argued by Witter and Melosh [31], the choice of the numerical tools is a challenge in itself since it will impact the results. The need to account for real 3D effects raises another challenge, which is the available panel of suitable tools: the numerical tools for solving large-scale models in reasonable computation times while offering the possibility to take into account geological structures are few. In the case of faulted settings, the presence of natural discontinuities that can result in heterogeneous stress fields and localized fluid flow pathways [32, 33] has promoted the development of fracture network models for the numerical simulation of fractured rocks [34]. Here, the key role of the regional fault network is taken into account due to a discrete numerical approach (DFN) for fault zones at the regional scale that remains a challenge alone. The model to represent the regional fault network is generated with the 3DEC™ software [35].
3.2. Geometry Building
The case study is a km2 area centred on the main thermal anomalies in the URG (Figure 1). The lateral east-west boundaries are chosen large enough to incorporate a large portion of the rift shoulders to limit the boundary effects in the subsequent computations (at least in this direction). The geometry building is focused on the conceptualization of the 3D fault zone network into a DFN at the studied scale; more specifically, the primary source of the data is the seismic fault map of the transnational database [36], where any valuable data (e.g., seismic reflection results, borehole data, outcrop traces, and scanline or window sampling) should be integrated into the structural interpretation to provide estimates on fault zone orientations, dip angles, and spatial extents [37–40]. Additional generic geological knowledge is required to complete the missing information (e.g., the dips of some deep fault zones) and to guide the conceptualization into a DFN (e.g., the fault zone priority rule dictating “who stops on who”). This includes the following nonexhaustive list of items: structural expertise on geological settings similar to the one studied (in our case, rifting systems), an understanding of the tectonic history, and the integration of generic structural tools (such as Riedel systems, see [41]). Furthermore, the numerical constraints and specificities of the codes used have to be taken into account.
Following the above-defined concept and using the geological knowledge of the URG, the discontinuities of the DFN model are the geological faults observed at different scales, from plurikilometric to kilometric scales. The largest discontinuities that cross the model are on the regional scale (primary faults of hundreds of kilometres known as major faults of the Variscan or tertiary rifting history). As discontinuities split the model blocks, the resulting block sizes decrease, allowing the integration of more local-scale discontinuities (faults of several kilometres from seismic profiles and lineament maps). The main faults integrated in the DFN model are the faults delimiting the URG and those inherited from the main Variscan NE-SW trending sinistral shearing (Lubine and Baden-Baden faults) ([7, 8, 42], among numerous others), as shown in Figure 2(a). The secondary fault network within the basement of the URG has been built from data of the GeORG project (within the basement and sedimentary layers, see [36]) or from data found in the literature [7, 27] (Figure 2(b)). The secondary fault network in the basement on both sides of the URG is from geological maps and digital terrain model (DTM) analysis (Figure 2(c)).
(a)
(b)
(c)
Note that the fault network is denser in the URG than on 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 the basement (estimated using the fault offset). It should be noted that additional geometric simplifications have been applied to fulfil the 3DEC™ topological requirements: (i)Planar surfaces are considered in 3DEC™, so the initial curved fault zones were flattened(ii)3DEC™ does not accept discontinuity tips within blocks, so the fault zones were extended until they stopped on a boundary or a previously created fault zone
The fault zone extension requires the definition of a priority rule to assess which fault zone stops on which one. The priority rule was defined according to the tectonic history first and then using a Riedel system for the contemporary fault zones ([43] evidenced a Riedel system for the Variscan fault zones).
Following these main steps based on the use of the geological structures and structural expertise, the final geometry of the DFN is composed of 351 fault zones in the region of interest (generated with the 3DEC™ software, Figure 3). The same 3D geometry will be used for the 3D flow and stress models.
4. 3D Flow Model at the URG Scale
Regionally moving groundwater is the basic and common cause of a wide variety of natural processes and phenomena that makes the gravity-driven groundwater flow a key process [19, 44]. The key reason for the geological scale effects of the groundwater flow is that the flow systems are ubiquitous and active over large spectra of time and space. Among the numerous processes associated with regional groundwater flow systems, the transport of heat is the most visible and best understood [45, 46]. Many geothermal systems are associated with high densities of faults and fractures [25, 47–50]. Geological discontinuities such as faults and stratigraphic layers play a key role in controlling the groundwater flow within faulted systems [51].
The numerical simulation of groundwater flows has become a common hydrogeological tool to investigate a variety of geological settings. Significant advances in regional flow system analysis have been made through the application of 3D groundwater flow models [52]. Considering that regional groundwater flows enable the approximation of the anomalies of the geothermal heat, we used a flow model to locate the preferential discharge areas. The flow model is built with the conviction that the patterns of basin-scale groundwater flow are mainly influenced by (1) the relief of the water table and (2) permeability heterogeneities of the basin’s rock framework (more specifically, permeability contrasts associated with faults). With that in mind, the geometry of the 3D fault network and its hydraulic connections with 3D elements (sedimentary units) need to be treated with special care. The modelling of a regional groundwater system requires large sets of data, stored in various forms and scales (maps, tables, spreadsheets, etc.). Thus, a GIS (geographical information system) software was used to improve the pre- and postprocessing of the data (ArcGIS). In the following subsections, we describe the main steps to build the flow model at the URG scale and the simulations realized to simulate the flow paths related to geothermal resources. The 3D flow model was built using the 3FLO software (Itasca®), which is designed to solve groundwater flow and transport equations in porous or fractured media (see the description of the code in Appendix A: 3FLO Numerical Code).
4.1. 3D Geometry: DFN and 3D Volumetric Elements
The starting point to build the 3D DFN flow model is the 3D regional fault network of the URG developed in Section 3. We used information provided by the 3DEC™ model (relevant data for each fault plane: strike, dip, coordinates of the mass centre, and intersections with other planes) to reproduce the same 3D network in 3FLO. In 3FLO (Appendix A: 3FLO Numerical Code), the flow within faults is supported by 1D elements (called pipes). These 1D elements are used to build a 2D grid, using two sets of pipes that are each regularly spaced. Both sets are oriented perpendicularly, with orientations of 0 and 90° (with respect to the plane dip), respectively, forming a grid, which is generated on each flow plane. According to the 3DEC structural model, the whole model domain of the URG for the flow model is a parallelepiped measuring 130 km in width (west to east), 150 km in length (north to south), and approximately 10 km in height (depending on the surface topography).
Based on an analysis of the fault sets (F1 to F5) in line with the knowledge of the regional tectonic history (Section 3.2 and Figure 4(d)), each fault has been assigned to a fault set as a function of its orientation.
Sedimentary deposits are modelled by 3D elements. We used the results provided by a transnational database [36] to build a geometric model of the basin. The thickness of the sedimentary layers ranges between 2 km and 6 km (in the eastern part), with an average value of 4 km. As a result, a 3D block domain (composed of tetrahedral elements) is built from digital elevation models (DEMs) of the basement top surface and the surface topography of the domain area. Note that the 3D block domain is discretized into more than 10,000 elements (Figure 4(e)).
4.2. Hydraulic Properties
The 3D fault network (Figure 4(d)) and 3D block domain (sedimentary cover, Figure 4(e)) have to be physically connected (to allow flow transfer). To do that, nodes are added at the intersection between the pipes and faces of the 3D elements. These two types of components (1D-3D elements) were joined in accordance with the tectonic history.
Following this, we obtain a 3D model of the URG. The rock matrix of the basement (between faults) is considered impermeable, and the fluid flow takes place within faults. As a first approximation, the hydraulic properties of the fault sets are uniform, with a transmissivity of 10-7 m2/s. The sedimentary cover (3D continuum elements) is considered permeable and hydraulically connected to faults (pipes) in accordance with the tectonic history. The sedimentary cover of the URG is composed of different lithology units. As a first approximation, we considered two main different units where the different units are combined according to their permeability values. The values of the physical parameters used in the model such as the permeability were estimated from geophysical investigations and hydraulic tests made on different geological units [53, 54] and, for instance, have been used in the hydrothermal model of [13]. The deepest one corresponds to Keuper, Muschelkalk, and Buntsandstein units with a permeability of 4.10-14 m2 and a porosity of 0.17 (see the blue layer in Figure 4(e)). The second one (above, see the gray layer in Figure 4(e)) corresponds to Jurassic, Oligocene, Miocene, Pliocene, and Quaternary units, considered less permeable with a permeability of 1.10-18 m2 and a porosity of 0.15 [13]. Note that the selection of 3D elements in order to apply the respective hydraulic properties is done using the DEMs provided by the transnational database [36], such as the position of the base and top of main geological units.
4.3. Initial and Flow Boundary Conditions
A hydraulic no-flow boundary condition is specified for the lateral (north, south, east, and west) and bottom ( km) model boundaries. For the top of the model, constant hydraulic heads are imposed (fixed). The relief of the water table at the regional scale can be either topography-controlled or recharge-controlled [55]. Here, we assume a topography-controlled water table system, which is justified by the fact that the system is in a humid region with a subdued topography and low hydraulic conductivity [56]. The current topographic elevation has been regionalized (averaged by zone), and the corresponding hydraulic head values have been applied as boundary conditions. As an initial condition, the hydraulic heads correspond to the minimum valley centre’s elevation (100 m).
4.4. Numerical Simulations to Assess Favourable Areas
4.4.1. Particle-Tracking Considerations
Once the flow model is built, the steady-state dynamic equilibrium of the groundwater flow system is simulated. Under steady-state conditions, the mathematical formulation of the model is highly simplified. The hydraulic conductivity is the only required parameter, contrary to other parameters such as the storability and specific yield related to transient responses, which are not necessary.
To study flow paths, a particle-tracking methodology that accounts for advective and dispersive transports is introduced (Figures 4(g) and 4(h)). Particles (over a thousand) are positioned close to the border faults of the URG at the top of the model that can be possibly identified as recharge areas [9, 13] (see Figure 4(g)). A series of simulations with different values for the random number generator is carried out to capture the different possible flow pathways. For the flow system modelled, several thousand flow paths are constructed. The simulated time is chosen large enough to allow particles to move out of the system or to be stuck in a zone. At the end of each simulation, an ensemble of flow paths (from recharge to discharge zones) is constructed (Figure 4(h)). For each particle, the following characteristics are extracted at each time step: position (, , and ), absolute time (), properties of the element containing the particle such as its type (1D or 3D for fault or sedimentary unit, respectively), its permeability, the name of the domain (fault sets or lithology units), etc.
4.4.2. Elaboration of Criteria and Indicators
To use the results of the particle-tracking method (flow paths) for geothermal resource characterization, we need to choose relevant hydraulic criteria that constitute a good signature of the geothermal anomalies. Then, based on these main characteristics of geothermal anomalies (hydraulic criteria), indicators are elaborated and extracted from each particle trajectory (Figure 4-Step 3) to use the results of the simulation that enable the identification of the geothermally favourable areas. The following assumptions are used: (1)Geothermal anomalies result from a state of imbalance over a long period where the only efficient mechanism (ubiquitous and active over large spectra of time and space) is the regional groundwater flow. Then, the flow path length (long enough to be regional) and the time required to cover this distance are used as the two main indicators to highlight regional groundwater loops(2)Geothermal anomalies reflect positive thermal anomalies resulting from upward flow in discharge areas due to ascending warm waters. Therefore, the particle must come from a deeper part and has to reach a sufficient depth along its trajectory to carry warm(er) water until the upflow area. The depth reached by the particle along its trajectory, the final position (or a position slightly upstream from it) of the particle along its track, and the flow direction are thus also used as main indicators
As mentioned just above, the indicators that are used are (1) the total distance travelled by the particle from the recharge area (starting point) to the discharge area (ending point) during time , (2) the time spent by the particle in the system to travel between the recharge and discharge areas, and (3) the depth reached by the particle throughout the trajectory. Note that the relevant information associated with each track is extracted to be analysed in a GIS software.
4.4.3. Flow Path Analysis
The main objective of the groundwater flow path analysis is to cross-correlate the properties of each particle trajectory (distance, time, depth, localization, and flow direction) to locate the preferential discharge areas (with specific features of geothermal anomalies) of regional groundwater loops.
After an individual flow path analysis, a spatial and qualitative analysis is realized to localize the preferential discharge areas of regional groundwater loops. A preferential discharge area corresponds to a zone with a high concentration of relevant particles (related to the deep regional flow circulation), which is likely to present thermal anomalies (Figure 4-Step 4). To select relevant particles, a preliminary data analysis is performed. The statistical distribution is observed to define threshold values that are used to assign weighting factors to particles. For example, to select particles that correspond to regional groundwater loops (related to a distance indicator), two threshold values are determined (equal to one-half and three-quarter of the maximal distance). The aim is to highlight the particles that cover greater distances. A weighting factor of 2 is attributed to each particle covering a distance of more than three-quarter of the maximal distance, which represents around 1% of the total particles. A weighting factor equal to 1 is attributed to each particle covering a distance of more than one-half and less than three-quarter of the maximal distance, which corresponds to 14.5% of total particles. By the same reasoning, two threshold values are calculated related to the maximal travel time; then, a weighting factor of 2 is attributed to each particle with a travel time of more than three-quarter of the maximal travel time, and a weighting factor equal to 1 is attributed to each particle having a travel time of more than half and less than three-quarter of the maximal travel time. For the indicator related to the depth reached by the particle along its trajectory, we considered that a particle, which has reached a sufficient depth of more than 2.5 km (related to the average geothermal gradient), is interesting and deserves a weighting factor. Then, a weighting factor of 1 is attributed to each particle that reached a depth of more than 2.5 km along its trajectory. To finish, the sum of weighting factors is calculated for each particle. The idea is that particles with the higher scoring are likely to produce thermal anomalies and correspond to the more relevant particles and inversely, the lower scoring to the less relevant particles. To eventually constitute a preferential target area (thermal anomaly), a final requirement has to be fulfilled: a sufficient density of relevant particles must be localized in a “restricted” spatial area.
4.5. Favourable Area Estimates: The Flow Model
Figure 5(a) presents the results of the flow path analysis at the URG scale where particles are plotted at a 2 km depth (equal to that of the available temperature map, see Section 6.2) as a function of the weighting factors. Particles with a longer distance and time (in accordance with regional groundwater loops) and reaching a depth over 2 km along the trajectory (in accordance with the temperature) are selected. The target areas are then delineated by defining the contours encompassing these relevant particles (Figure 5(a), green lines). Five different preferential areas are highlighted from the results of the 3D flow model, which are distributed as follows: two areas to the north of the URG (surrounding Speyer and between Landau and Mannheim), one area clearly located to the south of Soultz-sous-Forêts, and two areas located near and south of Strasbourg (Figure 5(a)). Figure 5(b) presents the trajectories associated with the most relevant particles providing the preferential areas. The northern areas are mainly related to particles from the eastern shoulder of the URG and along a SE-NW direction. The second area to the south of Soultz-sous-Forêts is mainly related to particles from the western edge of the URG and moving in a NW-SE direction. Note that few particles come from the eastern edge and in a SE-NW direction. The southernmost areas (surrounding Strasbourg) are related to particles from the eastern and western edges.
(a)
(b)
The results of the 3D flow model provide preferential discharge areas related to geothermal resources. Since the ultimate aim is to support expensive exploration campaigns, any information improving the confidence in the area location is valuable. When no in situ data is available, cross-validating the interpretation with results from other models is a way to increase the confidence level. In this study, we propose to build a 3D stress model to extrapolate preferential areas from the mechanical perspective and to finally compare them to the flow model estimates.
5. 3D Stress Model at the URG Scale
We propose to investigate the heterogeneities in the stress distribution that may favour the setup of thermal anomalies. The stress distribution in rock masses is significantly affected by structural and material heterogeneities (including mechanical behaviours/parameters and fault zones) and anthropic activity. Yale [57] gives a review of the in situ evidence of fault impact on the local deviation of stress from the far-field trend. The conclusions are that the main factors locally affecting the stress are the proximity to the fault zones ([29, 58, 59], among numerous others) and the segmentation of the medium (the more fault zones, the more deviation).
To build a 3D stress model of the URG, we use the code 3DEC™ [35].
5.1. Mechanical Properties
A full description of the code and its constitutive equations is provided in Appendix B: 3DEC Numerical Code, and only a short outline is given here. Since we focus on the presence of the fault zones, the matrix behaviour is kept simple and assumed to be elastic (Table 1). The fault zones, however, follow an elastoplastic law (Coulomb slip), with a possible dilation (Table 2). In addition, we assume that any large deformation occurring within the medium is localized into the joints and not into the rock matrix (see Appendix B: 3DEC Numerical Code for more details).
Due to the absence of data at the considered scale, the parameters are extrapolated from estimates obtained by fitting numerical models on in situ experiments [60, 61]. We extrapolated the parameters qualitatively considering fault zones with noncohesive infills (weaker elastic modulus, zero cohesion) and large wavelength profiles (smaller dilation angle, larger critical shear displacement) that are altered by numerous tectonic stages (smaller friction). The parameters are on the orders of magnitude of the parameters from models built at similar scales in the URG [62, 63].
5.2. Boundary Conditions
The displacement boundary conditions 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 presented in Section 3 is embedded in a continuous material having the dimensions of Buchmann and Connolly’s model (Figure 6). Following a procedure used by Peters [64], this embedding enables simulating the far-field nature of the boundary conditions. As a first approximation, the embedding material is assumed to have the same properties as the crystalline basement. In addition to these lateral boundary conditions, the normal displacement is fixed at the model base (roller-type boundary condition), and the top is left free.
5.3. Initial State
In addition to the present-day forces, the stress field is affected by past tectonic events that leave residual stresses through the rock mass. The prestressing method is commonly used to simulate the impact of past events [62, 65–67]. Assuming that the previous tectonic episode occurred a sufficiently long time ago to allow the relaxation of the medium, the initial state can be approximated using a gravitational loading. Gravity is activated within the blocks, and a uniaxial stress condition (equation (1)) is prescribed on the lateral boundaries of the model to account for the influence of the surrounding rock mass. Note that due to the high density of the fault zone network, the initial condition makes the stress distribution very heterogeneous through the rock mass. Once the initial equilibrium is reached, the conditions described in the previous section are applied to the model boundaries, while the gravity is still active within the model: where (Pa), (Pa), and (Pa) are the vertical, maximum horizontal, and minimum horizontal principal stresses (respectively); (m) is the depth; (m·s−2) is the vertical component of the gravity acceleration; (kg·m−3) is the rock mass density, and (−) is Poisson’s ratio.
5.4. Simulations to Assess Favourable Areas
Once the assembly of blocks and joints is balanced (see Appendix B: 3DEC Numerical Code for more details), the model gives an estimation of the mechanical equilibrium through the system. In particular, displacement and stress heterogeneities are obtained and can be inspected to highlight areas where the mechanical state can favour the presence of thermal anomalies. As previously, this requires the definition of an appropriate indicator. The definition of a mechanical indicator is difficult since stresses can act on flows in several ways. However, since stresses can act on the groundwater heads directly (impact on pressure gradients) or on the pore/fracture network (impact on medium effective permeability), most mechanical indicators can be categorized following the volumetric and deviatoric decomposition of the stress tensor: (i)The volumetric part of the stress tensor (e.g., mean stress) is directly linked with the groundwater heads: flow will happen from the more compressed areas to the less compressed zones [68]. “Volumetric” indicators could also relate to modifications of the pore network (e.g., fracturing), but since we consider hard rocks under compressive states, the pore network modification will mostly happen under the deviatoric influence(ii)The deviatoric part of the stress tensor does not affect the pressure heads but can affect the effective permeability of the medium. In the rock matrix, the stress deviator can be responsible for shear-induced fractures that greatly contribute to the permeability increase of the medium. For fault zones, deviatoric stresses result in shear displacements perpendicular to which upflow is favoured [24]. In addition, the shear stress maintains the fault zones in an active state, which is argued to reduce the fault sealing and hence favour the permeability [69]
In the first study, we selected a single scalar indicator to be chosen according to either the volumetric or the deviatoric part of the stress tensor. Denoting (Pa) the stress tensor (see Appendix B: 3DEC Numerical Code for more details), we chose the mean stress as the mechanical indicator: by releasing the efforts acting on the rock mass, the less compressed areas (i.e., low ) facilitate fluid discharges, including vertical upflows of fluids that were possibly heated up to the rock mass below. Low values of the indicators are thus used to identify favourable areas.
The different steps to delineate the favourable areas with the geomechanical model are summarized in Figure 7.
5.5. Favourable Area Estimates: The Stress Model
The equilibrium of the mechanical system is computed through 3DEC™. We obtain a stress map that can be used for delineating the preferential target areas. Prior to postprocessing the stress, the depth at which we search for favourable areas must be set. Rather than a fixed depth, we decided to interpret the results close to the sediment/basement interface. Since the DFN was built according to the fault zone traces observed at the top of the basement, the level of uncertainty in the geometry is more limited close to the interface. The mean stress distribution is presented in Figure 8(a) (the zones located in the horst are ignored).
Five different preferential areas are highlighted from the results of the 3D stress model. They are distributed as follows: two areas in the northwestern part of the URG (to the west of Speyer), one area located around Soultz-sous-Forêts, and two areas located in the southern part (one to the southwest of Strasbourg and one to the east) and to the south of Strasbourg (Figure 8).
6. From Regional-Scale Models towards the Prediction of Favourable Areas for Geothermal Resources
6.1. Cross-Analysis of Favourable Areas
The final step of the work is to estimate areas where the physical conditions 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 the geothermal characteristics and the physics considered. The most favourable areas are assumed to be at the intersection of flow-favourable and stress-favourable areas. They are estimated qualitatively by overlapping the different favourable areas obtained with each model (flow and stress models) (Figure 9(a)).
(a)
(b)
The cross-analysis highlights three main favourable areas: one is in the northern part of the URG, close to Speyer; a second one is to the south of Soultz-sous-Forêts, and a third one is to the south of Strasbourg (composed of three different subareas) (Figure 9(a)). The northern delineated area (west of Speyer) results from the overlapping of the two main favourable areas of the stress (blue circle orientated in the east-west direction) and flow (green circle orientated in the NS direction) models. Close to Soultz-sous-Forêts, the two distinct areas (blue circle orientated in a NE-SW direction and green circle orientated in a NW-SE) result in one main favourable area. Note that these two main favourable areas (west of Speyer and south of Soultz) are located in the western part of the rift (Figure 9(a)). For the four areas (two respective areas for each model corresponding to two green and two blue circles, respectively) located to the south of Strasbourg, the cross-analysis results in three subareas (including two small ones), where the largest one is orientated north-south in the western part of the URG (Figure 9(a)). To summarize, the cross-analysis highlights two main favourable areas that are located in the western part of the URG, to the south of Soultz-sous-Forêts, and to the west of Speyer. A third area can be distinguished to the south of Strasbourg, where three subareas are delineated.
6.2. Delineated Favourable Areas versus the Thermal Map
The main favourable areas from the cross-analysis are compared to the existing thermal map of the URG at a 2 km depth. From north to south, three main local thermal anomalies can be distinguished from the temperature map: one in the north, between Landau and Speyer; one in the centre, around Soultz-sous-Forêts; and one in the south, to the west of Strasbourg (Figure 9(b)). Then, from ~east to west, 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-Forêts in the centre of the URG, where the highest temperatures 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, 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 favourable areas: one in the north (west of Speyer), one in the centre (south of Soultz-sous-Forêts), and one in the south (south of Strasbourg). This distinct distribution (from north 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).
6.3. Discussion of the Cross-Analysis Results
6.3.1. Accuracy of the Areas Delineated through the Cross-Analysis
From a global perspective, the results obtained from the cross-analysis are very encouraging, with the delineation of three favourable areas distributed from north to south along the rift, which is in good agreement with the distribution of the main local thermal anomalies known in the URG. However, the accuracy of the results against the thermal map varies from one delineated area to another. From north to south: (i)the area between Landau and Speyer partially highlights the local thermal anomaly. The cross-analysis results delineate an area to the west of Speyer, approximately 20 kilometres away from the in situ area(ii)the area close to Soultz-sous-Forêts is the best correlated one by the model. Still, the final estimation is to the south of Soultz, while the actual area is most visible to its east(iii)the model results for the area close to Strasbourg are shifted eastwards compared to the thermal map. However, as opposed to the two other areas, the in situ data in this area is much more diffuse. The model results follow the natural tendency of the thermal anomalies, and in this regard, we conclude that this area is not the least precise one
The uneven precision of the model results steams from the cross-analysis of the stress and flow results. While the two hydraulically favourable areas tend to correlate with the actual thermal anomaly located between Landau and Speyer (Figure 9), the mechanical results are less precise and move the final area estimation northwards. In the end, the northern delineated area partially highlights the local thermal anomaly. Conversely, the mechanical results appear to better highlight the thermal anomaly located on the western edge of the rift than the flow model’s ones. According to the 3D flow model results, these hydraulically favourable areas are related to particles originating from the eastern border (see yellow circle, Figure 10(b)) and are characterized by regional flow pathways that mainly occur through NW-SE faults (F4 fault set) in a global direction SE to NW (Figure 5(b)). These few particles go through the centre of the URG (where the sedimentary basin reaches more than 4 km deep) and then migrate through the regional fault network up to the favourable areas. These deep regional flow paths appear consistent with the geochemical evidence [14]. Indeed, the geochemical signatures of geothermal brines collected from the granite basement (Landau, Rittershoffen, Soultz, and Insheim at the graben’s NW borders) suggest that they all reacted with deep sedimentary rock at temperatures close to 225°C at km [14]. This signature involves a migration of geothermal brines from the centre of the URG to the graben’s NW borders through a complex system of deep faults.
For the centre area, the delineated area (close to Soultz) that corresponds to the major thermal anomaly known is reasonably well estimated from a regional perspective. From the 3D flow model, this favourable area results from particles originating mainly from the western edge of the URG (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 (Figure 5(b)). The 3D cross-sectional view (Figure 10(b)) shows that particles move from the west (downward flow) along the URG east-dipping boundary fault to the east (upward flow), mainly along west-dipping faults and associated subvertical faults that are connected to these west-dipping faults (Figure 10(b)). These results suggest deep regional groundwater loops that occur through east-dipping (downward flow) and west-dipping (upward flow) faults that intersect at great depth and along the main NW-SE flow direction and where the main recharge area is located to the northwest of Soultz-sous-Forêts. These preferential flow pathways along the west-dipping fault zones characterize horst structures and are consistent with some previous works [9, 70]. Moreover, according to the results (particle flow paths), we can observe deep flow paths (few particles) from the SE (small blue circle in Figure 10(a)) and in a global SE-NW flow direction emerging close to the Soultz and Rittershoffen areas. These few particles go through the centre of the URG (where the depth of the sedimentary basin reaches more than 4 km) and then migrate through the regional fault network up to the favourable area (which could also be consistent with the results of [14]). The favourable area close to Soultz is the best located one by the cross-analysis because it presents the key characteristics from a hydraulic point of view and also from a mechanical point of view. Indeed, the results of the mechanical model show a less compressed region (Figure 8(b)) where the stress state can favour upflow. Thus, the major thermal anomaly known in the URG is the best located one by the numerical approach due to the combination of numerous key situations: structurally, hydraulically, and mechanically. Structurally, this region corresponds to a horst with a west-dipping fault (enabling the intersection at depth with east-dipping faults). Hydraulically, this area corresponds to a preferential discharge area of deep regional fluid circulation. Geomechanically, the mechanical stress in this region favours fluid circulation and then hot fluid upwelling in the rock matrix and fault.
To finish, for the southern area, the local thermal anomaly is particularly well delineated by the geomechanical model (less compressed area and faults exhibiting opening tendencies, see Figure 8). Second, preferential pathways (from the western border fault) that occur along the west-dipping fault zone also characterize the southern area (to the west of Strasbourg), and we can observe particles from the eastern border fault that emerge in this favourable area (see red circles, Figure 10(a)). All of this information seems consistent with the knowledge provided by the numerous studies realized within the URG (geophysics, geochemistry, etc.). The results demonstrate the potential of this type of numerical approach to provide constraints for the challenge of geothermal exploration and to move towards a better overall understanding of the URG geothermal systems.
6.3.2. Relevancy of the Model Cross-Analysis
(1) Choice of the Physics Described by the Models. The proposed approach is based on the choice of the single physics described by the models. Here, we propose an approach based on these main underlying hypotheses: the regional fault network affects the redistribution of stresses and regional groundwater fluid circulation. Additionally, the stresses significantly affect the flow properties, and there is a strong correlation between the regional groundwater flow and geothermal anomalies. The groundwater flow paths are 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 within the regional fault network, the stresses, and the distribution of the geothermal resources. These previous observations lead us to ask the following question: is it essential to consider the thermal physical processes to highlight the thermal anomalies in the context of geothermal exploration where the ultimate goal is to limit the area to explore and reduce the associated costs? Building a hydrothermal model at the regional scale could facilitate the interpretation of the 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 the definition of criteria and then have an impact on the delineation of the favourable areas. Thermal processes such as the heat transfer between faults and the surrounding matrix and at the intersection of the base of the thermal blanket and the fault (in which the upward transfer of heat occurs) result in redistributions of heat. These redistributions of heat could provide slightly shifted favourable areas in comparison with those currently delineated in this study. Thus, this type of numerical model could be relevant and, more specifically, useful to move towards a better overall understanding of the URG geothermal systems and should provide some constraints on the relative influence of each physics. In any case, building a hydrothermal model at the regional scale constitutes a real challenge and then requires further analysis. The promising results also demonstrate the key role of the regional fault network, which is at the basis of the two numerical physical models. Then, a better delineation of the favourable areas requires an improvement of each numerical model and thereby an improvement of the DFN geometry.
(2) Cross-Analysis versus Coupling. One of the main assumptions of the study presented in this paper is that the estimation of favourable areas can be obtained by cross-analysing results from separate single-physic models. Several authors use either one-way coupled or fully coupled models as an attempt to more accurately address the physical features involved in large-scale phenomena (>kilometric characteristic length) [71–75].
In our case, a solution for achieving one-way coupling would be to use mechanical results for parameterizing the groundwater flow model. Indeed, mechanical results predict huge heterogeneities in stress distribution, including in its deviatoric part (which can be related to second-order fracturation). As a result, the rock mass is expected to exhibit heterogeneous permeabilities (increase of permeability with fracturation). More complex effects could also occur in the fault zones such as permeability decreases in fault zones that are less shear active [69] or increases in permeability in the direction perpendicular to the shear displacement [22, 24].
Full coupling for the considered geometry (large scale and DFN) seems, however, to be unrealistic with the tools at hand today. The fully coupled models found in the literature consider much more reduced scales [74, 75]. Conversely, models exhibiting scales comparable to ours (~100 km characteristic length) incorporate either only a single physics and/or unfaulted media [26, 62, 65, 76]. CPU- or GPU-based parallel simulations show promising results for handling fully coupled models but need further progression to reach the scale and the medium segmentation considered in this study [71–73].
At this stage of development of the study, we anticipate the use of coupled models to introduce more uncertainties than improvements in the overall approach. Given the absence of information on coupled phenomena at large spatial and time scales, the coupling laws extrapolated from smaller scales might be erroneous, and the estimation of the coupling parameters would be very difficult (even qualitatively). In contrast, focusing on the prevailing processes related to the setup of geothermal resources for each model is expected to better constrain the results of the approach (i.e., positions of the favourable areas). In our opinion, the cross-correlation of relatively well-trusted models brings more objectivity in the results than estimation stemming from a coupled but more uncertain model. The case study presented in this paper proves that the cross-analysis tend to estimate reasonably well the thermal map obtained from in situ data and demonstrates how the overlapping results from uncoupled models can be used to refine the final estimates of the favourable areas.
6.4. Perspectives of Improvement for the Proposed Approach
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 with the temperature map at a 2 km depth, which underlines the need for improvements that will constitute the next steps of this work. Based on the results and discussion, we conclude that improvements can be made at three main levels: the geometry description, model precision, and result (cross-)analyses. From this perspective, the following paragraphs detail the opportunities for improvement and their relative importance.
6.4.1. Physical Model Formulations
Both the groundwater flow and mechanical models rely on several simplifying assumptions and could be refined in several ways (equations, parameters, boundary conditions, and initial states). However, building almost exact models is an unachievable task given the large spatial and time scales involved in light of the absence of parameter measurements and evolution of tectonic states and related mechanical boundary conditions [8]. Targeting a limited number of improvements could be carried out by hierarchizing them based on their impact on the results. Most notably, the boundary conditions and the use of nonheterogeneous parameters for the groundwater flow and mechanical models, respectively, are expected to be mostly responsible for the uneven accuracy of the areas they delineate (see Section 6.3.1.).
Nevertheless, the reasonable agreement of the model results with the thermal observations leads us to believe that the physical model accuracy, both in terms of state equations and parameters, is a secondary aspect in the approach. Rather, the geometry description seems to be of more importance, as is usually the case with models based on a dense DFN.
6.4.2. Geometry Description
The DFN geometry (regional fault network) requires a more detailed description. The geodynamic evolution of the rift system and existing conceptual models of the rifting zones are important sources of information to improve the DFN construction [77–79]. The interpretation of the two deep seismic lines illustrates the structural asymmetry of the Rhine Graben, which is characterized by thickness variations of the tertiary deposits and the presence of boundary faults in the north along the eastern flank and in the south along the western flank. The asymmetric evolution of the rift from north to south that affects the deep structural features has to be considered. Similarly, horizontal bedding planes and a simple conjugate fault pattern characterize the Cenozoic fill in the north, while in the south, a more complex structural pattern exists, with the presence of tilted fault blocks from west to east, conjugate normal faults, and a smooth anticlinal structure. Then, the steeply dipping normal fault that controls the western margin of the Cenozoic graben and merges at the midcrustal levels provides constraints on the fault system and its evolution with the depth. Building the geometry by integrating this type of constraint could directly affect the connectivity of the medium as well as the stress redistribution due to the model segmentation and consequently the groundwater flow paths and the geographical distribution of the discharge areas.
6.4.3. Result Analyses
The final estimations of the proposed approach depend highly on the way in which the model results are interpreted. First, considering the two models separately, the relation between the model results and the characteristics of the geothermal resource rely on the definition of the indicators for each model. Recalling that a particle-tracking method is used for the groundwater flow, the retained indicators are the residence time, the travelled distance, the depth reached by the particles, the flow direction, and the concentration of relevant particles. The first three indicators are linked with the temperature (fluid heating up the rock mass), the fourth one indicates resource upwelling, and the final one relates to the discharge areas. Because of the rather straightforward effect of flows on the thermal anomalies, we are rather confident in the chosen indicator. However, the relative interpretation of the particles could be improved through a rigorous statistical treatment. In this paper, and as a first approximation, we chose to set the weighting factors qualitatively to focus on the overall approach. Regarding the mechanical model, the impact on the thermal anomalies is only secondhand, through affecting the flow and/or thermal processes. In this paper, we retained the mean stress , arguing that less compressed areas favour fluid migration. The mean stress seemed to be a reliable indicator for identifying the true thermal anomalies. We remind here that, in accordance with the whole approach, the mechanical criterion can only be interpreted qualitatively, i.e., so as to highlight less compressed areas close to more compressed ones. In such cases, fluid migration is assumed to move towards less compressed areas, then move upwards due to stress reduction when getting closer to the Earth’s surface. In our opinion, a quantitative analysis is not realistic at this stage given the lack of information on the rock mass properties at the considered scale (law of behaviour and associated parameters, behaviour heterogeneities, second-order fracturation…). cannot be guaranteed to be the best one. Other indicators were tested such as the preferentially opened fault zones. Taking advantage of the discrete nature of the model, the opening tendency can be estimated through the fault zone normal displacement (m). More specifically, and to focus on the impact of the present-day stress state, the displacement increment can be plotted: , where (m) is the normal displacement due to the initial equilibrium. Figure 11 presents the contours of the fault zone averages ( (m2) being the area of the fault zone) in order to interpret the tendencies at the fault zone scale. NW-SE faults exhibit opening tendencies in contrast to EW and NE-SW faults that exhibit a tendency to close. Estimated favourable areas tend to be more widespread (light blue circles in Figure 11).
Interestingly, the results in the surroundings of Soultz-sous-Forêts correlate well with the interpretation in this region, tending to reinforce the idea that this region combines numerous key situations. Comparison of different mechanical results highlights the difficulty in choosing a mechanical criterion, which in the end might be a cross-correlation between rock matrix and fault zone criteria. On a broader perspective, these results question how the mechanical results should be interpreted: in a straightforward manner (as is done in this paper) or to parameterize the flow model? Thorough investigations must be carried out in the future to shed light on these interrogations.
Second, the cross-analysis of the model results itself could be improved. Here, as a first approximation, we proposed to set the favourable areas by choosing the overlapped regions that were qualitatively delineated. In the studies found in the literature (see, e.g., [80, 81]), the authors propose to quantitatively establish favourability maps by interpolating in situ data and using weighting factors according to the measurement confidence. In our case, because the physical models are built under several assumptions, it seems slightly too early to propose quantitative favourability maps. The more precise data at our disposal are the structural data. A first step could be to weight the results according to how close (and thus less extrapolated) they are to our confident structures.
Finally, depending on the ratio between the obtained and desired scales for the preferential target areas, additional systems might have to be set and run to progressively refine the estimation precision. If the proposed approach is used as a support for an exploration campaign whose aim is to site exploration wells, the results of a first region-scaled model might be too diffuse or uncertain by, e.g., providing several distinct areas in the same region. If more precise results are expected (e.g., if the range of the final estimates is too large or if there are concerns about the locations of the favourable areas), steps one to four could be repeated on progressively smaller-scaled geometries until the desired size is achieved or more confidence is obtained for the final estimates of the favourable areas. In this case, the results from larger-scaled models can provide boundary conditions for the smaller-scaled ones (see, e.g., [76]).
7. Conclusion
This paper is a current overview of our team’s efforts 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 geological system scale. The study relies on the building of the first 3D numerical flow (using the discrete-continuum method) and mechanical uncoupled models (using the distinct element method) at the URG scale (faulted geological setting). Studying the interactions of these main processes while coupling them with the effects of the regional structural elements (fault zones) helps us to understand their relative impacts on the distribution of geothermal resources. Here, we applied this approach in the URG context (Figure 12), where available data (such as seismic lines) allow the building of these types of models and a direct assessment of the results (using, for instance, the temperature map). Based on the fact that there exist strong links between the stress, fluid flow, and geothermal resources, we followed these main steps.
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 of the URG. This major step constitutes a challenge alone and marks the beginning of new opportunities to understand the key effect of the regional fault network on the geothermal resources. Indeed, this question remains an open question and still constitutes a challenge for geothermal engineers, who attempt to explain hydrothermal processes in this faulted geological context.
Then, we decline the DFN model of the regional fault network to build the first 3D flow and stress models of the URG using adapted numerical tools that explicitly consider the key role of the faults. The respective results issued from each numerical model are indirectly interpreted through the definition of criteria (using the link between the physics considered and the setting up of geothermal resources) that enable the elaboration of indicators to identify geothermally favourable areas related to the considered physics. For the flow model, the hydraulic criteria are related to the preferential discharge areas of the deep regional fluid circulation, and the indicators used are related to the properties of the particle trajectories (distance, depth, etc.). For the stress model, the stress criteria are related to the mechanical state that can favour the presence of thermal anomalies, and then, the related indicator chosen is the stress tensor (mean stress). As a final step, we proposed to compare the hydraulically and geomechanically favourable areas to locate thermal anomalies. A cross-analysis of the results of the 3D flow and stress models within the URG provides three preferential target areas for geothermal projects. These preferential areas are distributed along the rift as follows from north to south: to the west of Speyer, to the south of Soultz-sous-Forêts, and to the south of Strasbourg. The flow and stress models result on a common, rather condensed, overlap for two of the three areas (Speyer and Soultz-sous-Forêts) which are related to major local thermal anomalies currently identified in the western part of the rift. The overlapping is more diffuse for the third area where the in situ data is more scattered as well. Given that no parameter fitting was done, the delineated favourable areas and local known thermal anomalies tend to be in reasonable agreement, which underlines the potential of such a cross-numerical approach. These promising results suggest that the numerical approach (discrete approach to explicitly consider the effect of the fault network), the physics considered (stress and fluid flow) in the geothermal resource setup, and the method to interpret the results (criteria and associated indicators) are relevant to capturing the thermal processes and then can help to locate thermal anomalies. In other words, the geographical distribution of the geothermal resources is strongly correlated with the preferential discharge areas of the deep regional fluid circulation occurring within the regional fault network and related to the less compressed areas that favour these upwellings of deep hot fluids.
Moreover, a detailed analysis of each favourable area demonstrates the ability of this type of approach to move towards a better overall understanding of the URG geothermal systems. For example, the best favourable area delineated by the proposed approach corresponds to the major thermal anomaly known within the URG. A more detailed investigation of the delineated favourable area close to Soultz-sous-Forêts enables the highlighting of the key situations (structurally, geomechanically, and hydraulically) that make it a special location for geothermal resources. Structurally, the west-dipping fault resulting in the intersection at depth with east-dipping border faults forms a key structural feature. This V-shaped structure favours loops of deep fluid circulation. Hydraulically, this area corresponds to a preferential discharge area of deep regional fluid circulation. In addition, the mechanical stress in this region favours fluid circulation and then hot fluid upwelling in the rock matrix. In conclusion, the Soultz-sous-Forêts area is probably the major thermal anomaly because this area merges structurally, geomechanically, and hydraulically favourable conditions. This evidence is consistent with previous studies and demonstrates the primary role of the regional fault network. Then, the constraints provided by the respective models demonstrate what we can learn from these first models 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 very encouraging results underline the potential of such a cross-numerical approach to help locate geothermal resources. This type of approach could ensure the localization of four-fifths of all thermal anomalies. Indeed, the hydrothermal convection along faults (activated due to the tectonic context) may explain 75-85% of all temperature anomalies [21], as demonstrated in the URG context. This method proposes practical help in locating and estimating the extension of hidden thermal anomalies, with the prospects of reducing (1) the required number of exploration wells through the optimization of the exploration area and (2) the cost of each exploration well by minimizing the drilling depth needed to reach the target temperatures. Moreover, the opportunity for a better match between resources and needs offered by the 3D regional models will facilitate the deployment of deep geothermal energy in the territory and its effective integration in local and regional distribution networks while minimizing the environmental impacts. To conclude, the first 3D numerical flow (using the discrete-continuum method) and mechanical models (using the distinct element method) at the URG scale offer new research opportunities. The ultimate goal of our team’s work is to provide innovative exploration technologies/methodologies to limit the area to explore and reduce the risks associated with geothermal exploration.
Appendix
A. 3FLO Numerical Code
The 3FLO code is based on the finite element method. 3FLO can solve different types of problems, such as [82] (i)flow in fracture networks, represented by a 3D network of pipes. A network of channels (called pipes) is generated on each fracture. The connection of the channels from one plane to another through the fracture plane intersections (called tubes) results in a 3D network of 1D elements(ii)flow (saturated or unsaturated) in porous media by using mixed-hybrid 3D finite elements. These elements are more precise than classical (Galerkin) 3D elements and allow a proper computation of face fluxes. This helps minimize solute transport computation errors, as will be detailed later on. Solute transport equations can be used to assess the favourable geothermal areas(iii)flow interactions between fractures and porous media by allowing the superimposition of 1D elements generated on disk-shaped fracture planes (to simulate fractures) on 3D elements (to discretize the rock matrix)(iv)transport simulated by the particle-tracking method
The flow equation solved in the 1D classical elements is where (m3s−1) is the pipe conductivity; (m) the hydraulic head; (m) the abscissa along an element; (m2) the section of the element ; (m−1) the specific storage coefficient of the element ; and (s) the time. A pipe is formed by two infinite parallel ribbons, open enough so that a flow regime like the flow between two parallel plates can take place. To satisfy the laminar flow regime which dominates in rock masses, the pipe conductivity is correlated to a cubic law. The pipe section is then proportional to the cubic root of its conductivity.
For a flow plane with a regularly spaced grid of channels, the transmissivity of the fracture (m2/s) is related to the pipe conductivity and (m), the spaced grid size:
The diffusivity equation is solved in 3D elements: where (m·s−1) is the permeability of the porous media; (m) the hydraulic head; and a source term (s−1).
Solute transport is solved by using the random walk method. The flow is one-dimensional everywhere in pipes and tubes, except at intersections. The dispersion is therefore only longitudinal and is “completed” by the full mixing occurring at intersections. A pipe porosity (-) can be specified, as well as a pipe dispersivity (m). In a pipe, the particle displacement is rectilinear and uniform: where (m), (m·s−1), (s), and are the particle position along the pipe, the fluid speed, the time step, and the pipe porosity, respectively. The 3FLO code enables the simulation of the particle displacement in channels and in 3D elements. Indeed, in 3D elements, once the velocity is known, the movement of a particle can be separated into a longitudinal displacement along a flow line (corresponding to advection and longitudinal dispersion) and two orthogonal transverse displacements (simulating transverse dispersion). Moreover, direct exchanges are possible between 1D and 3D elements; particles can originate from a pipe and enter a 3D element and vice versa.
B. 3DEC Numerical Code
B.1. Overview of 3DEC Software
The mechanical software (3DEC™ for the three-dimensional Distinct Element Code, see Itasca 2016) relies on the distinct element method, where the contact forces between solids are explicitly incorporated. Initially proposed to study slopes in jointed rock masses [83], the DEM has been developed ever since, especially in recent years. Today, its scope of application has extended to rock mass scales [84–86] down to microscopic scales [87–89]. In addition to supporting discontinuous displacement and stress fields, the particularity of the DEM is to consider the full mechanical behaviours for the contacts. When the contacts depict fault zones, the law of behaviours embraces the behaviour of the infills and of the hanging and foot walls (which belong to the rock mass), and the law of behaviours can be chosen accordingly.
The DEM thus offers the possibility to (i)describe discontinuous displacement and stress fields(ii)explicitly account for mechanically active discontinuities and their impact on blocks(iii)use complex mechanical laws for the joints
The solution phase of the DEM is similar to that of any continuous numerical method: at each time step, the solver runs until mechanical equilibrium is achieved within the system. For the deformable blocks, the problem unknowns are the displacements (m) at the mesh grid points. is computed by solving the balance in where (m3) and (kg) are the volume and mass for the considered block, (Pa) the stress tensor, (N) the sum of external forces other than gravity forces, (m·s−2) the gravity vector, and (m) the displacement, and a top dot denotes a time derivative.
With the DEM, the external forces account for the interactions with the contiguous blocks. These interaction forces are obtained through the joint constitutive equations which, given a prescribed displacement, return the resulting forces. That is, in addition to the constitutive equations that must be given for the continuous methods, the DEM also requires joint constitutive equations to solve equation B.1. The complete solution scheme at each time step is then achieved by balancing the displacements resulting from the block deformations with the forces resulting from the block interactions (Figure 13).
At the whole-system scale (here, the URG model), the state equation solved in 3DEC™ is the mechanical balance of the assembly of blocks and joints. To be solved, the state equation must be completed with the constitutive equations, the boundary conditions, and the initial state.
B.2. Constitutive Equations
The mechanical constitutive equations, or laws of behaviour, give the relation between displacements and stresses. For the DEM, the constitutive equations must be provided for both the blocks and the joints. In this study, because we work in a highly segmented medium where the fault network is expected to have the greatest impact on the stress redistribution, the emphasis is set on the fault zone characterization.
B.2.1. Rock Matrix (Blocks)
A simple behaviour is considered for the rock matrix, and blocks are assumed to be linearly elastic, homogeneous, and isotropic materials: where (Pa) is the elastic stiffness tensor, (Pa) and (Pa) the current and initial stress tensors, respectively, and (−) the strain tensor. A convenient way of formulating the equation splits the strain tensor into volumetric (−) and deviatoric (−) parts: where (Pa) and (Pa) are the bulk and shear moduli, respectively. The volumetric strain is the trace of the strain tensor, and the deviatoric part is given by . The primary unknown in the mechanics is the displacement (m), and its relation with the strain is yielded by the compatibility condition:
B.2.2. Fault Zones (Joints)
A more sophisticated behaviour is considered for the fault zones. Fault zones are mechanical weaknesses in the rock mass and accommodate more deformation than the surrounding matrix. Because of the larger deformations involved, irreversible processes can more easily take place within the fault zones and should be considered in the model. Recalling that joints are flat surfaces in 3DEC™, the joint behaviour can be described along the normal and parallel directions to address opening/closing and shearing phenomena, respectively. With the Coulomb slip law, the shear behaviour is characterized by a bilinear relation: an elastic ramp and a plateau as soon as the onset of plasticity is reached: where (Pa·m−1) is the joint shear stiffness, (Pa) the shear stress, (Pa) the plastic yield stress, (m) the shear displacement, and (m) the plastic yield displacement. The plastic yield stress depends on the fault zone internal parameters and on the normal stress (Pa) acting on it: where (Pa) is the cohesion and (°) the angle of internal friction. When the onset of plasticity is reached , the exceeding displacement is permanent, hence the “plastic” denomination. When the joint lies in the plastic phase, another mechanical feature is triggered: dilation. Dilation is the opening of the fault zones under shearing and is related to the fault zones actually being irregular surfaces rather than flat planes. The normal behaviour of the joint is thus elastoplastic as well and reads where (Pa·m−1) is the joint normal stiffness, (m) the elastic part of the normal displacement, and (m) the irreversible part of the normal displacement due to dilation. The dilation relation introduces two additional material parameters, the dilation angle (°) and the critical shear displacement (m) above which dilation vanishes:
Both dilation parameters are conditioned by the geometry of the fault zone irregularities: relates to the shape, and depends on the maximum height of the irregularities.
The joint formulation remains valid for large displacements. The block equations however are restricted to infinitesimal strains (see equation (B.4)), implying that any large strain at the model scale would be localized into the joints.
B.2.3. Block and Joint Interactions
For each contact joint, the normal and shear stresses add up to the joint contact force (N), which depicts the interaction force between the two blocks sharing the joint: where (m2) is the contact area and (−) and (−) the unit normal vectors in the normal and shear directions, respectively, along the joint. For each block, the external force vector is expressed as the sum of the contact forces due to the joints surrounding the block: where (N) is the contact force of the -th joint surrounding the block.
Since the mechanical balance of the blocks explicitly depends on (see equation (B.1)), the impact of the joints on the blocks is directly expressed through this vector. Conversely, the impact of the blocks on the joints is obtained once the block mechanical balance is achieved: the balance modifies the displacement field through the model, which will affect the forces acting on the joints through equations (B.5) and (B.7).
Data Availability
The data used to build the models and support the findings of this study are included within the article or can be found within references cited (articles, database issued from the GeORG project results: http://www.geopotenziale.org/home?lang=2).
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was conducted in the framework of the IMAGE project (Integrated Methods for Advanced Geothermal Exploration, grant agreement No. 608553). This work was completed under the REFLET and TEMPERER projects (French National Research Agency grant agreement Nos. ANR_10-IEED-0802-02 and ANR-10-IEED-0801-02, respectively) which are supported by the French Government through the Investments for the Future programmes. This work also received internal funding from the French Geological Survey.