Sustained methane emissions from China after 2012 despite declining coal production and rice-cultivated area

We use the version 7.2 proxy nadir retrievals of GOSAT methane column data from the University of Leicester (Parker et al 2015) in our inverse analysis. GOSAT retrieves the atmospheric methane column by nadir measurements of solar backscatter (1.65 µm absorption band) (Kuze et al 2016). Observations are made at three circular pixels of 10 km diameter across the orbit track 260 km apart, separated by 260 km along the track. The same locations are sampled every three days. In China, GOSAT retrieves about 1000 observations per month, but there are more data in the west because of less cloud cover. The number of GOSAT retrievals over a given location does not vary significantly year-to-year and is similar for different seasons, though it is slightly larger during October–December in China (see figure S8 in SI). The ability of GOSAT to constrain regional methane emissions has been shown by previous high resolution inverse studies for India, Brazil, North America, North Africa, etc (Ganesan et al 2017, Tunnicliffe et al 2020, Maasakkers et al 2021, Western et al 2021).

We also include ground-based hourly measurements at (43.2^∘N, 145.5^∘E, 96 m above sea level) and Hateruma (24.1^∘N, 123.8^∘E, 46.5 m above sea level) stations in Japan operated by the Institute for Environmental Studies (NIES), which have been described in detail by Tohjima et al (2014). The two stations are relatively insensitive to China's emissions, but help to improve the constraints to boundary conditions.

We use the UK Met Office NAME (Numerical Atmospheric dispersion Modelling Environment) model, a Lagrangian particle dispersion model (Manning et al 2011), as the atmospheric transport model used to provide the relationship between emissions and concentrations in the atmosphere. The model domain is 54^∘E–170^∘W, 5^∘S–84^∘N. We derive an optimized estimate of spatially resolved methane emissions in China using Bayesian inverse analysis. The inversion minimizes the cost function $J(\mathbf{x})$ by solving $\nabla_\mathbf{x}{\,}J(\mathbf{x}) = 0$, with $J(\mathbf{x})$ defined as follows:

$$\begin{align} J(\mathbf{x}) &= \frac12(\mathbf{x} - \mathbf{x}_\mathrm{prior} )^T \mathbf{P}^{-1} (\mathbf{x} - \mathbf{x}_\mathrm{prior} ) \nonumber\\&\quad+ \frac{1}{2}(\mathbf{y} - \mathbf{H} \mathbf{x})^T \mathbf{R}^{-1} (\mathbf{y} - \mathbf{H} \mathbf{x}). \end{align} \tag{ 1 }$$

Here y is the vector of observations, H is the Jacobian matrix representing the sensitivities of observations to changes in the state vector x and $\mathbf{x}_\mathrm{prior}$ is the the prior value of x. P is the prior error covariance matrix and R is the observational error covariance matrix.

The state vector contains 250 elements for aggregated methane emissions per month (200 elements for China) within the model domain using a Gaussian Mixture Model (GMM) with radial basis functions based on spatial proximity and source type patterns (Turner and Jacob 2015). The resulting clustering is shown in figure S10 in SI. The use of the GMM enables us to retain high resolution for major sources up to native model resolution while coarsening resolution for weak or broadly distributed sources. The choice of the number of emission elements in China is based on the suggestions from Turner and Jacob (2015) that balance the aggregation and smoothing errors. The state vector also includes additional four elements for monthly boundary concentrations at the four edges of the model domain to reduce the influence of potential biases in boundary conditions on our derived fluxes (Lunt et al 2016, Tunnicliffe et al 2020) and one element for an offset parameter between satellite and surface data. The prior boundary concentrations at the domain edges for the Lagrangian particle dispersion model are from a GEOS-Chem $4^\circ\times5^\circ$ global simulation using methane emissions optimized with GOSAT satellite data Maasakkers et al (2019). Here our prior bottom-up inventory is not pre-optimized by the global inversion. The offset parameter (20 ppb a priori) accounts for any systematic differences between the GOSAT and the NIES ground-based measurements and their representations by the model Tunnicliffe et al (2020). Furthermore, we note that our inversion results are not significantly affected by changing the numbers of emission elements for China (within Turner et al suggested range) and the elements for the boundary concentrations (see figure S3 and S10 for in SI for details). We assume 100% prior error for each emission state vector element without correlation (i.e. off-diagonal elements in P are zero), 2% prior error for boundary conditions, and 50% prior error for the offset parameter. Previous studies (Maasakkers et al 2016, Bloom et al 2017, Zhang et al 2020) have suggested that error correlation is mainly related to natural wetlands (i.e. there are no large error correlations in anthropogenic emissions), while wetland emissions account for less than 5% of China's total emissions and are a minor source in our inversion. We also note that adding small error correlations (5%–20%) for emissions in P do not change our results (see figure S13 in the SI). Therefore taking a diagonal prior covariance matrix here is appropriate. Observational errors include GOSAT instrument, model transport, and representation errors. We use 10 ppb as the mean observational error (model error + instrument error) standard deviation derived by previous NAME inverse modeling against ground-based observations (Manning et al 2011). The GOSAT mean instrument error standard deviation is 11 ppb (Parker et al 2015), indicating that satellite instrument error dominates the observational error. For a given satellite observation, we take the maximum of the reported GOSAT instrument error and NAME model error as the observational error added in quadrature to the observational error covariance matrix.

The sensitivity matrix H is computed by NAME. We release particles at a rate of 2000 particles per hour for each GOSAT retrieval vertical level, over a one minute period centered around the retrieval time, and trace them back in time for 30 d, as done by previous studies (Ganesan et al 2017, Tunnicliffe et al 2020). For surface sites, we release 20 000 particles at the release height of the stations (100 m and 50 m above sea level for Cape Ochiishi and Hateruma, respectively). Particles are lost when they reach the surface or the four edges of the model domain. The sensitivity matrix H changes over time depending on different measurement locations and varying meteorology. The OH impact on regional inversions here is negligible. The majority of the sensitivity of a measurement to the surrounding emissions is within 100 km of the measurement location. This implies that air masses will typically travel for a few hours to a few days in this high sensitivity region, during which time and only 0.006%–0.2% of the methane (particles) will be lost due to reaction with OH. The model is driven by the Unified Model's Model meteorology (Walters et al 2014) with horizontal resolution spanning from $0.352^\circ \times 0.234^\circ$ to $0.141^\circ \times 0.094^\circ$ over 2010-2017. The model output is set to be $0.352^\circ \times 0.234^\circ$ for the inversion to be consistent with the lowest meteorological resolution used. NAME has been used extensively to calculate sensitivity matrices for inverses analyses of long-lived greenhouse gases and ozone-depleting gases (Manning et al 2011, Ganesan et al 2014, 2017, Rigby et al 2019, Tunnicliffe et al 2020).

We perform the inversion at monthly resolution. We optimize the logarithms of the emission to ensure positivity in all grid cells. Inversions under normal and lognormal assumptions for emissions are previously found to be consistent (Maasakkers et al 2019). We found our inversion results between the two assumptions using our baseline prior inventory differ less than 0.6% on sectoral and national scales. Under the lognormal assumption, the inverse problem is non-linear and can be solved numerically by the Levenberg-Marquardt iterative algorithm (Maasakkers et al 2019) with a convergence tolerance of 0.05% for maximum relative error:

$$\begin{align} \mathbf{x}_{n+1} &= \mathbf{x}_n + (\mathbf{R}^{-1} + \mathbf{H}_n^T\mathbf{P}^{-1}\mathbf{H}_n)^{-1} (\mathbf{H}_n^T\mathbf{P}^{-1}\mathbf{H}_n)^{-1}\nonumber\\&\quad\times(y-\mathbf{H}\mathbf{x}_n) + \mathbf{R}^{-1}(\mathbf{x}_n - \mathbf{x}_\mathrm{prior}) \end{align} \tag{ 2 }$$

where $\mathbf{H}_n = \partial \mathbf{y}/\partial (\ln \mathbf{x}$) is the Jacobian matrix at the nth iteration. Each individual element $\partial y_i/\partial \ln x_j$ of H_n can be recalculated in the iteration by $\partial y_i/\partial \ln x_j = x_j \partial y_i / \partial x_j$ with $\partial y_i / \partial x_j$ being individual elements of H for linear problems computed by the forward model. Validation of inversion results with observations is presented in figures S1 and S2 in SI.

To test the inversion sensitivity to the prior boundary conditions, we performed 100 inversions by randomly perturbing the boundary conditions on each grid cell at four domain edges (referred to as 'SENS$_\mathrm{BC}$'). The perturbations followed a uniform distribution within $\pm2\%$. We find that 95% (i.e. 2σ) of the ensemble 'SENS$_\mathrm{BC}$' results differ less than 0.8% for the national totals and 0.5%–2% for different source sectors (see figure S3). We also performed a sensitivity test on the model XCO₂ by using a different GOSAT proxy product from the RemoTeC v2.3.9, which used the CarbonTracker Model to model XCO₂ (Schepers et al 2012), while the University of Leicester GOSAT proxy product used in our inversion is based on an ensemble of model XCO₂ data (Parker et al 2015) (figure S4). These sensitivity tests show that inversion results are generally insensitive to different GOSAT proxy product and random perturbations on prior boundary conditions.

To evaluate the inversion sensitivity to the prior fractional information, we perform an ensemble of inversions using 1000 different prior estimates generated by perturbing the baseline prior. This ensemble is referred to as 'SENS$_\mathrm{Prior}$' and throughout the text the ensemble mean represents our posterior estimates. For each ensemble member, the relative fraction of a given source sector is perturbed in each grid cell. This is done by randomly selecting a source from coal, rice/aquaculture, livestock, waste, oil/gas, wetlands, or other sectors within the given grid cell and then increasing or decreasing the selected sector emissions following a uniform distribution ($\pm20\%$). This procedure will also change the total emissions in each grid cell and thus the total emissions for each ensemble member. We neglect grid cells that have small emissions less than 0.5 ton CH₄ per day. We do not perturb source sectors that have a fraction of 0 or 1 in a grid cell (i.e. have no source emissions or dominate the grid cell completely). Overall there are approximately 1500 gridded emissions in China that are perturbed. An example of perturbed sector fractions for the gridded emission inventory is shown in figure S9(a) in SI. We note that the effect of doing the emission perturbation on state vector elements instead of gridded emissions is small (as shown in figure S9 in SI). We prefer to do emission perturbation on the prior $0.352^\circ \times 0.234^\circ$ gridded inventory because it is independent of choice of state vector. This ensemble sensitivity test allows us to assess the effect on source attribution due to uncertainties in relative sector fractions within grid cells (in particular those that have mixed sources). Perturbed sector fractions at grid levels could change the patterns of state vector elements determined by the Gaussian mixture model we used, but we find this change is minimal to our inversion results (see figure S10 in SI).

We also perform two additional inversions using the EDGAR v4.2 inventory and a scaled EDGAR v4.2 inventory as the prior. The latter scales coal mining emissions in EDGAR v4.2 from 28 Tg CH₄ yr⁻¹ to 16.7 Tg CH₄ yr⁻¹ to match the total emissions of coal mining in the China-specific coal inventory (Sheng et al 2019). The purpose of these additional inversions is to demonstrate that large spatial errors in coal mining can lead to biased interpretation of inversion results for China as we discussed in the main text.

We estimate China's mean annual total methane emissions over 2010–2017 to be 57.6 ± 2.4 Tg CH₄ yr⁻¹ (of which 55.2 ± 2.3 Tg CH₄ yr⁻¹ are anthropogenic emissions) (figure 1). Our posterior anthropogenic estimates are consistent with the China's 2014 UNFCCC report (UNFCCC 2020) (55.3 ± 2.9 Tg CH₄ yr⁻¹) and with the previous top-down estimates over similar time periods (Maasakkers et al 2019, Miller et al 2019, Wang et al 2019, 2021). The year-to-year variability of our posterior emission estimates is smaller than that found by Wang et al (2021). This difference could be attributed to the fact that Wang et al refer to an East Asia region that includes many other countries outside of China. In particular, these non-China regions contain wetlands, which are one of the major drivers for year-to-year variability. While our posterior estimates are only 6% lower than the prior (61.2 Tg CH₄ yr⁻¹, figure 1), there are large positive and negative absolute differences on a sub-national/provincial scale that tend to cancel each other on the national scale (figure 2(c)). The relative differences are up to $50\%$ across the country and significant error reductions are mostly over large source regions (see figure S11 in SI).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find that China's emissions exhibit a statistically significant positive linear trend (0.36 ± 0.04 Tg CH₄ yr⁻²) for 2010–2017. This trend was slowing down at a statistically significant rate of -0.1 ± 0.04 Tg CH₄ yr⁻³ indicated by the negative coefficient for the third term in a Legendre polynomial regression (see SI for details). During 2010–2012 China's emissions increase at a rate of $0.77$ ±0.2 Tg CH₄ yr⁻²), then slow after 2012 to 0.3 ± 0.1 Tg CH₄ yr⁻² for 2012–2017. Our trend estimate for 2010–2015 (0.5 ± 0.2 Tg CH₄ yr⁻²) is smaller than 1.1 ± 0.4 Tg CH₄ yr⁻² suggested by a previous study (Miller et al 2019) for the same time period. This is discussed further below.

Fractions of each source in the prior emission inventories can be used to apportion emissions and emission trends to source sectors in countries such as China where the different sources have distinct spatial or temporal distributions (Sheng et al 2018, Tunnicliffe et al 2020). Coal mining, rice/aquaculture, livestock, and waste are dominant sources for China's anthropogenic emissions. CMM emissions are mainly in Shanxi province and Southwest China, whereas emissions from rice/aquaculture dominate in Southeast China, comparable to the spatial distribution of rice paddies retrieved from the MODIS satellite data (Zhang et al 2020). Livestock emissions are mainly concentrated in a number of provinces in north-central and Southwest China, which is consistent with major livestock provinces reported in China Statistical Yearbook (National Bureau of Statistics of China 2020). Waste emissions are localized in urban areas, correlated with population. Wetland emissions are small in China and only dominate a few places in Southeast China and North China. Emissions from other anthropogenic sources and natural sources are too small to dominate any grid cells. Regions having mixed sources (i.e. no source comprising >$50\%$ of a grid cell) are mainly in Sichuan and Shandong provinces. Overall, the locations dominated by these major source sectors are well separated at the $0.352^\circ\times0.234^\circ$ model spatial resolution and account for 70% of the total emissions from China (figure 2(d)).

Coal mining is the largest contributor to the total methane emissions in China with mean 2010–2017 emissions of 14.5 ± 1.4 Tg CH₄ yr⁻¹ (table S1). The second largest source in China is rice-aquaculture with mean emissions of 13.2 ± 1.0 Tg CH₄ yr⁻¹, followed by livestock at 11.2 ± 1.1 Tg CH₄ yr⁻¹ and waste at 10.2 ± 0.8 Tg CH₄ yr⁻¹ (table S1). Trends for these sectors are discussed below, with other minor source sectors showing no significant trends (see figure S5 in SI).

CMM emissions increased by 0.7 ±0.3 Tg CH₄ yr⁻² from 2010 to 2012, driving the national trend, but flattened afterwards with a smaller growth of 0.1 ± 0.06 Tg CH₄ yr⁻² for 2012–2017 (figure 3(a)). According to activity data, coal production peaked in 2013, and in 2016/2017 returned to levels similar to those of 2010 (National Bureau of Statistics of China 2020) (figure 3(b)). The derived emission trend is consistent with coal production activities for 2010–2012 (figure 3(b)), but deviates after that. This inconsistency between top-down estimates and production may be due to emissions from abandoned coal mines. Since 2010/2011, China has consolidated its coal industry to concentrate production in the existing larger and more efficient coal mines (typically state-run mines), and to gradually close a large number of small, village or town-owned coal mines (Zhu et al 2017, Sheng et al 2019) that are widely spread across the all coal mining fields in China (figure S12 in SI). Our CMM prior inventory includes 11 000 operating coal mines in 2010/2011, of which about 4000 mines have been closed and abandoned over 2010–2017 (Sheng et al 2019) (figure 3(b)) but are still allocated in the prior inventory for the inversions. Therefore the derived trend over coal mine fields comprises the contributions from abandoned mines since 2011. An actively venting abandoned mine can emit methane up to 40%–90% of its initial rate in the first 3–4 years, shrinking to 10% after 30 years (US EPA 2004). Increasing emission factors for active coal mines may also explain the inconsistency, but reported recovery rates of CMM have been increasing (Zhu et al 2017, Sheng et al 2019, Gao et al 2020). Thus emission factors are more likely to have decreased or remained stable, which would suggest about 1.7 Tg CH₄ yr⁻¹ coming from abandoned coal mines in 2017 assuming emissions in 2017 from active coal mining are similar to those in 2010.

Zoom In Zoom Out Reset image size

Our CMM emission trend differs from that derived in the previous inverse analysis (Miller et al 2019), which uses the EDGAR v4.2 as their prior and shows a continuous increase after 2012 1.0 ± 0.3 Tg CH₄ yr⁻² (figure 3(a)). We propose that this difference is due to two factors related to the spatial pattern of emissions in the prior. First, the trend in total emissions is higher after 2013. When we use the EDGAR v4.2 as the prior in our inversion rather than the China-specific coal inventory (Sheng et al 2019), we derive a similar larger trend after 2012 in the total emissions (figures 1). The magnitude of prior CMM emissions do not significantly influence the posterior trend as indicated by our inversion results using scaled-down EDGAR v4.2 coal emissions (figure 1). This finding indicates that inaccuracies in the spatial distribution of the coal sector in EDGAR can to some extent lead to errors in the derived trend in total emissions. Second, the EDGAR v4.2 inventory at coarse resolution ($2^\circ\times2.5^\circ$) that were used by the previous global inversion (Miller et al 2019) (figure S6 in SI), show that grid cells dominated by coal or mixed sources account for about 85% of the total emissions, in contrast to the high-resolution China-specific inventory used here which has 28% of total emissions from coal dominated grid cells. The strongly differing emission patterns between the two inventories remain even at higher spatial resolution (figure S7 in SI). This would lead to a higher apportionment of the total trend to the CMM sector and consequently affect the trends of other source sectors (figures 3(c), (e), (g)). These two factors together may explain the larger derived trend in CMM emissions in the previous work, compared to our estimates.

The emission estimates from rice-aquaculture have increased by 0.13 ± 0.05 Tg CH₄ yr⁻² from 2010 to 2017. Rice emissions are proportional to paddy area (Eggleston et al 2006), but the trend in derived emissions is opposite to that of reported paddy area after 2013 (2% decline from 2013 to 2017 (National Bureau of Statistics of China 2020)). The decline of rice paddy areas is largely related to ongoing conversion of paddy fields towards industrial-scale aquaculture (Yuan et al 2019). Freshwater aquaculture in China has emission factors about four times larger than rice cultivation and is estimated to emit 3.5 Tg CH₄ yr⁻¹ based on the total area of aquaculture ponds (Yuan et al 2019), but has been overlooked in previous gridded bottom-up inventories (such as EDGAR v4.3) and inverse analyses. The freshwater aquaculture industry is primarily located in the southeast of China (figure S12 in SI) (Wang et al 2020), where the majority of rice is grown, and more than half of aquaculture facilities been converted from rice paddies (Yuan et al 2019). Therefore the estimated emissions from rice-growing regions informed by the prior include contributions from rice paddies, aquaculture ponds converted from rice paddies, and to a large extent existing aquaculture ponds. China's total area of freshwater aquaculture ponds has increased more than 20% from 2010 to 2017 (National Bureau of Statistics of China 2020) (figure 3(d)). This could explain the opposite trends between emissions and rice paddy areas.

Waste emissions inferred by the inversion exhibit a positive trend of 0.06 ± 0.04 Tg CH₄ yr⁻² between 2010 and 2017 (figure 3(e)). Landfill methane emissions are partly related to the amount and type of solid waste and management at landfills (Cai et al 2018). Due to increased urbanization, the amount to solid waste added per year in China has increased by 50% from 2010 to 2017, and its trend accelerates after 2014 (National Bureau of Statistics of China 2020) (figure 3(f)). In addition, the number of wastewater facilities has increased by 40% from 2010 to 2014 based on the most recent available data (Zhao et al 2019) (figure 3(f)), which is potentially a significant source contributing to the waste trend. However we are not able to distinguish contributions from landfills and wastewater because they are both in urban areas and largely co-located at ∼30 km spatial resolution. Waste emissions are also reversely related to local management practices (e.g. recovery rate of waste methane for power generation). This could explain that waste emissions for 2010–2017 do not increase as fast as the amount of waste and the number of facilities.

China's livestock emissions do not exhibit a significant trend (0.01 ± 0.03 Tg CH₄ yr⁻²) over 2010–2017 (figure 3(g)). Livestock emissions include enteric fermentation and manure management. The former accounts for more than 90% of China's total livestock emissions and is linearly correlated to cattle population, which does not show a significant trend over 2010–2017 but suggests a small decrease between 2010 and 2012 (figure 3(h)). The emission trend from livestock is expected to remain stable in the future because the growing meat consumption in China is due to be met by increasing imports (National Bureau of Statistics of China 2020).