Model selection on solid ground: Rigorous comparison of nine ways to evaluate Bayesian model evidence.

Bayesian model selection or averaging objectively ranks a number of plausible, competing conceptual models based on Bayes' theorem. It implicitly performs an optimal trade-off between performance in fitting available data and minimum model complexity. The procedure requires determining Bayesian model evidence (BME), which is the likelihood of the observed data integrated over each model's parameter space. The computation of this integral is highly challenging because it is as high-dimensional as the number of model parameters. Three classes of techniques to compute BME are available, each with its own challenges and limitations: (1) Exact and fast analytical solutions are limited by strong assumptions. (2) Numerical evaluation quickly becomes unfeasible for expensive models. (3) Approximations known as information criteria (ICs) such as the AIC, BIC, or KIC (Akaike, Bayesian, or Kashyap information criterion, respectively) yield contradicting results with regard to model ranking. Our study features a theory-based intercomparison of these techniques. We further assess their accuracy in a simplistic synthetic example where for some scenarios an exact analytical solution exists. In more challenging scenarios, we use a brute-force Monte Carlo integration method as reference. We continue this analysis with a real-world application of hydrological model selection. This is a first-time benchmarking of the various methods for BME evaluation against true solutions. Results show that BME values from ICs are often heavily biased and that the choice of approximation method substantially influences the accuracy of model ranking. For reliable model selection, bias-free numerical methods should be preferred over ICs whenever computationally feasible.

A data-driven approach for modelling Karst spring discharge using transfer function noise models.

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 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., J Hydrol 600:126-508, 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 simulate spring discharge for a previously unseen testing 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 alternative to other approaches that should be considered in future studies.

Robust Data Worth Analysis with Surrogate Models.

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.

Eigenmodels to forecast groundwater levels in unconfined river-fed aquifers during flow recession.

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.

Evaluating Subsurface Parameterization to Simulate Hyporheic Exchange: The Steinlach River Test Site.

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.

An Electron-Balance Based Approach to Predict the Decreasing Denitrification Potential of an Aquifer.

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.

Sensitivity of Simulated Hyporheic Exchange to River Bathymetry: The Steinlach River Test Site.

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.

Quantifying River-Groundwater Interactions of New Zealand's Gravel-Bed Rivers: The Wairau Plain.

New Zealand's gravel-bed rivers have deposited coarse, highly conductive gravel aquifers that are predominantly fed by river water. Managing their groundwater resources is challenging because the recharge mechanisms in these rivers are poorly understood and recharge rates are difficult to predict, particularly under a more variable future climate. To understand the river-groundwater exchange processes in gravel-bed rivers, we investigate the Wairau Plain Aquifer using a three-dimensional groundwater flow model which was calibrated using targeted field observations, "soft" information from experts of the local water authority, parameter regularization techniques, and the model-independent parameter estimation software PEST. The uncertainty of simulated river-aquifer exchange flows, groundwater heads, spring flows, and mean transit times were evaluated using Null-space Monte-Carlo methods. Our analysis suggests that the river is hydraulically perched (losing) above the regional water table in its upper reaches and is gaining downstream where marine sediments overlay unconfined gravels. River recharge rates are on average 7.3 m3 /s, but are highly dynamic in time and variable in space. Although the river discharge regularly hits 1000 m3 /s, the net exchange flow rarely exceeds 12 m3 /s and seems to be limited by the physical constraints of unit-gradient flux under disconnected rivers. An important finding for the management of the aquifer is that changes in aquifer storage are mainly affected by the frequency and duration of low-flow periods in the river. We hypothesize that the new insights into the river-groundwater exchange mechanisms of the presented case study are transferable to other rivers with similar characteristics.

Optimal Design of Multitype Groundwater Monitoring Networks Using Easily Accessible Tools.

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.

Dual-tracer, non-equilibrium mixing cell modelling and uncertainty analysis for unsaturated bromide and chloride transport.

A model is presented for simulating one-dimensional advective dispersive solute transport in the vadose zone. The finite-volume, mixing-cell model uses drainage flux intervals as the index variable, which are calculated by a soil water balance model. The modelling approach considers solute transport from two different regions as well as a slow and a fast transport domain in each region as parallel transport processes. The model is applied to breakthrough curves of Cl(-) and Br(-) measured at different locations and different depths in the volcanic vadose zone of the Tutaeuaua subcatchment of Lake Taupo, New Zealand, following a dual tracer application. Estimates of transport parameter and model predictive uncertainty were derived using the differential evolution adaptive metropolis, DREAM(ZS) adaptive Markov chain Monte Carlo algorithm, a formal Bayesian likelihood function, observed leachate volumes, and Cl(-) breakthrough curves. The model was subsequently evaluated using Br(-) breakthrough curves from the dual tracer experiment and a previously conducted Br(-) tracer-only experiment. Uncertainty bounds derived by this MCMC method simultaneously capture the observed Br(-) and Cl(-) breakthrough curves and corresponding drainage volumes. Results suggest that the slow transport domain properties are relatively similar for different locations in the vadose zone and that the variability in contaminant fluxes is predominantly driven by structural variability of the vadose zone causing lateral flow.