2.1 Analytical solution to the inverse problem

The optimal estimate $\widehat{\mathit{x}}$ of a state vector x, given a prior estimate x_A, observation vector y, and normal error statistics given by prior and observational error covariance matrices S_A and S_O, respectively, is obtained by the minimization of the Bayesian scalar cost function 𝒥(x) given by the following (Brasseur and Jacob, 2017):

Here F(x) represents the forward model that simulates the observations y given x. In our application, the forward model is a CTM. The observational error covariance matrix S_O includes errors from both the measurement and the forward model, which collectively form the observing system. If the forward model is linear so that $\mathbf{F}\left(\mathit{x}\right)=\mathbf{K}\mathit{x}+\mathit{c}$, where $\mathbf{K}=\partial \mathit{y}/\partial \mathit{x}$ is the Jacobian matrix calculated by finite difference (see the Introduction) and c is a constant, then an analytical solution to the cost function minimum exists that yields both the posterior estimate $\widehat{\mathit{x}}$ and its error covariance matrix $\widehat{\mathbf{S}}$ as follows​​​​​​​:

Comparison of $\widehat{\mathbf{S}}$ and S_A defines the information content of the observing system, quantified by the averaging kernel matrix $\mathbf{A}=\partial \widehat{\mathit{x}}/\partial \mathit{x}$ that represents the sensitivity of the posterior emissions estimate $\widehat{\mathit{x}}$ to the true state x. A can be calculated as $\mathbf{A}=\mathbf{I}-\widehat{\mathbf{S}}{\mathbf{S}}_{\mathrm{A}}^{-\mathrm{1}}$ or equivalently as follows:

Equation (4) expresses the dependence of the averaging kernel matrix on the forward model and both error covariance matrices. The diagonal elements of A are commonly referred to as the averaging kernel sensitivities. They are highest in highly observed locations with uncertain, high emissions and lowest in poorly observed areas or in regions known to have low emissions. The sum of the sensitivities, or the trace of A, measures the number of pieces of information that can be independently quantified by the observing system, known as the degrees of freedom for signal or DOFS (Rodgers, 2000).

2.2 Optimal reductions in dimension and rank of inverse systems

We first consider the problem of optimally reducing the dimension and rank of an inverse system, as described in Sect. 2.1, with a known Jacobian matrix $\mathbf{K}\in {\mathbb{R}}^{m\times n}$. Figure 1 illustrates dimension and rank reductions for an emission grid over North America. The top left panel represents the original n-dimensional state space, i.e., the native-resolution grid. A linear transformation $\mathbf{\Gamma }\in {\mathbb{R}}^{k\times n}$ reduces the dimension of the state space from n to k. This transformation may reduce dimension discretely, as in the case of grid cell aggregation (top right panel), or non-discretely, in which case the k state vector components are themselves n-dimensional vectors (bottom right panel). A second linear transformation ${\mathbf{\Gamma }}^{\ast }\in {\mathbb{R}}^{n\times k}$ restores the dimension of the state space from k back to the original n. The resulting space, depicted in the bottom left, is a low-rank approximation of the original state space. The matrix $\mathbf{\Pi }={\mathbf{\Gamma }}^{\ast }\mathbf{\Gamma }$ transforms the original state space to the low-rank subspace. The inverse problem can be solved in any of these four spaces, although the eigenvector corrections generated in the non-discrete reduced-dimension space (bottom right panel) would be difficult to interpret.

Figure 1Dimension and rank reductions of a gridded emissions field. The linear transformation matrix Γ reduces the dimension of the original state space (a), either discretely by aggregating grid cells to generate a multiscale grid (b) or non-discretely by projecting along the patterns given by the rows of Γ (c, with positive values in red and negative in blue). The reverse transformation Γ^∗ restores the dimension but not the rank, producing a low-rank subspace of the original state space (d). The projection $\mathbf{\Pi }={\mathbf{\Gamma }}^{\ast }\mathbf{\Gamma }$ reduces the rank but not the dimension.

We wish to define matrices Γ and Γ^∗ that minimize the information loss associated with reducing the dimension or rank of the state vector. Bousserez and Henze (2018) show that the projection Π that maximizes the probability of restoring the original full dimension state vector x, given the reduced dimension state vector Γx, is given by $\mathbf{\Pi }={\mathbf{S}}_{\mathrm{A}}^{\mathrm{1}/\mathrm{2}}{\mathbf{UU}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{A}}^{-\mathrm{1}/\mathrm{2}}$, where $\mathbf{U}={\mathbf{S}}_{\mathrm{A}}^{\mathrm{1}/\mathrm{2}}\mathbf{\Gamma }{\left(\mathbf{\Gamma }{\mathbf{S}}_{\mathrm{A}}{\mathbf{\Gamma }}^{\mathrm{T}}\right)}^{-\mathrm{1}/\mathrm{2}}$. For a projection of this form, they show that information loss is minimized by maximizing ${\mathrm{DOFS}}_{\mathrm{\Pi }}=\mathrm{Tr}\left({\mathbf{A}}_{\mathrm{\Pi }}\right)=\mathrm{Tr}\left({\mathbf{U}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{A}}^{-\mathrm{1}/\mathrm{2}}{\mathbf{AS}}_{\mathrm{A}}^{\mathrm{1}/\mathrm{2}}\mathbf{U}\right)$, where A_Π and A are the reduced-rank and native-resolution averaging kernel matrices, respectively. They define as follows:

where the columns of W are the eigenvectors of Q, and Σ is a diagonal matrix of the corresponding eigenvalues ranked in descending order. The quantity Tr(A_Π) is maximized for a rank-k subspace when U=W_k, where W_k is the matrix of the first k columns of W. The corresponding optimal projection is then as follows:

This projection applies a dimension-reducing transformation Γ followed by a dimension-restoring transformation Γ^∗ as follows:

The columns of Γ^∗ give an eigenvector basis for the averaging kernel matrix, and the diagonal of Σ gives its eigenvalues, together defining the dominant patterns of information content. The fraction of information content explained by the first i columns of Γ^∗ is the sum of the i largest eigenvalues divided by the total DOFS (Bousserez and Henze, 2018). The eigenvalues can also be related to other measures of information content, including the Shannon and relative entropy differences (Rodgers, 2000; Xu, 2007). We will refer to the ordered list of the eigenvalues as the information content spectrum. On the basis of this spectrum, we can select k so that most of the information content is explained by the first k eigenvectors. Alternatively, we can select k so that all eigenvectors have a sufficiently large signal-to-noise ratio. The signal-to-noise ratio SNR of the ith eigenvector is given by the ith singular value of the pre-whitened Jacobian matrix $\stackrel{\mathrm{̃}}{\mathbf{K}}={\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}/\mathrm{2}}{\mathbf{KS}}_{\mathrm{A}}^{\mathrm{1}/\mathrm{2}}$ and is calculated as follows:

where σ_i is the ith ordered eigenvalue of Q (Rodgers, 2000).

2.3 Approximating the Jacobian matrix

Section 2.2 described optimal reductions in the dimension and rank of a state vector assuming knowledge of the native-resolution Jacobian matrix K. However, the n+1 forward model simulations needed to construct K may be prohibitively expensive. Here we present a two-step approach to construct a reduced-dimension or reduced-rank Jacobian matrix at much lower computational cost. We start from a low-cost, native-resolution estimate K⁽⁰⁾ (see below) and calculate the corresponding averaging kernel matrix A⁽⁰⁾. In the reduced-dimension method, we use A⁽⁰⁾ to construct a multiscale grid that maintains resolution in the areas of highest information content (top right panel of Fig. 1). We generate the updated, reduced-dimension Jacobian matrix ${\mathbf{K}}_{\mathrm{RD}}^{\left(\mathrm{1}\right)}\in {\mathbb{R}}^{m\times k}$ on the resulting grid using the forward model. In the reduced-rank method, we construct ${\mathbf{K}}_{\mathrm{\Pi }}^{\left(\mathrm{1}\right)}\in {\mathbb{R}}^{m\times n}$ on the basis of the k dominant eigenvectors of A⁽⁰⁾ by perturbing those patterns in the forward model, generating an approximation of the Jacobian matrix in a reduced-rank state space (bottom left panel of Fig. 1). In both methods, the updated Jacobian matrix improves the estimate of the averaging kernel matrix and its eigenvectors by incorporating information content from the forward model. We use either ${\mathbf{K}}_{\mathrm{RD}}^{\left(\mathrm{1}\right)}$ or ${\mathbf{K}}_{\mathrm{\Pi }}^{\left(\mathrm{1}\right)}$ to conduct a second update and construct the final Jacobian matrix.

In our demonstration case, we generate an initial estimate of the native-resolution Jacobian matrix K⁽⁰⁾ at no cost by assuming that a local perturbation of methane emissions Δx (kilograms per square meter per second; hereafter kg m⁻² s⁻¹) produces local dry column mixing ratio enhancements Δy (moles per mole; hereafter mol mol⁻¹), as determined by a simple column mass balance dependent on local wind speed and parameterized turbulent diffusion. We construct K⁽⁰⁾ column-wise by assuming that observation i responds to emissions in grid cell j as follows:

so that the elements $k_{ij}^{\left(\mathrm{0}\right)}$ of K⁽⁰⁾ are given by

where ${\mathit{\alpha }}_{ij}\in [\mathrm{0},\mathrm{1}]$ is a dimensionless coefficient providing a crude parameterization of turbulent diffusion, M_air and $M_{{\mathrm{CH}}_{\mathrm{4}}}$ are the molecular weights of dry air and methane, respectively, L is a ventilation length scale taken as the square root of the grid cell area, g is gravitational acceleration, U is the local wind speed taken here as 5 km h⁻¹, and p is the surface pressure. We assume α_ij=0.4 for observations in grid cell j and distribute the remaining mass over the three concentric rings surrounding that cell with α_ij=0.3/8, 0.2/16, and 0.1/24 from the inner to outer ring. Including a representation of turbulent diffusion increases the spatial coverage of the dominant patterns of information content. The exact form of the parameterization (e.g., the number of rings used or the values of α_ij) is unimportant.

The reduced-dimension and reduced-rank methods rely on characterizing the dominant patterns of information content of the observing system using the initial estimate of the averaging kernel matrix A⁽⁰⁾ corresponding to K⁽⁰⁾. A⁽⁰⁾ can provide a good approximation of A, even if the initial estimate of the Jacobian matrix K⁽⁰⁾ is crude, because the averaging kernel matrix depends strongly on the specified prior and observational error covariance matrices S_A and S_O (Eq. 4) and because, by assuming that observed concentrations are most sensitive to local emissions, K⁽⁰⁾ generates the highest information content where the observations are densest. This information content structure can then be refined by a two-step update.

2.4 Constructing the reduced-dimension Jacobian matrix

In an inverse system with a known native-resolution Jacobian matrix K, a reduced-dimension Jacobian matrix K_RD can be constructed on a multiscale grid that maintains the native resolution where the information content is highest and aggregates grid cells elsewhere (top right panel of Fig. 1). We refer to the state vector elements of this multiscale grid as clusters. An optimal multiscale grid maximizes the total DOFS and the averaging kernel sensitivities of each state vector element, referred to here as the DOFS per cluster. To construct this grid, we first define the state vector as a single element encompassing the inversion domain. We then add the native-resolution grid cells with the highest averaging kernel sensitivities to the state vector one by one, removing them from the original state vector element. For each new element x_i, we calculate the corresponding Jacobian matrix column $\partial \mathit{y}/\partial x_i$ and the resulting increase in DOFS. When the DOFS stabilize, we add instead clusters of two or more native-resolution grid cells and repeat this procedure. Clusters can be generated by, for example, K-means clustering, which aggregates spatially proximate grid cells. We repeat this process, increasing cluster size, until all native-resolution grid cells are allocated to the multiscale grid and the corresponding reduced-dimension Jacobian matrix K_RD is constructed. The DOFS convergence criteria and the sequence of cluster sizes can be selected to achieve the desired state vector dimension.

We apply this approach, beginning with our initial estimate K⁽⁰⁾ (Sect. 2.3), in a two-step update that iteratively improves the multiscale grid. Algorithm 1 describes this process in detail. Briefly, the information content for the initial multiscale grid is given by A⁽⁰⁾, which we use to identify the grid cells with the highest information content and construct a multiscale grid as described above. We compute the corresponding reduced-dimension Jacobian matrix ${\mathbf{K}}_{\mathrm{RD}}^{\left(\mathrm{1}\right)}$, introducing information content from the forward model. We identify the state vector elements where the forward model contributes the most information content by comparing the sensitivities given by the updated reduced-dimension averaging kernel matrix ${\mathbf{A}}_{\mathrm{RD}}^{\left(\mathrm{1}\right)}$ to the sensitivities given by A⁽⁰⁾. We disaggregate the clusters with the largest differences and update the reduced-dimension Jacobian, generating ${\mathbf{K}}_{\mathrm{RD}}^{\left(\mathrm{2}\right)}$. Convergence is rapid, and we find no need for further iteration. The analytical inversion can then be solved exactly on the multiscale grid with ${\mathbf{K}}_{\mathrm{RD}}^{\left(\mathrm{2}\right)}$.

2.5 Constructing the reduced-rank Jacobian matrix

In an inverse system with a known native-resolution Jacobian matrix K, a reduced-rank approximation of the Jacobian matrix K_Π can be constructed by calculating the linear relationship between emissions and observations for the most important patterns of information content rather than for individual or aggregate grid cells. A low-rank Jacobian corresponds to the state space shown in the bottom left panel of Fig. 1. We showed that the leading patterns of information content are given by the columns of the dimension-restoring transformation Γ^∗ (Eq. 8). For any selected value of k, the k leading patterns span a rank-k, dimension-n subspace of the original information content space. A Jacobian matrix can be constructed within this space by calculating the model response to perturbations of these patterns. The response of the forward model F to the jth normalized eigenvector ${\mathit{\gamma }}_j^{\ast }\in {\mathbb{R}}^n$, given by the jth column of Γ^∗, is as follows:

where β is any scalar sufficiently large to ensure numerical stability. The model responses ${\mathit{y}}_j,j\in \left\{\mathrm{1},\mathrm{\dots },k\right\}$, form the columns of the matrix ${\mathbf{K}}_{\mathit{\omega }}\in {\mathbb{R}}^{m\times k}$, which is the Jacobian matrix for an inverse system with a reduced-dimension state space spanned by the first k eigenvectors of the information content, illustrated by the bottom right panel of Fig. 1. This reduced-dimension Jacobian must be transformed to the original state dimension to enable physical interpretation of the posterior results. Bousserez and Henze (2018), following Bocquet et al. (2011), show that the reduced-dimension Jacobian matrix K_ω is given by ${\mathbf{K}}_{\mathit{\omega }}=\mathbf{K}{\mathbf{\Gamma }}^{\ast }$ and the reduced-rank Jacobian matrix K_Π by ${\mathbf{K}}_{\mathrm{\Pi }}=\mathbf{K}\mathbf{\Pi }=\mathbf{K}{\mathbf{\Gamma }}^{\ast }\mathbf{\Gamma }$. Thus, the reduced-rank Jacobian can be calculated from the reduced-dimension Jacobian by K_Π=K_ωΓ. The resulting Jacobian has dimension m×n and rank k.

In an inverse system without a known Jacobian matrix, the reduced-rank Jacobian matrix approximation can be constructed in a two-step update that iteratively improves the patterns of information content used as perturbations. Algorithm 2 describes this process in detail. Briefly, we use the initial estimate of the Jacobian matrix K⁽⁰⁾ (Sect. 2.3) to calculate the corresponding averaging kernel matrix A⁽⁰⁾ and the matrix of its eigenvectors ${{\mathbf{\Gamma }}^{\ast }}^{\left(\mathrm{0}\right)}$. When calculating ${{\mathbf{\Gamma }}^{\ast }}^{\left(\mathrm{0}\right)}$, we select the k⁽⁰⁾ eigenvectors that have a signal-to-noise ratio greater than some threshold. We use the signal-to-noise criterion, which is stricter than the information content criterion, to account for the errors in the initial estimate of the information content. We compute the forward model response to each of the eigenvectors using Eq. (12) and transform the resulting reduced-dimension Jacobian ${\mathbf{K}}_{\mathit{\omega }}^{\left(\mathrm{1}\right)}$ to the full-dimension state space with ${\mathbf{K}}_{\mathrm{\Pi }}^{\left(\mathrm{1}\right)}={\mathbf{K}}_{\mathit{\omega }}^{\left(\mathrm{1}\right)}{\mathbf{\Gamma }}^{\left(\mathrm{0}\right)}$. We calculate the associated averaging kernel matrix ${\mathbf{A}}_{\mathrm{\Pi }}^{\left(\mathrm{1}\right)}$ and the matrix of its eigenvectors ${{\mathbf{\Gamma }}_{\mathrm{\Pi }}^{\ast }}^{\left(\mathrm{1}\right)}$. Because ${\mathbf{K}}_{\mathrm{\Pi }}^{\left(\mathrm{1}\right)}$ is a reduced-rank approximation, its spectrum of information content is discontinuous at k⁽⁰⁾. We, therefore, use the spectrum of information content associated with the initial, full-rank estimate A⁽⁰⁾ to select the rank k⁽¹⁾ of the second update and calculate ${{\mathbf{\Gamma }}_{\mathrm{\Pi }}^{\ast }}^{\left(\mathrm{1}\right)}$. We use the k⁽¹⁾ eigenvectors that span most of the information content from the initial estimate and construct an updated reduced-rank Jacobian matrix approximation ${\mathbf{K}}_{\mathrm{\Pi }}^{\left(\mathrm{2}\right)}$ as above. The resulting Jacobian matrix ${\mathbf{K}}_{\mathrm{\Pi }}^{\left(\mathrm{2}\right)}$ is a rank ≈k⁽¹⁾ approximation that accurately quantifies the forward model where the observing system has high information content, and loses accuracy in areas with lower information content, where the observations are least able to constrain emissions. The resulting posterior solution is accurate in areas with high information content and defaults to the prior estimate elsewhere.