How Does it Work?

Geoscientific Principles & Underpinnings

Geomodeller draws on many aspects of Mathematics, Physics, Material Science, Geology, Geostatistics and Geophysics.

Perhaps the most fundamental improvement follows from observing how geological fabrics behave. This is a natural extension to the displine of geostatistics. Initially, numerical variabilities such as the grade of a gold vein, has driven geostatistical thinking over more than 30 years. It was inevitable then, that the form and shapes of naturally occurring geological bodies would also come under scrutiny. This work was conducted in France over several years (Chilès et.al. 2004).

The major outcome includes an insight into how the relationships between structural geology observations and contacts can be approximated by fitting a cubic function to the experimental variogram for the contacts and a linked differentiated cubic function to the parallel and perpendicular components of the dip vectors. Potential Field theory is suitable as it has close properties to the desired model.

1.1 Inference of the potential field

In usual geostatistical applications, the covariance or variogram of the variable under study is modelled from the sample variogram of the data. When modeling geology, we have no measurement of the potential T(x), and the potential increments used for the interpolation cannot be used for the inference of K since they all have a zero value. In its first implementation, the algorithm was used heuristically with a covariance model arbitrarily chosen by the user. That choice was more or less rationalised according to the following considerations:

• At the scale considered, geological interfaces are smooth rather than fractal surfaces. This implies that the covariance is twice differentiable. A cubic model is a good compromise among the various possible models, because it has the necessary regularity at the origin and has a scale parameter that can accommodate various situations.
• The scale parameter a and sill C of the covariance K(h) determine the sill of the variogram of the partial derivatives: it is equal to 14 C / a² in the case of an isotropic cubic covariance. When there is no drift and the geological body is isotropic (eg a granitic intrusion), the unit gradient vector can have any direction so that its variance is equal to one. The variance of each partial derivative is then equal to one third.
• A consistent choice for C once the scale parameter a is chosen, is therefore C = a² / 42. That value is an upper bound for C when the potential field has a drift. This is because the mean of the potential gradient is not equal to zero and its variance is less than one (its quadratic mean is zero by definition).
• Sensible measurement variances can also be defined (nugget effects).

The use of a heuristic model, however, implies two limitations:

1. The choice is usually not the best one; and
2. More importantly, this precludes any evaluation of the magnitude of the interpolation error.

A means to infer the covariance is thus a core issue of that approach. Since K cannot be inferred from the potential increments, it is inferred is from the gradient data. This is possible because the covariances of the partial derivatives can be derived from assuming a potential field.

In the case of an isotropic covariance K(h), which for simplicity will be denoted K(r) as a function of r = ||h||, the covariance of, say, δT(x) / δu and δT(x+h)/δu is:

• K"(||h||) when h is parallel to the u axis,
  
• K (||h||) / ||h|| when h is orthogonal to the u axis.

1.1.1 Isotropy vs anisotropy

The assumption of an isotropic covariance model is the standard starting position. It can become too restrictive and with V1.3 Geomodeller, options for anisotropy are more available. This allows you to model more reliably, thinner bodies such as dykes.

In practice the covariance K(h) is the sum of several cubic components Kp(h), each one possibly displaying a zonal or geometric anisotropy. With the initial formulation of this capability, the main anisotropy axes u, v, w, are common to all the components. Current development work allows for each geological series to have its own definition of anisotropy.

Using the formulae for the selected model the covariance parameters of K (nugget effect, scale parameter of each covariance component in the three main directions, sill of each component) are chosen so as to lead to a satisfactory global fit of the directional sample variograms of the three components of the gradient.

Figure 1 shows an example of such a fitting. 1485 structural data was sampled in an area of about 70 × 70 km² in the Limousin (Massif Central, France). The main (u, v, w) coordinates here coincide with the geographical (x, y, z) coordinates. Since the structural data is all located on the topographic surface, the variograms have been computed in the horizontal plane only. Note that the sill of the variogram of the vertical component is much lower than that of the horizontal components. This is due to the fact that the layers are subhorizontal so that the vertical component of the gradient displays limited variations around its non-zero mean. The model K includes three components, the second of which only depends on the horizontal component of h and the third one on the N-S component (zonal anisotropies).

1.1.2 Stationarity Property

An important assumption made is that any trend in the spatial variability of the geology can be removed or detrended. The aim in detrending the "geology", is to achieve a state where any remaining variability is essentially random and the "geology" is stationary.

This allows cokriging to be performed in the framework of a random function model. Formally, the mathematical function used to model geology T(x), is assumed to be a random function with a polynomial drift, and a stationary covariance K(h).

1.1.3 Detrending geology

Since the vertical usually plays a special role, the degree of the polynomial drift can be higher vertically than horizontally and the covariance can be anisotropic.

Not all geology is layer-cake. An intrusive geological body that has the shape of an ellipsoid, can be detrended assuming a quadratic drift, ie use ten coefficients for the drift function with degree less than or equal to two.

So, in general, a three dimensional quadratic drift function is the current practise. In theory, other detrending methods like sinusoidal terms could be used, but in usual applications geology is not regular enough for that.

It is also important to note that as no attempt is being made to follow the genesis of the geology, there is no need to try and "mass balance on a section by section basis". We are modeling what is observed, not simulating a folding and faulting progression.

1.1.4 Interpolating geology using the potential field method

In summary, the potential-field method defines a geological interface as an implicit surface, namely a particular isosurface of a scalar field defined in the 3D space - the potential field. The 3D interpolation of that potential field, based on universal cokriging, provides isosurfaces that honour all the data. Since no data measures the potential field itself, its covariance cannot be inferred directly. However, the covariance can be determined from the structural data, which makes it possible to associate sensible cokriging standard deviations to potential-field estimates and to translate them into uncertainties of the 3Dmodel.

The implementation of joint structural geology and contacts interpolation based upon experimentally derived methods, is a unique breakthrough in turning geology mapping and interpretation into a quantitative science. It is a fallacy to now claim there is no basis for how one should interpolate geology, so any / all methods are equally adaptable.

1.2 Radial Basis Functions

Universal kriging implies that all observations have an influence on every part of the model, no matter how far away an observation may be. However, it is now recognized that the majority of the influence is from those observations that are local to that part of your model. This turns out to be a local radial basis interpolator, using classical geostatistics.

A very similar outcome can also be arrived at with alternate mathematical thinking. One candidate is the use of the biharmonic equation and thin plate theory from engineering science, to interpolate the geology. This can be easily formuated as a radial basis function. In this case, there is no attempt to honour observed characteristics of how geology bodies should be modeled, just the condition that all observed contacts are honoured. Where there are many observations of geological contacts, say in an in-mine context with lots of bore-holes, this method will produce a satisfactory prediction of each geological contact surface.

1.2.1 Surface vs Volume

The important point to make here is that one should interpolate volumes rather than surfaces, and this is central to the Geomodeller potential field method. Thin plate surface interpolations do not naturally have this property. As the number of observations become sparser, the breakdown in being able to produce realistic geological bodies will become very pronounced for the thin plate spline methods. There is no natural constraint to have Top and Bottom surfaces for a unit follow the same trends, if using surface splining.

Not coincidently, well behaved geological volumes are also important for the geophysical modeling and inversion aspects of Geomodeller. Both the total mass and total magnetization of a unit is inherently tied to the volume of the body. It is now more generally recognized that it is the volume of the unit, not just individual bounding surfaces, that is important in producing believable 3D geology models. In later sections of this manual, many references to characterizing allowable volume changes while inversion is being pursued, will be made.

1.2.2 Faults

Several methods are used to handle faults. If faults delimit blocks and the geology is not correlated from one block to the other, it obviously suffices to process each block separately. These are termed "TERRAIN" faults.

For most other cases, the method used in Geomodeller is a transposition to 3D potential fields of the method proposed by Maréchal (1984). This handles faults in the 2D interpolation of the elevation of interfaces, where faults are entered as external drift functions. This method requires knowledge of the fault plane and also its zones of influence.

For the case where a normal fault intersects the whole study zone, the geology is divided into two sub zones D and D'. This fault induces a discontinuity of the potential field, whose amplitude is not known. Cokriging can accommodate that discontinuity whatever its amplitude by introducing a drift function complementing the L polynomial drift or detrending functions above, for example:

f L+1(x) = 1D(x).

If the polynomial drift functions include the first coefficient f 1(x) = x (first coordinate) due to the presence of a linear trend of the potential field, and we have good reasons to suspect not only a discontinuity but also a change of slope of the drift when crossing the fault, it is advisable to also introduce an additional drift function such as:

f L+2(x) = x 1D(x).

A finite fault can be modeled with a drift function with a bounded support, and whose value vanishes on the support boundaries; inside that support, the function takes on positive values on one side of the fault plane, with a maximum at the centre of the fault, and negative values on the other side. The fault plane is unlikely to be a planar surface. It is often only known by some points on its surface and unit vectors orthogonal to it. Its geometry can thus be modeled by a potential field too. Currently, faults are not honored in inversion.

1.3 Statistical Laws & Probability Distribution Functions

At the heart of much of the thinking in Geomodeller Inversion, are stochastic processes that respect statistical laws.

1. The main potential field interpolator in Geomodeller is based upon universal dual kriging, and as such represents a high level of geostatistics.
2. The laws that govern allowable changes in the volumes of geological bodies are constrained to follow Exponential and /or Log Normal distributions.
3. Finally, the lithological properties that you are required to specify for all your rock properties, also respect Normal, Lognormal etc. laws.

Each of these laws has a well defined and well understood mathematical basis.

For constraining the geometry of the geological bodies (their volumes), a Weibull probability distribution function is introduced to give access to more appropriate rules. The Weibull distribution has gained a large following in the automotive industry as a means of predicting mean time between failures for car component parts. The property sought is that the initial starting model should have a high probability of being right. This follows from the proposition that the geologist knows what he is doing and the reference model is close to the truth.

1.4 The need for a Close Starting Model

One of the central issues needed for a successful use of the fully constrained lithological inversion method, is to create a starting or reference model that is not only consistent with all the geological observations, but also "close" to honoring the observed geophysics signals for the majority of your study area.

The simple reason for this requirement follows from the use of Monte Carlo techniques. There is a very detail process of changing individual properties or geometries. A close staring model ensures you maximize the likelihood and minimize the time needed to explore probability space to firstly converge to an acceptable mis-fit, and then continue through gathering good statistical results.

Taking the example case of having an observed gravity grid and a forward model of the gravitational response from your model, you should immediately examine the two grids for a spatial comparison, and also an examination of the statistics, as shown in the above figure. In the situation shown, the likelihood of your starting position migrating to the required observed position, while also honouring the known constraints, is very low. It is for this reason, you are advised to continue working on improving your starting or reference model. This is done by looking at problem areas via sections, making changes and updating the 3D forward modelling.

To retain realistic geological bodies, you are advised to make use of volume and shape ratio constraints. Of course, if this constraint is not used, an apparent conversion to a good fit can be achieved automatically, but you are left with a fragmented predicted geology that bears little relation to your staring model.

Partners:

Affiliated Sites:

Intrepid Geophysics

GeoIntrepid

BRGM (France)

Communities:

3D Geology Forum

Recent News:

Website Translation:

If there is no flag for your preferred language:

1 Copy URL of webpage
2 Go to Babelfish Site
3 Paste URL of webpage