Search

Distribution forecasts Pover future quantities or events are routinely made in hydrology but usually traded for a (likelihood‐weighted) mean or median prediction to accommodate error measures or scoring functionssuch as the mean absolute error or mean squared error. Case in point is the so‐called KG efficiency (KGE) of Gupta et al. (2009, https://doi.org/10.1016/j.jhydrol.2009.08.003) and improvements thereof (Lamontagne et al., 2020, https://doi.org/10.1029/2020wr027101), which have rapidly gained popularity among hydrologists as alternative scoring functionsto the commonly used Nash and Sutcliffe (1970, https://doi.org/10.1016/0022‐1694(70)90255‐6) efficiency, but are equally exclusive in how they quantify model performance using only single‐valued output of the quantities of interest. This point‐valued mapping necessarily implies a loss of information about model performance. This paper advocates the use of probabilistic watershed model training, evaluation and diagnostics. Distribution evaluation opens a mature literature on scoring ruleswhose strong statistical underpinning provides, as we will demonstrate, the theory, context and guidelines necessary for the development of robust information‐theoretically principled metrics for watershed signatures. These so‐called hydrograph functionals are scalar‐valued mappings of major behavioral watershed functions embodied in a strictly properscoring rule. We discuss past developments that led to the current state‐of‐the‐art of distribution evaluation in hydrology and review scoring rules for dichotomous and categorical events, quantiles (intervals) and density forecasts. We are particularly concerned with elicitable functionals and scoring rule propriety, discuss the decomposition of scoring rules into a sharpness, reliability and entropy term and present diagnostically appealing strictly properdivergence scores of hydrograph functionals for flood frequency analysis, flow duration and recession curves. The usefulness and power of distribution‐based model evaluation and diagnostics by means of scoring rules is demonstrated on simple illustrative problems and discharge distributions simulated with watershed models using random sampling and Bayesian model averaging. The presented theory (a) enables a more complete evaluation of distribution forecasts, (b) offers a statistically principled means for watershed model training, evaluation, diagnostics and selection using hydrograph functionals and/or extreme events and (c) provides a universal framework for metric development of watershed signatures, promoting metric standardization and reproducibility. The past decades have witnessed an unbridled growth in goodness‐of‐fit metrics of hydrologic models. These metrics may satisfy the needs of hydrologists but lack conforming theory and principles. This state of affairs (a) elicits improper model training and evaluation, (b) provokes and supports misguided inferences, (c) impedes statistically‐principled uncertainty quantification, metric standardization and development of universal model benchmarks and (d) obfuscates determination of whether the model has finished learning. What is more, most hydrologic model evaluation metrics in use today are rather exclusive in how they quantify model performance using only single‐valued simulated output of the quantities of interest. Predictive distributions derived from (quasi)‐Bayesian methods or ensembles are usually traded for a (likelihood‐weighted) mean or median prediction to accommodate error measures (scoring functions) such as the mean absolute error. This implies a large loss of information. This paper develops a distribution‐based approach to hydrologic model evaluation and diagnostics. Distribution evaluation opens the necessary theory and guidelines for development of robust information‐theoretically principled metrics of watershed signatures. These so‐called hydrograph functionals are scalar‐valued mappings of major behavioral watershed functions embodied in a strictly properscoring rule. The hydrograph functionals offer a statistically principled means for hydrologic model evaluation, diagnostics and selection. Scoring rules of hydrograph functionals provide an information‐theoretically principled means for watershed model training, evaluation, and diagnosticsWe present strictly proper(divergence) scores for flood frequency analysis, flow duration, and recession curvesPropriety and elicitability offer useful working paradigms for metric development of hydrograph functionals Scoring rules of hydrograph functionals provide an information‐theoretically principled means for watershed model training, evaluation, and diagnostics We present strictly proper(divergence) scores for flood frequency analysis, flow duration, and recession curves Propriety and elicitability offer useful working paradigms for metric development of hydrograph functionals

In this Reply, we address the concerns of Knoben and Clark (2023, https://doi.org/10.1029/2022WR034294) or KC23 that “the assumptions needed to effectively use difference‐based variance estimation methods are not always met by hourly streamflow records.” There should be little doubt that the assumptions of our difference‐based estimator will sometimes be violated in hourly streamflow records but the results from de Oliveira and Vrugt (2022, https://doi.org/10.1029/2022wr032263) and confirmed by the findings in our Reply show that such violations are sporadic enough not to affect much the error variance estimates. Snowmelt as pointed out by KC23 (https://doi.org/10.1029/2022WR034294) may not have received sufficient attention in our paper, yet their 365‐day record is simply not long enough to demonstrate bias of our discharge error variance estimates (and their dependence on flow level). This would require analysis of a much longer, multi‐year, streamflow record of a snowmelt‐driven watershed. The snowmelt catchment analyzed in dOV22 (https://doi.org/10.1029/2022wr032263) did not demonstrate bias in the discharge error variance estimates. We also provide additional clarification on the interpretation of the variance estimates obtained with the nonparametric estimator, and discuss the main issues in the test case presented in Knoben and Clark (2023, https://doi.org/10.1029/2022WR034294)). We address the concerns of Knoben and Clark (2023) on the use of the nonparametric estimator with hourly discharge dataWe provide further clarification on the interpretation of the variance estimates obtained with the nonparametric estimator We address the concerns of Knoben and Clark (2023) on the use of the nonparametric estimator with hourly discharge data We provide further clarification on the interpretation of the variance estimates obtained with the nonparametric estimator

Nonstationary extreme value analysis (NEVA) can improve the statistical representation of observed flood peak distributions compared to stationary (ST) analysis, but management of flood risk relies on predictions of out‐of‐sample distributions for which NEVA has not been comprehensively evaluated. In this study, we apply split‐sample testing to 1250 annual maximum discharge records in the United States and compare the predictive capabilities of NEVA relative to ST extreme value analysis using a log‐Pearson Type III (LPIII) distribution. The parameters of the LPIII distribution in the ST and nonstationary (NS) models are estimated from the first half of each record using Bayesian inference. The second half of each record is reserved to evaluate the predictions under the ST and NS models. The NS model is applied for prediction by (1) extrapolating the trend of the NS model parameters throughout the evaluation period and (2) using the NS model parameter values at the end of the fitting period to predict with an updated ST model (uST). Our analysis shows that the ST predictions are preferred, overall. NS model parameter extrapolation is rarely preferred. However, if fitting period discharges are influenced by physical changes in the watershed, for example from anthropogenic activity, the uST model is strongly preferred relative to ST and NS predictions. The uST model is therefore recommended for evaluation of current flood risk in watersheds that have undergone physical changes. Supporting information includes a MATLAB®program that estimates the (ST/NS/uST) LPIII parameters from annual peak discharge data through Bayesian inference. Stationary predictions of flood peak distributions are preferred, overallExtrapolation of the nonstationary model parameter trend rarely improves the stationary prediction, even if an observed trend continuesUsing the most recent nonstationary parameters to predict with an updated stationary model is preferred for physically changing watersheds

What is the “best” model? The answer to this question lies in part in the eyes of the beholder, nevertheless a good model must blend rigorous theory with redeeming qualities such as parsimony and quality of fit. Model selection is used to make inferences, via weighted averaging, from a set of Kcandidate models, Mk;k=(1,…,K), and help identify which model is most supported by the observed data, Y∼=(y∼1,…,y∼n). Here, we introduce a new and robust estimator of the model evidence, p(Y∼|Mk), which acts as normalizing constant in the denominator of Bayes’ theorem and provides a single quantitative measure of relative support for each hypothesis that integrates model accuracy, uncertainty, and complexity. However, p(Y∼|Mk)is analytically intractable for most practical modeling problems. Our method, coined GAussian Mixture importancE (GAME) sampling, uses bridge sampling of a mixture distribution fitted to samples of the posterior model parameter distribution derived from MCMC simulation. We benchmark the accuracy and reliability of GAME sampling by application to a diverse set of multivariate target distributions (up to 100 dimensions) with known values of p(Y∼|Mk)and to hypothesis testing using numerical modeling of the rainfall‐runoff transformation of the Leaf River watershed in Mississippi, USA. These case studies demonstrate that GAME sampling provides robust and unbiased estimates of the evidence at a relatively small computational cost outperforming commonly used estimators. The GAME sampler is implemented in the MATLAB package of DREAM and simplifies considerably scientific inquiry through hypothesis testing and model selection. Science is an iterative process for learning and discovery in which competing ideas about how nature works are evaluated against observations. The translation of each hypothesis to a computational model requires specification of system boundaries, inputs and outputs, state variables, physical/behavioral laws, and material properties; this is difficult and subjective, particularly in the face of incomplete knowledge of the governing spatiotemporal processes and insufficient observed data. To guard against the use of an inadequate model, statisticians advise selecting the “best” model among a set of candidate ones where each might be equally plausible and justifiable a priori. Bayesian model selection uses probability theory to select among competing hypotheses; the key variable is the Bayesian model evidence, which provides a single quantitative measure of relative support for each hypothesis that integrates model accuracy, uncertainty, and complexity. Bayesian model selection has not entered into mainstream use in Earth systems modeling due to the lack of general‐purpose methods to reliably estimate the evidence. Here, we introduce a new method, called GAussian Mixture importancE (GAME) sampling. We demonstrate GAME power and usefulness for hypothesis testing using benchmark experiments with known target and numerical modeling of the rainfall‐runoff transformation of the Leaf River watershed (Mississippi, USA). We introduce a new estimator of Bayesian model evidence to select the best model among candidate onesGAME sampling combines MCMC simulation of the model posterior distribution and bridge samplingWe demonstrate that GAME sampling returns a robust estimate of the Bayesian model evidence under general conditions

The Community Land Model (CLM) contains many parameters whose values are uncertain and thus require careful estimation for model application at individual sites. Here we used Bayesian inference with the DiffeRential Evolution Adaptive Metropolis (DREAM(zs)) algorithm to estimate eight CLM v.4.5 ecosystem parameters using 1 year records of half‐hourly net ecosystem CO2exchange (NEE) observations of four central European sites with different plant functional types (PFTs). The posterior CLM parameter distributions of each site were estimated per individual season and on a yearly basis. These estimates were then evaluated using NEE data from an independent evaluation period and data from “nearby” FLUXNET sites at ~600 km distance to the original sites. Latent variables (multipliers) were used to treat explicitly uncertainty in the initial carbon‐nitrogen pools. The posterior parameter estimates were superior to their default values in their ability to track and explain the measured NEE data of each site. The seasonal parameter values reduced with more than 50% (averaged over all sites) the bias in the simulated NEE values. The most consistent performance of CLM during the evaluation period was found for the posterior parameter values of the forest PFTs, and contrary to the C3‐grass and C3‐crop sites, the latent variables of the initial pools further enhanced the quality‐of‐fit. The carbon sink function of the forest PFTs significantly increased with the posterior parameter estimates. We thus conclude that land surface model predictions of carbon stocks and fluxes require careful consideration of uncertain ecological parameters and initial states. The CLM parameters, estimated separately for four plant functional types, correlated with initial carbon‐nitrogen poolsParameter estimates improved model performance for an independent evaluation period and at independent sitesParameter estimates were more reliable for forest than for C3‐crop and C3‐grass PFTs and less dependent on the initial model states

Variance‐based analysis has emerged as method of choice for quantifying the sensitivity of the output, y, of a scalar‐valued square‐integrable function, f∈ L2(χ$\boldsymbol{\upchi }$), to its d≥ 1 input variables, x= {x1, …, xd}, with support x∈χ⊆Rd$\mathbf{x}\in \boldsymbol{\upchi }\subseteq {\mathbb{R}}^{d}$. The prototype of this approach, Sobol's method is a generalization of the analysis of variance (ANOVA) to d> 2 independent input variables and decomposes y, as sum of elementary functions of zeroth‐, first‐, second‐, up to dth‐order. This independence assumption is mathematically convenient but may not be borne out of the causal or correlational relationships between the x's. This paper is concerned with variance‐based sensitivity analysis (SA) for correlated input variables, for example, multivariate dependencies in a posterior parameter distribution. We use high‐dimensional model representation (HDMR) of Li et al. (2010, https://doi.org/10.1021/jp9096919), Li and Rabitz (2012, https://doi.org/10.1007/s10910-011-9898-0) and replace Sobol's elementary functions with so‐called component functions with unknown expansion coefficients to disentangle the structural, correlative and total contribution of input factors. We contrast the default HDMR methodology with cubic B‐splines and sequential coefficient estimation against its successor, HDMRextof Li and Rabitz (2012, https://doi.org/10.1007/s10910-011-9898-0), which uses polynomial component functions with an extended orthonormalized basis. Benchmark experiments confirm that HDMR and HDMRextparse out the structural and correlative contributions of input factors to the model output and infer an optimal experimental design with parameter correlation. Our last study applies HDMRextto probabilistic SA of a watershed model. The multivariate posterior parameter distribution supports model emulation and yields sensitivity indices that pertain to measured discharge data. We use high‐dimensional model representation for variance‐based sensitivity analysis with dependent input variablesWe use this method to disentangle the parameters' structural and correlative contributions across the posterior distributionCorrelation can profoundly affect input variable importance, optimal experimental design and physical underpinning of sensitivity indices We use high‐dimensional model representation for variance‐based sensitivity analysis with dependent input variables We use this method to disentangle the parameters' structural and correlative contributions across the posterior distribution Correlation can profoundly affect input variable importance, optimal experimental design and physical underpinning of sensitivity indices

The depth to bedrock controls a myriad of processes by influencing subsurface flow paths, erosion rates, soil moisture, and water uptake by plant roots. As hillslope interiors are very difficult and costly to illuminate and access, the topography of the bedrock surface is largely unknown. This essay is concerned with the prediction of spatial patterns in the depth to bedrock (DTB) using high-resolution topographic data, numerical modeling, and Bayesian analysis. Our DTB model builds on the bottom-up control on fresh-bedrock topography hypothesis of Rempe and Dietrich (2014) and includes a mass movement and bedrock-valley morphology term to extent the usefulness and general applicability of the model. We reconcile the DTB model with field observations using Bayesian analysis with the DREAM algorithm. We investigate explicitly the benefits of using spatially distributed parameter values to account implicitly, and in a relatively simple way, for rock mass heterogeneities that are very difficult, if not impossible, to characterize adequately in the field. We illustrate our method using an artificial data set of bedrock depth observations and then evaluate our DTB model with real-world data collected at the Papagaio river basin in Rio de Janeiro, Brazil. Our results demonstrate that the DTB model predicts accurately the observed bedrock depth data. The posterior mean DTB simulation is shown to be in good agreement with the measured data. The posterior prediction uncertainty of the DTB model can be propagated forward through hydromechanical models to derive probabilistic estimates of factors of safety. We introduce an analytic formulation for the spatial distribution of the bedrock depthBayesian analysis reconciles our model with field data and quantifies prediction and parameter uncertaintyThe use of a distributed parameterization recognizes geologic heterogeneities

In this paper, we introduce a relatively simple data‐driven method for the representation of the uncertainty in daily discharge records. The proposed method relies only on hourly discharge data and takes advantage of a nonparametric difference‐based estimator in the characterization of random errors in discharge time series. We illustrate with corrupted streamflow data that the nonparametric estimator provides an accurate characterization of the nature (homoscedastic or heteroscedastic) and magnitude of these errors. In addition, we demonstrate the practical usefulness of the estimator using discharge time series of 500+ watersheds of the Catchment Attributes and MEteorology for Large‐sample Studies data set. This analysis reveals that the magnitude of errors of aleatory nature in the investigated discharge records is rather small (less than 3% for 80% of the records). We then combine the effect of random errors and measurement frequency into a daily variance estimate, which serves as input to a streamflow generation approach. This procedure produces replicates of the discharge record which portray accurately the assigned streamflow uncertainty, preserve key statistical properties of the discharge record and are hydrologically realistic. The proposed method facilitates Bayesian analysis and supports tasks such as model diagnostics, data assimilation, uncertainty quantification and regionalization. We present a model‐agnostic method for the generation of replicates of discharge recordsThe nonparametric estimator provides an accurate estimate of random errors, which for the investigated discharge records are rather smallThe replicates portray accurately the assigned streamflow uncertainty and preserve statistical and hydrologic properties We present a model‐agnostic method for the generation of replicates of discharge records The nonparametric estimator provides an accurate estimate of random errors, which for the investigated discharge records are rather small The replicates portray accurately the assigned streamflow uncertainty and preserve statistical and hydrologic properties

Many watershed models used within the hydrologic research community assume (by default) stationary conditions, that is, the key watershed properties that control water flow are considered to be time invariant. This assumption is rather convenient and pragmatic and opens up the wide arsenal of (multivariate) statistical and nonlinear optimization methods for inference of the (temporally fixed) model parameters. Several contributions to the hydrologic literature have brought into question the continued usefulness of this stationary paradigm for hydrologic modeling. This paper builds on the likelihood‐free diagnostics approach of Vrugt and Sadegh ([Sadegh, M., 2013]) and uses a diverse set of hydrologic summary metrics to test the stationary hypothesis and detect changes in the watersheds response to hydroclimatic forcing. Models with fixed parameter values cannot simulate adequately temporal variations in the summary statistics of the observed catchment data, and consequently, the DREAM(ABC)algorithm cannot find solutions that sufficiently honor the observed metrics. We demonstrate that the presented methodology is able to differentiate successfully between watersheds that are classified as stationary and those that have undergone significant changes in land use, urbanization, and/or hydroclimatic conditions, and thus are deemed nonstationary. Stationarity paradigm is revisited using diagnostic model evaluation with DREAM(ABC)Nonstationarity is not readily apparent from statistical analysisNonstationarity is evident from analysis of summary metrics

Efficient two‐stage MCMC simulation for CPU‐demanding samplingReal‐world application to the inference of a multi‐layered 3D aquifer modelA 2‐5 times observed speed up in sampling efficiency This study reports on two strategies for accelerating posterior inference of a highly parameterized and CPU‐demanding groundwater flow model. Our method builds on previous stochastic collocation approaches, e.g., Marzouk and Xiu (2009) and Marzouk and Najm (2009), and uses generalized polynomial chaos (gPC) theory and dimensionality reduction to emulate the output of a large‐scale groundwater flow model. The resulting surrogate model is CPU efficient and serves to explore the posterior distribution at a much lower computational cost using two‐stage MCMC simulation. The case study reported in this paper demonstrates a two to five times speed‐up in sampling efficiency.

Model to predict Ks from the water retention parametersTreatment of connectivity and tortuosityExploitation of the pore‐size distribution through physically based concepts Prediction of flow through variably saturated porous media requires accurate knowledge of the soil hydraulic properties, namely the water retention function (WRF) and the hydraulic conductivity function (HCF). Unfortunately, direct measurement of the HCF is time consuming and expensive. In this study, we derive a simple closed‐form equation that predicts the saturated hydraulic conductivity, Ksfrom the WRF parameters of Brooks and Corey (1964). This physically based analytical expression uses an empirical tortuosity parameter (τ) and exploits the information embedded in the measured pore‐size distribution. Our proposed model is compared against the current state of the art using more than 250 soil samples from the Grenoble soil catalog (GRIZZLY) and hydraulic properties of European soils (HYPRES) databases. Results demonstrate that the proposed model provides better predictions of the saturated hydraulic conductivity values with reduced size of the 90% confidence intervals of about 3 orders of magnitude.

Nearly 75 years ago, Alfred C. Redfield observed a similarity between the elemental composition of marine plankton in the surface ocean and dissolved nutrients in the ocean interior. This stoichiometry, referred to as the Redfield ratio, continues to be a central tenet in ocean biogeochemistry, and is used to infer a variety of ecosystem processes, such as phytoplankton productivity and rates of nitrogen fixation and loss. Model, field and laboratory studies have shown that different mechanisms can explain both constant and variable ratios of carbon to nitrogen and phosphorus among ocean plankton communities. The range of C/N/P ratios in the ocean, and their predictability, are the subject of much active research. Here we assess global patterns in the elemental composition of phytoplankton and particulate organic matter in the upper ocean, using published and unpublished observations of particulate phosphorus, nitrogen and carbon from a broad latitudinal range, supplemented with elemental data for surface plankton populations. We show that the elemental ratios of marine organic matter exhibit large spatial variations, with a global average that differs substantially from the canonical Redfield ratio. However, elemental ratios exhibit a clear latitudinal trend. Specifically, we observed a ratio of 195:28:1 in the warm nutrient-depleted low-latitude gyres, 137:18:1 in warm, nutrient-rich upwelling zones, and 78:13:1 in cold, nutrient-rich high-latitude regions. We suggest that the coupling between oceanic carbon, nitrogen and phosphorus cycles may vary systematically by ecosystem.

Ground‐penetrating radar (GPR) data affected by waveguide dispersion are not straightforward to analyse. Therefore, waveguide dispersed common midpoint measurements are typically interpreted using so‐called dispersion curves, which describe the phase velocity as a function of frequency. These dispersion curves are typically evaluated with deterministic optimization algorithms that derive the dielectric properties of the subsurface as well as the location and depth of the respective layers. However, these methods do not provide estimates of the uncertainty of the inferred subsurface properties. Here, we applied a formal Bayesian inversion methodology using the recently developed DiffeRential Evolution Adaptive Metropolis DREAM(ZS)algorithm. This Markov Chain Monte Carlo simulation method rapidly estimates the (non‐linear) parameter uncertainty and helps to treat the measurement error explicitly. We found that the frequency range used in the inversion has an important influence on the posterior parameter estimates, essentially because parameter sensitivity varies with measurement frequency. Moreover, we established that the measurement error associated with the dispersion curve is frequency dependent and that the estimated model parameters become severely biased if this frequency‐dependent nature of the measurement error is not properly accounted for. We estimated these frequency‐dependent measurement errors together with the model parameters using the DREAM(ZS)algorithm. The posterior distribution of the model parameters derived in this way compared well with inversion results for a reduced frequency bandwidth which is an alternative, yet subjective method to reduce the bias introduced by this frequency‐dependent measurement error. Altogether, our inversion procedure provides an integrated and objective methodology for the analysis of dispersive GPR data and appropriately treats the measurement error and parameter uncertainty.

The dual‐probe heat‐pulse (DPHP) method is useful for measuring soil thermal properties; however, the probes of a DPHP sensor can deflect when inserted into the soil. Theoretical analysis has shown that measurements of thermal conductivity (λ) should be unaffected by deflection‐induced changes in probe spacing. To verify this result, the conductivities of water, dry sand, and saturated sand were measured using DPHP sensors with probes subject to inward deflection, no deflection, and outward deflection. No error in λ was detected when probes were deflected inward by an amount that caused a 14% reduction in probe spacing. Outward deflection (15% increase in spacing) caused error in λ estimates, but the errors were small (≤0.04 W m−1K−1) and likely to be of little consequence in most applications. We conclude that estimates of λ obtained with the DPHP method are largely unaffected by changes in probe spacing caused by deflection.

In catchments characterized by spatially varying hydrological processes and responses, the optimal parameter values or regions of attraction in parameter space may differ with location-specific characteristics and dominating processes. This paper evaluates the value of semi-distributed calibration parameters for large-scale streamflow simulation using the spatially distributed LISFLOOD model. We employ the Shuffled Complex Evolution Metropolis (SCEM-UA) global optimization algorithm to infer the calibration parameters using daily discharge observations. The resulting posterior parameter distribution reflects the uncertainty about the model parameters and forms the basis for making probabilistic flow predictions. We assess the value of semi-distributing the calibration parameters by comparing three different calibration strategies. In the first calibration strategy uniform values over the entire area of interest are adopted for the unknown parameters, which are calibrated against discharge observations at the downstream outlet of the catchment. In the second calibration strategy the parameters are also uniformly distributed, but they are calibrated against observed discharges at the catchment outlet and at internal stations. In the third strategy a semi-distributed approach is adopted. Starting from upstream, parameters in each subcatchment are calibrated against the observed discharges at the outlet of the subcatchment. In order not to propagate upstream errors in the calibration process, observed discharges at upstream catchment outlets are used as inflow when calibrating downstream subcatchments. As an illustrative example, we demonstrate the methodology for a part of the Morava catchment, covering an area of approximately 10 000 km2. The calibration results reveal that the additional value of the internal discharge stations is limited when applying a lumped parameter approach. Moving from a lumped to a semi-distributed parameter approach: (i) improves the accuracy of the flow predictions, especially in the upstream subcatchments; and (ii) results in a more correct representation of flow prediction uncertainty. The results show the clear need to distribute the calibration parameters, especially in large catchments characterized by spatially varying hydrological processes and responses.

Inverse modeling has become increasingly popular for estimating effective hydraulic properties across a range of spatial scales. In recent years, many different algorithms have been developed to solve complex multiobjective optimization problems. In this study, we compared the efficiency of the Nondominated Sorting Genetic Algorithm (NSGA‐II), the Multiobjective Shuffled Complex Evolution Metropolis algorithm (MOSCEM‐UA), and AMALGAM, a multialgorithm genetically adaptive search method for multiobjective estimation of soil hydraulic parameters. In our analyses, we implemented the HYDRUS‐1D model and used observed pressure head data at three different depths from the Spydia experimental field site in New Zealand. Our optimization problem was posed in a multiobjective context by simultaneously using three complementary RMSE criteria at each depth. We analyzed the trade‐off between these criteria and the adherent Pareto uncertainty. The results demonstrate that all three algorithms were able to find a good approximation of the Pareto set of solutions, but differed in the rate of convergence to this distribution. Small differences in performance of the various algorithms were observed because of the relative high dimension of the optimization problem in combination with the presence of multiple local optimal solutions within the three‐objective search space. The Pareto parameter sets yielded satisfactory results when simulating the transient tensiometric pressure at predetermined observation points in the investigated vadose zone profile. The overall best parameter set was found by AMALGAM with RMSE values of 0.14, 0.11, and 0.17 m at the 0.4‐, 1.0‐, and 2.6‐m depths, respectively. In contrast, the fit errors were substantially higher at these respective depths, with RMSE values ranging from 0.87 to 1.49 m, when using soil hydraulic parameters derived from laboratory analysis of small vadose zone cores.

Model nonlinearity and parameter interdependence violate the use of a first‐order approximation to obtain exact confidence intervals of parameters in soil hydrologic models. In this study, the posterior distribution of parameters in soil water retention and hydraulic conductivity functions is examined using observed water retention data and a laboratory transient multistep outflow experiment. Parameter uncertainties obtained with traditional first‐order approximations and uniform grid sampling strategies were compared with those obtained using the Metropolis algorithm, a Markov Chain Monte Carlo (MCMC) sampler. A diagnostic measure, based on multiple sequences generated in parallel, was used to test whether convergence of the Metropolis sampler to the posterior distribution had been achieved. Most significantly, as the Metropolis algorithm can cope with rough response surfaces generated by the objective function used, it not only successfully infers the multivariate posterior probability distribution of the model parameters, but also provides valuable insights in parameter interdependence in the full parameter space.