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

: Methane (CH 4 ) is one of the most important greenhouse gases causing the global warming effect. The mapping data of atmospheric CH 4 concentrations in space and time can help us better to understand the characteristics and driving factors of CH 4 variation as to support the actions of CH 4 emission reduction for preventing the continuous increase of atmospheric CH 4 concentrations. In this study, we applied a spatiotemporal geostatistical analysis and prediction to develop an approach to generate the mapping CH 4 dataset (Mapping-XCH 4 ) in 1 ◦ grid and three days globally using column averaged dry air mole fraction of CH 4 (XCH 4 ) 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 coefﬁcient 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 4 data with the ground measurements from the total carbon column observing network (TCCON). Moreover, we assessed the performance of this Mapping-XCH 4 dataset by comparing with the XCH 4 simulations from the CarbonTracker model and primarily investigating the variations of XCH 4 from April 2009 to April 2020. The results showed that the mean annual increase in XCH 4 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 4 is larger in South Asia and eastern China than in the other regions, which agrees with the XCH 4 simulations. The Mapping-XCH 4 shows a signiﬁcant linear relationship and a correlation coefﬁcient of determination (R 2 ) of 0.66, with EDGAR emission inventories over Monsoon Asia. Moreover, we found that Mapping-XCH 4 could detect the reduction of XCH 4 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 4 dataset has great potential for understanding the spatiotemporal variations of CH 4 concentrations induced by natural processes and anthropogenic emissions at a global and regional scale.


Introduction
Atmospheric methane (CH 4 ), as one of the most important greenhouse gases, is second only to carbon dioxide (CO 2 ) in contributing to global warming [1]. The global atmospheric CH 4 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 the spatiotemporal continuous XCH 4 on a global scale using these the available GOSAT long-term XCH 4 retrievals.
In this study, we aim to develop an approach to fill the gaps in XCH 4 retrievals and generate a global XCH 4 mapping dataset (Mapping-XCH 4 ) continuously in space and time based on spatiotemporal geostatistics using XCH 4 retrievals from GOSAT observations from 2009 to 2020. The data used in this study and the methodology of mapping XCH 4 based on spatiotemporal geostatistics are described in Section 2. The assessments of the generated Mapping-XCH 4 dataset, including validation of the Mapping-XCH 4 dataset, comparison with CH 4 simulations, and variation of XCH 4 revealed by the Mapping-XCH 4 dataset are shown in Section 3. The performance of the Mapping-XCH 4 dataset in revealing the spatiotemporal characteristics of XCH 4 corresponding to CH 4 emissions over Monsoon Asia is discussed in Section 4. Finally, the conclusion is described in Section 5.

GOSAT XCH 4 Retrievals
We collected the GOSAT XCH 4 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 4 data product using the "quality_flag" tag in the data product documents to get the valid datasets. It has been reported that these XCH 4 retrievals have a mean global bias of 1.9 ppb with a standard deviation of 13.4 ppb compared with groundbased observations of the Total Carbon Column Observing Network (TCCON) [32]. Figure 1 shows the density of the collected XCH 4 retrieval points from April 2009 to April 2020. As shown in Figure 1, there are many gaps in XCH 4 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. spatiotemporal continuous XCH4 on a global scale using these the available GOSAT longterm XCH4 retrievals.
In this study, we aim to develop an approach to fill the gaps in XCH4 retrievals and generate a global XCH4 mapping dataset (Mapping-XCH4) continuously in space and time based on spatiotemporal geostatistics using XCH4 retrievals from GOSAT observations from 2009 to 2020. The data used in this study and the methodology of mapping XCH4 based on spatiotemporal geostatistics are described in Section 2. The assessments of the generated Mapping-XCH4 dataset, including validation of the Mapping-XCH4 dataset, comparison with CH4 simulations, and variation of XCH4 revealed by the Mapping-XCH4 dataset are shown in Section 3. The performance of the Mapping-XCH4 dataset in revealing the spatiotemporal characteristics of XCH4 corresponding to CH4 emissions over Monsoon Asia is discussed in Section 4. Finally, the conclusion is described in Section 5.

GOSAT XCH4 Retrievals
We collected the GOSAT XCH4 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 XCH4 data product using the "quality_flag" tag in the data product documents to get the valid datasets. It has been reported that these XCH4 retrievals have a mean global bias of 1.9 ppb with a standard deviation of 13.4 ppb compared with groundbased observations of the Total Carbon Column Observing Network (TCCON) [32]. Figure 1 shows the density of the collected XCH4 retrieval points from April 2009 to April 2020. As shown in Figure 1, there are many gaps in XCH4 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.

Data for Validation and Analysis
We collected TCCON data and the CH4 simulations from Carbon Tracker-CH4 for the validation of the accuracy of Mapping-XCH4 dataset. Additionally, we collected the surface CH4 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 XCH4 retrievals. We collected the

Data for Validation and Analysis
We collected TCCON data and the CH 4 simulations from Carbon Tracker-CH 4 for the validation of the accuracy of Mapping-XCH 4 dataset. Additionally, we collected the surface CH 4 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 4 retrievals. We collected the available TCCON data released in the 2014 version over 22 sites shown in Figure 1a  CarbonTracker-CH 4 is released by the National Oceanic and Atmospheric Administration (NOAA). CarbonTracker introduces the assimilation of global atmospheric CH 4 observations from surface air samples and tall towers [34,35]. The CarbonTracker-CH 4 , CT2010, produces estimates of global atmospheric CH 4 mole fractions and surface-atmosphere fluxes released from 2000 to 2010. The model-simulated CH 4 concentration is gridded CH 4 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 4 mole fraction profile data from CarbonTracker-CH 4 are converted to XCH 4 (CT-XCH 4 ) 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 2 if averaging kernel smoothing is applied to model simulations or not [39].
To assess the performance of Mapping-XCH 4 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 2 , CH 4 , and N 2 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 4 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].

Mapping XCH 4 Based on Spatiotemporal Geostatistics
We developed an approach to generate global mapping XCH 4 data in a grid of 1 • by 1 • and intervals of 3 days from 2009 to 2020 by applying spatiotemporal geostatistics to XCH 4 retrievals derived from GOSAT observations. Figure 2 shows the main processing steps in the approach, including the analysis of spatiotemporal trend of XCH 4 , modeling of the correlation structure in space and time, and prediction based on space-time kriging method and generating mapping XCH 4 . The Mapping-XCH 4 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 XCH4 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 XCH4 on a global scale. In this study, we only used the XCH 4 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 4 on a global scale.
2.2.1. Modeling the Spatiotemporal Trend and Correlation Structure of XCH 4 XCH 4 can be represented by a variable Z = {Z(s, t)|s ∈ S, t ∈ T} that varies within a spatial domain S and time interval T. Its spatiotemporal variations can be separated into a deterministic trend component m(s, t) and a stochastic residual component R(s, t), as shown in Equation (1).
where w = (2 * π)/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, a 0 is the background XCH 4 at the starting time, and a 1-3 , β 1-4 and γ 1-4 are parameters to be estimated by the least-squares technique [26]. As shown in Equation (2), the deterministic spatial trend in m(s, t), depending on the spatial distribution of CH 4 sources and sinks [46], is modeled by a linear surface function [46]. The deterministic temporal trend in m(s, t), including the interannual XCH 4 increase and the inherent seasonal XCH 4 cycle, is determined by a set of annual harmonic functions with a simple linear function. The stochastic residual component R(s, t) in Equation (1), which represents the local variability not explained by the deterministic trend component, is derived by subtracting the trend components m(s, t) from the XCH 4 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(s 0 , t 0 ) at an unobserved position (s 0 , t 0 ) can be calculated as the linear weighted sum of the XCH 4 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 ∧ γ(h s , h t ) of XCH 4 data at spatial lag h s and temporal lag h t is given by: where Z(s i , t i ) is the observation at a spatial location and in a time (s i , t i ). N(h s , h t ) is the number of data pairs within a distance of (h s , h t ). Once the experimental variogram has been constructed, a spatiotemporal variogram model γ(h, u) 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 γ(h, u) 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): where γ S (h) and γ T (u) 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{sillγ S (h); sillγ T (u)}. sillγ S (h) and sillγ T (u) 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]. 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   , s t can be calculated as the linear weighted sum of the XCH4 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].
where      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 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.
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 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.
as a linear weighted sum of residual data within a kriging neighborhood [22,46] in space and time relative to the prediction location ( ) 0 0 , s t as shown in Equation (5).
where n is the number of observations to be used. i λ is the weighting factor assigned to a known observation ( ) , i i R s t so as to minimize the prediction error variance while maintaining unbiased prediction. The weighting factor is calculated by Equation (6).   Based on the modeled spatiotemporal variogram, space-time Kriging estimates an arbitrary target point at unobserved location (s 0 , t 0 ) from the stochastic residual component as a linear weighted sum of residual data within a kriging neighborhood [22,46] in space and time relative to the prediction location (s 0 , t 0 ) as shown in Equation (5).
where n is the number of observations to be used. λ i is the weighting factor assigned to a known observation R(s i , t i ) so as to minimize the prediction error variance while maintaining unbiased prediction. The weighting factor is calculated by Equation (6).
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) where σ 2 is the variance of kriging prediction, and 1 is the n × 1 unit vector. The deterministic trend component (the temporal trend of XCH 4 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.

Precision Evaluation of the Mapping-XCH 4 Dataset
We evaluated the performance of the mapping XCH 4 dataset using cross-validation and compared it with TCCON measurements and CH 4 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 ∧ Z s j , t j , j = 1, 2, . . . , N and the corresponding original dataset Z s j , t j , j = 1, 2, . . . , N , can be obtained. We selected four evaluation criteria, the correlation coefficient (R) between two datasets, the mean absolute prediction and the percentage of prediction error less than 15 ppb, and the population mean prediction error ME = 1 to evaluate the results [52]. Additionally, we also discussed the performance of the Mapping-XCH 4 dataset in detecting the spatial and timely variations from 2009 to 2020 over Monsoon Asia.  Figure 5, present a significant correlation coefficient (R 2 ) of 0.97 which indicate small prediction biases between the observed XCH 4 and the predicted XCH 4 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 4 tends to be larger than the predicted XCH 4 mostly at high concentrations, while it is smaller than the predicted XCH 4 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. 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 4 emissions in the two areas. statistical results of cross-validation for the five regions are show in Table 2. 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 XCH4 emissions in the two areas.  We can obtain the corresponding kriging variance as indicated in Equation (7) for each geostatistical prediction to evaluate the uncertainty of Mapping-XCH4 dataset. The kriging variance depends on both the density of XCH4 retrievals and the data homogeneity surrounding the prediction position. The denser the observations and more homogeneous the XCH4 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-XCH4 dataset from 2009 to 2020. Figure 6a shows that the uncertainty of Mapping-XCH4 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 XCH4 surrounding the prediction location induced by the effects of atmospheric conditions such as  We can obtain the corresponding kriging variance as indicated in Equation (7) for each geostatistical prediction to evaluate the uncertainty of Mapping-XCH 4 dataset. The kriging variance depends on both the density of XCH 4 retrievals and the data homogeneity surrounding the prediction position. The denser the observations and more homogeneous the XCH 4 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 4 dataset from 2009 to 2020. Figure 6a shows that the uncertainty of Mapping-XCH 4 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 4 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.

Evaluation of the
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 4 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 4 variation and atmospheric conditions in these areas. 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-XCH4 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 XCH4 variation and atmospheric conditions in these areas.

Comparison with TCCON Measurements
XCH4 from TCCON measurements has high accuracy and can be used for validation of the Mapping-XCH4 [54]. In this study, we selected 22 representative TCCON sites where XCH4 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-XCH4 by the comparison of the XCH4 measurements (X-axis) over TCCON sites plotted in different colors and the Mapping-XCH4 (Y-axis) that coincide with TCCON sites in geolocation and time. As seen from Figure 7, the correlation coefficient (CORR), indicating consistency between Mapping-XCH4 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.

Comparison with TCCON Measurements
XCH 4 from TCCON measurements has high accuracy and can be used for validation of the Mapping-XCH 4 [54]. In this study, we selected 22 representative TCCON sites where XCH 4 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 4 by the comparison of the XCH 4 measurements (X-axis) over TCCON sites plotted in different colors and the Mapping-XCH 4 (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 4 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) , 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-XCH4 dataset as described above.   A Table A1) , 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 4 dataset as described above. Figure 8 shows the spatial pattern of the annual average XCH 4 derived from the Mapping-XCH 4 dataset from 2009 to 2020. The local enlarged annual mean Mapping-XCH 4 for Monsoon Asia from 2010 to 2020 are shown in Figure A1 in the Appendix A. As shown in Figures 8 and A1 in the Appendix A, the highest XCH 4 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 4 generally varies from 1820 ppb to 1860 ppb from the northwest to southeast, and its difference is up to 60 ppb in winter. The XCH4 around tropical areas in South America and Africa shows higher XCH4 concentrations, which is likely related to wetland emissions. This may also be associated with the uncertainty of Mapping-XCH4 in the tropical area of South America (Figure 6b) caused by much smaller number of available GOSAT XCH4 retrievals (Figure 1a). Figure 9 shows the latitudinal and temporal variation in XCH4, which is computed by monthly averaged XCH4 within a 1° latitudinal band using the Mapping-XCH4 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-XCH4 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 XCH4 is. Mapping-XCH4 is higher in the Northern Hemisphere than in the Southern Hemisphere due to much higher anthropogenic CH4 emissions in the Northern Hemisphere. High XCH4 values are also present in the equator at approximately 15°N with high temperature and many wetlands, which largely impacts the enhancement of XCH4 surface emissions. The XCH 4 around tropical areas in South America and Africa shows higher XCH 4 concentrations, which is likely related to wetland emissions. This may also be associated with the uncertainty of Mapping-XCH 4 in the tropical area of South America (Figure 6b) caused by much smaller number of available GOSAT XCH 4 retrievals (Figure 1a). Figure 9 shows the latitudinal and temporal variation in XCH 4 , which is computed by monthly averaged XCH 4 within a 1 • latitudinal band using the Mapping-XCH 4 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 4 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 4 is. Mapping-XCH 4 is higher in the Northern Hemisphere than in the Southern Hemisphere due to much higher anthropogenic CH 4 emissions in the Northern Hemisphere. High XCH 4 values are also present in the equator at approximately 15 • N with high temperature and many wetlands, which largely impacts the enhancement of XCH 4 surface emissions. Figure 10 shows the timely variation in globally average XCH 4 from 2009 to 2020. The globally average XCH 4 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 4 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 4 is partly driven by the abundance of hydroxyl radicals (OH), the most important sink of CH 4 , in the atmosphere. As a result, a high abundance of OH in the summer drives down XCH 4 [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 4 emissions. Wetlands, such as swamps and lakes, have significant CH 4 emissions due to anaerobic environments [80,81].

Global XCH 4 Variations Revealed by the Mapping-XCH 4 Dataset XCH 4 in 1 • × 1 • grids from 2009 to 2020
from Figure 9 the Mapping-XCH4 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 XCH4 is. Mapping-XCH4 is higher in the Northern Hemisphere than in the Southern Hemisphere due to much higher anthropogenic CH4 emissions in the Northern Hemisphere. High XCH4 values are also present in the equator at approximately 15°N with high temperature and many wetlands, which largely impacts the enhancement of XCH4 surface emissions.  Figure 10 shows the timely variation in globally average XCH4 from 2009 to 2020. The globally average XCH4 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 XCH4 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 XCH4 is partly driven by the abundance of hydroxyl radicals (OH), the most important sink of CH4, in the atmosphere. As a result, a high abundance of OH in the summer drives down XCH4 [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 CH4 emissions. Wetlands, such as swamps and lakes, have significant CH4 emissions due to anaerobic environments [80,81]. The mean annual increase in XCH4 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]. Figure 11 presents the mean XCH4 in 2010 derived from the Mapping-XCH4 dataset and the model simulating XCH4 by CarbonTacker (CT-XCH4) in the same year, and Figure  12 shows their difference. The Mapping-XCH4 generally presents the similar spatial pattern to the CT-XCH4 (Figure 11). The Mapping-XCH4 in Southeast Asia, Northern India and the rainforest in Africa and Southern America, however, are 20-57 ppb larger than the CT-XCH4 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-XCH4 in central India, the central area of Northern Asia, North America, and South America is 10-37 ppb smaller than the modeled XCH4. The mean annual increase in XCH 4 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]. Figure 11 presents the mean XCH 4 in 2010 derived from the Mapping-XCH 4 dataset and the model simulating XCH 4 by CarbonTacker (CT-XCH 4 ) in the same year, and Figure 12 shows their difference. The Mapping-XCH 4 generally presents the similar spatial pattern to the CT-XCH 4 ( Figure 11). The Mapping-XCH 4 in Southeast Asia, Northern India and the rainforest in Africa and Southern America, however, are 20-57 ppb larger than the CT-XCH 4 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 4   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-XCH4 as the Mapping-XCH4 shows the lowest uncertainties in this area shown in Figure 6b. As a result, the large differences between Mapping-XCH4 and CT-XCH4 mostly arise in the wetland and the paddy fields, which imply that the CT-XCH4 likely have large uncertainty there due to the flaw of prior surface emission in simulating [35].

Spatiotemporal Characteristics of Mapping XCH4 Corresponding to the Surface Emissions
The variations in CH4 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-XCH4 and the annual surface CH4 emissions in 2010 from the EDGAR emission inventory was shown by Liu et al. [26]. Here, we primarily applied the Mapping-XCH4 dataset from 2009 to 2020 generated by the approach above to reveal the spatiotemporal characteristics of XCH4 corresponding to the surface emissions at the regional scale. Figure 13 demonstrates the spatial pattern of annual mean XCH4 in 2010 derived from the Mapping-XCH4 and CT-XCH4, respectively, the anthropogenic emissions from EDGAR and the uncertainty (Kstd) of Mapping-XCH4 in the same year. Comparing Mapping-XCH4 with CT-XCH4, as shown in Figure 13a,b), we can find that Mapping-XCH4 reveals a finer spatial variation than CT-XCH4. The spatial pattern of Mapping-XCH4,  The large difference in the middle of Africa is likely due to the uncertainty of m simulation in this tropical area rather than the Mapping-XCH4 as the Mapping-X shows the lowest uncertainties in this area shown in Figure 6b. As a result, the large ferences between Mapping-XCH4 and CT-XCH4 mostly arise in the wetland and paddy fields, which imply that the CT-XCH4 likely have large uncertainty there due to flaw of prior surface emission in simulating [35].

Spatiotemporal Characteristics of Mapping XCH4 Corresponding to the Surface Emissio
The variations in CH4 concentration primarily depend on surface emissions partly on sinks from chemical reactions in the atmosphere and soil, in addition to atm pheric transport by global and local circulations with seasonality. Surface emission clude anthropogenic emissions from agriculture and waste, fossil fuel production and natural emissions from wetlands, geological lakes, termites, biomass burning, etc. emissions are mostly from agriculture and waste (34%), wetlands (30%), and fossil production and use (19%) [82]. The significantly high correlation between Mapping-X and the annual surface CH4 emissions in 2010 from the EDGAR emission inventory shown by Liu et al. [26]. Here, we primarily applied the Mapping-XCH4 dataset from to 2020 generated by the approach above to reveal the spatiotemporal characteristi XCH4 corresponding to the surface emissions at the regional scale. Figure 13 demonstrates the spatial pattern of annual mean XCH4 in 2010 der from the Mapping-XCH4 and CT-XCH4, respectively, the anthropogenic emissions f EDGAR and the uncertainty (Kstd) of Mapping-XCH4 in the same year. Comparing M ping-XCH4 with CT-XCH4, as shown in Figure 13a,b), we can find that Mapping-X reveals a finer spatial variation than CT-XCH4. The spatial pattern of Mapping-X 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 4 as the Mapping-XCH 4 shows the lowest uncertainties in this area shown in Figure 6b. As a result, the large differences between Mapping-XCH 4 and CT-XCH 4 mostly arise in the wetland and the paddy fields, which imply that the CT-XCH 4 likely have large uncertainty there due to the flaw of prior surface emission in simulating [35].

Spatiotemporal Characteristics of Mapping XCH 4 Corresponding to the Surface Emissions
The variations in CH 4 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 4 and the annual surface CH 4 emissions in 2010 from the EDGAR emission inventory was shown by Liu et al. [26]. Here, we primarily applied the Mapping-XCH 4 dataset from 2009 to 2020 generated by the approach above to reveal the spatiotemporal characteristics of XCH 4 corresponding to the surface emissions at the regional scale. Figure 13 demonstrates the spatial pattern of annual mean XCH 4 in 2010 derived from the Mapping-XCH 4 and CT-XCH 4 , respectively, the anthropogenic emissions from EDGAR and the uncertainty (Kstd) of Mapping-XCH 4 in the same year. Comparing Mapping-XCH 4 with CT-XCH 4, as shown in Figure 13a,b, we can find that Mapping-XCH 4 reveals a finer spatial variation than CT-XCH 4 . The spatial pattern of Mapping-XCH 4 , moreover, is in better agreement with the CH 4 emissions of EDGAR ( Figure 14) than that of CT-XCH 4 . In particular, the Sichuan Basin in China, Bangladesh, and northern India with large emissions corresponds to high values in Mapping-XCH 4 . These large CH 4 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 2 columns indicated the large magnitude of industrial emissions over this region also [84]. The Mapping-XCH 4 does not respond to surface emissions only in the southern tip of the peninsula of South Asia where the large uncertainty of Mapping-XCH 4 (Figure 13d) due to a smaller number of available XCH 4 retrievals (Figure 1a). These results indicate that the Mapping-XCH 4 can reveal and capture the local CH 4 emissions in space and time as the Mapping-XCH 4 is generated by instantaneous satellite observations that respond to surface emissions [85,86].
Remote Sens. 2022, 14, x FOR PEER REVIEW 17 of 25 moreover, is in better agreement with the CH4 emissions of EDGAR ( Figure 14) than that of CT-XCH4. In particular, the Sichuan Basin in China, Bangladesh, and northern India with large emissions corresponds to high values in Mapping-XCH4. These large CH4 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 NO2 columns indicated the large magnitude of industrial emissions over this region also [84]. The Mapping-XCH4 does not respond to surface emissions only in the southern tip of the peninsula of South Asia where the large uncertainty of Mapping-XCH4 (Figure 13d) due to a smaller number of available XCH4 retrievals (Figure 1a). These results indicate that the Mapping-XCH4 can reveal and capture the local CH4 emissions in space and time as the Mapping-XCH4 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 4 in Monsoon Asia in 2010. Similarly, a significant linear relationship was also observed by the Mapping-XCH 4 data, with a coefficient of determination (R 2 ) of 0.66 and a p value less than 0.01. These results indicated that the variation in CH 4 concentration was greatly affected by the surface emissions.

Temporal Variations of XCH 4 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 4 revealed by the Mapping-XCH 4 datasets. Figure 15 presents a yearly increasing trend from 2009 to 2020 and the difference of annual enhancements for these sampling regions. Additionally, Figure 14 shows the correlation between surface emissions and Mapping-XCH4 in Monsoon Asia in 2010. Similarly, a significant linear relationship was also observed by the Mapping-XCH4 data, with a coefficient of determination (R 2 ) of 0.66 and a p value less than 0.01. These results indicated that the variation in CH4 concentration was greatly affected by the surface emissions.

Temporal Variations of XCH4 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 XCH4 revealed by the Mapping-XCH4 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 XCH4 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 XCH4, 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 XCH4 in Shanxi Province (S5) with  The relationship between CH4 surface emissions and the average XCH4 in 2018 in 1° by 1° grids derived from the Mapping-XCH4 dataset. The black-red line is derived from linear regression of average XCH4 data (Y-axis) and the base 10 logarithm of surface CH4 emissions (X-axis), which shows a significant linear relationship with R 2 equal to 0.66 (p value < 0.01).
Additionally, Figure 14 shows the correlation between surface emissions and Mapping-XCH4 in Monsoon Asia in 2010. Similarly, a significant linear relationship was also observed by the Mapping-XCH4 data, with a coefficient of determination (R 2 ) of 0.66 and a p value less than 0.01. These results indicated that the variation in CH4 concentration was greatly affected by the surface emissions.

Temporal Variations of XCH4 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 XCH4 revealed by the Mapping-XCH4 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 XCH4 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 XCH4, 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 XCH4 in Shanxi Province (S5) with 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 4 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 4 , 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 4 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 4 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 2 concentration in Shanxi Province during this lockdown period showed the similar change to XCH 4 [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 4 dataset for detecting the variations and the evidence of CH 4 above, find (1) the mean annual increase of XCH 4 and seasonal variation derived from the Mapping-XCH 4 dataset is agreement with that from in-situ observational network; (2) the Mapping-XCH 4 dataset could explain the variable of anthropogenic emissions with significant correlation between the Mapping-XCH 4 and EDGAR emission inventory (R 2 = 0.66) over Monsoon Asia; (3) the Mapping-XCH 4 can be used to detect spatiotemporal characteristics regionally, and the evidence of CH 4 variations induced by the special events, such as the decrease of XCH 4 during January-April in 2020 in China caused by reduction of anthropogenic emission due to the lock down of COVID-19 pandemic.

Conclusions
In this study, we proposed a data-driven approach based on a spatiotemporal geostatistical model to generate a global land Mapping-XCH 4 dataset in 1 • × 1 • grids and intervals of 3 days from 2009 to 2020 using XCH 4 retrievals derived by GOSAT observations. This Mapping-XCH 4 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 4 dataset by primarily investigating the spatial pattern and timely variation of XCH 4 and comparing it with the model simulations. The results show that the timely variations in XCH 4 characterized by the Mapping-XCH 4 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 4 revealed by the Mapping-XCH 4 dataset correspond to the distribution of surface CH 4 emissions from paddy fields, wetlands, and anthropogenic emissions. Moreover, the Mpping-XCH 4 dataset presents a finer spatial pattern than the model simulations. These results demonstrated that the Mapping-XCH 4 dataset could help us to investigate the spatiotemporal patterns of XCH 4 at global and regional scales. The Mapping-XCH 4 dataset has the advantage of being spatiotemporally continuous compared with the original XCH 4 retrievals with many gaps in space and time. The Mapping-XCH 4 dataset with long-term in grid, such as shown in Figure A1 in Appendix A, facilitates the detection of the driving factors of CH 4 variations combined with other satellite observation data, such as ecological parameters of vegetation and land cover related to CH 4 natural emissions, and inventory data of agriculture, wetlands, and anthropogenic emissions. This could be like the application of global mapping XCO 2 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 4 retrievals prevents us from resolving this problem. Therefore, the precision of Mapping-XCH 4 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 4 retrievals derived from satellite observations. Mapping-XCH 4 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 4 will hopefully be reduced along with increasing XCH 4 data available from multiple satellites to improve the geostatistical model.