Scientific Leader: Bibi Sarwat Naz
In hydropower and geothermal applications, modeling of shallow subsurface flow is of major importance in order to accurately simulate and predict the exchange of groundwater with streams under low-flow conditions, and the transport of energy. The major challenge is the representation of topographically driven groundwater convergence and streamflow generation, and of the geological heterogeneity across a number of space scales ranging from centimeters to thousands of kilometers in case of continental river systems. Constructing hydrological and geothermal models at this resolution over large spatial scales for scientific and operational applications constitutes a game changer, easily reaching up to 1012 degrees of freedom, where simulations must additionally assimilate observations to mitigate uncertainties in model data.
In EoCoE-II, the integration of hyper-resolved simulation of hydrological fluxes, routing along the river network, and management of storage reservoirs will be performed with a modernized version of ParFlow with adaptive mesh refinement (AMR) capability coupled with HYPERstreamHS model. The added values of these simulations will be shown by feeding ParFlow gridded runoff time series into the operational hydropower model embedded into HYPERstreamHS. The coupling will be specifically developed over the Italian Alpine region.
Previous work on geothermal reservoir characterization showed the successful application of optimal experimental design (OED) within the simulation code SHEMAT-Suite in order to identify optimal drilling locations for assessing uncertain reservoir parameters within a numerical reservoir model. However, the high computational cost has to date limited this approach to a numerical model with significantly reduced number of unknowns. Collaboration with experts in the EoCoE-II consortium will enable us to create a realistic geothermal reservoir model with vastly improved spatial resolution. Combining optimal experimental design for positioning boreholes with state-of-the art HPC techniques will improve the exploration and exploitation of geothermal reservoir systems, as it enables a sophisticated quantification of uncertainties in the subsurfac
ParFlow (https://parflow.org/) is a physics-based 3D parallel hydrologic model, which simulates surface and 3D subsurface flow, based on Richard’s and kinematic wave equations on a finite difference, finite volume grid with a preconditioned, implicit Newton-Krylov solver for the nonlinear PDEs. ParFlow will not be completely rewritten in EoCoE-II, but some code upscaling is planned as well as the activation of the Adaptive Mesh Refinement of the computing grids with the ad hoc numerical scheme and corresponding solver. It is also planned to have a strong effort on IO that will be coordinated with the effort on ensemble runs.
SHEMAT-Suite is a code for simulating single- or multi-phase heat and mass transport in porous media. It solves coupled problems comprising heat transfer, fluid flow, and species transport. SHEMAT-Suite can be applied to a range of hydrothermal or hydrogeological problems, be it forward or inverse problems.
HYPERstreamHS is distributed hydrological model based on the width function
instantaneous unit hydrograph (WFIUH) theory, which is specifically designed to
facilitate coupling with gridded climate datasets and climate model outputs.
HYPERstreamHS inherits the computational grid from the overlaying
meteorological forcing, still preserving geomorphological dispersion caused by
the river network irrespectively of the grid resolution. In addition
HYPERstreamHS is capable of simulating explicitly water transfers due to human
activities (e.g., diversion channels and storage reservoirs).
Pan-European spatially and temporally consistent 3 km soil moisture reanalysis using land surface data assimilation
fluxes between the land surface and atmosphere (Seneviratne et al., 2010). As a result, it plays an important role in many regional-scale applications, including meteorology, hydrology, flood forecasting, drought monitoring,
agriculture and climate change impact studies (Brocca et al., 2010). Because of its high spatiotemporal variability, it is difficult to monitor soil moisture at large spatial scales and remains the most difficult variable to obtain because there are no high-resolution soil moisture observations available at the continental scale. To produce pan-European spatially and temporally continuous information of soil moisture at 3 km resolution, the land surface data assimilation system CLM-PDAF
consisting of the Community Land Model (CLM; Oleson et al., 2008) and the Parallel Data Assimilation Framework (PDAF; Nerger and Hiller, 2013) coupled CLM-PDAF framework (Kurtz et al., 2016) was used. We implemented CLM-PDAF to
assimilate the satellite based soil moisture dataset ESA CCI (the European Space Agency Climate Change Initiative; Wagner et al., 2012) into CLM, producing a high-resolution European SSM reanalysis (called ESSMRA hereafter) dataset.
This product overcome the shortcomings of sparse spatial and temporal distributions in observations and provides a better estimate of SM than obtained only by modeling or by sparse observations alone. To evaluate the quality of ESSMRA, we compared the daily SM from ESSMRA against in-situ data using the SM measurement acquired from the International Soil Moisture Network (ISMN; Dorigo et al., 2011) across Europe. The ESSMRA dataset shows overall good agreement with in-situ observations at daily time scale as shown in Figure 2 for two networks REMEDHUS and SMOSMANIA located in Spain and France, respectively. This comparison shows that at daily scale the model is able to reproduce the daily variations in soil moisture fairly well, with overall correlation above 0.60 and RMSE ranges 0.04 to 0.19.
Overall, multiple validations revealed that ESSMRA data were consistent with the in-situ soil moisture data and other existing global reanalysis products. The relatively long time series and fine spatial resolution of this new European gridded ESSMRA dataset could provide a valuable data source for many hydrological applications. For example, it can be used as an initial input data for climate change analysis and for numerical weather prediction modelling to improve the model forecast in terms of location and amount of extreme precipitation events. This dataset will be also useful to understand the development and persistence of extreme weather events such as droughts, floods and heatwaves.
Experimental Design for Geothermal Modeling
Drilling boreholes during exploration and development of geothermal reservoirs not only involves high cost, but also bears significant risks of failure. In geothermal reservoir engineering, techniques of optimal experimental design (OED) can improve the decision making process. For instance, during exploration and production of geothermal fields, OED can be used to place additional slim holes, thus decreasing drilling costs and risk. Previous publications explained the formulation and implementation of this mathematical optimization problem and demonstrated its feasibility for finding borehole locations in two- and three-dimensional reservoir models that minimize the uncertainty of estimating hydraulic permeability of a model unit from temperature measurements (Seidler et al. 2014; Seidler et al. 2016). Subsequently, minimizing the uncertainty of the parameter estimation results in a more reliable parametrization of the reservoir simulation, improving the overall process in geothermal reservoir engineering.
OED is a mathematical optimization method. The general approach is to find optimalexperimental conditions for constraining model parameters. In other words, OED gives answers to the question: How do I have to design an experiment in order to collect data from which I can predict model parameters with least uncertainty? This is a question of the sensitivity of the model with respect to the unknown parameters. This sensitivity is mathematically described by the Fisher Matrix. It contains the first order derivative of the model output towards the parameters. For evaluating the information contained in the Fisher
matrix, OED criteria are formulated. One is the D-optimal design criterion which is based on the determinant of the Fisher Matrix (other criteria are based on the trace or eigenvalue). The minimum of the D-optimal criterion contains the maximum sensitivity.
Various OED techniques are implemented in the Environment for Combining Optimization and Simulation Software (EFCOSS) (Seidler et al. 2014). This software framework links mathematical optimization software with SHEMAT-Suite, our geothermal simulation code for fluid flow and heat transport through porous media, for addressing problems arising from geothermal modeling. Within EoCoE-II we will extend this OED approach with further functionalities for geothermal reservoir modeling and aim at increasing its performance in order to apply it to reservoir scale production runs. As a first step, synthetic test models have been set up for studying certain aspects, such as sensitivity of the OED result to a priori data.
Synthetic OED study: 2D geothermal reservoir
model above salt diapir[JF1]
Here we show first results of an OED study for a 2-dimensional synthetic model of a geothermal reservoir above a salt diapir (yellow unit), which is cross-cut by two high permeable faults (blue and red units) (modified after Rath et al. 2006). Each model unit is characterized by constant thermal and hydraulic properties (e.g., permeability, porosity, thermal conductivity). High heat conductivity of salt results in higher temperatures at shallower depth. This makes the sedimentary units above and close to salt diapirs interesting for direct use of geothermal energy. In addition, the permeable faults in our model provide pathways for advective heat transport, resulting in a heat transport towards the surface through the western fault (unit 11, blue). There is a borehole at location x0=7475 m with a certain depth z0=1475 m providing a temperature log Tlog0. The problem is to find the location x1 for an additional borehole of the same depth for measuring another temperature log Tlog1. These temperature data shall be used for estimating the fault permeability k10 and k11 with least uncertainty (i.e., optimally). The experimental condition is the horizontal position of the additional borehole. We assumed an a priori fault permeability k10=k11=5*10-15 m². The synthetic true fault permeability is 5*10-14 m². For this OED problem, the numerical simulation provides two possible optimal horizontal ranges for placing the additional borehole (Figure 1, top) : (i) from 1025 m to 1175 m and (ii) from 8225 m to 9975 m.
We will further use this model in the course of the project for simulating various OED problems and studying the impact of factors such as a priori assumptions, location of a priori data logs and their data quality.
Figure 2: Numerical
Figure 3: Numerical forward model of a synthetic geothermal reservoir above a salt diapir computed with SHEMAT-Suite Top: Reservoir structure in terms of geological model units, unit 9 in yellow is the salt and units 10 and 11 in red and blue are permeable faults. The remaining units are various sedimentary layers. X0 marks the location of the existing borehole and the black crosses depict the OED result in terms of normalized and binned D-optimality. Red arrows display two optimal ranges for an additional borehole. Middle: Darcy flow in terms of hydraulic reference head and darcy velocity (arrows) for the true reservoir properties. Bottom: Steady state temperature distribution for the true reservoir properties. Image credits: Johanna Bruckmann, M.Sc., E.ON Energy Research Center, RWTH Aachen University.
Impacts of groundwater representation on climate simulations
High-resolution large-scale predictions of hydrologic states and fluxes are important for many regional-scale applications and water resource management. However, because of uncertainties related to forcing data, model structural errors arising from simplified representations of hydrological processes or uncertain model parameters, model simulations remain uncertain. To quantify this uncertainty, multi-model simulations were performed at 3km resolution over the European continent using the Community Land Model (CLM3.5) and the ParFlow hydrologic model. While Parflow uses a similar approach as CLM in simulating the snow, vegetation and land-atmosphere exchange processes, it simulates three-dimensional variably saturated groundwater flow solving Richards equation and overland flow with a two-dimensional kinematic wave approximation. Both models were driven with the COSMO-REA6 reanalysis dataset at 6km resolution for the time period from 2000 to 2006 at an hourly time step, and both used the same datasets for the static input variables (such as topography, vegetation and soil properties). The performance of both models was analyzed through comparisons with independent observations including satellite-derived and in-situ soil moisture, evapotranspiration, river discharge, water table depth and total water storage datasets. GRACE satellite for total water storage (TWS) comparison over PRUDENCE regions (regions in black color in Figure 3) were used. The overall seasonal variation in TWS simulated by both models captured well when compared with the GRACE data (Figure 3b), However, both models show much stronger negative anomalies in summer over most regions. These simulations help to quantifying uncertainty associated with the representation of hydrological processes, such as groundwater flow and subsurface flow and its control on latent and sensible heat fluxes at the surface.

Geothermal potential for direct heating of a city quarter in the Lower-Rhine Embayment:
During the first phase of EoCoE researchers at the Institute for Applied Geophysics and Geothermal Energy (RWTH Aachen University) assessed the geothermal energy potential for a city quarter situated in the Lower Rhine Embayment. They had the opportunity to interact with stakeholders from the city of Geilenkirchen who are planning the reconstruction of a former NATO settlement next to the NATO-Airbase in Geilenkirchen-Teveren, north of Aachen (Germany). This way the simulations represent realistic scenarios for district heating and results may have a direct application in the renewed city quarter. Additionally, the study area in the Lower Rhine Embayment reveals results of more general interest for geothermal applications in sedimentary basins facing medium to high groundwater flow rates.
The integrated approach comprised geological modeling, numerical flow and heat transport model calibration, uncertainty analysis of flow model parameters and borehole heat exchanger (BHE) field simulations. The objective was to study BHE field performance in terms of long-term operation and efficiency, as well as its thermal impact on the ground.

The 3D coupled numerical flow and heat transport code SHEMAT-Suite includes a finite differences approach for calculating heat transfer processes within the borehole heat exchangers (Mottaghy & Dijkshoorn 2012). This allows for considering heterogeneous geological models and the influence of groundwater flow (i.e. advective heat transport) on the BHE performance. The approach is computationally feasible for whole BHE fields with hundreds of single BHEs and long-term simulations without the need for simplifications in operation time or temporal resolution of the BHEs. Additionally, code developments towards High-Performance Computing (HPC) and the use of modern HPC architectures enable simulating larger model areas or finer discretization in order to resolve heterogeneous or complex model geometries in more complex geological study areas.
More detailed information on the study results are presented in this poster: Infoposter_EoCoE Simulation of Borehole Heat Exchanger Fields based on an integrated model approach (open in a new page.
Naz, B. S., Kurtz, W., Montzka, C., Sharples, W., Goergen, K., Keune, J., Gao, H., Springer, A., Hendricks Franssen, H.-J., and Kollet, S. in Hydrol. Earth Syst. Sci., 23, 277-301, 2019.
Studying the influence of groundwater representations on land surface-atmosphere feedbacks during the European heat wave in 2003
Keune, J., F. Gasper, K. Goergen, A. Hense, P. Shrestha, M. Sulis and S. Kollet (2016) in J. of Geophys. Res. Atmos.. 121(13), 301-13,325.
A scale consistent terrestrial systems modeling platform based on COSMO, CLM, and ParFlow.
Shrestha, P., Sulis, M., Masbou, M., Dollet, S. And Simmer, C. (2014) in Monthly Weather Review, 142, 3466-3483.
Capturing the influence of groundwater dynamics on land surface processes using an integrated, distributed watershed model.
Kollet, S.J. and Maxwell, R.M. (2008) in Water Resources Research 44: W02402.
Proof-of-concept of regional scale hydrologic simulations at hydrologic resolution utilizing massively parallel computer resources.
Kollet, S.J., Maxwell, R.M., Woodward, C.S., Smith, S.G., Vanderborght, J., Vereecken, H., and Simmer, C. (2010) in Water Resources Research, 46, W04201.
Efficient Solution Techniques for Multi-phase Flow in Porous Media.
Büsing H (2018) in : Ivan Lirkov, Svetozar Margenov (Eds.), Large-Scale Scientific Computing : 11th International Conference, LSSC 2017: Sozopol, Bulgaria, June 5-9, 2017 : Revised Selected Papers
Verbundprojekt MeProRisk-II: Optimierungsstrategien und Risikoanalyse für tiefe geothermische Reservoire –Eine Machbarkeitsstudie.
Marquart G (2015) Teilprojekt A: Numerische Modellbildung und Optimierungsverfahren. Endbericht BMWi-Projekt 0325389A.
Implementing an effective finite difference formulation for borehole heat exchangers into a heat and mass transport code.
Mottaghy D, Dijkshoorn L (2012) in Renewable Energy, 45 :59-71.
Joint three-dimensional inversion of coupled groundwater flow and heat transfer based on automatic differentiation: sensitivity calculation, verification, and synthetic examples.
Rath V, Wolf A, Bücker HM (2006) in Geophys. J. Int., 167:453-466.
Optimal experimental design for reservoir property estimates in geothermal exploration. Comput Geosci, 20:375-383.
Seidler R, PadalkinaK, Bücker HM, Ebigbo A, HertyM, Marquart G, Niederau N (2016) in Computational Geosciences, 20:375-383.
On the Design of the EFCOSS Software Architecture When Using Parallel and Distributed Computing
Seidler R, Bücker H, Rostami M, Neuhäuser D (2014) in Proceedings of the 9th International Conference on Software Engineering and Applications (ICSOFT-EA-2014), pp 445-454
Ensemble-based stochastic
permeability and flow simulation of a sparsely sampled hard-rock aquifer.
Bruckmann J, Clauser C (2020) in Hydrogeology
Journal, 28(5): 1853-1869, doi: 10.1007/s10040-020-02163-5.
SHEMAT-Suite: An open-source code
for simulating flow, heat and species transport in porous media.
Keller J, Rath V, Bruckmann J, Mottaghy D,
Clauser C, Wolf A, Seidler R, Bücker H M, Klitzsch N (2020) in SoftwareX, 12:100533, doi: 10.1016/j.softx.2020.100533.
Work in progress