Numerical modeling (geology)

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. for example, folding and mantle convection. With advances in computer technology, the accuracy of numerical models has been improved. == Components ==

describe the flow of viscous fluid (the mantle rock). Second, since these equations are difficult to solve, discretization and numerical methods are chosen to make an approximation to the governing equations. Then, algorithms in the computer can calculate the approximated solutions. Finally, interpretation can be made from those solutions. For instance, in mantle convection modeling, the flow of mantle can first be visualized. Then, the relationship between the patterns of flow and the input parameters may be concluded. A general numerical model study usually consists of the following components: for example, the wave equation. • Discretization methods and numerical methods convert those governing equations in the mathematical models to discrete equations. These discrete equations can approximate the solution of the governing equations. Common methods include the finite element, finite difference, or finite volume method that subdivide the object of interest into smaller pieces (element) by mesh. These discrete equations can then be solved in each element numerically. The discrete element method uses another approach, this method reassembling the object of interest from numerous tiny particles. Simple governing equations are then applied to the interactions between particles. • Algorithms are computer programs that compute the solution using the idea of the above numerical methods. • Interpretations are made from the solutions given by the numerical models. == Properties ==

A good numerical model usually has some of the following properties: • Consistent: Numerical models often divide the object into smaller elements. If the model is consistent, the result of the numerical model is nearly the same as what the mathematical model predicts when the element size is nearly zero. In other words, the error between the discrete equations used in the numerical model and the governing equations in the mathematical model tends to zero when the space of the mesh (size of element) becomes close to zero. • Stable: In a stable numerical model, the error during the computation of the numerical methods does not amplify. The error of an unstable model will stack up quickly and lead to an incorrect result. A stable and consistent numerical model has the same output as the exact solution in the mathematical model when the spacing of the mesh (size of element) is extremely small. • Converging: The output of the numerical model is closer to the actual solution of the governing equations in the mathematical models when the spacing of mesh (size of element) reduces, which is usually checked by carrying out numerical experiments. • Conserved: The physical quantities in the models, such as mass and momentum, are conserved. Since the equations in the mathematical models are usually derived from various conservation laws, the model result should not violate these premises. • Bounded: The solution given by the numerical model has reasonable physical bounds with respect to the mathematical models, for instance mass and volume should be positive. • Accurate: The solution given by the numerical models is close to the real solution predicted by the mathematical model. == 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. Eulerian The Eulerian approach considers the changes of the physical quantities, such as mass and velocity, of a fixed location with time. 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. This equation is commonly used in numerical modeling in geology. The discrete element method uses another approach. The object is considered an assemblage of small particles. 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. As a result, the spectral method converges exponentially and is suitable for solving problems involving a high variability in time or space. 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. Consider a function f(x) with single-valued derivatives that are continuous and finite functions of x, according to Taylor's theorem: f(x+\Delta x) = f(x) + \Delta x f'(x) + \frac{1}{2}\Delta x ^2 f(x) + \frac{1}{6}\Delta x^3 f'(x)+\cdots and f(x-\Delta x) = f(x) - \Delta x f'(x) + \frac{1}{2}\Delta x ^2 f(x) - \frac{1}{6}\Delta x^3 f'(x)+\cdots Summing up the above expressions: 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, 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: and the space in the rock mass maybe filled with other substances like air and water. • Not elastic: Rock cannot perfectly revert to its original shape after stress is removed. Models that model rock as a discontinuum, using methods like discrete element and discrete fracture network methods, are also commonly employed. 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. Pecube Pecube is one of the numerical models developed to predict the thermochronological pattern. The first three terms on the right-hand side are the heat transferred by conduction in x , y and z directions while A is the advection. \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} + w \frac{\partial T}{\partial z} = \frac{\partial}{\partial x} \kappa \frac{\partial T}{\partial x} + \frac{\partial}{\partial y} \kappa \frac{\partial T}{\partial y} + \frac{\partial}{\partial z} \kappa \frac{\partial T}{\partial z} + A After the temperature field is constructed in the model, particle paths are traced and the cooling history of the particles can be obtained. The model is three-dimensional; and finite difference method. These two methods have been shown to produce similar results if the mesh is fine enough. grid used in MODFLOW|thumb 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 more features, including: • Interactions between groundwater and surface-water systems FLAC The Fast Lagrangian Analysis of Continua (FLAC) is one of the most popular approaches in modeling crustal dynamics. The approach has been used in 2D, 2.5D, and 3D studies of crustal dynamics, in which the 2.5D results were generated by combining multiple slices of two-dimensional results. This model uses the code called Flamar, which is a FLAC-like code that combines finite difference and finite element methods. The boundary condition used at the bottom is called "Winkler's pliable basement". It is at hydrostatic equilibrium and it allows the base to slip freely horizontally. In the model, the temperature of the core-mantle boundary (inner boundary) is a constant of 4273 K (about 4000°C), while that at the boundary between crust and mantle (outer boundary) is 973 K (about 700°C). finite volume, finite difference and spectral methods have all been used in modeling mantle convection, and almost every model used an Eulerian grid. Current approaches mostly uses a fixed and uniform grid. Finite difference approach In the 1960s to 1970s, mantle convection models using the finite difference approach usually used second-order finite differences. Finite volume approach Mantle convection modeled by finite volume approach is often based on the balance between pressure and momentum. The equations derived are the same as the finite difference approach using a grid with staggered velocity and pressure, in which the values of velocity and the pressure of each element are located at different points. In modeling three-dimensional geometry of the Earth, since the parameters of mantles vary at different scales, multigrid, which means using different grid sizes for different variables, is applied to overcome the difficulties. 'Yin-Yang' grid, and spiral grid. Finite element approach In the finite element approach, stream functions are also often used to reduce the complexity of the equations. modeling two-dimensional incompressible flow in the mantle, was one of the popular codes for modeling mantle convection in the 1990s. and 3D. Spectral method The spectral method in mantle convection breaks down the three-dimensional governing equation into several one-dimensional equations, which solves the equations much faster. It was one of the popular approaches in early models of mantle convection. However, the lateral changes of viscosity of mantle are difficult to manage in this approach, and other methods became more popular in the 2010s. The mantle convection model is fundamental in modeling the plates floating on it, and there are two major approaches to incorporate the plates into this model: rigid-block approach and rheological approach. Modeling of the flow of Earth's liquid outer core is difficult because: for example the Glatzmaier-Roberts model. Finite difference method has also been used in the model by Kageyama and Sato. Some study also tried other methods, like finite volume and finite element methods. . Seismology File:Global Seismic Wave Propagation Simulation.gif|thumb|299x299px|Simulation of seismic wave propagation through the Earth. However, due to limitations in computation power, in some models, the spacing of the mesh is too large (compared with the wavelength of the seismic waves) so that the results are inaccurate due to grid dispersion, in which the seismic waves with different frequencies separate. Some researchers suggest using the spectral method to model seismic wave propagation. == Errors and limitations ==

Sources of error While numerical modeling provides accurate quantitative estimation to geological problems, there is always a difference between the actual observation and the modeling results due to: Limitations Apart from the errors, there are some limitations in using numerical models: • Users of the models need a high level of knowledge and experience to prevent misuse and misinterpretation of results. == See also ==