24 Basic concepts of numerics
A weather model is a copy of the atmosphere. In it, all fields and differential operators are replaced by discretized versions, since you can only store a finite amount of information. One writes a substitution for a differential operator $D$
\[ \begin{align} D\to D', \end{align} \]
analogous for atmospheric conditions
\[ \begin{align} Z\to Z'. \end{align} \]
Then you define a dynamic, i.e. a procedure
\[ \begin{align} Z'\left(t\right)\stackrel{D'}{\to}Z'\left(t + \Delta t\right) \end{align} \]
to determine the time evolution of the model atmosphere with the aim of bringing the state trajectory of the model as close as possible to that of the real atmosphere. To do this, you define a distance
\[ \begin{align} \left|Z - Z'\right|\geq 0. \end{align} \]
Weather forecast has the following workflow:
Observation of the state of the atmosphere.
Determination of the initial state of the model.
Integration of the model.
Errors arise for the following reasons:
The initial state was not fully observed.
The initial state was incorrectly observed.
The governing equations are equations of statistical physics and therefore do not apply exactly.
The model can store less information than the real atmosphere.
The dynamics of the model are different from those of the real atmosphere.
There are rounding errors due to the finite computer precision.
The computer is making errors.
24.1 Basic concepts of discretizations
Let the solution of a differential equation be denoted by $F$ and the solution of a discretized version of the equation by $F_d$. A numerical procedure is called consistent if and only if the solution of the scheme converges to the solution of the differential equation as $\Delta x$ and $\Delta t$ become smaller, i.e. $\lim\limits_{\Delta\to 0}F_d = F$ in the above terms. A discretization is called self-consistent if all conserved quantities of the continuous system are also conserved quantities of the discretized system. This can also be transferred to approximations. A scheme is called stable if and only if the calculated solution is bounded for all times, i.e. if a $C>0$ exists with $|F_d|<C$ for all $t\geq t_0$ with $t_0$ as the starting time. Linear instability already arises with linear differential equations, whereas nonlinear instability can only arise with nonlinear differential equations.
24.2 Finite differences
Let $f:\mathbb{R}\to\mathbb{R}$ be an infinitely differentiable function. Using a Taylor series, one can write
\[ \begin{align} f\left(x + \Delta x\right) = f\left(x\right) + \Delta xf'\left(x\right) + \frac{f''\left(x\right)}{2!}\Delta x^2 + \frac{f'''\left(x\right)}{3!}\Delta x^3 + \dotsc\tag{24.5}\label{eq:taylor_links} \end{align} \]
and
\[ \begin{align} f\left(x - \Delta x\right) = f\left(x\right) - \Delta xf'\left(x\right) + \frac{f''\left(x\right)}{2!}\Delta x^2 - \frac{f'''\left(x\right)}{3!}\Delta x^3 + \dotsc\tag{24.6}\label{eq:taylor_rechts} \end{align} \]
Subtracting both equations yields
\[ \begin{align} f\left(x + \Delta x\right) - f\left(x - \Delta x\right) = 2\Delta xf'\left(x\right) + R', \end{align} \]
where the remainder $R'$ contains terms of third and higher order. Rearranging gives
\[ \begin{align} f'\left(x\right) = \frac{f\left(x + \Delta x\right) - f\left(x - \Delta x\right)}{2\Delta x} + R, \tag{24.8}\label{eq:approx_derivative_zentral} \end{align} \]
where $R$ is a polynomial in $\Delta x$ whose lowest power is 2; thus, Eq. (24.8) is a second-order approximation of the derivative, and one writes the symbol $\mathcal{O}\left(\Delta x^2\right)$ instead of $R$. This way of approximating derivatives is called a central difference quotient. A one-sided difference quotient, for example,
\[ \begin{align} f'\left(x\right) \approx\frac{f\left(x + \Delta x\right) - f\left(x\right)}{\Delta x} \end{align} \]
is also a suitable approximation, but it converges only to first order, as is easily seen by rearranging Eq. (24.5). The higher the order, the better, since the accuracy of a Taylor expansion increases with order. For this reason, central spatial difference quotients are used wherever possible. One can also derive an approximation for the second derivative of a function from Eqs. (24.5) and (24.6):
\[ \begin{align} f\left(x + \Delta x\right) + f\left(x - \Delta x\right)&= 2f\left(x\right) + f''\left(x\right)\Delta x^2 + \mathcal{O}\left(\Delta x^4\right)\nonumber\\ \Leftrightarrow f''\left(\Delta x\right) &= \frac{f\left(x + \Delta x\right) - 2f\left(x\right) + f\left(x - \Delta x\right)}{\Delta x^2} + \mathcal{O}\left(\Delta x^2\right) \end{align} \]
For example, for $f\left(x\right) = A\sin\left(kx\right)$ one has $f' = Ak\cos\left(kx\right)$, but Eq. (24.8) gives
\[ \begin{align} f' &\approx& \frac{A}{2\Delta x}\left(\sin\left(kx + k\Delta x\right) - \sin\left(kx - k\Delta x\right)\right) = \frac{A}{\Delta x}\cos\left(kx\right)\sin\left(k\Delta x\right) = \frac{\sin\left(k\Delta x\right)}{k\Delta x}f'. \end{align} \]
This shows that the approximation converges, as expected, to the correct derivative as $\Delta x$ approaches zero (since $\lim\limits_{\alpha\to 0}\frac{\sin\left(\alpha\right)}{\alpha} = 1$ by L'Hospital's rule). Writing $k = \frac{2\pi}{\lambda}$ with $\lambda$ as wavelength and setting $\Delta x = \lambda/2$, the central difference quotient always yields $f' = 0$. However, sufficiently long waves ($\lambda\gg 2\Delta x$) are well resolved.
24.3 Plane spectral method
The plane spectral method rests on Fourier-transforming the horizontal dependencies of the state variables. This transformation is set out in Section C.1. The quantities $x,y$ may be any horizontal coordinates, such as geographic longitude and latitude, the spherical coordinates of a rotated longitude-latitude grid, or coordinates defined within a map projection such as the Lambert projection. A well-known dynamical core that employs this method is ALADIN, which underlies the models AROME and HARMONIE.
24.4 Spherical spectral method
Let $\psi = \psi\left(\varphi, \lambda, z\right)$ be a scalar field. Since the spherical harmonics presented in Sect. C.5 form a complete orthonormal basis on spherical shells, the horizontal dependence of $\psi$ can be represented by this set of functions:
\[ \begin{align} \psi\left(\varphi, \lambda, z\right) = \sum_{l = 0}^\infty\sum_{m = -l}^l\newtilde{\psi}_{l, m}\left(z\right)Y_{l, m}\left(\varphi, \lambda\right) \end{align} \]
The definition of the spherical harmonics, Eq. (C.155), is written here in geographic coordinates:
\[ \begin{align} Y_{l, m}\left(\varphi, \lambda\right) = \sqrt{\frac{2l + 1}{4\pi}\frac{\left(l - m\right)!}{\left(l + m\right)!}}P_{l, m}\left(\sin\left(\varphi\right)\right)\exp\left(im\lambda\right) \end{align} \]
