# NUMERICAL MODELS | Methods

## NUMERICAL MODELS | Methods

Methods J Thuburn, University of Exeter, Exeter, UK Ó 2015 Elsevier Ltd. All rights reserved. Synopsis Numerical modeling of the atmosphere is challe...

Synopsis Numerical modeling of the atmosphere is challenging because of the enormous range of space and timescales involved, and because many qualitative features, such as balance and conservation, must be accurately captured. We discuss here a number of topics that have proved important historically in the development of numerical models of the atmosphere, including spatial discretization and its relation to representation of model ﬁelds, choice of vertical and horizontal grids, time integration schemes, advection schemes, and solution methods for elliptic problems.

temperature, etc.) must be stored in a computer as a ﬁnite set of values. There are several possible ways to do this.

Special Challenges in Numerical Modeling of the Atmosphere The atmosphere has a number of characteristics that make numerical solution of its governing equations especially challenging. Most obvious is the geometry. The atmosphere is highly anisotropic, with vertical length scales being typically much smaller than horizontal. Gravity acts in the vertical direction, and the atmosphere is strongly stratiﬁed and close to hydrostatic balance. The ratio of horizontal to vertical model resolution is commonly chosen to capture this anisotropy, and model grids and methods must be able to represent hydrostatic balance (see Section on Grids below). For global modeling, the grid and methods chosen must be able to represent the (approximately) spherical geometry of the Earth (again see Section on Grids). The atmosphere is strongly multiscale in both space and time. In terms of spatial energy spectra, the largest scales are energetically dominant, but spectra are rather shallow, implying that, whatever the resolution of an atmospheric model, there is signiﬁcant dynamical variability near the resolution limit, which requires careful handling by numerical methods, and on unresolved scales, which must be represented by ‘subgrid models.’ The atmosphere also supports dynamics with a huge range of timescales. Care is needed in modeling fast processes to ensure that the numerical solution is stable (see Section on Time Integration Schemes, Stability, Courant–Friedrichs–Lewy Criterion below). On the other hand, fast waves (acoustic waves, and to a lesser extent, inertio-gravity waves) tend to be energetically weak, implying that the dynamics is often slowly evolving and close to hydrostatic and geostrophic balance. The chosen numerical method must be able to represent this balance accurately. Also, certain properties of the atmosphere evolve slowly, either in a Lagrangian sense (e.g., the moisture content of an air parcel in the absence of condensation and evaporation) or in a global integral sense (e.g., the total mass or energy of the atmosphere). Numerical methods should be able to capture these conservation properties accurately (see Section on Conservation and Scale-Selective Dissipation below).

Representation of Data In a numerical model, the ﬁelds describing the threedimensional state of the atmosphere (density, velocity,

Encyclopedia of Atmospheric Sciences 2nd Edition, Volume 4

Gridpoint representation. In gridpoint methods, a ﬁeld is represented by its values sampled at a ﬁnite set of points called the grid or mesh points (Figure 1, top left). l Finite volume representation. Each value stored in the computer represents the value of a ﬁeld averaged over a ﬁnite volume or grid cell (Figure 1, top right). l Finite element representation. Each ﬁeld f(x) is represented as a series in terms of certain basis functions ji ðxÞ, and b i are stored in the computer (Figure 1, the coefﬁcients f bottom left): l

fðxÞ ¼

X

b i ji ðxÞ: f



i

Usually the ji are simple local functions such as piecewise polynomials that are nonzero only over a small number of grid cells. l Spectral representation. As for ﬁnite element methods, each ﬁeld is represented as a series in terms of certain basis functions, but here the basis functions are mutually orthogonal global functions. The simplest example in one dimension is a Fourier series representation, in which the basis functions are sines and cosines (Figure 1, bottom right). For modeling the atmosphere on a spherical planet, a spectral method based on spherical harmonics is often used. Some methods may be viewed in more than one way, for example as a ﬁnite volume representation for the purpose of looking at conservation properties but as a gridpoint representation for the purpose of calculating derivatives. In more than one dimension, combinations of the above methods can be used, e.g., spectral in the horizontal and gridpoint in the vertical. The number of data values used to represent a ﬁeld deﬁnes the model resolution: a greater number of data values corresponds to higher resolution and means that smaller-scale structures in the ﬁeld can be represented. For a current-day global weather forecast model, a typical horizontal resolution might be a few tens of kilometers, while a typical vertical resolution might stretch from a few meters or a few tens of meters near the surface to 2 or 3 km in the stratosphere. Climate models generally have coarser horizontal resolution to make long integrations affordable.

http://dx.doi.org/10.1016/B978-0-12-382225-3.00246-2

161

162

Numerical Models j Methods

Gridpoint

Finite volume

2

2

1

1

0

0

−1

0

0.2

0.4

0.6

0.8

1

−1

0

0.2

0.4

Finite element

0.6

0.8

1

0.6

0.8

1

Spectral

2

1 0.5

1

0 0

−1

−0.5

0

0.2

0.4

0.6

0.8

1

−1

0

0.2

0.4

Figure 1 Illustration of different discrete methods for representing the same one-dimensional function. Top left: representation by ﬁve gridpoint values (the sixth value is the same as the ﬁrst because the function is periodic). Top right: ﬁnite volume representation as ﬁve grid cell averages. Bottom left: ﬁnite element representation in terms of ﬁve piecewise linear ‘witches hat’ basis functions. Bottom right: the ﬁrst ﬁve Fourier components of the spectral representation of the function.

Round-Off Error, Truncation Error, Local Truncation Error, Consistency Numerical solution of differential equations involves approximations and hence, inevitably, errors. Round-off error arises because a ﬁnite number of digits, and hence a limited precision, are used in storing numbers in a computer. On a typical workstation using standard precision, the round-off error is about 3 in the eighth decimal digit. For most purposes this error is so small as to be insigniﬁcant. However, the effects of round-off error can be greatly ampliﬁed in certain calculations; for example, when adding small numbers to large ones or when taking differences of nearly equal numbers. The approximation of a continuous ﬁeld by a ﬁnite set of data values and the approximation of derivatives or integrals by algebraic operations such as differences or sums is called discretization. The errors made in these approximations are called truncation errors. They are a function of the resolution of the model, and are usually much more signiﬁcant than round-off errors. The local truncation error is the residual obtained when the true solution of a differential equation is substituted into the discrete approximation of that equation. It is a measure of how well the discretization approximates the differential equation. For example, consider the one-dimensional linear advection equation vf vf þv ¼ 0; vt vx



where the advecting velocity v is a constant. Let us discretize on a uniform grid with spatial interval Dx and time interval Dt so that fnj is the numerical approximation to f at x ¼ jDx, t ¼ nDt.

Approximating the space and time derivatives by simple difference formulas gives, for example, fnþ1  fnj j Dt

þv

fnj  fnj1 Dx

¼ 0:



If we substitute the true solution f(jDx, nDt) in place of fnj , etc., in eqn  and make use of Taylor series, we ﬁnd the residual, i.e., the local truncation error, as O(Dx) þ O(Dt). In this example the discretization is said to be ﬁrst-order accurate in space and time, because both Dx and Dt appear raised to the power one in the local truncation error. A discretization is said to be consistent if the local truncation error tends to zero as Dx / 0 and Dt / 0.

Spatial Discretization Methods: Finite Difference, Finite Volume, Finite Element, and Spectral A huge variety of methods are available for approximating the spatial derivatives that appear in the governing equations. The different methods are related to how the data are represented. Finite difference methods are perhaps conceptually the simplest. Rearranging Taylor series allows the construction of approximate expressions for derivatives in terms of point values of the data. Equation  shows a simple example in which vf/vx at (x ¼ jDx, t ¼ nDt) is approximated by ðfnj  fnj1 Þ=Dx. Higher-order accurate formulas and formulas for higher derivatives can also be derived. Finite volume methods allow one to deﬁne a local or global integral of a predicted ﬁeld, and so are natural when the emphasis is on integral conservation properties. The equations

Numerical Models j Methods are usually written in terms of the budget of each prognostic ﬁeld for each grid cell, and the method is deﬁned by constructing a suitable approximation for the ﬂux of each ﬁeld across each cell face. For ﬁnite element and spectral methods, the objective is b i evolve in time. One way to determine how the coefﬁcients f to do this is to require that the residual, that is, the error when the approximate solution  is substituted into the original differential equation, should be orthogonal to each of the basis functions. This approach typically leads to a set b i =dt multiplied by of simultaneous equations involving d f a sparse matrix called the mass matrix. Once these equations b i may be stepped forward using a b i =dt, f are solved for d f suitable time integration scheme. Discontinuous Galerkin and spectral element methods are closely related to the ﬁnite element method and have also been used in atmospheric models. For spectral methods, the basis functions are mutually orthogonal; therefore the mass matrix is diagonal and trivial to invert. On the other hand, nonlinear terms in the governing equations would be expensive to evaluate directly in terms of b i . Instead, it is simpler and more the spectral coefﬁcients f efﬁcient to transform the model ﬁelds from a spectral repreb i to a gridpoint representation in terms sentation in terms of f of fj to evaluate the nonlinear terms, and then to transform the result back into the spectral representation. The transforms can be carried out efﬁciently with the aid of a fast Fourier transform algorithm. This method is called the spectral transform method.

Grids Because of the strong anisotropy of the atmosphere, atmospheric models are almost invariably constructed with a certain number of levels or layers in the vertical, with essentially the same horizontal grid structure at each level. Thus it makes sense to discuss the vertical grid structure and horizontal grid structure separately. An obvious choice for the locations of the vertical levels is to place them at speciﬁed heights. However, in the presence of mountains, the bottom boundary condition is greatly simpliﬁed if the model levels are distorted to follow the terrain (Figure 2). It is common to deﬁne the levels to be more closely spaced near the bottom boundary, to better capture details of the atmospheric boundary layer.

Figure 2 Schematic showing three variants of a height-based coordinate and levels and how they vary in the vicinity of a mountain. Left: ﬂat, terrain-intersecting levels; middle: basic terrain-following levels; right: hybrid levels that are terrain-following near the surface but become ﬂat at high altitude.

163

The governing equations may be expressed in a variety of different vertical coordinates, and model levels chosen at speciﬁed values of the vertical coordinate. As well as height, example vertical coordinates include pressure, mass, potential temperature, and a Lagrangian vertical coordinate (i.e., one that moves with the ﬂow), and their corresponding terrain-following variants. Each has its own advantages and disadvantages compared with a height-based coordinate and levels. More accurate wave propagation and a more accurate representation of hydrostatic balance can be achieved by using a staggered vertical grid, in which different variables are interleaved rather than colocated. Two types of staggered vertical grid are widely used: those with the temperature variable at the same levels as the horizontal velocity components, commonly referred to as the Lorenz grid; and those with the temperature variable staggered with respect to the horizontal wind components, commonly referred to as the Charney–Phillips grid (Figure 3). The Charney–Phillips grid provides slightly more accurate wave propagation, while energy conservation is easier to achieve on the Lorenz grid. For a limited area model, the simplest choice for the horizontal grid is a rectangular Cartesian grid. For a global model, a longitude–latitude grid is the simplest choice. However, the clustering of resolution near the poles creates additional difﬁculties related to stability, convergence of elliptic solvers (see below), and scalability on massively parallel computers. This has motivated the development of models with a variety of more uniform horizontal grids. Some examples are shown in Figure 4. Staggered horizontal grids are possible, and can give increased accuracy in representing wave propagation and balance. Some examples of staggered rectangular Cartesian grids are shown in Figure 5. The commonest staggered grids are called the Arakawa A-, B-, C-, D-, and E-grids, after the work of Arakawa and colleagues who ﬁrst studied their wave dispersion properties. Analogous staggering of variables is possible on more general polygonal grids such as those shown in Figure 4.

Time Integration Schemes, Stability, Courant– Friedrichs–Lewy Criterion Given a solution at time t ¼ nDt (and perhaps at earlier time steps), we need a scheme to advance the solution to time

w

w, θ

u, v, ρ, θ

u, v, ρ

w

w, θ

u, v, ρ, θ

u, v, ρ

w

w, θ

Figure 3 Schematic showing the vertical placement of the prognostic ﬁelds for examples of the Lorenz vertical grid (left) and the Charney– Phillips vertical grid (right). u, v, and w are the velocity components, r is density, and q is potential temperature.

164

Numerical Models j Methods

Figure 4 Some examples of spherical grids used in atmospheric models. Top left to bottom right: Longitude–latitude grid, Yin-Yang grid, equal angle cubed sphere, conformal cubed sphere, triangular icosahedral grid, hexagonal–pentagonal icosahedral grid.

A

B

C u v

ρ , u, v

u v

ρ u v

v u

u v

ρ

u

v

Figure 5 Schematic showing the horizontal placement of prognostic ﬁelds for three of the rectangular grids studied by Arakawa and colleagues. Left: A-grid; middle: B-grid; right: C-grid. u, v are the horizontal velocity components and r is the density.

t ¼ (n þ 1)Dt. A great many such schemes are possible. Consider the following differential equation to illustrate ideas: df ¼ RðfÞ:  dt One of the simplest possible schemes approximates the time derivative by a forward difference: fnþ1  fn ¼ Rðfn Þ: Dt



Equation  uses this time integration scheme. This is an example of an explicit scheme, so-called because the formula

may be rearranged to give an explicit expression for fnþ1 in terms of known quantities. The scheme (5) is not very accurate. Greater accuracy can sometimes be obtained by using information from other time levels, giving a so-called multistep method. The leapfrog scheme fnþ1  fn1 ¼ Rðfn Þ:  2Dt is one example. Adams–Bashforth schemes also fall into this category. Another way to increase accuracy is via a multistage scheme. Here, intermediate solutions are computed at one or more

Numerical Models j Methods stages between step n and step n þ 1. The Matsuno or forward– backward scheme is one example: f  fn ¼ Rðfn Þ; Dt

fnþ1  fn ¼ Rðf Þ: Dt



Runge–Kutta schemes are also popular multistage schemes. An important consideration for time integration schemes is stability. There are several slightly different technical deﬁnitions of stability. For practical purposes, a scheme is unstable if it produces numerical solutions that grow in amplitude when the true solution does not grow. Explicit time integration schemes are typically stable only for sufﬁciently small Dt, that is, they are conditionally stable. Stable solutions for larger time steps can be obtained using implicit time integration schemes. A simple example is the trapezoidal or Crank–Nicolson scheme   fnþ1  fn 1 ¼ Rðfn Þ þ R fnþ1 : 2 Dt



Such schemes are called implicit because the unknowns fnþ1 appear in several places, and a system of (possibly nonlinear) simultaneous equations must be solved to step forward. See Section on Elliptic Problems below. Stability often depends on the coupling of the spatial discretization and the time integration scheme, in particular on the satisfaction of the Courant–Friedrichs–Lewy (CFL) criterion. The CFL criterion states that a necessary condition for stability is that the numerical domain of dependence of the solution at any point in space and time should contain the physical domain of dependence. For simple schemes, such as (3), the CFL criterion often reduces to the requirement that the Courant number cDt/Dx, where c is an appropriate wave or advective signal velocity, should be less than one. Various issues must be taken into consideration when choosing a time integration scheme (in combination with a spatial discretization). The scheme should, of course, be stable, but at the same time it should not introduce excessive dissipation or damping of the solution. Phase and group propagation of all waves of meteorological interest should be accurately captured. Even waves that are not of direct interest (often the case for acoustic waves and sometimes inertiogravity waves) must be captured well enough to represent adjustment toward balance. Multistep schemes may support computational modes, that is, numerical solutions that do not correspond to any solution of the original differential equation; some form of time ﬁltering may be necessary to control them. Finally, the overall cost should be acceptable. In practice, time integration schemes are often complex combinations of the kinds of scheme mentioned above. For example, early atmospheric models often used a leapfrog scheme (6) for the conservative wavelike terms combined with a forward scheme (like (5), but from step n  1 to step n þ 1) for dissipative terms, using each scheme on the terms for which it is conditionally stable. To avoid solving a global nonlinear system of equations, semi-implicit schemes treat a linearized subset of terms in the governing equations implicitly (like eqn ) and the remaining terms explicitly. Horizontally explicit, vertically implicit schemes treat only vertical derivative terms implicitly, so that only a relatively small system of equations needs to be solved for each vertical column. Finally, split explicit schemes are also popular. These avoid the need to solve

165

a system of simultaneous equations that arises with an implicit method while reducing computational cost by taking small time steps only for a subset of (relatively cheap) terms associated with fast waves and using longer time steps for the other (more expensive) terms.

Advection The transport of some ﬂuid property by ﬂuid parcels is called advection. The importance of advection for ﬂuid problems, and, in particular, the fact that many important quantities, such as potential temperature, potential vorticity, moisture, and long-lived chemical tracers often evolve slowly following ﬂuid parcels – that is, they are approximately conserved in a Lagrangian sense – have led to a great deal of effort to develop accurate numerical methods for advection. One approach to modeling advection is to combine an accurate approximation for the spatial derivatives (e.g., the vf/vx term in eqn ) with an accurate approximation for the time derivative term (e.g., the vf/vt term in eqn ). In this approach, better stability and accuracy is often obtained by using an upwind biased approximation for the spatial derivatives. An alternative approach is to re-express the advection equation as Df ¼ 0;  Dt where D/Dt indicates a derivative following a ﬂuid parcel rather than at a ﬁxed point in space. Semi-Lagrangian schemes follow this approach. For each gridpoint, xj say, they ﬁrst compute a departure point xd, which is the location at step n of the ﬂuid parcel arriving at xj at step n þ 1. The value of f at the departure point fnd is computed by interpolating the ﬁeld fn. Then the arrival point value is set equal to the departure point ¼ fnd . Semi-Lagrangian schemes are widely used value: fnþ1 j in atmospheric models because of their high accuracy, and because their stability is not limited by the advective Courant number, allowing longer time steps and making them very efﬁcient. It is often desirable to ensure that the numerical advection scheme does not lead to unphysical negative values of the advected ﬁeld, or to spurious ampliﬁcation of extrema (which could lead to spurious supersaturation in the moisture ﬁeld, for example). Such a desirable property goes by a variety of names, including monotone, monotonicity preserving, total variation diminishing, shape preserving, or locally bounded. It can be achieved by suitable modiﬁcations to interpolation schemes or to ﬂuxes, often referred to as limiters.

Elliptic Problems Elliptic problems arise in atmospheric models when some form of balance is assumed, e.g., when determining the stream function from the quasigeostrophic potential vorticity, or when determining the pressure ﬁeld for the Boussinesq or anelastic equations. They also arise when vertical vorticity and horizontal divergence are the predicted variables, and it is necessary to diagnose the velocity ﬁeld from them via a stream function and velocity potential. Finally, the system of equations that

166

Numerical Models j Methods

arises for the unknowns when using an implicit or semiimplicit time integration scheme typically takes the form of a spatial discretization of a Helmholtz problem   fnþ1  V\$ nVfnþ1 ¼ RHS;  where RHS indicates known terms. Spatial discretization of an elliptic problem leads to an equation of the form Mfnþ1 ¼ B;



where M is a sparse matrix and B is a known vector. Direct solution of eqn  is usually prohibitively expensive for three-dimensional problems. However, other methods are available. If a spectral representation of f was used and the basis functions could be chosen to be eigenfunctions of the Helmholtz operator, then M would be diagonal and the solution would be trivial. In practice the basis functions are usually chosen to be eigenfunctions of the horizontal part of the Helmholtz operator (for n constant), e.g., spherical harmonics for a global model. Then eqn  separates into a number of much smaller problems whose matrices correspond to the vertical part of the Helmholtz operator; these can be solved by direct matrix inversion. If a spectral method is not used, then eqn  may be solved iteratively. One type of iterative method improves (‘relaxes’) the solution at each gridpoint in turn, making a number of sweeps over the whole grid. Depending on the details of the relaxation, this gives rise to Jacobi, Gauss– Seidel, and successive over-relaxation methods. Convergence of these basic relaxation methods may be slow if the length scale n1/2 in eqn  is large compared with the grid length. In this case, convergence can be greatly accelerated by using a multigrid method, in which some iterations are carried out on coarser grids than that on which the ﬁnal solution is required. Another type of iterative method ﬁnds a sequence of successively better solutions over a growing hierarchy of vector spaces, called Krylov subspaces, whose bases are generated by iteratively applying M to some starting vector. Depending on the details, this gives rise to methods such as conjugate gradient, biconjugate gradient, minimum residual, generalized minimum residual, and generalized conjugate residual with k steps. The convergence of Krylov subspace methods can often be greatly improved by use of an appropriate preconditioner. That is, the problem  is rewritten as

or as

P1 Mfnþ1 ¼ P1 B



  MQ1 Qfnþ1 ¼ B



for some matrix P or Q whose inverse is easily calculated.

Conservation and Scale-Selective Dissipation Good Lagrangian conservation properties can often be achieved in a numerical model by the choice of an accurate

advection scheme, as discussed above. Depending on the application, integral conservation properties may also be important. Conservation of mass of air, mass of moisture and chemical tracers, energy, angular momentum, entropy, and potential enstrophy (among others) have all been addressed by model developers. Conservation of mass is straightforwardly achieved by the use of a ﬁnite volume method for the density or tracer density equation. Other conservation properties typically require some special measures in the numerical methods, and may only be achieved approximately. A variety of processes can lead to the growth of variability or ‘noise’ in a numerical model on marginally resolved scales. These include physically realistic processes such as the transfer of energy and potential enstrophy to different length scales through the nonlinearity of the governing equations, as well as numerical errors such as dispersion errors and aliasing (the misrepresentation of variability at one, unresolved, scale as variability at another, resolved, scale). Other contributors include the action of physical parameterization schemes, and small-scale forcing, for example, due to orography. To control such a build up of noise, all atmospheric models include some form of scale-selective dissipation that preferentially damps small scales. This may be in the form of extra dissipation terms explicitly added to the governing equations, or in the form of inherently dissipative numerical approximations. Care is needed to ensure that such dissipation is sufﬁcient to control the buildup of grid-scale noise without leading to excessive damping or adversely affecting conservation properties.

See also: Dynamical Meteorology: Acoustic Waves; Balanced Flow; Lagrangian Dynamics; Overview; Potential Vorticity; Primitive Equations. Gravity Waves: Buoyancy and Buoyancy Waves: Theory. Numerical Models: General Circulation Models; Mesoscale Atmospheric Modeling; Regional Prediction Models. Turbulence and Mixing: Turbulence, Two Dimensional.

Further Reading Durran, D.R., 1999. Numerical Methods for Wave Equations in Geophysical Fluid Dynamics. Springer-Verlag, New York. Haltiner, J.G., Williams, R.T., 1980. Numerical Prediction and Dynamic Meteorology. Wiley, Chichester. Lauritzen, P.H., Jablonowski, C., Taylor, M.A., Nair, R.D., 2011. Numerical Techniques for Global Atmospheric Models. Springer-Verlag, Berlin, Heidelberg. Pielke, R.A., 1984. Mesoscale Meteorological Modeling. Academic Press, London. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in FORTRAN 77: The Art of Scientiﬁc Computing, second ed. Cambridge University Press, Cambridge. Staniforth, A., Thuburn, J., 2012. Horizontal grids for global weather and climate prediction models: a review. Quarterly Journal of the Royal Meteorological Society 138, 1–26. Thuburn, J., 2008. Some conservation issues for the dynamical cores of NWP and climate models. Journal of Computational Physics 227, 3715–3730.