{'fraction_non_alphanumeric': 0.037791034402872695, 'fraction_numerical': 0.010540947526931542, 'mean_word_length': 4.472741679873217, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multivariate Gaussian Process Regression for Multiscale Data Assimilation and Uncertainty Reduction

We present a multivariate Gaussian process regression approach for parameter field reconstruction based on the field's measurements collected at two different scales, the coarse and fine scales. The proposed approach treats the parameter field defined at fine and coarse scales as a bivariate Gaussian process with a parameterized multiscale covariance model. We employ a full bivariate Mat\'{e}rn kernel as multiscale covariance model, with shape and smoothness hyperparameters that account for the coarsening relation between fine and coarse fields. In contrast to similar multiscale kriging approaches that assume a known coarsening relation between scales, the hyperparameters of the multiscale covariance model are estimated directly from data via pseudo-likelihood maximization. We illustrate the proposed approach with a predictive simulation application for saturated flow in porous media. Multiscale Gaussian process regression is employed to estimate two-dimensional log-saturated hydraulic conductivity distributions from synthetic multiscale measurements. The resulting stochastic model for coarse saturated conductivity is employed to quantify and reduce uncertainty in pressure predictions.


Introduction
Kriging is a widely used geostatistical tool for estimating the spatially distributed quantities from sparse measurements and for quantifying uncertainty of such estimates. When such estimates are used as parameter fields in a physics model, uncertainty in kriging estimates can be propagated throughout the model to quantify uncertainty in predictions. Observations of parameter fields are often gathered at different scales (often referred to as support scales). A common example is hydraulic permeability. Laboratory experiments performed on field samples (which are approximately 0.1 m large) provide values with the support equal to the sample size. On the other hand, indirect measurements such as cross-pumping tests have a support volume equaled to the distance between wells (10 −100 m), and correspond to local averages of the measured quantity over that support volume. Methodologies for incorporating multiscale measurements on stochastic models for spatial parameters are important in order to systematically reduce uncertainty using all available data sources.
In Li et al. [2009], the authors employ multivariate Gaussian process regression or cokriging, together with the properties of Gaussian processes, to formulate a multiscale simple kriging method for incorporating observations of the field at two scales under the assumption that coarsening is obtained by block averaging. This multiscale simple kriging method treats the fine and coarse scale fields as covariates in bivariate cokriging, with the coarse covariance and coarse-fine cross-covariance computed analytically from the fine scale covariance in terms of the block averaging operator [Vanmarcke, 2010]. While this approach to formulating multiscale simple kriging can be generalized to all linear coarsening operators (e.g. arithmetic moving averages, convolutions, etc.), it's not possible in general to obtain closed-form expressions for the coarse covariance and the coarsefine cross-covariance. Furthermore, this approach either presumes a priori knowledge of the structure of the coarsening operator and the coarse field support volume, which is in general application-specific and not available, or requires a sufficiently flexible parameterized linear coarsening operator with parameters to be estimated together with all other covariance hyperparameters.
As an alternative, in this work we propose to model the fine and coarse fields as components of a bivariate Gaussian process with parameterized bivariate covariance model and to estimate the shape and smoothness parameters of such model directly from multiscale data rather than deriving them from knowledge of the structure of the coarsening operator. Multivariate covariance models are commonly used for cokriging in geostatistics. We point the reader to Genton and Kleiber [2015] for a comprehensive review of multivariate models. For multiscale simple kriging, we propose employing the univari-ate Matérn kernel as the building block for multiscale covariance models as it provides sufficient flexibility for capturing the smoothness properties that arise from coarsening. In particular, we employ the full bivariate Matérn covariance kernel proposed in Gneiting et al. [2010]. Our approach can be extended to more than two support scales by employing the generalized construction proposed in Apanasovich et al. [2012]. We illustrate via numerical experiments that the full bivariate Matérn kernel captures the different degrees of smoothness of each scale and, via multiscale simple kriging, produces accurate estimates of the field at fine and coarse scales by assimilating data collected at fine and coarse scales.
This manuscript is structured as follows. In section 2, we introduce the multiscale simple kriging algorithm based on multivariate Gaussian process regression, and outline the model selection approach employed. In section 3, we summarize the block averagingbased covariance model, and introduce the bivariate Matérn model. Numerical experiments illustrating both inference from data and application to saturated flow modeling are presented in section 4. Finally, conclusions are given in section 5, together with possible avenues for future work.

Multiscale Simple Kriging
Let Ω ⊂ R d (d ∈ [1, 3]) be the subset of physical space of interest, and let Y : Ω → R denote a heterogeneous physical property to be estimated from discrete observations. We assume that measurements of Y are available at two different scales: (i) from soil samples (fine scale) with the measurement scale η f , and (ii) pumping-well tests (coarse scale) with the measurement scale η c . Our goal is to incorporate both sets of measurements to derive estimators of Y at either support scale.
We use y f to denote Y sampled at the support length-scale η f and y c for Y sampled at the support length-scale η c . Next, we assume that a "coarsening" operator R cf exists such that A total of N f and N c observations of y f and y f are available at locations , respectively. The observation locations are arranged into the N f × d and N c ×d matrices X f and X c , respectively. Similarly, we arrange the observation into the N f ×1 and N c ×1 column vectors y f s and y c s , respectively. Finally, we denote by D s ≡ (X f s , y f s , X c s , y c s ) the multiscale data set.
We assume that y f and y c are Gaussian processes. For univariate simple kriging, we employ the prior processes which are updated employing the single-scale data sets (X c s , y c s ) and (X f s , y f s ), respectively [Cressie, 2015]. Here, C c and C f denote parameterized covariance functions, with hyperparameter sets θ c and θ f , respectively, estimated from the corresponding data sets. For multiscale simple kriging [Li et al., 2009], we assume that the fine and coarse scale fields are components of a bivariate Gaussian process y( to be conditioned on the multiscale data set D s . Here, C denotes the parameterized multiscale covariance function, with hyperparameter set θ to be estimated from the multiscale data set. The multiscale simple kriging-based estimation procedure is as follows: First, a "valid" multiscale covariance model C(·, · | θ) is selected as explained later in this section. Second, an estimator θ * of the hyperparameters of the covariance model is computed from the multiscale data set D s (model selection stage). Finally, the estimator θ * is plugged into the covariance model, and the mean and covariance of the conditional bivariate GP y | D s , θ * are computed (plug-in conditioning stage).
The choice of multiscale covariance model and hyperparameters is not arbitrary. First, the model must be "valid", i.e. the block components C c , C f and C cf (≡ (C f c ) ) and the hyperparameters θ must be so that C is a symmetric positive definite (SPD) kernel. Second, in a weaker sense, covariance model and hyperparameters must reflect that y c is a coarsening of y f as given by eq. (1). The approach we pursue in this work is different from the approach of Li et al. [2009] and Cressie [2015], where block averaging with known support is employed as the coarsening operator, from which block components can be derived such that the previous conditions are satisfied. In contrast, in this work we propose a multiscale covariance model that can be used even when the support and/or the structure of the coarsening are unknown a priori.
For the remainder of this manuscript we employ the symbol C(·, · | θ) to denote both covariance matrices and functions, i.e., for N a × d and N b × d matrices X a and

Conditioning on Multiscale Data
For given hyperparameters θ, the posterior process conditional on the data D s is y | D s , θ ∼ GP(μ,Ĉ), with conditional mean and covariance given bŷ where y s = [(y c s ) , (y f s ) ] is the N × 1 vector of observations, N ≡ N f + N c , and C s (θ) ≡ C(X s , X s | θ) denotes the covariance matrix of the observations. From eq. (3) and eq. (4) we can find the conditional mean and covariance for each scale. For the coarse scale, we have y c | D s , θ ∼ GP(μ c ,Ĉ c ), witĥ Similar expressions can be found for the fine scale.

Model Selection
For the plug-in approach employed in this work we need to compute an estimate θ * of the covariance hyperparameters θ from the data set D s . Two approaches are commonly used in univariate geostatistics: variogram curve fitting, and (pseudo-log)likelihood maximization. In this work we employ the latter as experience shows that curve fitting is not adequate for estimating the smoothness of a differentiable processes [Stein, 1999], which plays an important role in estimating coarse fields (see section 3.1). Specifically, we take θ * ≡ arg max θ L(D s , θ), where L(D s , θ) is a (pseudo-)log likelihood function.
Two (pseudo-)log likelihood functions are employed in this work. The first one is the marginal likelihood (ML), commonly used in Bayesian model selection, and defined as the likelihood of the observations given the data and the hyperparameters: The computational cost of computing the ML log likelihood is O(N 3 ), chiefly driven by the inversion of the covariance matrix of the observations. We note that ML estimation is equivalent to the well-known residual maximum likelihood (REML) [Cressie, 2015] estimation for the case of zero mean priors, as in eq. (2).
The second function is the leave-one-out cross-validation (LOO-CV) pseudo-log likelihood, which consists of sum of the log predictive probabilities of each observation given all other observations and the hyperparameters. We denote by D −i = (X −i , y −i ) the data set excluding the ith observation. The prediction for the ith observation given D −i and θ has distribution p(y i |D s , θ) = N (y i |μ −i ,σ −i ), whereμ −i andσ −i are given by applying eq. (3) and eq. (4) to D −i , obtaininĝ The LOO-CV pseudo-log likelihood is therefore given by We note that computing the LOO-CV pseudo-likelihood does not require applying eq. (3) and eq. (4) N times, and instead can be computed with O(N 2 ) overhead from C −1 s (θ). Therefore, the computational cost of the LOO-CV pseudo-likelihood is O(N 3 ) [Rasmussen and Williams, 2005].

Multiscale Covariance Models
In this section we present bivariate covariance models to be used in multiscale simple kriging. The first model is the block-averaging-based model employed in Cressie [2015] and Li et al. [2009]. The second model is a bivariate Matérn covariance model. Both models require assuming a certain functional form for the fine-scale covariance (with hyperparameters estimated during model selection), but the block-averaging-based model requires a priori knowledge of the the structure of the coarsening, whereas the bivariate Matérn covariance model doesn't.

Block Averaging-Based Model
In this approach, we construct a multiscale covariance model directly from a model for the fine scale covariance C f and the relation eq. (1). We assume that R cf is the linear operator H ηc of the form where h ηc (x, z) is an averaging kernel with support η c satisfying Ω h(z, z ) dz = 1 for any z ∈ Ω. Linearity of the coarsening operator ensures that y c is a GP. If h ηc is known, the coarse-scale covariance and the fine-coarse-scale cross-covariance functions can be computed as In this model, the multiscale covariance hyperparameters are those of the fine scale covariance model. This approach was employed in Li et al.
[2009] to formulate the authors' multiscale simple kriging method, with the block averaging kernel for which explicit formulae for C c and C cf can be found in terms of C f [see Vanmarcke, 2010, chap. 5-7]. It must be noted that for the block averaging kernel eq. (13), y c is in general not isotropic, and only stationary away from ∂Ω.
This approach can be employed when the support length-scale is known and when eq. (10) accurately describes the averaging relation between y c and y f . The advantage of this approach is that the number of hyperparameters is lower than for a full multiscale model (as will be seen in section 3.2), thus simplifying the model selection process.
To illustrate this approach, we compute the coarse covariance and coarse-fine crosscovariance using the block averaging kernel eq. (13) for the fine scale isotropic covariance C f (r) = σ 2 f exp(−r/λ f ) in R 2 , and for various values of η c /λ f . Fine and coarse variograms and coarse-fine cross-variograms are presented in fig. 1. In this case, the coarse covariance is not isotropic, therefore we present in fig. 1 where e 1 is the unit vector along the 1st coordinate. A similar pseudo-isotropic variogram is computed from the coarse-fine cross-covariance.  As expected, the coarse covariance exhibit longer correlation length and smaller variance than the fine scale covariance. We also observe that the coarse correlation length increases and the coarse variance decreases with increasing η c /λ f , as noted in Tartakovsky et al. [2017]. Similar behavior is observed for the coarse-fine cross-covariance. Therefore, applying the multiscale simple kriging estimation method described in section 2 to the block-averaging model with fixed η c when the coarse length scale is not sufficiently well known may result in large estimation errors. This issue can be ameliorated by including η c as a hyperparameter of the multiscale covariance model to be estimated during the model selection stage. This requires assuming that the coarsening operator can be approximated by the block averaging operator with known averaging kernel.
Another effect of coarsening is in the change in smoothness from y f to y c , reflected by the order of mean square differentiability of the fields. This is seen visually in fig. 1 in the behavior of γ f andγ c for r → 0. For the exponential kernel, we have γ f = σ 2 f λ −1 f |r|+ O(r 2 ) for r → 0. According to theorem 2.7.2 in Stein [1999], under this condition, y f is not mean square differentiable. On the other hand, dγ c (r = 0)/dr ≈ 0 (indicated by the 'S'-shape of the pseudo-variogram), which indicates that y c is mean square differentiable. This change in smoothness is due to coarsening and correspond to high correlation between y c values at locations closely together, with the radius of such high correlation region increasing with increasing η c /λ f . Multiscale covariance models must capture this 'S'-shape in order to ensure that the coarse field estimator has the appropriate degree of smoothness due to coarsening.

Matérn Model
As an alternative to the block-averaging-based model presented in the previous section, we propose employing multivariate Matérn kernels [Gneiting et al., 2010] to model block components of the multiscale covariance function. The univariate Matérn kernel is given by where K ν is the modified Bessel function of the second kind, ν > 0 is a shape parameter, D is a d × d SPD matrix, and r D denotes the so-called Mahalanobis distance, given by The matrix D allows us to introduce anisotropy to the model. In the isotropic case, D ≡ λ −2 I d , where λ denotes the correlation length. The shape parameter ν controls the smoothness of the process. Specifically, a process with Matérn covariance is m times mean square differentiable if and only if ν > m. The Matérn family includes other well known kernels as special cases: For ν = 1/2, the Matérn kernel is equal to the exponential kernel exp(− r D /λ) (not differentiable), while for ν → ∞ we recover the Gaussian or squaredexponential kernel exp(− r 2 D /2λ 2 ) (infinitely differentiable).
In this work we propose employing the so-called full bivariate Matérn kernel as the multiscale covariance model, given by where ν c , ν f and ν cf are the shape parameters for the coarse, fine and coarse-fine components, respectively; λ c , λ f and λ cf are the correlation lengths of the coarse, fine, and coarse-fine components, respectively; σ c and σ f are the coarse and fine scale standard deviations, respectively, and ρ is the collocated correlation coefficient. The nugget terms σ 2 nc 1 r =0 and σ 2 nf 1 r =0 are added to the block diagonal components to capture Gaussian measurement noise. The full bivariate model, by including separate standard deviation, shape parameters and correlation lengths for the coarse and fine scale components, provides enough flexibility to capture the effects of coarsening discussed in the previous section, namely the changes in standard deviation, correlation length, and degrees of smoothness, from the fine to coarse scale.
We proceed to establish conditions for the hyperparameters that ensure that the multiscale covariance model is valid, that is, that the full bivariate Matérn function eq. (16) is SPD. To simplify the full bivariate model, we consider exclusively the case ν cf = 1 2 (ν c + ν f ). Numerical experiments show us that this case provides sufficient flexibility for modeling multiscale covariance functions for sufficient scale separation between y f and y c , i.e. for η c /λ f 1. For this case, sufficient conditions for validity are Gneiting et al. [2010] The number of hyperparameters to estimate during the model selection stage is 10 (λ c , λ f , λ cf , ν c , ν f , σ c , σ f , ρ, σ nc , and σ nf ). In general, we expect that σ c < σ f , ν c > ν f , λ c > λ f , and ρ > 0. These are not validity conditions and rather reflect the fact that y c is a coarsening of y f that is more smooth and has lower variance, and that y c and y f are by definition positively correlated.
The advantage of employing the full bivariate Matérn kernel as a multiscale covariance model is that it does not require a priori knowledge of the coarsening relation (e.g., whether is a linear coarsening, the analytic form of h ηc , the value of η c , etc.). Instead, the shape of the block variogram components is estimated from data, and reflected by the shape parameters ν c , ν f , and ν cf . The disadvantages are that the number of free hyperparameters is larger than for the block averaging model, and that model selection requires constrained minimization of the pseudo-likelihood subject to the nonlinear constraints eq. (17) and eq. (18), as opposed to unconstrained minimization.

Numerical Experiments
In this section we illustrate the proposed multiscale simple kriging approach that combines cokriging with the full bivariate Matérn kernel. We employ the proposed approach to estimate two-dimensional distribution of hydraulic conductivity K from synthetic multiscale hydraulic conductivity measurements. We model K at the fine scale as a stationary lognormal process, i.e., K = K G exp(Y ), where K G [L T −1 ] is the uniform geometric mean, assumed known, and Y , referred to as the normalized log-conductivity, is a zero-mean Gaussian process with isotropic exponential covariance.

Model Selection
We illustrate the model selection process by considering two sets of values of η f and η c , corresponding to two degrees of scale separation between y f and y c . For each test case, we obtain hyperparameter estimates using both the ML and LOO-CV approaches from multiscale observations. We compare the estimated coarse and fine variograms against the corresponding empirical variograms computed using the R package RandomFields [Schlather et al., 2017[Schlather et al., , 2015. We also compare against the true fine scale variogram and the coarse and coarse-fine pseudo-isotropic variograms.

For this experiment, the estimation domain is
is the domain length scale. A synthetic reference fine scale log-conductivity field is simulated on a 256×128 square grid. We take the fine support scale to be equal to the grid resolution ∆x, that is, η f = ∆x. The reference coarse scale log-conductivity field is then computed by block averaging with η c = m∆x, where m is a positive integer. At each scale, sampling locations are randomly chosen from the set of discretization centroids. White noise with variance σ 2 is added to the observations at both scales. All parameters and the number of measurements are listed in in table 1.
Constrained minimization of the negative log pseudo-likelihood is performed employing the MATLAB implementation of the interior point algorithm with a BFGS approximation of the Hessian. Correlation length, shape, and standard deviation hyperparameters are log-transformed to change their support from (0, ∞) to (−∞, ∞). Similarly, the collocated correlation coefficient is transformed as ξ = log[(1+ρ)/(1−ρ)] to change their support from (−1, 1) to (−∞, ∞). Initial values for the hyperparameters are chosen by inspection of the empirical variograms, and modified if the constrained minimization algorithm fails to converge. True parameters in italics (σc and ρ) are computed from eq. (11) and eq. (12) Test 1  fig. 2 shows both the reference fine and coarse scale fields and the observation locations. Similarly, fig. 3 presents the estimated, true, and empirical variograms. For this test we have significant scale separation as shown by the 46% variance decrease from y f to y c . From table 2 and fig. 3, we see that both ML and LOO-CV provide accurate estimates of the fine scale and coarse scale standard deviations and the collocated correlation coefficient, with LOO-CV overestimating the standard deviations. Both ML and LOO-CV overestimate λ f by approximately 30%, which indicates that more finescale observations are necessary to more accurately estimate the fine-scale correlation length. More importantly, both ML and LOO-CV correctly identify the difference in smoothness between scales. Even though ML and LOO-CV overestimate the fine scale shape parameter, both correctly estimate that the fine scale field is not mean-square differentiable. Furthermore, both methods estimate coarse correlation length and shape parameter larger than their fine scale counterparts, aligned with the observations in section 3.
We compute the mean and variance for the fine and coarse field conditioned on multiscale data by employing the Nyström method presented in appendix A: . The results are shown in fig. 4. It can be seen in fig. 4c and fig. 4d how multiscale measurements inform the conditioning for both scales. The measurements at a given scale reduce the conditional variance to the level of noise, while the measurements at the other scale reduce variance to a lesser but still noticeable degree.
For Test 2 we have less differentiation between scales compared to Test 1. fig. 5 shows both the reference fine and coarse scale fields. fig. 6 shows the estimated, true, and em-  pirical variograms. Additionally, fig. 7 shows the conditional mean and variance. Due to the lesser degree of separation between scales compared to Test 1, model selection is more challenging. Nevertheless, as for Test 1, we find that both ML and LOO-CV provide estimates for the full bivariate Matérn model that accurately resolve the structure of the multiscale reference model. We see that both ML and LOO-CV overestimate the shape parameter of the fine scale field, but correctly identify that the coarse field is smoother and has longer correlation length. As the degree of separation between scales is lesser, we can see in fig. 4c and fig. 4d that measurements on the fine scale more sharply reduce the conditional variance at the coarse scale, and that measurements on the coarse scale more sharply reduce the conditional variance at the fine scale.  We now illustrate how the estimated hyperparameters of the multivariate Matérn model vary for randomly chosen reference fields and observation locations. For this experiment we generate 500 draws of normalized fine scale log-conductivity fields from the same GP and compute for each one the corresponding coarse scale field via block aver- Results for Test 3 are summarized in fig. 8. These results, as those of Tests 1 and 2, indicate that the proposed approach is capable of identifying the different degrees of smoothness of each scale. We also see that the true parameters fall either inside or close the interquartile range of estimated hyperparameters with the exception of λ f . This indicates that more fine scale observations are generally required to accurately estimate the fine scale correlation length. We finally remark that fig. 8 indicates that both ML and LOO-CV result in similar estimate variability. Nevertheless, we recommend the use of LOO-CV estimation as it was more robust with respect to initial guess of hyperparameters.
To close this section, for Tests 1 and 2, we compare the fine and coarse scale fields estimated from multiscale data and estimated from the corresponding single scale datasets. For this purpose we compute the mean square error (MSE) of the predicted log-conductivity, defined for scale a = {c, f } as The results are summarized in table 3. It can be seen that, for both the fine and coarse scale fields and for both tests 1 and 2, going from single scale data to multiscale data results in a reduction in MSE. This reduction is more pronounced for the coarse scale field (24% for test 1 and 33% for test 2). In fact, we expect fine scale measurements to inform the coarse field estimate more than vice versa due to the indeterminacy of the inverse of the coarsening relation. Similarly, reduction in MSE is more pronounced for Test 2 than for Test 1, which is to be expected as we observe less scale separation between fine and coarse fields for Test 2.

Uncertainty Propagation
We present an application of multivariate GP regression to uncertainty propagation in physical systems. We consider Darcy flow in heterogeneous porous media, where we are interested in predicting spatial hydraulic pressure (h) profiles from multiscale observations of h and of the hydraulic conductivity K, and estimating the uncertainty in these predictions. We assume hydraulic head observations have a fine-scale support volume corresponding to piezometric observations, whereas hydraulic conductivity observations have fine-scale and coarse-scale support volumes corresponding to both laboratory experiments performed on field samples and cross-pumping tests.
In order to employ the Darcy flow model to propagate uncertainty in the spatial hydraulic conductivity distribution, both the hydraulic pressure and conductivity fields must be defined with the same support volume. As h observations are fine-scale, we require constructing a probabilistic model for the fine-scale hydraulic conductivity field K f from multiscale observations of K. We assume there's a data set D s of multiscale observations of the log hydraulic conductivity Y = log K. In particular we assume the data set of Test 1 in section 4.1.
The Darcy flow problem reads is a set of discrete nodes. Our goal is to estimate h f (x) given discrete observations of h f and the multiscale data set D s of observations of Y . We assume that observations h f s of the discrete hydraulic head field are taken as h f s = Hh f + h , where H ∈ R N hs ×N h is a sampling matrix with N hs N h , and h ∼ N (0, σ 2 eh I) is the normally distributed observation error. For this experiment, we take N hs = 20 head observations at observation locations randomly chosen from the set of discretization centroids. Furthermore, we assume an observation error We model the discrete log hydraulic conductivity y as a GP. Given the uncertainty in y and the observations h f s , the minimum mean square error estimator of the hydraulic head field,ĥ (predictive mean), and its covarianceĈ h (predictive covariance), are given whereh and C h are the prior mean and covariance of h f given the uncertainty in Y . We estimateh and C h from 1 × 10 3 Monte Carlo (MC) simulations of the GP y f | D s and the problem eqs. (20) to (22). Results for the profile x 2 = 0.5L are presented in fig. 9 and table 4. fig. 9a shows the mean and variance along the profile x 2 = 0.5L. It can be seen that the variance of h f due to uncertainty in Y is narrow because of conditioning y f on multiscale data. table 4 presents the L 2 -norm of the prior variance along the profile x 2 = 0.5L obtained when y f is conditioned on fine-scale data only and on multiscale data. Consistent with our observations in section 4.1, conditioning on multiscale data results in a reduction in uncertainty (∼ 30%), this time not only for the parameters but for the state itself.
Employing the MC estimate ofh and C h , we apply eqs. (23) and (24) to obtain the predictive mean and covarianceĥ andĈ h conditioned on h f observations. fig. 9b shows the mean and variance of the estimated hydraulic head along the profile x 2 = 0.5L, together with the reference profile. Good agreement of the predictive meanĥ with the reference profile is observed, with the reference profile falling inside the 95% confidence interval indicated by the predictive variance. table 4 shows the L 2 -norm of the conditioned variance when the GP model for y f is conditioned on fine-scale data and on multiscale data. It is again seen that by conditioning on multiscale data the uncertainty in the state is reduced.

Conclusions and Future Work
We have presented a multiscale simple kriging approach for parameter field reconstruction that incorporates measurements with two different support volumes (fine and  Table 4. L 2 -norm of the variance of normalized profile h f (x1, x2 = 0.5L) prior and after conditioning on h f observations. GP model for y f trained first on fine-scale data only, and then on multiscale data. h f (x 1 , x 2 = 0.5L) variance Fine-scale data only Multiscale data Prior 1.300 × 10 −5 9.488 × 10 −6 Conditioned 2.210 × 10 −6 1.840 × 10 −6 coarse scales). Our results show that treating each support scale as components of a bivariate Gaussian process with full bivariate Matérn covariance function is a viable strategy for modeling multiscale observations. We find that the full bivariate Matérn model with hyperparameters estimated from data accurately captures the stationary features of reference fields for scenarios exhibiting moderate to large separation between fine and coarse scales. In particular, we show that the model selection approaches employed in this work, marginal likelihood and leave-one-out cross-validation, both qualitatively identify the difference in smoothness and correlation lengths between scales. Nevertheless, we recommend the use of LOO-CV estimation as it was more robust with respect to initial guess of hyperparameters.
Furthermore, we show that conditioning on multiscale data can be used to reduce uncertainty in stochastic modeling of physical systems with heterogeneous model parameters sampled over multiple scales.
Future work will extend the proposed approach to more than two spatial scales by employing valid multivariate Matérn models [Apanasovich et al., 2012].

A: Nyström Method for Simulating GPs
In this section we present the Nyström method employed in this work for simulating GPs. Nyström methods are used to approximate the factorization of large covariance matrices in terms of the factorization of a smaller covariance matrices [Halko et al., 2011].
Let y(x) ∼ GP(0, C(x, x )). Our goal is to simulate realizations of y, y (k) * at the N nodes X * . Simulating y via factorization (e.g. Cholesky, eigendecomposition, etc.) is computationally prohibitive for large N as factorizations require O(N 3 ) operations and O(N 2 ) storage. Therefore, we aim to compute a low-rank approximation to the covariance matrix of the form C(X * , X * ) ≈ AB −1 A , where B is a rank-M spd matrix with M N .
By Mercer's theorem, the covariance kernel C(x, x ) allows the representation in terms of its eigenpairs {λ k , ψ k } ∞ k=1 , where the eigenpairs satisfy the relation λ i ψ i (x) = Ω C(x, x )ψ(x ) dx (A.2) and the orthogonality relation Ω ψ i (x)ψ j (x) dx = δ ij . We approximate the integral in eq. (A.2) using an M -node quadrature rule over Ω with nodesX k = (x 1 , . . . ,x M ) and uniform weight w (e.g., midpoint rule, MC sampling, etc.), obtaining Furthermore, letC ≡ C(X,X) denote the covariance matrix for the quadrature nodes, with eigendecompositionC =ŨLŨ ,Λ = diag(λ i , . . . ,λ M ). Approximating ψ i (x k ) and λ i in terms ofŨ andΛ, we obtain C(x,x k )Ũ ki , k = 1, . . . , M. (A.4) The relations in eq. (A.4) can be employed to construct an approximate K-L expansion, which then can be used to simulate the field. In this work we follow an alternative procedure, employed in section 4.2 for simulation and in section 4.1 to compute conditional mean and variances. Substituting eq. (A.4) into eq. (A.1), we obtain the following rank-M approximation to the covariance matrix C(X * , X * ), C(Y, Z) ≈ C(Y,X) ŨΛ −1Ũ C(X, Z) = C(Y,X)C −1 C(X, Z), (A.5) where Y and Z are two arbitrary sets of nodes.
For prediction, we employ eq. (A.5) to compute C c p (and C f p ) in eq. (7). For simulation, we note that the covariance matrixC admits the Cholesky factorizationC = LL . Substituting this factorization into eq. (A.5) for Y = Z = X * we obtain the rectangular factorization C(X * , X * ) = LL , L =L −1 C(X, X * ). (A.6) The realizations are simulated as y We employed the centroids of a 128 × 64 (M = 8192) uniform discretization of Ω as quadrature nodes for computing conditional mean and variance in section 4.1 and for MC simulations presented in section 4.2.