## Appendix A Methodological Detail
### Parameter Estimation: Fitting (Cross-) Semivariograms
The result in Equation (A3) follows because the two sets of standardized measurement errors, \(\{e_{i\uparrow}^{\ast}:l=1,\ldots,n_{i\downarrow}\}\), \(i=1,2\), are independent of one another; and they are each independent of the standardized small-scale-variation processes \(\eta_{1}^{\ast}(\cdot)\) and \(\eta_{2}^{\ast}(\cdot)\), and of the standardized micro-scale-variation processes \(\omega_{1}^{\ast}(\cdot)\) and \(\omega_{2}^{\ast}(\cdot)\). Weighted least squares is a popular approach for fitting univariate semivariogram models to their empirical estimates [47] since the method automatically favors lags where the (cross-) spatial dependence is strongest and downweights those lags associated with the fewest spatial pairs. An extension to fitting bivariate semivariogram models is now given. Let \(\Theta\) be a set containing the valid parameter space of the full bivariate Matern covariance model, as established in [41]. Taking the sum of the weighted squared differences between the empirical and model (cross-) semivariograms as the argument to be minimized, we minimize over \(\theta\in\Theta\). For a fixed set of spatial bins centered on \(\{h_{1},\ldots,h_{r}\}\) and a fixed tolerance \(\lambda\geq 0\), a multivariate weighted-least-squares estimate of \(\theta\) is given by  
\[\widehat{\theta}=\arg\min_{\theta\in\Theta}\left\{\sum_{i=1}^{2}\sum_{j=i}^{2 }\sum_{l=1}^{r}|N_{ij}(h_{l};\,\lambda)|\left(\frac{\widehat{\gamma}_{ij}(h_{ l})-\gamma_{ij}(h_{l};\,\,\theta)}{\gamma_{ij}(h_{l};\,\,\theta)}\right)^{2} \right\},\] (A4)  
where \(r\) is the number of spatial bins, \(\widehat{\gamma}_{ij}(\cdot)\) is given in Equation (A2), and \(\gamma_{ij}(\cdot\,\cdot\,\cdot)\) is given in Equation (A3). Note that Equation (A4) can be easily generalized from the bivariate to the multivariate setting. For our purposes, each set \(N_{ij}(\cdot\,\cdot\,\cdot)\) contains all pairs of points in one of the 30 equally spaced bins centered on \(\{h_{1},\ldots,h_{30}\}\), which span between 0 km and a maximum distance of 1000 km. The tolerance \(\lambda\) is determined such that the sets \(\{N_{ij}(h_{l};\,\lambda):l=1,\ldots,30\}\), \(i,j\in 1,2\), are contiguous and non-overlapping. This approach yields numbers of pairs that are well beyond the recommended minimum number of 30 pairs in each set \(N_{ij}(\cdot\,;\,\cdot)\); the importance of this is emphasized in [40], Section 2.4. To preserve total variation in prediction using Equations (11) and (12), each micro-scale-variation component, \(\sigma_{\omega,i}^{2}\), is estimated as  
\[\widehat{\sigma}_{\omega,i}^{2}=\left(\widehat{\tau}_{i}^{2}\widehat{S}_{R,i }^{2}-\widehat{\sigma}_{\omega,i}^{2}\right)\mathbf{1}\Big{\{}\widehat{\tau} _{i}^{2}\widehat{S}_{R,i}^{2}>\widehat{\sigma}_{\omega,i}^{2}\Big{\}};\quad i =1,2,\] (A5)  
where recall that \(\widehat{\tau}_{i}^{2}\) is an element of \(\widehat{\theta}\), \(\widehat{S}_{R,i}\) is given in Equation (6), and \(\mathbf{1}\{\cdot\}\) is the indicator function, which ensures that the estimated variance is non-negative. Recall that the values of \(\{v_{ij}\}_{ij}\), \(\{\ell_{ij}\}_{ij}\), and \(\{\rho_{ij}\}_{ij}\) that result in a positive-semidefinite bivariate covariance function belong to a constrained parameter space [41]. Accordingly, we check a necessary condition for validity of the spatial model given by Equation (8). At each prediction location, \(\mathbf{s}_{0}\), in the spatial domain, \(D\), we check that the joint covariance matrix  
\[\text{cov} \begin{pmatrix}\psi_{1}^{\ast}(\mathbf{s}_{0})\\ \mathbf{R}_{0}^{\ast}\end{pmatrix} = \begin{pmatrix}C_{11}(0)&\mathbf{c}_{0}^{ \ast}\\ \mathbf{c}_{0}&\boldsymbol{\Sigma}_{0}\end{pmatrix} \] (A6)  
is positive-semidefinite. Here, \(C_{11}(\cdot)\) is given by Equation (8); \(\boldsymbol{\Sigma}_{0}\) is a \(300\times 300\) covariance matrix given by Equation (9); and \(\mathbf{c}_{0}\) is a \(300\)-dimensional vector with elements given as described by Equation (10).