The aim of this white paper is to illustrate a workflow to build a model for a complex geological environment from boreholes data and outcrops. The approach is illustrated on a synthetic fluvio-glacial aquifer, in a case study inspired by a paper published by Comunian et al. (2011).
In our case study, the data set comprises 120 boreholes with four indicators (stratigraphic unit, lithology, dip and azimuth). As suggested by locally available outcrops, the geology is expected to present different kinds of structures depending on the stratigraphic unit. Therefore, building a geological model at once would probably not be satisfactory, whatever the method used. The aim of the simulation is to construct 3D models, as realistic as possible, based on the available data and the expected geological structure.
The general workflow proposed for this modeling consists of the following nested steps: (1) Delineate the stratigraphic units (layers); (2) Simulate the lithology (facies) within each layer and finally (3) Simulate physical properties such as permeability and porosity within each facies.
At each step, the most appropriate geostatistical method must be chosen while considering all the relevant data. For some steps, conventional geostatistical methods such as kriging or SGS (Sequential Gaussian Simulation) can be used.
We will focus on showing how DeeSse, a Direct Sampling Multiple Point Statistics (MPS) tool can be used at various steps of the workflow. In classical MPS implementations, the multiple-points statistics characterizing the training image (TI) are summarized into a tree or a list which is used later on, during the simulation process, to compute and sample a local probability density function. With DeeSse this initial step is short-cut: the training image is sampled directly during the simulation stage. This results in an improved efficiency, be it in terms of memory usage or computation times.
A unique feature of DeeSseis that it can be used for modeling categorical but also continuous properties. In a conventional way of doing, the lithology within each stratigraphic unit will be simulated using DeeSse(Step 2). In this step, it will be demonstrated how DeeSsecan take non-stationary orientations into account. In a more unusual mode, DeeSsewill also be used to simulate the elevation of the stratigraphic units’ bottom surfaces (step 1).
The paper first describes the data set and then focus on the modeling approach, detailing each step.
2. The data set
The data set is composed of borehole data and outcrops. The main features of each type of data are detailed in the next subsections.
2.1. Borehole data
The data set contains 120 boreholes. Each borehole was sampled every 20 centimeters along the vertical axis (z) between the topographic surface and the level z = -10 m. At each sample point, 5 variables were recorded:
- well_index = a unique identification number for each well.
- layer = identification of stratigraphic unit taking values from 1 to 4 and labeled from L1 to L4, from the bottom to the surface. This property has been interpreted along the boreholes.
- facies = identification of the lithologies taking values from 0 to 4 and coded by facies labeled F0 to F4.
- dip and azimuth angles in degrees: the structural orientation of the sedimentary deposits has been measured locally along the boreholes, but only in the stratigraphic units 2 and 3.
The spatial representation of the 4 last variables are given in Figure 1.
From the ‘layer’ property of the previous point set, the elevation of the bottom of the stratigraphic units L2 to L4 are retrieved and stored as three continuous properties in an additional 2D point set. Note that missing values can occur at locations where an interface has not been observed, most probably when a stratigraphic unit has been eroded.
An analysis of the borehole data shows that the geology presents different types of structures depending on the stratigraphic unit:
- In unit L1 (bottom): there are only a few thin horizontal structures (facies F1, F2, F3) within the matrix of facies F0. We may consider that these fine structures will not be important for flow simulations and completely neglect them assuming that L1 is homogeneous (filled exclusively with facies F0, see Figure 2).
- In unit L2 and L3: there are some cross stratifications made of thin structures of facies F2 and F3 within a matrix of facies F1. The orientations of the cross stratifications are turning, they are usually sub horizontal at the base of the stratigraphic unit and inclined on top.
- The unit L4 is homogeneous, it contains the facies F4 everywhere.
2.2. Outcrop observations
In addition to the borehole data, one outcrop is available, coming from a local quarry (Figures 2). This outcrop reveals some typical structures displaying cross stratifications, and often the alternation of three facies: a matrix of F1 and then an alternation of F2 and F3 forming some inclined structures.
3. Modeling strategy
To model such structure, the proposed workflow includes the following steps:
- Model the stratigraphic units
- Within each stratigraphic unit, model the lithology (facies)
- Within each lithology, model physical properties such as permeability and porosity.
Those steps are detailed hereafter.
Note that the workflow presented here represents only one realization. To obtain multiple realizations, it has to be put in a loop, generating alternative models for each realization at each of the three steps.
3.1. Step 1 – Modeling the stratigraphic units
In this step, the boundary surfaces of the various stratigraphic units are modeled. The sequence of tasks is as follows:
- Create a 3D grid containing the aquifer domain.
- Create a 2D grid with the same horizontal extension as the 3D grid.
- Simulate independently the elevation of the bottom of the stratigraphic units in the 2D grid by using DeeSse. The spatial structure of those surfaces is inspired by the (detrended) topography of a braided river (Figure 3). This continuous variable is used as training image. The corresponding elevations available in the 2D point set are used to condition the simulation.
- Use the following erosion rule to ensure coherency: more recent layers (above) erode older (underlying) ones.
- Create a region for each layer L1 to L4. Those regions will be used in the following steps.
The result is given in Figure 4.
3.2. Step 2 – Modeling the lithology within the stratigraphic units
Stratigraphic units L1 and L4 have been shown to be homogeneous or considered as such. Unit L1 is then filled with facies F0 and unit L4 with facies F4.
For units L2 and L3, the analysis of the borehole data has shown that the dip and azimuth are varying in space. The simulation of the facies within those layers must take this feature into account.
Three tasks will thus be necessary to achieve this step:
- Interpolate the orientation data in order to have orientation maps for both the dip and azimuth;
- Create a (set of) training image(s) representing the F2 and F3 objects, F1 being the matrix;
- Within each stratigraphic unit, use DeeSsewith the corresponding orientation data and TIs to simulate the lithology (facies);
Those tasks are detailed hereunder.
3.2.1. Interpolating the orientations within the stratigraphic units
The visualization of the data shows that the azimuth and dip are varying in a rather smooth manner but have sharp interfaces at the boundary between the different stratigraphic units. It means that it is reasonable to interpolate them using kriging. However, interpolation in each layer must be achieved independently from the other layers using only the data contained in that layer. The variogram models may also be adapted to each layer, if there is an evidence that they are different from one layer to the other.
Result for the dip is shown in Figure 5.
3.2.2. Creating the training image
The aim of this step is to build an elementary 3D training image for the simulation in stratigraphic units L2 and L3. We propose to consider thin ellipsoid objects “parallel” to the horizontal plane, resulting in a stationary TI. The local orientations of the structures will be accounted for during the simulation step by specifying local azimuth and dip angles in the input parameters of the simulation.
First, a 3D grid is created for the TI. Then, a Boolean object simulator is used with objects corresponding to facies F2 and F3, and a matrix being made of facies F1. The objects are made of the superposition of two half ellipsoids, one being rotated and translated with respect to the other.
A TI generated using this approach is shown in Figure 6.
3.2.3. Simulating the lithology
DeeSseis used at this step again to simulate the lithology within stratigraphic units L2 and L3 independently (Figure 7). Within each layer, the appropriate TI and orientation maps are used as parameters, along with the relevant conditioning data from the boreholes.
3.3. Step 3 – Modeling the porosity and permeability within the lithology
In this step, the various facies regions are filled with continuous properties: porosity and permeability. Sequential Gaussian simulation (SGS) has been chosen to simulate the porosity. The simulation is done with a varying anisotropy direction, using the orientations interpolated at step 3.2.1
The distribution of the porosity can be approximated by a Weibull distribution. Suggested values for the shape and scale parameters of this distribution for the various facies regions are given in Table 1.
Once porosity has been simulated, permeability can be directly computed from it using the relation of Kozeny-Carman:
where phi is the porosity and K, the permeability. Parameter A is a value that can vary depending on the rock type. Table 1 gives values for this parameter for each facies region.
Table 1: Parameter values for the distribution of porosity and for the computation of permeability using the Kozeny-Carman relation for each facies.
One realization of porosity is given in Figure 8.
The workflow presented in this paper is very generic and could be adapted to different modeling software. At each step, various methods/tools can be used depending on the amount of local data and conceptual knowledge available.
Our main point was to show how DeeSsecan be useful at the various stages of the modeling procedure, for simulating categorical but also continuous variables. DeeSseis particularly suited to accommodate non-stationary parameters such as orientations and trend maps. Consequently, DeeSseis able to generate realizations reflecting the potentially high complexity of geological domains.
Finally, note that the whole procedure can be semi-automated from the beginning (well data) to the final 3D model of the stratigraphic units, the lithology and permeability. Once the modeling decisions and parameters are fixed, all the steps can run automatically to generate multiple realizations efficiently.
The authors would like to thank Alexandre Boucher for his implication during the preparation of the synthetic case study used in this paper.
Comunian, A., Renard, P., Straubhaar, J., & Bayer, P. (2011). Three-dimensional high resolution fluvio-glacial aquifer analog–Part 2: Geostatistical modeling. Journal of hydrology, 405(1-2), 10-23.
Straubhaar J., Renard P., Mariethoz G.  Conditioning multiple-point statistics simulations to block data. Spatial Statistics 16 (2016)
Mariethoz G., Renard P., Straubhaar J.  The Direct Sampling method to perform multiple‐point geostatistical simulations. Water Resources Research, Vol. 46.