Spatio-Temporal Mapping of Multi-Satellite Observed Column Atmospheric CO 2 Using Precision-Weighted Kriging Method

: Column-averaged dry air mole fraction of atmospheric CO 2 (XCO 2 ), obtained by multiple satellite observations since 2003 such as ENVISAT/SCIAMACHY, GOSAT, and OCO-2 satellite, is valuable for understanding the spatio-temporal variations of atmospheric CO 2 concentrations which are related to carbon uptake and emissions. In order to construct long-term spatio-temporal continuous XCO 2 from multiple satellites with different temporal and spatial periods of observations, we developed a precision-weighted spatio-temporal kriging method for integrating and mapping multi-satellite observed XCO 2 . The approach integrated XCO 2 from different sensors considering differences in vertical sensitivity, overpass time, the field of view, repeat cycle and measurement precision. We produced globally mapped XCO 2 (GM-XCO 2 ) with spatial/temporal resolution of 1 × 1 degree every eight days from 2003 to 2016 with corresponding data precision and interpolation uncertainty in each grid. The predicted GM-XCO 2 precision improved in most grids compared with conventional spatio-temporal kriging results, especially during the satellites overlapping period (0.3–0.5 ppm). The method showed good reliability with R 2 of 0.97 from cross-validation. GM-XCO 2 showed good accuracy with a standard deviation of bias from total carbon column observing network (TCCON) measurements of 1.05 ppm. This method has potential applications for integrating and mapping XCO 2 or other similar datasets observed from multiple satellite sensors. The resulting GM-XCO 2 product may be also used in different carbon cycle research applications with different precision requirements.


Introduction
Spatio-temporal variation of atmospheric CO2 concentration reflects the balance between anthropogenic carbon emissions and terrestrial and oceanic carbon uptake or emissions [1]. Increased fossil fuel emissions after the start of the Industrial Revolution contribute to the continuous growth of atmospheric CO2 concentrations [2] from 277 parts per million (ppm) in 1750 [3] to 407.4±0.1 ppm in 2018 [4]. The growth rate of CO2 concentrations in the atmosphere is smaller than the rate of CO2 emitted by human activities because nearly 45% of the emissions are absorbed by oceans and the terrestrial biosphere each year [5]. The seasonal variations of terrestrial carbon uptake and emission contribute most to the seasonal cycle in atmospheric CO2 [6], which varies spatially due to nonuniform land-biosphere CO2 exchange [7]. In addition, there is spatial-temporal variability of atmospheric CO2 concentrations that can be used to study changes in regional land biosphere net CO2 fluxes, for example, seasonal cycle amplitude increase [8,9] and regional effects of extreme weather patterns like droughts [10,11]. Atmospheric CO2 concentrations have also been used for carbon flux or to estimate carbon uptake/emission changes using atmospheric inversion models [12,13] and a data-driven method [14]. The spatio-temporal variability of CO2 and related carbon sources/sinks distribution are still not fully understood [15,16]. A long time series of comparable global CO2 concentration datasets has the potential to improve our understanding of land-biosphere interactions and our ability to evaluate trends in regional terrestrial CO2 absorption capacity.
There are several methods to measure atmospheric CO2 concentrations, including surface measurements and satellite observations and approaches to estimate concentrations from model simulations. A network of surface CO2 monitoring station observations has been organized into the popular GLOBALVIEW-CO2 product and provides in situ measurements but is limited by station sparseness and the inherent spatial inhomogeneity of the surface atmosphere. Model simulations can provide continuous maps of CO2 using estimated surface fluxes and atmospheric mixing transport in addition to the previously noted sparse validation stations [17]. Satellite observations of atmospheric CO2 have the advantage of global coverage and high measurement density and can complement the surface network to advance our understandings of the carbon cycle and its changes [18,19]. With the development of remote sensing technology, there are several satellites used for atmospheric CO2 concentration observations [20,21]. Leveraging all available XCO2 datasets to construct a long time series, continuous, and comparable global CO2 concentrations would be useful.
Satellite-observed column-averaged dry air mole fraction of CO2 (XCO2) have been wildly used for carbon cycle studies, including CO2 enhancement detection induced by anthropogenic emissions [22,23], constraining model simulations of carbon fluxes [24,25], investigating carbon cycle responses to weather extremes [11,26], and improving understanding of vegetation uptake [27]. Satellites measuring XCO2 are the SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY (SCIAMACHY) onboard the Environmental Satellite (ENVISAT) [28], Greenhouse Gases Observing Satellite (GOSAT) [29], and Orbiting Carbon Observatory-2 (OCO-2) [30]. They observe XCO2 but with different spatial/temporal resolution, prior vertical profile estimates, local overpass time, data precision, and observing gaps. As a result, there is an opportunity to generate a long time series of global XCO2 dataset starting from 2003 using XCO2 retrievals from these satellites with careful integration and gap-filling.
Gap-filling satellite-observed XCO2 has been investigated in several studies from different perspectives [31][32][33]. Geostatistical approaches, especially kriging, were widely used for GOSAT observed XCO2 Level 3 product production [17,34,35]. Spatial-only geostatistical methods do not take into account the temporal correlation structure of CO2 data [19], which may provide extra information. In order to make full use of spatio-tempal correlation of atmospheric CO2, a new spatiotemporal kriging method was developed for the global mapping of XCO2 [19,36]. Because these methods were previously used for observations from a single satellite, measurement error could be assumed to be uniform and not interfere with the kriging approach. In order to produce high spatiotemporal resolution and a long time series of XCO2 from multiple satellite observations, the precision of different datasets should be considered in this geostatistical method.
In this study, to create the longest possible time series of XCO2 and leverage multiple measurements to improve precision when possible, we developed a precision-weighted spatiotemporal kriging method for gap filling of integrated XCO2 from multiple satellite observations. Datasets used in this study and data preprocessing are described in Section 2. XCO2 integration methods and global mapping can be found in Section 3. Results of global mapped XCO2 and its validation are shown in Section 4 and data quality considerations are discussed in Section 5. Finally, conclusions are presented in Section 6.

Dataset
In this study, we collected atmospheric CO2 concentration datasets from multiple satelliteobservations to produce a long time series of spatio-temporal continuous XCO2 from satellite observations, which was then evaluated by XCO2 from surface measurements and model simulations.

XCO2 from Multi-Satellite Observations
Atmospheric CO2 concentrations used here were the released Level 2 products that contain the full-physics retrievals of column-averaged CO2 in units of dry air mole fraction (XCO2). Satellite observations of XCO2 used here are from ENVISAT/SCIAMACHY, GOSAT, and OCO-2, which span from 2003 to 2016. XCO2 from SCIAMACHY onboard the European ENVISAT are obtained by the full-physics based Bremen Optimal Estimation-DOAS (BESD) algorithm (v02.01.01) [37], which span from January 2003 to March 2012 with a spatial/temporal resolution of 30x60 km every 35 days. XCO2 from GOSAT is produced by the algorithm of the Atmospheric CO2 Observations from Space (ACOS) team (v7.3 lite) [38], which span from June 2009 to May 2016 with spatial/temporal resolution of a diameter of 10.5 km every three days. XCO2 from OCO-2 is produced by the ACOS team (r9 lite) [39], which span from September 2014 to December 2016 with spatial/temporal resolution of 2.25 × 1.25 km every 16 days. Data quality was maximized by filtering XCO2 product by screening criteria specified by the corresponding user guides. Specifications of the three satellites that observed XCO2 are shown in Table 1. These satellites follow different orbits and have different gaps as shown in Figure 1. XCO2 retrievals from each satellite also have different data precision for different sensors, observation conditions, and retrieval methods. Therefore, we use different XCO2 precisions in time and space in the XCO2 integration and mapping with the spatio-temporal kriging method.

The Total Carbon Column Observing Network
The total carbon column observing network (TCCON) are upward-looking terrestrial Fourier transform spectrometers established for measuring atmospheric XCO2 and other trace gases from the surface [40]. The instruments have high accuracy with approximately 0.25% error in XCO2 retrievals, which has been extensively used for validation of satellite observations [21,40,41]. In this study, we selected 12 TCCON sites within the mapping area as shown in Figure 1, with at least five years of coincidental measurements in the period from January 2003 to December 2016 for validating combined XCO2 products.

XCO2 from Model Simulation
CarbonTracker (CT) is a modeling system that assimilates global atmospheric CO2 observations from the ground, tall tower, and aircraft, coupled with an atmospheric transport model for simulating global distributions of atmospheric CO2 and tracking CO2 sources and sinks. Model-simulated atmospheric CO2 from CarbonTracker 2017 [42] was used in three steps of the integration and mapping method, primarily to normalize differences in altitude sensitivity and overpass time. First, we used CarbonTracker CO2 profiles in a grid of 2° × 3° (latitude x longitude) in 3-hour intervals as the common profile to align the a priori CO2 profiles and averaging kernels of multi-satellite observations. Second, we adopted diurnal patterns of CO2 concentrations in different pressure layers from CarbonTracker to unify the CO2 from satellite observations to 13:00 local time. XCO2 from CarbonTracker (CT-XCO2) was calculated from the model CO2 profile data with 25 layers at the local time 13:00 by using a pressure weighting average method [43]. Third, we used the CT-XCO2 for comparing with our new globally mapped XCO2 product.

Method
In this study, we developed a method for spatio-temporal integration and mapping of multisatellite observed XCO2 considering variable data precision for producing globally mapped continuous XCO2 with spatial/temporal resolution of 1° × 1° every eight days from 2003 to 2016. We present the flow chart of global mapped XCO2 production and the precision weighted spatiotemporal kriging method in Figure 2. First, we adjusted a priori vertical CO2 profiles and averaging kernels of multiple satellite-observed XCO2 products to a common profile. Second, we corrected to a common local time and regularizing spatio-temporal scales of XCO2 from multiple satellites. Third, we used a modeled continuous XCO2 spatio-temporal random field for interpolation. Fourth, we developed a precision-weighted spatio-temporal kriging method for producing global maps of XCO2. Finally, we validated the new global mapped XCO2. We present details on developing the precisionweighted spatio-temporal kriging method including: 1) conventional spatio-temporal kriging method; 2) optimization of spatio-temporal correlation structure; 3) XCO2 prediction through integrating data precision, and; 4) uncertainty and precision of mapped XCO2.

Adjustment of a Priori Vertical Profiles and Averaging Kernels
Each satellite has different measurement sensitivity at different altitudes through the atmospheric column and, therefore, they use different averaging kernels based on a priori assumptions of vertical CO2 profiles to account for senor sensitivities in XCO2 retrieval algorithms. The corrections are based on prior vertical profile layers as shown in Table 1. The a priori profiles from different satellite retrievals should be adjusted to a common profile when comparing XCO2 from different instruments. Additionally, the smoothing effect of the retrievals should be considered by applying the averaging kernels [44] to reduce the effects from different instruments on XCO2 retrievals [37,45]. In this study, we introduce a common a priori XCO2 profile from CT to integrate XCO2 retrievals from ENVISAT/ SCIAMACHY, GOSAT, and OCO-2 by using Equation (1).
where XCO , is the adjusted XCO at observation time t, XCO , is the original XCO2 retrievals from satellites, h is a pressure-weighting vector, A is column-averaging kernels in the XCO2 retrieval algorithm, I is an identity matrix, X , is a set of common a priori CO2 profiles from CT, and X , is a set of a priori CO2 profiles used in XCO2 original satellite-specific retrieval algorithms. A priori CO2 profiles of each satellite, as shown in Table 1, were interpolated into the same 25 pressure layers of the CT model.

Unification of Observing Time and Spatio-Temporal Scales
The three satellites have different local overpass times: SCIAMACHY at 10:00, GOSAT at 13:00, and OCO-2 at 13:36 (Table 1). In order to reduce the effect of atmospheric CO2 concentrations diurnal variation [46,47], we introduce a correction coefficient to normalize the satellite observations local time to 13:00 based on diurnal variation of CT model simulations using Equation (2).
Where XCO , is the converted XCO at the reference time (rt, 13:00 local time); XCO , is the adjusted XCO2 derived from Equation (1) at satellite overpass time t; , and , are CO2 profiles from CT at times of rt and t, respectively; and is the pressure-weighting vector. Moreover, XCO2 is affected by the different fields of view and observing dates as shown in Table  1. In order to reduce this effect, we integrate spatial and temporal scales of XCO2 retrievals using precision weighted averaging of XCO , within 30 km by 30 km every 8 days using Equations (3a) and (3b). This unification also reduces computational complexity and preserves local spatiotemporal patterns. A temporal resolution of 8 days is also well-suited for biosphere-atmosphere interaction analysis with other 8-day resolution datasets like vegetation indices from the Moderate Resolution Imaging Spectroradiometer (MODIS). An integrated XCO2 dataset at 30 km resolution every 8 days from 2003 to 2016 (integrated-XCO2) was generated.
where XCO , is the integrated combination of XCO datasets, n is the number of observations in one unit, XCO , is converted satellite observed XCO2, is the weighting factor, which is determined by m (data precision of XCO , ). is an arbitrary constant used for normalizing the data precision weighting factor. We adopted the data precision from these satellites observed XCO2 Level 2 product. The data precision in ACOS-GOSAT and OCO-2 is XCO2 posterior error, and that in BESD-SCIAMACHY is 1-sigma uncertainty of the retrieved XCO2.

Modeling XCO2 Spatio-Temporal Random Field for use in Kriging
XCO2 increases from year-to-year varies by latitude and has a significant seasonal cycle in most locations [48,49]. In order to interpret spatio-temporal geostatistics, we need to construct a secondorder stationary random field represented by the stochastic residual component after removing interannual, latitudinal, and seasonal trends mentioned above, termed the deterministic mean. In this study, we adopted a fitting method for decomposing [50] the deterministic spatiotemporal mean and stochastic residual component of latitude-zonal XCO2. The fitting method is a combination of a linear function to fit the long-term increase and annual periodic function as shown in Equation (4).
where = 2 • π/T and T is the period of 46 time-units, t is the time in the time-unit, s represents the latitudinal zone, and a , and are parameters to be estimated. a is the cumulative annual increase for each time-unit determined by the Earth System Research Laboratory (ESRL) global annual CO2 growth rate [4]. The harmonic functions fit the annual cycle, semi-annual oscillation, seasonal variation and monthly variation of XCO2 [19]. ( , ) represents the spatio-temporal residual component of satellite observed integrated XCO2, which will be used for interpretation.

Precision Weighted Spatio-Temporal Kriging
In the ordinary kriging method, a predicted value at an arbitrary target point is estimated by considering the statistical properties of a set of observed data. As a result, the predicted value (Z(s )) at point s can be expressed as a weighted sum of the observational data as shown in Equation (5).
where is the weighting factor at the observation point , n is the total number of observational points to be used, and ( ) is the observed value at each point . The following subsections describe the precision-weighted spatio-temporal kriging method. Different kriging models were developed by adjusting the number of points 'n' and the weighting factor ' '.

Conventional Spatio-Temporal Kriging
In spatio-temporal geostatistical analysis of XCO2, kriging prediction of Z( , ) at a point ( , ) can be calculated as the linear weighted sum of the XCO2 values that minimizes the mean squared prediction error [19]. Weights of observations used for interpolation are determined by the geometry of observations and the spatio-temporal correlation structure of the data. Spatial and temporal information would be used for variogram modeling of the correlation structure [19,31] as shown in Equation (6).
where (ℎ , ℎ ) is an empirical variogram value at the lag (ℎ , ℎ ). ( , ) is the observational data. (ℎ , ℎ ) is the number of data pairs within a distance of (ℎ , ℎ ). Once the empirical variogram has been constructed, we need to select a spatio-temporal variogram model to fit it. As shown in Zeng et al. [19,31], the spatio-temporal variogram model adopted here (Equation (7)) is a combination of the product-sum model [51,52] and an extra global nugget model to capture the nugget effect [53] (last term in Equation (7)).
We selected the exponential model for the marginal variogram model γ (h ) and γ (h ) as shown in Equation (8).
where γ (h ; θ ) and γ (h ; θ ) are the marginal spatial and temporal variograms, [θ , θ , , N ] = [C , a , C , a , , N ] are parameters to be estimated, and N , C, and a are all greater than or equal to zero. In the exponential model, a, C and N represent the influence range, partial sill, and nugget effects, respectively. As a result, an arbitrary target point (Z (s , t )) to be estimated by using the spatio-temporal kriging method can be expressed as Equation (9).
where (s , t ) is the weight assigned to a known observation ( , ) so as to minimize the prediction error variance while maintaining an unbiased prediction. The prediction error variance, which is a measurement of prediction uncertainty, is given by where Γ(i, j) = γ(|s − s |, |t − t |), Υ (i, 1)= γ(|s − s |, |t − t |), and 1 is the n × 1 unit vector.

Optimization of Spatio-Temporal Correlation Structure
In a conventional spatio-temporal geostatistical analysis, all data pairs were used for spatiotemporal correlation structure using equal weight. XCO2 observations from different satellites/sensors, observing conditions and inversion methods have different data precision. So, the varying data precision should be considered in building the loss function for optimizing spatiotemporal correlation structure with precision weighting factor as shown in Equation (11).
where λ represents different weighting factors for data with different precisions. γ is the spatiotemporal variogram model shown in Equation (7). γ is the empirical variogram shown in Equation (6). m and represent data precision and a normalization term. The gradient descent method was used to calculate the optimal parameters. The partial derivative parameters (δ * ' , δ * ' , δ * ' , δ * ' , δ * ' , δ * ' ) that need to be estimated are shown in Equations where C * , a * , C * , a * , κ * and N * are parameters waiting to be optimized. Parameters, λ , γ , γ and γ , are the weighting factors (Equation (11-2)), exponential models (Equation (7) and Equation (8)) and empirical variogram values (Equation (6)), respectively. We optimized the structure through minimizing δ . Initial parameters β = (C , a , C , a , κ , N ) were obtained by using a least-squares approximation in the conventional spatio-temporal kriging method. Then, parameters were determined by a learning rate and partial derivative as shown in Equation (13).
where β represents ( C * , a * , C * , a * , κ * , N * ) . α is the learning rate, usually in the range of (10 ~ 10 ). is the partial derivative parameter. We set the operating condition as a current change of less than 1% of all the changes in this adjustment. As a result, one example of optimized spatio-temporal semi-variogram surface is shown in Figure 3.  55°N). Grey, black, and red points represent spatio-temporal semi-variogram that was calculated from experimental data, fitted models of the conventional and optimized correlation structure.

Integrating XCO2 Using Variable Data Precision
We estimated XCO2 in unobserved points using observational data and weighting factors from spatio-temporal correlation structure and data precision as shown in Equation (14).
where is the weighting factor from the spatio-temporal correlation structure and is from data precision. ( , ) is observed XCO2 used for Z(s , t ) estimation. These two weighting factors were calculated by Equations (15a) and (15b).
can be achieved by Equation (16).
where m , μ and s are data precision, arbitrary constant and the number of used observational data. We applied this method for global mapping of XCO2 (GM-XCO2), which provides a long time series of spatio-temporal continuous XCO2 dataset for global carbon cycle research.

Uncertainty and Precision of Mapped XCO2
We calculated data uncertainty and precision of the mapped XCO2 by using precision-weighted kriging. Uncertainty of predicted data from the spatio-temporal kriging method is shown in Equation where σ presents the prediction uncertainty, w , γ(h ) and are the estimated weighting factor, spatio-temporal variograms, and polynomial residuals, respectively as shown in Equation (15).
In addition, the precision of GM-XCO2 integrated from the observed data is shown in Equation (18), according to the error transfer equation.
As a result, these two sources of uncertainty can be used for data screening in GM-XCO2 application.

Validation of Global Mapped XCO2
We validated the mapped GM-XCO2 product using cross-validation, compared to TCCON measurements and model simulation. Cross-validation has been used in accuracy assessments for spatio-temporal kriging methods [19,54]. In this study, we adopted cross-validation based on the Monte Carlo sampling method used in Zeng et al. [19]. As above, we used the satellite-observed integrated XCO2 dataset (XCO ) for interpolation to GM-XCO2. The cross-validation was conducted by repeatedly (100 times) reserving 5% of the XCO data for validation. Predicted GM-XCO2 was compared with the 5% reserved XCO data in those corresponding spatio-temporal locations. We selected three statistical parameters for results evaluation: 1) the coefficient of determination (r 2 ), 2) the root mean square error (RMSE), and 3) the percentage of estimation bias less than 1 or 2 times of data precision.
We compared GM-XCO2 with the XCO2 measured from TCCON (TCCON-XCO2) in 12 sites with more than 5 years of observation. TCCON-XCO2 was calculated with the pressure-averaged method [43] and data observed at a local time of 11:00 to 15:00. In addition, we did a comparison between GM-XCO2 and XCO2 from the CT model simulation in the spatio-temporal change over the global area.

Integrated-XCO2 from Three Satellites
Here we presented latitudinal and temporal variability of the integrated XCO2 (top panel) and the difference between the integrated product and the original retrieval XCO2 values (bottom panel) from SCIAMACHY, GOSAT, and OCO-2 ( Figure. 4). The top panel shows XCO2 increased more than 30 ppm from 2003 to 2016 over most latitudes. Especially high XCO2 values occurred during the start of 2016 in the northern hemisphere (dark red). The bottom panel shows XCO2 adjustments, integrated XCO2 (XCO ) minus original XCO2 (XCO ) were mainly within −2.0 to 1.0 ppm and include seasonality, with high adjustment values in summer and low adjustment values in winter, except for some scattered grids in high or low latitudes. This could be caused by seasonal changes of the XCO2 averaging kernel (Appendix A: Figure A1) that was adopted for XCO2 adjustment in Equation (1). In addition, the adjustments decreased sharply in 2012 and 2016, with stability improvements of the newer satellite sensors.

Latitudinal and Temporal Variability of Globally Mapped XCO2
The latitudinal and temporal variability of globally mapped XCO2 (GM-XCO2) and its prediction uncertainty and precision are shown in Figure 5. Comparing Figure 5a to gaps in Figure 4a show missing observations have been reasonably filled. GM-XCO2 also shows the yearly increase and seasonal variation from 2003 to 2016. GM-XCO2 is higher in the northern hemisphere compared to the southern hemisphere, and the annual maximum GM-XCO2 lasts longer in mid-latitudes than high/low latitudes.
In the middle panel, kriging standard deviations (square root of σ in Equation 17) of the predicted GM-XCO2 represents the prediction uncertainty, which is determined by the available data around gaps. Higher uncertainty is shown in the tropic and high latitudes because of low numbers of robust observations in these latitudes corresponding to gaps in Figure 4a. The uncertainty also shows seasonal variation for different latitudes. In mid-high latitudes (35-60 °N), the uncertainty is high in winter for the observations affected by snow cover. In mid-low latitudes (extratropics within 35°N/°S), uncertainty is high in summer likely due to cloud contamination during the monsoon season.
In the bottom panel, GM-XCO2 precision is calculated from integrated XCO2 (XCO ) precision. The precision is significantly improved from the middle of 2009 and 2014. It was 1.5 to 2.5 ppm before 2009, 0.5 to 1.5 from 2009 to 2014, and below 0.5 after 2014. That is because of the observing precision improvement of the newer sensors. Generally, the observed GM-XCO2 precision is better in midlatitudes.

Comparison with Conventional Spatio-Temporal Kriging Results
The latitudinal and temporal differences between the results from the precision-weighted (this study) and conventional spatio-temporal kriging methods are shown in Figure 6. Precision ( Figure  6c

Spatial Distribution of GM-XCO2
GM-XCO2 provides important information about the mean spatial distribution and localized anomalies which could relate to local carbon uptake and emission. We present the spatial-temporal distribution of GM-XCO2 during different seasons in 2003, 2008, 2013, and 2015 ( Figure 7). Seasonal patterns of GM-XCO2 were similar with an annual increase of approximately around 2.0 ppm. In spring, high GM-XCO2 appeared in northern Canada, the North China Plain, and the Arabian Peninsula. In summer, extremely low GM-XCO2 occurred in mid-high latitudes in the northern hemisphere. In autumn, the north-south hemispheric GM-XCO2 gradient relaxes. In winter, high XCO2 returns over the North China Plain and Central Africa.  In addition, we give out the global spatial distribution of mean XCO2 in 2016 shown in Figure 8. We can find that the main high XCO2 of 2016 is distributed in Eastern China, Southeast Asian regions, Amazon forest regions, African forest regions, south of the United States, the Arabian Peninsula, and India. These regions could be related to three main conditions, including 1) an industrial development zone inhabited by human beings; 2) a large area of tropical/subtropical rainforest area; 3) large agricultural regions.

Evaluation Using Cross-Validation
Results of cross-validation using the precision-weighted spatio-temporal kriging method are shown in Figure 9. Predicted XCO2 (GM-XCO2) agreed well with integrated XCO2 (XCO ) with a high R 2 of 0.97, showing the interpolation method retains much of the input signal. Total RMSE (root mean square error) between predicted GM-XCO2 and integrated XCO2 (XCO ) is 1.36 ppm, which indicates good stability and the high precision of this method. Most of the predicted bias is within 1.0 ppm, except for some XCO2 data precision which is larger than 1.5 ppm (right panel). Specifically, 61% of the prediction bias is less than 1.0 ppm. Additionally, 70% and 80% of the predicted bias is within 1 and 2 times of XCO2 precision, respectively. Results from cross-validation suggest that this mapping method is effective and precise in gap-filling of multi-satellites observed XCO2, which could also be affected by original data precision.

Validation of GM-XCO2 with TCCON Measurements
High accuracy and continuous time series of XCO2 from TCCON measurements were used for validation of GM-XCO2 derived from satellite observations. In this study, we selected 12 sites with measurements for more than five years from 2003 to 2016 for comparison with GM-XCO2. XCO and GM-XCO2 data within 500 km of each TCCON site was used for this comparison (Figure 10). GM-XCO2 retained information from satellite observations, which captured the annual increase and seasonal variation of XCO2 well. The temporal variation of GM-XCO2 is consistent with TCCON XCO2. Comparison statistics between GM-XCO2 and TCCON XCO2 are shown in Table 2. If we assume TCCON measurements are accurate, the accuracy of GM-XCO2 performs well with the total averaged bias of 0.01 ppm across 303 data pairs. The precision of GM-XCO2 in most TCCON sites (9/12) is within 1.0 ppm with an averaged absolute bias of 0.92 ppm. The mean value of the standard deviation of the bias over 12 sites is 1.05 ppm.

Comparison between GM-XCO2 and CarbonTracker Simulated XCO2
GM-XCO2, the spatial-temporal continuous XCO2 from satellite observations, can provide a detailed distribution of XCO2 over global or regional land areas. In order to explore the advantages and disadvantages of GM-XCO2, we present the latitudinal-temporal change comparison between GM-XCO2 with CT-XCO2 and their local temporal change comparison in the northern hemisphere.

Comparison with Latitudinal and Temporal Variability of CT-XCO2
The latitudinal and temporal variability of the difference between GM-XCO2 ( Figure 5) and CT-XCO2 (Appendix A: Figure A2) and statistical summary are shown in Figure 11. The mean difference is 1.53 ± 0.80 ppm. XCO2 difference in the mid-latitudes showed better consistency with the mean value than that in low/high latitudes. The difference in high latitudes is smaller, especially for data before March 2012. In low latitudes, the difference varied with time, with lower values in summer and higher values in winter. The differences decreased with the satellite observations precision improvement from 2009 to 2016.
XCO2 from low/high latitudes, especially for data before March 2012, shows the largest differences between GM-XCO2 and CT-XCO2. The potential reasons for this are: 1) Sparse satellite observation in high/low latitudes limited the accuracy of GM-XCO2; 2) Limited precision of XCO2 from SCIAMCHY contributed to large differences compared to CT-XCO2; 3) Lacking or limited surface measurements in low/high latitudes constrained the CT model and affected the simulation of

Temporal Variability of GM-XCO2 and CT-XCO2 in Mid-Latitudes
GM-XCO2 may provide insights into local carbon uptake and emission. We present the temporal variability of XCO , GM-XCO2, and XCO2 from CarbonTracker (CT-XCO2) from 2003 to 2016 over mid-latitudes of North America and Eastern Asia ( Figure 12). Mean GM-XCO2 kept the original XCO2 information from the satellites observation by kriging the XCO . However, GM-XCO2 presented sharp changes in XCO2 during the peak and minimum of XCO2 in each year, which is not as smooth as that of CT-XCO2. In addition, the temporal variability of GM-XCO2 was more consistent with CT-XCO2 in Northern America than that of Eastern Asia. That could be attributed to the contribution from more surface measurements to constrain the models. As a result, GM-XCO2 might provide more information for areas with limited surface measurements, such as China.

Discussion
There are some advantages of GM-XCO2 used in global carbon cycle studies. GM-XCO2 provided spatial-temporal continuous XCO2 data from 2003 to 2016 which filled the gaps in observation as shown in Figure 1. It improved the data precision compared with XCO2 from conventional spatiotemporal kriging, especially for data in the satellites overlapping period. As a result, GM-XCO2 could help us understand the temporal and spatial changes in the global distribution of CO2 [17]. Because GM-XCO2 is from instantaneous satellite observations, it could capture detailed and abnormal XCO2 change which could be related to local carbon uptake and emission [10,55]. It could also be an important dataset for biosphere-atmosphere interactions by relating its changes to local biosphere parameter variations [7,11]. In addition, we provide XCO2 precision and interpolation uncertainty for each XCO2. Users can select different temporal and spatial segments of GM-XCO2 data for their specific application. As it was discussed that space-based XCO2 with precision (no bias) of 2.5 ppm could be used for matching ground-based network and precision within 1.0 ppm, it could help reduce inferred CO2 flux uncertainties significantly [56]. For anthropogenic CO2 monitoring, a precision requirement of GM-XCO2 must be better than 0.7 ppm [57].
However, there are also some challenges that need to be discussed. Prediction of GM-XCO2 in regions with little observations could have low precision and high interpolation uncertainty. GM-XCO2 in the tropical rainforest area and the winter in high latitudes should be used carefully as a few observations were limited by observing conditions [58,59]. The interpolation method might not describe the process of atmospheric transport perfectly, especially for regions with a complex climate like the Tibetan area. In addition, GM-XCO2 over a desert like the Sahara with an influence from high brightness reflection and complex dust for satellite observations inversion should also be used carefully [60].

Conclusion
In this study, we developed a precision-weighted spatio-temporal kriging interpolation method of multiple satellite observed XCO2. It not only used the spatio-temporal variability of XCO2, but also the precision of each observation for gap filling. The spatio-temporal correlation structure was optimized and the weighting of XCO2 with high precision was improved. The precision of predictions improved in most of the grids, especially for the satellites overlapping period (0.3 -0.5 ppm). It would be useful for gap-filling of increasing satellite observations not only in XCO2 but also for other data observed by multiply satellites with different precision.
We also produced a spatial-temporal continuous XCO2 product, which provides the data precision and interpolation uncertainty of each grid. It is spatio-temporal continuous XCO2 from satellite observations which could capture the detailed change of XCO2. It could be an available dataset for global carbon cycle studies like biosphere-atmosphere interaction.

Acronyms
Full Names