## Appendix A Methodological Detail
### Multivariate Predictive Covariance Matrix Approximation
, \tag{11}\]  
where each block \(\mathbf{C}_{ab,ji}\), \(i,j\in\{1,2\}\), is a \(150\times 150\) (cross-) covariance matrix having elements \(\{C_{ij}(h_{lk}):l=1,\ldots,150;\,k=1,\ldots,150\}\) for \(C_{ij}(\cdot)\) given in Equation (8) with estimated parameters substituted in, and \(h_{lk}=||\mathbf{s}_{l}-\mathbf{s}_{k}||\) the chordal distance between data locations \(\mathbf{s}_{l}\) and \(\mathbf{s}_{k}\). The locations \(\{\mathbf{s}_{l}:l=1,\ldots,150\}\) used in the predictor \(\widehat{\psi}_{1}^{*}(\mathbf{s}_{a})\) are associated with the values of sub-vector \(\mathbf{R}_{a1}^{*(-A)}\) (when \(i=1\)) or \(\mathbf{R}_{a2}^{*}\) (when \(i=2\)); and the locations \(\{\mathbf{s}_{k}:k=1,\ldots,150\}\) used in the predictor \(\widehat{\psi}_{1}^{*}(\mathbf{s}_{b})\) are associated with the values of sub-vector \(\mathbf{R}_{b1}^{*(-A)}\) (when \(j=1\)) or \(\mathbf{R}_{b2}^{*}\) (when \(j=2\)). Additionally, each block \(\mathbf{Q}_{ab,ii}\), \(i=1,2\), is a \(150\times 150\) matrix with elements  
\[\left\{\left(\frac{\sigma_{\omega,i}^{2}+\sigma_{e,all}^{2}}{\hat{S}_{R,i}^{2 }}\right)\mathbf{1}\{h_{lk}=0\}:\,l=1,\ldots,150;\,k=1,\ldots,150\right\}, \tag{12}\]  
where the measurement error variances \(\{\sigma_{e,all}^{2}\}_{l}\) correspond to locations associated with \(\mathbf{R}_{a1}^{*(-A)}\) (when \(i=1\)) or \(\mathbf{R}_{a2}^{*}\) (when \(i=2\)). The matrix \(\mathbf{\tilde{V}}_{1A}\) has elements given by Equation (10), which we approximate using the result in Equation (11). In fact, we need only compute the upper triangular elements since \(\mathbf{\tilde{V}}_{1A}\) is symmetric. Note that together, \(\mathbf{\tilde{Z}}_{1A}\) and \(\mathbf{\tilde{V}}_{1A}\) fully characterize \(\tilde{F}_{1A}\), the multivariate Gaussian predictive distribution we use for the cokriging prediction at each prediction location \(\mathbf{s}_{a}\in A\) based on the nearest \(150\) values of SIF (not in \(A\)) and the nearest \(150\) values of XCO\({}_{2}\). In the simpler case of kriging prediction, Equation (10) (and hence the weight vector \(\mathbf{\kappa}_{a}^{T}\equiv\mathbf{c}_{a}^{T}\mathbf{\Sigma}_{a}^{-1}\)) is updated using Equation (12). Then, to obtain the variance-covariance matrix for kriging prediction, \(\mathbf{\tilde{V}}_{1A}^{\text{tr}}\), the result in Equation (11) is updated by replacing \(\mathbf{c}_{ab}\) with \(\mathbf{c}_{ab,1}^{\text{tr}}\), \(\mathbf{c}_{ba}\) with \(\mathbf{c}_{ba,1}^{\text{tr}}\), and \(\mathbf{\Sigma}_{ab}\) with \(\mathbf{\Sigma}_{ab,1}^{\text{tr}}\equiv\mathbf{c}_{ab,1}^{\text{tr}}+\mathbf{Q}_ {ab,1}^{\text{tr}}\). The superscript (kr) is used to highlight that the elements of these components are obtained using the univariate covariance model given by Equation (10), rather than the bivariate covariance model given by Equation (10). In the simplest case, trend-surface-only prediction assumes no spatial dependence in \(\psi_{1}^{*}(\cdot)\), so the off-diagonal elements of Equation (11) are identically zero. Hence, the variance-covariance matrix for trend-surface-only prediction, \(\mathbf{\tilde{V}}_{1A}^{\text{tr}}\), is a diagonal matrix with diagonal elements given by the square of Equation (10).