Influence of the Main Border Faults on the 3 D Hydraulic Field of the Central Upper Rhine Graben

The Upper Rhine Graben (URG) is an active rift with a high geothermal potential. Despite being a well-studied area, the three-dimensional interaction of the main controlling factors of the thermal and hydraulic regime is still not fully understood. Therefore, we have used a data-based 3D structural model of the lithological configuration of the central URG for some conceptual numerical experiments of 3D coupled simulations of fluid and heat transport. To assess the influence of the main faults bordering the graben on the hydraulic and the deep thermal field, we carried out a sensitivity analysis on fault width and permeability. Depending on the assigned width and permeability of the main border faults, fluid velocity and temperatures are affected only in the direct proximity of the respective border faults. Hence, the hydraulic characteristics of these major faults do not significantly influence the graben-wide groundwater flow patterns. Instead, the different scenarios tested provide a consistent image of the main characteristics of fluid and heat transport as they have in common: (1) a topography-driven basin-wide fluid flow perpendicular to the rift axis from the graben shoulders to the rift center, (2) a N/NE-directed flow parallel to the rift axis in the center of the rift and, (3) a pronounced upflow of hot fluids along the rift central axis, where the streams from both sides of the rift merge. This upflow axis is predicted to occur predominantly in the center of the URG (northern and southern model area) and shifted towards the eastern boundary fault (central model area).

In general, 2D numerical models for the central URG (e.g., [13,14,[16][17][18]) predicted a topography-driven deep fluid flow characterized by a downflow at the borders of the URG and an upflow in the western part of the graben center.These simulations included one or more faults in the western part of the URG (mainly near to Soultz-sous-Forêts), which resulted in a preferred upflow of fluids along these faults.According to previous fluid chemistry characterization studies, it is likely that mixing of meteoric water and deep saline reservoir fluids occur at either both boundary faults (e.g., [25][26][27][28][29]) or at the W border fault [12,18].Fluid circulation is indicated by temperature profiles, for example, at Soultz-sous-Forêts and Rittershoffen (e.g., [28,30,31]).In summary, most studies proposed that deep convective circulation of hot fluids plays a major role in the hydraulic regime of the URG [8,12,14,18,28].Whether the deep circulation is occurring only in the fault zones or if cross-formation flow might exert additional influence on the graben hydraulics is to the author's knowledge still to be quantified (e.g., [29]).
Although the URG is a well-studied and utilized area, an integrated three-dimensional understanding of the main controlling factors for the coupled fluid and heat transport is still lacking.Freymark et al. [24] have assessed the lithosphere-scale conductive thermal field by means of 3D simulations.Despite that their conductive simulation reproduced the general trend of observed temperatures reasonably well, local deviations between observed and predicted temperatures occurred mainly at shallower depth levels (upper 2-3 km) of the URG.This latter aspect indicates that an additional component related to heat transported by moving fluids should be considered as well.
Following on these aspects, we present a detailed numerical investigation aiming at quantifying the role of fluid circulation on the three-dimensional thermal configuration inside the main sedimentary units of the URG.In general, deep subsurface fluid flow is controlled by (i) pressure gradients 2 Geofluids imposed on the system by variable recharge through the shallow subsurface, (ii) the hydraulic properties of the subsurface geological units as differentiated by the existence of lithological variations and structural discontinuities (faults and fractures), and (iii) variations in fluid density.To implement the hydrogeological configuration of the central URG into the 3D numerical simulations of coupled heat and fluid transport, we make use of the upper part of the data-based 3D structural model of the URG [24] complemented by the well-known geometry of the main border faults of the graben.By implementing general trends in surface pressure and temperature conditions affecting the underground, these numerical experiments provide new 3-dimensional insights into the regional patterns of groundwater flow and related thermal anomalies.Of special interest in the current study is (1) to quantify, by means of a sensitivity analysis, the influence of the main border faults on the basin-wide fluid flow in the central URG and (2) to evaluate the implications for the deep temperature distribution.

Method
2.1.Starting Model and Model Area.The 3D structural model of Freymark et al. [24] provides the basis for the geological configuration and thus for the lithology-dependent distribution of thermal and hydraulic properties in our numerical model.The structural model covers the URG, the western Molasse Basin, and the Hessian Depression and consists of 14 lithostratigraphic units, 7 upper crustal domains, a lower crystalline crust, and 2 domains of the lithospheric mantle.
The geometries of these model units were derived based on different observables (well data, reflection and refraction seismic data, geological maps, and existing 3D structural models; see references in Freymark et al. [24] for more information on the data base).In addition, 3D gravity modelling was performed to assess the internal structure of the crystalline crust while integrating all available seismic lines.Within this study, we focus on the subsurface above a depth of 8 km of the central URG, which hosts several sites currently used for geothermal energy production like Bruchsal, Landau, Insheim, Soultz-sous-Forêts, and Rittershoffen (Figure 1).In addition, the study area comprises the area of Strasbourg, which is currently under geothermal exploration.To define the model domain for the coupled thermo-hydraulic simulations, we have combined several lithostratigraphic units with comparable properties according to Freymark et al. [24].This results in 3 main sedimentary units within the URG that have distinct hydraulic and thermal rock properties (Table 1, Figure 2): (1) the permeable Cenozoic sands and marls; (2) the less permeable Keuper, Lias, and Dogger sediments; and (3) the permeable Muschelkalk, Buntsandstein, and Permo-Carboniferous carbonates and sandstones.Below the sedimentary units, an upper crystalline crustal layer of variable thickness has been included (at least 2 km thick).The laterally juxtaposed crystalline crustal units correspond to the upper crustal units of Freymark et al. ([24], Figure 2) and are parameterized accordingly (Table 1).In particular, the radiogenic heat production differs for the three crystalline crustal units.The base of the model is a flat surface at a depth of 8 km below sea level.
In addition to the 3D volumes representing the geological units, we have implemented the two main border faults of the URG as discrete model features, corresponding to those implemented in the model of GeORG [37].These faults are represented as slightly undulating surfaces dipping steeply towards the graben center (Figure 2).They extend from the topography to a model depth of 6 km b.s.l. as proposed by Baillieux et al. [30].
To adequately represent these structural attributes in the numerical simulation, a fully unstructured mesh is required.However, meshing complex geological structures are still a technical challenge [52].Each geological unit is defined as a closed 3D volume that shares the nodes with the neighboring geological units or faults.Faults are discretized as 2D elements embedded in the 3D finite element mesh such that each node of the fault is shared with the adjacent geological units.As layers and faults are defined by scattered points before meshing, their intersections have to be defined.A specific challenge is that the 3D volumes enclosed by the finite elements should not be too thin and that sharp angles should be avoided.Problems occur especially at thin, out-pinching layers.To have a good representation of the actual geometries, a large amount of elements are needed to mesh such layers, which increases the computational time for the simulations.Therefore, a balance has to be found between an adequate representation of the geology and an acceptable number of elements and thus computational time.
For our simulations, we have built a fully unstructured mesh consisting of about 1 million tetrahedral elements (Figure 2) using the software MeshIt [53].The mesh resolves the major structural features of the graben configuration in terms of geological units of varying permeability as well as thermal properties and in terms of the 3D configuration of the major border faults of the graben.
In addition, two main structural approximations have been implemented: (1) the thin out-pinching geometry of the Keuper/Lias/Dogger unit is approximated by extending this geological domain further to the North with a finite thickness of 50 m.Thus, the unit comprising Keuper, Lias, and Dogger is 50 m thicker than in reality in the northernmost model area (Figure 3(b)).( 2) The thin outcropping sediments on the graben shoulders were neglected.As a consequence, the graben shoulders are considered as upper crystalline crust without sedimentary cover.
For the numerical simulations, hydraulic and thermal properties are assigned to each model unit derived either from lab measurements or from former modelling studies (Table 1).Where lab measurements are available for subunits (e.g., Muschelkalk, Buntsandstein, and Permo-Carboniferous), we use the weighted mean considering the average thickness of each geological sublayer.Furthermore, bulk densities are calculated considering matrix densities derived by 3D gravity modelling [24] and measured porosities (Table 1).
The upper crystalline crust is differentiated into the main Variscan domains (Mid-German Crystalline High, Saxothuringian, and Moldanubian; Figure 2, Table 1).In the East of the central model area, the top of the crystalline basement is deepest reflecting the asymmetric rift geometry of the URG (Figure 4).However, in the northern model portion, this asymmetry is less pronounced.In contrast, in the South the basement is deeper at the western border of the URG (Figure 4).Since previous studies have shown that the locally highly fractured basement can contain deep saline groundwater (e.g., [29,55,56]), the basement is characterized as not completely impermeable in this study (3E − 18 m 2 [37]; Table 1).

Parameterization of the Boundary Faults.
To test the influence of the main border faults on the hydraulic system, we run a suite of simulations each different in terms of the hydraulic configuration of the faults.An overview of all the model runs is given in Table 2.We have tested configurations without faults (model A), with faults of different width (models B and C) and different permeability (models B/C, D, and E).We have refrained from testing scenarios with faults acting as a barrier to fluid flow since the border faults of the central URG are generally described as hydraulically active conduits (e.g., [13]).The permeability value of 5E − 14 m 2 was chosen as a "realistic" value, as it is based on an interpretation of pumping tests on the GRT-2 well in Rittershoffen [57].The other value of 1E − 12 m 2 was chosen to test the response of the system to a much higher permeability (model D).Hence, model A and model D represent end-member scenarios with the least and largest effects on fluid flow expected, respectively.With a much higher permeability at the eastern border fault (compared to the western fault; model E), we have tested if a one-sited increase in fluid velocities can cause a shift in fluid flow directions (and lead to an upflow of fluids in the western parts of the graben as proposed by various authors, e.g., [12,13,18,29]).Since we do not have graben-wide information on the internal structural variability of the faults (e.g., on their differentiation into fault cores and damage zones with finite widths), they are modelled with a horizontally and vertically homogeneous permeability value acting across a homogeneously wide domain.To further amplify the difference in the tested scenarios, the two permeability values have been coupled with two different values for fault width (1 and 10 m, respectively).With the different scenarios tested (Table 2), we attempt to systematically quantify the sensitivity of the hydraulic and thermal regime across the entire study area to variations in the hydraulic behavior of the faults.

Boundary Conditions.
The base and the lateral model boundaries are defined as no-flow boundaries in terms of fluid flow.Atmospheric pressure (0.1 MPa) is assigned to the top of the model as upper hydraulic boundary condition.Thus, the hydraulic head is fixed and it is equal to the topography (Figure 1; etopo1, [23]).The system is modelled as being fully water-saturated.Accordingly, largest pressure heads are prescribed in areas of highest elevations corresponding to up to 1070 m in the Black Forest and the Vosges Mountains in the SW and SE model area (Figure 1).Lowest pressure heads correspond to elevations of 100-200 m along the Rhine river bed (Figure 1).Despite being a simplification of the real surface hydraulic system, our choice of hydraulic boundary conditions provide a first-order approximation of the surface regional flow as believed to occur in the graben (e.g., [13]).
At the surface, the annual average surface temperature is assigned as upper thermal boundary condition (Figure 5(a)).High-resolution measurements of the annual average surface temperature in Germany [58] were complemented by a global data set for the French part of the model area [59].Accordingly, highest average surface temperatures of 11 °C are found in the URG, while coldest average surface temperatures of 5 °C characterize the topographic highs in the Vosges Mountains and the Black Forest.

Geofluids
The temperature distribution of the lithosphere-scale conductive thermal model of Freymark et al. [24] at 8 km depth below sea level was extracted and assigned as lower thermal boundary condition at the base of the coupled model (Figure 5(b)).This assures that the contribution to the global heat budget from the deep crustal domain is considered.Highest basal temperatures (up to 330 °C) occur below the URG in the northern central part of the model while coldest temperatures (up to 250 °C) are displayed in the north-western model area.The resulting pattern of basal temperature matches the lateral distribution of radiogenic heat production of the different Variscan crustal domains, superimposed by thermal blanketing from the Cenozoic sedimentary rocks [24].

Numerical Simulation.
To simulate the coupled transport of heat and fluid, we use the open-source software GOLEM [60].Golem is a flexible, parallel scalable finite elementbased simulator for thermo-hydro-mechanical process modelling in fractured porous rocks as based on the MOOSE framework [61].In a first step, the steady-state conductive thermal field and the steady-state pore pressure distribution are calculated.The results are used as initial conditions for the thermo-hydraulically coupled simulations.As proposed by Kaiser et al. [62], we assume that forced convection is largely suppressing free convection in this setting of high gradients in hydraulic head (here equal to topography).Thus, only pressure-driven advective heat transport is considered while fluid density and viscosity are assumed constant in our simulations.The coupled fluid and heat transport is calculated by solving the equations based on conservation of (1) fluid mass, (2) fluid momentum (here implemented as a linearized Darcy's law), and (3) internal energy as  7 Geofluids with temperature T, time t, porosity φ, fluid modulus k f , pore fluid pressure p f , Darcy velocity q D , permeability k, fluid viscosity μ f , fluid density ρ f , gravity acceleration g, bulk specific heat ρc b , bulk thermal conductivity λ b , specific heat capacity of the fluid c f , and radiogenic heat production H. Details of the numerical implementation are given by Cacace and Jacquey [60].
A pseudo-steady state was achieved by allowing the simulation to equilibrate for up to 2 Ma.In all runs presented in the study, we observed the reaching of pseudo-steady conditions after approximately 700 ka (Figure 6).These 700 ka represent a numerical simulation time, namely, the time after which quasi-steady-state conditions were achieved (in terms of modelled temperature).At this stage, it can be assumed that the interaction of all hydraulic and thermal controlling factors (implemented boundary conditions, physical properties, etc.) was effective throughout the complete model domain.Hence, allowing for steady-state conditions to appear is a prerequisite for fully assessing and understanding these interactions (in terms of fluid flow fields and related temperature variations).

Results
3.1.Hydraulic Field of the URG.The fluid flow field calculated based on the distribution of hydraulic properties and the hydraulic boundary conditions is exemplarily described for model scenario B (Table 2; Figures 7 and 8).
In general, all model scenarios show a regional fluid flow from the surrounding basement into the sedimentary infill of the URG (Figures 7, 8, and 9(b)).Along the main border faults, fluids flow downward as, for example, west of Soultz-sous-Forêts (Figures 7(b) and 9(b)), but locally also upward (Figure 7(b)).Inside the URG, most of the basin-wide fluid flow is perpendicular to the rift axis (Figure 7(a)).In the rift center, streams from NW and SE merge, further flow towards N/NE (Figure 7(a)), and thereby form a pronounced upflow axis (Figures 7(b), 8, and 9(b)).In the central model area (compare Figure 1), this upflow axis is located near to the E border fault (Figures 7(a), 8, and 9(b)), while in the northern model area it is located close to the center of the URG (Figures 7(a) and 8).
For the entire suite of model scenarios (Table 2), predicted fluid flow velocities range from a minimum of 2 7E − 8 m/year within the basement to a maximum of 5.2 m/year in the sedimentary infill of the URG (Table 3).Even higher fluid flow velocities are predicted for the fault planes (max.54 m/year; Table 3).
Compared to the fluid flow field of scenario B, the simulation without discrete faults (model scenario A; Table 2) shows no differences in the described basin-wide fluid flow trends, except for the fact that no fluid flow is predicted along the location of the main border faults (Figures 9 and 10(a)).
For model scenarios C, D, and E, differences in fluid flow velocity and directions are predicted locally along the main border faults (Tables 2 and 3; Figures 9 and 10).As the faults have a higher permeability than most of the surrounding geological units (models B-E), they provide preferential pathways for fluid flow.The higher the fault permeability, the higher the computed flow velocities inside the fault (Table 3).At the same time, the existence of such pathways generally enforces the downward directed flow, i.e., the infiltration rates in their direct vicinity (compare Figures 10(a) and 10(d)).This deep infiltration effect reaches even into the basement west of the fault for the scenario with an extremely conductive western fault (end-member model D).This model also shows that strong downward flow related to the fault hardly affects the basement east of the fault while it reduces the flow across the fault into the sediment fill of the graben.In general, fault-induced changes in the overall, graben-wide flow dynamics are hardly visible even when considering a very high permeability (1E − 12 m 2 ) at one (Figure 9(e)) or both border faults (Figure 9(d)).
In contrast to the sensitivity of the flow fields with regard to fault permeability, variations of the fault width as chosen in the range of 1 to 10 m have no significant effect on the hydraulic field (Figures 9(b), 9(c), 10(b), and 10(c)).
The comparison of model scenarios differing in the parametrization of the main border faults allows us also to uncover characteristics of the regional fluid flow field that are common to all models and thus independent of the fault behavior.Given the chosen setup of the models    2).On top, the topography and thus the assigned hydraulic head are shown.Rectangle in (a) shows the area zoomed into in Figure 10.
10 Geofluids   2, color code for geological units as in Figure 9).In contrast, the temperature distributions predicted by the coupled simulations of heat and fluid transport (Figures 11(d)-11(f)) reflect the fluid flow directions.Colder temperatures are predicted along the borders of the URG, where cold water is infiltrating the system (Figure 12).Higher temperatures are predicted around the center of the rift, where flows from the E and W borders merge and circulate upward due to forced convection.The most significant positive thermal anomaly is predicted at approximately the center of the URG caused by a pronounced upflow parallel to the graben axis (Figure 12).In addition, smaller anomalies are predicted that spatially correlate with the geothermal exploitation areas of Soultz-sous-Forêts (Figure 12(b)) and Bruchsal (Supplementary Material 1).Those thermal anomalies are caused by local forced convective upflow inside the Cenozoic sediments in addition to slightly upward-flowing fluids inside the basement (Figure 9).Such an upflow is also predicted for Rittershoffen (Figure 12(b)) and Landau (Supplementary Material 1).However, in these areas the upflow in the Cenozoic sediments is less pronounced than in the areas of Soultz-sous-Forêts and Bruchsal, causing thermal anomalies of smaller magnitudes.
Depending on the different fault configurations tested (Table 2), further differences emerge also for the temperature distribution in the respective models of coupled heat and fluid transport.We have analyzed predicted temperatures at a constant depth of 3 km below sea level, but at different distances with respect to the western border faults to compare the different models.As expected, the largest differences in modelled temperatures are observed between the end-member scenarios (models A and D) at smallest distances from the fault center (10 m): these scenarios differ by <9 °C in the absolute temperatures predicted.At a distance of 5 km from the border fault, however, the end-member scenarios differ only by 1-2 °C.All other scenarios tested differ by <2 °C in absolute temperatures, regardless of the distance to the fault center.We take these modelling results as an indication that the graben-bordering faults can be responsible for observed thermal anomalies in their direct proximity (at distances of <5 km).Thermal anomalies occurring in the central parts of the graben (at >5 km distance), however, require additional heterogeneities in hydraulic and thermal properties (as induced by lithological variabilities or additional faults and fracture zones not yet implemented in the model).While such border-fault-related variations in hydraulic properties influence temperatures only locally, it has to be emphasized that the coupled simulations in general reveal that the component of fluid flow results in first-order variations in the resulting 3D regional thermal configuration with respect to purely conductive heat transport (Figure 12).

Discussion
4.1.Observational Evidence.Because of its scale and due to a scarcity of relevant information, it is difficult at this stage to validate the presented regional models with observations.The modelled fluid dynamics, however, manifest themselves in local upflow, respectively, downflow, areas which can be compared to the locations of mineral and thermal springs (Figure 13; [69][70][71]) as well as observed artesian conditions (data only available east of the river Rhine; Figure 13; [72]).
Comparing those observations with the modelled Z -component of fluid flow velocity, a spatial correlation is visible between modelled upflow areas and observed springs (e.g., model scenario B in Figures 7(b) and 13).Most of the observed springs spatially coincide with predicted upflow areas, whereby the correspondence is more evident in the northern parts of the model area than in the south.
The model, however, does not predict upflow at all locations of mineral and thermal water springs or observed artesian conditions.Most deviations of the model from these observations are located on the graben shoulders, where modelled fluid flow velocities are significantly lower (down to 5E − 8 m/year; Figure 7(a)) and boundary effects near the closed model boundaries are likely to occur.Such discrepancies between the regional model and local observations can be related to the limited structural resolution of the model that, given the large lateral extent, does not integrate details of the local lithological differentiation and fracture and fault networks.To fully explain and reproduce these local observations, further studies resolving these local conditions will be required.
Another general feature of our findings is that the modelled upflow axis spatially coincides with the location of the river Rhine (Figure 13).This shows that not only the real river path follows the gradient of lowest topography, but that this area of maximum discharge is also the domain where uprise of warm fluids is strongest in response to the 3D variation of hydraulic pressure in the subsurface.
One of our major interests initiating this study was in the assessment of the influence of groundwater flow on the thermal field.We have already shown that the difference in temperatures predicted by the coupled thermo-hydraulic models with respect to a purely conductive thermal model is significant.For a comparison with observed thermal anomalies, there is a large number of temperature measurements available at boreholes spread over much of the central URG (314 temperature measurements from 75 wells [33,[63][64][65][66][67][68]; 13 Geofluids Figure 14).Only temperature data of best and good quality according to the respective author [33,[63][64][65][66][67][68] were used, while hydraulically disturbed temperature logs were not considered.Those temperature measurements are available with different spatial coverage for different wells and regionally show a large spatial variability in absolute values.For example, south of Strasbourg measured temperatures in 2 wells, which are located 1 km apart from each other, show a difference of 32 °C at the same depth level.Beside these horizontal variations, measured temperatures show large variations with depth.For example, in one of the wells south of Strasbourg, two measurements at a vertical distance of 23 m show a temperature difference of 15 °C.
One trend that can be derived from the calculated temperature differences (simulated minus observed temperature; Figure 14) is that close to the western border of the graben, modelled temperatures tend to be too low.However, also too high temperatures are predicted distributed over the whole model area.Strongest misfits (positive and negative) are located at a depth range between 700 and 2000 m below sea level.Comparing different well logs, several trends are evident.For some wells in the South, calculated temperatures underestimate observed values in the upper part, while the fit to observed temperatures in the deeper parts is better.In contrast, for some wells, as for example along the eastern border fault, predicted and observed temperatures fit very well in the upper parts, but the model is too cold in the deeper parts.
As, in general, the modelled temperatures are too low at the borders of the URG, the effect of the hydraulic upper boundary condition must be discussed.By assigning the hydraulic head to the topography, large hydraulic gradients are prescribed at the borders of the graben.These high gradients likely lead to an overestimated cooling effect.However, measured hydraulic head data show the same regional trends as the topography, which means that the modelled trends in imposed fluid direction can be regarded as meaningful.It is rather the absolute temperature values that are likely too low due to overestimated fluid velocities.
Figure 15 shows exemplarily the comparison of temperatures measured in two different wells with temperatures that were predicted for the same location by the conductive and the TH coupled model scenario B. It is evident that the temperature profile of Eschau 06 (Figure 15(a)) can be reproduced only in the shallow and in the very deep parts, while the small-scale variations between −400 and −1400 m depth cannot be resolved with the regional models.In the shallow parts of this well, the prediction of the TH coupled model is even better than the purely conductive model.In contrast, the temperature log of Bruchsal 2 in the northern model area (Figure 15(b)) is well reproduced by both models (coupled and purely conductive) in the upper 1300 m, but a discrepancy between measurement and simulations is obvious in the deeper parts.Whereas the TH coupled model is colder than the measured temperatures, the conductive model is too warm in the upper Keuper/Lias/Dogger unit and too cold below −1750 m depth.Both regional models cannot fully resolve the small-scale variability in measured temperature.
Thus, our comparison of modelled and measured temperatures reveals discrepancies that cannot be correlated with certain predefined components of the model.The spatial variability of temperature misfits is much smaller in scale than the spatial extent of model units and major trends in boundary conditions.Hence, changing the parametrization of the model in its current structural form and resolution would not improve the overall fit.Neither would a more detailed analysis of misfits lead to additional conclusions with respect to regional fluid flow trends (as controlled by the border faults in particular).Improving the fit of the model would require implementing local heterogeneities in hydraulic and thermal properties as could be related to lithological variabilities and the existence of fractures and faults.

Modelling Approach and Related
Insights into the Regional Fluid Flow Field.The models presented in this paper are the first 3-dimensional representations of deep regional fluid flow and related thermal anomalies for the area of the central URG.The 3D numerical models have been set up in a way to implement (i) major hydraulic and thermal rock property variations as can be derived from a data-driven geological model [24] and (ii) general trends in surface  temperature and pressure conditions assumed to control the subsurface thermo-hydraulic system.
The different model scenarios presented in this study have shown that major graben-bounding faults influence the fluid flow locally depending on their permeability, but do not significantly change the basin-wide hydraulic field.
We are able to identify aspects that are common to all model scenarios and thus define typical characteristics of the regional fluid flow.
In agreement with former studies (e.g., [13,18,29,73]), our simulations predict recharge dominating in the high-topography areas flanking the URG and discharge   prevailing in the center of the graben (Figures 7-9).Due to the high hydraulic gradient and the discontinuity between high-permeable sediments and the low-permeable basement, fluid flows downward along the borders of the rift into the sedimentary infill.In addition, fluids migrate slowly through the deep basement from the areas outside the URG thereby entering the sedimentary infill in areas characterized by low hydraulic potential (Figure 9).The most prominent difference of our 3D coupled simulations to former 2D studies is the location of the dominant upflow axis, which in our simulations locates in the (eastern) center of the URG.This result is in contrast to conclusions derived from previous studies predicting preferential upflow across the western center of the URG (e.g., [12,13,18,29]).The upflow axis in the presented study results from pressure-driven forced convective processes occurring throughout the basement and the overlying porous sediments due to the lowest hydraulic potential being located in the center of the URG (compare Figure 1).Accordingly, where the lowest hydraulic potential is imposed near the eastern border fault (as in the central model domain; Figure 1), the upflow axis is shifted towards the eastern center of the URG (Figures 1 and 7).
This difference to previous 2D numerical simulations (e.g., [12][13][14][16][17][18]) results from the consideration of the three-dimensional pressure variations in our numerical model, whereas former 2D models assumed cylindrical symmetry, and thus constant pressure conditions, perpendicular to the plane of the modelled section.With a limited 2D perspective on the topographical differences between the eastern (higher) and western (lower) graben shoulders (Figure 9), one might expect a prominent east-to-west-directed flow in the graben and upflow toward the west (cf.[13]).In contrast, due to its 3D nature, our model also captures the pressure-driven north/northeastward-directed component of the flow (Figure 7) caused by the overall south-north topographic gradient (i.e., the prevailing trend from the Alps to the Rhenish Massif).The resulting south-north potential is most significant inside the graben, where it also interacts with thickness maxima of the more permeable sedimentary units that tend to be located east of graben axis (Figure 3).Thus, we propose this interaction (of north-directed gradients and accumulation of sediments) to suppress E-W-directed flow due to the differences in hydraulic potentials as induced by the differently high graben shoulders.
Given that the setup of the presented regional models implies a number of simplifications, it is worth discussing how the main findings of this study concerning (i) the impact of border faults and (2) the main characteristics of regional fluid flow would be altered by changing the model input parameters.

Vertical Resolution and Parameterization.
To test the influence of a vertically larger number of hydraulic units on the hydraulic field, we have performed the simulation presented in Supplementary Material 2. By vertically differentiating three additional sedimentary layers inside the central URG, however, we find the same main characteristics of the hydraulic field as in the series of the simplified models presented above.As mentioned in Section 4.1, the differentiation of hydraulic and thermal properties on the basis of regionally traceable geological units is not suitable to correctly reproduce local observations.The heat production of the Saxothuringian upper crust, for instance, is modelled as an assumed average value on a regional scale of 2.5 μW/m 3 [47], but locally increases up to 7 μW/m 3 [74,75].An important example for large regional variability in hydraulic parameters would be the facies variations inside the sedimentary units not accounted for in this study.Hence, addressing questions concerning the fluid dynamics of a specific site in the graben system requires setting up local models describing property variability with a critical precision.Such regional models presented here would, however, be suitable to providing the required boundary conditions for local models.4.2.2.Fault Parameterization.Another simplification of the presented models concerns the structural and parametric representation of the main border faults.Cacace et al. [3], for example, studied the influence of the inner structure of a fault zone on the hydraulic field and demonstrated that a tight fault core inside a permeable damage zone can cause local differences compared to simulations implementing a homogenous fault.Moreover, the depths of the border faults in our sensitivity analysis represent a compromise between the proposed depths of former studies.Although some models of the area of Soultz-sous-Forêts (e.g., [30,73]) propose a depth of −6 km for the main border faults as well, most modelling studies integrate the border faults to an average depth of around −3 km (e.g., [13,16]).In contrast, interpretations of deep seismic lines propose a depth of −15 km and even deeper for the main border faults (e.g., [76][77][78]).Our choice to limit the faults to −6 km depth is motivated by the assumption that permeability generally decreases with depth so that the fault loses its role as a fluid pathway at larger depths.In general, the way we have modelled the main border faults characterizes the presented series of models as a conceptual approach towards a better understanding of the overall impact of fault properties on the graben-wide fluid flow.The models thus are not suitable for correctly reproducing local observations and predicting processes at specific locations near the fault (which would require additional local information).4.2.3.Advective vs. Convective Heat Transport.One more simplification of the presented TH coupled simulations is that only pressure-driven advective heat transport is considered while fluid density and viscosity are assumed constant.Previous studies proposed free convection (1) in the crystalline crust below the URG (e.g., [14]), (2) at the basement/sediment-boundary (e.g., [12,18]), and/or (3) in the fault zones inside the URG sediments (e.g., [8,30]).First of all, even without allowing for density-driven flow to occur in our models, there is a large number of smaller-scale convection systems predicted for the domains of the sediments and the basement as well as across the basement/sediment-boundary (Figure 9).Hence, no thermal and density-driven flow is required to explain such observations described by previous studies.In addition, it is worth mentioning that the conductive thermal field used as initial model predicts the highest deep temperatures in the eastern part of the URG (Figures 5(b) and 11(c)).In consequence, buoyancy-driven upflow, if existing, should occur in the eastern center of the URG and downflow at the borders and shoulders of the rift.Thus, the main characteristics of the flow regime would not be changed by considering temperaturedependent fluid density and viscosity.
Beside these arguments inherent in the modelled situation, the formation of density-driven convection cells in general requires very specific conditions with respect to the structure of an aquifer and the interaction with hydraulic-head imposed pressure gradients so that previous studies concluded that free convection is less likely to occur on the regional scale in nature than forced convective flow [62].A differentiation of the URG sediments into a larger number of geological units as would be advisable according to local observations (e.g., [32]) would result in thinner aquifers divided by additional aquitards.As a result of this increased vertical heterogeneity, free convection inside the sedimentary infill of the URG would be even less likely to occur.Bjørlykke et al. [79], for instance, have shown that impermeable layers of only <1 m thickness already can split or even inhibit the formation of convection cells.
Beside temperature dependency, also variations in fluid salinity can change the density of the fluid.However, it is shown that salinity is quite heterogeneously distributed across the URG (e.g., 4-20 g/l at the western border reaches up to 120-200 g/l in Bruchsal and Bühl [12]).Thus, we assumed a constant fluid density on a pure water approach in our study, as we did not have an appropriate data base for implementing such complexity into the model.4.2.4.Inner-Rift Faults.Whatever the physical process controlling fluid flow (free or forced convection) in specific parts of the graben system, it is the distributed occurrence of smaller-scale faults and fractures that seems to be decisive for providing the hydraulic pathways (e.g., [8,18,30]).For example, in Soultz-sous-Forêts most fluid flow is observed in the vicinity of deep fracture zones [80].Accordingly, a 1-2 km thick hydrothermal alteration zone in the uppermost parts of the highly fractured basement is described [81,82].However, this thick hydrothermal alteration zone seems to be a local phenomenon [81].In Rittershoffen, the hydrothermally altered granite has a thickness of approximately 200 m only, while no such zone is found in Landau and Insheim [81].Given these large differences in the vertical extent of a strongly fractured basement domain and the lack of a regional coverage of such information, we had to refrain from implementing a corresponding additional unit into the model.Instead, we have chosen a constant permeability value for the crystalline crust down to a depth of −8 km that considers the existence of secondary permeability; i.e., representing a value still allows for fluid flow in the basement instead of regarding the crust completely impermeable.
Such inner-inner faults are observed across wide parts of the central URG, in particular within Triassic and deeper geological domains where convection is described to occur inside fracture systems connecting the basement with the main reservoirs, the latter being bounded by the impervious Muschelkalk (e.g., [81]).Based on our sensitivity analysis revealing that even the largest faults at the graben borders do affect the flow fields only in their immediate vicinity (as described in Section 3), we conclude that these inner-rift faults would also impose only localized modifications of the regional fluid flow (and heat transport).On the other hand, if pervading the geological units with a high spatial frequency, they might impose to the units a spatially continuous (and probably anisotropic) increase in hydraulic conductivities (secondary permeability).As mentioned above, in general more observations (such as from hydraulic leak tests) would be required to validate the hydraulic parametrization of the different geological unitswith respect to both the consideration of structurally controlled pathways (faults and fractures) as well as the upscaling of average matrix permeabilities as derived from laboratory measurements on rock samples.
Whether fluid flow is occurring parallel to the URG inside the inner-rift faults (e.g., [8]) or if the faults are connected by convection perpendicular to the URG (e.g., [18]) is still debated.A larger impact on the basin-wide fluid field would anyway only be expected if there are large impermeable fault zones that act as hydraulic barriers.Peters [83] has performed a slip tendency analysis for a large number of faults with known geometries in the URG.This analysis has shown that there are only a few faults that are almost perpendicular to the maximum horizontal stress and thus have a low slip tendency and most likely a very low permeability.These faults strike predominantly in the ESE-WNW direction, i.e., almost parallel to fluid flow direction proposed by our simulations, and thus they would not significantly change the graben-wide fluid flow.
In general, previous studies (e.g., [13,18,29,30,73]) explain major thermal anomalies, such as the one at Soultz-sous-Forêts, as the result of a high basal heat flow and upward convection of deep and hot fluids within inner-rift faults.Our model demonstrates that the phenomenon of small-scale thermal anomalies may occur even without the existence of inner-rift faults.Thereby, our models do allow for upflow from the basement up into the Cenozoic (Figures 9  and 10), the latter being modelled with the highest permeability (according to laboratory derived measurements; Table 1) although locally being known as lower permeable and not influenced by fluid flow (e.g., as in Soultz-sous-Forêts [11]).It remains to be investigated why despite the overall high permeabilities supporting the upflow of warm fluids, the modelled temperatures along the western parts of the graben still tend to be too low (Figure 14).4.2.5.Boundary Conditions.Finally, we have taken strong assumptions concerning the boundary conditions applied to the thermo-hydraulic simulations.Beside the hydraulic contrasts between the basement and the sedimentary graben fill, it is the hydraulic upper boundary condition that has a major effect on the calculated fluid flow fields (e.g., Figure 9).By assigning the hydraulic head to the topography, we assume the underground to be completely filled with groundwater and we likely overestimate the hydraulic gradients and thus fluid flow velocities and their effect on the modelled temperatures.This may also affect the time to reach steady-state conditions in the simulations.Extensive infiltration of cold water forced in areas of high hydraulic gradients could lead to an overestimation of deep cooling and a prolonged time of equilibration.
An alternative upper boundary condition would be interpolated hydraulic head data.Those data exist, however, only at certain measuring stations (e.g., data from the federal state office of the environment Baden-Württemberg), while using the completely known topographical variability to mimic hydraulic head gradients is possible for the entire model domain.Setting the hydraulic head equal to topography does not mean that we deny the existence of an unsaturated zone in the study area.Moreover, we are aware that the modelled fluid dynamics thus may locally not reproduce reality in terms of absolute valueswhich anyway is difficult to test given the scarcity of observations in terms of absolute fluid velocity and direction in the deep subsurface.With our modelling approach, we assume that the continuously available topographical trends mimic the shallow pressure conditions in terms of the locations of highs and lows.Compared to hydraulic head data that typically represent a subdued replica of topography, we present models with more extreme hydraulic gradients, in particular between the graben shoulders and the graben center.This imposes a maximum hydraulic effect to the main border faults that are located where the gradients tend to be largest.Hence, our conclusion that their effects on the regional fluid flow are minimal would be confirmed by applying hydraulic head data as upper boundary condition.Similarly, our main findings about regional flow directions seem to be confirmed when speculating about more realistic upper boundary conditions.Since the graben center hosts a major river system, the differences between topography and hydraulic heads there should be minimal while maximizing towards the graben shoulders.The related reduction of recharge rates from the graben shoulders towards the graben center would then be associated with a relative increase in the role of graben-parallel flow and low potentials related to larger thickness of more permeable sediments (Figure 3).Hence, setting up our models with observed hydraulic heads most probably will also lead to a pronounced north-and northeast-ward flow and the formation of an upflow axis in the graben center.
The northward flow predicted by our models casts some doubts also on the southern and northern lateral boundary conditions that are modelled as being closed to fluid flow.At the southern boundary, one would expect a stronger inflow into the modelled domain and thus even a reinforcement of the northward flow.It is more difficult, however, to speculate about the impacts of opening the northern boundary to fluid flow.The currently closed boundary does not seem to be associated with a significant upflow as could be expected at least (Figure 7(b)).A solution for the problem of appropriately setting lateral boundary conditions (given the lack of corresponding fluid pressure data) would be to set up hydraulic models that are even larger than the investigated central URG and derive locally predicted pressure conditions from these models.

Conclusion
In this conceptual study, we performed first 3-dimensional numerical simulations of coupled fluid and heat transport in the URG.With our focus on the influence of the main border faults on the 3D hydraulic field of the URG, we gained valuable new insights into the hydraulic system of the URG: (1) A general northward flow is predicted for the subsurface of the URG by our model, which indicates the importance of 3D effects (2) The main border faults are less important than the hydraulic head and the permeability contrast between sediments and basement (3) Different permeabilities and widths of the main border faults in a geologically reasonable range have no significant effect on the graben-wide hydraulic and thermal field (4) Upflow (forced convective) is predicted even without considering density-driven flow (5) The regional deep fluid flow is directed from the margins (recharge at topographic highs) towards the center of the URG (discharge at topographic low) (6) Previously proposed upflow in the western parts of the graben is difficult to explain by this model, but could be related to the difference between 2D and 3D (7) Inner faults might play an important rolealso for the thermal anomalies observed (8) More realistic hydraulic boundary conditions would probably lead to more realistic resultshowever, the previously mentioned main conclusions stay valid Finally, we strongly encourage further studies to test the influence of, e.g., the parameterization of the geological layers, density-driven flow, and more realistic boundary conditions on the hydraulic and thermal field of the URG.

Figure 1 :
Figure 1: Topography (etopo1 [23]) of the model area in the central Upper Rhine Graben (URG) with selected areas of geothermal exploration (S: Strasbourg) and exploitation (L: Landau Pfalz; I: Insheim; B: Bruchsal; SsF: Soultz-sous-Forêts; R: Rittershoffen).The rectangular model area is elongated parallel to the main rift axis (NNE-SSW) and has a horizontal extent of 87 × 153 km.The red rectangle in the map on the top right indicates the area of the larger Rhine Graben model of Freymark et al. [24] based on which the geological configuration of the current model (blue rectangle) is taken.The map on the left is shown in UTM32N and rotated counter-clockwise by approximately 20 °.

Figure 4 :
Figure 4: Depth of top basement.The map is shown in UTM32N and rotated counter-clockwise by approximately 20 °.

Figure 5 :
Figure 5: (a) Annual average surface temperature [58, 59] as upper thermal boundary condition and (b) temperature distribution at 8 km depth below sea level from the conductive thermal model of Freymark et al. [24] as lower thermal boundary condition.Maps are shown in UTM32N and rotated counter-clockwise by approximately 20 °.

Figure 6 :
Figure 6: Temperature evolution with ongoing simulation time for model scenario B. Each line represents a point inside the model at which a temperature measurement is available.In total, 314 temperature measurements were available for this study[33, 63- 68].Exemplarily, two different trends are marked in blue and red.For more details, refer to the main text.

Figure 7 :
Figure 7: Hydraulic field (stream lines) of model scenario B, taken as exemplary for all scenarios where border faults have been implemented.(a) Fluid velocities and (b) z-component of the fluid velocity that differentiates the flow into upward (red colors) and downward (blue colors) flow.For both (a) and (b), the colored flow lines are representative for the entire depth range of the model.Maps are shown in UTM32N and rotated counter-clockwise by approximately 20 °.

Figure 8 :Figure 9 :
Figure 8: 3D fluid flow illustrated exemplarily for model scenario B. Locations of the cross sections were chosen to run through active geothermal power plants and geothermal exploration areas (L: Landau; B: Bruchsal; SsF: Soultz-sous-Forêts; R: Rittershoffen; S: Strasbourg; vertical exaggeration ×2).A zoom into the central part of the middle cross section (between vertical dashed lines) is presented in Figures 9 and 10.Details of the northern and southern cross sections are presented in the Supplementary Material 1.

4. 8 Figure 10 :
Figure 10: Fluid flow at the W border fault in the central model area for different model scenarios (for location refer to Figure 9; for model explanation refer to Table2, color code for geological units as in Figure9).
contrast, in the central model area the upflow is predicted near the eastern border fault (4) a fluid flow which is partitioned between the upper Cenozoic aquifer and the deeper Permo-Mesozoic aquifer by the Mesozoic aquitard with highest flow velocities in the upper aquifer.Fluid flow is slowest in the Mesozoic aquitard and the crystalline basement (Table3)3.2.Thermal Field of the URG.Considering the additional heat transport by fluid flow results in clear differences in predicted temperature distribution when compared to a purely conductive model ([24], initial thermal conditions in this study).Figure11illustrates these differences exemplarily by comparing the temperature distributions at 1, 2, and 3 km depth below sea level predicted by the conductive simulation and by the simulation of coupled heat and fluid transport (exemplary model scenario B).It is evident that the temperature range at the same depth level is different for the two types of models and that the wavelength of temperature variations is considerably smaller in the model of coupled heat and fluid transport.The initial thermal field as derived from a purely conductive thermal model (Figures11(a)-11(c)) is characterized by a long-wavelength thermal anomaly in the URG, caused by (1) the higher basal heat input in response to the high radiogenic heat production of the Saxothuringian upper crust and (2) the thick thermally low conductive sediments[24].Temperature maxima shift from the southern model area at shallow depth (Figure11(a)) to the northern model area at greater depth (Figures 11(b) and 11(c)), which is related to the thickness distribution of the insulating Cenozoic sediments.At a depth of −3 km (Figure 11(c)), the thermal blanketing effect of the thick sediments in the North has a much stronger effect than at −1 km.

Figure 11 :
Figure 11: Temperature distribution at 1, 2, and 3 km b.s.l.depth below the URG predicted by the simulation of purely conductive heat transport (a, b, c) and by the simulation of coupled heat and fluid transport (exemplary model scenario B; d, e, f).The black line indicates location of the cross section shown in Figure 12.Maps are shown in UTM32N and rotated counter-clockwise by approximately 20 °.

Figure 12 :
Figure 12: Cross sections through the (a) initial conductive thermal model and (b) coupled thermo-hydraulic model scenario B (location of cross section as for Figure 9; SsF: Soultz-sous-Forêts; R: Rittershoffen).

Figure 13 :
Figure 13: Observed mineral and thermal water springs (green; [69-71]), observed artesian conditions (red; only available for East of the Rhine; [72]) and the river Rhine (light blue line) on top of the Z-component of the fluid flow velocity predicting upflow areas in red/orange colors (model scenario B, compare Figure 7(b)).The map is shown in UTM32N and rotated counter-clockwise by approximately 20 °.
data points Simulated temperature (conductive model) Simulated temperature (TH coupled model B)

Figure 14 :
Figure 14: 3D perspective view on the border faults (grey shaded areas) and temperature measurement points [33, 63-68].Points are color-coded according to the temperature differences (with model scenario B simulated minus observed temperature) (vertical exaggeration ×3; maximum depth of the faults 6 km b.s.l.).

Table 2 :
Fault parameterization in the different model scenarios.

Table 3 :
Simulated fluid velocity ranges for each model scenario and geological unit in m/year.