# 20.4: Numerical Errors and Instability

- Page ID
- 9662

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)

\( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)

( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)

\( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)

\( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\)

\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)

\( \newcommand{\Span}{\mathrm{span}}\)

\( \newcommand{\id}{\mathrm{id}}\)

\( \newcommand{\Span}{\mathrm{span}}\)

\( \newcommand{\kernel}{\mathrm{null}\,}\)

\( \newcommand{\range}{\mathrm{range}\,}\)

\( \newcommand{\RealPart}{\mathrm{Re}}\)

\( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)

\( \newcommand{\Argument}{\mathrm{Arg}}\)

\( \newcommand{\norm}[1]{\| #1 \|}\)

\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)

\( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)

\( \newcommand{\vectorA}[1]{\vec{#1}} % arrow\)

\( \newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow\)

\( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

\( \newcommand{\vectorC}[1]{\textbf{#1}} \)

\( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)

\( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)

\( \newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)

Causes of NWP errors include **round-off error**,** truncation error**, **numerical instability**, and **dynamical instability**. Dynamical instability related to initial-condition errors will be discussed later in the section on chaos. Additional errors not considered in this section are coding bugs, computer viruses, user errors, and numerical or physical approximations (simplifications & parameterizations).

# 20.4.1. Round-off Error

**Round-off error** exists because computers represent numbers by a limited number of binary bits (e.g., 32, 64, 128 bits). As a result, some real decimal numbers can be only approximately represented in the computer. For example, a 32-bit computer can resolve real numbers that are different from each other by about 3x10^{–8} or greater. Any finer differences are missed.

To demonstrate, I wrote a computer program to start with x = 0.0, and then repeatedly add 0.1 to x (printing x at each step) until it reaches x = 3.0, at which point I programmed it to stop. When I used single precision (32-bits), my program never stopped. After 30 additions it had found x = 2.9999993, but since this was not exactly equal to 3.0, the program kept adding 0.1 in an **infinite loop** (i.e., ran forever). When I tried it again using double precision (64 bits) it also never stopped, getting only as close to 3.0 as x = 3.0000000000000013.

Namely, the slight error between decimal and binary representations of a number can accumulate and cause unexpected outcomes of conditional tests (“if” statements). Most modern computers use many bits to represent numbers. Nonetheless, always consider round-off errors when you write programs.

The first equations of fluid mechanics were formulated by Leonhard Euler in 1755, using the differential calculus invented by Isaac Newton in 1665 and Gottfried Leibniz in 1675, and using partial derivatives as devised by Jean le Rond d’Alembert in 1746.

Terms for molecular viscosity were added by Claude-Louis Navier in 1827 and George Stokes in 1845. The equations describing fluid motion are often called the **Navier-Stokes equations**. These primitive equations for fluid mechanics were refined by Herman von Helmholtz in 1888.

About a decade later Vilhelm Bjerknes in Norway suggested that these same equations could be used for the atmosphere. He was a very strong proponent of using physics, rather than empirical rules, for making weather forecasts.

In 1922, Lewis Fry Richardson in England published a book describing the first experimental numerical weather forecast — which he made by solving the primitive equations with mechanical desk calculators. His book was very highly regarded and well received as one of the first works that combined physics and dynamics in a thorough, interactive way.

It took him 6 weeks to make a 6 h forecast. Unfortunately, his forecast of surface pressure was off by an order of magnitude compared to the real weather. Because of the great care that Richardson took in producing these forecasts, most of his peers concluded that NWP was not feasible. This discouraged further work on NWP until two decades later.

John von Neumann, a physicist at Princeton University’s Institute for Advanced Studies, and Vladimir Zworykin, an electronics scientist at RCA’s Princeton Laboratories and key inventor of television, proposed in 1945 to initiate NWP as a way to demonstrate the potential of the recently-invented electronic computers. Their goal was to simulate the global circulation. During the first few years they could not agree on how to approach the problem.

Von Neumann formed a team of theoretical meteorologists including Carl-Gustav Rossby, Arnt Eliassen, Jule Charney, and George Platzman. They realized the need to simplify the full primitive equations in order to focus their limited computer power on the long waves of the global circulation. So Charney and von Neumann developed a simple one-layer barotropic model (conservation of absolute vorticity).

Their first electronic computer, the ENIAC, filled a large room at the Univ. of Pennsylvania. It used 18,000 vacuum tubes that generated tremendous heat and frequently burned out. The research team had to translate the differential equations into discrete form, write the code in machine language (FORTRAN, C and python had not yet been invented), decide the forecast domain size, and do many preliminary calculations using slide rules and mechanical calculators.

Their first ENIAC forecasts were made in MarchApril 1950, for three weather case studies over North America. This was the start of modern NWP.

# 20.4.2. Truncation Error

**Truncation error** was already discussed, and refers to the neglect (i.e., truncation) of higher-order terms in a Taylor series approximation to local gradients. If we retain more terms in the Taylor series, then the result is a higher-order solution that is more accurate, but which takes longer to run because there are more terms to compute. If we truncate the series at lower order, the numerical solution is faster but less accurate. In NWP, time and space difference schemes are chosen as a compromise between accuracy and speed.

# 20.4.3. Numerical Instability

**Numerical instability** causes forecasts to **blow up**. Namely, the numerical solution rapidly diverges from the true solution, can have incorrect sign, and can approach unphysical values (±∞). Truncation error is one cause of numerical instability.

Numerical instability can also occur if the wind speeds are large, the grid size is small, and the time step is too large. For example, eq. (20.16) models advection by using temperature in neighboring grid cells. But what happens if the wind speed is so strong that temperature from a more distant location in the real atmosphere (beyond the neighboring cell) can arrive during the time step ∆t? Such a physical situation is not accounted for in the numerical approximation of eq. (20.16). This can create numerical errors that amplify, causing the model to blow up (see Fig. 20.11).

Such errors can be minimized by taking a small enough time step. The specific requirement for stability of advection processes in one dimension is

\begin{align}\Delta t \leq \frac{\Delta X}{|U|}\tag{20.18}\end{align}

with similar requirements in the y and z directions. This is known as the **Courant-FriedrichsLewy (CFL) stability criterion**, or the Courant condition. When modelers use finer mesh grids with smaller ∆X values, they must also reduce ∆t to preserve numerical stability. The combined effect greatly increases model run time on the computer. For example, if ∆X and ∆Y are reduced by half, then so must ∆t, thereby requiring 8 times as many computations to complete the forecast.

For other physical processes such as diffusion and wave propagation, there are other requirements for numerical stability. To preserve overall stability in the model, one must satisfy the most stringent condition; that is, the one requiring the smallest time step. Some high-resolution NWP models use time steps of ∆t = 5 s or less.

For advection, one way to avoid the time-step limitation above is to use a **semi-Lagrangian** method. This scheme uses the wind at each grid point to calculate a **backward** **trajectory**. The backward trajectory indicates the source location for air blowing into the grid cell of interest. This source location need not be adjacent to the grid-cell of interest. By carrying the values of meteorological variables from the source to the destination during the time step, advection can be successfully modeled (i.e., be numerically stable) even for long time steps.

Theoretically, the smallest horizontal wavelength you can resolve with data at discrete grid points is 2·∆X. However, the finite-difference equations that are used to describe advection and other dynamics in NWP models are unable to handle 2·∆X waves. Namely, these waves either do not advect at all (Fig. 20.11d), or they are numerically unstable.

To avoid such unphysical behavior, small wavelength waves are numerically filtered out of the model. As a result, the smallest waves that are usually retained in NWP models are about 5 to 7·∆X.

Hence, the actual **resolution** (i.e., the smallest weather features that can be modeled) are about 7 times the **grid spacing**. Stated another way, if you know the size of the smallest weather system or terrain-related flow that you want to be able to forecast, then you need to design your NWP model with horizontal grid spacing ∆X smaller than 1/7 of that size.

**Sample Application**

What grid size, domain size, number of grid points, and time steps would you use for a numerical model of a hurricane, and how many computations would be needed to make a 3-day forecast? How fast should your computer be? [Hint: Use info from the Hurricane chapter.]

**Find the Answer**

This is an example of how you **design an NWP system**, including both the software and hardware. Assume the smallest feature you want to resolve is a thunderstorm in the eyewall. If tropical thunderstorms are about 14 km in diameter, then you would want __ ∆X__ = (14 km)/7 =

__to horizontally resolve it.__

**2 km**Hurricanes can be 300 km in diameter. To model the whole hurricane and a bit of its surrounding environment, you might want a horizontal __ domain of 500 km by 500 km__. This works out to (500 km / 2 km) = 250 grid points in each of the x and y directions, giving (250)2 = 62,500 grid points in the horizontal. If you want a model with 50 vertical levels, then you need (50) · (62,500) =

__.__

**3,125,000 grid points total**If you want to be able to forecast hurricanes up through category 5 (wind speed > 69 m s^{–1}), then design for a max wind of 80 m s^{–1}. The CFL criterion (eq. 20.18) gives ∆t = (2000 m)/(80 m s^{–1}) = 25 s. Thus, a 3-day forecast would require (72 h) · (3600 s h^{–1}) / (25 s) = __ 10,368 time steps__.

The temperature forecast eq (20.17) has about 43 arithmetic operations (adds, subtracts, multiplies, divides). We have 7 equations of motion, so this gives about (7 ·43 ≈) 300 operations. You must do these operations at each grid point for each time step, giving a total = 3,125,000 x 10,368 x 300 ≈ 10^{13} operations.

But we haven’t included the calculations for all the other physics (clouds, turbulence, precipitation, radiation, etc.) that must be done at each grid point. As a quick estimate, round up to 10^{15} **floating-point** (real-number) operations.

But you need to complete all these calculations quickly, in order to be useful as a forecast to warn people to evacuate. Suppose you design the model to finish within 3 h (=10,800 s) of computer run time. Thus, your computer must be powerful enough to compute at the rate of (10^{15} operations)/(10,800 s) ≈ 10^{11} = __ 100 gigaflops__ (where 1

**flops**= 1

**floating-point operation per second**).

**Check:** Units OK. Physics OK.

**Exposition:** The number of calculations needed to make a hurricane forecast is tremendous, and requires a powerful computer. As computer power increases, NWP modelers strive for finer horizontal and vertical resolution spanning larger domains, and including more accurate (and complex) physics parameterizations.

Namely, NWP modelers always tend to fully use all the computer power available, and dream of even more powerful computers.