A comparison of approaches to include outcrop information in overburden thickness estimation

Abstract

The estimation of overburden sediment thickness is important in hydrogeology, geotechnics and geophysics. Usually, thickness is known precisely at a few sparse borehole data. To improve precision of estimation, one useful complementary information is the known position of outcrops. One intuitive approach is to code the outcrops as zero thickness data. A problem with this approach is that the outcrops are preferentially observed compared to other thickness information. This introduces a strong bias in the thickness estimation that kriging is not able to remove. We consider a new approach to incorporate point or surface outcrop information based on the use of a non-stationary covariance model in kriging. The non-stationary model is defined so as to restrict the distance of influence of the outcrops. Within this distance of influence, covariance parameters are assumed simple regular functions of the distance to the nearest outcrop. Outside the distance of influence of the outcrops, the thickness covariance is assumed stationary. The distance of influence is obtained thru a cross-validation. Compared to kriging based on a stationary model with or without zero thickness at outcrop locations, the non-stationary model provides more precise estimation, especially at points close to an outcrop. Moreover, the thickness map obtained with the non-stationary covariance model is more realistic since it forces the estimates to zero close to outcrops without the bias incurred when outcrops are simply treated as zero thickness in a stationary model.

Appendix 1

The computation of the non-stationary covariance is illustrated along a profile in 2D. Three borehole data, one outcrop point and one estimation point are respectively located at x ₁(100,0), x ₂(600,0), x ₃(1000,0), x _out(500,0) and x ₀(450,0). The points x ₂, x _out and x ₀ are within the distance of influence of the outcrop here selected as a _out=120 m. Following Eqs. 5–8, the covariance parameters on these points can be determined. Table 2 gives the resulting covariance parameters associated to each point.

Table 2 Non-stationary covariance model parameters for each point

The same stationary and isotropic covariance model as for the case study is used (exponential model with C0 _S  = 9.5 m², C = 78 m² and a _S  = 6686 m). Then data-to-data and data-to-estimation point covariance are calculated following Eqs. 3 and 4. In the equation, x _i and x _j represent the 2 × 1 vector of coordinates at points i and j. In the model, the NS covariance at point x _i is considered isotropic. Hence, the kernel matrix $\Upsigma_{i}$ is a diagonal matrix with a ² _i repeated along the diagonal, and similarly for $\Upsigma_{j}$ (with a _i obtained from Table 2). As an example, for the pair of points x ₂ − x _out, one computes:

$$ \begin{aligned} Q_{2,{\rm out}}&=\left[100, 0\right] \left(\frac{ \left[ \begin{array}{ll} 6246.2^2 & 0 \\ 0 & 6246.2^2 \end{array} \right] +\left[\begin{array}{ll} 120^2 & 0 \\ 0 & 120^2 \end{array} \right]}{2} \right)^{-1}\left[ \begin{array}{l} 100 \\ 0 \end{array} \right] \\ &= 5.1243 \times 10^{-4} \end{aligned} $$

(13)

Hence, one gets using Eq. 3:

$$ \begin{aligned} R^{\rm NS}_{2,{\rm out}}&=6246.2 \times 120 \times (5.1243 \times 10^{-8}) \times \exp(-\sqrt{5.1243 \times 10^{-4}} \\ &= 0.0375 \end{aligned} $$

(14)

Finally, as the point variance was assumed stationary at 78 m²

$$ C^{\rm NS}_{2,{\rm out}}=0.0375*\sqrt{(78)(78)}=2.9289 \hbox { m}^2 $$

The resulting data-to-data and data-to-estimation point covariances are listed in Table 3 for K-NS and in Table 4 for K-SO. For that particular example, the outcrop kriging weight for K-NS is less than with K-SO. However, as the estimation point gets closer to the outcrop, the weight assigned to the outcrop point increases (see Fig. 8).

Table 3 Non-stationary covariances for K-NS and kriging weights

Table 4 Staionary covariances for K-SO and kriging weights

Fig. 8

Kriging weight assigned to outcrop as a function of the estimation point for K-SO and K-NS

The model is more general than the above example suggests as the kernel matrix $\Upsigma_i$ can be any symmetric positive-definite matrix instead of a diagonal matrix. This would enable to model spatially varying orientation and anisotropy ratio of the ellipses of ranges (i.e. spatially varying geometric anisotropy). As an example, at point x _i with correlation ranges of 1,000 and 200 m along azimuths 30° and 120° respectively, the associated kernel matrix is:

$$ \begin{aligned} \Upsigma_i &= U\Uplambda U^T\\ &= \left[\begin{array}{ll} \cos(30) & \sin(30) \\ -\sin(30) &\cos(30) \end{array} \right] \left[\begin{array}{ll} 200^2 & 0 \\ 0 & 1000^2 \end{array} \right] \left[\begin{array}{ll} \cos(30) & -\sin(30) \\ \sin(30) &\cos(30) \end{array} \right] \\ &= \left[\begin{array}{ll} 280 000 & 415 692 \\ 415 692 & 760 000 \end{array} \right] \end{aligned} $$

(15)

where U is a rotation matrix.

Similarly, the variance can be modeled as a spatially varying function like here for the range and the nugget effect. Rivest and Marcotte (2012) shows an example where the variance of a contaminant is described by a decreasing function of the distance to the source of the contaminant.