2.1 The Bayesian 4D-Var setup

The following equation describes the relationship between a parameter vector, c∈R^m, that we wish to optimize and an observation vector, y∈Rⁿ:

where $f(\cdot ;\cdot )$ is the forward model mapping from parameter space to observation space, μ∈R^p is a control vector of parameters that remain fixed during the optimization, and $\mathbf{R}\in {\mathbf{R}}^{n\times n}$ is the observation covariance matrix. To streamline notation, we denote henceforth the forward model using only the parameter vector as input, i.e., f(c). To regularize the problem and provide uncertainty quantification of the estimated parameter vector, c in Eq. (1) is given a Gaussian prior distribution, yielding the following Bayesian generative model:

where c^b∈R^m is the prior mean and B is the prior covariance matrix. Finding the mode of the posterior c∣y defined by Eqs. (2) and (3) is the objective of 4D-Var data assimilation, a method simultaneously optimizing the parameter vector over all time steps, which has been applied to a wide range of applications in Earth sciences (Deng et al., 2014; Kalmikov and Heimbach, 2014; Liu et al., 2016). The 4D-Var method can be regarded as a least-squares optimization with an ℓ₂ regularizer (ridge regression), or, equivalently, as maximum a posteriori (MAP) estimation in the Bayesian paradigm. This connection means that the 4D-Var optimization is connected to the posterior resulting from the prior and likelihood in Eqs. (2) and (3). From the Bayesian perspective, the 4D-Var cost function F(c) is the negative log-posterior density of the scaling factors given the observations:

where C∈R is a normalizing constant for the posterior distribution and π(c∣y) denotes the posterior density. Thus, finding the MAP estimator, i.e., the c that maximizes the posterior density, is equivalent to finding the vector c that minimizes the 4D-Var cost function. As such, uncertainty quantification can be handled through the covariance matrix of the posterior, which we denote as $\mathbf{\Sigma }\mathbf{=}\text{Cov}\left(\mathit{c}\mathbf{\mid }\mathit{y}\right)$.

To facilitate exposition, we assume the forward model is linear, i.e., f(c)=Ac, though we address strategies to apply the Monte Carlo procedure described in Sect. 2.2 to nonlinear forward models in Sect. 2.2.1. As such, we assume the linear Bayesian model,

on which the primary analysis of this technical note is performed. Not only is the linearity assumption helpful to expose the validity of the Monte Carlo algorithm under consideration, but it is also valid in carbon flux inversion, our primary application of interest. We discuss the necessary adjustments to apply this approach to carbon flux inversion in Sect. 3.2. Assuming linearity, Gaussian prior, and Gaussian likelihood ensures the equivalence between the posterior mode and expectation, and hence optimizing Eq. (4) is equivalent to obtaining the posterior expectation. We emphasize that the matrix A is not explicitly available to us but is instead implicitly defined by a computational model. Thus, even though the linear assumption allows explicit analytical interpretation, we are not able to interact with the matrix in arbitrary ways in practice.

Assuming linearity of the forward model, the posterior mean (mode) and covariance of c are analytically tractable. Hence, $\mathit{c}\mid \mathit{y}\sim \mathcal{N}(\mathit{\alpha },\mathbf{\Sigma })$, where, by Eqs. 4.3 and 4.7 in Rodgers (2000), the posterior mean and covariance equations, we have

Note that α is also the MAP estimator of c. Since the forward model A is only known implicitly via a computer simulator, making direct use of Eqs. (7) and (8) is numerically intractable. Instead, the 4D-Var cost function in Eq. (4) is typically minimized using the L-BFGS-B algorithm (Byrd et al., 1995) with the cost function gradient computed numerically using the adjoint of f (Henze et al., 2007). L-BFGS-B can often find a reasonable approximation of the posterior mean and mode α in a handful of iterations. We now describe the Monte Carlo method, which provides an approach for uncertainty quantification of α despite the analytical intractability of the posterior covariance Σ.

2.2 The Monte Carlo procedure

To execute the Monte Carlo procedure introduced in Chevallier et al. (2007), we generate M ensemble members. For each $k=\mathrm{1},\mathrm{2},\mathrm{\dots },M$, we sample a new prior mean c_k independently and identically distributed (i.i.d.) and new observation y_k i.i.d. as follows:

where c_e is a chosen prior expectation to perturb and y_e∈Rⁿ is a chosen control observation to perturb. Chevallier et al. (2007) operated in a simulation setting where they knew the true parameter, which we call c_clim following their notation. As such, they set c_e=c_clim and y_e=Ac_clim. When using real data, c_clim is not known, and so it is sensible to set this expectation using a physically reasonable value (for example, see Sect. 3.2 for how to set these expectations in more recent carbon flux inversion settings). However, as we show below, the choice of the expectations in Eqs. (9) and (10) does not particularly matter under the linear Gaussian assumption (see Eq. 12), as long as $\text{Cov}({\mathbf{B}}^{-\mathrm{1}}{\mathit{c}}_k,{\mathbf{A}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathit{y}}_k)=\mathrm{0}$. Therefore, each Monte Carlo iteration involves sampling a pair $({\mathit{c}}_k,{\mathit{y}}_k)\in {\mathbf{R}}^m\times {\mathbf{R}}^n$ according to Eqs. (9) and (10).

The MAP estimator from Eq. (8) corresponding to prior mean c_k and observation y_k is analytically tractable for each ensemble member when A is explicitly known:

Since ${\mathit{c}}_{\mathrm{MAP}}^k$ is a linear function of two Gaussian random variables, it is also Gaussian. Therefore, linear combinations of it are also Gaussian, which will be a useful fact in the analysis performed in Sect. 2.3. The covariance matrix of this MAP estimator, henceforth denoted as ${\mathbf{\Sigma }}_{{\mathit{c}}_{\mathrm{MAP}}^k}$, is

where the covariance decomposition in line two follows from $\text{Cov}({\mathbf{B}}^{-\mathrm{1}}{\mathit{c}}_k,{\mathbf{A}}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathit{y}}_k)=\mathrm{0}$ and the final equality follows from ${\mathbf{\Sigma }}^{-\mathrm{1}}=\left(A^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathbf{A}+{\mathbf{B}}^{-\mathrm{1}}\right)$. Thus, the covariance matrix of the Monte Carlo ensemble MAP estimators is equal to the desired posterior covariance Σ. We also note that using the expectation choices from Chevallier et al. (2007) of c_e=c_clim and y_e=Ac_clim means that $\mathbb{E}\left[{\mathit{c}}_{\mathrm{MAP}}^k\right]={\mathit{c}}_{\mathrm{clim}}$.

The above covariance equality is the key fact allowing this method to work, as it allows us to compute an empirical estimate of the posterior covariance by sampling from two unconditional distributions and solving the 4D-Var objective. To the best of our knowledge, proof of this equality has not appeared in previous literature on this method. However, there is a similar method used in the spatial statistics literature to sample from conditional random fields as shown in Chap. 3, Sect. 3.6.2, of Cressie (1993) and Chap. 7, Sect. 7.3.1, of Chiles and Delfiner (2012). As described in Chiles and Delfiner (2012), this sampling method is able to sample from complex conditional random fields if one can sample from unconditional distributions with the appropriate distributions. This requirement more broadly matches the covariance equality requirement of this Monte Carlo procedure and suggests a potential connection between these two sampling procedures.

Since the linear forward model is not explicitly available in most 4D-Var scenarios, each ensemble member MAP estimator ${\mathit{c}}_{\mathrm{MAP}}^k$ must be obtained with an iterative optimization algorithm minimizing Eq. (4) with c_k and y_k as the prior mean and observation vectors. Once these ensemble members are obtained, we could in principle estimate the posterior covariance matrix Σ with the empirical covariance estimator $\widehat{\mathrm{\Sigma }}$ based on the Monte Carlo ensemble as follows:

where $\stackrel{\mathrm{‾}}{\mathit{c}}=\frac{\mathrm{1}}{M}{\sum }_{k=\mathrm{1}}^M{\mathit{c}}_{\mathrm{MAP}}^k$. However, in practice, most data assimilation applications, like carbon flux inversion, are high dimensional, making direct interaction with these covariance matrices difficult. Indeed, accurate estimation of Σ using Eq. (13) would require an enormously large Monte Carlo ensemble and would require storing and working with an m×m matrix, where m∼10⁵ or larger. Fortunately, we often care about the variance of one-dimensional summaries of c, such as the posterior variance of some linear combination of c, as opposed to the full posterior covariance matrix. For instance, we might wish to estimate North American carbon fluxes over a particular month.

Obtaining quantities of the above type is mathematically implemented using a linear functional of the underlying high-dimensional parameter. That is, we wish to characterize the posterior of $\mathit{\phi }\left(\mathit{c}\right)={\mathit{h}}^{\top }\mathit{c}$, where h∈R^m contains weights necessary to aggregate the desired vector components. Hence, building off Eqs. (7) and (8), we obtain the posterior distribution for the functional of interest:

We wish to obtain the posterior variance of this functional. Define ${\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}=\text{Var}\left(\mathit{\phi }\right(\mathit{c})\mid \mathit{y})={\mathit{h}}^{\top }\mathbf{\Sigma }\mathit{h}$. We could inefficiently estimate this using ${\widehat{\mathit{\sigma }}}_{\mathit{\phi }}^{\mathrm{2}}={\mathit{h}}^{\top }\widehat{\mathbf{\Sigma }}\mathit{h}$, but we wish to avoid working directly with the full empirical covariance matrix. The following algebraic steps provide a better alternative:

where ${\mathit{\phi }}_k={\mathit{h}}^{\top }{\mathit{c}}_{\mathrm{MAP}}^k$ and $\stackrel{\mathrm{‾}}{\mathit{\phi }}={\mathit{h}}^{\top }\stackrel{\mathrm{‾}}{\mathit{c}}=\frac{\mathrm{1}}{M}{\sum }_{k=\mathrm{1}}^M{\mathit{\phi }}_k$. The above algebra shows that the posterior variance of the functional can be computed using the functionals of the Monte Carlo samples without having to form the full empirical covariance matrix. See Algorithm 1 for a succinct exposition of the above procedure. Notice that the functional does not need to be specified when creating the Monte Carlo ensemble. As long as the ensemble $\mathit{\{}{\mathit{c}}_{\mathrm{MAP}}^k{\mathit{\}}}_{k=\mathrm{1}}^M$ is stored and made available to the end users, they may evaluate post hoc the uncertainty in any functional that is of interest in their specific use case.

2.3 Quantifying the Monte Carlo uncertainty

Although Sect. 2.2 establishes the equality between the Monte Carlo MAP estimator ensemble member covariance and the posterior covariance (and therefore between the equality of the Monte Carlo ensemble member functional variance and the posterior functional variance), we have not yet established that the empirical covariance matrix (and functional variance) converges in terms of probability to the true posterior covariance matrix (and functional variance). There are consistency results showing the empirical covariance matrix converging in terms of probability to the true covariance matrix (see, for instance, Chap. 6 in Wainwright, 2019). However, since this application is primarily concerned with linear functionals of the form $\mathit{\phi }\left(\mathit{c}\right)={\mathit{h}}^{\top }\mathit{c}$ as described in Sect. 2.2, we can appeal directly to the consistency of the sample variance as shown, for example, in Chap. 5 in Casella and Berger (2002).

Additionally, using the above algorithm, we would like to know either the uncertainty in the variance estimate given the number of Monte Carlo samples or the number of samples required to obtain a particular level of Monte Carlo uncertainty in the variance. In essence, we would like to quantify the uncertainty in our uncertainty. To do so, we take a frequentist approach and construct confidence intervals on ${\widehat{\mathit{\sigma }}}_{\mathit{\phi }}^{\mathrm{2}}$. The confidence intervals can be constructed by recognizing that the ratio of the Monte Carlo functional empirical posterior variance to the true functional posterior variance scaled by (M−1) follows a ${\mathit{\chi }}_{M-\mathrm{1}}^{\mathrm{2}}$ distribution.

Since ${\mathit{c}}_{\mathrm{MAP}}^k$ is Gaussian, the random variables ${\mathit{\phi }}_k={\mathit{h}}^{\top }{\mathit{c}}_{\mathrm{MAP}}^k$ ($k=\mathrm{1},\mathrm{\dots },M)$ are sampled independently and identically from a Gaussian distribution with some mean and variance ${\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}={\mathit{h}}^{\top }\mathbf{\Sigma }\mathit{h}$. By Theorem 5.3.1 of Casella and Berger (2002), we have the following distributional result:

Thus, for $\mathit{\alpha }\in (\mathrm{0},\mathrm{1})$, the distribution in Eq. (20) enables the creation of a 1−α confidence interval for the true posterior variance, ${\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}$, as a function of the empirical posterior variance, ${\widehat{\mathit{\sigma }}}_{\mathit{\phi }}^{\mathrm{2}}$. Using the exact distribution in Eq. (20), we can create either one- or two-sided confidence intervals. Focusing on the two-sided case, we have

where ${\mathit{\chi }}_{M-\mathrm{1},\mathit{\alpha }/\mathrm{2}}^{\mathrm{2}}$ is the $\mathit{\alpha }/\mathrm{2}$ quantile of a chi-squared distribution with M−1 degrees of freedom. Hence, with some algebraic manipulation we arrive at the confidence interval of the posterior variance,

Since in practice we would like to characterize uncertainty in the same units as the flux estimate, we can provide an analogous confidence interval for the posterior standard deviation by taking square roots of all the terms within the probability statement in Eq. (22), giving

Equation (23) facilitates the computation of a $(\mathrm{1}-\mathit{\alpha })\times \mathrm{100}$ % frequentist interval estimator of the Bayesian credible interval for the functional of interest φ. For each endpoint of the true Bayesian credible interval, we find a confidence interval such that the probability that both endpoint confidence intervals simultaneously cover the true credible interval endpoints is 1−α. Let $\mathit{\gamma }\in (\mathrm{0},\mathrm{1})$ and ${\mathit{\phi }}_{\mathrm{MAP}}={\mathit{h}}^{\top }{\mathit{c}}_{\mathrm{MAP}}$ be the functional MAP estimator as described in Sect. 2.1. Because the posterior is Gaussian and φ_MAP is a linear functional, it is a one-dimensional Gaussian random variable. Hence, the Bayesian $(\mathrm{1}-\mathit{\gamma })\times \mathrm{100}$ % credible interval is computed as follows:

where $z_{\mathrm{1}-\mathit{\gamma }/\mathrm{2}}$ is the $\mathrm{1}-\mathit{\gamma }/\mathrm{2}$ quantile of a standard Gaussian distribution. Equation (23) allows us to construct the aforementioned endpoint confidence intervals as follows. For readability, define $L^{\mathrm{2}}:=\frac{M-\mathrm{1}}{{\mathit{\chi }}_{M-\mathrm{1},\mathrm{1}-\mathit{\alpha }/\mathrm{2}}^{\mathrm{2}}}$ and $R^{\mathrm{2}}:=\frac{M-\mathrm{1}}{{\mathit{\chi }}_{M-\mathrm{1},\mathit{\alpha }/\mathrm{2}}^{\mathrm{2}}}$. Thus, we have the following,

More concisely, defining

it follows that

The intervals $\underset{\mathrm{¯}}{I}$ and $\stackrel{\mathrm{‾}}{I}$ quantify uncertainty in uncertainty and provide a rigorous probabilistic characterization of the Monte Carlo procedure's uncertainty.

In practice, the original Bayesian credible interval in Eq. (24) can be modified to account for the Monte Carlo uncertainty. To obtain an upper bound on the Bayesian credible interval, we apply an inflation factor (defined above as R) and thus obtain the interval $\left[{\underset{\mathrm{¯}}{\mathit{\phi }}}_u,{\stackrel{\mathrm{‾}}{\mathit{\phi }}}_u\right]=\left[{\mathit{\phi }}_{\mathrm{MAP}}-z_{\mathrm{1}-\mathit{\gamma }/\mathrm{2}}{\widehat{\mathit{\sigma }}}_{\mathit{\phi }}R,{\mathit{\phi }}_{\mathrm{MAP}}+z_{\mathrm{1}-\mathit{\gamma }/\mathrm{2}}{\widehat{\mathit{\sigma }}}_{\mathit{\phi }}R\right]$ such that $\mathbb{P}\left\{\left[{\underset{\mathrm{¯}}{\mathit{\phi }}}_u,{\stackrel{\mathrm{‾}}{\mathit{\phi }}}_u\right]\supset \left[{\underset{\mathrm{¯}}{\mathit{\phi }}}^{*},{\stackrel{\mathrm{‾}}{\mathit{\phi }}}^{*}\right]\right\}=\mathrm{1}-\mathit{\alpha }/\mathrm{2}$. This probability is $(\mathrm{1}-\mathit{\alpha }/\mathrm{2})$ instead of (1−α), since the probability in Eq. (21) evaluated with only the lower bound yields a probability of $(\mathrm{1}-\mathit{\alpha }/\mathrm{2})$. Following the same steps as above, we obtain lower and upper bounds on the lower and upper endpoints of the credible interval of Eq. (24), respectively, holding with probability of exactly $(\mathrm{1}-\mathit{\alpha }/\mathrm{2})$. Similarly, to obtain a lower bound on the Bayesian credible interval, we apply a deflation factor (defined above as L) and thus obtain the interval $\left[{\underset{\mathrm{¯}}{\mathit{\phi }}}_l,{\stackrel{\mathrm{‾}}{\mathit{\phi }}}_l\right]=\left[{\mathit{\phi }}_{\mathrm{MAP}}-z_{\mathrm{1}-\mathit{\gamma }/\mathrm{2}}{\widehat{\mathit{\sigma }}}_{\mathit{\phi }}L,{\mathit{\phi }}_{\mathrm{MAP}}+z_{\mathrm{1}-\mathit{\gamma }/\mathrm{2}}{\widehat{\mathit{\sigma }}}_{\mathit{\phi }}L\right]$ such that $\mathbb{P}\left\{\left[{\underset{\mathrm{¯}}{\mathit{\phi }}}_l,{\stackrel{\mathrm{‾}}{\mathit{\phi }}}_l\right]\subset \left[{\underset{\mathrm{¯}}{\mathit{\phi }}}^{*},{\stackrel{\mathrm{‾}}{\mathit{\phi }}}^{*}\right]\right\}=\mathrm{1}-\mathit{\alpha }/\mathrm{2}$, holding with equality by the same logic as that used for the inflation factor.

Observing that the aforementioned inflation and deflation factors monotonically approach 1 as the number of Monte Carlo samples M gets large, the effect of the Monte Carlo procedure wanes as the number of samples grows. As shown in Table 2, the deflation factor monotonically approaches 1 from below as M gets large while the inflation factor monotonically approaches 1 from above as M gets large. As is characteristic of DA methods, each Monte Carlo iteration requires a non-trivial amount of computation, which practically restricts the number of Monte Carlo samples that can be obtained. As such, the inflated interval protects against underestimating the uncertainty, while the deflated interval provides a lower bound or “best-case” scenario for the uncertainty in the Bayesian procedure.

Table 2Inflation and deflation factors for Monte Carlo (MC) estimated posterior standard deviation with α=0.05. When M=100, by inflating the MC-estimated posterior standard deviation by a factor of 1.1607 (inflating by 16.07 %), the extra uncertainty resulting from the MC procedure is accounted for with 97.5 % confidence. Similarly, when M=100, deflating the MC-estimated posterior standard deviation by a factor of 0.8785 provides a lower bound on the true underlying Bayesian uncertainty with 97.5 % confidence. When considered simultaneously, the inflation and deflation factors bracket the true uncertainty with 95 % confidence.

Download Print Version | Download XLSX

3.1 Low-dimensional example

We construct a two-dimensional example to provide a numerical demonstration that this MC procedure computes a consistent estimate of the posterior covariance and is numerically close in practice. Define a linear forward model by

where ϵ>0. Let c_true∈R² be the true value of some physical parameter, and suppose we set up Eqs. (5) and (6) as

We use the values in Table 3 to demonstrate the covariance agreement between the analytical equations for the posterior and the Monte Carlo ensemble member in addition to showing their agreement with the empirical covariance computed from the Monte Carlo ensemble members. We also provide an empirical demonstration that the empirical covariance matrix converges to the posterior covariance in the Frobenius norm as the number of ensemble elements M gets large.

Table 3Parameter settings for the low-dimensional example.

Download Print Version | Download XLSX

Using the settings in Table 3, we obtain the following posterior covariance matrix by Eq. (7):

For the analytical covariance of the MAP estimator, we obtain the following matrix using Eq. (12):

Indeed, these matrices are expected to be the same. Using simulated ensemble members, we obtain the following empirical covariance matrix using Eq. (13):

This empirical covariance matrix is qualitatively close to the analytical matrices, even for only M=10². To show how the empirical covariance computed from the Monte Carlo ensemble elements converges to the true posterior covariance, we repeat the above implementation of the Monte Carlo procedure for $M=\mathrm{100},\mathrm{200},\mathrm{\dots },{\mathrm{10}}^{\mathrm{4}}$. For each M value, we compute the Frobenius norm of the difference between the true posterior covariance and the empirical covariance to measure the estimation error, i.e., $\Vert \mathbf{\Sigma }-\widehat{\mathrm{\Sigma }}{\Vert }_F$. Since this estimation error is itself random, for each setting of M, we sample 100 realizations of the Monte Carlo procedure to obtain a mean norm difference. The results of this procedure are presented in Fig. 1 and show that the estimation error exponentially decreases as the ensemble size gets large. By regressing ${\log }_{\mathrm{10}}\Vert \mathbf{\Sigma }-\widehat{\mathbf{\Sigma }}{\Vert }_F$ on log ₁₀M, we estimate the error convergence rate to be CM^−0.49, where C≈10^0.22, as shown by the fitted line in Fig. 1.

Figure 1The Monte Carlo procedure is implemented on the low-dimensional example for a variety of ensemble sizes ($M=\mathrm{100},\mathrm{200},\mathrm{\dots },{\mathrm{10}}^{\mathrm{4}}$) to show that the empirical covariance resulting from the Monte Carlo procedure quickly gets close, in terms of the Frobenius norm, to the posterior covariance as the number of ensemble members gets large. For each ensemble size M, we sample 100 realizations of the Monte Carlo procedure (shown in the gray dots) and find the average over these 100 realizations to see the expected behavior at each ensemble size. The orange regression line results from regressing ${\log }_{\mathrm{10}}\Vert \mathbf{\Sigma }-\widehat{\mathbf{\Sigma }}{\Vert }_F$ on log ₁₀M, where the slope estimates the exponential error convergence rate.

Download

3.2 Carbon flux inversion OSSE

We show an example of this Monte Carlo procedure being used to compute posterior uncertainties for global carbon fluxes along with the adjusted uncertainties from Sect. 2.3. We first detail adjustments to the 4D-Var setup and the Monte Carlo procedure for carbon flux inversion, followed by the particulars of our simulation study and the Monte Carlo procedure results.

3.2.1 Applying 4D-Var and the Monte Carlo procedure to carbon flux inversion

Following along with the mathematical setup of Henze et al. (2007), the prior and posterior distributions are defined in a scaling-factor space, and hence the prior and posterior distributions of the physical quantity of interest are obtained by multiplying each respective scaling factor by a control quantity. In carbon flux estimation, the control quantity is a control flux, typically an ansatz CO₂ flux between the Earth's surface and the atmosphere. Note that if the prior distribution mean in scaling-factor space is unity, then the control flux is also the prior mean in the physical quantity of interest space. Mathematically, let c∈R^m denote a scaling-factor vector; $\stackrel{\mathrm{̃}}{\mathit{y}}\in {\mathbf{R}}^n$ be the observation vector; μ∈R^m be the control physical quantity; and $\mathit{\theta }=\mathit{c}\cdot \mathit{\mu }\in {\mathbf{R}}^m$ be the physical quantity implied by c and μ, where ⋅ denotes element-wise multiplication. We substitute these elements into the original mathematical model stated in Eq. (1). The observation vector $\stackrel{\mathrm{̃}}{\mathit{y}}$ is a sequence of X${}_{{\mathrm{CO}}_{\mathrm{2}}}$ observations produced by a remote sensing satellite, e.g., GOSAT or OCO-2 (O'Dell et al., 2012). The forward model, f, is a composition of an atmospheric transport model mapping scaling factors to atmospheric CO₂ concentrations with a remote sensing observation operator mapping CO₂ concentrations to X${}_{{\mathrm{CO}}_{\mathrm{2}}}$ scalar values. The atmospheric transport model is known to be affine in its scaling factors due to the physics of CO₂ atmospheric transport. The exact mapping from atmospheric CO₂ concentrations to X${}_{{\mathrm{CO}}_{\mathrm{2}}}$ is nonlinear, but in line with Liu et al. (2016), we use an affine approximation involving the known GOSAT averaging kernel. As such, the affine composed function f is of the form $f(\mathit{c};\mathit{\mu })=\mathbf{A}(\mathit{c}\cdot \mathit{\mu })+z$, where $\mathbf{A}\in {\mathbf{R}}^{n\times m}$ is the linear forward model matrix, c⋅μ denotes the component-wise multiplication of μ by the scaling factors c, and z is comprised of the non-biospheric CO₂ contribution to the observations along with the prior mean of the X${}_{{\mathrm{CO}}_{\mathrm{2}}}$ retrieval algorithm. As such, we define $\mathit{y}=\stackrel{\mathrm{̃}}{\mathit{y}}-z$, giving the linear model

$$\begin{array}{}\text{(34)}& \mathit{y}=\mathbf{A}(\mathit{c}\cdot \mathit{\mu })+\mathit{ϵ},\mathit{ϵ}\sim \mathcal{N}(\mathrm{0},\mathbf{R}),\end{array}$$

to which the previous analysis can nearly be applied. These definitions lead to the following Bayesian generative model:

$$\begin{array}{}\text{(35)}& {\displaystyle }\mathit{c}\sim \mathcal{N}\left({\mathit{c}}^b,\mathbf{B}\right),\text{(36)}& {\displaystyle }\mathit{y}\mid \mathit{c}\sim \mathcal{N}\left(\mathbf{A}\right(\mathit{c}\cdot \mathit{\mu }),\mathbf{R}),\end{array}$$

where c^b∈R^m is the scaling factor prior mean and B is the prior covariance matrix.

In this study, the prior covariance B is parameterized with a single real value, B:=b²I_m, where $b\in {\mathbf{R}}_{+}$. This is in line with several published studies (Deng et al., 2014; Liu et al., 2016) and implies that all prior spatiotemporal indices are statistically independent. Similarly, the noise covariance R is assumed to be a diagonal matrix where each diagonal element is simply the variance of the corresponding X${}_{{\mathrm{CO}}_{\mathrm{2}}}$ observation. As such, each diagonal element depends on the uncertainty in its corresponding X${}_{{\mathrm{CO}}_{\mathrm{2}}}$ retrieval, and the observations are assumed statistically independent given the scaling factors.

The uncertainty quantification objective is then to find the posterior uncertainty in a linear functional, defined by vector h∈R^m, of the physical quantity, ${\mathit{h}}^{\top }\mathit{\theta }={\mathit{h}}^{\top }\left(\mathit{c}\cdot \mathit{\mu }\right)$, given the observations. As shown in Algorithm 1, we can obtain an estimate of this uncertainty without computing the entire posterior covariance matrix as long as the covariance of the Monte Carlo ensemble elements is equal to the posterior covariance. We show this equivalence in the following argument.

For ease of notation, we rewrite Eq. (36) using a shorthand notation for the ⋅ operation (see Appendix A), i.e., $\mathbf{A}(\mathit{c}\cdot \mathit{\mu })={\mathbf{A}}_{\mathit{\mu }}\mathit{c}$, so that

$$\begin{array}{}\text{(37)}& \mathit{y}\mid \mathit{c}\sim \mathcal{N}({\mathbf{A}}_{\mathit{\mu }}\mathit{c},\mathbf{R}).\end{array}$$

Hence, applying Eqs. (7) and (8) and accounting for the control flux we have

$$\begin{array}{}\text{(38)}& {\displaystyle }\begin{array}{rl}\mathbf{\Sigma }& ={\left({\displaystyle \frac{\mathrm{1}}{b^{\mathrm{2}}}}{\mathbf{I}}_m+{\mathbf{A}}_{\mathit{\mu }}^{\top }{\mathbf{R}}^{-\mathrm{1}}{\mathbf{A}}_{\mathit{\mu }}\right)}^{-\mathrm{1}}\\ & ={\left(\left({\mathbf{A}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathbf{A}\right)\cdot \mathit{\mu }{\mathit{\mu }}^{\top }+{\displaystyle \frac{\mathrm{1}}{b^{\mathrm{2}}}}{\mathbf{I}}_m\right)}^{-\mathrm{1}},\end{array}\text{(39)}& {\displaystyle }\begin{array}{rl}\mathit{\alpha }& =\mathbf{\Sigma }\left({\mathbf{A}}_{\mathit{\mu }}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathit{y}+{\displaystyle \frac{\mathrm{1}}{b^{\mathrm{2}}}}{\mathit{c}}^b\right)\\ & =\mathbf{\Sigma }\left(\left({\mathbf{A}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathit{y}\right)\cdot \mathit{\mu }+{\displaystyle \frac{\mathrm{1}}{b^{\mathrm{2}}}}{\mathit{c}}^b\right),\end{array}\end{array}$$

where Eq. (38) follows from Corollary 4 and Eq. (39) follows from Lemma 5 in Appendix A. These covariance and expectation equations define the Gaussian posterior for c. The posterior distribution for the physical quantity, $\mathit{\theta }=\mathit{c}\cdot \mathit{\mu }$, is $\mathit{\theta }\mid \mathit{y}\sim \mathcal{N}(\mathit{\delta },\mathrm{\Gamma })$, where

$$\begin{array}{}\text{(40)}& {\displaystyle }\mathit{\delta }=\mathbb{E}[\mathit{\theta }\mid \mathit{y}]=\mathbb{E}[\mathit{c}\cdot \mathit{\mu }\mid \mathit{y}]=\mathbb{E}[\mathit{c}\mid \mathit{y}]\cdot \mathit{\mu }=\mathit{\alpha }\cdot \mathit{\mu },\text{(41)}& {\displaystyle }\begin{array}{rl}\mathrm{\Gamma }& =\text{Cov}[\mathit{\theta }\mid \mathit{y}]=\text{Cov}[\mathit{c}\cdot \mathit{\mu }\mid \mathit{y}]=\text{Cov}[\mathit{c}\mid \mathit{y}]\cdot \mathit{\mu }{\mathit{\mu }}^{\top }\\ & =\mathbf{\Sigma }\cdot \mathit{\mu }{\mathit{\mu }}^{\top }\end{array}\end{array}$$

and we have used Lemmas 1 and 2 from Appendix A.

The argument used in Eq. (12) to establish the equivalence between the posterior covariance and the Monte Carlo ensemble element covariance also holds under the scaling-factor and control flux setup in carbon flux inversion. One simply needs to apply the lemmas in the Appendix to account for the component-wise relationship between the scaling factors and the control flux to establish covariance equality in the scaling-factor space. This equality further extends to the physical quantity space, i.e., carbon flux space. The estimator of the physical quantity corresponding to ${\mathit{c}}_{\mathrm{MAP}}^k$ is

$$\begin{array}{}\text{(42)}& {\mathit{\theta }}_k={\mathit{c}}_{\mathrm{MAP}}^k\cdot \mathit{\mu }.\end{array}$$

Using the result from Lemma 2 in Appendix A, the covariance matrix Cov[θ_k] of this estimator is

$$\begin{array}{}\text{(43)}& \begin{array}{rl}\text{Cov}\left[{\mathit{\theta }}_k\right]& =\text{Cov}[{\mathit{c}}_{\mathrm{MAP}}^k\cdot \mathit{\mu }]=\text{Cov}\left[{\mathit{c}}_{\mathrm{MAP}}^k\right]\cdot \mathit{\mu }{\mathit{\mu }}^{\top }\\ & =\mathbf{\Sigma }\cdot \mathit{\mu }{\mathit{\mu }}^{\top }=\mathrm{\Gamma }.\end{array}\end{array}$$

Hence, the covariance matrix of the Monte Carlo physical quantity ensemble element is equal to the posterior covariance matrix of that physical quantity. Because of this equality, the empirical variance of the functional ensemble can be used to estimate the posterior variance of the functional of interest defined by h, as shown in Step 3 of Algorithm 1. Note that because the validity of Algorithm 1 is based upon the covariance equivalence and not the functional of interest, once a Monte Carlo ensemble has been obtained, one can obtain posterior uncertainties for any collection of desired functionals without having to generate a new Monte Carlo ensemble. This feature is particularly appealing in an application like carbon flux inversion, since scientists are usually interested in estimating posterior uncertainties for many regions (each of which can be encoded by a functional), which would be computationally cumbersome if each functional required a new Monte Carlo ensemble. It also enables the obtaining of uncertainties in new functionals that become the subject of scientific interest only after the generation of the Monte Carlo ensemble.

Figure 2(a) Estimated posterior $(\mathrm{1}-\mathit{\gamma })\times \mathrm{100}\mathit{\%}=\mathrm{95}\mathit{\%}$ credible intervals around the monthly global flux functionals show markedly improved uncertainty over the prior during boreal summer months. The three interval types shown are the unchanged MC-estimated intervals (red), the inflated MC-estimated intervals (gray), and the deflated MC-estimated intervals (orange). As described in Sect. 2.3, for each month, the true upper and lower credible interval endpoints are contained within the inflated and deflated endpoints with probability $\mathrm{1}-\mathit{\alpha }=\mathrm{0.95}$. Note that $\widehat{\mathit{\sigma }}$ is a shorthand notation for the empirical functional standard deviation as defined in Step 3 of Algorithm 1; $\stackrel{\mathrm{‾}}{\mathit{\mu }}$ is the globally averaged control flux to which the prior uncertainty, b, is applied for each month; and R and L are the inflation and deflation factors, respectively, as defined in Table 2. We observe that with even as few as M=60 ensemble members, at the monthly/global scale, the magnitude of the Monte Carlo sampling uncertainty is small in comparison to the posterior uncertainty. (b) The percent reduction in uncertainty from prior to posterior for the monthly global fluxes is most significant during the boreal summer. The light-blue curve shows the percent reduction estimated with the unchanged MC-estimated posterior standard deviation, while the magenta curve shows the percent reduction estimated with the inflated MC-estimated posterior standard deviation. The true reduction is larger than the reduction shown by the latter curve with 97.5 % confidence.

Download