## 2 Materials and Methods
### Statistical Validation of the coSIF Data Product
Because we assume the data are Gaussian, each \(F_{1a}\) is a Gaussian distribution with mean \(\operatorname{E}\left(Z_{1a}\mid\mathbf{Z}^{-A}\right)\approx\widehat{Z}_{1a}\) and variance \(\operatorname{var}\left(Z_{1a}\mid\mathbf{Z}^{-A}\right)\approx RMSPE\left( \widehat{Z}_{1a}\right)^{2}\); see Equations (13) and (15), respectively. Recall that the approximations are due to the predictor \(\widehat{Z}_{1a}\) using only the nearest 150 data points for SIF and for XCO\({}_{2}\), and that the 150 data points for SIF are outside of the validation block \(A\). Through the predictive distribution, \(F_{1a}\), we can quantify the uncertainty in our Level 3 coSIF data product. We compare these quantified uncertainties to analogous ones obtained from the univariate kriging model and the trend-surface-only model using scoring rules, which measure how well a probabilistic prediction aligns with what is actually observed. Here, we consider the interval score which assesses calibration and sharpness by favoring narrow prediction intervals and penalizing those that do not contain the true observation. With a 95% two-sided prediction interval given by \([p,q]\subset\mathbb{R}\) for \(Z_{1a}\) at validation location \(\mathbf{s}_{a}\in A\), the interval score (INT) is given by  
\[INT(p,q;\;Z_{1a})\equiv(q-p)+\frac{2}{\alpha}(p-Z_{1a})\mathbf{1}\{Z_{1a}<p\} +\frac{2}{\alpha}(Z_{1a}-q)\mathbf{1}\{Z_{1l}>q\}, \tag{18}\]  
for \(a=1,\dots,n_{A}\), where \(\alpha=0.05\) and \(\mathbf{1}\{\cdot\}\) is the indicator function [51]. Another scoring rule that assesses both the first and second moments of a predictive distribution is the Dawid-Sebastiani score (DSS; [52]) defined as \[DSS(F_{1a},\,Z_{1a})\equiv\frac{\big{(}Z_{1a}-\operatorname{E}\big{(}Z_{1a}\mid \mathbf{Z}^{-A}\big{)}\big{)}^{2}}{\operatorname{var}(Z_{1a}\mid\mathbf{Z}^{-A} )}+2\log\sqrt{\operatorname{var}(Z_{1a}\mid\mathbf{Z}^{-A})}\, \tag{19}\]  
for \(a=1,\ldots,n_{A}\), where recall that \(F_{1a}\) is a distribution with mean \(\operatorname{E}\big{(}Z_{1a}\mid\mathbf{Z}^{-A}\big{)}\approx\widehat{Z}_{1a}\) and variance \(\operatorname{var}\big{(}Z_{1a}\mid\mathbf{Z}^{-A}\big{)}\approx RMSPE \Big{(}\widehat{Z}_{1a}\Big{)}^{2}\). Because \(F_{1a}\) is Gaussian, Equation (19) is functionally equivalent to the logarithmic score, which is one of the most-used proper scoring rules [53]. The INT and DSS are computed pointwise, for each validation location, and typically the average of each collection of scores is used to compare the different prediction frameworks, without any regard to the spatial dependence between them. The validation summaries introduced above focus on an average of pointwise predictions, but since the predictions are spatially correlated within a given block \(A\), it is also important to consider multivariate scoring rules. These scoring rules assess \(\widehat{\mathbf{Z}}_{1A}\) as _a single multivariate prediction_ of \(\mathbf{Z}_{1A}\) through the corresponding multivariate predictive distribution denoted as \(F_{1A}\). Again, because the process model in Equation (2) is Gaussian, \(F_{1A}\) is a multivariate Gaussian distribution with conditional mean vector \(\operatorname{E}\big{(}\mathbf{Z}_{1A}\mid\mathbf{Z}^{-A}\big{)}\) and conditional variance-covariance matrix \(\operatorname{var}\big{(}\mathbf{Z}_{1A}\mid\mathbf{Z}^{-A}\big{)}\equiv \mathbf{V}_{1A}\). Hence, \(F_{1A}\) can be assessed using the multivariate DSS:  
\[MDSS(F_{1A},\,\mathbf{Z}_{1A})\equiv\Big{(}\mathbf{Z}_{1A}-\operatorname{E} \Big{(}\mathbf{Z}_{1A}\mid\mathbf{Z}^{-A}\Big{)}\Big{)}^{T}\mathbf{V}_{1A}^{- 1}\Big{(}\mathbf{Z}_{1A}-\operatorname{E}\Big{(}\mathbf{Z}_{1A}\mid\mathbf{Z }^{-A}\Big{)}\Big{)}+\log\det\mathbf{V}_{1A}. \tag{20}\]  
The MDSS provides a scoring metric to assess multivariate predictive distributions (e.g., [51], Equation (25)).