Inverse Modeling of Hydrologic Parameters in CLM4 via Generalized Polynomial Chaos in the Bayesian Framework

Abstract

In this work, generalized polynomial chaos (gPC) expansion for land surface model parameter estimation is evaluated. We perform inverse modeling and compute the posterior distribution of the critical hydrological parameters that are subject to great uncertainty in the Community Land Model (CLM) for a given value of the output LH. The unknown parameters include those that have been identified as the most influential factors on the simulations of surface and subsurface runoff, latent and sensible heat fluxes, and soil moisture in CLM4.0. We set up the inversion problem in the Bayesian framework in two steps: (i) building a surrogate model expressing the input–output mapping, and (ii) performing inverse modeling and computing the posterior distributions of the input parameters using observation data for a given value of the output LH. The development of the surrogate model is carried out with a Bayesian procedure based on the variable selection methods that use gPC expansions. Our approach accounts for bases selection uncertainty and quantifies the importance of the gPC terms, and, hence, all of the input parameters, via the associated posterior probabilities.

1. Introduction

Land surface models (LSMs) are used by scientists to quantitatively simulate the exchange of water and energy fluxes at the Earth–atmosphere interface. LSMs have progressed from simplistic schemes defining simply the surface boundary conditions for general circulation models (GCMs) to complicated models incorporating modules representing biogeochemical, hydrological, and energy cycles at the surface–atmosphere interface over the last few decades (Pitman [1]). The model parameters are frequently related with definite physical meaning and have an influence on the primary model outputs, such as water and energy fluxes. They are built on mathematical formulations of the laws of physics. The assumption for Inter-comparison of Land Surface Parameterization Schemes is that the parameters are quantifiable and transferable between sites with similar physical properties or site conditions (PILPS) (Bastidas et al. [2], Henderson-Sellers et al. [3,4]). However, default assignment of parameter values are actually inappropriate according to (Bastidas et al. [2], Rosero et al. [5]). In the meantime, given the high dimensionality of the parameter space and the complexity of the land surface system, more research is needed to determine which parameters are more uncertain and the potential for using observations to constrain or calibrate the uncertain parameters in order to better capture uncertainty in the resulting land surface states (Hou et al. [6], Huang et al. [7]). The goal of this research is to quantify the uncertainty in a subset of parameters in a community LSM called the Community Land Model (CLM), which is the land component of the Community Earth System Model (CESM) (formerly known as the Community Climate System Model (CCSM) (Collins et al. [8], Gent et al. [9], Lawrence et al. [10]). It is possible that model structural uncertainty in LSMs is attributable to simplified assumptions or representations of actual processes or events. Many of these assumptions hold true only in certain situations. Furthermore, LSMs are prone to ambiguity in terms of input parameter values, owing to the fact that many input characteristics, such as those related to land cover and land use, as well as soil qualities, are difficult to assess at the scales at which they are utilized. A common technique in land surface modeling is to define a set of universally applicable default parameter values. Over the last two decades, the land surface modeling community has worked hard to resolve uncertainty in model parameters, data, and model structure. The purpose of this study is to reduce model parameter uncertainty using generalized polynomial expansion and Bayesian inversion.

In a high-dimensional parameter space, the computational cost of stochastic inversion is significant. Surrogate models can be utilized to solve this problem instead of numerical simulators. Ensemble simulations, which are required to build surrogate models, can be efficiently executed in a task-parallel manner on supercomputing facilities. Surrogate creation, on the other hand, is a lengthy procedure. Surrogates are not utilized very often in climate model or LSM calibration. Using observations of latent heat fluxes, the authors used different surrogates (e.g., polynomials and/or universal kriging) to calibrate hydrological parameters of CLM 4.0 in (Ray et al. [11], Huang et al. [12]).

For the model–data mismatch, two competing models were utilized to construct a composite of measurement error and (a crude estimate of) CLM structural error. In (Gong et al. [13]), the authors employed adaptive surrogate-based optimization to estimate the parameters of the Common Land Model employing six observables at once; 12 independent parameters were (deterministically) calibrated.

Using Bayesian compressive sensing (BCS) and polynomial chaos expansions, (Sargsyan et al. [14]) sought to create surrogates for five variables of interest from CLM4 with prognostic carbon and nitrogen modules turned on (i.e., CLM4-CN) (PCEs). They discovered that the input–output relationship in CLM4-CN could be made up of qualitatively diverse regimes (i.e., living or dead vegetation regimes associated with various locations in the parameter space), necessitating the use of clustering and classification to generate piecewise PCEs.

The applicability of employing gPC for CLM4 hydrological model calibration is investigated in this study. Based on [15], we provide a completely Bayesian technique that combines fully Bayesian statistics, variable selection, and generalized polynomial chaos surrogate models to handle the uncertainty quantification and model inversion problem in CLM4. The procedure produces a cheap mathematical/statistical approximation of the model output (latent heat flux) as a function of a set of model parameters. Bayesian inversion of the model parameters’ given observations is performed by using the produced cheap gPC surrogate model instead of the expensive computer model output in the likelihood function, and then by performing Bayesian parametric inference facilitated by Markov chain Monte Carlo methods. The method allows for dimension reduction and the selection of the important model parameters that significantly influence the output parameter by computing inclusion posterior probabilities.

The layout of the paper is as follows: in Section 2, we describe the study site, input data, and the conducted numerical simulations, and present the parameters of interest that we calibrate; in Section 3, we present the inversion methodology using gPC in the Bayesian framework; in Section 4, we evaluate the inversion results; in Section 5, we draw our conclusions.

2. Dataset and Parameterization

The FLUXNET database (www.fluxdata.org, accessed on 20 March 2022) contains half-hourly observations of ecosystem CO2, heat fluxes, and meteorological data of more than 250 sites worldwide and for a total of 960 site-years. The study site in this study is US-ARM (ARM Southern Great Plains site, Lamont, Oklahoma) (Fischer et al. [16], Torn [17]). It has a vegetation type of croplands, covered with temporary crops followed by harvest and a bare soil period (e.g., single and multiple cropping systems), a humid subtropical climate, and clay-type soil texture.

Observational data used in parameter estimation are observed latent heat fluxes and runoff measurements, which are processed and gap-filled to obtain daily and monthly averaged data. We consider 10 critical hydrological parameters that are likely to have dominant impacts on the simulation of surface and subsurface runoff, latent and sensible fluxes, and soil moisture, as suggested in existing literature (Hou et al. [6], Niu et al. [18,19], Oleson et al. [20,21]). The selected parameters are fmax, Cs, fover, fdrai, qdrai,max (denoted as Qdm hereinafter), Sy, b, Ys, Ks, and $\theta$s. Explanations of the 10 parameters and their prior information are shown in Table 1 in Hou et al. [6]. Prior distributions of the parameters are specified as independent truncated normal distributions with location, scale, and range parameters specified according to the values in Table 1 in Hou et al. [6]; in the cases that prior variance was not available, we considered uniform prior distributions. We recall that such values for the hyper-parameters are acceptable according to existing literature Oleson et al. [20,21]. It is important to highlight that, although truncated normals, these priors are essentially flat due to their large variances. Moreover, 256 samples were generated by using quasi Monte Carlo (QMC) sampling on the basis of the prior pdfs; QMC provides a good dispersion between sample points and can achieve a good uniformity in dimensionality at around 10 Caflisch [22], Wang and Sloan [23]. Numerical simulations corresponding to sampled parameter sets were conducted, which yield the data matrix of inputs (i.e., realizations of the 10 parameters) and outputs (i.e., latent heat fluxes), which enables the development of response surfaces or surrogates that can be used for sensitivity analysis, parameter ranking, and model calibration.

3. Bayesian Methodology

We describe a synergy of Bayesian methods aiming to quantify the importance of input CLM4 model parameters and calibrate these parameters against real measurements, as well as build a surrogate model describing the input–output relation in the CLM4 model, in the Bayesian framework.

3.1. Bayesian Inverse Problem Setup

Bayesian inverse methods allow for the uncertainty quantification of input parameters of a computer model from observations [24]. Bayesian inference is performed through the posterior distribution, which is derived according to the Bayes theorem. We need to specify two components: the likelihood function representing experimental information from the measurements, and the prior distribution representing the researcher’s prior information about the uncertain parameters.

We consider that the output value observed $u^{\mathrm{f}}\in {\mathbb{R}}^{d_u}$ is associated to some unknown input $\xi$ via a forward model (e.g., CLM4) $u(·):{\mathbb{R}}^{d_{\xi }}\to \mathbb{R}$, and possibly contaminated by some additive observational noise ${\epsilon }^{\mathrm{f}}$ (residuals); namely $u^{\mathrm{f}}=u\left(\xi \right)+{\epsilon }^{\mathrm{f}}$. A reasonable assumption is to model ${\epsilon }^{\mathrm{f}}$ as a normally distributed variable with mean zero and variance ${\sigma }^2$. This is justified by central limit theorem arguments, and the fact that, in the noise term, there are several insignificant random measurement errors accumulated from observations and modeling errors due to model parameterization. Thus, the likelihood is:

$$\mathcal{L}\left(u^{\mathrm{f}}\right|\xi )={\mathrm{N}}_{d_u}\left(u^{\mathrm{f}}\right|u\left(\xi \right),\mathrm{diag}({\sigma }_1^2,\dots ,{\sigma }_{d_u}^2)).$$

To account for the uncertainty about the input parameters we assigned prior distribution $\xi \sim {\pi }_{\xi }(·)={\prod }_{i=1}^{d_{\xi }}{\pi }_{{\xi }_i}(·)$, where ${\pi }_{{\xi }_i}\left({\xi }_i\right)$ are considered to be shifted beta distributions $\mathrm{Beta}(a_{{\xi }_i},b_{{\xi }_i},{\xi }_{i,min},{\xi }_{i,max})$, with known ${\xi }_i\in [{\xi }_{i,min},{\xi }_{i,max}]$ and $a_{{\xi }_i}>0$, $b_{{\xi }_i}>0$. Here, ${\xi }_{i,min}$, and ${\xi }_{i,max}$ can be specified because a reasonable range of the model parameter is often a priori known. In addition, $a_{{\xi }_i}$ and $b_{{\xi }_i}$ can be specified by using the method of moments or the maximum entropy method [25] because, often, reasonable values for the mean and variance of the model parameter are a priori known. A computationally convenient choice for the priors of the variance $\left\{{\sigma }_j^2\right\}$ is the inverse gamma prior distribution $\mathrm{IG}(a_{{\sigma }^2},b_{{\sigma }^2})$ with parameters $a_{{\sigma }^2}>0$, and $b_{{\sigma }^2}>0$, as it is a semi-conjugate prior to the likelihood. By considering the likelihood variance ${\sigma }^2$ as random unknown parameters and treating the problem in the Bayesian framework, we let the data decide a proper value for ${\sigma }^2$ through the posterior distribution. The prior hyper-parameters are considered to be fixed values and are pre-defined by the researcher.

The posterior distribution of $\xi ,{\sigma }^2$ given the observations $u^{\mathrm{f}}$ is

$$\begin{array}{cc}\hfill {\pi }_{\xi ,{\sigma }^2}(\xi ,{\sigma }^2|u^{\mathrm{f}})& =\frac{\mathcal{L}\left(u^{\mathrm{f}}\right|\xi ,{\sigma }^2){\pi }_{\xi }\left(\xi \right){\pi }_{{\sigma }^2}\left({\sigma }^2\right)}{\int \mathcal{L}\left(u^{\mathrm{f}}\right|\xi ,{\sigma }^2){\pi }_{\xi }\left(\xi \right){\pi }_{{\sigma }^2}\left({\sigma }^2\right)\mathrm{d}(\xi ,{\sigma }^2)},\hfill \end{array}$$

(1)

by using the Bayes’ theorem, and marginal distribution ${\pi }_{\xi }\left(\xi \right|u^{\mathrm{f}})$ quantifies the uncertainty about the model parameters $\xi$.

The posterior distribution density (1) is usually intractable. If $u(·)$ was known or was cheap to compute, inference on $\xi$ could be obtained using standard Markov chain Monte Carlo (MCMC) methodology [26]. In Algorithm 1, we present a simple MCMC sampler that updates $\xi$ and ${\sigma }^2$ iteratively in two stages. The output sample ${\left\{({\xi }^{\left(t\right)},{\sigma }^{2,\left(t\right)})\right\}}_{t=1}^T$ of Algorithm 1 can be used to perform inference on $\xi$ and ${\sigma }^2$, e.g., expectation of any function $h(·)$ of $\xi$ can be approximated as

$$\mathrm{E}\left(h\left(\xi \right)\right)\approx \frac{1}{T}\sum _{t=1}^{T}h\left({\xi }^{\left(t\right)}\right)$$

for large T.

In the present application of the CLM4 model, $u(·)$ is prohibitively expensive to compute iteratively because the forward model (CLM4) is too expensive to run. This prevents us from directly using this method, and, particularly, Algorithm 1. To overcome this issue, we built a cheap but accurate proxy (called surrogate model) for $u(·)$, and plugged it into (1). The construction of such a surrogate model in the Bayesian framework is discussed in Section 3.2.

Algorithm 1 Updates of a single MCMC swap.

Update $\xi$: Simulate $\xi$ from a Metropolis–Hastings transition probability that targets ${\pi }_{\xi |{\sigma }^2}\left(\xi \right|u^{\mathrm{f}},{\sigma }^2)$, where ${\pi }_{\xi ,{\sigma }^2}\left(\xi \right|u^{\mathrm{f}},{\sigma }^2)\propto \mathrm{N}\left(u^{\mathrm{f}}\right|u\left(\xi \right),{\sigma }^2)\mathrm{N}\left(\xi \right|{\mu }_{\xi },{\sigma }_{\xi }^2)$. Update ${\sigma }^2$: Draw ${\sigma }^2$ from $\mathrm{IG}\left(\frac{1}{2}+a_{{\sigma }^2},\frac{1}{2}{\left|u^{\mathrm{f}}-u\left(\xi \right)\right|}_2^2+b_{{\sigma }^2}\right)$.

3.2. Surrogate Model Specification

We describe a fully Bayesian procedure for building a surrogate model to be used as a cheap but accurate proxy of $u(·)$ in (1) based on gPC expansions and MCMC methods. The highlight is that, apart from evaluating a surrogate model, the procedure is able to quantify the importance of each PC basis via inclusion posterior probabilities, which allows for the selection of the important input parameters: Fmax, Cs, Fover, Fdrai, Qdm, Sy, B, Psis, Ks, thetas.

3.2.1. Generalized Polynomial Chaos Expansion

We consider the output parameter $u\left(\xi \right)$ as a function of the vector of random input variables $\xi \in \mathsf{\Xi }\subseteq {\mathbb{R}}^{d_{\xi }}$ that follow distribution $f(·)$.

The output parameter $u\left(\xi \right)$ can be modeled as an expansion of PC bases $\left\{{\psi }_{\alpha }(·)\right\}$ and PC coefficients $\left\{c_{\alpha }\right\}$:

$$u\left(\xi \right)=\sum _{\alpha \in {\mathbb{N}}_0^{d_{\xi }}}{c}_{\alpha }{\psi }_{\alpha }\left(\xi \right),$$

(2)

for $\xi \sim f(·)$ [27]. We denote multi-indices $\alpha :=({\alpha }_1,\dots ,{\alpha }_{d_{\xi }})$ of size $d_{\xi }$ that are defined on a set of non-negative integers ${\mathbb{N}}_0^{d_{\xi }}:=\{({\alpha }_1,\dots ,{\alpha }_{d_{\xi }}):{\alpha }_j\in \mathbb{N}\cup \left\{0\right\}\}$. The family of polynomial bases $\{{\psi }_{\alpha }(·);\alpha \in {\mathbb{N}}_0^{d_{\xi }}\}$ contains multidimensional orthogonal polynomial bases with respect to the probability measure $f(·)$ of $\xi$. Each multidimensional PC basis ${\psi }_{\alpha }(·)$ results as a tensor product of univariate orthogonal polynomial bases ${\psi }_{{\alpha }_j}(·)$ of degree ${\alpha }_j\in {\mathbb{N}}_0^1$, namely:

$${\psi }_{\alpha }\left(\xi \right)=\prod _{j=1}^{d_{\xi }}{\psi }_{{\alpha }_j}\left({\xi }_j\right),{\alpha }_j\in {\mathbb{N}}_0^1,$$

(3)

where ${\mathrm{E}}_f\left({\psi }_{{\alpha }_j}\left(\xi \right){\psi }_{{\alpha }_{j^{\prime }}}\left(\xi \right)\right)=Z_j{\delta }_0(j-j^{\prime })$, for $j,j^{\prime }=1,\dots ,d_{\xi }$ and $Z_j={\mathrm{E}}_f\left({\psi }_{{\alpha }_j}^2\left(\xi \right)\right)$.

Often, the PC bases $\left\{{\psi }_{\alpha }(·)\right\}$ are chosen so that they are orthogonal with respect to the distribution $f(·)$. A quite general family of polynomials is the Askey family [28], whose members are associated with standard distributions. In this work, we focus on the use of Jacobi polynomial bases [28], which can be defined recursively as:

$$\begin{array}{cc}\hfill {\psi }_0\left(z^{\left(\xi \right)}\right)& =1;\hfill \\ \hfill {\psi }_1\left(z^{\left(\xi \right)}\right)& =\frac{1}{2}[a-b+(a+b-2)z_{\xi }];\hfill \\ \hfill {\psi }_j\left(z^{\left(\xi \right)}\right)& =\frac{(2j+a+b-1)[(2j+a+b)(2j+a+b-2)z^{\left(\xi \right)}+a^2+b^2]}{2j(j+a+b)(2j+a+b-2)}{\psi }_{j-1}\left(z^{\left(\xi \right)}\right)\hfill \\ \hfill & -\frac{2(j+a+b-1)(j+b-1)(2j+a+b)}{2j(j+a+b)(2j+a+b-2)}{\psi }_{j-2}\left(z^{\left(\xi \right)}\right),\hfill \\ \hfill & j=2,\dots ,p_{\xi }-1,\hfill \end{array}$$

where $z^{\left(\xi \right)}$ is a linear transformation $z^{\left(\xi \right)}:[{\xi }_{min},{\xi }_{max}]\to [-1,1]$, and ${\xi }_{min}$, ${\xi }_{max}$ are the minimum and maximum of $\xi$.

In practice, for computational reasons, we truncate (2) by considering a finite sub-set of PC bases based on the total truncated rule. Hence, the expansion form becomes:

$$u_{p_{\xi }}\left(\xi \right):=\sum _{\alpha \in A}{c}_{\alpha }{\psi }_{\alpha }\left(\xi \right),$$

(4)

which uses, for only a finite set of multi-indices, A, such that $A=\{\alpha \in {\mathbb{N}}^{d_{\xi }}:{\sum }_{i=1}^{d_{\xi }}{\alpha }_i\le p_{\xi }\}$ with cardinality $m_{\xi }:=\frac{(p_{\xi }+d_{\xi })!}{p_{\xi }!d_{\xi }!}$. Other truncation rules can be adopted [14,29].

The PC coefficients $\{c_{\alpha };\alpha \in {\mathbb{N}}_0^{d_{\xi }}\}$ are equal to $c_{\alpha }={\mathrm{E}}_f\left(u\left(\xi \right){\psi }_{\alpha }\left(\xi \right)\right)/Z_{\alpha }$, where $Z_{\alpha }={\mathrm{E}}_f\left({\left({\psi }_{\alpha }\left(\xi \right)\right)}^2\right)$, for $a\in {\mathbb{N}}_0^{d_{\xi }}$ [27]; however, they are intractable. Moreover, the number of unknown PC coefficients is of order $d_{\xi }^{p_{\xi }}$ and grows rapidly with the dimension $d_{\xi }$ and PC degree $p_{\xi }$. This causes computational problems such as over-fitting [15,30], especially when the dimensionality of $\xi$ is high. If we consider a smaller PC degree or carelessly omit PC bases, it may lead to a significant increase in bias and hence produce inaccurate surrogate models. Hence, there is a particular interest in keeping only the inputs or bases that significantly affect the output in the gPC. The Bayesian procedure in Section 3.2.2 effectively addresses these matters.

3.2.2. Bayesian Training Procedure

We describe a stochastic and automatic Bayesian procedure that accurately evaluates the PC coefficients and the gPC surrogate model, while it allows for the selection of the the significant PC bases, and hence input model parameters. This procedure is able to trade off efficiently between the bias (caused by omitting bases) and the over-fitting, which are both required. Furthermore, it can select the significant PC bases and estimate the PC coefficients simultaneously, while providing credible intervals.

We assume that there is a training dataset $\mathcal{D}=\left\{(u_j,{\xi }_j)\right\}{}_{j=1}^{n_{\xi }}$ available, where $n_{\xi }$ is the size of the dataset, ${\xi }_j$ denotes the random input value, and $u_j:=u\left({\xi }_j\right)$ denotes the output value corresponding to the j-th input value ${\xi }_j$. Then, we model

$$\begin{array}{cc}\hfill u_j& =u_{p_{\xi }}\left({\xi }_j\right)+{ϵ}_j,\mathrm{for}j=1,\dots ,n_{\xi },\hfill \end{array}$$

(5)

where ${ϵ}_j\in \mathbb{R}$ is the error term, associated with the j-th observation. Equation (5) can be written in matrix form:

$$u=Xc+ϵ,$$

(6)

where $u:={(u_j;j=1:n_{\xi })}^{\top }$, $ϵ:={({ϵ}_j;j=1:n_{\xi })}^{\top }$, $X_a:=({\psi }_a\left({\xi }_j\right);j=1:n_{\xi })$, and $X:=(X_a;a\in A)$ is an $n_{\xi }\times m_{\xi }$ dimensional matrix of basis functions.

The likelihood $\mathcal{L}\left(u\right|c,{\sigma }^2):=\mathcal{L}\left(\left\{u_j\right\}\right|\left\{{\xi }_j\right\},c,{\sigma }^2)$ is:

$$\begin{array}{cc}\hfill \mathcal{L}\left(u\right|c,{\sigma }^2)& =\prod _{j=1}^{n_{\xi }}\mathsf{N}\left(u_j|\psi {\left({\xi }_j\right)}^{\top }c,{\sigma }^2\right);\hfill \\ \hfill & =\mathsf{N}\left(u|X^{\top }c,{\mathbb{I}}_m{\sigma }^2\right),\hfill \end{array}$$

(7)

where $\mathsf{N}(·|\mu ,{\sigma }^2)$ denotes the Gaussian probability density function with mean $\mu$ and variance ${\sigma }^2$. The likelihood function is considered as a proximity measure of the truncated gPC expansion to the observations. We specify a hierarchical prior model $\pi (c,\gamma ,{\sigma }^2,\lambda ,\rho )$:

$$\begin{array}{cc}\hfill c_a|{\gamma }_a,{\sigma }^2,\lambda & \sim {\gamma }_a\mathsf{N}\left(c_a\right|0,{\sigma }^2/\lambda )+(1-{\gamma }_a){\delta }_0\left(c_a\right),a\in A;\hfill \\ \hfill {\gamma }_a& \sim \mathrm{Bernoulli}\left(\rho \right),a\in A;\hfill \\ \hfill {\sigma }^2|a_{\sigma },b_{\sigma }& \sim \mathrm{IG}(a_{\sigma },b_{\sigma });\hfill \\ \hfill \lambda |a_{\lambda },b_{\lambda }& \sim \mathrm{G}(a_{\lambda },b_{\lambda });\hfill \\ \hfill \rho |a_{\rho },b_{\rho }& \sim \mathsf{Beta}(a_{\rho },b_{\rho }),\hfill \end{array}$$

(8)

where $a_{\lambda }$, $b_{\lambda }$, $a_{\rho }$, and $b_{\rho }$ are fixed prior hyper parameters, and predetermined by the researcher. In the Bayesian framework, inference on the uncertain model parameters can be performed based on the posterior distribution

$$\pi (c,\gamma ,{\sigma }^2,\lambda ,\rho |\mathcal{D})=\frac{\mathcal{L}\left(u\right|c,{\sigma }^2)\pi (c,\gamma ,{\sigma }^2,\lambda ,\rho )}{\int \mathcal{L}\left(u\right|c,{\sigma }^2)\pi (c,\gamma ,{\sigma }^2,\lambda ,\rho )\mathrm{d}(c,\gamma ,{\sigma }^2,\lambda ,\rho )}.$$

(9)

Particular interest lies in the computation of the inclusion probabilities $\pi \left({\gamma }_a\right|\mathcal{D})$ that refer to the marginal posterior probability that the a-th basis is important, and the $\pi \left(c_a\right|\mathcal{D},{\gamma }_a)$ that refers to the posterior density of the a-th PC coefficient.

To fit the Bayesian model, we resort to Markov chain Monte Carlo (MCMC) methods because the posterior distribution (9) is intractable and cannot be sampled directly. The conjugate prior model (8) allows for the design of a Gibbs sampler [31,32] whose updates involve sampling from the full conditional distributions of the parameter updated. In Algorithm 2, we represent a pseudo-code of one sweep of the Gibbs samples, along with the associated full conditional distributions. We highlight that the procedure is fully automatic because there is no need to tune the algorithmic parameters involved in Algorithm 2. The notation $c_{-a}$ ( and $X_{-a}$) refers to the vector c (and matrix X) excluding the a-th element (and column). Moreover, we denote $m_{\xi ,\gamma }={\sum }_{a\in A}{\gamma }_a$.

Algorithm 2 Updates of the Gibbs swap.

Update $\left\{({\gamma }_a,c_a)\right\}$: For $a\in A$, Compute ${\mu }_a$, and $s_a^2$, where $$\begin{array}{cc}\hfill {\mu }_a& ={\left(X_a^{\top }{X}_a+\lambda \right)}^{-1}{X}_a^{\top }\left(u-X_{-a}{c}_{-a}\right),\hfill \\ \hfill s_a^2& ={\sigma }^2{(X_a^{\top }{X}_a+\lambda )}^{-1}\hfill \end{array}$$ Update ${\gamma }_a$: draw ${\gamma }_a$ from $\mathrm{Bernoulli}\left(P_a^{\left(\gamma \right)}\right)$, where $$P_a^{\left(\gamma \right)}={\left[1+\frac{1-\rho }{\rho }\sqrt{\frac{2\pi {\sigma }^2}{\lambda }}\mathsf{N}\left(0\right|{\mu }_a,s_a^2)\right]}^{-1},$$ Update $c_a$: draw $c_a$ from $\pi \left(c_a\right|u,X,\gamma ,c_{-a},\rho ,{\sigma }^2,\lambda )$, where $$\pi \left(c_a\right|u,X,\gamma ,c_{-a},\rho ,{\sigma }^2,\lambda )=\left\{\begin{array}{cc}{\delta }_0\left(c_a\right)\hfill & ,\mathrm{if}{\gamma }_a=0\hfill \\ \mathsf{N}\left(c_a\right|{\mu }_a,s_a^2)\hfill & ,\mathrm{if}{\gamma }_a=1\hfill \end{array}\right..$$ Update ${\sigma }^2$:   Sample ${\sigma }^2$ from $\mathrm{IG}\left(\frac{n_{\xi }}{2}+\frac{m_{\xi ,\gamma }}{2}+a_{\sigma },b_{\sigma }+\frac{n_{\xi }}{2}{\left|u-Xc\right|}_2+\frac{1}{2}\lambda {\left|c\right|}_2^2\right)$. Update $\lambda$: , Sample $\lambda$ from $\mathsf{G}\left(\frac{m_{\xi ,\gamma }}{2}+a_{\lambda },\frac{1}{2{\sigma }^2}{\left|c_a\right|}_2^2+b_{\lambda }\right)$. Update $\rho$:   Sample $\rho$ from $\mathrm{Beta}(m_{\xi ,\gamma }+a_{\rho },m_{\xi }-m_{\xi ,\gamma }+b_{\rho })$.

In order to evaluate the surrogate model, as well as quantify the importance of the input model parameters, we consider two Bayesian procedures. One is the Bayesian model averaging [33] and the other is the median probability model.

Bayesian Model Averaging

Bayesian model averaging is most suitable when we are interested in the predictive ability of the surrogate model.

Let ${\left\{\left({\gamma }^{\left(t\right)},c^{\left(t\right)},{\sigma }^{2,\left(t\right)},{\lambda }^{\left(t\right)},{\rho }^{\left(t\right)}\right)\right\}}_{t=1}^T$ be a Gibbs sample from Algorithm 2. Estimates and associated standard errors of $\left\{c_a\right\}$ can be computed as ${\widehat{c}}_a=\frac{1}{T}{\sum }_{t=1}^T{c}_a^{\left(t\right)}$ and $\mathrm{s}.\mathrm{e}.\left({\widehat{c}}_a\right)=s_a^c\sqrt{\frac{{\varrho }_a^c}{T}}$, where $s_a^c$ is the sample standard deviation and ${\varrho }_a^c$ is the sample integrated autocorrelation time of $\left\{c_a^{\left(t\right)};t=1:T\right\}$ for $a\in A$ [26,33]. Estimates for $\mu$, v, and $u\left(\xi \right)$ can be computed by Monte Carlo integration using the ergodic average of quantities in Gibbs sample; for instance:

$$\begin{array}{cc}\hfill \widehat{u}\left(\xi \right)& =\frac{1}{T}\sum _{t=1}^T\left(\sum _{a\in A}{c}_a^{\left(t\right)}{\psi }_a\left(\xi \right)\right)=\sum _{a\in A}{\widehat{c}}_a{\psi }_a\left(\xi \right);\hfill \\ \hfill \widehat{P_a}& =\frac{1}{T}\sum _{t=1}^T𝟙\left({\gamma }_a\right);\hfill \\ \hfill \widehat{\mu }& =\frac{1}{T}\sum _{t=1}^{T}c{}_0^{\left(t\right)}={\widehat{c}}_0,\hfill & \hfill \widehat{v}& =\frac{1}{T}\sum _{t=1}^T\sum _{\alpha \in A-\left\{0\right\}}{\left(c_{\alpha }^{\left(t\right)}\right)}^2{Z}_{\alpha };\hfill \\ \hfill {\widehat{\sigma }}^2& =\frac{1}{T}\sum _{t=1}^T{\sigma }^{2,\left(t\right)},\hfill & \hfill \widehat{\lambda }& =\frac{1}{T}\sum _{t=1}^T{\lambda }^{\left(t\right)};\hfill \\ \hfill {\widehat{c}}_a& =\frac{1}{T}\sum _{t=1}^T{c}_a^{\left(t\right)},\hfill & \hfill \widehat{\rho }& =\frac{1}{T}\sum _{t=1}^T{\rho }^{\left(t\right)}.\hfill \end{array}$$

(10)

Median Probability Model Based Evaluation

A parsimonious (or sparse) surrogate model involving only significant basis functions and important input model parameters can be obtained by examining the estimated inclusion probabilities $\pi \left({\gamma }_a\right|\mathcal{D})$. A suitable probabilistic basis selection mechanism is the median probability model (MPM) [34].

Let ${\left\{\left({\gamma }^{\left(t\right)},c^{\left(t\right)},{\sigma }^{2,\left(t\right)},{\lambda }^{\left(t\right)},{\rho }^{\left(t\right)}\right)\right\}}_{t=1}^T$ be a Gibbs sample from Algorithm 2. The posterior distribution that the PC basis is significant can be estimated as ${\widehat{P}}_a=\widehat{\pi }\left({\gamma }_a\right|\mathcal{D})=\frac{1}{T}{\sum }_{t=1}^T{\gamma }_a^{\left(t\right)}$. We call $\left\{P_a\right\}$ the marginal inclusion posterior probabilities. The inclusion parameters are estimated as ${\widehat{\gamma }}_a^{\left(\mathrm{MPM}\right)}=𝟙({\widehat{P}}_a\in (0.5.1))$, for $a\in A$. MPM consists of PC coefficients $\left\{c_a\right\}$, and PC bases $\left\{{\psi }_a(·)\right\}$, whose marginal inclusion probabilities are such that ${\widehat{P}}_a>0.5$ (and hence ${\widehat{\gamma }}_a^{\left(\mathrm{MPM}\right)}=1$). Hereafter, we will refer to them as significant PC coefficients and PC bases.

After the selection of the significant PC bases, inference about the uncertain parameters can be made by re-running Algorithm 2 for fixed $\gamma$ equal to ${\widehat{\gamma }}^{\left(\mathrm{MPM}\right)}$ in order to obtain a a larger Gibbs sample. The new Gibbs sample can be used to perform inference and estimation. For instance, estimates for $\mu$, v, and $u\left(\xi \right)$ can be computed by Monte Carlo integration and using the equations of the estimators in (10). Note that a number of the coefficients $\left\{c^{\prime }{}_a^{\left(t\right)}\right\}$ will be constantly zero for all $t=1,\dots ,T$. The reason is because, unlike in the BMA approach, here, we consider only a single subset of PC bases, and hence the inclusion parameter is a fixed parameter equal to ${\gamma }^{\left(t\right)}={\gamma }^{\left(\mathrm{MPM}\right)}$, for $t=1,\dots ,T$.

The MPM approach allows the selection of the important input model parameters that significantly affect the output. If an input parameter ${\xi }_j$ is not represented by any significant PC basis in the gPC expansion, it would be reasonable to consider that input parameter ${\xi }_j$ does not significantly affect the output model and, hence, is omitted from the analysis.

4. Analysis of the US-ARM Data-Set

Here, we consider the US-ARM data set. The main interest lies in computing the posterior distributions of the 10 (random) input parameters of CLM4 $\xi =\left(\mathrm{Fmax},\mathrm{Cs},\mathrm{Fover},\mathrm{Fdrai},\mathrm{Qdm},\mathrm{Sy},\mathrm{B},\mathrm{Psis},\mathrm{Ks},\theta \mathrm{s}\right)$, given an observed latent heat flux (LH) measurement $u^{\left(\mathrm{f}\right)}$ as shown in

Table 1. Observed value $u^{\left(\mathrm{f}\right)}$ for the US-ARM data-set.

Months	Measurement
January-April	$15.542$	$22.017$	$41.365$	$59.095$
May-July	$58.377$	$58.813$	$45.107$	$41.362$
August-December	$31.250$	$28.645$	$17.635$	$12.778$

.

We apply the above methodology, which involves two stages: (i) building a surrogate model to replace the accurate but expensive forward model CLM4 according to the methodology in Section 3.2, and (ii) conducting inversion of the 10 inputs according to the theory in Section 3.1.

Surrogate Model Building Step

For each month, we built a surrogate model that maps the input of CLM4 $\xi$ to the output LH u. For this type of dataset, Hou et al. [6] and Sun et al. [35] observed that the dependency between the output parameters LH corresponding to different months is weak and, hence, can be neglected. Therefore, here, we built surrogate models for each month independently. An advantage of assuming independence is that it leads to a simpler parameterization for the statistical model, which is easier to treat and interpret.

To build the gPC expansion, we followed the procedure in Section 3.2. For the design of the gPC expansion, we considered PC bases from the Jacobi polynomial family, and the total truncated rule with polynomial degree 4. The bases in the gPC expansion had the following arrangement: first, the main effects as Fmax, Cs, Fover, Fdrai, Qdm, Sy, B, Psis, Ks, $\theta$, and, then, the interactions by increasing order, while interactions of the same order were arranged so that their interacting bases correspond to main effects in the aforesaid order. This rule identifies the order of the columns of the design matrix X, as well as that of the elements of the vector C in (6). The parameters of the Jacobi PC bases were set according to the prior information of the input of CLM4 in ([6], Table 1) by matching the moments of the corresponding shifted beta distribution, to which, they are orthogonal. We considered the prior model (8) with hyper-parameters $a_{\lambda }={10}^{-3}$, $b_{\lambda }={10}^{-3}$, $a_{\sigma }={10}^{-3}$, $b_{\sigma }={10}^{-3}$, $a_{\rho }=1$, and $b_{\rho }=1$. This choice of hyper-parameters leads to weakly informative priors. This is a reasonable choice because there is lack of prior information about the parameterization of the surrogate model. For training the gPC expansion, we used the training data set US-ARM. We ran Algorithm 2 for $2·{10}^5$ iterations, where the first ${10}^5$ were discarded as burn-in.

In Figure 1 and Figure 2, we present the marginal posterior probabilities $\{Pr\left({\gamma }_j\right|\mathcal{D})\}$ computed by the ergodic average of the occurrences of the corresponding PC bases in the Gibbs sample. We observe that, during the period May-August, the marginal inclusion probabilities are higher than those of the rest months. This indicates that the input parameters of CLM4 may have a larger impact on the output LH during those months. From the hydrology point of view, this is expected because LH is higher on average and has a larger variability during these months, and the effects of hydrological parameters are expected to be more pronounced in the summer months.

We can infer that the input parameters Fmax, Cs, Fdrai, Qdm, and B are significant according to the MPM rule. This is because these input parameters are represented by significant PC basis functions; namely, the corresponding marginal inclusion probabilities in Figure 1 and Figure 2 are greater than $0.5$. The order of the indices in Figure 1 and Figure 2 is the same as that mentioned above in the surrogate model building step section; namely, first, the main effects are followed by their interactions in the order of Fmax, Cs, Fover, Fdrai, Qdm, Sy, B, Psis, Ks, $\theta$. The first 11 indexes correspond to the main effects, indices 12–66 correspond to the second-order interactions, indices 67–286 correspond to third-order interactions, and the rest of the indices correspond to fourth-order interactions. From the hydrology perspective, this is reasonable, because these parameters are major factors controlling the drainage and runoff generation, which, in turn, impact heat fluxes. The results are also consistent with the previous work in (Hou et al. [6]).

In Figure 3 and Figure 4, we use box-plots to represent the posterior density estimates of the PC coefficients generated by the Gibbs sampler. In these figures, we consider only the coefficients corresponding to bases characterized as significant according to the MPM rule; namely, those with $Pr\left({\gamma }_j\right|\mathcal{D})<0.5$. We observe that significant coefficients do not include zero in the range where the main posterior mass is concentrated. This shows that the method is consistent. Moreover, we observe that the significant PC coefficients that correspond to the period May-August have, in general, larger absolute values. This indicates that the variance of the output LH during those months is larger compared to that of the rest months, which is as expected. Finally, we can possibly infer that different input parameters $\xi$ may have a different effect of LH at different times, since we can observe different significant coefficients, or bases, in different months.

The gPC expansion of LH as a function of the input parameters $\xi$ can be evaluated according to the estimator in (10), which is the ergodic average of the Gibbs sample.

4.1. Inversion Step

We calibrate the 10 CLM parameters $\xi =\left(\mathrm{Fmax},\mathrm{Cs},\mathrm{Fover},\mathrm{Fdrai},\mathrm{Qdm},\mathrm{Sy},\mathrm{B},\mathrm{Psis},\mathrm{Ks},\theta \right)$ against the measurement of the parameter LH $u^{\left(\mathrm{f}\right)}$ in Table 1, as in Section 3.1. We considered beta priors on the CLM parameters whose hyper-parameters are specified by using the method of moments and based on the prior information in Hou et al. [6].

Calibration was performed by running the MCMC sampler (Algorithm 1) and evaluating the posterior distributions according to the procedure in Section 3.1. Even though, in the previous step, we detected that the input model parameters Fmax, Cs, Fdrai, Qdm, and B are the significant ones, we also considered them for calibration to obtain information about them as well. In order to make the MCMC sampler tractable, we replaced the forward model CLM4 $u(·)$, in Algorithm 1, with the estimated gPC expansion that serves as a surrogate model. We ran the MCMC sampler (Algorithm 1) for $2·{10}^4$ iterations and discarded the first ${10}^4$ as burn-in.

In Figure 5 and Figure 6, we present the estimated posterior densities of the input parameters of CLM4, as generated by the MCMC sampler. The blue bars correspond to the histogram estimate, whereas the red line corresponds to the kernel density estimate. These posterior distributions allow us to find a reasonable range of input values that correspond to the given value of output $u^{\mathrm{f}}$. The associated box-plots of the marginal posteriors of the model parameters at each individual figure indicate the range of the main posterior density. We can see that we have successfully managed to shorten the ranges of the possible values for most of the input model parameters. For instance, we observe that the main density on the marginal posterior distribution of the density of Qdm is around the area $[-4,-2.5]$ in ${log}_{10}$ scale.

4.2. Model Validation

We validated the effectiveness of the Bayesian inversion procedure. We evaluated the predictive distribution of the output parameter LH by using the derived gPC surrogate model and the MCMC sample of the input parameters generated by Algorithm 1. In Figure 7 and Figure 8, we present the resultant predictive distributions of the output LH for the 12 months. The blue bars are the histogram estimate, the red line is the kernel density estimate of the predictive distribution, and the green arrow represents the observed output value of LH. The plots show that, for each month, the observed output value $u^{\mathrm{f}}$ for the LH lies below the modes of each of the marginal predictive distributions. This implies that the proposed methodology is valid and that the surrogate model derived from the method is able to produce accurate predictions.

5. Conclusions

We focused on evaluating the uncertainties associated with the hydrologic parameters in CLM4 in the Bayesian framework. We presented a Bayesian methodology for the uncertainty quantification (UQ) framework that couples a generalized polynomial chaos model with the Bayesian variable selection methods. We presented a fully Bayesian methodology that involves two steps: building a surrogate step, which is the construction of a surrogate model to express the input–output mapping; and the inverse modeling step, which is the evaluation of the posterior distribution of the input parameters for a given value of the output LH.

For the construction of the surrogate model, we propose a Bayesian procedure, based on variable selection methods, that uses gPC expansions and accounts for bases selection uncertainty. The advantage of this approach is that it can quantify the significance of the gPC terms, and, hence, the importance of the input parameters, in a probabilistic manner. The input posterior distributions were evaluated according to Bayesian inverse modeling. Our numerical experiments suggested that the method presented is suitable for performing the inverse modeling of hydrologic parameters in CLM4, and able to effectively describe the uncertainty related to these parameters.

Our future work involves the comparison of the proposed method against other methods based on neural networks, and generalized linear models, on which, we are currently working.