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 hydrologic and geothermal models at this resolution over large spatial scales for scientific and operational applications constitutes a game changer, easily reaching up to 10^{12} 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. The added values of these simulations will be shown by feeding ParFlow gridded runoff time series into the operational hydropower model that 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 subsurface.

**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.

**ExaTerr** is a new development that aims at building a common software platform for both SHEMAT-Suite and ParFlow. This platform will be based on Kokkos, a software technology strongly pushed by the US DoE which holds in his heart performance portability. The idea pursued here is to have a single code base that could be executed on both CPU- and GPU-based platforms.

### Success story 1:

**Continental-scale high-resolution land surface data assimilation system**

The land surface data assimilation system CLM-PDAF (consisting of the Community Land Model (CLM 3.5) (Oleson et al., 2004) and the Parallel Data Assimilation Framework (PDAF)(Kurtz et al., 2016) was used to update the soil moisture estimates from the land surface model utilizing the coarse resolution satellite soil moisture data. A key capability of CLM-PDAF is the support for data assimilation that combines land surface processes with satellite and in-situ observations for the estimation of optimal land surface states. The data assimilation structure in CLM-PDAF allows to directly ingest remotely sensed observations of land surface conditions to produce accurate, spatially and temporally consistent fields of land surface states, with reduced associated error. The CLM-PDAF uses an Ensemble Kalman Filter (EnKF) algorithm (Burgers et al., 1998) to generate assimilated or reanalysis products. To effectively simulate the background-error covariances, a large enough ensemble size needs to be maintained in the data assimilation process, which linearly increases the computational resource requirements. The CLM-PDAF is designed for high-performance computing infrastructures and can efficiently cope with the high computational burden of ensemble-based data assimilation.

Hence, we implemented the CLM-PDAF over Europe to provide downscaled estimates of the soil moisture with complete spatiotemporal coverage by combining historical satellite SM observations with a high resolution LSM using data assimilation techniques. Using the CLM-PDAF, the satellite based soil moisture dataset ESA CCI (the European Space Agency Climate Change Initiative (Wagner et al., 2012) was assimilated into CLM using the EnKF producing a high-resolution European SSM reanalysis (called ESSMRA hereafter) dataset. This product overcomes the shortcomings of sparse spatial and temporal datasets and provides a better estimate of SM than obtained only by modeling or by sparse observations alone. The 3 km ESSMRA is generated by first implementing the regional land surface model setup coupled with the data assimilation framework as shown in Figure 1 (a). In the second step the ESA CCI satellite-based data is assimilated into the CLM-PDAF setup to generate the daily 3km ESSMRA product over Europe for the 2000 – 2015 time period Figure 1 (b).

References:

Burgers, G., Leeuwen, P. J. v., and Evensen, G. Analysis scheme in the ensemble Kalman filter , Monthly weather review, 126 , no. 6, 1719–1724, 1998.

Kurtz, W. , He, G., Kollet, S.J., Maxwell, R.M. , Vereecken, H. and Franssen, H.J. H., TerrSysMPPDAF (version 1.0): a modular high-performance dataassimilation framework for an integrated land surfacesubsurface model, Geoscien-tific Model Development,9, no. 4, 1341–1360, 2016.

Oleson, K., Dai, Y., Bonan, G. B., Bosilovichm, M., Dickinson, R., Dirmeyer, P., … Zeng, X. (2004). Technical Description of the Community Land Model (CLM) (No. NCAR/TN-461+STR). University Corporation for Atmospheric Research. doi:10.5065/D6N877R0.

Wagner, W. , Dorigo, W. , de Jeu, R. , Fernandez, D. , Benveniste, J. , Haas, E., and Ertl, M. Fusion of active and passive microwave observations to create an essential climatevariable data record on soil moisture, ISPRS Annals of the Photogrammetry, RemoteSensing and Spatial Information Sciences (ISPRS Annals),7, 315–321, 2012.

### Success story 2:

**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 optimal experimental 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**

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 x_{0}=7475 m with a certain depth z_{0}=1475 m providing a temperature log Tlog_{0}. The problem is to find the location x_{1} for an additional borehole of the same depth for measuring another temperature log Tlog_{1}. These temperature data shall be used for estimating the fault permeability k_{10} and k_{11} with least uncertainty (i.e., optimally). The experimental condition is the horizontal position of the additional borehole. We assumed an a priori fault permeability k_{10}=k_{11}=5*10^{-15} m^{2}. The synthetic true fault permeability is 5*10^{-14} m^{2}. For this OED problem, the numerical simulation provides two possible optimal horizontal ranges for placing the additional borehole (Figure 2, 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.

### Success story 3

**Impacts of groundwater representation on climate simulations**

By coupling ParFlow to land and meteorological models, the fully coupled water cycle from groundwater to atmosphere can be simulated, a novel exploration in that atmospheric models rarely incorporate a dynamic water table and three dimensional subsurface flow. This study couples ParFlow to the meteorological model over Europe during the 2003 heat wave in order to investigate the effects of various lower boundary conditions and configurations on land-atmosphere moisture exchange and thermal energy (Keune et al., 2016).

Improved soil moisture estimates by assimilating satellite observations: This study investigates the value of assimilating coarse-resolution remotely sensed soil moisture data into high-resolution land surface models for improving soil moisture and runoff modeling. The soil moisture estimates in this study, with complete spatio-temporal coverage and improved spatial resolution from the assimilation, offer a new reanalysis product for the monitoring of surface soil water content and other hydrological fluxes at 3 km resolution over Europe (Naz et al., 2019).

**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).

D4.2 – Analysis of temperature ranges and fluid flow rates in a 3D subsurface model

D4.3 – Layout of an energy supply system for Teveren based on geothermal direct-heat

Improving soil moisture and runoff simulations at 3 km over Europe using land surface data assimilation

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