1
Preliminaries

This chapter provides a short introduction to mathematical models consisting of systems of partial differential equations (PDEs) along with auxiliary (boundary and initial) conditions. We discuss how these equations can be solved, either exactly or using numerical methods. We also briefly consider the important issues of precision and stability of a numerical solution. A Matlab script is provided at the end of the chapter to enable readers to compare an analytical solution with its corresponding numerical approximation.

1.1 Mathematical Models

The application of the principles of conservation of mass, momentum, and energy combined with experimentally derived laws produces sets of PDEs that describe variations in velocity (or displacement), pressure, and temperature in space and time. When combined with boundary and initial conditions, these equations constitute mathematical models that can be solved and studied in a way somewhat similar to performing experiments in a laboratory. Whether a model is mathematical or analogue, both are simplified abstractions of reality. However, such models are useful because they can help isolate the influence of certain parameters or scenarios, study complex system interactions, and make predictions.

An example of a mathematical model that has important application in Earth science is the heat conduction equation, often more generally referred to as the diffusion equation. A complete derivation of the heat conduction equation is given in Appendix A. In one dimension (1D), the heat conduction equation can be written as follows:

Here, T is the temperature (K), x is the distance (m), t is the time (s), ? is the rock density (kg m), c is the specific heat capacity (J kg K), k is the thermal conductivity (W m), and A is the rate of internal heat production per unit volume (J s m). In Equation 1.1, the temperature (the unknown) is referred to as the dependent variable, while t and x are known as independent variables. This type of equation is called a "partial differential equation" since the dependent variable depends on more than one independent variable. The physical parameters ?, c, k, and A are assumed to be known. Obtaining a solution to the equation means finding the function T(x, t) (i.e., T as a function of x and t) that satisfies the PDE.

More generally, the heat equation just introduced is also referred to mathematically as a parabolic (initial value) problem, which are typically of the form

Parabolic equations involve time-dependent behavior (term 1) and dissipation (terms 2 and 3), together which tend to smooth the solution with increasing time (at least for linear problems). Note that the signs in front of the second-order spatial derivatives on the right-hand side of 2 are necessarily positive; otherwise, the solutions grow rather than decay in time. Note also that the solution to parabolic equations depends on the initial value of the solution at (hence the name initial value problems). The other two major classes of PDEs are elliptic (boundary value) problems and hyperbolic. Elliptic equations are typically associated with steady-state problems. Examples of elliptic equations are Poisson's equation,

and Laplace's equation

which govern incompressible potential flow and steady heat transfer. Note that these equations don't involve any time derivatives and so their solutions depend only on the boundary conditions (hence the name boundary value problems) and any source (if present). Hyperbolic (initial value) PDEs involve time-dependent wave-like solutions. An example of a hyperbolic equation is the first-order wave equation

Here, the first term accounts for time-dependent behavior, while the second and third terms translate the solution laterally without any dissipation. Hyperbolic equations are common in problems involving flowing fluids.

1.2 Boundary and Initial Conditions

The solution to a PDE is not unique until boundary conditions are imposed. Boundary conditions essentially "ground" the solution to some specific physical scenario. There are four types of boundary conditions commonly encountered in the solution of PDEs:

1. Dirichlet, where the value of the solution is imposed on the boundary
2. Neumann boundary conditions, where the derivative of the solution is imposed on the boundary
3. Robin boundary conditions, where one specifies some linear combination of the solution and its derivative
4. Periodic (or repeating) boundary conditions, where one assumes that the solution at one end of the model domain is equal to the solution at the other end

The number of boundary conditions necessary to determine a solution to a differential equation matches the order of the highest spatial derivative in the differential equation. For example, Equation 1.1 contains a second-order spatial derivative and so two boundary conditions must be specified, one at each end of the domain. The equation also contains a first-order time derivative, so we must also provide an initial condition. This means we must define the value of T everywhere (over the entire domain) at . Equation 1.5 has only first-order spatial derivatives and so requires only one boundary condition in each direction. In this case, the boundary condition should be imposed at the end of the domain from where flow arrives, whereas the downstream end should be left unconstrained so that the flow can exit uninhibited.

1.3 Analytical Solutions

For relatively simple PDEs and for certain boundary conditions and initial conditions, it may be possible to find an exact (also known as a closed-form or analytical) solution. As an example, consider 1D heat transfer about a steadily creeping, narrow, planar, vertical fault. In this case, Equation 1.1 needs to be solved with A given by (e.g., see McKenzie and Brune, 1972)

where t is the (constant) shear stress (Pa) resolved on the fault plane, v is the fault slip rate (m ), and d(x0) is the Dirac function, that is, when , 0 when 0, and . The initial temperature at is assumed to be 0C everywhere. The spatial domain extends horizontally from to on either side of the fault located at . The boundary conditions are that the first derivative of the temperature vanishes at . The exact solution to Equation 1.1 combined with 1.6 can be written down directly using the Green's function for this equation (Morse and Feshbach, 1953, p. 981). The solution is

where ? ()) is the thermal diffusivity (m2) and erf is the error function . This solution can easily be evaluated exactly at any desired x and t once the values for the various physical parameters are specified (as done in the following).

1.4 Numerical Solutions

Although it is normally always desirable to obtain exact solutions to the PDE(s) being investigated, in practice this is often not possible. A closed-form solution may either not exist, or it may be too complicated to be of practical use. This may be due a number of factors, including nonlinearities in the governing equation, variable material properties, complicated geometries or boundary conditions, and so on. In such cases, one must resort to numerical methods that provide an approximate solution to the governing differential equation(s). Today, with powerful computers, many complicated problems can be solved quickly using numerical techniques.

The process of obtaining a computational solution consists of two stages shown schematically in Figure 1.1. The first stage converts the continuous PDE and auxiliary conditions (boundary and initial conditions) into a discrete system of algebraic equations. This first stage is called "discretization" and may be performed using various methods (one of which is the finite element method or FEM). The second stage involves solving the system of algebraic equations (normally performed on a computer, see Appendix B) to obtain an approximate solution to the original PDE. This second stage typically will involve some standard mathematical method such as Gaussian elimination.

Figure 1.1 Major steps involved in obtaining a numerical solution to a PDE.

Two important issues that must be considered when obtaining a numerical solution to PDEs are error and stability. All numerical methods introduce discretization errors, which in principle can be reduced by increasing the spatial and temporal resolution. This can be achieved by increasing the number of nodes (in time or space) where the solution is computed, or equivalently, by decreasing the spacing between nodes. It both cases, this should be performed without changing the total spatial or temporal extent of the model domain. Ideally, a numerical solution will converge to the exact solution as the resolution is increased. Even if an exact solution doesn't exist, one should always check that the numerical solution doesn't change significantly as the numerical resolution is changed, indicating that convergence has been achieved. Other errors may also...