## 2 Materials and Methods
### Methods
#### 2.2.2 Spatiotemporal Kriging  
It is a common mapping method for spatiotemporal discrete points to a continuum. The process of calculating weight coefficients and estimating weighted averages is basically the same as that of spatial Kriging interpolation, which extends the computational dimension from space to space-time. Suppose \(\left\{Z(s,t),(s,t)\in D\times T\subseteq\mathfrak{R}^{d+1}\right\}\) is a space-time random field, \(s\ =\ (s_{1},s_{2},\cdots,s_{d})\) (\(d\leq 3\), generally) is the space coordinate and \(t\leq T\) is the time coordinate. The interpolation formula of the spatiotemporal Kriging is as follows [38] (p. 39):  
\[\widehat{Z(s,t)}_{0}\ =\ \sum_{i\ =\ 1}^{n}\lambda_{i}Z(s,t)_{i} \tag{2}\]  
where \(\widehat{Z(s,t)}_{0}\) is the predicted result of point \(s\) at \(t\) time. \(Z(s,t)_{i}\) denotes the observation; \(n\) is the number of observations; and \(\lambda_{i}\) is the weight coefficient which constitutes the optimal set of coefficients that satisfy the minimum difference between the estimation and the observation at the point \((s,t)_{0}\), i.e.,: \(\underset{\lambda_{i}}{\min}Var\left(\widehat{Z(s,t)}_{0}-Z(s,t)_{0}\right)\). Simultaneously, satisfying the conditions for unbiased estimation (Equation (3)) [38] (pp. 31-32):  
\[E\left(\widehat{Z(s,t)}_{0}-Z(s,t)_{0}\right)\ =\ 0 \tag{3}\]  
For the residuals \(\hat{R}(s,t)\) obtained in this study, we need to test its stationarity before constructing the spatiotemporal Kriging model. After passing the test, we first calculated the spatiotemporal empirical semi-variogram function and then fit the spatiotemporal theoretical variogram function model. In this step, we chose the spatiotemporal product-sum model as the theoretical variogram model to fit the spatiotemporal variational structure of the spatiotemporal geographic data [18; 36], which is obtained by transforming the known pure spatial covariance and pure temporal covariance functions by addition and multiplication, mixing, and integrals. The covariance function and variational function are as follows (Equations (4) and (5), [38], pp. 25-27):  
\[C(h_{s},h_{t})\ =\ k_{1}C_{s}(h_{s})C_{t}(h_{t})+k_{2}C_{s}(h_{s})+k_{3}C_{t}(h_{t}) \tag{4}\]  
\[\gamma(h_{s},h_{t})\ =\ (k_{1}C_{t}(0)+k_{2})\gamma_{s}(h_{s})+(k_{1}C_{s}(0)+k_{3 })\gamma_{t}(h_{t})-k_{1}\gamma_{s}(h_{s})\gamma_{t}(h_{t}) \tag{5}\]  
In Equation (4), \(\gamma(h_{s},h_{t}),\gamma_{s}(h_{s}),\gamma_{t}(h_{t})\) are the corresponding spatiotemporal, spatial, and temporal variogram functions, respectively, while \(C(0,0),C_{s}(0),C_{t}(0)\) are the corresponding sill values, respectively. Generally, the maximum value approaching stability is taken when computing the spatiotemporal variogram function for the sample. This study implemented spatiotemporal Kriging interpolation based on the \"gstat\" package in the R programming language.  
Finally, the constructed spatiotemporal model was used to interpolate the residuals. The trend term was added to the residual term so as to obtain the final XCO2 estimates. For a more comprehensive understanding, please refer to the previous works [39; 40].  
Figure 2: CNN structure.