Search

Knowledge of the recurrence intervals of precipitation extremes is vital for infrastructure design, risk assessment, and insurance planning. However, trends and shifts in rainfall patterns globally pose challenges to the application of extreme value analysis (EVA) which relies critically on the assumption of stationarity. In this paper, we explore: (1) the suitability of nonstationary (NS) models in the presence of statistically significant trends, and (2) their potential in modeling out-of-sample data to improve frequency analysis of extreme precipitation. We analyze the benefits of using a nonstationary Generalized Extreme Value (GEV) model for annual extreme precipitation records from Southern Brazil. The location of the GEV distribution is allowed to change with time. The unknown GEV model parameters are estimated using Bayesian techniques coupled with Markov chain Monte Carlo simulation. Next, we use GAME sampling to compute the evidence (and their ratios, the so-called Bayes factors) for stationary and nonstationary models of annual maximum precipitation. Our results show that the presence of a statistically significant trend in annual maximum precipitation alone does not justify the use of a NS model. The location parameter of the GEV distribution must also be well defined, otherwise, stationary models of annual maximum precipitation receive more support by the data. These findings reiterate the importance of accounting for GEV model parameters and predictive uncertainty in frequency analysis and hypothesis testing of annual maximum precipitation data records. Furthermore, a meaningful EVA demands detailed knowledge about the origin and persistence of observed changes.

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

The two-term infiltration equation (Formula presented.) is commonly used to determine the sorptivity, (Formula presented.) (Formula presented.), and product, (Formula presented.) (Formula presented.), of the dimensionless multiple (Formula presented.) and saturated soil hydraulic conductivity (Formula presented.) (Formula presented.) from cumulative vertical infiltration measurements (Formula presented.) (L) at times (Formula presented.) (T). This reduced form of the quasi-analytical power series solution of Richardson's equation of Philip enjoys a solid physical underpinning but at the expense of a limited time validity. Using simulated infiltration data, Jaiswal et al. have shown this time validity to equal about 2.5 cm of cumulative infiltration. The goals of this work are twofold. First, we investigate the extent to which cumulative infiltration measurements larger than 2.5 cm bias the estimates of (Formula presented.) and (Formula presented.). Second, we investigate the impact of epistemic errors on the inferred time validities and parameters. Partial infiltration curves up to 2.5 cm of cumulative vertical infiltration improve substantially the agreement between actual and least squares estimates of (Formula presented.) and (Formula presented.). But this only holds if the data generating infiltration process follows Richardson's equation and experimental conditions satisfy assumptions of soil homogeneity and a uniform initial water content. Otherwise, autocorrelated cumulative infiltration residuals will bias the least squares estimates of (Formula presented.) and (Formula presented.). Our findings reiterate and reinvigorate earlier conclusions of Haverkamp et al. and show that epistemic errors deteriorate the physical significance of the coefficients of infiltration functions. As a result, the parameters of infiltration functions cannot simply be used in storm water and vadose zone flow models to forecast runoff and recharge at field and landscape scales unless these predictions are accompanied by realistic uncertainty bounds. We conclude that the time validity of Philip's two-term equation is an elusive theoretical quantity with arbitrary physical meaning.

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((Formula presented.)), to its d ≥ 1 input variables, x = {x1, …, xd}, with support (Formula presented.). 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, HDMRext of 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 HDMRext parse 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 HDMRext to probabilistic SA of a watershed model. The multivariate posterior parameter distribution supports model emulation and yields sensitivity indices that pertain to measured discharge data.

This paper is concerned with the formulation of an adequate likelihood function in the application of Bayesian epistemology to uncertainty quantification of hydrologic models. We focus our attention on a special class of likelihood functions (hereinafter referred to as distribution-adaptive likelihood functions), which do not require prior assumptions about the expected distribution of the residuals, rather inference takes place over the hypotheses (model parameters) and space of distribution functions. Our goals are threefold. First, we present theory of a revised implementation of the generalized likelihood (GL) function of Schoups and Vrugt (2010) wherein residual standardization precedes the treatment of serial correlation. This so-called GL+ function, enjoys a solid statistical underpinning and guarantees a more robust joint inference of the autoregressive coefficients and residual properties. Then, as secondary goal, we present a further generalization of the GL+ function, coined the universal likelihood (UL) function, which extends applicability to highly asymmetrical lepto- and platy-kurtic residual distributions. The UL function builds on the 5-parameter skewed generalized Student's t distribution of Theodossiou (2015) which makes up a large family of continuous probability distributions including (but not limited to) symmetric and skewed forms of the generalized normal, generalized t, Laplace, normal, Student's t, and Cauchy-Lorentz distributions. As our third and last goal, we present the use of strictly proper scoring rules to evaluate, compare and rank likelihood functions. These scoring rules condense the accuracy of a distribution forecast to a single value while retaining attractive statistical properties. The GL+ and UL functions are illustrated using data of a simple autoregressive scheme and benchmarked against the GL function, Student t likelihood (SL) of Scharnagl et al. (2015) and normal likelihood (NL) for a conceptual hydrologic model using measured streamflow data. Our results show that, (i) the GL+ function is superior to the GL function, (ii) the active set of nuisance variables exerts a large control on the performance of the GL+, SL and UL functions, (iii) the treatment of autocorrelation deteriorates the scoring rules and performance metrics of the forecast distribution, (iv) a leptokurtic distribution is favored for discharge residuals, (v) scoring rules are indispensable in our search for the true forecast distribution, and (vi) the use of multiple strictly proper scoring rules turns the selection of an adequate likelihood function into a multi-criteria problem.

Evolutionary Polynomial Regression (EPR) has found widespread application and use for model structure development in engineering and science. This hybrid evolutionary approach merges real world data and explanatory variables to generate well-structured models in the form of polynomial equations. The simple and transparent models produced by this technique enable us to explore, via sensitivity analysis, the robustness of the derived models. Yet, existing EPR frameworks do not make explicit use of sensitivity analysis in the selection of robust and high-fidelity model structures. In this paper, we develop a multi-step sensitivity-driven method which combines the strengths of differential evolution and model selection via Monte Carlo simulation to explore the input–output relationships of model structures. In the first step, our hybrid approach automatically determines the optimum number of terms of the polynomial equations. In a subsequent step, our algorithm explores the mean parametric response of each explanatory variable used in the mathematical formulation to select a final model structure. Finally, in our selection of the most robust mathematical structure, we take explicit consideration of the prediction uncertainty of the simulated output. We illustrate and evaluate our EPR method for different engineering problems involving modeling and prediction of the moisture content and creep index of soils. Altogether, our results demonstrate that the use of sensitivity analysis as an integral part of model structure search and selection will lead to robust models with high predictive ability.

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.

The Kling-Gupta efficiency, hereafter referred to as KG efficiency rather than its common abbreviation KGE, proposed by Gupta et al. (2009) has become a widely used metric for evaluating the goodness-of-fit of n-vectors of observations, ỹ=[ỹ1ỹ2…ỹn]⊤, and corresponding model simulations, yθ=y1θy2θ…ynθ⊤. This metric rectifies some of the shortcomings of the coefficient of determination, R2, also known to hydrologists as the efficiency of Nash and Sutcliffe (1970), by using a Euclidean-distance based weighting of the correlation, bias and temporal variability of the observed, ỹ, and simulated, yθ, data. But as the KG efficiency is not borne out of assumptions with respect to the statistical distribution of the residuals, eθ=ỹ-yθ, we cannot formally characterize its uncertainty. The NS efficiency suffers a similar problem, yet, statistical theory postulates that its confidence intervals should follow a beta distribution in certain special cases. Without a formal description of the confidence intervals of the KG efficiency, we cannot (amongst others) quantify parameter uncertainty, compute confidence and prediction limits on simulated model responses, inform decision makers about critical modeling uncertainties, evaluate model adequacy and assess the information content of calibration data. More fundamentally, without confidence intervals we cannot establish whether the KG efficiency is a consistent, efficient and unbiased estimator. In this paper we present an empirical description of the confidence intervals of the KG efficiency. We relate the unknown probability distribution of the KG efficiency to the measurement errors of the training data record, ỹ, and use the bootstrap method to carry out statistical inference. We illustrate our method by application to a simple linear regression function for which the least squares parameter confidence regions are exactly known and two hydrologic models of contrasting complexity. The empirical parameter confidence regions and/or intervals of the KG efficiency are compared to those derived from generalized least squares, objective function contouring and Bayesian analysis using Markov chain Monte Carlo simulation. The marginal parameter distributions of the KG efficiency are generally well described by a normal distribution. Results further confirm that the distribution of the KG efficiency is a complex function of data length and the magnitude, distribution and structure of the measurement errors. This prohibits an analytic description of the empirical confidence regions and/or intervals of the KG efficiency and reiterates the need for the bootstrap method.

Stormwater capture systems have the potential to address many urban stormwater management challenges, particularly in water-scarce regions like Southern California. Here, we investigate the potential best-case limits of water supply and stormwater retention benefits delivered by a 10,000 m3 stormwater capture system equipped with real-time control (RTC) on a university campus in Southern California. Using a copula-based conditional probability analysis, two performance metrics (percent of water demand satisfied and the percent of stormwater runoff captured) are benchmarked relative to (1) precipitation seasonality (historical rainfall and a counterfactual in which the same average annual rainfall is distributed evenly over the year); (2) annual precipitation (dry, median, and wet years); and (3) three RTC algorithms (no knowledge of future rainfall or perfect knowledge of future rainfall 1 or 2 days in advance). RTC improves stormwater retention, particularly for the highly seasonal rainfall patterns in Southern California, but not water supply. Improvements to the latter will likely require implementing stormwater capture RTC in conjunction with other stormwater infrastructure innovations, such as spreading basins for groundwater recharge and widespread adoption of green stormwater infrastructure.

Soil water infiltration is central to hydrologic studies and lends itself for detailed experimentation and mathematical–physical modeling. The most rigorous of these approaches numerically solve Richards’ equation, possibly coupled through a heterogeneous soil to a surface water routine and groundwater model. This approach is computationally expensive and prone to mass balance errors and overparameterization. Analytic solutions of the infiltration process obviate the need for specification of hydraulic functions and simplify computation and inverse determination of soil properties. This paper investigates the usefulness of Parlange's three-parameter infiltration equation for forward and inverse modeling of vertical infiltration experiments. This quasi-exact implicit solution of Richards’ equation, also credited to Haverkamp and coworkers in recent literature, is valid for the entire infiltration event, matches cumulative infiltration data from different soils and its (super-)parameters, S, Ks, and β, exhibit a solid mathematical–physical underpinning. Nonetheless, Parlange's equation has not entered mainstream use for infiltration simulation and data analysis in the absence of a robust, exact and efficient numerical solution. This paper builds upon the recent work of Jaiswal et al. and presents theory, algorithms, and source codes of two numerical procedures for forward and inverse modeling of Parlange's infiltration equation. We illustrate the procedures using measured infiltration data from the Soil Water Infiltration Global (SWIG) database. Our findings highlight the potential of Parlange's equation for infiltration modeling and hydraulic characterization of the soil sorptivity, S [L T−1/2], saturated hydraulic conductivity, Ks [L T−1], and unitless coefficient, β. Parlange's infiltration equation provides a powerful alternative to mathematically more convenient explicit infiltration equations that suffer physical underpinning and/or a limited time validity.

Many different equations have been proposed to describe quantitatively one-dimensional soil water infiltration. The unknown coefficients of these equations characterize soil hydraulic properties and may be estimated from a n record, (Formula presented.), of cumulative infiltration measurements using curve fitting techniques. The two-term infiltration equation, (Formula presented.), of Philip has been widely used to describe measured infiltration data. This function enjoys a solid mathematical–physical underpinning and admits a closed-form solution for the soil sorptivity, S [L T−1/2], and multiple, c [−], of the saturated hydraulic conductivity, Ks [L T−1]. However, Philip's two-term equation has a limited time validity, tvalid [T], and thus cumulative infiltration data, (Formula presented.), beyond (Formula presented.) will corrupt the estimates of S and Ks. This paper introduces a novel method for estimating S, c, Ks, and tvalid of Philip's two-term infiltration equation. This method, coined parasite inversion, use as vehicle Parlange's three-parameter infiltration equation. As prerequisite to our method, we present as secondary contribution an exact, robust and efficient numerical solution of Parlange's infiltration equation. This solution admits Bayesian parameter estimation with the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm and yields as byproduct the marginal distribution of Parlange's β parameter. We evaluate our method for 12 USDA soil types using synthetic infiltration data simulated with HYDRUS-1D. An excellent match is observed between the inferred values of S and Ks and their “true” values known beforehand. Furthermore, our estimates of c and tvalid correlate well with soil texture, corroborate linearity of the (Formula presented.) relationship for (Formula presented.), and fall within reported ranges. A cumulative vertical infiltration of about 2.5 cm may serve as guideline for the time-validity of Philip's two-term infiltration equation.

Strength properties of most sedimentary and metamorphic rocks are known to vary with direction. Knowledge of this so-called rock anisotropy is of utmost importance for reliability analysis and engineering design. The purpose of this paper is twofold. First, we propose a formulation of the Hoek–Brown (HB) failure criterion, which calculates strength anisotropy using a non-uniform scaling of the stress tensor. We use two scaling factors, CN and CS, to link the orientation of the anisotropy planes with the loading direction. As we assume isotropic parameters for intact rock, our HB model formulation is relatively easy to use and has the additional advantage that it does not demand any modifications to the HB failure criterion. Second, we embed our HB model formulation in a Bayesian framework and illustrate its power and usefulness using experimental data of anisotropic rock samples published in the literature. Results demonstrate that our HB model formulation predicts accurately measured peak strengths of rocks with different degrees of anisotropy, confining stresses and anisotropy orientations. The uncertainty in peak strength of anisotropic rocks can be quite large, reiterating the need for an explicit treatment of strength anisotropy uncertainty in rock mechanics studies. The Bayesian methodology is general-purpose, and, as such, can help better inform geotechnical engineers, contractors and other professionals about rock conditions and design reliability and assist decision makers in determining the overall risks of engineering structures.

Many different equations ranging from simple empirical to semi-analytical solutions of the Richards equation have been proposed for quantitative description of water infiltration into variably saturated soils. The sorptivity, S, and the saturated hydraulic conductivity, Ks, in these equations are typically unknown and have to be estimated from measured data. In this paper, we use so-called characteristic time (tchar) to design a new method, referred to as the characteristic time method (CTM) that estimates S, and Ks, from one-dimensional (1D) cumulative infiltration data. We demonstrate the usefulness and power of the CTM by comparing it with a suite of existing methods using synthetic cumulative infiltration data simulated by HYDRUS-1D for 12 synthetic soils reflecting different USDA textural classes, as well as experimental data selected from the Soil Water Infiltration Global (SWIG) database. Results demonstrate that the inferred values of S and Ks are in excellent agreement with their theoretical values used in the synthetically simulated infiltration experiments with Nash–Sutcliffe criterion close to unity and RMSE values of 0.04 cm h−1/2 and 0.05 cm h−1, respectively. The CTM also showed very high accuracy when applied on synthetic data with added measurement noise, as well as robustness when applied to experimental data. Unlike previously published methods, the CTM does not require knowledge of the time validity of the applied semi-analytical solution for infiltration and, therefore, is applicable to infiltrations with durations from 5 min to several days. A script written in Python of the CTM method is provided in the supplemental material.

In the past decades, many different approaches have been developed in the literature to quantify the load-carrying capacity and geotechnical stability (or the Factor of Safety, F s) of variably saturated hillslopes. Much of this work has focused on a deterministic characterization of hillslope stability. Yet, simulated F s values are subject to considerable uncertainty due to our inability to characterize accurately the soil mantle's properties (hydraulic, geotechnical and geomorphologic) and spatiotemporal variability of the moisture content of the hillslope interior. This is particularly true at larger spatial scales. Thus, uncertainty-incorporating analyses of physically based models of rain-induced landslides are rare in the literature. Such landslide modeling is typically conducted at the hillslope scale using gauge-based rainfall forcing data with rather poor spatiotemporal coverage. For regional landslide modeling, the specific advantages and/or disadvantages of gauge-only, radar-merged and satellite-based rainfall products are not clearly established. Here, we compare and evaluate the performance of the Transient Rainfall Infiltration and Grid-based Regional Slope-stability analysis (TRIGRS) model for three different rainfall products using 112 observed landslides in the period between 2004 and 2011 from the North Carolina Geological Survey database. Our study includes the Tropical Rainfall Measuring Mission (TRMM) Multi-satellite Precipitation Analysis Version 7 (TMPA V7), the North American Land Data Assimilation System Phase 2 (NLDAS-2) analysis, and the reference 'truth' Stage IV precipitation. TRIGRS model performance was rather inferior with the use of literature values of the geotechnical parameters and soil hydraulic properties from ROSETTA using soil textural and bulk density data from SSURGO (Soil Survey Geographic database). The performance of TRIGRS improved considerably after Bayesian estimation of the parameters with the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm using Stage IV precipitation data. Hereto, we use a likelihood function that combines binary slope failure information from landslide event and 'null' periods using multivariate frequency distribution-based metrics such as the False Discovery and False Omission Rates. Our results demonstrate that the Stage IV-inferred TRIGRS parameter distributions generalize well to TMPA and NLDAS-2 precipitation data, particularly at sites with considerably larger TMPA and NLDAS-2 rainfall amounts during landslide events than null periods. TRIGRS model performance is then rather similar for all three rainfall products. At higher elevations, however, the TMPA and NLDAS-2 precipitation volumes are insufficient and their performance with the Stage IV-derived parameter distributions indicate their inability to accurately characterize hillslope stability.

Hydroelectric power is a relatively cheap, reliable, sustainable, and renewable source of energy that can be generated without toxic waste and considerably lower emissions of greenhouse gases than fossil fuel energy plants. Conventional hydroelectric plants produce energy by the controlled release of dammed reservoir water to one or more turbines via a penstock. The kinetic energy of the falling water produces a rotational motion of the turbine shaft and this mechanical energy is converted into electricity via a power generator. Dam-based plants are among the largest and most flexible power producing facilities in the world, yet their construction and operation is costly and can damage and disrupt upstream and downstream ecosystems and have catastrophic effects on downriver settlements and infrastructure. Run-of-the-river (RoR) hydroelectric stations are an attractive and environmentally friendly alternative to dam-based facilities. These plants divert water from a flowing river to a turbine and do not require the formation of a reservoir. Despite their minimal impact on the surrounding environment and communities, the potential of RoR plants has not been fully explored and exploited. For example, in the United States it is estimated that RoR plants could annually produce 60,000 MW, or about 13% of the total electricity consumption in 2016. Here, we introduce a numerical model, called HYdroPowER or HYPER, which uses a daily time step to simulate the technical performance, energy production, maintenance and operational costs, and economic profit of a RoR plant in response to a suite of different design and construction variables and record of river flows. The model is coded in MATLAB and includes a built-in evolutionary algorithm that enables the user to maximize the RoR plant's power production or net economic profit by optimizing (among others) the penstock diameter, and the type (Kaplan, Francis, Pelton and Crossflow) design flow, and configuration (single/parallel) of the turbine system. Unlike other published models, this module of HYPER carefully considers each turbine's design flow, admissible suction head, specific and rotational speed in evaluating the technical performance, cost and economic profit of a RoR plant. Two case studies illustrate the power and practical applicability of HYPER. Some of their results confirm earlier literature findings, that (I) the optimum capacity and design flow of a RoR plant is controlled by the river's flow duration curve, (II) a highly variable turbine inflow compromises energy production, and (III) a side-by-side dual turbine system enhances considerably the range of workable flows, operational flexibility and energy production of a RoR plant. HYPER includes a GUI and is available upon request from the authors.

Vegetation plays an important role in regulating global carbon cycles and is a key component of the Earth system models (ESMs) that aim to project Earth's future climate. In the last decade, the vegetation component within ESMs has witnessed great progress from simple "big-leaf" approaches to demographically structured approaches, which have a better representation of plant size, canopy structure, and disturbances. These demographically structured vegetation models typically have a large number of input parameters, and sensitivity analysis is needed to quantify the impact of each parameter on the model outputs for a better understanding of model behavior. In this study, we conducted a comprehensive sensitivity analysis to diagnose the Community Land Model coupled to the Functionally Assembled Terrestrial Simulator, or CLM4.5(FATES). Specifically, we quantified the first- and second-order sensitivities of the model parameters to outputs that represent simulated growth and mortality as well as carbon fluxes and stocks for a tropical site with an extent of 1×1°. While the photosynthetic capacity parameter (Vc;max25) is found to be important for simulated carbon stocks and fluxes, we also show the importance of carbon storage and allometry parameters, which determine survival and growth strategies within the model. The parameter sensitivity changes with different sizes of trees and climate conditions. The results of this study highlight the importance of understanding the dynamics of the next generation of demographically enabled vegetation models within ESMs to improve model parameterization and structure for better model fidelity.

Low frequency loop-loop electromagnetic induction (EMI) is a widely-used geophysical measurement method to rapidly measure in situ the apparent electrical conductivity (ECa) of variably-saturated soils. Here, we couple Bayesian inversion of a quasi-two-dimensional electromagnetic (EM) model with image compression via the discrete cosine transform (DCT) for subsurface electrical conductivity (EC) imaging. The subsurface EC distributions are obtained from multi-configuration EMI data measured with a CMD-Explorer sensor along two transects in the Chicken Creek catchment (Brandenburg, Germany). Dipole-dipole electrical resistivity tomography (ERT) data are used to benchmark the inferred EC fields of both transects. We are especially concerned with the impact of the DCT truncation method on the accuracy and reliability of the inversely-estimated EC images. We contrast the results of two different truncation approaches for model parametrization. The first scenario considers an arbitrary selection of the dominant DCT coefficients and their prior distributions (a commonly-used approach), while the second methodology benefits from geostatistical simulation of the EMI data pseudosection. This study demonstrates that DCT truncation based on geostatistical simulations facilitates a robust selection of the dominant DCT coefficients and their prior ranges, resulting in more accurate subsurface EC imaging from multi-configuration EMI data. Results based on geostatistical prior modeling present an excellent agreement between the EMI- and ERT-derived EC fields of the Chicken Creek catchment.

Crosshole ground-penetrating radar (GPR) is a widely used measurement technique to help inspect the structural integrity of man-made underground structures. In a previous paper, we have introduced a Bayesian framework for inversion of crosshole GPR experiments to help back out defects in concrete underground structures. Here, we evaluate the practical usefulness of our inversion framework by application to waveform data from a real-world GPR survey of a diaphragm wall panel with two embedded structure defects. We also use this case study to further refine our methodology by introducing the elements of a two-stage inversion method to help delineate the exact location and shape of small structure defects. Herein, a low-resolution inversion composed of relatively few inversion coefficients (stage-1) is used to determine roughly the presence of structure defects, followed by a second inversion (stage-2) with much enhanced spatial resolution in those areas classified with anomalous or suspicious permittivity values. This two-stage inversion approach uses more wisely CPU-resources by focusing primarily on those areas of the concrete structure that have been classified as anomalies. We investigate the benefits of this two-stage inversion scheme using a synthetic and real-world case study involving waveform data of a diaphragm wall panel measured with crosshole GPR. Our results demonstrate that the proposed two-stage inversion method recovers successfully the location and shape of structure defects, at a computational cost that is considerably lower than the original inversion framework.

Land transactions between farmers are responsible for landscape changes in rural areas. The price a farmer is willing to pay (WTP) for vacant land depends on the distance of the parcel to the farmstead. Detailed quantitative knowledge of this WTP– distance relationship is of utmost importance for accurate modelling of land markets, and for the design and implementation of effective and robust land consolidation schemes. Practical experience suggests, however that it is not particularly easy to back out the WTP–distance relationship from empirical transaction data. Here, we present a novel statistical framework to help quantify the relationship between a farmer's WTP and the distance of his/her farmstead to the vacant parcel. We describe a land market with a simple statistical model and simulate an artificial archive of land transactions via Monte Carlo sampling. The parameters of our virtual market are estimated from a historical archive of land transactions in the Province of Gelderland The Netherlands, using minimization of the divergence (relative entropy) between the observed and simulated joint distributions of distance and transaction price. A reasonable agreement was observed between the observed and simulated bivariate distributions of distance and transaction price. Our results demonstrate that for short distances (500–1000m) any additional metre distance reduces the WTP by about 60 € ha-1. The impact of distance on WTP gradually levels off with larger distance: beyond 5 km the effect has reduced to less than 0.5 € ha-1

This essay illustrates some recent developments to the DiffeRential Evolution Adaptive Metropolis (DREAM) MATLAB toolbox of Vrugt (2016) to delineate and sample the behavioural solution space of set-theoretic likelihood functions used within the GLUE (Limits of Acceptability) framework (Beven and Binley, 1992, 2014; Beven and Freer, 2001; Beven, 2006). This work builds on the DREAM(ABC) algorithm of Sadegh and Vrugt (2014) and enhances significantly the accuracy and CPU-efficiency of Bayesian inference with GLUE. In particular it is shown how lack of adequate sampling in the model space might lead to unjustified model rejection.

Summary:Biological models contain many parameters whose values are difficult to measure directly via experimentation and therefore require calibration against experimental data. Markov chain Monte Carlo (MCMC) methods are suitable to estimate multivariate posterior model parameter distributions, but these methods may exhibit slow or premature convergence in high-dimensional search spaces. Here, we present PyDREAM, a Python implementation of the (Multiple-Try) Differential Evolution Adaptive Metropolis [DREAM(ZS)] algorithm developed by Vrugt and ter Braak (2008) and Laloy and Vrugt (2012). PyDREAM achieves excellent performance for complex, parameter-rich models and takes full advantage of distributed computing resources, facilitating parameter inference and uncertainty estimation of CPU-intensive biological models. Availability and implementation:PyDREAM is freely available under the GNU GPLv3 license from the Lopez lab GitHub repository at http://github.com/LoLab-VU/PyDREAM. Contact:c.lopez@vanderbilt.edu. Supplementary information:Supplementary data are available at Bioinformatics online.

The daunting complexity of ecosystems has led ecologists to use mathematical modelling to gain understanding of ecological relationships, processes and dynamics. In pursuit of mathematical tractability, these models use simplified descriptions of key patterns, processes and relationships observed in nature. In contrast, ecological data are often complex, scale-dependent, space-time correlated, and governed by nonlinear relations between organisms and their environment. This disparity in complexity between ecosystem models and data has created a large gap in ecology between model and data-driven approaches. Here, we explore data assimilation (DA) with the Ensemble Kalman filter to fuse a two-predator-two-prey model with abundance data from a 2600+ day experiment of a plankton community. We analyse how frequently we must assimilate measured abundances to predict accurately population dynamics, and benchmark our population model's forecast horizon against a simple null model. Results demonstrate that DA enhances the predictability and forecast horizon of complex community dynamics.

Abstract. Land surface models (LSMs) use a large cohort of parameters and state variables to simulate the water and energy balance at the soil–atmosphere interface. Many of these model parameters cannot be measured directly in the field, and require calibration against measured fluxes of carbon dioxide, sensible and/or latent heat, and/or observations of the thermal and/or moisture state of the soil. Here, we evaluate the usefulness and applicability of four different data assimilation methods for joint parameter and state estimation of the Variable Infiltration Capacity Model (VIC-3L) and the Community Land Model (CLM) using a 5-month calibration (assimilation) period (March–July 2012) of areal-averaged SPADE soil moisture measurements at 5, 20, and 50 cm depths in the Rollesbroich experimental test site in the Eifel mountain range in western Germany. We used the EnKF with state augmentation or dual estimation, respectively, and the residual resampling PF with a simple, statistically deficient, or more sophisticated, MCMC-based parameter resampling method. The performance of the calibrated LSM models was investigated using SPADE water content measurements of a 5-month evaluation period (August–December 2012). As expected, all DA methods enhance the ability of the VIC and CLM models to describe spatiotemporal patterns of moisture storage within the vadose zone of the Rollesbroich site, particularly if the maximum baseflow velocity (VIC) or fractions of sand, clay, and organic matter of each layer (CLM) are estimated jointly with the model states of each soil layer. The differences between the soil moisture simulations of VIC-3L and CLM are much larger than the discrepancies among the four data assimilation methods. The EnKF with state augmentation or dual estimation yields the best performance of VIC-3L and CLM during the calibration and evaluation period, yet results are in close agreement with the PF using MCMC resampling. Overall, CLM demonstrated the best performance for the Rollesbroich site. The large systematic underestimation of water storage at 50 cm depth by VIC-3L during the first few months of the evaluation period questions, in part, the validity of its fixed water table depth at the bottom of the modeled soil domain.

We investigate the uncertainty in bedrock depth and soil hydraulic parameters on the stability of a variably-saturated slope in Rio de Janeiro, Brazil. We couple Monte Carlo simulation of a three-dimensional flow model with numerical limit analysis to calculate confidence intervals of the safety factor using a 22-day rainfall record. We evaluate the marginal and joint impact of bedrock depth and soil hydraulic uncertainty. The mean safety factor and its 95% confidence interval evolve rapidly in response to the storm events. Explicit recognition of uncertainty in the hydraulic properties and depth to bedrock increases significantly the probability of failure.

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.

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 K candidate models, (Formula presented.), and help identify which model is most supported by the observed data, (Formula presented.). Here, we introduce a new and robust estimator of the model evidence, (Formula presented.), 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, (Formula presented.) 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 (Formula presented.) 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.

Geophysical data can help to discriminate among multiple competing subsurface hypotheses (conceptual models). Here, we explore the merits of Bayesian model selection in hydrogeophysics using crosshole ground-penetrating radar data from the South Oyster Bacterial Transport Site in Virginia, USA. Implementation of Bayesian model selection requires computation of the marginal likelihood of the measured data, or evidence, for each conceptual model being used. In this paper, we compare three different evidence estimators, including (1) the brute force Monte Carlo method, (2) the Laplace-Metropolis method, and (3) the numerical integration method proposed by Volpi et al. (2016). The three types of subsurface models that we consider differ in their treatment of the porosity distribution and use (a) horizontal layering with fixed layer thicknesses, (b) vertical layering with fixed layer thicknesses and (c) a multi-Gaussian field. Our results demonstrate that all three estimators provide equivalent results in low parameter dimensions, yet in higher dimensions the brute force Monte Carlo method is inefficient. The isotropic multi-Gaussian model is most supported by the travel time data with Bayes factors that are larger than 10100 compared to conceptual models that assume horizontal or vertical layering of the porosity field.

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 CO2 exchange (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.

Distribution forecasts P over 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 functions such 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 functions to 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 rules whose 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 proper scoring 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 proper divergence 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. Plain Language Summary: 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 proper scoring rule. The hydrograph functionals offer a statistically principled means for hydrologic model evaluation, diagnostics and selection. Key Points: 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 [ABSTRACT FROM AUTHOR]

Crosshole ground-penetrating radar (GPR) is a widely used measurement technique to help inspect the structural integrity of man-made underground structures, yet the resulting waveform and travel-time data can be difficult, complex and challenging to interpret. Here, we introduce the elements of a Bayesian inversion method for analyzing crosshole GPR data to guide detection of defects (weakness zones) in underground concrete structures. This framework uses as main building blocks the two-dimensional finite-difference time-domain (FDTD) simulator, the discrete cosine transform (DCT) method, and the DREAM(ZS) algorithm to reconstruct the relative permittivity field of an underground concrete structure from full-waveform GPR inversion. The FDTD simulator solves numerically Maxwell's equations in the time and space domain of the crosshole GPR experiment and simulates iteratively the electromagnetic (EM) waveforms. The DCT algorithm transforms the Cartesian parameterization to the frequency domain and reduces drastically the dimensionality of the parameter space by retaining only the low-frequency DCT-coefficients. Markov chain Monte Carlo (MCMC) simulation with the DREAM(ZS) algorithm is used to estimate the posterior distribution of the DCT-coefficients. The usefulness and applicability of the FDTD-DCT-DREAM(ZS) framework is demonstrated on a synthetic test example involving a unit square concrete structure with a small defect. Our results demonstrate that the proposed method successfully detects and locates defects in concrete structures. The inversion results appear rather insensitive to the noise level of the measured GPR waveforms, and the amount of data used (number of receiving antenna positions), as long as a sufficient number of measurements is available. The more DCT-coefficients that are used to characterize the concrete structure, the more accurate the results, yet the larger the posterior uncertainty of the reconstructed permittivity fields.

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.

Bayesian inference has found widespread application and use in science and engineering to reconcile Earth system models with data, including prediction in space (interpolation), prediction in time (forecasting), assimilation of observations and deterministic/stochastic model output, and inference of the model parameters. Bayes theorem states that the posterior probability, p(H|Y~) of a hypothesis, H is proportional to the product of the prior probability, p( H) of this hypothesis and the likelihood, L(H|Y~) of the same hypothesis given the new observations, Y~, or p(H|Y~)∝p(H)L(H|Y~). In science and engineering, H often constitutes some numerical model, ℱ(x) which summarizes, in algebraic and differential equations, state variables and fluxes, all knowledge of the system of interest, and the unknown parameter values, x are subject to inference using the data Y~. Unfortunately, for complex system models the posterior distribution is often high dimensional and analytically intractable, and sampling methods are required to approximate the target. In this paper I review the basic theory of Markov chain Monte Carlo (MCMC) simulation and introduce a MATLAB toolbox of the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm developed by Vrugt et al. (2008a, 2009a) and used for Bayesian inference in fields ranging from physics, chemistry and engineering, to ecology, hydrology, and geophysics. This MATLAB toolbox provides scientists and engineers with an arsenal of options and utilities to solve posterior sampling problems involving (among others) bimodality, high-dimensionality, summary statistics, bounded parameter spaces, dynamic simulation models, formal/informal likelihood functions (GLUE), diagnostic model evaluation, data assimilation, Bayesian model averaging, distributed computation, and informative/noninformative prior distributions. The DREAM toolbox supports parallel computing and includes tools for convergence analysis of the sampled chain trajectories and post-processing of the results. Seven different case studies illustrate the main capabilities and functionalities of the MATLAB toolbox.

Photosynthetic capacity, determined by light harvesting and carboxylation reactions, is a key plant trait that determines the rate of photosynthesis; however, in Earth System Models (ESMs) at a reference temperature, it is either a fixed value for a given plant functional type or derived from a linear function of leaf nitrogen content. In this study, we conducted a comprehensive analysis that considered correlations of environmental factors with photosynthetic capacity as determined by maximum carboxylation (V(cm)) rate scaled to 25 degrees C (i.e., V(c),25; μmol CO2 x m(-2)x s(-1)) and maximum electron transport rate (J(max)) scaled to 25 degrees C (i.e., J25; μmol electron x m(-2) x s(-1)) at the global scale. Our results showed that the percentage of variation in observed V(c),25 and J25 explained jointly by the environmental factors (i.e., day length, radiation, temperature, and humidity) were 2-2.5 times and 6-9 times of that explained by area-based leaf nitrogen content, respectively. Environmental factors influenced photosynthetic capacity mainly through photosynthetic nitrogen use efficiency, rather than through leaf nitrogen content. The combination of leaf nitrogen content and environmental factors was able to explain -56% and -66% of the variation in V(c),25 and J25 at the global scale, respectively. Our analyses suggest that model projections of plant photosynthetic capacity and hence land-atmosphere exchange under changing climatic conditions could be substantially improved if environmental factors are incorporated into algorithms used to parameterize photosynthetic capacity in ESMs.

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

In the past decades significant progress has been made in the fitting of hydrologic models to data. Most of this work has focused on simple, CPU-efficient, lumped hydrologic models using discharge, water table depth, soil moisture, or tracer data from relatively small river basins. In this paper, we focus on large-scale hydrologic modeling and analyze the effect of parameter and rainfall data uncertainty on simulated discharge dynamics with the global hydrologic model PCR-GLOBWB. We use three rainfall data products; the CFSR reanalysis, the ERA-Interim reanalysis, and a combined ERA-40 reanalysis and CRU dataset. Parameter uncertainty is derived from Latin Hypercube Sampling (LHS) using monthly discharge data from five of the largest river systems in the world. Our results demonstrate that the default parameterization of PCR-GLOBWB, derived from global datasets, can be improved by calibrating the model against monthly discharge observations. Yet, it is difficult to find a single parameterization of PCR-GLOBWB that works well for all of the five river basins considered herein and shows consistent performance during both the calibration and evaluation period. Still there may be possibilities for regionalization based on catchment similarities. Our simulations illustrate that parameter uncertainty constitutes only a minor part of predictive uncertainty. Thus, the apparent dichotomy between simulations of global-scale hydrologic behavior and actual data cannot be resolved by simply increasing the model complexity of PCR-GLOBWB and resolving sub-grid processes. Instead, it would be more productive to improve the characterization of global rainfall amounts at spatial resolutions of 0.5° and smaller.

Catchment urbanization perturbs the water and sediment budgets of streams, degrades stream health and function, and causes a constellation of flow, water quality, and ecological symptoms collectively known as the urban stream syndrome. Low-impact development (LID) technologies address the hydrologic symptoms of the urban stream syndrome by mimicking natural flow paths and restoring a natural water balance. Over annual time scales, the volumes of stormwater that should be infiltrated and harvested can be estimated from a catchment-scale water-balance given local climate conditions and preurban land cover. For all but the wettest regions of the world, a much larger volume of stormwater runoff should be harvested than infiltrated to maintain stream hydrology in a preurban state. Efforts to prevent or reverse hydrologic symptoms associated with the urban stream syndrome will therefore require: (1) selecting the right mix of LID technologies that provide regionally tailored ratios of stormwater harvesting and infiltration; (2) integrating these LID technologies into next-generation drainage systems; (3) maximizing potential cobenefits including water supply augmentation, flood protection, improved water quality, and urban amenities; and (4) long-term hydrologic monitoring to evaluate the efficacy of LID interventions.

We present a Bayesian inversion method for the joint inference of high-dimensional multi-Gaussian hydraulic conductivity fields and associated geostatistical parameters from indirect hydrological data. We combine Gaussian process generation via circulant embedding to decouple the variogram from grid cell specific values, with dimensionality reduction by interpolation to enable Markov chain Monte Carlo (MCMC) simulation. Using the Matérn variogram model, this formulation allows inferring the conductivity values simultaneously with the field smoothness (also called Matérn shape parameter) and other geostatistical parameters such as the mean, sill, integral scales and anisotropy direction(s) and ratio(s). The proposed dimensionality reduction method systematically honors the underlying variogram and is demonstrated to achieve better performance than the Karhunen-Loève expansion. We illustrate our inversion approach using synthetic (error corrupted) data from a tracer experiment in a fairly heterogeneous 10,000-dimensional 2-D conductivity field. A 40-times reduction of the size of the parameter space did not prevent the posterior simulations to appropriately fit the measurement data and the posterior parameter distributions to include the true geostatistical parameter values. Overall, the posterior field realizations covered a wide range of geostatistical models, questioning the common practice of assuming a fixed variogram prior to inference of the hydraulic conductivity values. Our method is shown to be more efficient than sequential Gibbs sampling (SGS) for the considered case study, particularly when implemented on a distributed computing cluster. It is also found to outperform the method of anchored distributions (MAD) for the same computational budget. Key Points: Joint Bayesian inference of Gaussian conductivity fields and their variograms A dimensionality reduction that systematically honors the underlying variogram Distributed multiprocessor implementation is straightforward