Introduction to numerical simulations

Introduction to numerical simulations

Over the past decades, a series of sophisticated multidimensional (2D, 3D, dD) simulation technologies have been developed to address the need for a precise definition of the earth's (sub)surface conditions, geometries and their development (Turner, 1991, 2000). The models have to have the capability to effectively explain the geometry of rocks and minerals, their spatial and temporary relationships, their displacements by tectonic forces and/or the fluid flow through a rock (Kelk, 1991). A model according to the definition above is not necessarily a computer model. The terms model and simulation are regularly used as synonyms for each other although they describe different concepts.

A model is a logical representation of one or more processes that simultaneously or sequentially affect a material (or a material property) and therewith induce some kind of change of the material or its properties.

A simulation is the actual application of a model to a specific case, often over time. This can, but not necessarily has to be, a computer simulation.

Models are typically used when it is impossible or impracticable to create experimental conditions that enable the direct observation of processes. As such, a model can be seen as a simplified view of the complex reality. In our everyday work, geologists (or in more general terms scientists) constantly invent, improve or dismiss models.

The short definition of the terms model and simulation already hint at one of the main advantages of using models, and subsequently simulations derived from these models, in geology: most of the geological processes active in nature are simply to slow for direct observations (e.g. mountain building or subduction), are inaccessible (e.g. emplacement of plutons/batholiths deep in the earth's crust, subduction) and/or are unobservable “in situ” (e.g. metamorphic reactions between minerals, subduction).

Computer simulations developed hand-in-hand with the rapid growth of the computer. More and more complex simulations were performed using the increasing computational power. Usually, computer simulations are performed for problems where no closed form analytical solution can be obtained.

Computer simulations can be classified following different criteria such as:

• stochastic or deterministic models

• steady-state or dynamic models

• continuous or discrete models

• etc.

Computer simulations in geology are always steady-state or dynamic models. In steady-state models, a set of parameters and their relationship is defined. The simulation then tries to find a state where all these parameters are in equilibrium. Countless applications of steady-state models to geological problems have been published and the range of problems analyzed using this kind of simulation ranges from exhumation of the European Alps (Bernet et al., 2001), climatic controls on hillslope angles in the Himalayas (Gabet et al., 2004) to analyzing rheological and porosity equations during compaction (Boudreau & Bennett, 1999) or predicting the freshwater depth in islands (Bobba, 1998).

In dynamic models, the set of parameters that is initially defined changes during the simulation. In geological simulations, this usually is the change over time but may also be an increased or decreased flux of e.g. a liquid or a change of concentration. Examples for dynamic models are e.g. the simulation of time scales of magmatic processes (Edwards & Russell, 1998), simulations of subduction zones (e.g. Billen et al., 2003, Yoshioka & Wortel, 1995, Shijie & Gurnis, 1995), growth of fibres (Köhn et al., 2000), development of crystal morphology during crack opening (Hilgers et al., 2001) to bathymetric limits of nautiloids (Hewitt et al., 1989).

It is therefore more practical to classify geological simulations according to the numerical methods that are used:

• Stochastic models (e.g. Monte Carlo models, Ising models, Potts models etc.)

Stochastic models are a class of computational algorithms that rely on chance. They are most often used to simulate physical systems with a large number of coupled degrees of freedom such as fluids or cellular structures. Examples include the evaluation of ground conditions (Rosenbaum & Stevens, 1991), reservoir modelling (Tamhane et al., 1999) or grain growth (Enting, 1977; Yu et al., 2008).

• Finite difference models (FDM)

The finite difference method is well suited to solve ordinary and partial differential equations. The area (or volume) to be simulated is first divided into a finite amount of nodes. At each node, a finite difference between neighboring nodes (in space and/or time) can be calculated that is used to replace the derivatives appearing in the differential equation. Finite difference models are very often used to solve problems or fluid flow or thermodynamics.

One example for a finite difference simulation would be the reconstruction of the thermal history of variscian crust (Kosakowski et al., 1998). FDM are also frequently used to solve the Navier-Stokes equation (e.g. Cottet, 1990).

• Finite element models (FEM)

Finite element models are also well suited to solve partial differential equations, but the area (or volume) to be simulated is divided into elements (usually simple geometric objects). Each element is assigned a set of basis functions that describe its properties and relation to the surrounding elements. FEMs are most often used in simulations that include the movement or displacement of material, have a varying precision of the simulation domain or where the simulation domain is very complex.

The are several differences between the finite element and finite difference models, the two main differences probably are:

• FDMs are an approximation to the differential equation, FEMs are an approximation to the solution.

• FEMs are able to handle very complex geometries and/or boundary conditions with varying precision while the FDM (in its basic form) is restricted to a regularly spaced grid of nodes.

A typical example for a finite element simulation would be the behavior of rocks at various deformation conditions (e.g. Paterson, 2007, Li et al., 2007).

• Finite Volume models (FVM)

This method is a method to represent and evaluate partial differential equations as algebraic equations. It is similar to the FDM but each point in the grid is thought to have a small volume. Each flux entering a given volume is thought to be identical to a flux leaving neighboring volumes. This method therefore is conservative and very often used in computational fluid dynamics (e.g. Yoshida & Ogawa, 2005).

• Boundary models

Boundary models are well suitable to simulate the movement of distinct boundaries due to driving forces. The driving forces may for example be the surface energy of grains or a chemical potential. Usually, only the change of parameters in the immediate surroundings of the boundary are important for this kind of simulation (Bons et al., 2007 and references therein).

• Phase field models

Phase field models are mathematical models to solve interfacial problems. The major difference to the boundary models is that the boundaries are not represented by sharp change over an infinitesimal thin boundary but rather over a larger area (the phase field). Pure phases have distinct values (e.g. +1 or -1) while within boundaries intermediate values do occur. They are most often used in material science but have recently been successfully applied to geological problems. (e.g. Kazaryan et al., 2002, Wendler et al.,accepted)

To comprehensively explain the various different methods would go far beyond the scope of this thesis. For a detailed description of e.g. Monte Carlo models, finite difference models or the finite element method, the book Microdynamics Simulation (Bons et al., 2007) is a very good starting point.