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

Abstract

Column-averaged dry air mole fraction of atmospheric CO₂ (XCO₂), 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₂ concentrations which are related to carbon uptake and emissions. In order to construct long-term spatio-temporal continuous XCO₂ 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₂. The approach integrated XCO₂ from different sensors considering differences in vertical sensitivity, overpass time, the field of view, repeat cycle and measurement precision. We produced globally mapped XCO₂ (GM-XCO₂) 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₂ 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² of 0.97 from cross-validation. GM-XCO₂ 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₂ or other similar datasets observed from multiple satellite sensors. The resulting GM-XCO₂ product may be also used in different carbon cycle research applications with different precision requirements.

1. Introduction

Spatio-temporal variation of atmospheric CO₂ 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 CO₂ concentrations [2] from 277 parts per million (ppm) in 1750 [3] to 407.4 ± 0.1 ppm in 2018 [4]. The growth rate of CO₂ concentrations in the atmosphere is smaller than the rate of CO₂ 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 CO₂ [6], which varies spatially due to non-uniform land-biosphere CO₂ exchange [7]. In addition, there is spatial-temporal variability of atmospheric CO₂ concentrations that can be used to study changes in regional land biosphere net CO₂ fluxes, for example, seasonal cycle amplitude increase [8,9] and regional effects of extreme weather patterns like droughts [10,11]. Atmospheric CO₂ 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 CO₂ and related carbon sources/sinks distribution are still not fully understood [15,16]. A long time series of comparable global CO₂ concentration datasets has the potential to improve our understanding of land-biosphere interactions and our ability to evaluate trends in regional terrestrial CO₂ absorption capacity.

There are several methods to measure atmospheric CO₂ concentrations, including surface measurements and satellite observations and approaches to estimate concentrations from model simulations. A network of surface CO₂ monitoring station observations has been organized into the popular GLOBALVIEW-CO₂ 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 CO₂ using estimated surface fluxes and atmospheric mixing transport in addition to the previously noted sparse validation stations [17]. Satellite observations of atmospheric CO₂ 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 CO₂ concentration observations [20,21]. Leveraging all available XCO₂ datasets to construct a long time series, continuous, and comparable global CO₂ concentrations would be useful.

Satellite-observed column-averaged dry air mole fraction of CO₂ (XCO₂) have been wildly used for carbon cycle studies, including CO₂ 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 XCO₂ 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 XCO₂ 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 XCO₂ dataset starting from 2003 using XCO₂ retrievals from these satellites with careful integration and gap-filling.

Gap-filling satellite-observed XCO₂ has been investigated in several studies from different perspectives [31,32,33]. Geostatistical approaches, especially kriging, were widely used for GOSAT observed XCO₂ Level 3 product production [17,34,35]. Spatial-only geostatistical methods do not take into account the temporal correlation structure of CO₂ data [19], which may provide extra information. In order to make full use of spatio-tempal correlation of atmospheric CO₂, a new spatio-temporal kriging method was developed for the global mapping of XCO₂ [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 spatio-temporal resolution and a long time series of XCO₂ 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 XCO₂ and leverage multiple measurements to improve precision when possible, we developed a precision-weighted spatio-temporal kriging method for gap filling of integrated XCO₂ from multiple satellite observations. Datasets used in this study and data preprocessing are described in Section 2. XCO₂ integration methods and global mapping can be found in Section 3. Results of global mapped XCO₂ and its validation are shown in Section 4 and data quality considerations are discussed in Section 5. Finally, conclusions are presented in Section 6.

2. Dataset

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

2.1. XCO₂ from Multi-Satellite Observations

Atmospheric CO₂ concentrations used here were the released Level 2 products that contain the full-physics retrievals of column-averaged CO₂ in units of dry air mole fraction (XCO₂). Satellite observations of XCO₂ used here are from ENVISAT/SCIAMACHY, GOSAT, and OCO-2, which span from 2003 to 2016. XCO₂ 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. XCO₂ from GOSAT is produced by the algorithm of the Atmospheric CO₂ 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. XCO₂ 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 XCO₂ product by screening criteria specified by the corresponding user guides. Specifications of the three satellites that observed XCO₂ are shown in

Table 1. Specifications of multiple satellites that observed XCO_2.

Attributes\Satellites	ENVISAT/SCIAMACHY	GOSAT	OCO-2
Period of selected data	January 2003–March 2012	June 2009–May 2016	September 2014–December 2016
Repeat cycle (days)	35	3	16
Field of view (km)	30 × 60	Diameter of 10.5	2.25×1.25
Overpass local time	10:00	13:00	13:36
version	BESD v02.01.01	ACOS v7.3	OCO2 r9
Profile layers number	10	20	20
Criteria of data screening	XCO2_quality_flag=0;	XCO2_quality_flag=0; gain=H; land_fraction>90; warn_level<10	XCO2_quality_flag=0; gain=H; land_fraction>90
Name referred hereafter	SCI-XCO2	GOS-XCO2	OCO-XCO2
Reference	[37]	[38]	[39]

. These satellites follow different orbits and have different gaps as shown in Figure 1. XCO₂ retrievals from each satellite also have different data precision for different sensors, observation conditions, and retrieval methods. Therefore, we use different XCO₂ precisions in time and space in the XCO₂ integration and mapping with the spatio-temporal kriging method.

2.2. The Total Carbon Column Observing Network

The total carbon column observing network (TCCON) are upward-looking terrestrial Fourier transform spectrometers established for measuring atmospheric XCO₂ and other trace gases from the surface [40]. The instruments have high accuracy with approximately 0.25% error in XCO₂ 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 XCO₂ products.

2.3. XCO₂ from Model Simulation

CarbonTracker (CT) is a modeling system that assimilates global atmospheric CO₂ observations from the ground, tall tower, and aircraft, coupled with an atmospheric transport model for simulating global distributions of atmospheric CO₂ and tracking CO₂ sources and sinks. Model-simulated atmospheric CO₂ 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 CO₂ profiles in a grid of 2° × 3° (latitude x longitude) in 3-hour intervals as the common profile to align the a priori CO₂ profiles and averaging kernels of multi-satellite observations. Second, we adopted diurnal patterns of CO₂ concentrations in different pressure layers from CarbonTracker to unify the CO₂ from satellite observations to 13:00 local time. XCO₂ from CarbonTracker (CT-XCO₂) was calculated from the model CO₂ profile data with 25 layers at the local time 13:00 by using a pressure weighting average method [43]. Third, we used the CT-XCO₂ for comparing with our new globally mapped XCO₂ product.

3. Method

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

3.1. Preprocessing

3.1.1. 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 CO₂ profiles to account for senor sensitivities in XCO₂ 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 XCO₂ 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 XCO₂ retrievals [37,45]. In this study, we introduce a common a priori XCO₂ profile from CT to integrate XCO₂ retrievals from ENVISAT/ SCIAMACHY, GOSAT, and OCO-2 by using Equation (1).

$${\mathrm{XCO}}_{2_{\mathrm{adj},\mathrm{t}}}={\mathrm{XCO}}_{2_{\mathrm{ret},\mathrm{t}}}+{\mathit{h}}^T·\left(\mathit{I}-\mathit{A}\right)·\left({\mathrm{X}}_{\mathrm{M},\mathrm{t}}-{\mathrm{X}}_{\mathrm{a},\mathrm{t}}\right)$$

(1)

where ${\mathrm{XCO}}_{2_{\mathrm{adj},\mathrm{t}}}$ is the adjusted ${\mathrm{XCO}}_{2_{\mathrm{ret}}}$ at observation time t, ${\mathrm{XCO}}_{2_{\mathrm{ret},\mathrm{t}}}$ is the original XCO₂ retrievals from satellites, h is a pressure-weighting vector, A is column-averaging kernels in the XCO₂ retrieval algorithm, I is an identity matrix, ${\mathrm{X}}_{\mathrm{M},\mathrm{t}}$ is a set of common a priori CO₂ profiles from CT, and ${\mathrm{X}}_{\mathrm{a},\mathrm{t}}$ is a set of a priori CO₂ profiles used in XCO₂ original satellite-specific retrieval algorithms. A priori CO₂ profiles of each satellite, as shown in Table 1, were interpolated into the same 25 pressure layers of the CT model.

3.1.2. 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 CO₂ 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).

$${\mathrm{XCO}}_{2_{\mathrm{con},\mathrm{rt}}}=\frac{{\mathit{h}}^T{\mathrm{X}}_{\mathrm{M},\mathrm{rt}}}{{\mathit{h}}^T{\mathrm{X}}_{\mathrm{M},\mathrm{t}}}·{\mathrm{XCO}}_{2_{\mathrm{adj},\mathrm{t}}}$$

(2)

where ${\mathrm{XCO}}_{2_{\mathrm{con},\mathrm{rt}}}$ is the converted ${\mathrm{XCO}}_{2_{\mathrm{adj}}}$ at the reference time (rt, 13:00 local time); ${\mathrm{XCO}}_{2_{\mathrm{adj},\mathrm{t}}}$ is the adjusted XCO₂ derived from Equation (1) at satellite overpass time t; $X_{\mathrm{M},\mathrm{rt}}$ and $X_{\mathrm{M},\mathrm{t}}$ are CO₂ profiles from CT at times of rt and t, respectively; and $\mathit{h}$ is the pressure-weighting vector.

Moreover, XCO₂ 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 XCO₂ retrievals using precision weighted averaging of ${\mathrm{XCO}}_{2_{\mathrm{con},\mathrm{rt}}}$ within 30 km by 30 km every 8 days using Equations (3-1) and (3-2). 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 XCO₂ dataset at 30 km resolution every 8 days from 2003 to 2016 (integrated-XCO₂) was generated.

$${\mathrm{XCO}}_{2_{\mathrm{int},\mathrm{rt}}}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{p}}_{\mathrm{i}}·{\mathrm{XCO}}_{2_{\mathrm{con},{\mathrm{rt}}_{\mathrm{i}}}}$$

(3-1)

$${\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{p}}_{\mathrm{i}}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\frac{{\mathsf{\mu }}^2}{{\mathrm{m}}_{\mathrm{i}}^2}=1$$

(3-2)

where ${\mathrm{XCO}}_{2_{\mathrm{int},\mathrm{rt}}}$ is the integrated combination of ${\mathrm{XCO}}_{2_{\mathrm{con}}}$ datasets, n is the number of observations in one unit, ${\mathrm{XCO}}_{2_{\mathrm{con},{\mathrm{rt}}_{\mathrm{i}}}}$ is converted satellite observed XCO₂, $p_i$ is the weighting factor, which is determined by ${\mathrm{m}}_{\mathrm{i}}$ (data precision of ${\mathrm{XCO}}_{2_{\mathrm{con},{\mathrm{rt}}_{\mathrm{i}}}}$). $\mu$ is an arbitrary constant used for normalizing the data precision weighting factor. We adopted the data precision from these satellites observed XCO₂ Level 2 product. The data precision in ACOS-GOSAT and OCO-2 is XCO₂ posterior error, and that in BESD-SCIAMACHY is 1-sigma uncertainty of the retrieved XCO₂.

3.2. Modeling XCO₂ Spatio-Temporal Random Field for use in Kriging

XCO₂ 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 second-order stationary random field represented by the stochastic residual component after removing inter-annual, 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 XCO₂. 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).

$$\mathrm{m}\left(\mathrm{s},\mathrm{t}\right)={\mathrm{a}}_0\left(s\right)+{\mathrm{a}}_t\left(s\right)+{\displaystyle \sum }_{i=1}^4\left({\beta }_i\left(s\right)·\sin \left(i\omega t\right)+{\gamma }_i\left(s\right)·\cos \left(i\omega t\right)\right)+R\left(s,t\right)$$

(4)

where $\omega =2·\mathsf{\pi }/\mathrm{T}$ and T is the period of 46 time-units, t is the time in the time-unit, s represents the latitudinal zone, and ${\mathrm{a}}_0$, ${\beta }_{1-4}$ and ${\gamma }_{1-4}$ are parameters to be estimated. ${\mathrm{a}}_t$ is the cumulative annual increase for each time-unit determined by the Earth System Research Laboratory (ESRL) global annual CO₂ growth rate [4]. The harmonic functions fit the annual cycle, semi-annual oscillation, seasonal variation and monthly variation of XCO₂ [19]. $R\left(s,t\right)$ represents the spatio-temporal residual component of satellite observed integrated XCO₂, which will be used for interpretation.

3.3. 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 ($\mathrm{Z}\left({\mathrm{s}}_0\right)$) at point ${\mathrm{s}}_0$ can be expressed as a weighted sum of the observational data as shown in Equation (5).

$$\mathrm{Z}\left({\mathrm{s}}_0\right)={\displaystyle \sum }_{i=1}^n{\omega }_i·Z\left(s_i\right)$$

(5)

where ${\omega }_i$ is the weighting factor at the observation point $s_i$, n is the total number of observational points to be used, and $Z\left(s_i\right)$ is the observed value at each point $s_i$. 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 ‘${\omega }_i$’.

3.3.1. Conventional Spatio-Temporal Kriging

In spatio-temporal geostatistical analysis of XCO₂, kriging prediction of Z($s_0$,$t_0$) at a point ($s_0$,$t_0$) can be calculated as the linear weighted sum of the XCO₂ 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).

$${\widehat{\mathsf{\gamma }}}_{ST}\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$$

(6)

where ${\widehat{\gamma }}_{ST}\left(h_s,h_t\right)$ is an empirical variogram value at the lag ($h_s,h_t$). $Z\left(s_i,t_i\right)$ is the observational data. $N\left(h_s,h_t\right)$ is the number of data pairs within a distance of ($h_s,h_t$). 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)).

$${\mathsf{\gamma }}_{\mathrm{ST}}\left({\mathrm{h}}_{\mathrm{s}},{ \mathrm{h}}_{\mathrm{t}};{ \mathsf{\theta }}_{\mathrm{s}},{ \mathsf{\theta }}_{\mathrm{t}}, \mathsf{\kappa },{ \mathrm{N}}_{\mathrm{ST}}\right)={\mathsf{\gamma }}_S\left({\mathrm{h}}_{\mathrm{s}};{\mathsf{\theta }}_{\mathrm{s}}\right)+{\mathsf{\gamma }}_T\left({\mathrm{h}}_{\mathrm{t}};{\mathsf{\theta }}_{\mathrm{t}}\right)-\kappa ·{\mathsf{\gamma }}_S\left({\mathrm{h}}_{\mathrm{s}};{\mathsf{\theta }}_{\mathrm{s}}\right)·{\mathsf{\gamma }}_T\left({\mathrm{h}}_{\mathrm{t}};{\mathsf{\theta }}_{\mathrm{t}}\right)+{ \mathrm{N}}_{\mathrm{ST}}·{\mathsf{\gamma }}_0\left({\mathrm{h}}_{\mathrm{s}};{\mathrm{h}}_{\mathrm{t}}\right)$$

(7)

We selected the exponential model for the marginal variogram model ${\mathsf{\gamma }}_S\left({\mathrm{h}}_{\mathrm{s}}\right)$ and ${\mathsf{\gamma }}_T\left({\mathrm{h}}_{\mathrm{t}}\right)$ as shown in Equation (8).

$$\mathsf{\gamma }\left(\mathrm{h};\mathsf{\theta }=\left[\mathrm{C}, \mathrm{a}\right]\right)=\left\{\begin{array}{c}0,\hspace{1em}\hspace{1em} h=0\\ C·\left(1-e^{-\frac{h}{a}}\right),\hspace{1em}\hspace{1em} h\ne 0\end{array}\right.$$

(8)

where ${\mathsf{\gamma }}_S\left({\mathrm{h}}_{\mathrm{s}};{\mathsf{\theta }}_{\mathrm{s}}\right)$ and ${\mathsf{\gamma }}_T\left({\mathrm{h}}_{\mathrm{t}};{\mathsf{\theta }}_{\mathrm{t}}\right)$ are the marginal spatial and temporal variograms, [${\mathsf{\theta }}_{\mathrm{s}}$,${ \mathsf{\theta }}_{\mathrm{t}}$,$\kappa ,{ \mathrm{N}}_{\mathrm{ST}}$] = [${\mathrm{C}}_{\mathrm{s}}$, ${\mathrm{a}}_{\mathrm{s}}$, ${\mathrm{C}}_{\mathrm{t}}$, ${\mathrm{a}}_{\mathrm{t}}$, $\kappa$, ${\mathrm{N}}_{\mathrm{ST}}$] are parameters to be estimated, and ${\mathrm{N}}_{\mathrm{ST}}$, C, and a are all greater than or equal to zero. In the exponential model, a, C and ${\mathrm{N}}_{\mathrm{ST}}$ represent the influence range, partial sill, and nugget effects, respectively.

As a result, an arbitrary target point ($\hat{\mathrm{Z}}\left({\mathrm{s}}_0,{\mathrm{t}}_0\right)$) to be estimated by using the spatio–temporal kriging method can be expressed as Equation (9).

$$\hat{\mathrm{Z}}\left({\mathrm{s}}_0,{\mathrm{t}}_0\right)={{\displaystyle \sum }}_{i=1}^n{\omega }_i\left({\mathrm{s}}_0,{\mathrm{t}}_0\right)·Z\left(s_i,t_i\right) with {{\displaystyle \sum }}_{i=1}^n{\omega }_i\left({\mathrm{s}}_0,{\mathrm{t}}_0\right)=1$$

(9)

where ${\omega }_i\left({\mathrm{s}}_0,{\mathrm{t}}_0\right)$ is the weight assigned to a known observation $Z\left(s_i,t_i\right)$ 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

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

(10)

where $\mathsf{\Gamma }\left(\mathrm{i},\mathrm{j}\right)=\mathsf{\gamma }\left(\left|{\mathrm{s}}_i-{\mathrm{s}}_j\left|,\right|{\mathrm{t}}_i-{\mathrm{t}}_j\right|\right)$, ${\mathsf{{\rm Y}}}_0\left(\mathrm{i}, 1\right)$=$\mathsf{\gamma }\left(\left|{\mathrm{s}}_i-{\mathrm{s}}_0\left|,\right|{\mathrm{t}}_i-{\mathrm{t}}_0\right|\right)$, and 1 is the n × 1 unit vector.

3.3.2. Optimization of Spatio-Temporal Correlation Structure

In a conventional spatio-temporal geostatistical analysis, all data pairs were used for spatio-temporal correlation structure using equal weight. XCO₂ 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 spatio-temporal correlation structure with precision weighting factor as shown in Equation (11).

$$\mathsf{\delta }={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}{\mathsf{\lambda }}_{\mathrm{i}}·{\left({\mathsf{\gamma }}_{\mathrm{ST}}-{\hat{\mathsf{\gamma }}}_{\mathrm{ST}}\right)}^2={\displaystyle \sum }_{i=1}^m{\mathsf{\lambda }}_{\mathrm{i}}·{\left({\mathsf{\gamma }}_s\left(h_s\right)+{\mathsf{\gamma }}_T\left(h_t\right)-\kappa ·{\mathsf{\gamma }}_s\left(h_s\right) ·{\mathsf{\gamma }}_T\left(h_t\right)+N_{ST}-\hat{\mathsf{\gamma }}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)\right)}^2$$

(11-1)

$${\displaystyle \sum }_{i=1}^n{\mathsf{\lambda }}_{\mathrm{i}}={\displaystyle \sum }_{i=1}^n\frac{{\mu }_1{}^2}{m_i^2}=1$$

(11-2)

where ${\mathsf{\lambda }}_{\mathrm{i}}$ represents different weighting factors for data with different precisions. ${\mathsf{\gamma }}_{\mathrm{ST}}$ is the spatio-temporal variogram model shown in Equation (7). ${\hat{\mathsf{\gamma }}}_{\mathrm{ST}}$ is the empirical variogram shown in Equation (6). ${\mathrm{m}}_{\mathrm{i}}$ and ${\mu }_1$ represent data precision and a normalization term.

The gradient descent method was used to calculate the optimal parameters. The partial derivative parameters (${\mathsf{\delta }}_{{\mathrm{C}}_{\mathrm{s}}^{\ast }}^{‘}$,${ \mathsf{\delta }}_{{\mathrm{a}}_{\mathrm{s}}^{\ast }}^{‘}$, ${\mathsf{\delta }}_{{\mathrm{C}}_{\mathrm{t}}^{\ast }}^{‘}$,${ \mathsf{\delta }}_{{\mathrm{a}}_{\mathrm{t}}^{\ast }}^{‘}$, ${\mathsf{\delta }}_{{\mathsf{\kappa }}^{\ast }}^{‘}$, ${\mathsf{\delta }}_{{\mathrm{N}}_{\mathrm{ST}}^{\ast }}^{‘}$) that need to be estimated are shown in Equations (12-1)–(12-6).

$${\mathsf{\delta }}_{{\mathrm{C}}_{\mathrm{s}}^{\ast }}^{‘}=2·{\displaystyle \sum }_{i=1}^m\left\{{\mathsf{\lambda }}_{\mathrm{i}}·\left({\mathsf{\gamma }}_{\mathrm{ST}}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)-\hat{\mathsf{\gamma }}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)\right)\right\}·\left(1-e^{-\frac{\left|\right|h_s\left|\right|}{a_s^{\ast }}}\right)$$

(12-1)

$${\mathsf{\delta }}_{{\mathrm{a}}_{\mathrm{s}}^{\ast }}^{‘}=2·{\displaystyle \sum }_{i=1}^m\left\{{\mathsf{\lambda }}_{\mathrm{i}}·\left({\mathsf{\gamma }}_{\mathrm{ST}}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)-\hat{\mathsf{\gamma }}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)\right)\right\}·\left(-{\mathrm{C}}_{\mathrm{s}}^{\ast }·e^{-\frac{\left|\right|h_s\left|\right|}{a_s^{\ast }}}\right)·(||h_s||·a_s^{\ast }{}^{-2})$$

(12-2)

$${\mathsf{\delta }}_{{\mathrm{C}}_{\mathrm{t}}^{\ast }}^{‘}=2·{\displaystyle \sum }_{i=1}^m\left\{{\mathsf{\lambda }}_{\mathrm{i}}·\left({\mathsf{\gamma }}_{\mathrm{ST}}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)-\hat{\mathsf{\gamma }}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)\right)\right\}·\left(1-e^{-\frac{\left|\right|h_t\left|\right|}{a_t^{\ast }}}\right)$$

(12-3)

$${\mathsf{\delta }}_{{\mathrm{a}}_{\mathrm{t}}^{\ast }}^{‘}=2·{\displaystyle \sum }_{i=1}^m\left\{{\mathsf{\lambda }}_{\mathrm{i}}·\left({\mathsf{\gamma }}_{\mathrm{ST}}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)-\hat{\mathsf{\gamma }}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)\right)\right\}·\left(-{\mathrm{C}}_{\mathrm{t}}^{\ast }·e^{-\frac{\left|\right|h_t\left|\right|}{a_t^{\ast }}}\right)·(||h_t||·a_t^{\ast }{}^{-2})$$

(12-4)

$${\mathsf{\delta }}_{{\mathsf{\kappa }}^{\ast }}^{‘}=2·{\displaystyle \sum }_{i=1}^m\left\{{\mathsf{\lambda }}_{\mathrm{i}}·\left({\mathsf{\gamma }}_{\mathrm{ST}}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)-\hat{\mathsf{\gamma }}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)\right)\right\}·\left({\mathsf{\gamma }}_{\mathrm{S}}\left({\mathrm{h}}_{\mathrm{s}};{ \mathsf{\theta }}_{\mathrm{s}}^{\ast }\right)·{\mathsf{\gamma }}_{\mathrm{T}}\left({\mathrm{h}}_{\mathrm{t}};{ \mathsf{\theta }}_{\mathrm{t}}^{\ast }\right)\right)$$

(12-5)

$${\mathsf{\delta }}_{{\mathrm{N}}_{\mathrm{ST}}^{\ast }}^{‘}=2·{\displaystyle \sum }_{i=1}^m\left\{{\mathsf{\lambda }}_{\mathrm{i}}·\left({\mathsf{\gamma }}_{\mathrm{ST}}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)-\hat{\mathsf{\gamma }}\left({\mathrm{h}}_s,{\mathrm{h}}_t\right)\right)\right\}$$

(12-6)

where ${\mathrm{C}}_{\mathrm{s}}^{\ast },{\mathrm{a}}_{\mathrm{s}}^{\ast },{\mathrm{C}}_{\mathrm{t}}^{\ast },{\mathrm{a}}_{\mathrm{t}}^{\ast },{\mathsf{\kappa }}^{\ast }{ \mathrm{and} \mathrm{N}}_{\mathrm{ST}}^{\ast }$ are parameters waiting to be optimized. Parameters, ${\mathsf{\lambda }}_{\mathrm{i}}$,${ \mathsf{\gamma }}_{\mathrm{ST}}$, $\mathsf{\gamma }$ and $\hat{\mathsf{\gamma }}$, 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 $\mathsf{\delta }$. Initial parameters ${\mathsf{\beta }}_0=\left({\mathrm{C}}_{\mathrm{s}0},{\mathrm{a}}_{\mathrm{s}0},{\mathrm{C}}_{\mathrm{t}0},{\mathrm{a}}_{\mathrm{t}0},{\mathsf{\kappa }}_0,{\mathrm{N}}_{\mathrm{ST}0}\right)$ 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).

$${\mathsf{\beta }}_i={\mathsf{\beta }}_{i-1}+\mathsf{\alpha }·{\delta }_{i-1}^{\prime }$$

(13)

where $\mathsf{\beta }$ represents (${\mathrm{C}}_{\mathrm{s}}^{\ast },{\mathrm{a}}_{\mathrm{s}}^{\ast },{\mathrm{C}}_{\mathrm{t}}^{\ast },{\mathrm{a}}_{\mathrm{t}}^{\ast },{\mathsf{\kappa }}^{\ast },{ \mathrm{N}}_{\mathrm{ST}}^{\ast })$. $\mathsf{\alpha }$ is the learning rate, usually in the range of $({10}^{-4} ~ {10}^{-2})$. ${\delta }^{\prime }$ 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.

3.3.3. Integrating XCO₂ Using Variable Data Precision

We estimated XCO₂ in unobserved points using observational data and weighting factors from spatio-temporal correlation structure and data precision as shown in Equation (14).

$$\mathrm{Z}\left({\mathrm{s}}_0,{\mathrm{t}}_0\right)={\displaystyle \sum }_{i=1}^n{\omega }_{1i}·{\omega }_{2i}·Z\left(s_i,t_i\right)$$

(14)

where ${\omega }_{i1}$ is the weighting factor from the spatio-temporal correlation structure and ${\omega }_{i2}$ is from data precision.$Z\left(s_i,t_i\right)$ is observed XCO₂ used for $\mathrm{Z}\left({\mathrm{s}}_0,{\mathrm{t}}_0\right)$ estimation. These two weighting factors were calculated by Equations (15-1) and (15-2).

$${\displaystyle \sum }_{i=1}^s{\omega }_{1i}·\gamma \left(h_{1i}\right)+\epsilon =\gamma \left(h_{10}\right)$$

$${\displaystyle \sum }_{i=1}^s{\omega }_{1i}·\gamma \left(h_{2i}\right)+\epsilon =\gamma \left(h_{20}\right)$$

$$\dots \dots$$

$${\displaystyle \sum }_{i=1}^s{\omega }_{1i}·\gamma \left(h_{si}\right)+\epsilon =\gamma \left(h_{s0}\right)$$

(15-1)

$${\displaystyle \sum }_{i=1}^s{\omega }_{1i}·{\omega }_{2i}=1$$

(15-2)

where ${\omega }_{1i}$ and ${\omega }_{2i}$ are weighting factors in Equation (14). Equation (15-2) was used to control unbiased estimation. $\gamma \left(h_{1i}\right)\dots$ $\gamma \left(h_{si}\right)$ and $\gamma \left(h_{10}\right)$ … $\gamma \left(h_{s0}\right)$ are spatio-temporal variograms for observations used for estimation of $\mathrm{Z}\left({\mathrm{s}}_0,{\mathrm{t}}_0\right)$. $\epsilon$ is the polynomial residuals. ${\omega }_{2i}$ can be achieved by Equation (16).

$${\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{s}}{\omega }_{2i}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{s}}\frac{{\mathsf{\mu }}_2{}^2}{{\mathrm{m}}_{\mathrm{i}}^2}=\mathrm{s}$$

(16)

where ${\mathrm{m}}_{\mathrm{i}}$, ${\mathsf{\mu }}_2$ and $\mathrm{s}$ are data precision, arbitrary constant and the number of used observational data.

We applied this method for global mapping of XCO₂ (GM-XCO₂), which provides a long time series of spatio-temporal continuous XCO₂ dataset for global carbon cycle research.

3.3.4. Uncertainty and Precision of Mapped XCO₂

We calculated data uncertainty and precision of the mapped XCO₂ by using precision-weighted kriging. Uncertainty of predicted data from the spatio-temporal kriging method is shown in Equation (17)

$${\mathsf{\sigma }}_0{}^2={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{s}}{\mathrm{w}}_{1\mathrm{i}}·\mathsf{\gamma }\left({\mathrm{h}}_{\mathrm{i}0}\right)+\epsilon$$

(17)

where ${\mathsf{\sigma }}_0{}^2$ presents the prediction uncertainty, ${\mathrm{w}}_{1\mathrm{i}}$, $\mathsf{\gamma }\left({\mathrm{h}}_{\mathrm{i}0}\right)$ and $\epsilon$ are the estimated weighting factor, spatio-temporal variograms, and polynomial residuals, respectively as shown in Equation (15).

In addition, the precision of GM-XCO₂ integrated from the observed data is shown in Equation (18), according to the error transfer equation.

$${\mathsf{\epsilon }}_0=\mathrm{sqrt}({\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{s}}{\mathsf{\epsilon }}_{\mathrm{i}}{}^2·{\mathrm{w}}_{1\mathrm{i}}{}^2·{\mathrm{w}}_{2\mathrm{i}}{}^2)$$

(18)

As a result, these two sources of uncertainty can be used for data screening in GM-XCO₂ application.

3.4. Validation of Global Mapped XCO₂

We validated the mapped GM-XCO₂ 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 XCO₂ dataset (${\mathrm{XCO}}_{2_{\mathrm{int}}}$) for interpolation to GM-XCO₂. The cross-validation was conducted by repeatedly (100 times) reserving 5% of the ${\mathrm{XCO}}_{2_{\mathrm{int}}}$ data for validation. Predicted GM-XCO₂ was compared with the 5% reserved ${\mathrm{XCO}}_{2_{\mathrm{int}}}$ data in those corresponding spatio-temporal locations. We selected three statistical parameters for results evaluation: (1) the coefficient of determination (r²), (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-XCO₂ with the XCO₂ measured from TCCON (TCCON-XCO₂) in 12 sites with more than 5 years of observation. TCCON-XCO₂ 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-XCO₂ and XCO₂ from the CT model simulation in the spatio-temporal change over the global area.

4. Results

4.1. Integrated-XCO₂ from Three Satellites

Here we presented latitudinal and temporal variability of the integrated XCO₂ (top panel) and the difference between the integrated product and the original retrieval XCO₂ values (bottom panel) from SCIAMACHY, GOSAT, and OCO-2 (Figure 4). The top panel shows XCO₂ increased more than 30 ppm from 2003 to 2016 over most latitudes. Especially high XCO₂ values occurred during the start of 2016 in the northern hemisphere (dark red). The bottom panel shows XCO₂ adjustments, integrated XCO₂ (${\mathrm{XCO}}_{2_{\mathrm{int}}}$) minus original XCO₂ (${\mathrm{XCO}}_{2_{\mathrm{ret}}}$) 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 XCO₂ averaging kernel (Appendix A: Figure A1) that was adopted for XCO₂ adjustment in Equation (1). In addition, the adjustments decreased sharply in 2012 and 2016, with stability improvements of the newer satellite sensors.

4.2. Globally-Mapped XCO₂

4.2.1. Latitudinal and Temporal Variability of Globally Mapped XCO₂

The latitudinal and temporal variability of globally mapped XCO₂ (GM-XCO₂) 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-XCO₂ also shows the yearly increase and seasonal variation from 2003 to 2016. GM-XCO₂ is higher in the northern hemisphere compared to the southern hemisphere, and the annual maximum GM-XCO₂ lasts longer in mid-latitudes than high/low latitudes.

In the middle panel, kriging standard deviations (square root of ${\mathsf{\sigma }}_0{}^2$ in Equation 17) of the predicted GM-XCO₂ 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-XCO₂ precision is calculated from integrated XCO₂ (${\mathrm{XCO}}_{2_{\mathrm{int}}}$) 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-XCO₂ precision is better in mid-latitudes.

4.2.2. 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) of GM-XCO₂ from precision weighted spatio-temporal kriging is improved for most of the predictions, especially for the SCIAMACHY and GOSAT overlapping periods (June 2009–March 2012) or GOSAT and OCO-2 overlaps (September 2014–May 2016). The precision improved by 0.3–0.5 ppm over latitudes with overlapping observations. GM-XCO₂ (Figure 6a) is enhanced by 0.1–0.2 ppm for the summer of 2009 to 2011 and reduced by 0.2–0.3 for the summer of 2014 to 2016. Differences in the uncertainty of prediction (Figure 6b) are small for both kriging results (conventional and precision-weighted) because they are based on the same observations (${\mathrm{XCO}}_{2_{\mathrm{int}}}$), which results in a similar spatio-temporal correlation structure as shown in Figure 3.

4.2.3. Spatial Distribution of GM-XCO₂

GM-XCO₂ 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-XCO₂ during different seasons in 2003, 2008, 2013, and 2015 (Figure 7). Seasonal patterns of GM-XCO₂ were similar with an annual increase of approximately around 2.0 ppm. In spring, high GM-XCO₂ appeared in northern Canada, the North China Plain, and the Arabian Peninsula. In summer, extremely low GM-XCO₂ occurred in mid-high latitudes in the northern hemisphere. In autumn, the north-south hemispheric GM-XCO₂ gradient relaxes. In winter, high XCO₂ returns over the North China Plain and Central Africa.

In addition, we give out the global spatial distribution of mean XCO₂ in 2016 shown in Figure 8. We can find that the main high XCO₂ 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.

4.3. GM-XCO₂ Validation

4.3.1. Evaluation Using Cross-Validation

Results of cross-validation using the precision-weighted spatio-temporal kriging method are shown in Figure 9. Predicted XCO₂ (GM-XCO₂) agreed well with integrated XCO₂ (${\mathrm{XCO}}_{2_{\mathrm{int}}}$) with a high R² of 0.97, showing the interpolation method retains much of the input signal. Total RMSE (root mean square error) between predicted GM-XCO₂ and integrated XCO₂ (${\mathrm{XCO}}_{2_{\mathrm{int}}}$) 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 XCO₂ 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 XCO₂ precision, respectively. Results from cross-validation suggest that this mapping method is effective and precise in gap-filling of multi-satellites observed XCO₂, which could also be affected by original data precision.

4.3.2. Validation of GM-XCO₂ with TCCON Measurements

High accuracy and continuous time series of XCO₂ from TCCON measurements were used for validation of GM-XCO₂ 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-XCO₂. ${\mathrm{XCO}}_{2_{\mathrm{int}}}$ and GM-XCO₂ data within 500 km of each TCCON site was used for this comparison (Figure 10). GM-XCO₂ retained information from satellite observations, which captured the annual increase and seasonal variation of XCO₂ well. The temporal variation of GM-XCO₂ is consistent with TCCON XCO₂. Comparison statistics between GM-XCO₂ and TCCON XCO₂ are shown in

Table 2. Comparison statistics between GM-XCO₂ and TCCON XCO₂. Bias is calculated using GM-XCO₂ minus TCCON XCO₂ for each coincident data pair and averaged for each site.

Sites	Location (Latitude, Longitude)	Coincident Data Pairs	Averaged Bias (ppm)	Averaged Absolute Bias (ppm)	Standard Deviation (ppm)
Bialystok	(53.23°N, 23.02°E)	249	−0.19	0.73	0.92
Bremen	(53.10°N, 8.85°E)	260	0.21	0.97	1.25
Karlsruhe	(49.10°N, 8.44°E)	226	0.51	0.90	0.98
Orleans	(47.97°N, 2.11°E)	232	0.34	0.71	0.85
Garmisch	(47.48°N, 11.06°E)	361	0.62	1.05	1.16
Park Falls	(45.94°N, 90.27°W)	499	0.00	0.74	0.96
Lamont	(36.60°N, 97.49°W)	381	−0.45	0.77	0.92
Tsukuba	(36.05°N, 140.12°E)	210	0.73	1.70	1.89
JPL/Caltech	(34.20°N, 118.18°W)	243	−1.06	1.19	0.97
Saga	(33.24°N, 130.29°E)	204	−0.33	0.76	0.91
Darwin	(12.43°S, 130.89°E)	434	−0.47	0.89	1.00
Wollongong	(34.41°S, 150.88°E)	341	0.09	0.58	0.75
Overall	-	303	0.01	0.92	1.05

. If we assume TCCON measurements are accurate, the accuracy of GM-XCO₂ performs well with the total averaged bias of 0.01 ppm across 303 data pairs. The precision of GM-XCO₂ 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.

4.4. Comparison between GM-XCO₂ and CarbonTracker Simulated XCO₂

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

4.4.1. Comparison with Latitudinal and Temporal Variability of CT-XCO₂

The latitudinal and temporal variability of the difference between GM-XCO₂ (Figure 5) and CT-XCO₂ (Appendix A: Figure A2) and statistical summary are shown in Figure 11. The mean difference is 1.53 ± 0.80 ppm. XCO₂ 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.

XCO₂ from low/high latitudes, especially for data before March 2012, shows the largest differences between GM-XCO₂ and CT-XCO₂. The potential reasons for this are: (1) Sparse satellite observation in high/low latitudes limited the accuracy of GM-XCO₂; (2) Limited precision of XCO₂ from SCIAMCHY contributed to large differences compared to CT-XCO₂; (3) Lacking or limited surface measurements in low/high latitudes constrained the CT model and affected the simulation of accurate XCO₂. GM-XCO₂ provided spatial-temporal continuous XCO₂ based on satellite observations that are different from the model simulation.

4.4.2. Temporal Variability of GM-XCO₂ and CT-XCO₂ in Mid-Latitudes

GM-XCO₂ may provide insights into local carbon uptake and emission. We present the temporal variability of ${\mathrm{XCO}}_{2_{\mathrm{int}}}$, GM-XCO_2, and XCO₂ from CarbonTracker (CT-XCO₂) from 2003 to 2016 over mid-latitudes of North America and Eastern Asia (Figure 12). Mean GM-XCO₂ kept the original XCO₂ information from the satellites observation by kriging the ${\mathrm{XCO}}_{2_{\mathrm{int}}}$. However, GM-XCO₂ presented sharp changes in XCO₂ during the peak and minimum of XCO₂ in each year, which is not as smooth as that of CT-XCO₂. In addition, the temporal variability of GM-XCO₂ was more consistent with CT-XCO₂ 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-XCO₂ might provide more information for areas with limited surface measurements, such as China.

5. Discussion

There are some advantages of GM-XCO₂ used in global carbon cycle studies. GM-XCO₂ provided spatial-temporal continuous XCO₂ data from 2003 to 2016 which filled the gaps in observation as shown in Figure 1. It improved the data precision compared with XCO₂ from conventional spatio-temporal kriging, especially for data in the satellites overlapping period. As a result, GM-XCO₂ could help us understand the temporal and spatial changes in the global distribution of CO₂ [17]. Because GM-XCO₂ is from instantaneous satellite observations, it could capture detailed and abnormal XCO₂ 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 XCO₂ precision and interpolation uncertainty for each XCO₂. Users can select different temporal and spatial segments of GM-XCO₂ data for their specific application. As it was discussed that space-based XCO₂ 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 CO₂ flux uncertainties significantly [56]. For anthropogenic CO₂ monitoring, a precision requirement of GM-XCO₂ must be better than 0.7 ppm [57].

However, there are also some challenges that need to be discussed. Prediction of GM-XCO₂ in regions with little observations could have low precision and high interpolation uncertainty. GM-XCO₂ 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-XCO₂ 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].

6. Conclusions

In this study, we developed a precision-weighted spatio-temporal kriging interpolation method of multiple satellite observed XCO₂. It not only used the spatio-temporal variability of XCO₂, but also the precision of each observation for gap filling. The spatio-temporal correlation structure was optimized and the weighting of XCO₂ 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 XCO₂ but also for other data observed by multiply satellites with different precision.

We also produced a spatial-temporal continuous XCO₂ product, which provides the data precision and interpolation uncertainty of each grid. It is spatio-temporal continuous XCO₂ from satellite observations which could capture the detailed change of XCO₂. It could be an available dataset for global carbon cycle studies like biosphere-atmosphere interaction.