## 2 Materials and Methods
### Modeling Multivariate Spatial Dependence between SIF and XCO\({}_{2}\)  
In what follows, we discuss the modeling approach for a single pair of months (e.g., SIF in July and XCO\({}_{2}\) in August). Let \(D\) be the spatial region of interest shown in Figure 3 and described in Section 2.2, tessellated by 0.05-degree-resolution grid cells. Then, \(\mathbf{s}\in D\) refers to the centroid \(\mathbf{s}\) of a grid cell in \(D\). The true SIF process averaged over a given month is \(\{Y_{1}(\mathbf{s}):\mathbf{s}\in D\}\), and the true XCO\({}_{2}\) process averaged over the following month is \(\{Y_{2}(\mathbf{s}):\mathbf{s}\in D\}\), but the OCO-2 satellite observations are imperfect and incomplete on \(D\). We denote the gridded and monthly averaged SIF and XCO\({}_{2}\) data (subsequently simply referred to as \"data\") by the vectors \(\mathbf{Z}_{1}\) and \(\mathbf{Z}_{2}\), respectively, with \(\mathbf{Z}_{i}\equiv(Z_{i1},\ldots,Z_{in_{i}})^{T}\), \(i=1,2\), where there are \(n_{1}\) grid cells that contain SIF data and \(n_{2}\) grid cells that contain XCO\({}_{2}\) data. In Figure 3 (Top), \(n_{1}=40\),780, and in Figure 3 (Bottom), \(n_{2}=15\),514; see also Figure 2. In the spatial domain \(D\), we model the noisy data value \(Z_{il}\) at grid cell location \(\mathbf{s}_{il}\in D\) as  
\[Z_{il}=Y_{i}(\mathbf{s}_{il})+\epsilon_{il};\quad l=1,\ldots,n_{i},\ i=1,2, \tag{1}\]  
where \(\{\epsilon_{il}\}_{il}\) is a set of (non-stationary) independent measurement errors, with mean zero and variances \(\{\sigma_{\varepsilon,l}^{2}\}_{il}\) obtained from the Lite files and averaged as described in Section 2.2. The latent processes \(Y_{1}(\cdot)\) and \(Y_{2}(\cdot)\) are defined on all \(\mathbf{s}\in D\) and are modeled as  
\[Y_{i}(\mathbf{s})=\mu_{i}(\mathbf{s})+\psi_{i}(\mathbf{s});\quad i=1,2. \tag{2}\]  
Here, \(\mu_{i}(\cdot)\equiv\mathrm{E}(Y_{i}(\cdot))\) represents the deterministic trend component, which can be thought of as _large-scale variation_, while \(\psi_{i}(\cdot)\) is a random process with zero mean that accounts for any remaining spatial variation. For \(\mathbf{s}\in D\), we write  
\[\psi_{i}(\mathbf{s})=\eta_{i}(\mathbf{s})+\omega_{i}(\mathbf{s});\quad i=1,2, \tag{3}\]  
where \(\eta_{i}(\cdot)\) is a mean-zero Gaussian process that captures the _smooth small-scale variation_, and \(\omega_{i}(\cdot)\) is a process uncorrelated at the grid scale with mean zero and constant variance \(\sigma_{\omega,j}^{2}\) that accounts for _micro-scale variation_. Because our objective is spatial prediction, we take a relatively parsimonious approach for modeling \(\mu_{i}(\cdot)\) to avoid over-fitting the data. At each location \(\mathbf{s}\in D\), estimates for the deterministic trend components, \(\widehat{\mu}_{i}(\mathbf{s})\), \(i=1,2\), are obtained through regression on a vector of covariates, \(\mathbf{x}_{i}(\mathbf{s})\), as follows:  
\[\widehat{\mu}_{i}(\mathbf{s})=\mathbf{x}_{i}(\mathbf{s})^{T}\widehat{\mathbf{ \beta}}_{i};\quad i=1,2, \tag{4}\]  
where \(\widehat{\mathbf{\beta}}_{i}\) is a vector of ordinary-least-squares-estimated regression coefficients obtained from fitting \(\{\mathbf{x}_{i}(\mathbf{s}_{il})^{T}\widehat{\mathbf{\beta}}_{i}\}_{l}\) to \(\{Z_{il}\}_{l}\) separately for \(i=1,2\). For each \(i=1,2\), the elements of \(\mathbf{x}_{i}(\cdot)\) include an intercept term, and \(k_{i}\) spatial bisquare basis functions (e.g., [38]), here used to express a spatial trend. We obtained the bisquare basis functions from the package \(FRK\) in R [39] and chose them such that their spatial supports cover the data region with centers arranged in a regular, \(6\times 10\) grid. Some basis functions over water were not appropriate for the land-based SIF data; three were identified and withheld from \(\mathbf{x}_{1}(\cdot)\), the vector of covariates for SIF. Hence, \(k_{1}=57\) and \(k_{2}=60\).