Accounting for Local Geological Variability in Sequential Simulations—Concept and Application

Heterogeneity-preserving property models of subsurface regions are commonly constructed by means of sequential simulations. Sequential Gaussian simulation (SGS) and direct sequential simulation (DSS) draw values from a local probability density function that is described by the simple kriging estimate and the local simple kriging variance at unsampled locations. The local simple kriging variance, however, does not necessarily reflect the geological variability being present at subsets of the target domain. In order to address that issue, we propose a new workflow that implements two modified versions of the popular SGS and DSS algorithms. Both modifications, namely, LVM-DSS and LVM-SGS, aim at simulating values by means of introducing a local variance model (LVM). The LVM is a measurement-constrained and geology-driven global representation of the locally observable variance of a property. The proposed modified algorithms construct the local probability density function with the LVM instead of using the simple kriging variance, while still using the simple kriging estimate as the best linear unbiased estimator. In an outcrop analog study, we can demonstrate that the local simple kriging variance in sequential simulations tends to underestimate the locally observed geological variability in the target domain and certainly does not account for the spatial distribution of the geological heterogeneity. The proposed simulation algorithms reproduce the global histogram, the global heterogeneity, and the considered variogram model in the range of ergodic fluctuations. LVM-SGS outperforms the other algorithms regarding the reproduction of the variogram model. While DSS and SGS generate a randomly distributed heterogeneity, the modified algorithms reproduce a geologically reasonable spatial distribution of heterogeneity instead. The new workflow allows for the integration of continuous geological trends into sequential simulations rather than using class-based approaches such as the indicator simulation technique.


Introduction
Drawing conclusions from uncertain data in Earth sciences is rather usual than unusual. Each measurement in geoscientific studies is affected by measurement errors and represents only a subset of the natural variability of geological media. The natural variability is a substantial business-critical controlling factor of different types of subsurface utilization such as mining, hydrocarbon and geothermal exploitation, carbon capture and storage, or nuclear waste disposal. The physical variability of rocks is defined as the complexity or heterogeneity of a system within time and space [1]. Even marginal discrepancies from the predicted property distributions in the subsurface can lead to inaccurate simulations of a quarrie's production potential or a reservoir's recovery and life-time [2,3]. Especially the small-scale variability of rock physical properties makes field-sized predictions still challenging.
Natural heterogeneity and the corresponding property distribution in time and space can be modeled through interpolation, statistical regression, machine learning or stochastic simulation [4,5] by using a number of observations or training data. Due to technical, economic or temporal limitations, geoscientific sampling campaigns practically always end up in scarce data sets within a target domain Ω. Accordingly, estimates of properties often do not account for or misfit the observed geological structures in the field and especially conventional interpolation techniques such as kriging produce smooth transitions at sharp geological boundaries. Moreover, they may fail to reproduce the global statistics appropriately. Conventional interpolations tend to underestimate the presence of values in the upper tail of a distribution and likewise in the lower tail, too [5]. Consequently, major geological heterogeneities, such as faults, major bounding surfaces, or physicochemical anomalies, are very likely not to be reproduced appropriately by a continuous random function (RF) [6].
In contrast to conventional interpolation techniques, stochastic simulations aim to reproduce the variance and the histogram observed in the global data [7,8]. Based on either being constrained or not, stochastic simulations split up into unconditional and conditional simulations [9]. Unconditional Monte Carlo-based simulations reproduce the original histogram without spatial constraints. The realizations produced by those methods, however, are regularly far away from representing the true spatial distribution and constitute "most likely" cases at the best. Conditional simulations, in contrast, aim to reproduce the original property distribution by means of discretely sampled points together with spatial characteristics such as the observed variogram model [10].
One type of conventional simulation algorithms is represented by the sequential Gaussian simulation (SGS) in which the local variability is simulated by sampling the local probability density function (PDF) derived from the local simple kriging variance σ 2 SK . This parameter results from the previously performed interpolation of the standard normally distributed data set [11]. Early field studies have proven the potential of this method to predict rock properties at unknown locations and to assess the uncertainty that can be expected in the area of interest [12][13][14][15]. More recent approaches lead to modifications of the SGS algorithm without the need to transform the original variable into standard normal space. That technique-better known as direct sequential simulations (DSS)-may, for example, sample from the global histogram rather than from the local PDF [9] or perform a quantile-quantile back-transformation into the original variable's space after the simulation. Those approaches can reproduce both the original histogram and experimental semivariogram model as well [10]. The local PDF derived from σ 2 SK , however, mainly reflects the degree of uncertainty induced by the interpolation method itself and does not necessarily reflect the local variability observed on a smaller scale than Ω.
In order to enhance the accuracy of sequential simulations, we propose a new workflow, which incorporates the local variability derived from measurements on a subset of Ω into SGS and DSS under the consideration of measurement errors. The modified SGS and DSS algorithms utilize a global representation of the locally observable variance, named local variance model (LVM), in order to draw a value at an unsampled location. Accordingly, the algorithms are called LVM-SGS and LVM-DSS. Before simulation, an integer programming optimization analysis is performed in order to optimize the robustness of the underlying interpolation function. Instead of sampling from the local PDF, which is generated by means of σ 2 SK , or by solving a global optimization problem, our parametric approach simulates a local PDF based on a measurement-constrained and geology-driven variance extracted from the LVM. The local PDF hereby is simulated with a Box-Muller transform [16].
The method was tested and validated in a case study, which has been conducted in a potential geothermal reservoir formation in southwestern Germany. Therefore, we measured the intrinsic permeability, representing a key parameter in many types of subsurface utilization, on a set of samples taken from an active quarry. Ω is represented by a 3-D outcrop model, which is constructed by means of photogrammetric outcrop wall reconstruction. The model covers a volume of 9000 m 3 . Small-scale variability is derived from rock samples, which are taken from two representative rock cubes. Those are regarded as Ω b and cover a volume of 0.0156 m 3 and 0.008 m 3 , respectively. The rock cubes are taken from the same outcrop, from which the global samples are taken from. Eventually, our approach is compared to the conventional SGS and DSS algorithms and assessed by its ability to reproduce the global variogram model and the geological heterogeneity.

Spatial Variability
In order to reduce the probability of economic failure in mining industries, the concept of the regionalized variable had been developed by Matheron [17] in the 1960s. The regionalized variable is a function that takes a definite value at each point of space. In geological media that regionalized variable often proved to be too complex to be expressed by mathematical functions. A regionalized variable is assumed to show a more or less steady continuity in space accompanied by local fluctuations (Figure 1). In geological media, those fluctuations usually result from the physical variability observed at smaller scales. Lithological and physical variability is subject of numerous geoscientific studies [18][19][20][21] and is commonly termed heterogeneity. In the Oxford Dictionary [22] the word heterogeneity is defined as a Difference or diversity in kind from other things or a Composition from diverse elements or parts; multifarious composition. In most works, this term is used to describe that an object consists of multiple subsets being different to one another in one or more attributes. Li and Reynolds [1] restrict the term to be the variability of a system property in three-dimensional space. Fitch et al. [23] provide a set of methods to quantify heterogeneity within a sample of observations including the coefficient of variation (c v ), where σ is the standard deviation and µ is the arithmetic mean and the Dykstra-Parsons coefficient (c dp ) where x n is the nth percentile of a distribution. The continuity of a regionalized variable is thus dependent on the continuity of the geological media and may or may not provide continuity in a mathematical sense. In this work, we will use the term property for a regionalized variable, the term field for the (quasi-)continuous spatial distribution of a property, and the term target domain Ω for an area of interest. When we mention global and local characteristics, we refer to characteristics of Ω and its subsets Ω b , respectively.

Geostatistical Interpolation
Geostatistical interpolation techniques aim to estimate a value at unsampled locations of a property in Ω and build the base for sequential simulations. The most popular geostatistical interpolation technique is kriging. In the following subsections, we will briefly describe the theory behind kriging and focus on its variety simple kriging (SK). Moreover, we will discuss practical computational aspects such as neighborhoods.

Spatial Neighborhood
As the system of linear equations for geostatistical estimations might grow very large, those algorithms require subset-sampling in order to perform reasonably. Therefore, a 3-D search ellipsoid can be used to find the neighbors of a point in a mesh [24]. This ellipsoid can be defined by six properties: azimuth α; dip β; plunge γ together with the radius in X r x , Y r y , and Z direction r z of the ellipsoid. α, β, and γ define the ellipsoid's clockwise rotation around the Z, X, and Y axes in this exact order. Accordingly, the rotation matrix T can be defined as After translating the mesh such that x x = x y = x z = 0 and rotating it according to Equation (3), Equation (4) can be used to determine, whether a point x with the transformed coordinates x x , x y and x z is located inside or on the boundary of the search ellipsoid (≤1) or not (>1).

Variography
The variographic analysis is a crucial prerequisite for numerous geostatistical interpolation techniques. Hereby, the experimental semivariogram represents the cumulative dissimilarity of a discrete set of point-pairs with n c representing the count of point-pairs within the distance classes h of identical distance increments (Equation (5)).
The continuous counterpart, represented by the theoretical semivariogram γ theo , is an approximation of the experimental semivariogram assuming z(x) to be a stationary random field [25]. γ theo is used to fit the experimental variogram. The spherical variogram model γ sph with a nugget effect is a popular nested model used to fit the experimental semivariogram [26,27], which is calculated by with n being the nugget, b the sill and a the range [6]. The variogram model represents a covariance function c with the relationship γ(h) theo = c(0) − c(h), where c is a positive definite, even function and c(0) = n + b in case of a spherical variogram model with nugget effect. Semivariograms can be used to quantify the spatial or time correlation of a random variable [27][28][29]. c and γ theo are input variables for geostatistical interpolation algorithms.

Simple Kriging
Kriging is a commonly used stochastic technique to interpolate geological rock properties in space and time [30]. The kriging estimator is the best linear unbiased estimator (BLUE) of a property as it minimizes the error variance. It incorporates the covariance structure of the globally sampled values into the weights for predicting the value z(x 0 ) at an unsampled location x 0 [31]. Therefore, z(x 0 ) is calculated by weighting the values of the sampled locations and building a linear combination of those what gives where w k is the weight of the sampled point x k with the value z(x k ). The kriging types primarily differ by their derivation of the weight vector. For all kriging systems, a system of linear equations must be solved as it is outlined in the following paragraphs, in which we will consider the simple kriging (SK) technique [32] and expand it by the integration of a locally varying mean [33]. Therefore, we modify Equation (7) into in which the known stationary mean µ has been added [6]. While SK assumes that µ is globally constant and known, SK with locally varying mean assumes µ to be constant only in the neighborhood of x 0 . In order to obtain the SK weights, a system of n linear equations must be solved in which n stands for the number of considered neighbors. This system of equations is defined as which corresponds to with c as covariance function and x n as point with known value [25]. In SK each interpolated point provides a simple kriging variance σ 2 SK [5], which we can calculate by means of the formula The quality of a kriging interpolation is dependent on the variogram model and its goodness of fit to the experimental semivariogram.

Consideration of Measurement Error Variance
We already saw that kriging induces a local interpolation error by itself, namely, σ 2 SK . There are, however, also other components which bias the interpolation result. Besides σ 2 SK , the local and unknown variability of z(x) in Ω b as well as the measurement error variance σ 2 m might play an important role ( Figure 2). Integrating σ 2 m into an interpolation can be achieved by estimating the measurement error precision σ m with a variance of σ 2 m and incorporating it into the kriging system of linear equations giving In contrast to the conventional formula, σ 2 m with regard to the considered known value at x k is added in the diagonal of the matrix [25]. This accounts for the heteroscedastic nature of geological parameters as they commonly show a higher variability for high values and a lower variability for low values.

Sequential Simulation
In contrast to geostatistical interpolation techniques, sequential simulations aim to reproduce the global statistics in form of the considered variogram model and the global histogram. Therefore, in order to account for the spatial heterogeneity of a rock property, the sequential Gaussian simulation (SGS) and the direct sequential simulation (DSS) algorithm can be utilized for univariate simulation. SGS is based on the multi-Gaussian approach [33], which assumes that the kriging error is standard normally distributed with µ = 0 and σ 2 SK = 1. This requires that each one-point cumulative density function (CDF) of any linear combination of the RV is normally distributed, that all subsets of the RF are multivariate normal, that the two-point distribution is normal and that all conditional distributions of subsets of the RF are normal [33]. If the RF fulfills the requirements, then the simple kriging estimate and variance characterize the posterior cumulative CDF under consideration of the normal score variogram model. Thus, we need to transform the original distribution's CDF into standard normal space for SGS. In order to transform any point in the CDF (F(Z(u))) of any random variable Z(u) to a random function Y(u) and vice versa the following equation can be applied, where G −1 is the inverse Gaussian CDF of Y(u), which is also named quantile function [34], and φ is the inverse Gaussian CDF of F(Z(u)). Thus, z and y correspond to the same probabilities. For each previously interpolated point x j now a random value of the normal distribution N µ SK , σ 2 SK , whose PDF defines as is drawn as z(x 0 ) using the Box-Muller transform [35]. We can perform this transform by applying the equation with u 1 and u 2 as random numbers ∈ [0, 1], σ as the standard deviation, and µ as the mean of the original distribution. The simulation is eventually back-transformed into the original space using a quantile-quantile back-transformation mapping technique. The reproduction of the covariance model, however, does not require the multi-Gaussian approach as long as the estimate and variance are derived from the SK estimation [9,10]. Thus, the conditional distribution type, which is chosen in order to simulate the variability at each point, does not necessarily need to be Gaussian. With this in mind, it is evident that a normal score transform is not needed before performing a sequential simulation. This results in the DSS approach, which commonly samples from the global PDF by determining the sampling interval from the local PDF [9].

Cross-Validation
In order to assess the quality of a realization, models, which are constructed by means of interpolation or simulation techniques, should be validated. Commonly, interpolations are validated by cross-validation. This technique is usually performed by using point removal procedures called leave-p-out cross-validation (LpO CV). For the LpO CV, p randomly selected samples are removed from the input data set of size n with 0 < p < n and the interpolation or simulation is performed without these samples [36]. As measures of goodness of fit, the mean-square error (MSE, Equation (16)), the root-mean-square error (RMSE, Equation (17)), and the mean-absolute error (MAE, Equation (18)) of the realization can be calculated as and whereẑ(x k ) are the simulated points. While Willmott et al. [37] question the status of the triangle inequality for the RMSE, which is required for a distance function metric, Chai and Draxler [38] show that the RMSE in fact fulfills this condition. Thus, if the model errors follow a normal distribution, the RMSE is to favor over the MAE [38].

Ergodic Fluctuations
The minimum requirement for geostatistical simulations is their ability to reproduce the original data, the global summary statistics and the global variogram model [8,39]. Erdogic fluctuations refer to the discrepancy between the model parameters and the realizations' statistics [6]. In the case of the variogram model, the discrepancy of a realization to the variogram model is related to the limitation of the integrated constraints to a limited neighborhood. This, in fact, leads to higher errors at far ranges within the simulation. In this study, we quantified the ergodic fluctuation of a realization by estimating the average MSE between the experimental semivariogram and the variogram model. If a realization's discrepancy among the experimental variogram and variogram model does not exceed the original values discrepancy, the variogram reproduction is said to be within the range of ergodic fluctuations.

Sequential Simulation using a Local Variance Model
In this section, we will describe how the SGS and DSS algorithms need to be modified in order to sample from a local variance model (LVM). The LVM can be described as a global representation of the locally observable variance σ 2 LV M in one mesh cell. Thus, the LVM can be referred to as the local geological heterogeneity. The LVM is constructed using a mapping technique in which the value of the mapped variances is constrained by a set of measurements. Those are intended to represent the small-scale variability present at the mapped position. Subsequently, the variance is interpolated onto Ω. The basic concept of interpolating a distribution in space is illustrated in Figure 3c. The sequential simulations are performed on the nodes of Ω using a modification of the SGS and DSS algorithms, namely, the LVM-SGS and LVM-DSS. Our basic idea is that, if and only if the geological heterogeneity is exceeding σ 2 SK at x k , we will sample from the LVM-constructed PDF instead of from the kriging-derived PDF. Otherwise, if the interpolation error is greater than the expectable geological heterogeneity, we will sample from the kriging-derived PDF. The generalized algorithm is displayed in Algorithm 1. All analyses have been conducted with the open-source software GeoReVi [41] in which the new algorithms have been implemented as extensions in the C# programming language (Appendix A).  (13), Back-transform the simulated values into the original space

Case Study
In order to test and evaluate the new workflow with the modified algorithms, we conducted an outcrop analogue study in a quarry in Germany. In the following subsections, we will outline the object of investigation, the sampling strategy and the modeling techniques used to implement the LVM-SGS and LVM-DSS algorithms. We decided to use the intrinsic permeability k for the implementation as that property plays a critical role in numerous types of subsurface utilization-especially with regard to subsurface reservoirs.

Object of Investigation
An actively quarried sandstone outcrop (long. 7.647546, lat. 49.523821) in Obersulzbach, which is located in the Saar-Nahe basin in southwestern Germany, has been selected as object of investigation (Figure 3a). The outcrop exposes the Disibodenberg Formation of the innervariscan Rotliegend Group, which constitutes a deeply buried [42] potential hydrothermal reservoir unit [43] in the northern Upper Rhine Graben. The Disibodenberg Formation in the quarry is composed of two Bouma sequences ( Figure 3b) from a lacustrine delta, which deposited during Permian times. There were two selection criteria being decisive for selecting the quarry. On the one hand side, the sedimentary beds are ≥2 m thick and laterally continuous. Moreover, the outcrop is actively mined, which reduces the impact of recent weathering onto the permeability. Moreover, it was possible to extract both rock samples from the outcrop wall as well as oriented rock cubes from different representative lithofacies types in order to conduct multi-scale three-dimensional investigations. The outcrop measures 50 × 15 × 10 m and thus owns the size of a typical cell in common static and dynamic reservoir models (see, e.g., in [44]).

Sampling Strategy
Numerous studies showed that the physical variability in geological media must be integrated as a function of measurement volume, also known as the representative elementary volume (REV) [45].
The REV denotes a volume, at which a representative amount of heterogeneity is captured by one measurement [46] minimizing the smaller-scale fluctuations. Therefore, a multi-scale approach based on the concept of the REV has been implemented. Accordingly, 39 cylindrical rock samples with diameters and lengths of four centimeters were extracted from the outcrop wall. The samples were taken from six 1-D profiles covering the entire quarry area (Figure 4a). More information regarding the sample positions and orientations can be found in Linsel [41]. Those samples were used for the global field simulations. The quarry contains sequences from a prodelta mouthbar deposited as turbiditic densites. The sequences graduate from a high-energetic depositional environment at the base to a low-energetic environment at the top as the flow velocity is steadily declining. The sequences consist of heterogeneous, intraclast-rich sandstones at the base and of trough cross-bedded, ripple cross-bedded and homogeneous sandstones at the top. Consequently, the sequences can be declared as Bouma sequences containing the Bouma A to Bouma E intervals in a fluvial-dominated lacustrine-deltaic depositional environment [47]. Based on that, we assumed that the variability within one Bouma sequence is highest at the base and lowest at the top (Figure 3b). Accordingly, two rock cubes of 0.2 × 0.2 × 0.2 m (OSB2_c) and 0.25 × 0.25 × 0.25 m (OSB1_c) were taken-one from the top (Bouma E) and one from the base (Bouma A) of one sequence-in order to capture both the highest and the lowest variability. The locations of the cubes within the quarry and the strata are shown Figure 3a,b.
We selected two types of lithofacies: OSB1_c, a discontinuously cross-bedded, intraclast-rich lithofacies type and OSB2_c, a homogeneous lithofacies type without macroscopically observable internal bounding surfaces. In total, 79 rock cylinders were extracted from rock cube OSB1_c and 29 from OSB2_c. More information regarding the sampling process can be found in Linsel et al. [40]. Those samples were used for constraining the LVM.

Laboratory Measurements
The cylinder samples were cut, oven-dried at 105 • C and measured in the laboratory for determining the intrinsic gas permeability k at unsaturated conditions. k can be considered one of the key parameters of geothermal reservoir rocks with regard to hydrothermal systems in porous aquifers [48]. k was measured with the Hassler cell Darmstadt permeameter. The device's functionality is described in detail in Filomena et al. [49]. The permeability is provided in the industry-standard unit millidarcy (mD), where 1 mD corresponds to 9.869 × 10 −16 m 2 . The permeability measurement provides an error variance between 0 and 0.15 mD 2 in the range of the observed values [50].

Mesh Generation
In order to construct Ω, the outcrop wall is modeled using a photogrammetric representation that was downsampled into 40 × 20 faces and subsequently interpolated using Shepard's p-value IDW interpolation, which we can write as where d is the Euclidean distance between the the point with the known value x k and the point with the unknown value x 0 and p is an exponent factor to influence the weights non-linearly. IDW has been applied with a short search radius of five meters and a power parameter of four. The interpolation result has an RMSE of 0.024 m, which can be considered low for the surface interpolation. The resulting outcrop surface is used as a bounding surface for a hexahedral mesh, which represents Ω, that is composed of 75,240 cells ( Figure 4b, Table 1). The rock cubes, which represent Ω b , are constructed by an orthogonal, hexahedral mesh containing 25,230 (OSB2_c) and 64,000 cells (OSB1_c), respectively. The volume of an average cell of the outcrop mesh is roughly eight times the volume of OSB1_c and 15 times the volume of OSB2_c (Table 1).  The variance σ 2 c derived from the measurements conducted on the samples from the rock cubes is assumed to represent the variance σ 2 Ω b that can be expected in one cell of the outcrop mesh so that with σ 2 LV M being the local sample variance, which we can calculate by means of the formula where n is the total number of samples, µ is the mean and x i is the sample at the ith location.

Spatial Variability
The variogram analysis reveals a range of 0.3 m and 0.2 m for the rock cube samples OSB1_c and OSB2_c, respectively, and a range of 18 m for the outcrop samples (Figure 5a,d,g). The sill is slightly higher in the outcrop region as it is in the rock cubes. Moreover, the outcrop samples show a weak nugget effect. Generally, a scale effect can be observed in which the variance increases with the considered volume. This effect is also present in the descriptive statistics (Figure 5c,f,i).
The measurements from the outcrop region show a c v of 0.28 and a c dp of 0.31. The histogram indicates a normal distribution of k ranging from 0.7 mD to 4.6 mD (Figure 5b). A two-sided Kolmogorov-Smirnov test [51], which is based on an implementation of Simard and Ecuyer [52], confirmed the hypothesis that all samples come from a normal distribution. The application of Tukey's outlier detection method [53] reveals no statistical outliers in the sample. By applying the classification scheme of Corbett and Jensen [54], the sample can be classified as being very homogeneous.
The local histogram of k from OSB1_c shows a bimodality in the distrubtion ranging from 0.7 to 3.9 mD (Figure 5e). OSB2_c's histogram shows a unimodal range from 0.8 to 1.5 mD (Figure 5h). Again, no statistical outliers can be detected. The local variability of OSB1_c is significantly higher than that of OSB2_c. k of OSB1_c provides a c v of 0.3 and a c dp of 0.4 while measurements from OSB2_c show values of 0.2 for c v and 0.18 for c dp . c v and c dp of OSB1_c tend to cover the variability of the global data. This result is in good agreement with the REV theory from Nordahl and Ringrose [45]. Both rock cubes can be classified being very homogeneous as well. Thus, we can observe a significant small-scale variability. The bedding structures in OSB1_c are well preserved in the permeability field of the k interpolation, which gradually increases from low values between 0.7 and 2 mD in the lower beds to higher values between 2 and 4 mD in the upper beds (Figure 6a). In OSB2_c the trend is running diagonally through the rock cube ( Figure 6b); however, no macroscopic bounding surfaces are visible, which could have had a control on the field of k here. It should be noted, however, that the range of k is significantly smaller here compared to OSB1_c.

Constructing the LVM
The LVM is constructed by means of a 3-D architectural element mapping of both Bouma sequences in the quarry. The base and the top of the sequences are mapped which are being used to constrain the LVM by the locally observable variance σ 2 LV M . The exploratory data analysis reveals that the variance of k in OSB1_c is five times larger than that of OSB2_c. This is in accordance with the sedimentological mapping, which indicates a higher heterogeneity at the base of the Bouma sequence.
It is assumed that OSB1_c represents the most heterogeneous and OSB2_c the most homogeneous lithofacies type in the Bouma sequences as it is illustrated in Figure 3b. Accordingly, the positions of those lithofacies types are mapped throughout the quarry area and parameterized with σ 2 LV M , which has been determined by the k measurements of OSB1_c and OSB2_c. Thus, we use σ 2 LV M = 0.43 for mapping the base boundaries of the sequences throughout the outcrop area. Likewise, σ 2 LV M = 0.07 is used as a local variance for the topmost boundary of the single sequences. The mapping locations of LV M are shown in Figure 7a. The points mapped onto the outcrop model are subsequently interpolated onto Ω by using a SK-based interpolation procedure for parametric PDFs (Figure 3c). The interpolation is conducted using 5 neighbors, a range of five meters, a sill of 0.005, a nugget of 0 and a plunge of 10 • as the strata gently dip towards south. Figure 7b shows the constructed LVM which is being used by the sequential simulation algorithms. It should be noted that we have a decent offset in the LVM in the area of the central fault zones.

Optimizing the BLUE for Sequential Simulation
Prior to sequential simulation, the optimal SK conditions with regard to the integrated measurement error variance σ 2 m and the selected neighborhood are determined. Therefore, a simple integer programming optimization is performed using varying measurement error variances (0.0 mD 2 ≤ σ 2 m ≤ 0.15 mD 2 ) and a varying number of neighbors (10 ≤ n n ≤ 20) as inequality constraints. We can express the optimization problem as min σ 2 m ∈R, n n ∈N SK (σ 2 m , n n ) (22) subject to 0 ≤ σ 2 m ≤ 0.15 10 ≤ n n ≤ 20, in which the SK error SK in form of the RMSE and MAE must be minimized. The response surface of the numerical optimization process indicates that the SK error is generally declining when σ 2 m is increasing. The lowest errors are produced with an n n of 10, 11, and 20. This sensitivity of the SK error on the number of neighboring points is not unusual. The numerical optimization reveals that the optimal conditions for SK are met at n n = 20 and σ 2 m = 0.15 which yields a RMSE of 0.708 mD (Figure 8). The interpolation error can be reduced by 16.5% for the RMSE and by 18.5% for the MAE. The final SK realization and the spatial distribution of σ 2 SK for that exact model is illustrated in Figure 9. It should be noted that the spatial distribution of σ 2 SK in a sequential simulation is different as previously simulated locations are considered as well. The interpolation error RMSE is minimized using the inequality constraints given in Equation (22). (a) RMSE response surface with regard to the incorporated measurement error variance σ 2 m and the maximum number of neighbors n n using a leave-one-out cross-validation. (b) Cross sections through the response surface of (a).
The final modeling variables for the sequential simulations are given in Table 2  The statistical and spatial characteristics of σ 2 SK and σ 2 LV M differ tremendously. σ 2 SK is unimodally distributed, whereas σ 2 LV M provides a bimodal distribution (Figure 10a). It is evident that σ 2 SK covers the total range of the considered covariance model while σ 2 LV M 's range is more limited. The probability of simulating variances between 0.2 and 0.43 mD 2 is higher when sampling from the LVM instead of the local SK variance (Figure 10b). The median between σ 2 LV M and σ 2 SK differs by ≈ 0.08 mD 2 , which indicates that the variability simulated in a realization of conventional sequential simulation algorithms is systematically underestimated.
With regard to the variogram model, σ 2 LV M has a range of 5 m and a sill of 0.36 mD 2 , and σ 2 SK has a range of 0.3 m and a sill of 0.1 mD 2 . Thus, σ 2 SK seems to be spatially uncorrelated and random. However, the grade of variability in the eastern part of the outcrop is slightly higher than in the western part. Therefore, in contrast to σ 2 LV M , σ 2 SK obviously does not provide the simulation algorithm with a spatial trend when simulating the local variability.

Model Validation
All algorithms reproduce the considered variogram model within the range of ergodic fluctuations after back-transformation (Figure 11a-d). The quality of variogram reproduction has been evaluated by calculating the average mean square error MSE of all realizations between the experimental variogram and the variogram model. The best reproduction is produced by the LVM-DSS and LVM-SGS algorithms, while the latter one provides the lowest degree of ergodic fluctuations with MSE = 0.066 mD 2 . All realizations reproduce short-range dissimilarities well but slightly underestimate the dissimilarity at medium ranges. DSS and SGS tend to gentle underestimation at far ranges which is a drawback of limited neighborhoods. This effect, however, is less expressed in the LVM-based algorithms. For both types of sequential simulation, the LVM-based algorithm outperforms the conventional conditional simulation approaches.

Visual Outputs
It is evident that all simulation algorithms provide visually comparable results ( Figure 12). It should be noted that the quadrilaterals of the 3-D models are subdivided using the Catmull-Clark scheme [55] for visualization. Within this scheme, a new point in a quadrilateral is calculated by with x k+1 j as the new point at subdivision step k + 1 in the center of the element j with n vertices at the subdivision step k. This technique smooths the observable patterns in the models. There is an obvious trend in all realizations, which indicates that the highest values are located in the eastern part of the quarry and the lowest values in the western part. Having in mind that the applied algorithms are conditional, this trend is in well accordance with the constraints as given by the global measurements, which also provide the highest values in the eastern part of the quarry and the lowest values in the western part (Figure 4a). The trend is most clearly depicted in the DSS and LVM-DSS realizations ( Figure 12). SGS and DSS tend to construct homogeneous regions more likely than their LVM equivalents. Thus, those algorithms might indicate a homogeneity, which is likely not present in the strata. Moreover, the heterogeneity of the LVM equivalents is more realistically oriented along the bounding surfaces in the quarry than the models produced by the conventional algorithms.

Discussion
In this study, we present a workflow that accounts for the locally observable geological variability in modified versions of conventional sequential simulation algorithms. Our approaches produce similar outputs as the conventional algorithms and reproduce the global variance model together with the global summary statistics, which are important criteria for the validity of a statistical simulation [8,10,39]. Our results are confirming the concept of the REV [45], in which the complexity of a continuous random variable is increasing with reducing the scale of observation. Moreover, we can confirm that σ 2 SK constitutes no measure for the local estimation accuracy [56] as it is only reflecting the spatial configuration of the constraining data points being simultaneously independent on the constraints' values [6]. There are, however, two points which must be raised in order to discuss the benefits as well as the drawbacks of our approach.

Construction of the LVM
The main source of errors in the proposed workflow is based on the construction of the LVM. The LVM has been derived by an integrated approach of measuring the local variability in the most homogeneous and most heterogeneous lithofacies types in the sedimentary succession. The statistical analysis revealed that this assumption proved to be true as the heterogeneity measures in OSB1_c indicate a way higher variability as is present in OSB2_c. This, in fact, is building the basis for this study. The variance has been assumed to be constant at the base and at the top of a Bouma sequence. This assumption is limited by the number of samples taken within this case study. By constructing the LVM with an SK interpolation, we assume that the variance in one sedimentary Bouma sequence is continuous in a mathematical sense. This assumption might be proved to be too simple in future studies. In order to validate those results, more local samples would be necessary to constrain the LVM. This is a drawback in comparison to conventional SGS and DSS algorithms as those are not dependent on estimating a global variability model. Figure 13a,b illustrates the relationship between the LVM and a DSS realization (a) and an LVM-DSS realization, respectively (b). Although the overall trend remains identical among both types of simulation, the spatial distribution of local variability is uncorrelated and inherently different. In the DSS realization, the heterogeneity within the region is randomly distributed. The most heterogeneous areas in the LVM-DSS realization reside in the light areas-in which σ 2 LV M is high-whereas the most homogeneous regions reside in the dark ones-where σ 2 LV M is low. As the spatial distribution of σ 2 SK is primarily dependent on the distance to the constraining neighbors, the SGS and DSS algorithms, in contrast to their LVM-based modifications, cannot account for a realistic spatial distribution of the local geological variability. This observation is conceptually illustrated in Figure 13c, which shows the spatial relationship between σ 2 SK and σ 2 LV M , as implied by the results of our study. It is evident that the conventional algorithms underestimate the local geological variability in close ranges to conditional data. It is also evident that σ 2 SK systematically underestimates the natural variability present in the geological medium, which is investigated in this study (Figure 10a). Therefore, SGS and DSS might not be able reproduce the total geological variability as shown in this study, which is an advantage of the proposed algorithms instead. (c) Conceptual illustration showing the spatial distribution of the constraining measurements k(x j ) and the spatial relationship between the simple kriging estimate µ SK with the measurement error m and the two parameters used to simulate k in this study namely σ 2 SK and σ 2 LV M . Pr stands for the probability of k under the condition that k belongs to the Gaussian distribution described by µ SK together with either σ 2 SK or σ 2 LV M .

Conclusions
In this study, we propose a new workflow, which incorporates the locally observed variability σ 2 LV M into sequential simulations. We could demonstrate that the local simple kriging variance σ 2 SK differs from σ 2 LV M in local volumes of the target region. Therefore, the DSS and SGS algorithms have been modified by the replacement of σ 2 SK through the measurement-derived σ 2 LV M within one mesh cell. This replacement has been done if and only if σ 2 LV M ≥ σ 2 SK . The LVM has been constructed by means of geological mapping and the assumption that the variability is highest in the most heterogeneous lithology and lowest in the least heterogeneous lithology in a Bouma sequence. The proposed approach can be used in any type of spatial property simulation but is especially tailored for geological media.
The LVM-DSS and LVM-SGS approaches reproduce the observed variability in the sedimentary succession adequately yet reproducing the minimum required statistical measures of a valid simulation including the global histogram, the global heterogeneity, and the variogram model. Moreover, in contrast to their conventional representatives, the LVM-based algorithms account for the spatial distribution of the expected local variance adequately. Once the LVM is derived, it may be integrated into other geostatistical simulation algorithms such as the turning bands method [57][58][59][60].
From our results we conclude the following.
1. The distance metrics RMSE and MAE in spatial interpolations can be optimized with regard to the measurement error variance and the optimal neighborhood. 2. Geological samples always represent a small subset of the local variability, which should be accounted for by high-resolution sampling at a random basis at the least. 3. The simple kriging variance does not necessarily account for the magnitude of local variability in geological media and definitely does not account for its spatial distribution. 4. The fact that the local simple kriging variance does not reflect a geological trend might lead to unforeseen problems when using sequential simulation-derived models as a basis for subsurface utilization processes because the full geological heterogeneity might not have been taken into account properly. 5. By introducing a measurement-constrained, geology-driven local variance model, the spatial distribution of the variance that is expected in the investigated quarry can be integrated into sequential simulations. This allows to simulate the geological variability, which might be greater than the simulated variability in conventional sequential simulation algorithms.
Future research should focus on comparing σ 2 SK and σ 2 LV M under the consideration of other physicochemical properties, other geological settings, and other scales. This might require adapting the assumptions on the spatial continuity of the variability which should, however, always be based on reliable geological analyses. Funding: A.L. has received financial support by a PhD scholarship from the Friedrich-Ebert-Stiftung, Germany, which is gratefully acknowledged.