Spatiotemporal Geostatistical Analysis and Global Mapping of CH₄ Columns from GOSAT Observations

Abstract

Methane (CH₄) is one of the most important greenhouse gases causing the global warming effect. The mapping data of atmospheric CH₄ concentrations in space and time can help us better to understand the characteristics and driving factors of CH₄ variation as to support the actions of CH₄ emission reduction for preventing the continuous increase of atmospheric CH₄ concentrations. In this study, we applied a spatiotemporal geostatistical analysis and prediction to develop an approach to generate the mapping CH₄ dataset (Mapping-XCH₄) in 1° grid and three days globally using column averaged dry air mole fraction of CH₄ (XCH₄) data derived from observations of the Greenhouse Gases Observing Satellite (GOSAT) from April 2009 to April 2020. Cross-validation for the spatiotemporal geostatistical predictions showed better correlation coefficient of 0.97 and a mean absolute prediction error of 7.66 ppb. The standard deviation is 11.42 ppb when comparing the Mapping-XCH₄ data with the ground measurements from the total carbon column observing network (TCCON). Moreover, we assessed the performance of this Mapping-XCH₄ dataset by comparing with the XCH₄ simulations from the CarbonTracker model and primarily investigating the variations of XCH₄ from April 2009 to April 2020. The results showed that the mean annual increase in XCH₄ was 7.5 ppb/yr derived from Mapping-XCH_4, which was slightly greater than 7.3 ppb/yr from the ground observational network during the past 10 years from 2010. XCH₄ is larger in South Asia and eastern China than in the other regions, which agrees with the XCH₄ simulations. The Mapping-XCH₄ shows a significant linear relationship and a correlation coefficient of determination (R²) of 0.66, with EDGAR emission inventories over Monsoon Asia. Moreover, we found that Mapping-XCH₄ could detect the reduction of XCH₄ in the period of lockdown from January to April 2020 in China, likely due to the COVID-19 pandemic. In conclusion, we can apply GOSAT observations over a long period from 2009 to 2020 to generate a spatiotemporally continuous dataset globally using geostatistical analysis. This long-term Mpping-XCH₄ dataset has great potential for understanding the spatiotemporal variations of CH₄ concentrations induced by natural processes and anthropogenic emissions at a global and regional scale.

1. Introduction

Atmospheric methane (CH₄), as one of the most important greenhouse gases, is second only to carbon dioxide (CO₂) in contributing to global warming [1]. The global atmospheric CH₄ concentration has increased from a preindustrial level of 722 ppb to approximately 1750 ppb in 2000 and continues to increase [2,3,4,5] due to the influence of anthropogenic emissions such as fossil fuel combustion, agricultural planting, livestock husbandry, biological emissions, etc. [6,7]. Countries all over the world have begun to take actions to reduce CH₄ emission to prevent the continuous increase of atmospheric CH₄ concentrations [2]. It is necessary to develop a spatiotemporally continuous CH₄ data on a global scale to aid in understanding the spatial and temporal variabilities in CH₄ emissions so as to support the actions of CH₄ emission reduction [8,9].

Ground-based measurement networks have provided long-term and high-precision greenhouse gas data [10,11]. These stations, however, are sparsely and unevenly distributed and highly costly [12]. Satellite observations, which can obtain the column averaged dry-air mole fraction of CH₄ (XCH₄), have the advantages of global coverage and high observation density. It has become an important way to obtain global and regional atmospheric CH₄ [13]. Satellite measurements of atmospheric CH₄ mainly include the thermal and near-infrared sensor for carbon observation Fourier transform spectrometer (TANSO-FTS) operated on the Greenhouse Gases Observing Satellite (GOSAT) [14], the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) onboard the Environmental Satellite (ENVISAT) [15], the tropospheric monitoring instrument (TROPOMI) onboard the European Space Agency’s Sentinel-5 Precursor [16], and the infrared atmospheric sounding interferometer (IASI) instrument flying onboard the MetOp platform [17]. GOSAT, launched in 2009, is a first spacecraft specially measuring the atmospheric CO₂ and CH₄ column abundances [18].

GOSAT observations, which are sensitivity in detecting the CH₄ variation in space and time, has been providing long-term XCH₄ retrievals starting from 2009 as an important data source. The XCH₄ retrievals from GOSAT observations have been produced and released as L2 products by the National Institute for Environmental Studies (NIES) team [18]. XCH₄ are retrieved from an optimal estimation inversion algorithm using the shortwave infrared spectrum observed by TANSO-FTS. These XCH₄ retrievals, however, irregularly distribute globally with many gaps in space and time due to GOSAT observing mode [19] such as return period, sampling mode of footprints etc., and the effects of cloud and aerosols [20]. It is difficult to reveal the fine variations of XCH₄ in space and time by these XCH₄ retrievals with many gaps and in irregular distribution.

Several previous studies have indicated that geostatistical approaches [21,22] can effectively fill in the gaps in satellite observations utilizing correlation structures within greenhouse gas observations and linear unbiased prediction capabilities [23,24,25]. Spatial-only kriging, a conventional geostatistical method, was widely adopted to generate the XCH₄ mapping data product. Liu et al. [26] applied spatial kriging interpolation method to obtain the daily mean value of greenhouse gas concentration in East Asia using the L2 product data of XCO₂ and XCH₄ retrievals by GOSAT satellite. They only used a single year of the XCH₄ retrievals to model the spatial variation to fill those gaps. This conventional spatial-only geostatistical method, which only makes use of spatial correlation, does not consider the temporal correlation structure of the CH₄ variations, including the annual increase and seasonal cycles. A geostatistical analysis has been extended from space to time [27,28,29] with the accumulation of available data and the development of spatiotemporal variation function. Liu et al. [30] developed a the space-time kriging interpolation approach based on the geostatistical method to fill gaps and generated the spatiotemporal continues XCH₄ data in Monsoon Asia [30]. This research indicated that the spatiotemporal geostatistical method showed better performance than the spatial-only geostatistical method. The spatiotemporal geostatistical approach not only provides a probabilistic framework for data analysis and prediction [21] by utilizing the joint spatial and temporal dependences between observations, but also applies a larger dataset in both space and time to support stable parameter estimation and prediction to ensure high accuracy prediction [31]. Previous these studies, however, mostly focused on the regional scale and short periods. A mapping XCH₄ data doesn’t developed on a global scale and in long-term yet currently. GOSAT observations has been providing the XCH₄ retrievals in a long period from 2009 to now globally. We need to develop an approach to generate the spatiotemporal continuous XCH₄ on a global scale using these the available GOSAT long-term XCH₄ retrievals.

In this study, we aim to develop an approach to fill the gaps in XCH₄ retrievals and generate a global XCH₄ mapping dataset (Mapping-XCH₄) continuously in space and time based on spatiotemporal geostatistics using XCH₄ retrievals from GOSAT observations from 2009 to 2020. The data used in this study and the methodology of mapping XCH₄ based on spatiotemporal geostatistics are described in Section 2. The assessments of the generated Mapping-XCH₄ dataset, including validation of the Mapping-XCH₄ dataset, comparison with CH₄ simulations, and variation of XCH₄ revealed by the Mapping-XCH₄ dataset are shown in Section 3. The performance of the Mapping-XCH₄ dataset in revealing the spatiotemporal characteristics of XCH₄ corresponding to CH₄ emissions over Monsoon Asia is discussed in Section 4. Finally, the conclusion is described in Section 5.

2. Materials and Methods

2.1. Data Used

2.1.1. GOSAT XCH₄ Retrievals

We collected the GOSAT XCH₄ retrievals data as Level 2 product (v2.81) released for general users spanning from April 2009 to April 2020, which covers the global land area. We filtered these XCH₄ data product using the “quality_flag” tag in the data product documents to get the valid datasets. It has been reported that these XCH₄ retrievals have a mean global bias of 1.9 ppb with a standard deviation of 13.4 ppb compared with ground-based observations of the Total Carbon Column Observing Network (TCCON) [32]. Figure 1 shows the density of the collected XCH₄ retrieval points from April 2009 to April 2020. As shown in Figure 1, there are many gaps in XCH₄ retrievals and the density of data over high latitude and low latitude regions is smaller than that in other regions (Figure 1a). Additionally, the density of retrievals in summer is less than that in winter because of limitations from clear sky conditions and solar zenith angles.

2.1.2. Data for Validation and Analysis

We collected TCCON data and the CH₄ simulations from Carbon Tracker-CH₄ for the validation of the accuracy of Mapping-XCH₄ dataset. Additionally, we collected the surface CH₄ emissions data in 2010 from the Emission Database for Global Atmospheric Research (EDGAR).

TCCON, a global network of ground-based Fourier transform spectrometers (FTS) established for the validation of near-infrared total column measurements which has been extensively used for validation of satellite observations [11,33]. The instruments have high accuracy with approximately 0.5% (5 ppb) error in XCH₄ retrievals. We collected the available TCCON data released in the 2014 version over 22 sites shown in Figure 1a from 2015 to 2019 to validate the Mapping-XCH₄ data.

CarbonTracker-CH₄ is released by the National Oceanic and Atmospheric Administration (NOAA). CarbonTracker introduces the assimilation of global atmospheric CH₄ observations from surface air samples and tall towers [34,35]. The CarbonTracker-CH₄, CT2010, produces estimates of global atmospheric CH₄ mole fractions and surface-atmosphere fluxes released from 2000 to 2010. The model-simulated CH₄ concentration is gridded CH₄ at 4° (latitude) and 6° (longitude) with 25 vertical layers before 2006 and 34 vertical layers after 2006 [36] and with a temporal resolution of 3 h. In comparison with ground-based observations, the model outputs show a mean bias of −10.4 ppb [35]. The CH₄ mole fraction profile data from CarbonTracker-CH₄ are converted to XCH₄ (CT-XCH₄) by using the pressure-average method described by Connor et al. [37]. The averaging kernel effect [38] is not considered in the conversion since it has been indicated that the difference is less than 0.1% for XCO₂ if averaging kernel smoothing is applied to model simulations or not [39].

To assess the performance of Mapping-XCH₄ further, we collected anthropogenic emission data in a 0.1 grid from EDGAR (version 6.0, globally) from 2010 to 2018, which are released by the European Commission’s Joint Research Centre and Netherlands Environmental Assessment Agency [40,41]. EDGAR provides emissions of the three main greenhouse gases (CO₂, CH₄, and N₂O) and fluorinated gases per sector and country. The data are mainly from point source emissions and the global energy statistics database of the International Energy Agency (IEA). The EDGAR CH₄ emission data include the emissions from the energy industry, fuel transformation/nonenergy, agriculture, solid waste disposal, fossil fuel fires, large-scale biomass burning, etc. The EDGAR data is a fundamental reference for many studies of surface emissions and have been widely used [42].

2.2. Mapping XCH₄ Based on Spatiotemporal Geostatistics

We developed an approach to generate global mapping XCH₄ data in a grid of 1° by 1° and intervals of 3 days from 2009 to 2020 by applying spatiotemporal geostatistics to XCH₄ retrievals derived from GOSAT observations. Figure 2 shows the main processing steps in the approach, including the analysis of spatiotemporal trend of XCH₄, modeling of the correlation structure in space and time, and prediction based on space-time kriging method and generating mapping XCH₄. The Mapping-XCH₄ data are generated by the optimal prediction of the variable at an unsampled location and time of satellite observation [43,44].

In this study, we only used the XCH₄ retrievals with land fractions larger than 90 because the GOSAT retrievals over the ocean from glint mode observations have not been fully validated [45]. We implemented this processing for five mainland subregions in a global land, including Eurasia, Africa, North America, South America, and Oceania. This partition facilitates the processing and geostatistical mapping of XCH₄ on a global scale.

2.2.1. Modeling the Spatiotemporal Trend and Correlation Structure of XCH₄

XCH₄ can be represented by a variable $Z=\left\{Z\left(s,t\right)\right|s\in S,t\in T\}$ that varies within a spatial domain S and time interval T. Its spatiotemporal variations can be separated into a deterministic trend component $m\left(s,t\right)$ and a stochastic residual component $R\left(s,t\right)$, as shown in Equation (1).

$$Z\left(s,t\right)=m\left(s,t\right)+R\left(s,t\right)$$

(1)

Equations (2) and (3).

$$m\left(s,t\right)=\left[{\mathrm{a}}_0+{\mathrm{a}}_1\ast s\left(1\right)+{\mathrm{a}}_2\ast s\left(2\right)\right]+\left[{\mathrm{a}}_3\ast t+{\displaystyle {\sum }_{i=1}^4\left({\beta }_i\sin \left(i\omega t\right)+{\gamma }_i\cos \left(i\omega t\right)\right)}\right]$$

(2)

where $w=(2\ast \pi )/T$, $s(1)$ and $s(2)$ are the latitude and longitude of the spatial location, respectively. $t$ is the marking sequence number of the corresponding time unit, ${\mathrm{a}}_0$ is the background XCH₄ at the starting time, and ${\mathrm{a}}_{1–3}$, ${\beta }_{1–4}$ and ${\gamma }_{1–4}$ are parameters to be estimated by the least-squares technique [26]. As shown in Equation (2), the deterministic spatial trend in $m\left(s,t\right)$, depending on the spatial distribution of CH₄ sources and sinks [46], is modeled by a linear surface function [46]. The deterministic temporal trend in $m\left(s,t\right)$, including the interannual XCH₄ increase and the inherent seasonal XCH₄ cycle, is determined by a set of annual harmonic functions with a simple linear function.

The stochastic residual component $R\left(s,t\right)$ in Equation (1), which represents the local variability not explained by the deterministic trend component, is derived by subtracting the trend components $m\left(s,t\right)$ from the XCH₄ data $Z$.

According to the period of GOSAT observation with three days, we set every three days as a time unit, and the annual period $T$ is 122 time-units. Figure 3 shows the trend components modeled in space and time. We can find an increasing trend in space from north to south and a temporal trend with a clear annual increase and a seasonal cycle from Figure 3.

In spatiotemporal geostatistical analysis, the optimal kriging prediction of $Z\left(s_0,t_0\right)$ at an unobserved position $\left(s_0,t_0\right)$ can be calculated as the linear weighted sum of the XCH₄ values that minimizes the mean squared prediction error [22]. The weights for the observations are determined by the distribution of observations and the variogram model, which characterizes the spatiotemporal correlation structure of the data. Therefore, estimating and modeling the variogram are crucial steps in kriging prediction for the global land mapping [46].

The experimental variogram value $\stackrel{\wedge }{\gamma }\left(h_s,h_t\right)$ of XCH₄ data at spatial lag $h_s$ and temporal lag $h_t$ is given by:

$$\stackrel{\wedge }{\gamma }\left(h_s,h_t\right)=\frac{1}{2N\left(h_s,h_t\right)}{{\displaystyle {\sum }_{i=1}^{N\left(h_s,h_t\right)}\left[Z\left(s_i+h_s,t_i+h_t\right)-Z\left(s_i,t_i\right)\right]}}^2$$

(3)

where $Z\left(s_i,t_i\right)$ is the observation at a spatial location and in a time $\left(s_i,t_i\right)$. $N\left(h_s,h_t\right)$ is the number of data pairs within a distance of $\left(h_s,h_t\right)$. Once the experimental variogram has been constructed, a spatiotemporal variogram model $\gamma \left(h,u\right)$ to fit it. The reason for such a model fit is to make sure the variance model is positive definitive which is required in the calculation of kriging prediction, as is described in Section 2.2.2 below. The spatiotemporal variogram model $\gamma \left(h,u\right)$ adopted here is an inseparable combination of the product-sum model [47] and an extra global nugget to capture the nugget effect [48], as given by Equation (4):

$$\gamma \left(h,u\right)={\gamma }_S\left(h\right)+{\gamma }_T\left(u\right)-K{\gamma }_S\left(h\right){\gamma }_T\left(u\right)+N_{ST}$$

(4)

where ${\gamma }_S\left(h\right)$ and ${\gamma }_T\left(u\right)$ are the exponential marginal variograms in space and time, respectively, $K$ is the fixed spatiotemporal binding parameter to be estimated, and $N_{ST}$ is the global base station value. The admissible value for this $K$ is dependent on the sill values of the marginal variograms, namely, $0<K<1/\max \left\{sill{\gamma }_S\left(h\right);sill{\gamma }_T\left(u\right)\right\}$. $sill{\gamma }_S\left(h\right)$ and $sill{\gamma }_T\left(u\right)$ represent the partial sills of the exponential marginal variograms in space and time, respectively. We used the iterative nonlinear weighted least-square method to estimate all parameters [49].

Figure 4 presents the different experimental variogram (left) for the five regions calculated from all GOSAT retrievals by Equation (4) and its fitted variogram model (right) by the product-sum model for five regions. In the variograms shown in Figure 4, the semi-variance values are called nugget and sill when spatial and temporal lag are closet to 0 and infinity, respectively. The ratio of nugget to sill, expressed as a percentage, can be used as an indicator to classify data dependence [50]. In general, a ratio larger than 0.75 indicates a strong spatial or temporal dependence.

Table 1. The Parameters of the Spatiotemporal Variogram Model and the Corresponding Nugget/Sill Ratios.

XCH4 Retrievals	Spatial Lag (km)	Temporal Lag (Days)	Nugget/Sill
Eurasia	300	28	0.43
Africa	600	54	0.53
North America	900	45	0.66
South America	700	36	0.43
Oceania	1800	24	0.63

presents the parameters of spatiotemporal variation in five regions. The ratio values of nugget to sill for all regions shown in Table 1 are less than 0.75, which indicates an overall spatial and temporal dependence. The ratio in Eurasia and South America are the lowest which indicates stronger spatial and temporal dependence than the other three regions.

2.2.2. Generating the Mapping-XCH₄ Dataset Using Space-Time Kriging

Based on the modeled spatiotemporal variogram, space-time Kriging estimates an arbitrary target point at unobserved location $\left(s_0,t_0\right)$ from the stochastic residual component $R=\left\{R\left(s,t\right)\right|s\in S,t\in T\}$. Supposing that $\stackrel{\wedge }{R}\left(s_0,t_0\right)$ is a prediction of $R\left(s_0,t_0\right)$, we estimated $\stackrel{\wedge }{R}\left(s_o,t_o\right)$ as a linear weighted sum of residual data within a kriging neighborhood [22,46] in space and time relative to the prediction location $\left(s_0,t_0\right)$ as shown in Equation (5).

$$\stackrel{\wedge }{R}\left(s_0,t_0\right)={\displaystyle {\sum }_{i=1}^n\left[{\lambda }_{i}R\left(s_i,t_i\right)\right]}$$

(5)

where $n$ is the number of observations to be used. ${\lambda }_i$ is the weighting factor assigned to a known observation $R\left(s_i,t_i\right)$ so as to minimize the prediction error variance while maintaining unbiased prediction. The weighting factor is calculated by Equation (6).

$$\lambda ={\left({\gamma }_0+\mathbf{1}\frac{\mathbf{1}-{\mathbf{1}}^T{\mathrm{\Gamma }}^{-1}{\gamma }_0}{{\mathbf{1}}^T{\mathrm{\Gamma }}^{-1}\mathbf{1}}\right)}^T{\mathrm{\Gamma }}^{-1}, {\displaystyle {\sum }_{i=1}^n{\lambda }_i}=1$$

(6)

where $\lambda ={[{\lambda }_1,{\lambda }_2,{\lambda }_3,\dots ,{\lambda }_n]}^T$, $T$ is the transpose of matrix, $\mathrm{\Gamma }\left(i,j\right)=\gamma \left(\left|s_i-s_j\right|,\left|t_i-t_j\right|\right)$ represents the matrix of variogram values between observations, and ${\gamma }_0\left(i,1\right)=\gamma \left(\left|s_i-s_0\right|,\left|t_i-t_0\right|\right)$ represents the column vector of variogram values between observations and predictions.

The variance of prediction error, which is a measurement of prediction uncertainty, is shown in Equation (7)

$${\sigma }^2={\gamma }_0{}^T{\mathrm{\Gamma }}^{-1}{\gamma }_0-\frac{{\left({\mathbf{1}}^T{\mathrm{\Gamma }}^{-1}{\gamma }_0-\mathbf{1}\right)}^2}{{\mathbf{1}}^T{\mathrm{\Gamma }}^{-1}\mathbf{1}}$$

(7)

where ${\sigma }^2$ is the variance of kriging prediction, and 1 is the n × 1 unit vector. The deterministic trend component (the temporal trend of XCH₄ as shown in Equation (2) and Figure 3) is tightly constrained since a lot of data available. As indicated by the fitted parameters, the uncertainties are very small compared to the predictive variance in the residual component. Therefore, the uncertainty in the regression mean term is not considered here.

To reduce the computational complexity and maintain local variability, data used in the prediction were searched within an appropriate space-time kriging neighborhood centered on the predicting point [22]. The kriging neighborhood used is a moving cylinder in space and time as described in Zeng et al. [44]. The radii of the search range are initially set to 300 km in space and 20 time-units. The 10 km and 1 time unit for each search process is used as the increment lags if the number of observations is less than 20 within the cylinder neighborhood. At the same time, we set the search range radii limits to 500 km and 40 time-units.

2.3. Precision Evaluation of the Mapping-XCH₄ Dataset

We evaluated the performance of the mapping XCH₄ dataset using cross-validation and compared it with TCCON measurements and CH₄ data simulated by the atmospheric transport model.

Cross-validation, which is a widely used method for evaluating the prediction accuracy of statistical models [51], can be used to evaluate the prediction accuracy of the spatiotemporal geostatistical method. Cross-validation is implemented by first removing the observation data and then making its prediction using the remaining data. As a result, two datasets, the predicted dataset $\left\{\stackrel{\wedge }{Z}\left(s_j,t_j\right),j=1,2,\dots ,N\right\}$ and the corresponding original dataset $\left\{Z\left(s_j,t_j\right),j=1,2,\dots ,N\right\}$, can be obtained. We selected four evaluation criteria, the correlation coefficient (R) between two datasets, the mean absolute prediction error $\left\{MAE=\frac{1}{N}{\displaystyle \sum _{j=1}^N\left|Z\left(s_j,t_j\right)-\stackrel{\wedge }{Z}\left(s_j,t_j\right)\right|}\right\}$ and the percentage of prediction error less than 15 ppb, and the population mean prediction error $\left\{ME=\frac{1}{N}\left({\displaystyle \sum _{j=1}^{N}Z\left(s_j,t_j\right)-\stackrel{\wedge }{Z}\left(s_j,t_j\right)}\right)\right\}$ to evaluate the results [52].

Additionally, we also discussed the performance of the Mapping-XCH₄ dataset in detecting the spatial and timely variations from 2009 to 2020 over Monsoon Asia.

3. Evaluation of the Mapping-XCH₄ Dataset from 2009 to 2020

3.1. Precision of the Mapping-XCH₄ Dataset

3.1.1. Results of Cross-Validation

Figure 5 shows the biases of the predicted values comparing with the corresponding observed values by the cross-validation. The results of cross-validation, as shown in Figure 5, present a significant correlation coefficient (R²) of 0.97 which indicate small prediction biases between the observed XCH₄ and the predicted XCH₄ data, where the mean MAE is generally 7.76 ppb and the standard deviation (STD) is 6.91 ppb. Figure 5a demonstrates a systematic bias where the observed XCH₄ tends to be larger than the predicted XCH₄ mostly at high concentrations, while it is smaller than the predicted XCH₄ mostly at lower concentrations. This is likely related to the effects of the sensor capability and sensitivity in GOSAT observations for areas with extremely high and low concentrations. The statistical results of cross-validation for the five regions are show in

Table 2. Results of Cross-Validation for Mapping-XCH₄. N is the Number of Validated Sample Pairs.

Area	N (*104)	R	MAE (ppb)	Percent (MAE < 15 ppb)	ME
Eurasia	44.28	0.9408	8.6319	84	0.0146
Africa	40.76	0.9636	6.9318	91	0.0021
North America	13.26	0.9314	8.8291	83	−0.0072
South America	13.26	0.9604	7.7473	88	−0.1034
Oceania	11.91	0.9456	6.2717	93	0.0188

. It can be seen from Table 2 that the number of predicted data with MAE less than 15 ppb accounts for greater than 83 percent of all validated samples (N). The biases in Eurasia and North America are larger than those in the other regions likely due to the diverse and active surface XCH₄ emissions in the two areas.

3.1.2. Uncertainty of Mapping-XCH₄

We can obtain the corresponding kriging variance as indicated in Equation (7) for each geostatistical prediction to evaluate the uncertainty of Mapping-XCH₄ dataset. The kriging variance depends on both the density of XCH₄ retrievals and the data homogeneity surrounding the prediction position. The denser the observations and more homogeneous the XCH₄ variation surrounding the prediction position are, the lower the prediction uncertainty [53].

Figure 6a demonstrates the spatiotemporal variation in the kriging standard deviations (Kstd), which is the root of the kriging variance of the Mapping-XCH₄ dataset from 2009 to 2020. Figure 6a shows that the uncertainty of Mapping-XCH₄ is larger in mid-high latitudes, 35°N–65°N, and around the tropical region, 10°S, where Kstd ranges from 14 ppb to 17 ppb, and seasonally presents a maximum in December or January each year in these areas. In mid-low latitudes, the uncertainty is lower, and Kstd presents a maximum in July or August each year. These high uncertainties are likely due to fewer available observations, which can be seen in Figure 1, less homogeneous variation in XCH₄ surrounding the prediction location induced by the effects of atmospheric conditions such as clouds, water vapor, aerosols, etc., and observation of the geometry of the sun-target-sensor geometry.

The magnitude of annual average Kstd, which is shown in Figure 6b, generally presents larger and strip variation in Eurasia compared with the other regions. In particular, Mapping-XCH₄ shows higher uncertainties in the southern area of Asia, where Kstd is up to 16 ppb. These results are likely induced by the available GOSAT observations, as shown Figure 1, and inhomogeneous XCH₄ variation and atmospheric conditions in these areas.

3.1.3. Comparison with TCCON Measurements

XCH₄ from TCCON measurements has high accuracy and can be used for validation of the Mapping-XCH₄ [54]. In this study, we selected 22 representative TCCON sites where XCH₄ measurements were available from 2015 to 2019. We took 1° × 1° grid and one month as a space-time unit to calculate the mean values of TCCON. Figure 7 shows the biases of Mapping-XCH₄ by the comparison of the XCH₄ measurements (X-axis) over TCCON sites plotted in different colors and the Mapping-XCH₄ (Y-axis) that coincide with TCCON sites in geolocation and time. As seen from Figure 7, the correlation coefficient (CORR), indicating consistency between Mapping-XCH₄ and TCCON data, is as high as 0.92, and the mean absolute error (MAE), mean error (ME), and standard deviation (STD) are 8.1 ppb, 1.07 ppb, and 11.42 ppb, respectively.

Overall verification results by TCCON sites (Appendix A 

Table A1. Bias Statistics between the Monthly Average Mapping-XCH₄ and TCCON XCH₄, which are Calculated as Their Differences for each Coincident Data Pair and Averaged for each Site.

Sites	Location (Latitude, Longitude)	Coincident Data Pairs	MAE (ppb)	ME (ppb)	Reference
Lauder	(−45.04°N, 169.68.5°E)	53	5.93	−3.92	[55]
Wollongong	(−34.41°N, 150.89°E)	55	6.33	−5.57	[56]
Darwin	(−12.43°N, 130.89°E)	44	3.07	−1.85	[57]
Hefei	(31.90°N, 117.17°E)	16	7.47	−4.45	[58]
Saga	(33.24°N, 130.29°E)	57	7.48	5.59	[59]
Caltech	(34.14°N, 118.13°E)	35	4.57	1.75	[60]
Jet Propulsion Laboratory	(34.20°N, 118.18°E)	37	5.01	1.04	[61]
Edwards	(34.96°N, 117.88°E)	44	6.47	2.10	[62]
Nicosia	(35.14°N, 33.38°E)	44	7.71	3.53	[63]
Tsukuba	(36.05°N, 140.12°E)	57	6.68	4.90	[64]
Anmeyondo	(36.54°N, 126.33°E)	26	7.36	1.67	[65]
Lamont	(36.60°N, 97.49°E)	60	7.05	−6.10	[66]
Rikubetsu	(43.46°N, 143.77°E)	57	6.92	4.33	[67]
Park Falls	(45.94°N, 90.27°E)	60	4.17	−0.30	[68]
Zugspitze	(47.42°N, 10.98°E)	51	14.60	10.30	[69]
Garmisch	(47.48°N, 11.06°E)	58	4.38	−0.41	[70]
Orléans	(47.97°N, 2.11°E)	55	6.58	−1.68	[71]
Paris	(48.85°N, 2.36°E)	48	10.61	−9.42	[72]
Karlsruhe	(49.1°N, 8.44°E)	59	6.40	−2.00	[73]
Bremen	(53.1°N, 8.85°E)	52	8.27	−1.17	[74]
Bialystok	(53.23°N, 23.02°E)	44	7.77	−3.22	[75]
East Trout Lake	(54.36°N, 104.99°E)	36	6.87	−4.22	[76]
Overall	-	1048	8.10	1.07	-

) [55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76], the MAE is generally 8.10 ppb, while the largest MAEs are over Paris and Zugspitze located in Eurasia, which are 10.61 ppb and 14.60 ppb, respectively. The reasons for these larger biases are likely related to the uncertainty of the Mapping-XCH₄ dataset as described above.

3.2. Global XCH₄ Variations Revealed by the Mapping-XCH₄ Dataset XCH₄ in 1° × 1° grids from 2009 to 2020

Figure 8 shows the spatial pattern of the annual average XCH₄ derived from the Mapping-XCH₄ dataset from 2009 to 2020. The local enlarged annual mean Mapping-XCH₄ for Monsoon Asia from 2010 to 2020 are shown in Figure A1 in the Appendix A. As shown in Figure 8 and Figure A1 in the Appendix A, the highest XCH₄ presents in eastern China and Southeast Asia, which is related to high emissions induced by human activities related to fossil fuel usage and paddy agriculture. It is known that the area of paddy fields in Southeast Asia, including eastern China, accounts for more than 50 percent of the global paddy fields. Spatially, the largest varying amplitude shows in China, where XCH₄ generally varies from 1820 ppb to 1860 ppb from the northwest to southeast, and its difference is up to 60 ppb in winter.

The XCH₄ around tropical areas in South America and Africa shows higher XCH₄ concentrations, which is likely related to wetland emissions. This may also be associated with the uncertainty of Mapping-XCH₄ in the tropical area of South America (Figure 6b) caused by much smaller number of available GOSAT XCH₄ retrievals (Figure 1a).

Figure 9 shows the latitudinal and temporal variation in XCH₄, which is computed by monthly averaged XCH₄ within a 1° latitudinal band using the Mapping-XCH₄ dataset in 1° × 1° grids and intervals of three days from April 2009 to April 2020. It can be seen from Figure 9 the Mapping-XCH₄ demonstrates an obvious spatial change depending on latitude and a temporally increasing trend in all latitudinal bands. The higher the latitude is, the lower the XCH₄ is. Mapping-XCH₄ is higher in the Northern Hemisphere than in the Southern Hemisphere due to much higher anthropogenic CH₄ emissions in the Northern Hemisphere. High XCH₄ values are also present in the equator at approximately 15°N with high temperature and many wetlands, which largely impacts the enhancement of XCH₄ surface emissions.

Figure 10 shows the timely variation in globally average XCH₄ from 2009 to 2020. The globally average XCH₄ presents an annual increase from 2009 to 2020 and seasonal variation which are in agreement with the Global Atmosphere Watch (GAW) report [77]. The maximum XCH₄ appeared in November and December, and the minimum appeared in June and July. The seasonal amplitude was up to 11.4 ppb. The seasonal variation in XCH₄ is partly driven by the abundance of hydroxyl radicals (OH), the most important sink of CH₄, in the atmosphere. As a result, a high abundance of OH in the summer drives down XCH₄ [78,79]. Regionally, it could be related to the growth cycle of paddies and vegetation, as well as the change in temperature and humidity caused by the influence of monsoons [80]. Higher temperature and humidity will increase CH₄ emissions. Wetlands, such as swamps and lakes, have significant CH₄ emissions due to anaerobic environments [80,81].

The mean annual increase in XCH₄ was 7.5 ppb/yr during the eleven years from April 2009 to April 2020, which was slightly greater than the mean growth rate of 7.3 ppb/yr during the ten years from Jan 2010 to Dec 2019 reported by the GAW report based on an in situ observational network [77].

3.3. Comparison with Model Simulations of XCH₄

Figure 11 presents the mean XCH₄ in 2010 derived from the Mapping-XCH₄ dataset and the model simulating XCH₄ by CarbonTacker (CT-XCH₄) in the same year, and Figure 12 shows their difference. The Mapping-XCH₄ generally presents the similar spatial pattern to the CT-XCH₄ (Figure 11). The Mapping-XCH₄ in Southeast Asia, Northern India and the rainforest in Africa and Southern America, however, are 20–57 ppb larger than the CT-XCH₄ as shown in Figure 12. These regions are associated with strong effects of surface emissions, mostly from paddy fields and wetlands in southeastern Asia and tropical areas, and some fossil fuel exploitation and landfills. Mapping-XCH₄ in central India, the central area of Northern Asia, North America, and South America is 10–37 ppb smaller than the modeled XCH₄.

The large difference in the middle of Africa is likely due to the uncertainty of model simulation in this tropical area rather than the Mapping-XCH₄ as the Mapping-XCH₄ shows the lowest uncertainties in this area shown in Figure 6b. As a result, the large differences between Mapping-XCH₄ and CT-XCH₄ mostly arise in the wetland and the paddy fields, which imply that the CT-XCH₄ likely have large uncertainty there due to the flaw of prior surface emission in simulating [35].

4. Discussion

4.1. Spatiotemporal Characteristics of Mapping XCH₄ Corresponding to the Surface Emissions

The variations in CH₄ concentration primarily depend on surface emissions and partly on sinks from chemical reactions in the atmosphere and soil, in addition to atmospheric transport by global and local circulations with seasonality. Surface emissions include anthropogenic emissions from agriculture and waste, fossil fuel production and use, natural emissions from wetlands, geological lakes, termites, biomass burning, etc. The emissions are mostly from agriculture and waste (34%), wetlands (30%), and fossil fuel production and use (19%) [82]. The significantly high correlation between Mapping-XCH₄ and the annual surface CH₄ emissions in 2010 from the EDGAR emission inventory was shown by Liu et al. [26]. Here, we primarily applied the Mapping-XCH₄ dataset from 2009 to 2020 generated by the approach above to reveal the spatiotemporal characteristics of XCH₄ corresponding to the surface emissions at the regional scale.

Figure 13 demonstrates the spatial pattern of annual mean XCH₄ in 2010 derived from the Mapping-XCH₄ and CT-XCH₄, respectively, the anthropogenic emissions from EDGAR and the uncertainty (Kstd) of Mapping-XCH₄ in the same year. Comparing Mapping-XCH₄ with CT-XCH_4, as shown in Figure 13a,b, we can find that Mapping-XCH₄ reveals a finer spatial variation than CT-XCH₄. The spatial pattern of Mapping-XCH₄, moreover, is in better agreement with the CH₄ emissions of EDGAR (Figure 14) than that of CT-XCH₄. In particular, the Sichuan Basin in China, Bangladesh, and northern India with large emissions corresponds to high values in Mapping-XCH₄. These large CH₄ emissions are mostly related to paddy fields, fossil fuel production and use, and the stagnant effect of the basin terrain [83]. Satellite observations of CO columns and NO₂ columns indicated the large magnitude of industrial emissions over this region also [84]. The Mapping-XCH₄ does not respond to surface emissions only in the southern tip of the peninsula of South Asia where the large uncertainty of Mapping-XCH₄ (Figure 13d) due to a smaller number of available XCH₄ retrievals (Figure 1a). These results indicate that the Mapping-XCH₄ can reveal and capture the local CH₄ emissions in space and time as the Mapping-XCH₄ is generated by instantaneous satellite observations that respond to surface emissions [85,86].

Additionally, Figure 14 shows the correlation between surface emissions and Mapping-XCH₄ in Monsoon Asia in 2010. Similarly, a significant linear relationship was also observed by the Mapping-XCH₄ data, with a coefficient of determination (R²) of 0.66 and a p value less than 0.01. These results indicated that the variation in CH₄ concentration was greatly affected by the surface emissions.

4.2. Temporal Variations of XCH₄ for Various Surface Emissions

We selected the regions with typical surface emissions which locations are marked as S1–S7 in Figure 13a,c, to investigate the variation of XCH₄ revealed by the Mapping-XCH₄ datasets. Figure 15 presents a yearly increasing trend from 2009 to 2020 and the difference of annual enhancements for these sampling regions.

It can be seen from Figure 15a that the Tibetan Province (S7) and Taklimakan Deserts (S6) with the lowest emissions (Figure 14), present the minimum XCH₄ values of 1790 ± 6.1 ppb and 1808 ± 2.9 ppb, respectively, which are less than the average for all of Chinese land (1818 ± 23.6 ppb). Chongqing (S1), Hunan province (S2), Shanghai city (S3), and the Sichuan Basin (S4) nearly show the highest XCH₄, 1858 ± 4.0 ppb, 1857 ± 3.0 ppb, 1853 ± 2.4 ppb, and 1851 ± 8.1 ppb, respectively, which is mostly related to the emissions from the dense oil and gas fields, and paddy fields in addition to the atmospheric accumulation effects due to the geographic basin structure [85]. The XCH₄ in Shanxi Province (S5) with high emissions shows the highest standard deviation (9.8 ppb), which is likely induced by many point sources of small or mesoscale coal mines scattered in this region.

The Taklimakan Desert (S6) has the smallest amplitude in seasonal cycle due to less human activity and without vegetation. The Tibetan Province (S7) has the largest amplitude, which is likely induced by the XCH₄ emissions from wetlands and meadow grasslands tending to be sensitive to seasonal variations in temperature on the Tibetan Plateau above 4000 m above sea level, the reason of which needs further investigation.

Figure 15b shows the difference of annual enhancement between years during the three months from January to April for the sampling regions. It can be found from Figure 15b that the annual enhancement in 2020 is lower than that in 2019 for the regions except the Shanxi Province (S5) and Taklimakan Deserts (S6), which is likely related with the reduction of anthropogenic emissions caused by the lockdown during this period from January to April in 2020 due to the coronavirus disease 2019 (COVID-19) [87,88]. The NO₂ concentration in Shanxi Province during this lockdown period showed the similar change to XCH₄ [89], which is likely because the needs of heating in winter did not decrease; thus, there was no significant reduction of emissions from the fossil fuel production and use.

As a result, primarily evaluating the reasonability and potential of the mapping-XCH₄ dataset for detecting the variations and the evidence of CH₄ above, find (1) the mean annual increase of XCH₄ and seasonal variation derived from the Mapping-XCH₄ dataset is agreement with that from in-situ observational network; (2) the Mapping-XCH₄ dataset could explain the variable of anthropogenic emissions with significant correlation between the Mapping-XCH₄ and EDGAR emission inventory (R² = 0.66) over Monsoon Asia; (3) the Mapping-XCH₄ can be used to detect spatiotemporal characteristics regionally, and the evidence of CH₄ variations induced by the special events, such as the decrease of XCH₄ during January–April in 2020 in China caused by reduction of anthropogenic emission due to the lock down of COVID-19 pandemic.

5. Conclusions

In this study, we proposed a data-driven approach based on a spatiotemporal geostatistical model to generate a global land Mapping-XCH₄ dataset in 1° × 1° grids and intervals of 3 days from 2009 to 2020 using XCH₄ retrievals derived by GOSAT observations. This Mapping-XCH₄ dataset shows better precision with a high correlation coefficient, 0.97, a small mean absolute prediction error of 7.66 ppb in the cross-validation, and good agreement with TCCON sites with a standard deviation of 11.42 ppb. We evaluated the performance of the Mapping-XCH₄ dataset by primarily investigating the spatial pattern and timely variation of XCH₄ and comparing it with the model simulations. The results show that the timely variations in XCH₄ characterized by the Mapping-XCH₄ dataset are generally in agreement with the characteristics based on the ground measurements and significantly correlate with the EDGAR emission inventory. The spatial patterns of XCH₄ revealed by the Mapping-XCH₄ dataset correspond to the distribution of surface CH₄ emissions from paddy fields, wetlands, and anthropogenic emissions. Moreover, the Mpping-XCH₄ dataset presents a finer spatial pattern than the model simulations. These results demonstrated that the Mapping-XCH₄ dataset could help us to investigate the spatiotemporal patterns of XCH₄ at global and regional scales. The Mapping-XCH₄ dataset has the advantage of being spatiotemporally continuous compared with the original XCH₄ retrievals with many gaps in space and time. The Mapping-XCH₄ dataset with long-term in grid, such as shown in Figure A1 in Appendix A, facilitates the detection of the driving factors of CH₄ variations combined with other satellite observation data, such as ecological parameters of vegetation and land cover related to CH₄ natural emissions, and inventory data of agriculture, wetlands, and anthropogenic emissions. This could be like the application of global mapping XCO₂ data which are generated by GOSAT observations using geostatistical analysis as well [90,91].

This study assumed that the spatiotemporal correlation structure is similar over the entire processing area. However, this is not always true, as the spatiotemporal variations in different locations are different. The irregular distribution of the original XCH₄ retrievals prevents us from resolving this problem. Therefore, the precision of Mapping-XCH₄ is mostly due to the number of available GOSAT observations where the more observations there are, the better the geostatistical modeling. Moreover, it also depends on the accuracy of the original XCH₄ retrievals derived from satellite observations. Mapping-XCH₄ in tropical rainforest areas and high latitudes should receive more attention due to the few observations and limitations of observation conditions there. The uncertainty of Mapping-XCH₄ will hopefully be reduced along with increasing XCH₄ data available from multiple satellites to improve the geostatistical model.