## Appendix A Methodological Detail  
### Parameter Estimation: Fitting (Cross-) Semivariograms  
Let \(\psi_{i}^{*}(\cdot)=(Y_{i}(\cdot)-\widehat{\mu}_{i}(\cdot)-\hat{\mu}_{R,i})/\hat{S }_{R,i}\), \(i=1,2\), be the standardized version of the de-trended random process \(\psi_{i}(\cdot)\) defined in Equation (3). The empirical estimator of a stationary (cross-) covariance function, \(\text{cov}(\psi_{i}^{*}(\mathbf{s}+\mathbf{h}),\ \psi_{i}^{*}(\mathbf{s}))\), for \(i,j\in\{1,2\}\), is known as the covariogram for \(i=j\) and the cross-covariogram for \(i\
eq j\). The (cross-) covariogram is a method-of-moments estimator that is highly susceptible to bias and trend contamination ([40], Section 2.4.1). Instead, another measure of multivariate spatial dependence is advocated in [32], namely the (cross-) semivariogram function given by  
\[\gamma_{ij}^{*}(\mathbf{h})\equiv\frac{1}{2}\text{var}(\psi_{i}^{*}(\mathbf{s} +\mathbf{h})-\psi_{j}^{*}(\mathbf{s}));\quad i,j\in\{1,2\}, \tag{1}\]  
for any \(\mathbf{s}\), \(\mathbf{s}+\mathbf{h}\in D\). This is sometimes called the pseudo cross-semivariogram when \(i\
eq j\)[31]. To study spatial variation in \(\psi_{1}^{*}(\cdot)\) and \(\psi_{2}^{*}(\cdot)\), the empirical (cross-) semivariogram can be found in terms of residuals calculated from the data as in Equations (5) and (6) ([40], Section 3.4.3). Under isotropy, the empirical (cross-) semivariogram is given by  
\[\widehat{\gamma}_{ij}(h)=\frac{1}{2|N_{ij}(h;\ \lambda)|}\sum_{l,k(\xi_{il}, \xi_{jk})\in N_{ij}(l_{k}\,\lambda)}\Bigl{(}\Bigl{(}R_{il}^{*}-\overline{R}_{i }^{*}\Bigr{)}-\Bigl{(}R_{jk}^{*}-\overline{R}_{j}^{*}\Bigr{)}\Bigr{)}^{2}, \tag{2}\]  
where \(h=||\mathbf{h}||\geq 0\), \(l=1,\ldots,n_{i}\), \(k=1,\ldots,n_{j}\), and \(i,j\in\{1,2\}\). Here, subtracting \(\overline{R}_{i}^{*}=\frac{1}{n_{i}}\sum_{l=1}^{n_{i}}R_{il}^{*}\) removes any remaining constant-mean variation, and \(|N_{ij}(h;\ \lambda)|\) is the cardinality of \(N_{ij}(h;\ \lambda)\), for \(\lambda\geq 0\). The set \(N_{ij}(h;\ \lambda)\) represents all pairs of points \((\mathbf{s}_{il},\ \mathbf{s}_{jk})\) whose chordal distance, \(h=||\mathbf{s}_{il}-\mathbf{s}_{jk}||\), is in the spatial bin \([h-\lambda,\ h+\lambda]\). We obtain a bivariate semivariogram model to fit to the empirical estimator given in Equation (2) from a convenient relationship between semivariogram functions and covariance functions given in Equation (3) below. Under second-order stationarity assumptions, it is shown in [32] that (cross-) semivariogram functions at spatial lag \(h\) can be obtained from (cross-) covariance functions. The relationship is modified here to account for information about the stochastic micro-scale variation and measurement error that are also contained in the residuals:  
\[\gamma_{ij}(h;\ \theta)=\frac{1}{2}\Bigl{(}\Bigl{(}C_{ii}(0;\ \theta^{\circ})+\tau_{i}^{2}\Bigr{)}+\Bigl{(}C_{jj}(0;\ \theta^{\circ})+\tau_{j}^{2}\Bigr{)}\Bigr{)}-C_{ij}(h;\ \theta^{\circ});\quad i,j\in\{1,2\}, \tag{3}\]  
for \(h\geq 0\), where \(\tau_{i}^{2}=(\sigma_{\omega,i}^{2}+\bar{\sigma}_{e,i}^{2})/\hat{S}_{R,i}^{2}\) is the variability due to micro-scale variation and measurement error, with \(\hat{S}_{R,i}\) given in Equation (6); \(\sigma_{\omega,i}^{2}\) is as described for Equation (3); and \(\bar{\sigma}_{e,i}^{2}\equiv\text{med}(\{\sigma_{e,i}^{2};\ l=1,\ldots,n_{i}\})\) is the median of the measurement error variances given as described for Equation (1). Here, \(\theta^{\circ}\equiv\{\rho_{12},\sigma_{11},\rho_{22},\
u_{11},\
u_{12},\
u_{ 22},\ell_{11},\ell_{12},\ell_{22}\}\) is the set of parameters that completely characterizes \(C_{ij}(\cdot)\), \(i,j\in\{1,2\}\), the full bivariate Matern covariance model given in Equation (8), and \(\theta\equiv\{\theta^{\circ},\ \tau_{1}^{2},\ \tau_{2}^{2}\}\) is the set of parameters thatcompletely characterizes the spatial dependence model.