Influences of hydraulic boundary conditions on the deep fluid flow in a 3D regional model (central Upper Rhine Graben)

Knowledge of groundwater flow is of high relevance for groundwater management or the planning of different subsurface utilizations such as deep geothermal facilities. While numerical models can help to understand the hydrodynamics of the targeted reservoir, their predictive capabilities are limited by the assumptions made in their setup. Among others, the choice of appropriate hydraulic boundary conditions, adopted to represent the regional to local flow dynamics in the simulation run, is of crucial importance for the final modelling result. In this work, we systematically address this problematic in the area of the central part of the Upper Rhine Graben. We quantify how and to which degree different upper boundary conditions and vertical cross-boundary fluid movement influence the calculated deep fluid flow conditions in the area under study. Robust results, which are insensitive to the choice of boundary condition, are: (i) a regional groundwater flow component descending from the graben shoulders to rise at its centre and (ii) the presence of heterogeneous hydraulic potentials at the rift shoulders. Contrarily, results affected by the chosen boundary conditions are: (i) calculated flow velocities, (ii) the absolute position of the upflow axis, and (iii) the evolving local flow dynamics. If, in general, the investigated area is part of a supra-regional flow system—like the central Upper Rhine Graben is part of the entire Upper Rhine Graben—the inflow and outflow across vertical model boundaries need to be considered.


Introduction
The Upper Rhine Graben (URG, Fig. 1) is a 300 km long continental rift system, flanked by strongly uplifted rift shoulders and extending along parts of the central European Alpine foreland (Lopes Cardozo and Behrmann 2006). It is part of the European Cenozoic Rift System (Lopes Cardozo and Behrmann 2006) extending from the Mediterranean Sea to the North Sea (Dèzes et al. 2004, Fig. 2a). The URG stretches in SSW -to NNE direction from the Jura Mountains in the south to the Rhenish Massif in the north, and its width varies between 30 and 40 km (Fig. 1). It is flanked by outcropping crystalline rocks of the Black Forest and the Vosges Mountains, and its main depocentre is filled with permeable sediments.
The geometry and morphology of the URG, in combination with the permeability of its sedimentary infill, has led to the establishment of a large-scale fluid circulation system (Lampe and Person 2002). Prominent subsurface thermal anomalies are present in the URG, which are currently exploited for geothermal energy production (Agemar et al. 2013).
Previous studies agreed on the presence of a regional groundwater flow from the Alps to the Rhenish Massif mainly driven by differences in the water table topology and locally enhanced by local recharge systems at the graben flanks. Understanding the hydrodynamics within the URG has been the subject of numerous investigations (Aquilina et al. 2000;Clauser and Villinger 1990;Freymark et al. 2019Freymark et al. , 2017Koltzer et al. 2019;Lampe and Person 2000; Länderübergreifende Organisation für Grundwasserschutz am Rhein 2013; Person and Garven 1992;Stober and Bucher 2015;Thierion et al. 2012). It is worth mentioning that the previous models rarely cover the whole extent of the URG (Fig. 1), and are rather limited both on the horizontal and vertical scales, and that the different studies have lead to partly contradicting conceptual models.
As an example, Aquilina et al. (2000) and Stober and Bucher (2015) investigated the fluid flow in the URG along a 2D profile extending from the Vosges Mountains (west of URG) through the URG to the Black Forest (east of URG). The flow regime was inferred based on chemical analysis of fluid samples. Both models predicted the main recharge area to be located at the graben shoulders (Black Forest and Vosges Mountains). From there, meteoric water was thought to flow into the graben sediments following the main hydraulic gradients. Aquilina et al. (2000) inferred ascending of these waters across the graben faults, preferentially in the form of convective cells. Stober and Bucher (2015) predicted hot deep fluids to ascend rather vertically through the Buntsandstein and Muschelkalk, and, though to a lesser extent, along the major fault systems. Both models, possibly because of their 2D character, proposed only graben-perpendicular flow. This last conclusion is in contrast to the main results from the study by Bächler et al. (2003) that indicated with 3D simulations a dominant graben-parallel flow direction, while neglecting any graben-perpendicular flow. The main limitations from this first family of models was to be limited to a two-dimensional representation of the regional hydrodynamics. More recent numerical models that considered the full three-dimensional characteristics of the natural groundwater flow (e.g., Freymark et al. 2019;Koltzer et al. 2019) have illustrated the variability of deep groundwater flow with respect to its causative forces (advective versus free convective fluid flow). Based on 3D-coupled thermal and hydraulic numerical simulations of the URG, Freymark et al. (2019) addressed the influence of the border faults' width and permeability on the hydraulic and thermal field. To parameterize the surface hydraulic boundary forcing, they relied on the commonly used assumption of representing the groundwater surface potential by the model topography. As main outcome, these authors suggest a basinwide graben-perpendicular flow from the graben shoulders to the central URG, with an upflow axis approximately below the river Rhine and close to the eastern border fault in the central part of their model area. This is partly in contrast to the two previously described models by Aquilina et al. (2000) and Stober and Bucher (2015).
Despite all past and recent efforts, the causative processes responsible for the flow dynamics in the URG are yet to be fully understood. In this study, we test the reliability of the model of Freymark et al. (2019) and analyze its sensitivity with respect to different boundary forcings. In doing so, we rely on numerical modelling techniques, more precisely, we use the software FEFLOW © (Diersch 2014) to test different boundary condition scenarios. It is of interest to quantify the influence of hydraulic boundary conditions on the resulting modelling outcomes as the latter are critical in controlling the characteristics of a groundwater system (De Filippis et al. 2017;Franke et al. 1987) and can influence the calculated hydrothermal field (Frick et al. 2015. When cross validating their modelling results against available borehole temperature measurements, Freymark et al. (2019) noted a systematic misfit in their model, which they related to their simplistic hydraulic boundary settings, leading to an overestimated hydraulic potential at the rift shoulders.
Inspired by this latter conclusion, we identified the following two open questions for the present study: (a) Does a water table as a subdued replica of the topography used as an upper hydraulic boundary condition represent a reliable alternative in the URG? (b) What role do open/closed vertical/lateral boundaries play for hydraulic models of the URG?
To answer these questions, we systematically investigate the influence of different approximations of the hydraulic boundary settings on the regional groundwater dynamics in the URG. We make use of the 3D structural model of Freymark et al. (2019) as it represents the most up-to-date integrative model of the subsurface. In a first step, we characterize the first order features and related driving forces of the hydraulic flow processes under a fixed boundary condition setup. This model realization will serve as a reference comparative model. In a second step, the influence of different hydraulic boundary conditions is tested: a potentiometric head surface based on available measurements of the  groundwater table and open vertical boundaries along the  model borders, where inflow and outflow are considered. We build on previous studies to properly characterize the structural configuration of the study area (GeORG-Projektteam 2013a, b, c, d), and the resulting hydraulic regime, the reliability of the latter being constrained by available hydraulic measurement data provided by geological state agencies and societies in France and Germany, respectively.

Hydrogeological setting and structural model
Due to the availability of numerous geological and geophysical data, the structure and evolution of the URG is well known. Figure 2b sketches the URG with the graben's main border faults, the course of the river Rhine and different geological units. The rocks encountered in the URG are the products of the geological evolution from Paleozoic to recent times. For our simulations, we make use of the structural model of Freymark et al. (2019). Relying on this study, we shortly introduce the different geological units in chronological order together with their hydrogeological characteristics and refer to Freymark et al. (2019) for more information on the structural model.
The pre-Permian basement as lowermost model unit consists of the crystalline crust consolidated during the Variscan Orogeny (Kroner et al. 2011). From north to south, five main segments of the Variscides below the URG are named by Kossmat (1927) and Scholtz (1930) as Rhenohercynian (RH), Northern Phyllite Zone (NPZ), Mid-German Crystalline High (MGCH), Saxothuringian (ST) and Moldanubian (MD, Fig. 2b). The crystalline crust is in the simulations characterized by a homogeneous and low hydraulic conductivity (Table 1). In the area of the rift shoulders, this unit is confined between the topography and the lower model boundary. Its thickness reaches up to 9 km in the area of the Black Forest and the Vosges Mountains in the southern part of the rift shoulders (Fig. 3).
The model unit above the crystalline basement consists of Rotliegend, Buntsandstein, and Muschelkalk. Of these, the volcano-sedimentary successions of the Rotliegend are restricted to several SW-NE striking intramontane basins (GeORG-Projektteam 2013b). Following upward, the fractured continental sandstone aquifer of the Buntsandstein is distributed below most of the graben area (Stober and Bucher 2015). The marls and fossil-rich limestones of the Muschelkalk (Walter and Dorn 1995) cover the  Keuper and Jurassic (Lias/Dogger) sediments are overlying the Muschelkalk with a cumulative thickness decreasing toward the north, where they pinch out (Fig. 3). This is approximated in the model by a minimum thickness of 50 m. Of these, the Keuper deposits of alternating sandstones, clays, marls, and evaporites form an aquitard (Länderübergreifende Organisation für Grundwasserschutz am Rhein 2006). Accordingly, this unit is considered as low permeable.
Following upward, all Cenozoic sediments compose the uppermost aquifer. It is thicker in the north than in south in response to two separate locations of subsidence in Eocene times, of which the northern continued to subside until Pleistocene times (Reicherter et al. 2008). The Cenozoic fluvial and lacustrine deposits (Ellwanger et al. 2003) are the thickest sedimentary layer reaching a thickness of up to 3 km northward (Fig. 3).
The Cenozoic sediments of the graben are flanked by the uplifted rift shoulders, where the sedimentary cover was eroded and the crystalline basement as well as parts of the Mesozoic sediments were exposed. This asymmetric rift shoulder uplift due to the isostatic response to extension and subsidence (Illies and Greiner 2011) was weaker in the north than in the south and stronger in the east than in the west. Thus, the Black Forest is on average 200 m higher than the Vosges Mountains (Pflug 1982). The present 3D model is a geological, rectangular model of the central URG parallel to the main rift axis (Fig. 4). It extends about 150 km in N-S direction and about 90 km in E-W direction and to a depth of 8 km below mean sea  level. The structural model has an overall horizontal and vertical grid resolution of ≤ 1 km. Besides the above introduced model units, it encompasses the western and eastern main border faults. Each stratigraphic unit is parametrized as being isotropic and homogeneous as indicated in Table 1. Likewise, physical properties are assigned to the two major faults following Baujard et al. (2017), who investigated the properties of the major fault in Rittershofen.
In the model, a given stratigraphic layer can be found at very different depths due to post-and syn-depositional deformation during the rifting. In general, the faulted structure results in a complex hydrogeological situation within the graben, providing local hydraulic connections between the stratigraphic layers (Stober and Bucher 2015). Together with the generally high hydraulic conductivity of the sedimentary rocks in the Rhine valley, cross-formational flow of pore fluids may occur (Stober and Bucher 2015). Such local features are not resolved in the model.

Physical background
To simulate the fluid flow and evaluate the influence of the different boundary conditions, we used the commercial software package FEFLOW © (Diersch 2014). The governing (partial differential) equations for fluid flow in porous and fractured media are solved numerically in FEFLOW © by a multidimensional finite-element method (Diersch 2014). The solver in FEFLOW © respects the conservation laws of mass, energy, and linear momentum. We assume steady-state conditions, and solve for the following form of the fluid mass balance: Where q f is derived from the fluid momentum balance under laminar flow conditions, Darcy's law as: where K is the hydraulic conductivity tensor, p is the pressure, is the density of water, g is the gravitational acceleration, and z is the elevation of the water table. As we solve the problem under steady-state conditions, we are not accounting for variations of the fluid properties. For this reason, the effects of freshwater recharge and saltwater intrusion on density-driven fluid flow is not further discussed in this study.

Numerical model and boundary conditions
For this study, we needed a simulation software that allows to change the boundary conditions easily. With FEFLOW © , it is possible to set and change boundary conditions at any node in the model very efficiently. Moreover, we used FEFLOW © , to be consistent with the software framework as we used it to simulate the fluid flow of the entire URG (Koltzer et al. 2019), from which we used the pressure distribution to include the vertical inflow and outflow across the southern and northern model boundaries. To use FEFLOW © as modelling software, it was necessary to convert the structural model from Freymark et al. (2019), who simulated the coupled fluid and heat transport with the software GOLEM (Cacace and Jacquey 2017). We therefore first created a reference model using the same hydraulic settings applied by Freymark et al. (2019). Figure 4 shows the numerical model including the sedimentary layers, the Variscan domains, and the two main border faults. The latter are normal faults dipping with 55-85° toward the graben centre (Reicherter et al. 2008). In the model, the faults were prolonged downward to 6 km depth below sea level, and are considered to be 10 m wide and approximately 30 km apart from each other on average. The two main border faults are implemented as discrete features in FEFLOW © as geometric representations of lower spatial dimension (here a 2D fault in a 3D model) and can be assigned a fluid conductance different from the porous medium (Diersch 2014).
We used an unstructured finite-element mesh generated by Freymark et al. (2019) with the software MeshIt (Cacace and Blöcher 2015). Therefore, Delauney triangulation was applied to the faults and tetrahedrization (3D Delaunay triangulation) to the sedimentary layers and the crust. The mesh consists of 796260 elements and 155601 nodes.
With respect to the two main questions described in the introduction to this paper, to investigate the effects of different hydraulic boundary conditions on the deep fluid flow in the central URG, we compare the reference model (after For the two other model scenarios, the settings differ only in terms of the imposed boundary condition at the surface and at lateral northern and southern boundaries. In all model runs, the hydraulic boundary condition imposed along the base of the model and the western and eastern lateral boundaries (crystalline basement as no-flow boundary) were kept the same (Table 2). For the reference model, a pressure of 0 kPa is assigned along the topmost surface, that is, the resulting hydraulic head was considered equal to the topography, as applied by Freymark et al. (2019) (Fig. 5b).
In the first scenario (S1), we test a more realistic representation of the upper boundary condition. Therefore, we make use of interpolated groundwater surface pressures from 288 measured groundwater head data points: for the French part of the model provided by the Bureau de recherches géologiques et minières (2016) and for the    (Fig. 5a). Subsequently, the hydraulic head surface is cropped to the model area; thus, data points lying outside the model area are also taken into account to avoid boundary effects. Despite considering these information, the interpolation comes with some uncertainties. As data density and coverage are unevenly distributed, interpolation effects could be left in the final grid for areas devoid of data. Furthermore, a correction to the potentiometric surface due to the presence of surface water bodies (e.g., rivers) is not applied. This latter aspect is particularly relevant for the mountain valleys, where the main rivers are located, and water levels of up to -330 m above m.s.l. are encountered (Fig. 5c). To minimize this interpolation error, the water level in these specific domains is approximated to correspond to the topography (Fig. 5d) and was manually smoothed to avoid numerical artefacts. In contrast, in the Rhine lowland under wetlands (location 2 in Fig. 5c) and close to the rift shoulders (location 3 in Fig. 5c), the negative water table values, representing heads above land surface, can be assumed as realistic, because of the proximity to upland areas and presence of confined or semi-confined conditions (Ad-hoc Arbeitsgruppe Hydrogeologie 2016). S2 considers across boundary fluid flow along the southern and northern lateral boundaries of the model in addition to the newly derived water table surface. This way, the regional south-to-north flow component due to high hydraulic gradients along the Alps in the south compared to lower gradients in foreland in the north is integrated. For this purpose, the hydraulic head distributions with depth along the southern and northern vertical model boundaries are extracted from a larger model (Koltzer et al. (2019); location in Fig. 1) covering the entire URG from the Alps in the south to central Germany in the north. The simulations based on this larger model (Koltzer et al. 2019) addressed the effects of forced convective fluid flow, but do not consider the border faults as hydraulic pathways.
To transfer the hydraulic head distribution from the larger regional model to the current model along the northern and southern vertical boundaries, a workflow with four steps is adopted ("Appendix: four-step workflow"; Fig. 10).

Results
The results for the reference model are briefly summarized here; for more details, we refer to Freymark et al. (2019). Given the topographic elevation, the graben's rift shoulders represent the main recharge areas for meteoric water to infiltrate into the basement. From here, fluid flows either directly toward the border faults, where it is either channelled upward, or migrates to areas of lower hydraulic potential in the centre of the graben. Reaching the lowest hydraulic potential, groundwater originating from both rift shoulders merge and ascend to the surface. Within a 3D perspective (Fig. 6a), ascending fluids form a south-north-oriented upflow axis that superposes the perpendicular flow toward the graben centre. This upflow axis develops close to the eastern graben border in the central part of the model, where the cumulative thicknesses of the upper and lower aquifers are largest and the base of the graben sediments is deepest (up to 7000 m b.s.l.). In the northern and southern model area, where the sediment thickness is more evenly distributed, the upflow axis shifts to the west, approximately under the river Rhine in the centre of the graben (Fig. 6a).
Changing the hydraulic boundary conditions leads to differences in the simulated groundwater flow pattern, as exemplified in Fig. 6. The upflow axis is discontinuous and changes location (black stippled lines in Fig. 6b, c). Moreover, a detailed analysis of the results revealed that the flow direction both in the upper and lower aquifer at graben size, as well as the location and extent of secondary local flow systems, change. These results are important for other application studies in this area, where upflow of heated fluids may be of relevance (i.e., geothermal application sites) or for studies in 2D that include for example temperature/concentration effects and related density changes of the fluid.
By visually inspecting Fig. 6, the main variation among the three model realizations is found in the location and geometry of the preferential upflow axis (marked by the black lines in all three figures). In the reference model, this upflow axis is continuous and crosses the entire model area from south to north, being spatially located under the river Rhine (Fig. 6a). Both scenarios S1 and S2 are characterized by a more complex and discontinuous preferential upflow configuration in the central model area (stippled lines in Fig. 6b, c). Additionally, the upflow axis is shifted to the west by up to 7 km in the southern model area. In the northernmost model area, preferential upflow cannot be distinguished. These variations in predicted regional flow dynamics are the results of the different hydraulic boundary conditions at the top of the models. As we find no difference in the upflow location between S1 and S2, we can exclude that the open vertical southern and northern boundaries influence this process. Accordingly, the location of the predicted upflow axis or discharge area directly depends on the assigned upper hydraulic boundary condition.
The second aspect to point out is the difference in the direction of fluid streamlines between the three models, resulting from a different hydraulic pressure distribution. In the reference model, where the upper hydraulic boundary condition overestimates the hydraulic head, several local flow systems evolve forming subsidiary local scale recharge and discharge regions. In both S1 and S2, the streamlines in the crystalline basement show a more homogeneous trend than in the reference model, being mainly parallel. This effect is most pronounced along and across the graben shoulders, where the differences between the adopted boundary conditions are largest, though it can also be observed in the sediments at some distance from these domains (Fig. 6b, c).
The results of simulations S1 and S2 permit to investigate the fluid flow in the deeper aquifer system with respect to other controlling factors than the imposed hydraulic pressure of the topography. In response to the smoother hydraulic head surface adopted in S1 and S2, the system is less forced to follow the pressure-driven flow. This is particularly evident if comparing the upflow axis in of the reference model (coincides with the river Rhine) and S1 and S2 (more complex geometry of upflow axis).
For more detail, we show the topography (Fig. 3a), the thickness of the aquitard (Fig. 3d) and the top of the crystalline crust (Fig. 3g) in comparison to the upflow axis of the reference model (black lines) and S1/S2 (black dotted lines). While the upflow axis correlates with the path of the river Rhine (or with the topography) in the reference model, a spatial correlation between the geometry of the deep aquitard (Keuper/Lias/Dogger) is found in scenarios S1 and S2. In areas where the aquitard is thicker than 1 km (red in Fig. 3d), the fluid flow is bounded and the upflow axis is going around these thicker areas or even splitting up in their southern part.
The direction of modelled streamlines is also influenced by considering or neglecting fluid flow across the northern and southern boundaries. Closed vertical boundaries in the reference model and in S1 force the streamlines to be parallel at these boundaries. Open vertical boundaries lead to streamlines which are rather perpendicular to the boundary and thus allow a graben-parallel flow, which is not interrupted at the model boundaries.
The influence of open model boundaries can be tracked northward from the southern boundary and southward from the northern boundary. At the graben shoulders where the crystalline crust crops out, the effect of open or closed vertical boundaries on the simulated flow field can be traced to up to 10 km distance from the model boundaries. In the more permeable layers of the sedimentary graben fill, the effect reaches less far (up to 5 km, Fig. 6c). The influence of the boundary effect correlates spatially with the height of hydraulic gradients imposed at the upper boundary: boundary effects are largest in the Black Forest (southeast), medium in the Vosges Mountains (southwest), and smallest at the northern boundary.
Changing the perspective to a 3D view from the south reveals a consistency among all three models with respect to a regional flow pattern from the graben shoulders to the rift sediments (Fig. 7). In areas of maximum elevation, recharge Upflow axis in S1 and S2 c b a (R) develops, while neighbouring domains of lower elevations act as discharge areas (D). On a more local scale, the reference model is characterized by an increase in the number of secondary flow systems mainly at the rift shoulders and in the upper aquifer. This observation confirms the conclusion that even minimum differences in the topology of the adopted potentiometric surface have a first-order effect on the development of flow systems and the consequent evolution of recharge and discharge areas. The depth influence of these secondary flow systems in the graben sediments does not exceed the base of the first aquitard. This indicates that the shallower and deeper aquifers are hydraulically decoupled. In the deeper aquifer, the flow field has a more regional character being directed toward the rift centre in all three model scenarios. The hydraulic role of the aquitard disconnecting the deep from the shallow flow regimes is supported by the results obtained in areas where the aquitard is discontinuous. These areas correlate spatially with the main upflow axis, where both aquifers are hydraulically connected, thereby producing similar groundwater dynamics in all model realizations. Figure 8 shows the results of a quantitative analysis comparing the magnitudes of simulated Darcy velocities in the different models. Therefore, we plot the minimum, median, and maximum values as well as the 25 and 75 % quantile. The values are determined in the centre of each element and separately for each model unit. The results of this analysis demonstrate highest Darcy velocities in the reference model (black in Fig. 8), lower velocities in S2, and lowest velocities in S1 (blue and red, respectively, in Fig. 8). This is consistent  Fig. 1 . a The boundary condition of hydraulic head distribution in the reference model (black) and scenario S1 (red). The cross section is vertically 20 times exaggerated. b Simulated streamlines of the reference model and c of the scenario S1. R recharge area; D discharge area with the hydraulic pressure set with the upper and lateral boundary conditions that is highest in the reference model (equal to the topography) and lower in S1 and S2 (groundwater surface). The higher velocities in S2 compared to S1 are related to the imposed N-S gradient in S2 that adds a forcing to the flow that is missing in S1. The differences between the three models are however limited to less than one order of magnitude. Groundwater flow velocities are not directly measurable, but have been deduced from chemical data, mineral reactions, and numerical simulations. Comparing simulated Darcy velocities with published values highlights that simulations in this study underestimate the velocities independent of the applied boundary condition configuration. Aquilina et al. (2000) proposed fluid velocities of 1 m a -1 in the Buntsandstein based on 14 C data at the western border and Person and Garven (1992) predicted a maximum Darcy velocity of 0.13 m a -1 in the Buntsandstein/ Muschelkalk based on numerical models. In this study, the flow velocities in the Permo-Carboniferous/Buntsandstein/ Muschelkalk layer rarely reach meters or centimeters per year.

Discussion
The results of the three simulations have shown that specific processes in the URG are predicted independently of the chosen set of boundary conditions. Therefore, we can interpret those flow characteristics as generic and realistic features of the deep fluid flow in the URG. In contrast, varying the imposed boundary condition results in different flow behaviors caused either by the specific geometry of the upper boundary condition or by the hydraulic potential for specific areas, as especially is the case at the rift shoulders. Comparing our model results for all three scenarios, we are able to provide additional support for the major findings of the 3D numerical studies of Freymark et al. (2019): (1) a regional topography-driven fluid flow is confirmed where main recharge occurs at the rift shoulders. (2) Meteoritic waters infiltrating at both rift shoulders merge in the centre of the rift where they form a major north/northeast directed upflow axis, parallel to the graben axis. Consequently, the direction of flow strongly depends on the location within the rift system (Fig. 6a, b, and c). This upflow axis is the result of the applied hydraulic gradients due to pressure-driven advective flow and evolves even when lowering the absolute magnitude of the hydraulic pressure at the upper boundary (scenarios S1 and S2).
To assess if this flow pattern is also consistent for the case that hydraulic conductivity is anisotropic in the sediments, we run a test simulation for the reference model implementing an anisotropic hydraulic conductivity of sediments with ten times lower hydraulic conductivity in vertical direction than in the two horizontal directions (kx:ky:kz = 10:10:1). The results revealed that the upflow axis in the reference model is still evolving approximately under the river Rhine, but it is less pronounced and more discontinuous.
In general, our results show that the apparent contradictions of previous conceptual models can be reconciled if the full 3D flow field is considered. All the three simulations  Person & Garven (1992) predict vertical fluid ascent through the Buntsandstein and Muschelkalk (Stober and Bucher 2015) in concert with graben-parallel flow (Bächler et al. 2003) and confirm the main results of Freymark et al. (2019). At this stage, a further validation of the model scenarios presented against observations is difficult, since relevant information such as measured groundwater flow velocities are not available for the study area. However, the local upflow and downflow areas (expressed by the z-component of the Darcy velocity vector) of the simulated groundwater flow field can be compared to the locations of mineral and thermal springs and to artesian conditions (Fig. 6). A correlation between modelled upflow areas and observed discharge areas in the three models can be taken as a first qualitative reliability index. As already noted by Freymark et al. (2019), the reference model (topography equal to the groundwater surface) predicts the upflow zones very well, though not in all areas. In S1 and S2, there is no improvement in fitting the locations of springs and artesian conditions (Fig. 6). A possible reason for the remaining misfit of modelled and observed zones of rising fluids is the limited structural resolution of our model realizations. According to Loges et al. (2012), the thermal springs in the northern URG partly result from deep artesian aquifers intersected by faults, which might also be the case for other thermal springs in the URG. As faults and fractures in the URG apart from the main border faults are not differentiated in our current model realizations, fault-bound springs cannot be predicted (Fig. 9). Though implementing more structural detail still represents a numerical challenge, future work should aim on improving in this aspect.

Effects of the upper boundary condition on the flow systems
For the upper boundary condition, we tested two extreme settings. While the water table equal to the topography (reference model) leads to a larger number of spatially limited local flow systems below largest gradient differences, a water table as a subdued replica of the topography based on measured groundwater levels (S1 and S2) predicts more homogeneous flow lines. The higher frequency of hydraulic head variations in the reference model leads to the evolution of secondary flow systems with a small spatial extent in comparison to S1 and S2 (Fig. 7). These flow systems evolve primarily in the uppermost layer which is the upper aquifer in the Rhine valley and the basement at the rift shoulders. Even small local differences in the hydraulic potential imposed by the hydraulic boundary condition have an effect on the development of these shallow flow systems and, therefore, impact the evolution of predicted local recharge and discharge areas. In contrast, the lower aquifer is in the largest part of the URG disconnected from these shallow influences of topography as it is separated from the upper aquifer through the aquitard, which is in most of the model area around 500 m thick and locally thicker than 1 km (Fig. 3d). Only in the northernmost part of the model area, it is thinned out and, therefore, not decoupling upper and lower aquifer from each other. Because of this disconnection, the deep aquifer is not as much affected by local shallow topography-induced hydraulic head undulations as the upper aquifer. Conversely, we can say, the fluid flow in the deep aquifer is controlled by the regional component (from graben flanks toward graben centre and from south to north) and the hydraulic configuration at the top of the deep aquifer. There is one gap in the disconnection between deep and shallow aquifer, which is at the main upflow axis. The potentiometric surface realized by interpolation of available hydraulic head data shows lower hydraulic potentials at the rift shoulders. This, in turn, results in lower graben-perpendicular hydraulic gradients than in the reference model. The graben-parallel flow is not affected by the different upper boundary conditions as the depth to the water table is close to zero in the Rhine valley. As a main consequence, the graben-perpendicular flow becomes less vigorous in counteracting the graben-parallel flow component in S1 and S2. Without additional constraining local observations on groundwater movement, we cannot differentiate between the two upper boundary conditions concerning their reliability. Nevertheless, we consider that scenarios S1 and S2 to predict a more realistic hydraulic behaviour as the westward shift and a more heterogeneous main upflow axis in S1 and S2 than in the reference model are more consistent with the data base on springs and artesian conditions.
The strong sensitivity of the flow models with respect to the upper boundary condition underlines the great importance of appropriate constraints on the shallow hydraulic pressure variations realized in the upper hydraulic boundary condition.

Effects of open and closed vertical model borders
As noted by Koltzer et al. (2019), the central part of the URG cannot be considered as a hydraulically closed system. Indeed, the URG is characterized by inflow from the south and outflow to the north, whereas the influence of east-west flow across the rift shoulders is considered negligible. This is related to the fact that the main aquifers are located in the N-S-oriented sediment-filled graben with a hydraulic gradient from south to north, whereas the eastern and western vertical boundaries intersect almost impermeable basement rocks. In contrast, within the graben fill, the graben-perpendicular flow from the flanks to the centre is larger than the graben-parallel flow (Koltzer et al. 2019). Moreover, the top basement surface is characterized by a northwardinclined slope within the graben. This additionally favours north-directed, graben-parallel flow. Accordingly, different hydraulic boundary conditions at the southern and northern vertical model boundaries will influence the graben-parallel flow. A quantitative comparison of the effects induced by open versus closed vertical boundary conditions at the southern/northern model margins shows that such effects evolve only in the direct vicinity of the border regions (up to 10 km; Fig. 6c). Nevertheless, open vertical southern and northern model boundaries implemented in S2 lead to overall higher Darcy velocities between the model margins (Fig. 7).
In summary, the effect of open or closed vertical boundaries on the general fluid flow in the URG is rather small. The streamlines are shaped parallel to the model boundaries in an up to 10 km-wide stripe in the vicinity of the boundaries when those are considered to be closed (Fig. 6c). Apart from these portions next to the model boundaries, the boundary effects of closed lateral boundaries are negligibly small for our study, where the model size is 150 km in N-S direction. In contrast, for modelling studies on a more local scale, it can be crucial to consider appropriate hydraulic boundaries. Another point is that the overall predicted Darcy velocities are smaller with closed boundaries than with open boundaries (Fig. 9). This is an effect of the boundary, present in the whole model domain and not spatially restricted to the boundary nearfield. Moreover, boundary effects may be considerably larger if the slope of the aquifer base is steeper and/ or the difference in hydraulic gradient between the model boundaries is larger. Accordingly, the potential effects of lateral boundary conditions should always be evaluated for every specific hydrogeological setting.

Combined effects of upper and vertical boundary conditions
The superposed effects of a measured hydraulic head as upper hydraulic boundary condition and a hydraulic gradient between the southern and northern model boundaries can be evaluated by incorporating both in the model scenario S2. The main differences to the other scenarios consist in (1) a general increase of velocities in the entire model domain, (2) a shift of the central upflow axis, and (3) an overall decrease in the number of small local flow systems.
While the results of scenario S2 qualitatively appear more realistic than those of the other tested scenarios, a true validation remains a challenge. The reason for this lies in the lack of validation data compilations that consistently consider both surface (precipitation) and subsurface (deep recharge) effects. The recharge rates, which are rather representing the hydrological setting in the soil or at most down to the depth of the uppermost aquifer, are not representative for deep subsurface infiltration as predicted by our simulations. For areas where low permeable rocks are at the surface, high recharge in response to high precipitation is predicted that is in conflict with the real hydrogeology, where small permeabilities lead to small recharge rates. In addition, the role of the unsaturated zone is largely uncertain and resolving the processes related to the latter still requires simulation meshes of significantly higher resolution than those for fluid flow models of the saturated zone as different physical processes are involved.
In the Black Forest and the Vosges Mountains, high precipitation rates result in likewise high estimated groundwater recharge rates. In contrast, in the area west of the river Rhine, locally negative groundwater recharge values are reported (Länderübergreifende Organisation für Grundwasserschutz am Rhein 2013). In further modelling studies, a recharge data set applied as an upper flux boundary condition is conceivable. Here, average annual data from Länderübergreifende Organisation für Grundwasserschutz am Rhein (2013) for the French study part and from Kopp et al. (2018) for the German part of the model are available. However, a hydraulic conductivity of 10 × 10 −11 m s -1 at the graben shoulders will not allow for infiltration of meteoric water. This is inconsistent with relatively high predicted groundwater recharge rates for the graben shoulders (e.g., up to 400 mm a -1 recharge at the Black Forest after Kopp et al. 2018) which are an example for the above-mentioned problem. Alternatively, Herrmann (2010) provides a complex method to build regional groundwater pressure surfaces for the area of the Federal State of Hesse. The study of Goderniaux et al. (2013) in northern France proposes to use the recharge as function of hydraulic head. In this approach the correlation between the altitude of the hydraulic head with the amount of recharge is used and tested. Following this method that would mean for our simulations, that the topography as upper boundary condition for the hydraulic head (reference model) would represent the highest amount of recharge and the lower hydraulic head boundary condition in S1 and S2 could be compared to low recharge rates.

Model limitations
As the focus in this study is on quantifying the effects of different hydraulic boundary conditions, only hydraulic simulations are carried out. So far, only advective (pressuredriven) flow of a fluid with constant properties is considered, whereas possible effects of fluid density due to temperature and/or salinity are neglected. The effects of gypsum and halite dissolution, which are important for the southeastern adjoining area of the URG in Switzerland, as studied in Huggenberger et al. (2015) and Zidane et al. (2014aZidane et al. ( , 2014b are neglected in these model realizations. Moreover, in the model, the four geological layers and two faults are parameterized as being homogeneous and isotropic, which is a clear simplification. Effects on the flow regime related to additional (hydro-) geological and thermal complexity may superpose the effects related to the choice of hydraulic boundary conditions and need to be addressed separately in future studies.
Finally, using the proposed four-step method (see "Appendix: four-step workflow") also implies some shortcomings: (1) Larger regional models may not always be available. (2) In case of irregular model boundaries, it may be very timeconsuming to apply a vertical hydraulic boundary condition.
(3) Numerical instabilities may arise due to possible inconsistencies between the applied upper hydraulic boundary condition from groundwater measurements and the vertical hydraulic boundary condition extracted from a larger model. For the URG, the instabilities encountered in this context were larger at the southern than at the northern boundary. To avoid these instabilities completely, one suggestion for further modelling studies would be to use the same input data and methods for the preparation of the upper hydraulic boundary condition. Using the groundwater fluxes instead of extracting the hydraulic head distribution as a vertical hydraulic boundary condition could be a valid alternative to quantitatively consider the direction and magnitude of flow into and out of the model area.

Conclusions
In this conceptual study, we performed three 3D numerical simulations of fluid flow in the URG. By choosing different hydraulic upper and lateral boundary conditions, we gained valuable new insights on effects of the chosen boundary conditions on the simulated 3D hydraulic field of the URG.
-The general fluid flow pattern is not affected by the different chosen boundary conditions. It is characterized (i) by recharge areas at the rift shoulders, where descending waters flow through the graben fill perpendicular to the graben to rise along an axis parallel to graben in the central part, and (ii) a northward directed graben-parallel flow in the graben sediments. -The strong sensitivity of the flow models with respect to the upper boundary condition underlines the great importance of appropriate constraints on the shallow hydraulic pressure variations realized in the upper hydraulic boundary condition. -The specific hydraulic potential at the rift shoulders influences the direction of flow, the Darcy velocities, and the location of the upflow axis. -The geometry of the upper boundary condition influences the amount of locally restricted flow systems, even minimal changes in the groundwater pressure surface are causing differences.
-For other model areas worldwide, the decision to use a replica of the topography or hydraulic head data, closed or open model boundaries would depend on the (hydro-) geological situation and the size of the model area, modelling aim, and the amount of available hydraulic head data. A case-specific evaluation is recommended. For this case study, we can conclude that the simulation results are significantly more sensitive to the upper than to the vertical boundary conditions.

Appendix: four step workflow
We established a four-step method to transfer the hydraulic head distribution along vertical model boundaries from a larger regional model to this model (Fig. 10). In a first step, the results (pressure with depth) of the steady-state simulations from Koltzer et al. (2019) were imported to the open source software Paraview (Squillacote et al. 2007). In a second step, the vertical slices within the supra-regional model at the location of the northern and southern border of the model domain were generated with Paraview. These slices were further cropped to the size of the model domain borders. The pressure values along the grid nodes of the two slices were exported and recalculated to hydraulic heads (h) according to Eq. 3 (step 3) Finally, the hydraulic head values have been applied as a Dirichlet boundary condition of constant hydraulic head along vertical borders in FEFLOW © (step 4).
Data availability statement All data to reproduce the modelling results of all three described models are published in GFZ Data Services (Koltzer et al. 2021