## Appendix A Methodological Detail
### Kriging and Trend-Surface-Only Prediction for Validation Comparisons  
The cokriging method described in Sections 2.3 and 2.4 is able to leverage dependence between the latent SIF and XCO\({}_{2}\) processes, respectively, \(Y_{1}(\cdot)\) and \(Y_{2}(\cdot)\), for prediction and uncertainty quantification using the noisy data \(\mathbf{Z}_{1}\) and \(\mathbf{Z}_{2}\). Now, through two successive simplifications of cokriging, we describe two special cases for producing predictions and corresponding prediction standard errors. In the first simplification, we assume there is no cross-dependence between \(Y_{1}(\cdot)\) and \(Y_{2}(\cdot)\), and hence we seek to model \(Y_{1}(\cdot)\) usingalone; we refer to this model as a univariate spatial dependence model. Recall that at each data location \(\mathbf{s}_{l}\in D\), we model the elements of \(\mathbf{Z}_{1}\) as  
\[Z_{1l}=Y_{1}(\mathbf{s}_{l})+\epsilon_{1l}=\mu_{1}(\mathbf{s}_{l})+\psi_{1}( \mathbf{s}_{l})+\epsilon_{1l};\quad l=1,\dots,n_{1}, \tag{10}\]  
where \(\{\epsilon_{1l}\}_{l}\) is a set of (non-stationary) independent measurement errors, assumed to have mean zero and variances \(\{\sigma_{e,1l}^{2}\}_{l}\) given as described for Equation (1). Estimates for the large-scale trend-surface, \(\mu_{1}(\cdot)\), are obtained identically to Equation (4) when \(i=1\). The process \(\psi_{1}^{*}(\cdot)=\eta_{1}^{*}(\cdot)+\omega_{1}^{*}(\cdot)\) is a combination of smooth small-scale variation and micro-scale variation _on the standardized scale_ and is the focus of the univariate spatial dependence model, from which the optimal linear predictor known as the kriging predictor is derived. As in the cokriging method, the micro-scale-variation term, \(\omega_{1}^{*}(\cdot)\), is a process uncorrelated at the grid scale with mean zero and variance \(\sigma_{\omega,1}^{2}/\hat{S}_{R,1}^{2}\), where \(\hat{S}_{R,1}\) is given in Equation (6) when \(i=1\). To model small-scale spatial variation in \(\eta_{1}^{*}(\cdot)\) only, the spatial dependence captured in Equation (8) can be simplified to  
\[C_{1}(h)=\sigma_{1}^{2}M(h;\,\
u_{1},\ell_{1}), \tag{11}\]  
for chordal distance \(h\geq 0\), scale parameter \(\sigma_{1}>0\), and \(M(\cdot\,;\
u_{1},\ell_{1})\) the Matern correlation function [43; 44] given in Equation (8) with smoothness parameter \(\
u_{1}>0\) and correlation length \(\ell_{1}>0\). For inference, we obtain an empirical semivariogram using the standardized spatial residuals, \(\mathbf{R}_{1}^{*}=(R_{11}^{*},\dots,R_{1n_{1}}^{*})^{T}\), given by Equation (6) when \(i=1\). Using the (univariate) relationship between semivariograms and covariance functions given in Equation (10) when \(i=j=1\), we obtain an estimate, \(\widehat{\theta}\), for \(\theta=\{\sigma_{1},\
u_{1},\ell_{1},\tau_{1}^{2}\}\) through fitting the empirical semivariogram by weighted least squares [47]. As a simplification of cokriging, the kriging prediction at validation location \(\mathbf{s}_{a}\in A\), \(\hat{Z}_{1a}^{\text{kr}}\), and corresponding root-mean-squared prediction error, \(RMSPE\left(\hat{Z}_{1a}^{\text{kr}}\right)\), use only the noisy SIF data \(\mathbf{Z}_{1}\). Let \(\mathbf{R}_{a1}^{*(-A)}\) be the 150-dimensional subset of \(\mathbf{R}_{1}^{*}\) whose values correspond to the 150 data locations found in \(D\) but not in \(A\) that are closest in chordal distance to \(\mathbf{s}_{a}\).