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.

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.

Lateral subsurface flow (LSF) is a phenomenon frequently occurring in the field induced by local water saturation along horizon boundaries under nonequilibrium conditions. However, observations of LSF in undisturbed soils under controlled irrigation in the laboratory are limited but needed for model improvement, prediction, and quantification of LSF. We present a method for extracting an undisturbed soil monolith along a soil horizon boundary and introduce an experimental setup for the measurement of LSF and an irrigation device for simulating rainfall. An experimental test run was simulated using HYDRUS 2D. Water infiltrating into the monolith and flowing either laterally along the horizon boundary or vertically through the bottom horizon could be separately captured by suction discs at the side and the bottom. Thus, a clear distinction between lateral and vertical flow was possible. Pressure heads and water contents were recorded by tensiometers and frequency domain reflectometry (FDR) sensors distributed across the monolith in a regular two-dimensional, vertical, cross-sectional pattern. Sensor readings indicated the presence of nonequilibrium conditions within the monolith. Modeling results could reproduce the lateral and vertical outflow of the monolith under constant irrigation, thus showing that water flow within the monolith under steady-state conditions can be explained by the Richards equation and the van Genuchten–Mualem model. The presented method can be used to improve and verify models designed for the prediction of the onset of LSF including that induced by local nonequilibrium conditions.

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.

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.

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.

Summary Understanding the hydrological and hydrogeochemical responses of hillslopes and other small scale groundwater systems requires mapping the velocity and direction of groundwater flow relative to the controlling subsurface material features. Since point observations of subsurface materials and groundwater head are often the basis for modelling these complex, dynamic, three-dimensional systems, considerable uncertainties are inevitable, but are rarely assessed. This study explored whether piezometric head data measured at high spatial and temporal resolution over six years at a hillslope research site provided sufficient information to determine the flow paths that transfer nitrate leached from the soil zone through the shallow saturated zone into a nearby wetland and stream. Transient groundwater flow paths were modelled using MODFLOW and MODPATH, with spatial patterns of hydraulic conductivity in the three material layers at the site being estimated by regularised pilot point calibration using PEST, constrained by slug test estimates of saturated hydraulic conductivity at several locations. Subsequent Null Space Monte Carlo uncertainty analysis showed that this data was not sufficient to definitively determine the spatial pattern of hydraulic conductivity at the site, although modelled water table dynamics matched the measured heads with acceptable accuracy in space and time. Particle tracking analysis predicted that the saturated flow direction was similar throughout the year as the water table rose and fell, but was not aligned with either the ground surface or subsurface material contours; indeed the subsurface material layers, having relatively similar hydraulic properties, appeared to have little effect on saturated water flow at the site. Flow path uncertainty analysis showed that, while accurate flow path direction or velocity could not be determined on the basis of the available head and slug test data alone, the origin of well water samples relative to the material layers and site contour could still be broadly deduced. This study highlights both the challenge of collecting suitably informative field data with which to characterise subsurface hydrology, and the power of modern calibration and uncertainty modelling techniques to assess flow path uncertainty in hillslopes and other small scale systems.

Low-land alluvial gravel aquifers are formed from, and tend to be recharged, by rivers. These interconnected river - groundwater systems can be highly dynamic with groundwater levels following the seasonality of the hydrological regime of the river. The associated groundwater resources are regularly under stress during summer periods when abstractive demand is high and recharge is low. Predicting lead-times for critical groundwater levels allows for a more flexible and adaptive groundwater management. An eigenmodel approach is proposed here as a way of making such predictions, fast and efficiently. The eigenmodel is a mathematical concept that represents the hydraulic function of a groundwater aquifer as a set of conceptual linear reservoirs, arranged in-series. River recharge, land surface recharge, and groundwater abstraction for irrigation are considered as model forcings. The eigenmodel approach is demonstrated on three wells of the unconfined Wairau Aquifer in the Marlborough District of New Zealand, which are used for water resources management. Individual eigenmodels were calibrated to historic data and predictive uncertainty bounds were determined by Markov chain Monte Carlo sampling. Hindcasting of past recession periods showed a low predictive error of the models and a good coverage of the predictive uncertainty bounds. The main advantage of the approach is a 4-orders of magnitude higher computational efficiency compared to a numerical benchmark model. This allows for probabilistic simulation in operational forecasting of groundwater levels. The framework is implemented as a web application for 30-day operational forecasts that comprises automatic data downloads and model input generation, stochastic simulation, uncertainty estimation, visualization, and daily updates on a website.

Summary Groundwater modelers face the challenge of how to assign representative parameter values to the studied aquifer. Several approaches are available to parameterize spatial heterogeneity in aquifer parameters. They differ in their conceptualization and complexity, ranging from homogeneous models to heterogeneous random fields. While it is common practice to invest more effort into data collection for models with a finer resolution of heterogeneities, there is a lack of advice which amount of data is required to justify a certain level of model complexity. In this study, we propose to use concepts related to Bayesian model selection to identify this balance. We demonstrate our approach on the characterization of a heterogeneous aquifer via hydraulic tomography in a sandbox experiment (Illman et al., 2010). We consider four increasingly complex parameterizations of hydraulic conductivity: (1) Effective homogeneous medium, (2) geology-based zonation, (3) interpolation by pilot points, and (4) geostatistical random fields. First, we investigate the shift in justified complexity with increasing amount of available data by constructing a model confusion matrix . This matrix indicates the maximum level of complexity that can be justified given a specific experimental setup. Second, we determine which parameterization is most adequate given the observed drawdown data. Third, we test how the different parameterizations perform in a validation setup. The results of our test case indicate that aquifer characterization via hydraulic tomography does not necessarily require (or justify) a geostatistical description. Instead, a zonation-based model might be a more robust choice, but only if the zonation is geologically adequate.

Bayesian model averaging (BMA) ranks the plausibility of alternative conceptual models according to Bayes' theorem. A prior belief about each model's adequacy is updated to a posterior model probability based on the skill to reproduce observed data and on the principle of parsimony. The posterior model probabilities are then used as model weights for model ranking, selection, or averaging. Despite the statistically rigorous BMA procedure, model weights can become uncertain quantities due to measurement noise in the calibration data set or due to uncertainty in model input. Uncertain weights may in turn compromise the reliability of BMA results. We present a new statistical concept to investigate this weighting uncertainty, and thus, to assess the significance of model weights and the confidence in model ranking. Our concept is to resample the uncertain input or output data and then to analyze the induced variability in model weights. In the special case of weighting uncertainty due to measurement noise in the calibration data set, we interpret statistics of Bayesian model evidence to assess the distance of a model's performance from the theoretical upper limit. To illustrate our suggested approach, we investigate the reliability of soil-plant model selection following up on a study by Wohling et al. (2015). Results show that the BMA routine should be equipped with our suggested upgrade to (1) reveal the significant but otherwise undetected impact of measurement noise on model ranking results and (2) to decide whether the considered set of models should be extended with better performing alternatives.

Study region The Lerma catchment, a small (7.3 km 2 ) sub-catchment of the Ebro Basin in northern Spain. Study focus The Lerma catchment underwent a monitored transition to irrigated agriculture, using water from outside the catchment, between 2006 and 2008. This transition has successfully been simulated using the partial-differential-equation-based model HydroGeoSphere, simulating coupled evapotranspiration, surface water, and groundwater flow in the catchment. We use the calibrated model to study how irrigation practices influence the response of the Lerma catchment to the climate change projected for northern Spain. We consider four different irrigation scenarios: no irrigation, present irrigation, climate-adapted irrigation with current crops, and adapted irrigation for crops requiring less water. The climate scenarios are based on four regional climate models and two downscaling methods. New hydrological insight The simulated catchment responses to climate change show clear differences between the irrigation scenarios. In future climate, groundwater levels and base flows decrease more when irrigation is present than without irrigation, because groundwater levels and base flow in present climate are already at low levels without irrigation. In contrast, annual peak discharges increase more in non-irrigated cases than in irrigated cases. Irrigation increases water availability and an associated rise in potential evapotranspiration results in higher actual evapotranspiration during summer. In non-irrigated scenarios, by contrast, actual evapotranspiration in summer is controlled by precipitation and thus decreases in future climate.

A Bayesian model averaging (BMA) framework is presented to evaluate the worth of different observation types and experimental design options for (1) more confidence in model selection and (2) for increased predictive reliability. These two modeling tasks are handled separately because model selection aims at identifying the most appropriate model with respect to a given calibration data set, while predictive reliability aims at reducing uncertainty in model predictions through constraining the plausible range of both models and model parameters. For that purpose, we pursue an optimal design of measurement framework that is based on BMA and that considers uncertainty in parameters, measurements, and model structures. We apply this framework to select between four crop models (the vegetation components of CERES, SUCROS, GECROS, and SPASS), which are coupled to identical routines for simulating soil carbon and nitrogen turnover, soil heat and nitrogen transport, and soil water movement. An ensemble of parameter realizations was generated for each model using Monte-Carlo simulation. We assess each model's plausibility by determining its posterior weight, which signifies the probability to have generated a given experimental data set. Several BMA analyses were conducted for different data packages with measurements of soil moisture, evapotranspiration (ETa), and leaf area index (LAI). The posterior weights resulting from the different BMA runs were compared to the weight distribution of a reference run with all data types to investigate the utility of different data packages and monitoring design options in identifying the most appropriate model in the ensemble. We found that different (combinations of) data types support different models and none of the four crop models outperforms all others under all data scenarios. The best model discrimination was observed for those data where the competing models disagree the most. The data worth for reducing prediction uncertainty depends on the prediction to be made. LAI data have the highest utility for predicting ETa, while soil moisture data are better for predicting soil water drainage. Our study illustrates, that BMA provides an objective framework for data worth analysis with respect to both model discrimination and model calibration for a wide range of applications.

This study determines the aspects of river bathymetry that have the greatest influence on the predictive biases when simulating hyporheic exchange. To investigate this, we build a highly parameterized HydroGeoSphere model of the Steinlach River Test Site in southwest Germany as a reference. This model is then modified with simpler bathymetries, evaluating the changes to hyporheic exchange fluxes and transit time distributions. Results indicate that simulating hyporheic exchange with a high-resolution detailed bathymetry using a three-dimensional fully coupled model leads to nested multiscale hyporheic exchange systems. A poorly resolved bathymetry will underestimate the small-scale hyporheic exchange, biasing the simulated hyporheic exchange towards larger scales, thus leading to overestimates of hyporheic exchange residence times. This can lead to gross biases in the estimation of a catchment's capacity to attenuate pollutants when extrapolated to account for all meanders along an entire river within a watershed. The detailed river slope alone is not enough to accurately simulate the locations and magnitudes of losing and gaining river reaches. Thus, local bedforms in terms of bathymetric highs and lows within the river are required. Bathymetry surveying campaigns can be more effective by prioritizing bathymetry measurements along the thalweg and gegenweg of a meandering channel. We define the gegenweg as the line that connects the shallowest points in successive cross-sections along a river opposite to the thalweg under average flow conditions. Incorporating local bedforms will likely capture the nested nature of hyporheic exchange, leading to more physically meaningful simulations of hyporheic exchange fluxes and transit times.