## 2 Materials and Methods
### Modeling Multivariate Spatial Dependence between SIF and XCO\({}_{2}\)
The bivariate spatial variation across the locations associated with \(\mathbf{R}_{0i}^{*}\) and \(\mathbf{R}_{0j}^{*}\), \(i,j\in\{1,2\}\), is described by the variance-covariance matrix,  
\[\mathrm{cov}(\mathbf{R}_{0}^{*})\equiv\mathbf{\Sigma}_{0}= \begin{pmatrix} \mathbf{C}_{0,11}&\mathbf{C}_{0,12}\\ \mathbf{C}_{0,21}&\mathbf{C}_{0,22}\end{pmatrix} + \begin{pmatrix}\frac{\sigma _{\omega 1}^{2}}{\hat{S}_{R,1}^{*}}\mathbf{I}+\mathrm{diag}(\sigma_{\epsilon,01}^{*})& \mathbf{0}\\ \mathbf{0}&\frac{\sigma_{\omega 2}^{2}}{\hat{S}_{R,2}^{*}}\mathbf{I}+\mathrm{diag}( \sigma_{\epsilon,02}^{*})\end{pmatrix} , \tag{9}\]  
where each block \(\mathbf{C}_{0,ij}\) is a 150 \(\times\) 150 (cross-) covariance matrix having elements given by \(C_{ij}(\cdot)\) from Equation (8) with estimated parameters substituted in; \(\mathbf{I}\) is the 150 \(\times\) 150 identity matrix; and \(\sigma^{*}_{e,0i}\) has elements \(\{\sigma^{2}_{e,0i}/\delta^{2}_{R,i}:l=1,\ldots,150\}\) from the definition of \(\epsilon^{*}_{il}\) in Equation (7), corresponding to locations associated with \(\mathbf{R}^{*}_{0i}\). Similarly, we obtain the covariance vector \(\mathbf{c}_{0}=(\mathbf{c}^{T}_{01},\ \mathbf{c}^{T}_{02})^{T}\), with elements of the \(j\)-th sub-vector (\(j=1,2\)) given by \(\mathbf{c}_{0j}\equiv(\text{cov}(\psi^{*}_{1}(\mathbf{s}_{0}),\ R^{*}_{0 jl}):l =1,\ldots,150)^{T}\), where  
\[\text{cov}(\psi^{*}_{1}(\mathbf{s}_{0}),\ R^{*}_{0 jl})=C_{ij}(h_{0l})+\frac{ \sigma^{2}_{\omega,1}}{\delta^{2}_{R,1}}\mathbf{1}\{j=1,\ h_{0l}=0\}; \tag{10}\]  
\(h_{0l}=||\mathbf{s}_{0}-\mathbf{s}_{l}||;\mathbf{s}_{l}\) is the spatial location associated with \(R^{*}_{0 jl};||\cdot||\) denotes chordal distance; and \(\mathbf{1}\{\cdot\}\) is the indicator function. Notice the absence of measurement error variances in Equation (10), since we are predicting the (noiseless) residual process. Then, by minimizing the mean squared prediction error (MSPE) over all linear predictors of \(\psi^{*}_{1}(\mathbf{s}_{0})\) given \(\mathbf{R}^{*}_{0}\), we obtain the cokriging prediction and minimized MSPE (Equations (6) and (7) in [32]):  
\[\widehat{\psi}^{*}_{1}(\mathbf{s}_{0})=\mathbf{c}^{T}_{0}\mathbf{\Sigma}^{-1 }_{0}\mathbf{R}^{*}_{0}, \tag{11}\]  
\[MSPE(\widehat{\psi}^{*}_{1}(\mathbf{s}_{0}))=\frac{\sigma^{2}_{\omega,1}}{ \delta^{2}_{R,1}}+C_{11}(0)-\mathbf{c}^{T}_{0}\mathbf{\Sigma}^{-1}_{0}\mathbf{ c}_{0}. \tag{12}\]  
The estimator given in Equation (11) is an unbiased linear predictor for \(\psi^{*}_{1}(\mathbf{s}_{0})\) under squared error loss, and the uncertainty measure we use is the square root of the MSPE given in Equation (12); notice that \(\widehat{\psi}^{*}_{1}(\mathbf{s}_{0})\) is an approximation to the optimal predictor, \(\mathrm{E}(\psi^{*}_{1}(\mathbf{s}_{0})\ |\ \mathbf{R}^{*})\), since it is based on a subset of the available data. To obtain predictions and uncertainty quantification on the original scale of the latent SIF process, we reverse the standardization to yield  
\[\widehat{Y}_{1}(\mathbf{s}_{0}) =\hat{\mu}_{R,1}+\widehat{\mu}_{1}(\mathbf{s}_{0})+\widehat{\psi }^{*}_{1}(\mathbf{s}_{0})\hat{S}_{R,1}, \tag{13}\] \[RMSPE\Big{(}\widehat{Y}_{1}(\mathbf{s}_{0})\Big{)} =\sqrt{MSPE(\widehat{\psi}^{*}_{1}(\mathbf{s}_{0}))\delta^{2}_{R, 1}}\,. \tag{14}\]  
Operationally, these predictions and their root-mean-squared prediction errors, given by Equations (13) and (14), respectively, can be produced at each location in a highly parallelizable manner. Hence, for each month, our 0.05-degree resolution coSIF data product consists of \(\Big{\{}\widehat{Y}_{1}(\mathbf{s}_{0}):\mathbf{s}_{0}\in D\Big{\}}\) and \(\Big{\{}RMSPE\Big{(}\widehat{Y}_{1}(\mathbf{s}_{0})\Big{)}:\mathbf{s}_{0}\in D \Big{\}}\).