## Appendix A Methodological Detail
### Multivariate Predictive Covariance Matrix Approximation  
Recall that for spatial locations \(\{\mathbf{s}_{1},\ldots,\mathbf{s}_{n_{A}}\}\) in a validation block \(A\), \(F_{1A}\) is the multivariate Gaussian predictive distribution with conditional mean vector \(\mathrm{E}\big{(}\mathbf{Z}_{1A}\mid\mathbf{Z}^{-A}\big{)}\) and conditional variance-covariance matrix \(\mathrm{var}\big{(}\mathbf{Z}_{1A}\mid\mathbf{Z}^{-A}\big{)}\equiv\mathbf{V} _{1A}\) (Section 2.5). In what follows, we approximate \(\mathrm{E}\big{(}\mathbf{Z}_{1A}\mid\mathbf{Z}^{-A}\big{)}\) by the localized predictor \(\widehat{\mathbf{Z}}_{1A}\equiv\Big{(}\widehat{Z}_{11},\ldots,\widehat{Z}_{1 n_{A}}\Big{)}^{T}\), with elements given by the right-hand-side of Equation (13), and \(\mathbf{V}_{1A}\) by the matrix \(\widehat{\mathbf{V}}_{1A}\) whose elements we now derive. The diagonal elements of \(\widehat{\mathbf{V}}_{1A}\) can be obtained via the square of Equation (15); here, we derive the off-diagonal elements of \(\widehat{\mathbf{V}}_{1A}\). Each element is given by  
\[\widetilde{V}_{1A,ab}\equiv\mathrm{cov}\Big{(}Z_{1a},\,Z_{1b}\mid\mathbf{Z}_{ ab}^{-A}\Big{)};\quad a,b\in\{1,\ldots,n_{A}\}, \tag{15}\]  
where the vector \(\mathbf{Z}_{ab}^{-A}\) is composed of sub-vectors \(\mathbf{Z}_{ab,1}^{-A}\) and \(\mathbf{Z}_{ab,2}\). Here, the sub-vector \(\mathbf{Z}_{ab,1}^{-A}\) is a subset of the SIF values in \(\mathbf{Z}_{1}\)_corresponding to the 150 locations not in \(A\) that are closest to \(\mathbf{s}_{a}\)_, and the 150 locations not in \(A\) that are closest to \(\mathbf{s}_{b}\)_. Note that there may be overlap in each set of 150 locations, resulting in a vector of dimension smaller than 300. The sub-vector \(\mathbf{Z}_{ab,2}\) is defined similarly, but it does _not_ exclude the locations in the validation region because its entries are \(\mathrm{XCO}_{2}\) data, not SIF. Using the definition of the noisy data value, \(Z_{1a}\), from Equation (A7), Equation (A15) becomes  
\[\widetilde{V}_{1A,ab}=\hat{S}_{R,1}^{2}\mathrm{cov}\Big{(}\psi_{1}^{*}(\mathbf{s }_{a}),\;\psi_{1}^{*}(\mathbf{s}_{b})\;|\;\mathbf{R}_{ab}^{*(-A)}\Big{)}+ \sigma_{\epsilon,1a}^{2}\mathbf{1}\{a=b\};\quad a,b\in\{1,\ldots,n_{A}\},\] (A16)  
where \(\mathbf{1}\{\cdot\}\) is the indicator function and the variances \(\{\sigma_{\epsilon,1a}^{2}\}_{a}\) are given as described for Equation (1). The vector \(\mathbf{R}_{ab}^{*(-A)}\) has sub-vectors \(\mathbf{R}_{ab,1}^{*(-A)}\) and \(\mathbf{R}_{ab,2}^{*}\), with elements obtained, respectively, from \(\mathbf{Z}_{ab,1}^{-A}\) and \(\mathbf{Z}_{ab,2}\), by the standardization process involving \(\hat{\mu}_{R,1}\), \(\hat{\mu}_{R,2}\), \(\hat{S}_{R,1}\), and \(\hat{S}_{R,2}\) as described in Equation (6). Because the joint distribution of \(\mathbf{R}\) is Gaussian, \(\mathrm{cov}\Big{(}\psi_{1}^{*}(\mathbf{s}_{a}),\,\psi_{1}^{*}(\mathbf{s}_{b} )\;|\;\mathbf{R}_{ab}^{*(-A)}\Big{)}=\mathrm{E}\Big{(}\mathrm{cov}\Big{(}\psi_ {1}^{*}(\mathbf{s}_{a}),\;\psi_{1}^{*}(\mathbf{s}_{b})\;|\;\mathbf{R}_{ab}^{* (-A)}\Big{)}\Big{)}\). Then, by the law of total expectation,  
\[ \begin{split}\mathrm{cov}&\Big{(}\psi_{1}^{*}( \mathbf{s}_{a}),\,\psi_{1}^{*}(\mathbf{s}_{b})\;|\;\mathbf{R}_{ab}^{*(-A)} \Big{)}\\ &=\mathrm{E}\Big{(}\mathrm{cov}\Big{(}\psi_{1}^{*}(\mathbf{s}_{a }),\,\psi_{1}^{*}(\mathbf{s}_{b})\;|\;\mathbf{R}_{ab}^{*(-A)}\Big{)}\Big{)}\\ &=\mathrm{E}\Big{(}\mathrm{E}\Big{\{}\Big{(}\psi_{1}^{*}( \mathbf{s}_{a})-\mathrm{E}\Big{\{}\psi_{1}^{*}(\mathbf{s}_{a})\;|\;\mathbf{R} _{ab}^{*(-A)}\Big{\}}\Big{)}\Big{(}\psi_{1}^{*}(\mathbf{s}_{b})-\mathrm{E} \Big{\{}\psi_{1}^{*}(\mathbf{s}_{b})\;|\;\mathbf{R}_{ab}^{*(-A)}\Big{\}} \Big{)}\;|\;\mathbf{R}_{ab}^{*(-A)}\Big{\}}\Big{)}\\ &=\mathrm{cov}\Big{(}\psi_{1}^{*}(\mathbf{s}_{a})-\mathrm{E} \Big{\{}\psi_{1}^{*}(\mathbf{s}_{a})\;|\;\mathbf{R}_{ab}^{*(-A)}\Big{\}},\; \psi_{1}^{*}(\mathbf{s}_{b})-\mathrm{E}\Big{\{}\psi_{1}^{*}(\mathbf{s}_{b}) \;|\;\mathbf{R}_{ab}^{*(-A)}\Big{\}}\Big{)},\end{split}