## Appendix A Methodological Detail
### Multivariate Predictive Covariance Matrix Approximation
\] (A17)  
where the final equality uses that \(\mathrm{E}\Big{(}\psi_{1}^{*}(\mathbf{s}_{a})-\mathrm{E}\Big{\{}\psi_{1}^{*}( \mathbf{s}_{a})\;|\;\mathbf{R}_{ab}^{*(-A)}\Big{\}}\Big{)}=0\), for \(a=1,\ldots,n_{A}\). By the properties of the multivariate Gaussian distribution (e.g., [55], Chapter 4),  
\[\mathrm{E}\Big{(}\psi_{1}^{*}(\mathbf{s}_{a})\;|\;\mathbf{R}_{ab}^{*(-A)} \Big{)}\approx\mathbf{c}_{a}^{\tau}\mathbf{\Sigma}_{a}^{-1}\mathbf{R}_{a}^{*(- A)}=\widehat{\psi}_{1}^{*}(\mathbf{s}_{a}),\] (A18)  
with components as defined for Equation (11). The vector \(\mathbf{R}_{a}^{*(-A)}\) has 150-dimensional sub-vectors \(\mathbf{R}_{a1}^{*(-A)}\) (for SIF residuals) and \(\mathbf{R}_{a2}^{*}\) (for \(\mathrm{XCO}_{2}\) residuals), and since each sub-vector is determined with respect to \(\mathbf{s}_{a}\) only, \(\mathbf{R}_{a}^{*(-A)}\) is potentially a subset of \(\mathbf{R}_{ab}^{*(-A)}\), in which case Equation (A18) is an approximation. Denote the cokriging weights in Equation (A18) by \(\mathbf{\alpha}_{a}^{\tau}\equiv\mathbf{c}_{a}^{\tau}\mathbf{\Sigma}_{a}^{-1}\). Then, from Equation (A17) and the approximation in Equation (A18), we obtain  
\[ \begin{split}\mathrm{cov}&\Big{(}\psi_{1}^{*}( \mathbf{s}_{a}),\;\psi_{1}^{*}(\mathbf{s}_{b})\;|\;\mathbf{R}_{ab}^{*(-A)} \Big{)}\\ &\approx\mathrm{cov}\big{(}\psi_{1}^{*}(\mathbf{s}_{a})-\widehat{ \psi}_{1}^{*}(\mathbf{s}_{a}),\;\psi_{1}^{*}(\mathbf{s}_{b})-\widehat{\psi}_{1 }^{*}(\mathbf{s}_{b})\big{)}\\ &=\mathrm{cov}\Big{(}\psi_{1}^{*}(\mathbf{s}_{a})-\mathbf{ \alpha}_{a}^{\tau}\mathbf{R}_{a}^{*(-A)},\;\psi_{1}^{*}(\mathbf{s}_{b})- \mathbf{\alpha}_{b}^{\tau}\mathbf{R}_{b}^{*(-A)}\Big{)}\\ &=\mathrm{cov}(\psi_{1}^{*}(\mathbf{s}_{a}),\;\psi_{1}^{*}( \mathbf{s}_{b}))-\mathrm{cov}\Big{(}\psi_{1}^{*}(\mathbf{s}_{a}),\;\mathbf{R}_{ b}^{*(-A)}\Big{)}\mathbf{\alpha}_{b}\\ &\quad\quad-\mathbf{\alpha}_{a}^{\tau}\mathrm{cov}\Big{(}\mathbf{ R}_{a}^{*(-A)},\;\psi_{1}^{*}(\mathbf{s}_{b})\Big{)}+\mathbf{\alpha}_{a}^{\tau} \mathrm{cov}\Big{(}\mathbf{R}_{a}^{*(-A)},\;\mathbf{R}_{b}^{*(-A)}\Big{)} \mathbf{\alpha}_{b}\\ &=C_{11}(h_{ab})+\frac{\sigma_{\epsilon,1}^{2}}{\hat{S}_{R,1}^{2}} \mathbf{1}\{a=b\}-\mathbf{\alpha}_{ab}^{\tau}\mathbf{\alpha}_{b}-\mathbf{ \alpha}_{a}^{\tau}\mathbf{\alpha}_{b}+\mathbf{\alpha}_{a}^{\tau}\mathbf{\Sigma}_{ ab}\mathbf{\alpha}_{b},\end{split} \] (A19)  
where \(C_{11}(\cdot)\) is the covariance function given in Equation (8); \(h_{ab}=||\mathbf{s}_{a}-\mathbf{s}_{b}||\) is the chordal distance between validation locations \(\mathbf{s}_{a}\) and \(\mathbf{s}_{b}\); \(\sigma_{\epsilon,1}^{2}\) is the variance of the micro-scale process \(\omega_{1}(\cdot)\); \(\hat{S}_{R,1}\) is given in Equation (6) when \(i=1\); and \(\mathbf{1}\{\cdot\}\) is the indicator function. The vector \(\mathbf{c}_{ab}\) is made up of sub-vectors \(\mathbf{c}_{ab,1}\equiv\Big{(}\mathrm{cov}\Big{(}\psi_{1}^{*}(\mathbf{s}_{a}),\;R_{ b1I}^{*(-A)}\Big{)}:l=1,\ldots,150\Big{)}^{\top}\) and \(\mathbf{c}_{ab,2}\equiv\big{(}\mathrm{cov}\big{(}\psi_{1}^{*}(\mathbf{s}_{a}),\;R_{ bI}^{*}\big{)}:l=1,\ldots,150\big{)}^{\top}\), where \[\text{cov}\Big{(}\psi_{1}^{*}(\mathbf{s}_{a}),\,R_{b1I}^{*(-A)}\Big{)} =C_{11}(h_{ab})+\frac{\sigma_{\omega,1}^{2}}{\hat{S}_{R,1}^{2}} \mathbf{1}\{h_{ab}=0\}, \tag{10}\] \[\text{cov}(\psi_{1}^{*}(\mathbf{s}_{a}),\,R_{b2I}^{*}) =C_{12}(h_{ab}).\]  
We obtain the vector \(\mathbf{c}_{ba}=(\mathbf{c}_{ba,1}^{T},\,\mathbf{c}_{ba,2}^{T})^{T}\) similarly. The covariance matrix \(\mathbf{\Sigma}_{ab}\) in Equation (11) is given by  
\[\mathbf{\Sigma}_{ab}\equiv\text{cov}\Big{(}\mathbf{R}_{a}^{*(-A)},\,\mathbf{R}_{b }^{*(-A)}\Big{)}= \begin{pmatrix}\mathbf{C}_{ab,11}&\mathbf{C}_{ab,12}\\ \mathbf{C}_{ab,21}&\mathbf{C}_{ab,22}\end{pmatrix} + \begin{pmatrix}\mathbf{Q}_ {ab,11}&\mathbf{0}\\ \mathbf{0}&\mathbf{Q}_{ab,22}\end{pmatrix}