Search

Bayesian inference has become the most popular approach to uncertainty assessment in vadose zone hydrological modeling. By combining prior information with observations and model predictions, it became popular among hydrologists as it enables them to infer parameter posterior distributions, verify model adequacy, and assess the model's predictive uncertainty. In particular, the posterior distribution is frequently the variable of interest for modelers as it describes the epistemic uncertainty of model parameters conditioned on measurements. Gradient-free Markov-Chain Monte Carlo (MCMC) ensemble samplers based on Differential Evolution (DE) or Affine Invariant (AI) strategies have been used to approximate the posterior distribution, which is frequently anisotropic and correlated in vadose zone-related problems. However, a rigorous benchmark of different MCMC algorithms to provide guidelines for their application in vadose zone hydrological model calibration is still missing. In this study, we elucidate the behavior of MCMC ensemble samplers by performing an in-depth comparison of four samplers that use AI moves or DE-based strategies to approximate the target density. Two Rosenbrock distributions, and one synthetic and one actual case study focusing on the inverse estimation of soil hydraulic parameters using HYDRUS-1D, are used to compare algorithms in different dimensions. The analysis reveals that AI-based samplers are immune to affine transformations of the target density, which instead double the autocorrelation time for DE-based samplers. This behavior is reiterated in the synthetic scenario, for which AI-based algorithms outperform DE-based strategies. However, this performance gain disappears when the number of soil parameters increases from 7 to 16, with both samplers exhibiting poor acceptance rates, which are not improved by increasing the number of chains from 50 to 200 or by mixing different strategies.

The Wairau Plain in Marlborough, New Zealand, host a shallow, highly conductive gravel aquifer which is mainly recharged by the braided Wairau River. The groundwater of the Wairau Aquifer is an important source of drinking water in the region and provides irrigation water for viticulture. Managing the groundwater is increasingly challenging because of a decline in levels and spring flows combined with a natural seasonality and an expected high vulnerability to climate change and overexploitation.A prerequisite for groundwater managment and limit-setting is the knowledge of the water supply (recharge) as well as the water discharge and takes from the aquifer. Both are spatialy and temporally highly variable and difficult to measure. In addition, the aquifer is influenced by ephemeral streams which could potentially lead to unknown groundwater sinks and sources and periodical shifts of boundary locations.In order to investigate these aspects and to estimate the water balance components and fluxes in the associated groundwater system, a coupled surface water – groundwater flow model (MODFLOW) has been set up. A variety of observations, ranging from groundwater heads to spring/river flows and exchange rates, were built into the model. Metered groundwater abstraction data was used to test and adapt a spatio-temporally distributed soil water balance model that is coupled to the MODFLOW model.The complex surface flow network of the Wairau Plain and strong contrasts of hydraulic conductivity (up to 4 orders of magnitude) between older hydrogeological units and more recent, interwoven fluvial sediments have been a major challenge in the model setup. An iterative procedure with step-wise increasing complexity of parameterization was adopted to derive a plausible model structure, that is commensurable with the various observation types. Model parameters were calibrated and (linear) uncertainty bounds estimated using the PEST software.The model performs well with respect to the plausibility of the groundwater flow field and the overall water balance. Groundwater heads, river-groundwater exchange flows and spring flows are generally well covered by the predictive uncertainty bounds during the evaluation period (data not utilized for calibration). The model provides insights into the relative contributions and the seasonality of the various water balance components of the Wairau Aquifer. It has been confirmed that the Wairau River is the main contributor to groundwater recharge, but also that recharge in a given year is matched by equal levels of discharge to low-land springs and off-shore. Although large river flood events during the wet period lead to interim excess groundwater storage, the recession to pre-flood groundwater levels is surprisingly fast (less than one year). On the other hand, unremarkable (high return period) flood events in the summer period, which keep river flows at elevated levels, have a noticeable effect on groundwater storage. This effect can last up until the following dry season.The model proved useful for assessing the status quo of the Wairau Aquifer groundwater resources. This is a prerequisite for the search and testing of alternative management options that are imminent due to the observed trends in groundwater levels.

Braided rivers are an important source of groundwater recharge in New Zealand. They consist of multiple temporary channels in a gravel environment and, as a consequence, their interactions with groundwater are complex and highly variable in space and time at different scales. Recently, the gravels of the contemporary braidplain of these rivers have been described and referred to as the ‘braidplain aquifer’. It is within this aquifer that hyporheic and parafluvial flows occur. In these systems, the groundwater recharge to the deeper regional aquifer is actually the water exchange between the braidplain and the regional aquifers. Some of these braided rivers are perched above the regional water level in their main losing section, which means that an unsaturated zone exists between the braidplain and regional aquifer. This complexity calls for the use of 3D fully integrated hydrological models to represent groundwater – surface water interactions in these environments. However, these complex models are very computationally intensive, which strongly limits their use in parameter inference and uncertainty quantification schemes as well as their applicability to regional scale problems.We present a modelling framework that includes a 3D fully coupled HydroGeoSphere (HGS) model and several 2D cross-sectional HYDRUS-2D models (with 1, 2 and 3 layers). This framework aims at simplifying the model while ensuring the appropriate simulation of the groundwater recharge. We demonstrate our modelling approach on the relatively small Selwyn River. Piezometric data and groundwater recharge estimates derived from satellite photography were available for this river. First, stochastic simulations were performed using the 2D cross-sectional models and compared to observations in order to explore the validity of different subsurface conceptualizations and parameter values. Second, the selected conceptualization and parameter values were used to parameterize the 3D fully coupled HGS model. Third, the groundwater recharge simulated by the 3D and the 2D models were compared. Our results demonstrate that the observations can only be reproduced with a minimum of 3 distinct layers, with a lower permeability layer in the middle. Furthermore, this modelling exercise revealed the primary importance of the width and thickness of the braidplain aquifer as they determine the infiltration front width and the pressure head applied to the braidplain aquifer bottom, respectively. This shows that the properties, dimensions and water level in the subsurface are controlling the groundwater recharge from the perched braided river rather than the river characteristics. Moreover, we show that a 2D cross-sectional model can effectively replace the 3D fully coupled model to simulate groundwater recharge from the perched braided river and that this reduces the model run time by 3 orders of magnitude. Finally, some analytical equations, which can be easily implemented in regional groundwater models, were tested as a further simplification of the 2D model.

With the increasing threat of climate change, projections of its impact on the availability of natural resources, such as groundwater, are becoming significantly more important. While extrapolation or prediction is a valid application for environmental system models, their assumptions and limitations are often not clearly communicated or investigated. This is especially true for data-driven models, which have been applied more frequently, and often with great success, to groundwater problems in the last decade. But are these techniques applicable to long-term future predictions of the impact of climate scenarios on groundwater resources, and if so, under which conditions and limitations?Within the context of KlimaKonform (Technische Universität Dresden, 2023), a research project studying the impact of climate change on the Central German Uplands, ensembles of multi-layer perceptron (MLP) models have been derived to estimate groundwater levels for a variety of wells in a region in Saxony-Anhalt in Germany. Once trained to replicate historical time series with climatological input data such as precipitation and temperature, these model ensembles were then tasked to predict the groundwater levels up until the end of the current century under different climate scenarios.We analyse the plausibility of the model ensemble predictions to shed light on the above-proposed question: what are factors (and metrics) of success, as well as limitations and possible failures, of data-driven methods (MLPs in this case) for long-term prediction? First, we propose that using ensembles of models, rather than “the single-best” trained model, is a necessity for such applications. We identify different methods of pre-processing of input and target data, structural model setup as well as target-oriented post-processing of ensemble simulations as vital factors for the successful application of MLPs to long-term prediction of groundwater levels under climate change scenarios. We further highlight remaining limitations and pose the question of how they could potentially be overcome. Technische Universität Dresden, 2023. KlimaKonform: Forschungsprojekt KlimaKonform. https://klimakonform.uw.tu-dresden.de/

Characterization of karst systems and forecast of their state variables are essential for groundwater management and engineering in karst regions. These objectives can be met by the use of process-based discrete-continuum models (DCMs). However, results of DCMs may suffer from inversion nonuniqueness. It has been demonstrated that the joint inversion of observations regulated by different natural processes can tackle the nonuniqueness issue in groundwater modeling. However, this has not been tested for DCMs thus far. This research proposes a methodology for the joint inversion of hydro-thermo-chemo-graphs, applying to two small-scale sink-to-spring experiments at Freiheit Spring, Minnesota, USA. In order to address conceptual uncertainty, a multimodel approach was implemented, featuring seven mutually exclusive variants. Spring hydro-thermo-chemo-graphs, for all the variants simulated by MODFLOW-CFPv2, were jointly inverted using a weighted least squares algorithm. Subsequently, models were compared in terms of inversion and forecast performances, as well as parameter uncertainties. Results reveal the suitability of the DCM approach for simultaneous inversion and forecast of hydro-physico-chemical behavior of karst systems, even at a scale of meters and seconds. The estimated volume of the tracer conduit passage ranges from approximately 46–51 m3, which is comparable to the estimate from the flood-pulse method. Moreover, it was demonstrated that the thermograph and hydrograph contain more information about aquifer characteristics than the chemograph. However, this finding can be site-specific and should depend on the analysis scale, the considered conceptual models, and the hydrological state, which are potentially affected by minor unaccountable processes and features.

Transmission losses are the loss in the flow volume of a river as water moves downstream. These losses provide crucial ecosystem services, particularly in ephemeral and intermittent river systems. Transmission losses can be quantified at many scales using different measurement techniques. One of the most common methods is differential gauging of river flow at two locations. An alternative method for non-perennial rivers is to replace the downstream gauging location by visual assessments of the wetted river length on satellite images. The transmission losses are then calculated as the flow gauged at the upstream location divided by the wetted river length. We used this approach to estimate the transmission losses in the Selwyn River (Canterbury, New Zealand) using 147 satellite images collected between March 2020 and May 2021. The location of the river drying front was verified in the field on six occasions and seven differential gauging campaigns were conducted to ground-truth the losses estimated from the satellite images. The transmission loss point data obtained using the wetted river lengths and differential gauging campaigns were used to train an ensemble of random forest models to predict the continuous hourly time series of transmission losses and their uncertainties. Our results show that the Selwyn River transmission losses ranged between 0.25 and 0.65 m3s-1km-1 during most of the 1-year study period. However, shortly after a flood peak the losses could reach up to 1.5 m3s-1km-1. These results enabled us to improve our understanding of the Selwyn River groundwater–surface water interactions and provide valuable data to support water management. We argue that our framework can easily be adapted to other ephemeral rivers and to longer time series.

Karst aquifers are important sources of fresh water on a global scale. The hydrological modelling of karst spring discharge, however, still poses a challenge. In this study we apply a transfer function noise (TFN) model in combination with a bucket-type recharge model to simulate and predict karst spring discharge. The application of the noise model for the residual series has the advantage that it is more consistent with assumptions for optimization such as homoscedasticity and independence. In an earlier hydrological modeling study, named Karst Modeling Challenge (KMC; Jeannin et al (2021)), several modelling approaches were compared for the Milandre Karst System in Switzerland. This serves as a benchmark and we apply the TFN model to KMC data, subsequently comparing the results to other models. Using different data-model combinations, the most promising data-model combination is identified in a three-step least-squares calibration. To quantify uncertainty, the Bayesian approach of Markov-chain Monte Carlo (MCMC) sampling is subsequently used with uniform priors for the previously identified best data-model combination. The MCMC maximum likelihood solution is used to predict spring discharge for the evaluation period, indicating a superior performance compared to all other models in the KMC. It is found that the model gives a physically feasible representation of the system, which is supported by field measurements. While the TFN model simulated rising limbs and flood recession especially well, medium and baseflow conditions were not represented as accurately. The TFN approach poses a well-performing data-driven grey-box alternative to other approaches that should be considered in future studies.

Highly detailed physically based groundwater models are often applied to make predictions of system states under unknown forcing. The required analysis of uncertainty is often unfeasible due to the high computational demand. We combine two possible solution strategies: (1) the use of faster surrogate models; and (2) a robust data worth analysis combining quick first-order second-moment uncertainty quantification with null-space Monte Carlo techniques to account for parametric uncertainty. A structurally and parametrically simplified model and a proper orthogonal decomposition (POD) surrogate are investigated. Data worth estimations by both surrogates are compared against estimates by a complex MODFLOW benchmark model of an aquifer in New Zealand. Data worth is defined as the change in post-calibration predictive uncertainty of groundwater head, river-groundwater exchange flux, and drain flux data, compared to the calibrated model. It incorporates existing observations, potential new measurements of system states ("additional" data) as well as knowledge of model parameters ("parametric" data). The data worth analysis is extended to account for non-uniqueness of model parameters by null-space Monte Carlo sampling. Data worth estimates of the surrogates and the benchmark suggest good agreement for both surrogates in estimating worth of existing data. The structural simplification surrogate only partially reproduces the worth of "additional" data and is unable to estimate "parametric" data, while the POD model is in agreement with the complex benchmark for both "additional" and "parametric" data. The variance of the POD data worth estimates suggests the need to account for parameter non-uniqueness, like presented here, for robust results.

A major effect of environment on crops is through crop phenology, and therefore, the capacity to predict phenology for new environments is important. Mechanistic crop models are a major tool for such predictions, but calibration of crop phenology models is difficult and there is no consensus on the best approach. Here we propose an original, detailed approach, a protocol, for calibration of such models. The protocol covers all the steps in the calibration work-flow, namely choice of default parameter values, choice of objective function, choice of parameters to estimate from the data, calculation of optimal parameter values and diagnostics. The major innovation is in the choice of which parameters to estimate from the data, which combines expert knowledge and data-based model selection. First, almost additive parameters are identified and estimated. This should make bias (average difference between observed and simulated values) nearly zero. These are “obligatory” parameters, that will definitely be estimated. Then candidate parameters are identified, which are parameters likely to explain the remaining discrepancies between simulated and observed values. A candidate is only added to the list of parameters to estimate if it leads to a reduction in BIC (Bayesian Information Criterion), which is a model selection criterion. A second original aspect of the protocol is the specification of documentation for each stage of the protocol. The protocol was applied by 19 modeling teams to three data sets for wheat phenology. All teams first calibrated their model using their “usual” calibration approach, so it was possible to compare usual and protocol calibration. Evaluation of prediction error was based on data from sites and years not represented in the training data. Compared to usual calibration, calibration following the new protocol reduced the variability between modeling teams by 22% and significantly reduced prediction error.

Bayesian mechanistic modeling often suffers from overconfident and biased posterior distributions for parameters and predictions. This phenomenon arises because the fundamental assumption of Bayesian Model Analysis is violated: the underlying model is assumed to be true, but in fact, it is a simplification of reality with structural errors that show at least during some periods of the modeled time span. As a result, a compromise solution in parameter space is identified that can formally fit the full data set best, but this parameter set will not be representative of the true system state. Neither will it be representative of the “compensation mode” in which the model is whenever structural error kicks in. As a logical consequence, predictions will be biased and their intervals too narrow. The longer the data set used for calibration, the stronger the misleading effect. Typical sources of severe structural deficits that produce dynamically occurring errors are missing or misspecified processes in the model. We propose a formal time-windowed Bayesian analysis to overcome this general problem. When performing Bayesian updating on shorter time windows, the assumption of a (quasi-) true model becomes more plausible, and by sliding this window through the calibration time series, we let the model adjust its posterior parameter distributions according to the current strength of error. These time-shifting parameter distributions allow us to (1) identify periods of statistically significant model error occurrence via measuring time-varying Bayesian model evidence, (2) diagnose potential sources of model error by understanding the time-varying parameter compensation mechanisms, and (3) predict with more realistic uncertainty intervals by distribution averaging. We demonstrate the proposed method on a set of synthetic and real-world scenarios of soil moisture modeling. With this example, we also highlight its usefulness to analyze dynamic systems in a wide range of disciplines, such as water quality modeling, decision support, and risk assessment. Results show that the time sequence of posterior parameter distributions (and dependent model mechanisms such as water retention curves and unsaturated hydraulic conductivity functions) provides valuable insights into the model’s weaknesses, and it also provides guidance for model improvement.

Efficient 1.5-D hillslope model using 1D Richards' and 1D Boussinesq equations coupled by the Method Of Lines Robert Mietrach, Thomas Wöhling Abstract Modeling lateral subsurface flow at the hill slope scale by considering the full 2D or 3D spacial domain can be quite complex from a computational point of view. To reduce simulation times a new prototype model consisting of a set of 1D Richards' equations linked with the Boussinesq equation is presented. The Method Of Lines approach to solving the Richards' equation has already been shown to be an efficient and stable alternative to established solution methods, such as low-order finite difference and finite element methods applied to the mixed form of Richards' equation. Besides its beneficial properties in numerical challenging scenarios, the Method Of Lines approach allows for easier integration of additional differential equations which proves advantageous here, when integrating the Boussinesq equation into the combined model. In the combined model it is assumed that lateral flow primarily happens at soil layer interfaces. The Boussinesq equation is than used to link the 1D columns at these soil layer interface nodes in the lateral direction. Thus enabling water transport between adjacent columns and therefor along the hillslope. In an analog procedure it would be possible to extend the presented model to also simulate for surface runoff. Simulations for several synthetic setups have been carried out and compared to solutions to the full 2D problem from the software Hydrus. The results show good agreement between the two approaches, with the benefit of reduced simulation times and increased numerical stability of the presented model.

The origins of the Richard’s equation date back more than a century but it is still the most commonly used model to describe variably saturated soil water movement. It requires the specification of the characteristic relationships between soil water content, tensiometric pressure and hydraulic conductivity, which can be subsumed as soil hydraulic models. Numerous soil hydraulic models and variants thereof have been developed to mimic the behaviour of natural soils, including hysteresis, macro-pore flow, and dynamic non-equilibrium flow. Despite the progress being made, it is often difficult to accurately simulate retention data derived from undisturbed field soils.This work presents a benchmark inverse modelling study for 1D soil water movement and field retention data from a wetting-drying cycle using state-of-the-art soil hydraulic models. The main aim is to test the ability of the different models to reproduce the field data. The soil hydraulic models tested are, among others, the van Genuchten-Mualem model (VGM, 1976, 1980), VGM with hysteresis (Kool and Parker, 1987), Brooks-Corey (1966), Dual porosity (Durner, 1994) and the non-equilibrium flow model by Diamantopoulus (2015).In our study, we used an implementation of the Richards equation with the highly efficient and numerically stable Methods-of-Line scheme. Best-fit estimates and parameter posterior distributions were derived using the Markov-Chain Monte Carlo sampling algorithm DREAM_ZS and time series of soil water content and tensiometric pressure. The field data shows clear signs of non-equilibrium flow. It originates from an intensively studied, inverted-lysimeter site with Pumice soils under grass from the central part of the North Island of New Zealand.Results demonstrate that none of the models was able to accurately mimic soil water content and tensiometric pressure data simultaneously at all times. Model deficiencies were identified particularly for the two wetting events, where all models underestimated soil water content while tensiometric pressure matched the data closely. We hypothesise that at least part of the discrepancies relate to an oversimplification of the hydraulic conductivity function for non-equilibrium flow.This study is limited to a single data set and by several assumptions that are commonly made in inverse parameter estimation studies. The better assessment and implementation of measurement error (structures) might alleviate (or mask) some of the discrepancies between model simulations and data. However, this is apparently not the solution to the problem. Dynamic non-equilibrium flow has been observed in natural soils in several well-conducted field experiments. Our results are just one example that demonstrates the need to improve soil hydraulic modelling by revisiting the physics of fundamental processes in natural soils.

In the coastal plains of New Zealand, braided rivers lose a considerable amount of water to the groundwater system, and for many aquifers are the largest source of recharge. Quantifying the recharge rates/transmission losses and how they relate to river stage and flow is particularly challenging. A commonly used approach to estimating recharge rates/transmission losses is differential flow gauging, where river discharges are measured simultaneously at multiple locations along a given reach. However, differential flow gauging is labour-intensive and limited by the accuracy of river discharge measurements, particularly in braided river systems.We have developed an alternative method for ephemeral rivers using river stage monitoring and widely available satellite photography. The method was applied to the upstream part of the Selwyn River (Canterbury, New Zealand), which is perennial in its mountainous environment, but becomes ephemeral once it crosses its alluvial plain. The river stage near the downstream boundary of the perennial reach was monitored for the period of March 2020 to May 2021 and a stage-discharge rating curve was developed. Downstream of the monitoring station, the river becomes ephemeral with the drying front location changing over time. On a number of 146 suitable satellite photographs, taken within the stage recording period, and retrieved from the Planet application program interface, we identified the position of the drying front, and used this to determine the length of the active (wet) river channel. This enabled us to calculate the average river transmission losses by dividing the river discharge at the monitoring station by the downstream active river length. The transmission losses estimated using the satellite photography correspond well with the losses estimated using seven sets of independent differential flow gauging surveys, given the respective uncertainties of both methods. The average estimated transmission losses range from 0.2 to 1 m3/s/km. Most of the estimated losses are below 0.4 m3/s/km and correspond to baseflow periods. The highest losses occur shortly after peak flows and decrease exponentially with time after the peak.We hypothesize that the high losses, occurring shortly after peak flows, are due to the replenishment of the shallow braid plain aquifer associated with the river. Lower losses, occurring during baseflow periods, represent groundwater recharge to the deeper regional aquifer. Groundwater recharge to the deeper regional aquifer appears to be linearly correlated with the groundwater head in the shallow aquifer. This response is consistent with the presence of an unsaturated zone that has been identified between the shallow (riverine) and deeper (regional) aquifers. Furthermore, we have successfully trained a random forest regression model to reconstruct the transmission losses for every day of the study period. The daily transmission loss dataset can now be used to evaluate our physically-based groundwater – surface water interaction models, currently under development, as well as support water management in the Selwyn basin.

Globally around half of the population resides in the tropics. Understanding the regional trends and physical controls of streamflow drought in the present day helps us to recognize the future changes in the hydrological cycle impacting the freshwater supplies. Few studies have examined the climatic and catchment controls on the propagation of streamflow droughts in tropical catchments. This study unveils the regional pattern of drought onset and deficit volume and identifies their evolution by analyzing the dominant climatic and physiographic controls. Two crucial streamflow drought signatures, the time of onset and deficit volume, are extracted from daily observed streamflow records of 82 rain-dominated catchments from peninsular India (8-24º N, 72-87º E). A daily variable threshold approach with an 80% exceedance probability of the flow record is used to identify drought events. Despite a decreasing trend in deficit volume, we show a delayed shift in drought onset of about a week for more than 50% of catchments. However, droughts with mean onset time clustered around the summer (Mar-May) season show a sharp rise in deficit volume trends. We show that while dynamic (precipitation and soil moisture) factors influence the onset of droughts, both static (catchment, topographic, and soil properties) and dynamic properties play a significant role in deciding the deficit volume. Based on a statistical approach (Taylor skill score and non-parametric dependence metrics), we identify static features, surface (30 cm) and sub-surface (up to 1m) soil organic stock and cation exchange capacity (CEC) to be dominant soil controls, also topographic ruggedness index, slope, topographic wetness index, and curvature are found to be the main controlling catchment attributes. The derived insights add new avenues in understanding the causal chain of physical processes linking climatic and physiographic controls on streamflow drought mechanisms, elucidating tropical climate response to water availability in a changing climate.

Estimating the effects of climate change on water resources is an important topic for researchers in hydrological sciences. Climate model outputs can be used as inputs for calibrated numerical groundwater models to predict the effect of different climate scenarios on groundwater levels. Unfortunately, spatially explicit numerical models are spatially constrained, data-hungry, and difficult to set-up and to calibrate. Furthermore, the involved model run-times make proper uncertainty analysis of the model predictions computationally expensive.Machine-learning (ML) tools have become recognized as a powerful alternative for numerical models for different applications in hydrological science, but come with their own challenges. They have been successfully used for the prediction of groundwater levels in single bores based on historical data, but their application for estimating groundwater levels for climate scenarios is still a matter of active research. To identify the potential of ML techniques for this application, two different ML algorithms have been applied to predict groundwater levels for several climate scenarios at a number of groundwater wells with long-term historical data time series. The two ML algorithms are (i) multi-layer perceptrons (MLP) with a closed feedback loop and (ii) long short-term memory (LSTM) networks. For each observation well, several thousand versions of these ML models with differing setups are trained to historical, monthly data time series up to 2015 and then applied for the estimation of monthly groundwater levels. These estimations are computed by using climate model outputs for different scenarios, as well as artificially generated scenarios, as drivers for the ML models to test their sensitivity and plausibility to those input data series. The high quantity of model versions for each bore are utilized to generate mean groundwater level estimates and accompanying uncertainty bands via Bayesian Model Averaging (BMA).Both ML techniques are able to match the historical data time series for the different bores with small uncertainty, but differ in their ability to predict long-term groundwater levels from climate change and artificial scenarios. The long-term simulations of the MLP models show believable trends and appropriate uncertainty bands, while the LSTM networks seem to underestimate the uncertainty of their future predictions.

Deterministic hydrological models with uncertain, but inferred‐to‐be‐time‐invariant parameters typically show time‐dependent model errors. Such errors can occur if a hydrological process is active in certain time periods in nature, but is not resolved by the model or by its input. Such missing processes could become visible during calibration as time‐dependent best‐fit values of model parameters. We propose a formal time‐windowed Bayesian analysis to diagnose this type of model error, formalizing the question “In which period of the calibration time‐series does the model statistically disqualify itself as quasi‐true?” Using Bayesian model evidence (BME) as model performance metric, we determine how much the data in time windows of the calibration time‐series support or refute the model. Then, we track BME over sliding time windows to obtain a dynamic, time‐windowed BME (tBME) and search for sudden decreases that indicate an onset of model error. tBME also allows us to perform a formal, sliding likelihood‐ratio test of the model against the data. Our proposed approach is designed to detect error occurrence on various temporal scales, which is especially useful in hydrological modeling. We illustrate this by applying our proposed method to soil moisture modeling. We test tBME as model error indicator on several synthetic and real‐world test cases that we designed to vary in error sources (structure and input) and error time scales. Results prove the successful detection errors in dynamic models. Moreover, the time sequence of posterior parameter distributions helps to investigate the reasons for model error and provide guidance for model improvement., Key Points: We propose a data‐driven method for model‐structural error detection. Our method rests on a statistically rigorous Bayesian framework without prior assumptions about error sources or patterns. We confirm successful error detection on various temporal scales in synthetic test cases and present insights from a real‐world case study., German Research Foundation (DFG), Cluster of Excellence, University of Stuttgart, https://doi.org/10.18419/darus-1836

Funding Information: This work was in part supported by the Collaborative Research Center 1253 CAMPOS (Project 7: Stochastic Modeling Framework), funded by the German Research Foundation ( DFG , Grant Agreement SFB 1253/1 2017 ), the Academy of Finland through projects AICropPro ( 316172 ) and DivCSA ( 316215 ) and Natural Resources Institute Finland (Luke) through a strategic project BoostIA, the BonaRes project ''Soil3'' ( BOMA 03037514 ) of the Federal Ministry of Education and Research ( BMBF ), Germany, the Deutsche Forschungsgemeinschaft ( DFG , German Research Foundation ) under Germany's Excellence Strategy - EXC 2070–390732324 EXC (PhenoRob), the project BiomassWeb of the GlobeE programme (Grant number: FKZ031A258B ) funded by the Federal Ministry of Education and Research ( BMBF , Germany), the INRA ACCAF meta-programme, the German Federal Ministry of Education and Research ( BMBF ) in the framework of the funding measure “Soil as a Sustainable Resource for the Bioeconomy – BonaRes”, project “BonaRes (Module B): BonaRes Centre for Soil Research , subproject B” (grant 031B0511B ), the National Key Research and Development Program of China ( 2017YFD0300205 ), the National Science Foundation for Distinguished Young Scholars ( 31725020 ), the Priority Academic Program Development of Jiangsu Higher Education Institutions ( PAPD ), the 111 Project ( B16026 ), and China Scholarship Council , the Agriculture and Agri-Food Canada's Project 1387 under he Canadian Agricultural Partnership, the DFG Research Unit FOR 1695 ‘Agricultural Landscapes under Global Climate Change – Processes and Feedbacks on a Regional Scale, the U.S. Department of Agriculture National Institute of Food and Agriculture (award no. 2015-68007-23133 ) and USDA / NIFA HATCH grant N. MCL02368 , the National Key Research and Development Program of China ( 2016YFD0300105 ), The Broadacre Agriculture Initiative, a research partnership between University of Southern Queensland and the Queensland Department of Agriculture and Fisheries , the Academy of Finland through project AI-CropPro ( 315896 ), the JPI FACCE MACSUR2 project, funded by the Italian Ministry for Agricultural, Food, and Forestry Policies ( D.M. 24064/7303/15 of 6/Nov/2015), the Ministry of Education, Youth and Sports of Czech Republic through SustES - Adaption strategies for sustainable ecosystem services and food security under adverse environmental conditions (project no. CZ.02.1.01/0.0/0.0/16_019/000797 ). The order in which the donors are listed is arbitrary. Funding Information: This work was in part supported by the Collaborative Research Center 1253 CAMPOS (Project 7: Stochastic Modeling Framework), funded by the German Research Foundation (DFG, Grant Agreement SFB 1253/1 2017), the Academy of Finland through projects AICropPro (316172) and DivCSA (316215) and Natural Resources Institute Finland (Luke) through a strategic project BoostIA, the BonaRes project ''Soil3'' (BOMA 03037514) of the Federal Ministry of Education and Research (BMBF), Germany, the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC 2070?390732324 EXC (PhenoRob), the project BiomassWeb of the GlobeE programme (Grant number: FKZ031A258B) funded by the Federal Ministry of Education and Research (BMBF, Germany), the INRA ACCAF meta-programme, the German Federal Ministry of Education and Research (BMBF) in the framework of the funding measure ?Soil as a Sustainable Resource for the Bioeconomy ? BonaRes?, project ?BonaRes (Module B): BonaRes Centre for Soil Research, subproject B? (grant 031B0511B), the National Key Research and Development Program of China (2017YFD0300205), the National Science Foundation for Distinguished Young Scholars (31725020), the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD), the 111 Project (B16026), and China Scholarship Council, the Agriculture and Agri-Food Canada's Project 1387 under he Canadian Agricultural Partnership, the DFG Research Unit FOR 1695 ?Agricultural Landscapes under Global Climate Change ? Processes and Feedbacks on a Regional Scale, the U.S. Department of Agriculture National Institute of Food and Agriculture (award no. 2015-68007-23133) and USDA/NIFA HATCH grant N. MCL02368, the National Key Research and Development Program of China (2016YFD0300105), The Broadacre Agriculture Initiative, a research partnership between University of Southern Queensland and the Queensland Department of Agriculture and Fisheries, the Academy of Finland through project AI-CropPro (315896), the JPI FACCE MACSUR2 project, funded by the Italian Ministry for Agricultural, Food, and Forestry Policies (D.M. 24064/7303/15 of 6/Nov/2015), the Ministry of Education, Youth and Sports of Czech Republic through SustES - Adaption strategies for sustainable ecosystem services and food security under adverse environmental conditions (project no. CZ.02.1.01/0.0/0.0/16_019/000797). The order in which the donors are listed is arbitrary. Publisher Copyright: © 2021 Elsevier Ltd Calibration, the estimation of model parameters based on fitting the model to experimental data, is among the first steps in many applications of process-based models and has an important impact on simulated values. We propose a novel method of developing guidelines for calibration of process-based models, based on development of recommendations for calibration of the phenology component of crop models. The approach was based on a multi-model study, where all teams were provided with the same data and asked to return simulations for the same conditions. All teams were asked to document in detail their calibration approach, including choices with respect to criteria for best parameters, choice of parameters to estimate and software. Based on an analysis of the advantages and disadvantages of the various choices, we propose calibration recommendations that cover a comprehensive list of decisions and that are based on actual practices.

International audience; The complexity of karst groundwater flow modelling is reflected by the amount of simulation approaches. The goal of the Karst Modelling Challenge (KMC) is comparing different approaches on one single system using the same data set. Thirteen teams with different computational models for simulating discharge variations at karst springs have applied their respective models on one single data set coming from the Milandre Karst Hydrogeological System (MKHS). The approaches include neural networks, reservoir models, semi-distributed models and fully distributed groundwater models. Four and a half years of hourly or daily meteorological input and hourly discharge data were provided for model calibration. The validation comprised forecasting one year of discharge, without the observed discharge data. The model performance was evaluated using the volume conservation, Nash-Sutcliffe efficiency (NSE) and the Kling-Gupta efficiency (KGE) applied on the total discharge and individual flow components. As a result, the comparison of model performances is a challenging task due to the differences in the model architecture but also required time steps: some of the models require aggregated daily steps while others could be run using hourly data, which provided some interesting differences depending on how the data was transformed. The use of instantaneous data (e.g. value at noon) produces less bias that averaging hourly data over one day. The transformation of hourly into daily data produces a decrease of Nash and KGE of 0.05 to 0.08 (i.e. from 1 to ~0.93). The resulting simulations (forecasted values for year 2016) produced KGEs ranging between 0.83 and 0.37 (0.83 to −0.24 for NSE). Although the simulations matched the monitored flows reasonably well, most models struggled to simulate baseflow conditions accurately. In general, the models that performed the best for this exercise were the global ones (Gardenia and Varkarst), with a limited number of parameters, which can be calibrated using automatic calibration procedures. The neural network models also showed a fair potential, with one providing reasonable results despite the relatively short dataset available for warming-up (4.5 years). Semi-and fully distributed models also suggested that with some more effort they could perform well. The accuracy of model predictions does not seem to increase by using models with more than 9–12 calibration parameters. An evaluation of the relative errors between the forecasted and the observed values revealed that for most models, 50% of the forecasted values contained more than 50% of difference against the observed discharge rate, with 25% having a difference larger than 100%. A significant part of the poorly forecasted values corresponded to base-flow which was surprising given that as base-flow is generally much easier to predict than peak flow. Hence, this shows that modelling approaches and criteria for the calibration are too oriented towards peak-flow sections of the hydrographs, and that improvements could be gained by more focus on the base-flow.

Hyporheic exchange is the interaction of river water and groundwater, and is difficult to predict. One of the largest contributions to predictive uncertainty for hyporheic exchange has been attributed to the representation of heterogeneous subsurface properties. Our study evaluates the trade-offs between intrinsic (irreducible) and epistemic (reducible) model errors when choosing between homogeneous and highly complex subsurface parameter structures. We modeled the Steinlach River Test Site in Southwest Germany using a fully coupled surface water-groundwater model to simulate hyporheic exchange and to assess the predictive errors and uncertainties of transit time distributions. A highly parameterized model was built, treated as a "virtual reality" and used as a reference. We found that if the parameter structure is too simple, it will be limited by intrinsic model errors. By increasing subsurface complexity through the addition of zones or heterogeneity, we can begin to exchange intrinsic for epistemic errors. Thus, the appropriate level of detail to represent the subsurface depends on the acceptable range of intrinsic structural errors for the given modeling objectives and the available site data. We found that a zonated model is capable of reproducing the transit time distributions of a more detailed model, but only if the geological structures are known. An interpolated heterogeneous parameter field (cf. pilot points) showed the best trade-offs between the two errors, indicating fitness for practical applications. Parameter fields generated by multiple-point geostatistics (MPS) produce transit time distributions with the largest uncertainties, however, these are reducible by additional hydrogeological data, particularly flux measurements.

Numerical models for reactive transport can be used to estimate the breakthrough of a contaminant in a pumping well or at other receptors. However, as natural aquifers are highly heterogeneous with unknown spatial details, reactive transport predictions on the aquifer scale require a stochastic framework for uncertainty analysis. The high computational demand of spatially explicit reactive-transport models hampers such analysis, thus motivating the search for simplified estimation tools. We suggest performing an electron balance between the reactants in the infiltrating solution and in the aquifer matrix to obtain the hypothetical time of dissolved-reactant breakthrough at a receptor if the reaction with the matrix was instantaneous. This time we denote as the advective breakthrough time for instantaneous reaction (τinst ). It depends on the amount of the reaction partner present in the matrix, the mass flux of the dissolved reactant, and the stoichiometry. While the shape of the reactive-species breakthrough curve depends on various kinetic parameters, the overall timing scales with τinst . We calculate the latter by particle tracking. The effort of computing τinst is so low that stochastic calculations become feasible. We apply the concept to a two-dimensional test case of aerobic respiration and denitrification. A detailed spatially explicit reactive-transport model includes microbial dynamics. Scaling the time of local breakthrough curves observed at individual points by τinst decreased the variability of electron-donor breakthrough curves significantly. We conclude that the advective breakthrough time for instantaneous reaction is efficient in estimating the time over which an aquifer retains its degradation potential.

Groundwater resource management often requires detailed numerical models that make calibration and predictive uncertainty analysis computationally challenging. Reduced order models (ROMs) alleviate the computational burden but can potentially lead to predictive bias and underestimation of uncertainty. A paired model approach has previously been proposed to estimate the predictive uncertainty of models compared to highly complex, synthetic realities. This approach is modified to analyze and compare the simplification error for groundwater predictions of a real-world numerical MODFLOW model of the Wairau Plains Aquifer in the Marlborough Region of New Zealand. Two different ROM types were applied in this study to predict groundwater heads, spring discharge and river–groundwater exchange fluxes: (1) a drastically simplified MODFLOW model, and (2) artificial neural networks (ANNs). The different ROMs exhibit very different patterns of bias and magnitude of model simplification error. The method accurately captures the simplification error for most predictions by the MODFLOW model, but underestimates the error for predictions highly dependent on the variability of the complex model. The simplified MODFLOW model shows significant parameter surrogacy and non-optimality of simplification, thus questioning the adherence to expert-knowledge based parameter limits. For predictions where historic data sets are available, ANNs provide superior predictive power. However, they cannot be applied to predictions of data types and locations not contained in the calibration data set. For those predictions, simplified numerical models can be applied with varying degree of accuracy.

The classical formulation of Richards' equation is relying on a unique functional relationship between water content, conductivity and pressure head. Some phenomena like hystersis effects in the water content during wetting and drying cycles and hydraulic non-equillibrium cannot be accounted for with this formulation. Therefor it has been extended in different ways in the past to be able to include these effects in the simulation. Each modification comes with its own challenges regarding implementation and numerical stability.The Method Of Lines approach to solving the Richards' equation has already be shown to be an efficient and stable alternative to established solution methods, such as low-order finite difference and finite element methods applied to the mixed form of Richards' equation.In this work a slightly modified Method Of Lines approach is used to solve the pressure based 1D Richards' equation. A finite differencing scheme is applied to the spatial derivative and the resulting system of ordinary differential equations is reformulated as differential-algebraic system of equations. The open-source code IDAS from the Sundials suite is used to solve the DAE system. Different extensions to Richards' equation have been incorporated into the model to address the shortcomings mentioned above. These extensions are a model able to simulate preferential flow using a coupled two domain approach, a simple hysteretic model to account for hysteresis in the water retention curve and also two models to either fully or partially calculate hydraulic non-equillibrium effects. To verify the numerical robustness of the extended model, stochastic parameterizations were generated that represent the full range of all soil types. Simulations were carried out using these parameter sets and real-world meteorological boundary conditions at 10 minutes time intervals, that exhibit drastic flux changes and poses numerical challenges for classical solution methods.The results show that not only does the extended model converge for all parameterizations, but that numerical robustness and performance is maintained. Where applicable the results have been verified against solutions from the software Hydrus and show good agreement with those.

For decades, lumped rainfall-runoff models have been used for hydrological analysis and forecasting such as operational flood forecasting. However, the accuracy of model forecasts depends on the ability of the model to approximate the dominant hydrological processes of the catchment under consideration. These processes are site specific and, therefore, the choice of a particular model is challenging. A large number of hydrological models has been developed and applied in various regions of the world. Model choice has often been hampered in the past by technical problems such as different programming languages, different software platforms, and different input formats and requirements.The Modular Assessment of Rainfall Runoff Model Toolbox (MARRMoT) unites 46 lumped models from around the world within the same Matlab® framework with standardized inputs. The model equations have been simplified and adapted for this purpose. As a result, it is possible to test a large number of different models with comparatively little effort. The models implemented in MARRMoT vary in their structural complexity and have between 1-24 parameters and between 1-8 storages.Here MARRMoT was used in order to find a model or model ensemble suitable for the simulation of precipitation-runoff relationships in the Wairau River catchment, New Zealand. The catchment area is assumed to have predominantly homogeneous runoff-generating properties. Model input data (precipitation and potential evapotranspiration) was derived from the Virtual Climate Station Network by National Institute of Water and Atmospheric Research, NZ.In a first scenario, 42 selected models from MARRMoT were calibrated for the Wairau River catchment using 45 years of Wairau River flow data, an in-built nonlinear unconstrained optimization algorithm and the model fitness criteria Kling-Gupta-Efficiency (KGE). In two further scenarios, calibrations using the KGE with inversely transformed flows (KGEi) as well as a mixed form of the two criteria (KGEm) were realized. Model performance was further evaluated based on different performance criteria such as NSE, RMSE and R². It was demonstrated that the model ranking depends on the choice of the performance.Evaluating the model performance for the different calibration scenarios showed that a few models with very different structures performed well to reproduce the flow data. No decisive structural feature could be identified which all models have in common and which led to a good representation of the rainfall-runoff processes in the Wairau River catchment. However, the differentiated consideration of flow routing and a high degree of flexibility seem to benefit model performance. Deficits in the modeling can be seen in the discharge peaks, which are not correctly simulated by many models. The simulation of fast direct runoff with lumped models seems to be less accurate for the relatively large catchment area of ​​the Wairau River (3430 km²).Eventually, three models (GR4J, FLEX-I and HBV-96) demonstrated a high performance in all three calibration scenarios and were identified as suitable for further use in the Wairau River catchment.

Groundwater is a major resource for drinking water supply and irrigation of crops in many parts of the world. Many groundwater aquifers are already fully allocated, but the demand is projected to increase further while, concurrently, climate change may cause more variability in the natural supply. This poses enormous challenges for the future management of groundwater resources and a paradigm shift from traditional, threshold-based management strategies to more flexible, adaptive management strategies.For that purpose, operational forecasting tools are required that predict future states of groundwater aquifers under various scenarios and to predict the risk of critical states which would have adverse effects either for the environment or for water users, or both.The mathematical description of the complex interactions particularly of shallow, unconfined river-fed aquifers typically requires the use of spatially explicit numerical models. These are, however, not suitable for operational forecasting due to lengthy run-times and extensive data requirements. This also poses strong limitations with respect to predictive uncertainty analysis – which should be an integral part of predictive management tools. Model simplification or model surrogates are the method of choice to circumvent the problem.An operational forecasting tool is presented here to predict groundwater heads and groundwater storage in the unconfined Wairau Plain Aquifer in Marlborough, New Zealand, during flow recession times. The tool uses low-complexity “eigenmodels” to describe groundwater flow and to provide an early warning for critical groundwater storage levels to the Marlborough District Council, which manages the groundwater resource. These critical levels have been approached more frequently during the past years when the natural recession of groundwater storage in summer is exacerbated by groundwater abstraction to satisfy the irrigation water demand of the Plain’s viticulture.The forecasting tool requires, amongst others, daily forecasts of Wairau River flows because the river is the major recharge source for the aquifer. Flow forecasts and their uncertainty are computed i) by using a master recession curve for predictions during flow recession times and ii) by a lumped rainfall-runoff model for times of aquifer storage recovery. This allows a broad evaluation of forecasting scenarios. The tool has been tested and is operational for recession times (worst-case scenario predictions; Wöhling & Burbery, 2020). The rainfall-runoff model performs reasonably well in predicting river flows and correspondingly in predicting groundwater storage recovery for historic data (hindcasting). The 30-day predictive uncertainty bounds generally cover the observations of river flows, groundwater levels and aquifer storage. The predictive accuracy of the tool largely depends on the predictive accuracy of the drivers, particularly of the areal estimates of precipitation that drives the rainfall-runoff model and the river-groundwater exchange function that describes aquifer recharge rates. References Wöhling T, Burbery L (2020). Eigenmodels to forecast groundwater levels in unconfined river-fed aquifers during flow recession. Science of the Total Environment, 747, 141220, doi: 10.1016/j.scitotenv.2020.141220.

Rising concentrations of dissolved organic carbon (DOC) in inland waters are observed and investigated intensely in the last decades. The development of adaptive measures requires the forecasting of DOC-exports from catchments. Since DOC is exported from river catchments along hydrological pathways it is evident that the investigation of runoff generation, retention and travel times along flow paths are important to quantify DOC-loads and to develop a forecast model.To gain comprehensive insights in runoff formation and DOC export in a small forested catchment in the Bavarian Forest National Park we apply a nested multi-tracer approach, combining experimental and analytical methods with the aim to develop a hydrological forecast model which is able to reproduce the dominant mobilization- and export processes of DOC in forested mountain catchments. The use of multiple tracers combines different approaches to determine source areas, flow paths and retention times of runoff water in catchments. Stable isotopes (d2H, d18O) are suitable as natural tracers to estimate contributions from precipitation to stream discharge. With the additional use of geochemical tracers (e.g. DOC, SiO2) contributions from groundwater and the organic and mineral soil horizons can be estimated. Combined with a nested approach these analyses can be conducted on different spatial scales, enabling the development of scalable prognostic models of runoff formation in catchments.To complement the limited information from historic data sets we instrumented two hill transects to observe lateral contributions from hill slopes and to investigate potential preferential flow paths. Water samples from stream-, soil-, ground- and precipitation water were collected during two flood events and analysed for stable isotopes and chemical compositions. To support the nested approach, the sampling sites were chosen at strategical sites within the catchment, including the instrumented hill transects and the stream network from the creek to the catchment outlet.Preliminary results of stable isotope analysis show, that after dry periods nearly no event water seems to contribute to runoff formation, whereas after wet periods the proportions can be up to 40 %. A strongly delayed reaction of the groundwater was observed which suggests that deep groundwater is not contributing to stream flow, but a possible mobilization of pre-event water in the riparian zone was observed as a response to precipitation events.A likely major source of DOC is in the organic soil horizons due to storage and degradation of organic material. This is supported by higher DOC-concentrations in the soil water from these horizons. In how far residence times, precipitation intensities and flow paths activation from different source areas influence concentration peaks of DOC in the stream will be analysed in the next steps.The results of the recent field campaign help to identify the dominant processes of runoff generation and DOC mobilization on different temporal and spatial scales and for different antecedent system states. The data and insights gained from the field campaign will be used to develop and calibrate process models for hypothesis testing and further analyses to eventually develop a forecast model for DOC mobilization.

Plant phenology, which describes the timing of plant development, is a major aspect of plant response to environment and for crops, a major determinant of yield. Since climate change is projected to alter crop phenology worldwide, there is a large effort to predict phenology as a function of environment. Many studies have focused on comparing model equations for describing how phenology responds to weather but the effect of crop model calibration, also expected to be important, has received much less attention. The objective here was to obtain a rigorous evaluation of prediction capability of wheat crop phenology models, and to analyze the role of calibration. The 27 participants in this multi-model study were provided experimental data for calibration and asked to submit predictions for sites and years not represented in those data. Participants were instructed to use and document their 99usual99 calibration approach. Overall, the models provided quite good predictions of phenology (median of mean absolute error of 6.1 days) and did much better than simply using the average of observed values as predictor. Calibration was found to compensate to some extent for differences between models, specifically for differences in simulated time to emergence and differences in the choice of input variables. Conversely, different calibration approaches led to major differences in prediction error between models with the same structure. Given the large diversity of calibration approaches and the importance of calibration, there is a clear need for guidelines and tools to aid with calibration. Arguably the most important and difficult choice for calibration is the choice of parameters to estimate. Several recommendations for calibration practices are proposed. Model applications, including model studies of climate change impact, should focus more on the data used for calibration and on the calibration methods employed.

Predicting wheat phenology is important for cultivar selection, for effective crop management and provides a baseline for evaluating the effects of global change. Evaluating how well crop phenology can be predicted is therefore of major interest. Twenty-eight wheat modeling groups participated in this evaluation. Model predictions depend not only on model structure but also on the parameter values. This study is thus an evaluation of modeling groups, which choose the structure and fix or estimate the parameters, rather than an evaluation just of model structures. Our target population was wheat fields in the major wheat growing regions of Australia under current climatic conditions and with current local management practices. The environments used for calibration and for evaluation were both sampled from this same target population. The calibration and evaluation environments had neither sites nor years in common, so this is a rigorous evaluation of the ability of modeling groups to predict phenology for new sites and weather conditions. Mean absolute error (MAE) for the evaluation environments, averaged over predictions of three phenological stages and over modeling groups, was 9 days, with a range from 6 to 20 days. Predictions using the multi-modeling group mean and median had prediction errors nearly as small as the best modeling group. For a given modeling group, MAE for the evaluation environments was significantly correlated with MAE for the calibration environments, which suggests that it would be of interest to test ensemble predictors that weight individual modeling groups based on performance for the calibration data. About two thirds of the modeling groups performed better than a simple but relevant benchmark, which predicts phenology by assuming a constant temperature sum for each development stage. The added complexity of crop models beyond just the effect of temperature was thus justified in most cases. Finally, there was substantial variability between modeling groups using the same model structure, which implies that model improvement could be achieved not only by improving model structure, but also by improving parameter values, and in particular by improving calibration techniques.

A robust solution to Richards' equation for complex soil hydraulic models using the Method Of LinesRobust numerical solutions are required for automatic parameter estimation, uncertainty analysis of soil hydraulic models, but also to quantify and legitimate model complexity.The Method Of Lines approach to solve Richards' equation has already be shown to be an efficient and stable alternative to established methods, namely low-order finite difference and finite element methods applied to the mixed form of Richards' equation. Besides its beneficial properties in numerical challenging scenarios, the Method Of Lines approach allows for easier integration of additional differential equations which proves advantageous where further processes should be included in the modeling.In this work a slightly modified Method Of Lines approach is used to solve the pressure based Richards' equation. A finite differencing scheme is applied to the spatial derivative and the resulting system of ordinary differential equations is reformulated as differential-algebraic system of equations (DAE). The open-source code IDAS from the Sundials suite is used to solve the DAE system. This solution has been extended to include hydraulic models that can account for hysteresis, dual-permeability and non-equilibrium effects. The different hydraulic model implementations have been verified against results from the software Hydrus and show good agreement with those.Bayesian model selection techniques and the concept of the model confusion matrix can be used to examine the legitimacy of a given model's complexity with regards to available input data.To generate the necessary data a Monte Carlo Sampling over a range of soil parameters was carried out for the models of different complexity. The computations were performed at the high-performance computing facilities at TU Dresden using the developed code.The results of the analysis show the identifiability of the models, i.e. how well a model recognizes itself through Bayesian model selection when it was the one that has generated the data. This is a useful technique when building model ensembles for diagnostic or predictive purposes.

The identification of origin, flow paths and transit times of water in catchments is an important component for process-based model development for runoff prediction. Hydrological studies offer, combined with isotope data, the possibility to quantify interactions between different compartments in catchments. In the context of this work it is examined to what extent event sampling of precipitation, streamflow, soil water and groundwater and the evaluation of their isotopic ratios δ2H and δ18O enable complex hydrological process investigations in the small forested mountain catchment of the river Große Ohe in the Bavarian Forest National Park. Within this study process analyses are carried out on small scales, e.g. runoff formation on hill slopes and on catchment scale as integrative process analysis. The water samples were collected during a small flood event and analysed for the isotope ratios δ2H and δ18O using a Picarro. A hydrograph separation was carried out through a comprehensive evaluation of the concentration profiles during the event. In combination with further hydrological and soil hydrological observations possible areas of origin and retention times of the water were determined. A strongly delayed reaction of the groundwater was observed which suggests that groundwater is not contributing to stream flow during a flood event, but a possible mobilization of pre-event water in the riparian zone can be observed as a response to precipitation events. The knowledge gained hereby is the basis for further process analysis and model development.

The Upper Wairau Plain Aquifer serves as the major resource for drinking water and irrigation in the region of Blenheim in Marlborough, New Zealand. Natural recession of groundwater levels and storage occurs annually during the summer months in the upper part of the highly conductive gravel aquifer. Due to a number of particularly dry summers, aquifer storage has reached critical levels several times in the past which could lead to future restrictions on groundwater abstraction to be imposed by the Marlborough District Council (MDC) who manages the resource. The MDC requires at any given time an early warning whether or not these critical levels are likely to be reached. Correspondingly, an operational framework was developed to forecast Wairau Plains Aquifer groundwater levels and storage. The tool creates more lead time for the MDC in decision making and adaptive management of the Wairau Plain groundwater resources.A numerical groundwater flow model of the Wairau Plain was previously set up to understand the main drivers of aquifer storage (Wöhling et al. 2018). Since that model posed practical restrictions to be of use for operational management purposes, we tested several data-driven surrogate models with easily attainable inputs that could be derived directly by an automated database query of the MDC monitoring network. Here, a tailor-made version of the Eigenmodel approach (Sahuquillo, 1983) is used to predict Wairau Aquifer groundwater heads and coupled with Markov chain Monte Carlo (MCMC) sampling for model calibration and parameter uncertainty analysis. Several Eigenmodels are embedded in a modular prediction framework that allows for a flexible description of critical model inputs depending on different states of knowledge and on different purposes of the analysis. A Wairau River flow master-recession curve has been derived from historic (observed) time series data to provide the boundary condition for the major recharge source of the aquifer.The Eigenmodels perform very well in hind-casting the recessions of historic groundwater levels at selected locations of the Wairau Plain Aquifer. Periods with critical groundwater levels were successfully detected and accurately reproduced. The models are efficient and fast, which is a prerequisite for the operational management support tool. The way how aquifer recharge is described as a function of river discharge proved to be very sensitive to the accuracy of the results. Future plans are investigations on an improved knowledge of this relationship, but also to implement the propagation of input uncertainty through the framework, in addition to the treatment of parametric and predictive uncertainty which is implemented already. ReferencesSahuquillo, A. (1983). An eigenvalue numerical technique for solving unsteady linear groundwater models continuously in time. Water Resources Research 19(1): 87-93.Wöhling T., Gosses, M. Wilson, S., Davidson, P. (2018). Quantifying river-groundwater interactions of New Zealand's gravel-bed rivers: The Wairau Plain. Groundwater, 56(4), 647-666.

Physically-based groundwater models allow highly detailed spatial resolution, parameterization and process representation, among other advantages. Unfortunately, their size and complexity make many model applications computationally demanding. This is especially problematic for uncertainty and data worth analysis methods, which often require many model runs.To alleviate the problem of high computational demand for the application of groundwater models for data worth analysis, we combine two different solutions:a) the use of surrogate models as faster alternatives to a complex model, and b) a robust data worth analysis method that is based on linear predictive uncertainty estimation, coupled with highly efficient null-space Monte Carlo techniques. We compare the performance of a complex benchmark model of a real-world aquifer in New Zealand to two different surrogate models: a spatially and parametrically simplified version of the complex model, and a projection-based surrogate model created with proper orthogonal decomposition (POD). We generate predictive uncertainty estimates with all three models using linearization techniques implemented in the PEST Toolbox (Doherty 2016) and calculate the worth of existing, “future” and “parametric” data in relation to predictive uncertainty. To somewhat account for non-uniqueness of the model parameters, we use null-space Monte Carlo methods (Doherty 2016) to efficiently generate a multitude of calibrated model parameter sets. These are used to compute the variability of the data worth estimates generated by the three models.Comparison between the results of the complex benchmark model and the two surrogates show good agreement for both surrogates in estimating the worth of the existing data sets for various model predictions. The simplified surrogate model shows difficulties in estimating worth of “future” data and is unable to reproduce “parametric” data worth due to its simplification in parameter representation. The POD model was able to successfully reproduce both “future” and “parametric” data worth for different predictions. Many of its data worth estimates exhibit a high variance, though, demonstrating the need of robust data worth methods as presented here which (to some degree) can account for parameter non-uniqueness. Literature:Doherty, J., 2016. PEST: Model-Independent Parameter Estimation - User Manual. Watermark Numerical Computing, 6th Edition.

Excessive groundwater pumping due to immense agricultural, industrial and municipal demand poses a major threat of aquifer depletion in many areas around the world. The impact of climate change on the global hydrological cycle has further exacerbated the situation. Accurate and reliable prediction of long-term aquifer balance terms is a key prerequisite to manage groundwater sustainably. To deal with uncertainties of such predictions, lumped (conceptual) hydrological models could help with their computational speed that allows for Monte-Carlo simulation. Compared to more complex models, lumped models are fast, lean on data requirement and capable to quantify uncertainty. However, lumped models are mainly designed to simulate river discharge only, not aquifer storage. Even the standard practice for calibrating lumped hydrological models only includes river discharge, as data on groundwater storage is not directly accessible. In this study, we hypothesize that we can extend the HBV model by additional water budget and groundwater storage terms, and calibrated it on both groundwater storage data and discharge data. Then, we test whether its predictions of groundwater storage levels withstand validation tests. To avoid problems with unavailability of data for calibration and validation in a first proof of concept, we build a virtual reality with a MODFLOW-based model, driven with synthetic weather data over a period of more than 50 years. For rigorous testing, we cast calibration into the framework of Bayesian parameter inference, and validate with metrics that assess the appropriateness of the Bayesian prediction distribution of groundwater storage. We test our idea in the Wairau Plain aquifer, New Zealand. Poor understanding of recharge mechanisms and hence declining groundwater levels are the major hindrance for sustainable groundwater management in our study area. We pay specific attention to river-groundwater exchange processes, to the forecast of aquifer storage dynamics, and to groundwater depletion in a hypothetical, persistent draught. The purpose is to provide a proof of concept whether lumped models can be adapted and made suitable to predict declining groundwater resources up to full depletion, as an uncertainty-aware decision support system for sustainable management.

In recent years, proper orthogonal decomposition (POD) has become a popular model reduction method in the field of groundwater modeling. It is used to mitigate the problem of long run times that are often associated with physically-based modeling of natural systems, especially for parameter estimation and uncertainty analysis. POD-based techniques reproduce groundwater head fields sufficiently accurate for a variety of applications. However, no study has investigated how POD techniques affect the accuracy of different boundary conditions found in groundwater models. We show that the current treatment of boundary conditions in POD causes inaccuracies for these boundaries in the reduced models. We provide an improved method that splits the POD projection space into a subspace orthogonal to the boundary conditions and a separate subspace that enforces the boundary conditions. To test the method for Dirichlet, Neumann and Cauchy boundary conditions, four simple transient 1D-groundwater models, as well as a more complex 3D model, are set up and reduced both by standard POD and POD with the new extension. We show that, in contrast to standard POD, the new method satisfies both Dirichlet and Neumann boundary conditions. It can also be applied to Cauchy boundaries, where the flux error of standard POD is reduced by its head-independent contribution. The extension essentially shifts the focus of the projection towards the boundary conditions. Therefore, we see a slight trade-off between errors at model boundaries and overall accuracy of the reduced model. The proposed POD extension is recommended where exact treatment of boundary conditions is required.

Selecting a “best” model among several competing candidate models poses an often encountered problem in water resources modeling (and other disciplines which employ models). For a modeler, the best model fulfills a certain purpose best (e.g., flood prediction), which is typically assessed by comparing model simulations to data (e.g., stream flow). Model selection methods find the “best” trade-off between good fit with data and model complexity. In this context, the interpretations of model complexity implied by different model selection methods are crucial, because they represent different underlying goals of modeling. Over the last decades, numerous model selection criteria have been proposed, but modelers who primarily want to apply a model selection criterion often face a lack of guidance for choosing the right criterion that matches their goal. We propose a classification scheme for model selection criteria that helps to find the right criterion for a specific goal, i.e., which employs the correct complexity interpretation. We identify four model selection classes which seek to achieve high predictive density, low predictive error, high model probability, or shortest compression of data. These goals can be achieved by following either nonconsistent or consistent model selection and by either incorporating a Bayesian parameter prior or not. We allocate commonly used criteria to these four classes, analyze how they represent model complexity and what this means for the model selection task. Finally, we provide guidance on choosing the right type of criteria for specific model selection tasks. (A quick guide through all key points is given at the end of the introduction.)

The common practice of infrequent (e.g., monthly) stream water quality sampling for state of the environment monitoring may, when combined with high resolution stream flow data, provide sufficient information to accurately characterise the dominant nutrient transfer pathways and predict annual catchment yields. In the proposed approach, we use the spatially lumped catchment model StreamGEM to predict daily stream flow and nitrate concentration (mg L −1 NO 3 -N) in four contrasting mesoscale headwater catchments based on four years of daily rainfall, potential evapotranspiration, and stream flow measurements, and monthly or daily nitrate concentrations. Posterior model parameter distributions were estimated using the Markov Chain Monte Carlo sampling code DREAM ZS and a log-likelihood function assuming heteroscedastic, t-distributed residuals. Despite high uncertainty in some model parameters, the flow and nitrate calibration data was well reproduced across all catchments (Nash-Sutcliffe efficiency against Log transformed data, NSL, in the range 0.62–0.83 for daily flow and 0.17–0.88 for nitrate concentration). The slight increase in the size of the residuals for a separate validation period was considered acceptable (NSL in the range 0.60–0.89 for daily flow and 0.10–0.74 for nitrate concentration, excluding one data set with limited validation data). Proportions of flow and nitrate discharge attributed to near-surface, fast seasonal groundwater and slow deeper groundwater were consistent with expectations based on catchment geology. The results for the Weida Stream in Thuringia, Germany, using monthly as opposed to daily nitrate data were, for all intents and purposes, identical, suggesting that four years of monthly nitrate sampling provides sufficient information for calibration of the StreamGEM model and prediction of catchment dynamics. This study highlights the remarkable effectiveness of process based, spatially lumped modelling with commonly available monthly stream sample data, to elucidate high resolution catchment function, when appropriate calibration methods are used that correctly handle the inherent uncertainties.

We simulate aquifer-scale reactive transport using an approach based on travel times and relative reactivity. The latter quantifies the intensity of the chemical reaction relative to a reference reaction rate with identical concentrations and can be interpreted as the strength of electron-donor (or electron-acceptor) release by the matrix, scaled by a reference release. In general, the relative reactivity is a spatially variable property reflecting the geology of the formation. In the proposed approach, we track the path of individual water parcels through the aquifer and evaluate the age of the water parcels and the relative reactivity integrated along their trajectories. By switching from spatial discretization to cumulative relative reactivity, advective-reactive transport can be simulated by solving a single system of ordinary differential equations for each combination of concentrations in the inflow. We test the validity of the approach in a two-dimensional test case of steady-state groundwater flow and reactive transport involving aerobic respiration and denitrification. Here we compare steady-state concentration distributions of the spatially explicit virtual truth, accounting for dispersive mixing, with the approximation based on cumulative relative reactivity and show that the errors introduced by neglecting dispersive mixing are minor if the target quantities are the mass fluxes crossing a control plane or being collected by a well. We further demonstrate the efficiency of the approach in a synthetic three-dimensional case study. The proposed approach is computationally so efficient, that ensemble runs to assess statistical distributions of concentration time series of reactive solutes become feasible, which is not practical with a spatially explicit model. This article is protected by copyright. All rights reserved.

Monitoring networks are expensive to establish and to maintain. In this paper, we extend an existing data-worth estimation method from the suite of PEST utilities with a global optimization method for optimal sensor placement (called optimal design) in groundwater monitoring networks. Design optimization can include multiple simultaneous sensor locations and multiple sensor types. Both location and sensor type are treated simultaneously as decision variables. Our method combines linear uncertainty quantification and a modified genetic algorithm for discrete multilocation, multitype search. The efficiency of the global optimization is enhanced by an archive of past samples and parallel computing. We demonstrate our methodology for a groundwater monitoring network at the Steinlach experimental site, south-western Germany, which has been established to monitor river-groundwater exchange processes. The target of optimization is the best possible exploration for minimum variance in predicting the mean travel time of the hyporheic exchange. Our results demonstrate that the information gain of monitoring network designs can be explored efficiently and with easily accessible tools prior to taking new field measurements or installing additional measurement points. The proposed methods proved to be efficient and can be applied for model-based optimal design of any type of monitoring network in approximately linear systems. Our key contributions are (1) the use of easy-to-implement tools for an otherwise complex task and (2) yet to consider data-worth interdependencies in simultaneous optimization of multiple sensor locations and sensor types.