In the past decades, Bayesian methods have found widespread application and use in environmental systems modeling. Bayes theorem states that the posterior probability of a hypothesis is proportional to the product of the prior probability of this hypothesis and the likelihood of the hypothesis given the new/incoming observations. In science and engineering, this hypothesis often constitutes some numerical simulation model, which summarizes using algebraic, empirical, and differential equations, state variables and fluxes, all our theoretical and/or practical knowledge of the system of interest, and unknown model parameters which are subject to inference using some data of the observed system response. The Bayesian approach is intimately related to the scientific method and uses an iterative cycle of hypothesis formulation (model), experimentation and data collection, and theory/hypothesis refinement to elucidate the rules that govern the natural world. Unfortunately, model refinement has proven to be very difficult in large part because of the poor diagnostic power of residual based likelihood functions.

In this thesis, I present seven different chapters (publications) that introduce the theory, concepts and implementation of a diagnostic approach to model identification and evaluation. This approach, coined approximate Bayesian computation (ABC), relaxes the need for an explicit likelihood function in favor of one or more summary statistics, which when rooted in the relevant environmental theory have a much stronger and compelling diagnostic power than some average measure of the size of the error residuals. The proposed methodology is statistically coherent, and provides new insights into and guidance on model structural (epistemic) errors, model (hypothesis) refinement, system nonstationarity, and the information content of experimental data. What is more, the proposed diagnostic approach provides a much-needed statistical underpinning of the popular generalized likelihood uncertainty framework (GLUE) of Beven and co-workers, and presents a powerful alternative to subjective regularized inversion methods used in geophysical inversion. Finally, the DREAM(ABC) algorithm, developed to solve the diagnostic model evaluation problem, is orders of magnitude more efficient than commonly used ABC sampling methods, thereby permitting inference of parameter-rich system models.

The topic of infiltration of water into variably saturated soils has received much

attention in the soil physics literature in the past decades. Many different equations

have been proposed to describe quantitatively the infiltration process. These equations

range from simple empirical equations to more advanced deterministic descriptions

of the infiltration process and semi-analytical solutions of Richards’ equation. The

unknown coefficients in these infiltration functions signify hydraulic properties and

must be estimated by curve fitting to measured cumulative infiltration data, Ie(t).

2

From all available infiltration functions, the two-term equation, I(t) = S

√

t + cKst of

Philip [1957] has found most widespread application and use. This popularity has not

only been cultivated by detailed physical and mathematical analysis, the two-term

infiltration equation is also easy to implement and admits a closed-form solution

for the soil sorptivity, S (L T1/2

), and multiple c (−) of the saturated hydraulic

conductivity, Ks (L T−1

). Yet, Philip’s two-term infiltration function has a limited

time validity, tvalid (T), and consequently, measured cumulative infiltration data, Ie(t),

beyond t = tvalid (T) should not be used to estimate S and Ks (among others). The

theoretical treatise in Philip [1957] provides a closed-form solution for the maximum

time validity, t

+

valid, of the two-term infiltration equation. It is not particularly easy to

experimentally corroborate these theoretical findings as this demands prior knowledge

of c, S and Ks

. What is more, the maximum time validity, t

+

valid may not characterize

properly the actual time validity, tvalid. In this paper, we introduce a new method to

determine simultaneously the values of the coefficient c, hydraulic parameters, S and

Ks

, and time validity, tvalid, of Philip’s two-term infiltration equation. Our method is

comprised of two main steps. First, we determine independently the soil sorptivity,

S, and saturated hydraulic conductivity, Ks by fitting the semi-implicit infiltration

equation of Haverkamp [Haverkamp et al., 1994] to measured cumulative infiltration

data. This step uses the DiffeRential Evolution Adaptive Metropolis (DREAM)

algorithm of Vrugt [2016] and returns as byproduct the marginal distribution of the

3

parameter β in Haverkamp’s infiltration equation. In the second step, the maximum

likelihood values of S and Ks are used in Philip’s two-term infiltration equation, and

used to determine the optimal values of c and tvalid via model selection using the

Bayesian information criterion. To benchmark, test and evaluate our approach we

use cumulative infiltration data simulated by HYDRUS-1D [Simunek et al., 2008]

for twelve different USDA soil types with contrasting textures. This allows us to

determine whether our procedure is unbiased as the inferred S and Ks of the synthetic

data are known before hand. Results demonstrate that the estimated values of S and

Ks are in excellent agreement with their ”true” values used to create the artificial

infiltration data. Furthermore, our estimates of c and tvalid are dependent on soil

texture and fall within the ranges stipulated in the literature.

In 2013, the World Meteorological Organization (WMO) urged the global community for coordinated international action against accelerating and potentially devastating climate change. Preliminary data indicated that CO2 levels increased more between 2012 and 2013 than during any other year since 1984, and this was possibly related to reduced uptake by the Earth's biosphere in addition to the steadily increasing emissions from the Earth's surface. In the upcoming decades, it will be critical for scientists and policy makers to not only resolve the problem of carbon emissions by assessing human behavior, but also to understand as thoroughly as possible the underlying coupled processes of the Earth's atmosphere and biosphere in order to adequately measure and estimate the fluxes of carbon, water, and energy that are dictating the climatic trends we observe today. Fortunately, our ability to understand Earth's processes and predict climate change is improving.

This thesis covers a suite of environmental models and numerical methods to disentangle information found both in observed data as well as model simulations. Various methods are applied such as parameter estimation with Markov Chain Monte Carlo (MCMC), state estimation with data assimilation using the Ensemble Kalman Filter (EnKF), and sensitivity analysis of model parameters using the Fourier Amplitude Sensitivity Test (FAST), which all in one way or another offer treatments to predictive uncertainty. Furthermore, applying these methods on more sophisticated and complex models can be impossible sometimes due to their high CPU costs; in this thesis model emulators are built using Polynomial Chaos Expansion (PCE) to reduce the computational burden for various environmental models. Overall, our goal in this dissertation is to present what tools are currently available for making predictions of environmental systems, with emphasis on maintaining accuracy of model simulations when compared to observed data, optimizing the efficiency of computationally heavy models to minimize their run time costs, and obtaining fidelity of model structures to properly represent the underlying hydrologic, biophysical, and biogeochemical processes occurring on our Earth.