Technique to solve geological problems by computational simulation
In geology, numerical modeling is a widely applied technique to tackle complex geological problems by computational simulation of geological scenarios.
Numerical modeling uses mathematical models to describe the physical conditions of geological scenarios using numbers and equations. Nevertheless, some of their equations are difficult to solve directly, such as partial differential equations. With numerical models, geologists can use methods, such as finite difference methods, to approximate the solutions of these equations. Numerical experiments can then be performed in these models, yielding the results that can be interpreted in the context of geological process. Both qualitative and quantitative understanding of a variety of geological processes can be developed via these experiments.
Numerical modelling has been used to assist in the study of rock mechanics, thermal history of rocks, movements of tectonic plates and the Earth's mantle. Flow of fluids is simulated using numerical methods, and this shows how groundwater moves, or how motions of the molten outer core yields the geomagnetic field.
History.
Prior to the development of numerical modeling, analog modeling, which simulates nature with reduced scales in mass, length, and time, was one of the major ways to tackle geological problems, for instance, to model the formation of thrust belts. Simple analytic or semi-analytic mathematical models were also used to deal with relatively simple geological problems quantitatively.
In the late 1960s to 1970s, following the development of finite-element methods in solving continuum mechanics problems for civil engineering, numerical methods were adapted for modeling complex geological phenomena, for example, folding and mantle convection. With advances in computer technology, the accuracy of numerical models has been improved. Numerical modeling has become an important tool for tackling geological problems, especially for the parts of the Earth that are difficult to observe directly, such as the mantle and core. Yet analog modeling is still useful in modeling geological scenarios that are difficult to capture in numerical models, and the combination of analog and numerical modeling can be useful to improve understanding of the Earth's processes.
Components.
A general numerical model study usually consists of the following components:
Properties.
A good numerical model usually has some of the following properties:
Computation.
The following are some key aspects of ideas in developing numerical models in geology. First, the way to describe the object and motion should be decided (kinematic description). Then, governing equations that describe the geological problems are written, for example, the heat equations describe the flow of heat in a system. Since some of these equations cannot be solved directly, numerical methods are used to approximate the solution of the governing equations.
Kinematic descriptions.
In numerical models and mathematical models, there are two different approaches to describe the motion of matter: Eulerian and Lagrangian. In geology, both approaches are commonly used to model fluid flow like mantle convection, where an Eulerian grid is used for computation and Lagrangian markers are used to visualize the motion. Recently, there have been models that try to describe different parts using different approaches to combine the advantages of these two approaches. This combined approach is called the arbitrary Lagrangian-Eulerian approach.
Eulerian.
The Eulerian approach considers the changes of the physical quantities, such as mass and velocity, of a "fixed location" with time. It is similar to looking at how river water flows past a bridge. Mathematically, the physical quantities can be expressed as a function of location and time. This approach is useful for fluid and homogeneous (uniform) materials that have no natural boundary.
Lagrangian.
The Lagrangian approach, on the other hand, considers the change of physical quantities, such as the volume, of "fixed elements" of matter over time. It is similar to looking at a certain collection of water molecules as they flow downstream in a river. Using the Lagrangian approach, it is easier to follow solid objects which have natural boundary to separate them from the surrounding.
Governing equations.
Following are some basic equations that are commonly used to describe physical phenomena, for example, how the matter in a geologic system moves or flows and how heat energy is distributed in a system. These equations are usually the core of the mathematical model.
Continuity equation.
The continuity equation is a mathematical version of stating that the geologic object or medium is continuous, which means no empty space can be found in the object. This equation is commonly used in numerical modeling in geology.
One example is the continuity equation of mass of fluid. Based on the law of "conservation of mass", for a fluid with density formula_1 at position formula_2 in a fixed volume formula_3 of fluid, the rate of change of mass is equal to the outward fluid flow across the boundary formula_4:
formula_5
where formula_3 is the volume element and formula_7 is the velocity at formula_2.
In Lagrangian form:
formula_9
In Eulerian form:
formula_10
This equation is useful when the model involves continuous fluid flow, like the mantle is over geological time scales.
Momentum equation.
The momentum equation describes how matter moves in response to force applied. It is an expression of Newton's second law of motion.
Consider a fixed volume formula_3 of matter. By the law of "conservation of momentum", the rate of change of volume is equal to: 
formula_14
where formula_3 is the volume element, formula_16 is the velocity.
After simplifications and integrations, for any volume formula_3, the Eulerian form of this equation is:
formula_18
Heat equation.
The heat equations describe how heat energy flows in a system.
From the law of conservation of energy, the rate of change of energy formula_19 of a fixed volume formula_3 of mass is equal to:
Mathematically:
formula_26
where formula_3 is the volume element, formula_16 is the velocity, formula_29 is the temperature, formula_30 is the conduction coefficient and formula_31 is the rate of heat production.
Numerical methods.
Numerical methods are techniques to approximate the governing equations in the mathematical models.
Common numerical methods include finite element method, spectral method, finite difference method, and finite volume method. These methods are used to approximate the solution of governing differential equations in the mathematical model by dissecting the domain into meshes or grids and applying simpler equations to individual elements or nodes in the mesh. 
The discrete element method uses another approach. The object is considered an assemblage of small particles.
Finite element method.
The finite element method subdivides the object (or domain) into smaller, non-overlapping elements (or subdomains) and these elements are connected at the nodes. The solution for the partial differential equations are then approximated by simpler element equations, usually polynomials. Then these element equations are combined into equations for the entire object, i.e. the contribution of each element is summed up to model the response of the whole object. This method is commonly used to solve mechanical problems. The following are the general steps of using the finite element method:
Spectral method.
The spectral method is similar to the finite element method. The major difference is that spectral method uses basis functions, possibly by using a Fast Fourier Transformation (FFT) that approximates the function by the sum of numerous simple functions. These kinds of basis functions can then be applied to the whole domain and approximate the governing partial differential equations. Therefore, each calculation takes the information from the whole domain into account while the finite element method only takes the information from the neighborhood. As a result, the spectral method converges exponentially and is suitable for solving problems involving a high variability in time or space.
Finite volume method.
The finite volume method is also similar to the finite element method. It also subdivides the object of interest into smaller volumes (or elements), then the physical quantities are solved over the control volume as fluxes of these quantities across the different faces. The equations used are usually based on the conservation or balance of physical quantities, like mass and energy.
The finite volume method can be applied on irregular meshes like the finite element method. The element equations are still physically meaningful. However, it is difficult to get better accuracy, as the higher order version of element equations are not well-defined.
Finite difference method.
The finite difference method approximates differential equations by approximating the derivative with a difference equation, which is the major method to solve partial differential equations.
 Consider a function formula_32 with single-valued derivatives that are continuous and finite functions of formula_33, according to Taylor's theorem:
formula_34
and
formula_35
Summing up the above expressions:
formula_36
Ignore the terms with higher than 4th power of formula_37, then:
formula_38
formula_39
The above is the "central-difference" approximation of the derivatives, which can also be approximated by "forward-difference":
formula_40
or "backward-difference":
formula_41
The accuracy of the finite differences can be improved when more higher order terms are used.
Discrete element method.
The discrete element method, sometimes called distinct element method, is usually used to model discontinuous materials, such as rocks with fractures like joints and bedding, since it can explicitly model the properties of discontinuities. This method was developed to simulate rock mechanics problems at the beginning.
The main idea of this method is to model the objects as an assemblage of smaller particles, which is similar to building a castle out of sand. These particles are of simple geometry, such as a sphere. The physical quantities of each particle, such as velocity, are continuously updated at the contacts between them. This model is relatively computationally intensive, as a large quantity of particles needs to be used, especially for large-scale models, like a slope. Therefore, this model is usually applied to small-scale objects.
Bonded-particle model.
There are objects that are not composed of granular materials, such as crystalline rocks composed of mineral grains that stick to each other or interlock with each other. Some bonding between particles is added to model this cohesion or cementation between particles. This kind of model is also called a bonded-particle model.
Applications.
Numerical modeling can be used to model problems in different fields of geology at various scales, such as engineering geology, geophysics, geomechanics, geodynamics, rock mechanics, hydrogeology, and stratigraphy. The following are some examples of applications of numerical modeling in geology.
Specimen to outcrop scale.
Rock mechanics.
Numerical modeling has been widely applied in different fields of rock mechanics. Rock is a material that is difficult to model because rock are usually:
In order to model the behaviors of rock, a complex model that takes all the above characteristics into account is needed. There are many models modeling rock as a continuum using methods like finite difference, finite element, and boundary element methods. One of the disadvantages is that the ability of modeling cracks and other discontinuities is usually limited in these models. Models that model rock as a discontinuum, using methods like discrete element and discrete fracture network methods, are also commonly employed. Combinations of both methods have also been developed.
Numerical modeling enhances the understanding of mechanical processes in rock by conducting numerical experiments, and is useful for design and construction works.
Regional-scale.
Thermochronology.
Numerical modeling has been used to predict and describe the thermal history of the Earth's crust, which allows geologists to improve their interpretation of thermochronological data. Thermochronology can indicate the time at which a rock cooled below a particular temperature. Geologic events, like the development of a faults and surface erosion, can change the thermochronological pattern of samples collected on the surface, and it is possible to constrain the geologic events by these data. Numerical modeling can be used to predict the pattern.
The difficulties of thermal modeling of the Earth's crust mainly involve the irregularity and the changes of the Earth's surface (mainly erosion) through time. Therefore, in order to model the morphological changes of the Earth's surface, the models need to solve heat equations with boundary conditions that change with time and have irregular meshes.
Pecube.
Pecube is one of the numerical models developed to predict the thermochronological pattern. It solves the following generalized heat transfer equation with advection using finite element method. The first three terms on the right-hand side are the heat transferred by conduction in formula_42, formula_43 and formula_44 directions while formula_45 is the advection.
formula_46
After the temperature field is constructed in the model, particle paths are traced and the cooling history of the particles can be obtained. The pattern of thermochronological age can then be computed. 
Hydrogeology.
In hydrogeology, groundwater flow is often modeled numerically by the finite element method and finite difference method. These two methods have been shown to produce similar results if the mesh is fine enough.
MODFLOW.
One of the well-known programs in modeling groundwater flow is MODFLOW, developed by the United States Geological Survey. It is a free and open-source program that uses the finite difference method as the framework to model groundwater conditions. The recent development of related programs offers