2.1 Model framework

The aim is to infer some parameters of interest, here a spatial field of a priori emissions scaled by some factor, x, from some measurement. The a priori emissions at a given location are generally informed by spatially resolved bottom-up inventories (as in Sect. 3.1.3) or extrapolation from some reported emissions value. The value x is a multiplicative scaling of this a priori value for emissions, most generally expressed as a quantity of gas per unit time per unit area. The approach taken in this work uses a Gaussian Markov random field for fast and efficient calculation of spatial correlation for x. This means that the emissions field is required to be a latent Gaussian field, which will be discussed further in Sect. 2.2. Net surface fluxes of many greenhouse gases are positive at the scales resolved by the model. In this work, due to the usage of a latent Gaussian field, which must be defined over both positive and negative values, we choose to look at the deviations of emissions from the prior mean emissions field. This is in an effort to fit the physical model to the imposed statistical model that fast computation requires. Taking the approach that, for many regional inverse problems involving longer-lived trace gases, there is a linear relationship between emissions which are constant in time and observed atmospheric concentration, the relationship between measurements and emissions is

where y is a vector of the residual between the measured and a priori predicted measurement; H is a Jacobian (or sensitivity) matrix, which maps the surface emissions to the measurements; u is a vector of independent and identically distributed variables containing the contribution to the measurement of mole fractions at the boundary of the domain minus the prior mean contribution, with an associated sensitivity matrix K; and ϵ is some stochastic error. The variable x has a Gaussian prior probability with zero expectation and a covariance described by a Gaussian random field. This will be solved using a hierarchical framework (Sect. 2.4).

2.2 Gaussian Markov random fields

The emissions scaling from its a priori value x is spatially continuous over the domain of interest. We assume that it exhibits a spatial correlation structure because emissions at one location are generally not independent from all other locations in the same field. We choose to model the covariance in this field using a Matérn covariance function, which Stein (1999) shows is well-suited to natural systems due to its flexibility, with other correlation structures (e.g. exponential) being special cases of the Matérn family (Guttorp and Gneiting, 2006). A stationary Gaussian random field with a Matérn covariance is the solution to a particular stochastic partial differential equation (Whittle, 1954, 1963), given by

where κ is the spatial scale parameter, τ influences the variance, Δ is the Laplacian operator, Ω is the spatial domain and 𝒲(s) is a Gaussian stochastic process for locations s, noting that the dependence is only on the Euclidean distance; the process is therefore isotropic. For an example of non-stationary and anisotropic fields see Marques et al. (2019). The smoothness parameter α gives a continuous domain Markov field for integer values and is set to α=2 (see Whittle, 1954). Smaller values will give more short-scale variability and can be difficult to differentiate from noise. Lindgren et al. (2011) show that if the field is represented using a Gaussian Markov random field (see Rue and Held, 2005), then it is very efficient to directly construct the precision, or inverse covariance, matrix using a basis function representation

where ψ_k(⋅) is a piecewise linear function in each triangle of a mesh construction of the spatial domain, where ψ_k(⋅) is 1 at vertex k and 0 at all other vertices. This mesh is created using constrained refined Delaunay triangulation (Shewchuk, 2002), placing nodes at the main points of interests and infilling the rest of the space using some condition of minimum and maximum length of the vertices. In this work we choose to represent the UK and Ireland with an evenly spaced denser mesh with a coarser mesh outside of this region. A mesh could be further refined, for example by creating an even denser mesh close to the measurement site where sensitivity is higher. It is important to extend the mesh beyond the region where the measurements are sensitive to emissions as the mesh is constructed using Neumann boundary conditions, which trigger reflection and therefore overestimation close to the boundaries. This has no effect on the inferred emissions, provided that the mesh is extended far enough around the domain of interest. Figure 1 shows the mesh used in this work, excluding the extended outer mesh region.

Figure 1A mesh constructed using constrained refined Delaunay triangulation, where the distribution of nodes is denser around the United Kingdom and Ireland.

2.3 Extension to spatio-temporal correlation

A spatio-temporal extension to the forward model (Eq. 1) is possible by including the spatial correlation structure introduced in Sect. 2.2 in a temporal framework (Cameletti et al., 2013). Similar to Eq. (1), the deviation from the prior mean measurement at site l made at time t in its simplest form is

where i represents each of the total n nodes and j represents the edge of the domain at each of the four cardinal directions. We make the assumption that measurements made at a given time are independent, giving a vectorised observation vector y_t at each time. This can be generalised to include correlated measurements, where the covariance function should result in a sparse matrix to retain computational efficiency. Following this, the matrix H from Eq. (1) becomes the sparse block diagonal matrix

which operates on the vectorised spatio-temporal scaling of the emissions field $\mathit{x}=[{\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{2}},\mathrm{\dots }{\mathit{x}}_m]$ to model the observations $\mathit{y}=[{\mathit{y}}_{\mathrm{1}},{\mathit{y}}_{\mathrm{2}},\mathrm{\dots }{\mathit{y}}_m]$. The time varying structure of H_t and x applies also to K_t, and u_t. We impose the temporal correlation structure between emissions over time to be a first-order autoregressive model, where

and

where ϕ is the temporal correlation and ${\mathbf{Q}}_S^{-\mathrm{1}}$ is the spatial correlation structure described by a Matérn field using the stochastic partial differential equation approach for a Gaussian Markov random field. We have chosen a first-order autoregressive model as we believe that emissions are generally similar to those at the previous time step. In the given model, ϕ is the correlation between the previous time step and the current time step. This means that the emissions at time t have a similarity of ϕ to emissions at t−1 plus some spatially correlated random effect δ_t. The vectorisation of x allows a separable covariance structure for the temporal and spatial covariances, which means that the spatio-temporal precision matrix can be expressed using a Kronecker product as follows (Mardia et al., 1979):

Estimating hourly emissions at each time t soon makes inference prohibitive due to burden. Instead we make the assumption that emissions are constant over a 3-month period to reduce the computational size. We continue with the assumption that measurements within a single time period are independent, although this can be generalised if required. In this experiment we make the assumption that emissions are constant over a 3-month period and that the correlation between these 3-month periods is autoregressive of the order of 1.

2.4 Hierarchical model

Inferring the emissions and the related uncertainties requires a hierarchical model to infer the quantities of interest from measurements while estimating some unknown parameters which are necessary for inference. The main focus of this work is to estimate the posterior distribution of the emissions field x based on observations y. We follow a typical Bayesian hierarchical framework

where θ is vector of hyperparameters describing the variances and covariances in x, u and y, noting that $\mathit{x}⟂⟂\mathit{y}\mid \mathit{\theta }$. We assume

where the precision matrix of the model–measurement uncertainty Q_y contains the hyperparameter for the standard deviation of the model–measurement standard deviation σ_y, σ_BC is the standard deviation of the prior for u and I is the identity matrix. Together the hyperparameters make the vector $\mathit{\theta }=(\mathit{\rho },\mathit{\sigma },\mathit{\varphi },{\mathit{\sigma }}_{\text{BC}},{\mathit{\sigma }}_y)$, which have independent prior distributions. The hyperparameters for the spatial precision matrix Q(θ) are transformations of the variables in Eq. (2),

where ρ is the range parameter and σ is the marginal standard deviation of the latent field (Lindgren et al., 2011). We use penalised complexity priors to define the prior probabilities for these parameters (Simpson et al., 2017; Fuglstad et al., 2018). Penalised complexity priors allow the formation of priors when there is only a vague understanding of their true values while enforcing more constraint than using a broad uniform or a Jeffreys prior. This uses the information loss of deviating from some baseline estimate of the parameter. Penalised priors do not increase the computational speed for the given case and are chosen instead for their intuitiveness. Penalised complexity priors are specified by the probability of the parameters being less than or greater than some baseline value

where ρ₀ and σ₀ are the baseline values and p_ρ and p_σ are the associated probabilities defined by the user. The hyperparameter ϕ controls the temporal correlation between latent variables, and its prior probability is defined on a beta distribution scaled between −1 and 1 as

where a and b are assumed coefficients. The matrix Q_y is the precision matrix for the combined measurement and model errors. The diagonal of Q_y contains the square of the hyperparameter σ_y, where the prior probability is defined on $\log \frac{\mathrm{1}}{{\mathit{\sigma }}_y^{\mathrm{2}}}$ as

The marginal standard deviation of the a priori boundary conditions σ_BC is also constructed this way, giving

A Bayesian hierarchical model requires a method of inference to estimate the parameters of interest and any parameters that are not of direct interest but required, and uncertain, in order to infer the other parameters in the hierarchy. Many methods of inference exist and have been applied to the problem of estimating emissions of trace gases (see references in Sect. 1). A promising approach using the correlation structure in Sect. 2.2 and 2.3 is INLA (Lindgren and Rue, 2015). Section 2.5 outlines this approach to inference.

2.5 Inference using an integrated nested Laplacian approximation

An integrated nested Laplacian approximation (Rue et al., 2009) provides a fast and efficient framework to infer the latent variable x and hyperparameters θ from measurements y. The calculation of the INLA is possible using the R-INLA package (Lindgren and Rue, 2015). In this work we use R-INLA version 17.06.20. The speed in this approach comes from solving the marginal posteriors for x_i, i.e. each element i of the latent field, through the numerical integration

where j is the jth element in θ and −j indicates all by the element j. Equations (11) and (12) make use of a Laplace approximation by approximating p(θ∣y) by

where $p_{\text{G}}(\mathit{x}\mid \mathit{\theta },\mathit{y})$ is a Gaussian approximation to the full conditional of x and the right-hand side is evaluated at x^∗(θ), which is the modal probability of x for a given θ. This approximation is exact if p(θ∣y) is Gaussian and gives a good approximation for log concave problems (Tierney and Kadane, 1986). Then, it is possible to approximate $p(x_i\mid \mathit{y},\mathit{\theta })$ using another Laplace approximation

where $p_{\text{G}}({\mathit{x}}_{-i}\mid x_i,\mathit{\theta },\mathit{y})$ is the Gaussian approximation to ${\mathit{x}}_{-i}\mid x_i,\mathit{\theta },\mathit{y}$ evaluated at the mode ${\mathit{x}}_{-i}^{\ast }(x_i,\mathit{\theta })$. See Rue et al. (2009) and Martins et al. (2013) for a more in-depth description.

While this method relies on calculating only the marginal posterior distributions of x_i, it is still possible to predict a linear combination of the field to provide regional emissions totals (e.g. country totals). We define a linear predictor of emissions for a given region η^*, defined using the basis function representation of the mesh where each node contains information on the spatial area represented by that node and its connecting vertices and zeros for all nodes that are not in the region of interest. Then we can approximate the linear combination of the parameters of interest, giving a combined emissions total η^*, using Eqs. (11) and (12) and transforming the predicted latent field by the weightings containing the area information.

3.1 Transport model and data

3.1.1 Measurement data

Figure 2The case studies use the four UK DECC network measurement stations located on this map. RGL is Ridge Hill, TAC is Tacolneston, TTA is Angus and MHD is Mace Head station.

The case study observations are from four measurement sites: three in the UK and one in Ireland, which are part of the UK Deriving Emissions related to Climate Change (DECC) network (Stanley et al., 2018). Figure 2 shows the location of these four measurement stations: Ridge Hill in the west of England, Angus on the east coast of Scotland, Tacolneston on the east coast of England and Mace Head on the west coast of Ireland. Measurements are made quasi-continuously throughout this network but are averaged into hourly values here, consistent with the time step of the meteorology that drives the atmospheric transport model (Sect. 3.1.2). The data set contains ∼ 10 000 measurements to demonstrate the capabilities of the method at handling moderate data volumes. We consider the scalability of this method in the discussion (Sect. 4).

3.2 Pseudo-data experiment

We test the method by performing an inversion using pseudo-data for four consecutive time periods of 1 month. By creating a known-emissions field we are able to validate the method through comparing the inferred emissions to the known emissions, which is not possible in the real world. We form a synthetic-emissions field by allowing the emissions to deviate from the prior mean emissions according to a Matérn field (see Sect. 2.2). We choose to simulate the data using σ=0.5 as the uncertainty in the EDGAR v4.3.2. inventory is around 50 % for methane emissions and ρ=3.25, which is similar to the correlation length scale in UK emissions in the EDGAR v4.3.2. inventory estimated using a variogram (Cressie, 2015), although the correlation length scale in the uncertainty is unknown. We use ϕ=0.8 as expert experience suggests that UK emissions of methane are generally highly correlated in time.

To create the synthetic-emissions field, the NAME sensitivities for each measurement at each grid cell, detailed in Sect. 3.1.1, are multiplied by the corresponding EDGAR inventory emissions in that grid cell, which are then transformed into the triangulation nodes in Fig. 1. This forms the matrix H. We randomly sample the full spatio-temporal precision matrix for the latent field (Sect. 2.3) using the GMRFLib library (Rue and Follestad, 2001) to generate the latent field x. In this experiment we treat the boundary conditions as known as in practice these are generally well-constrained (or treated as known) during an inversion. The observations are simulated using the simulated latent field and sensitivities following Sect. 2.1 with additive Gaussian noise with a standard deviation of 15 % of y. The synthetic observations contain a total of 11 520 measurement points, which are used to infer 646 emissions nodes for each time period, i.e. the nodes of the mesh in Fig. 1. Figure 3a–d show the synthetic-emissions field in terms of deviation from some prior emissions field. The synthetic total deviation in emissions from the prior mean for the UK for each period is 0.07, 1.29, 1.01 and −0.33 Tg yr⁻¹.

Figure 3Time-correlated deviation of emissions from the prior mean, simulated using a Matérn field for (a) the first time period, (b) the second time period, (c) the third time period and (d) the fourth time period.

The inference needs prior probabilities for the hyperparameters, which are known exactly here, but we set them to be deliberately incorrect, but feasible based on true prior knowledge, to check that the inversion method can still recover the correct emissions. For the inversion we assign a prior probability for ϕ using a=6.5 and b=0.1, for σ using σ₀=0.1 and p_σ=0.01 and for ρ using ρ₀=5 and p_ρ=0.5. We base the constraint on the spatio-temporal emissions on the assumption that methane emissions in the UK are likely to be strongly correlated between time periods, that they are unlikely to vary by more than 10 % of the a priori emissions from the previous time step and that there is little knowledge of the spatial correlation structure. For the model measurement error we assign the prior probability on the log precision as $\log \frac{\mathrm{1}}{{\mathit{\sigma }}_y^{\mathrm{2}}}\sim \mathcal{N}(-\mathrm{5},\mathrm{1})$, which represents approximately 68 % probability that the error falls between 8 and 20 ppb, in line with previous experience.

Figure 4The inferred mean deviation of emissions from the prior mean for (a) the first time period, (b) the second time period, (c) the third time period and (d) the fourth time period.

Figure 4 shows the resultant inferred mean changes in the emissions field for each time period using INLA (Sect. 2.5) and the hierarchical model in Sect. 2.4. The change in the emissions field is the inferred latent field scaling parameter multiplied by the inventory emissions in Sect. 3.1.3. Figure S1 shows maps of the 95 % uncertainty in the difference in emissions. A linear combination of the posterior latent field multiplied by the inventory emissions field gives total emissions for the UK. The inferred posterior mean change in total emissions for the UK from the prior mean with their associated 95 % uncertainties for the four time periods are 0.01 [−0.10, 0.12], 1.23 [1.10, 1.36], 1.07 [0.97, 1.18] and −0.30 [−0.41, −0.18] Tg yr⁻¹, which are shown in Fig. 5. The inferred emissions for the UK as a whole agree well with the synthetic emissions, with all true synthetic emissions totals falling within the 95 % uncertainty of the inferred total emissions. These synthetic tests cover both small and large differences between the a priori estimate of methane emissions in the UK and their true value. The larger differences are most likely an overestimate of what would be expected in practice, although this provides a good test case. The posterior mean spatial distribution of emissions in Fig. 4 is qualitatively similar to those in Fig. 3, although we avoid reading heavily into comparisons of the mean estimates plotted on spatial maps (see Gelman and Price, 1999).

Hyperparameter estimation is less accurate than for the latent field. The estimation of the noise σ_y generally captures the imposed noise well, which had a mean value of 5.5 ppm and an estimated value of 6.6 [6.5, 6.8] ppm. The correlation structures are less well-captured. The temporal correlation ϕ was estimated with a mean value of 0.7 [0.6, 0.8]; the range ρ has a value of 1.7 [1.3, 2.1] and the marginal standard deviation of the latent field σ has a value of 0.8 [0.7,0.9]. It is promising that the posterior mean estimates of all hyperparameters show an improvement on their prior mean or baseline values, although only the true value for ρ falls within the estimated 95 % uncertainty.

Figure 5The inferred and true difference between UK emissions and a synthetic inventory value for each of the four time periods. The crosses show the inferred mean difference in total emissions with their associated 95 % uncertainty. The red line indicates one-to-one agreement.

Download

3.3 Real-data experiment

This section presents methane emissions estimates for the UK in 2014. The year is split into four time periods: January to March, April to June, July to September and October to December. Based on the synthetic data set in Sect. 3.2, we use the same prior probabilities for the hyperparameters as the inversion. The prior boundary conditions are distributed with ${\mathit{\mu }}_{{\mathit{\sigma }}_{\text{BC}}}=\mathrm{3.2}$ and ${\mathit{\sigma }}_{{\mathit{\sigma }}_{\text{BC}}}=\mathrm{0.4}$, based on expert judgement. The a priori values for these boundary conditions come from the MOZART-4 model as in Sect. 3.1.3. The linear mapping from the latent field to the measurements is generated using the NAME-derived sensitivities described in Sect. 3.1.2 multiplied by the inventory emissions detailed in Sect. 3.1.3.

Figure 6Inferred mean difference in methane emissions for the UK in 2014 compared to the EDGAR inventory for (a) January to March, (b) April to June, (c) July to September and (d) October to December.

Figure 6 shows maps of the inferred mean difference in emissions from the a priori inventory for the four time periods using INLA as in Sect. 2.5. This result is the mean posterior scaling for the latent field multiplied by the inventory value (see Sect. 3.1.3). Figure S2 shows maps of the 95 % uncertainty in the difference in emissions. A linear combination of the posterior latent field multiplied by the inventory emissions field gives total emissions for the UK. We estimate total UK methane emissions in 2014, with their associated 95 % uncertainty, of 2.14 [1.90, 2.38], 2.28 [2.06, 2.51], 2.62 [2.39, 2.84] and 2.09 [1.87, 2.30] Tg yr⁻¹ for the periods January to March, April to June, July to September and October to December, respectively. These are plotted along with the inventory estimate in Fig. 7. This emissions trend suggests that, for 2014, there was an increase in methane emissions in the UK during the summer months compared to the winter months. The uncertainties, however, are large, meaning that this increase may not be as stark as suggested by the mean estimates. Combining these emissions into a mean annual emission for 2014 with its associated 2-standard-deviation uncertainty assuming that the time periods are correlated with the modal posterior value of ϕ gives 2.28±0.33 Tg yr⁻¹. The emissions are similar to mean UK estimates from previous hierarchical inversions using NAME of 2.09 [1.65, 2.67] Tg yr⁻¹ by Ganesan et al. (2015) and 2.28 [2.04, 2.52] Tg yr⁻¹ by Lunt et al. (2016). All of these estimates are broadly in line with the United Kingdom's 2014 emissions estimate of 2.13 Tg yr⁻¹ reported in its national inventory (Department of Business, Energy and Industrial Strategy, 2019).

Figure 7Estimated total UK methane emissions for 2014. The prior mean comes from EDGAR v4.32 inventory data (red dashed). For comparison, the figure also shows the posterior means for UK methane emissions from Ganesan et al. (2015) (black dotted) and Lunt et al. (2016) (black dash-dotted).

Download