25 Time step method

Let $\psi: A \subseteq \mathbb{R}^3 \times \mathbb{R} \to \mathbb{R}^n$ be the atmospheric state vector. All fundamental physical laws are first order in time, therefore one can write

\[ \begin{align} \frac{d\psi}{dt} = F\left(\psi\right), \end{align} \]

where $F$ contains the physics. Discretizing this onto a temporal grid $n\Delta t$ with $n \in \mathbb{Z}$ and interpolating between two time steps $n$, $n + l$, one obtains

\[ \begin{align} \psi_{n + l} - \psi_{n} = \int_{n\Delta t}^{\left(n + l\right)\Delta t}F\left(\psi\right)dt. \end{align} \]

The number of time steps involved is $l + 1$. The integral on the right-hand side can now be discretized in a variety of ways. The explicit Euler method is obtained, for $l = 1$, via the approximation

\[ \begin{align} \int_{n\Delta t}^{\left(n + 1\right)\Delta t}F\left(\psi\right)dt \approx F\left(\psi_{n}\right)\Delta t \end{align} \]

with $\psi_{n} \coloneqq \psi\left(n\Delta t\right)$. The implicit Euler method, on the other hand, is obtained via

\[ \begin{align} \int_{n\Delta t}^{\left(n + 1\right)\Delta t}F\left(\psi\right)dt \approx F\left(\psi_{n + 1}\right)\Delta t \end{align} \]

This chapter aims to explain how the time discretization of a linear differential equation can affect the dispersion relation and the stability of a wave solution. The linear advection equation

\[ \begin{align} \frac{\partial f}{\partial t} + U\frac{\partial f}{\partial x} = 0\tag{25.5}\label{eq:test_lineare_instab} \end{align} \]

considered with $U\not = 0$ as a homogeneous advection velocity and $f$ as some particle property. In complex notation, the analytical solution reads

\[ \begin{align} f = \exp\left(i\left(kx - \omega t\right)\right), \end{align} \]

Inserting this yields $-i\omega + Uik = 0$, i.e. as dispersion relation

\[ \begin{align} \omega = Uk. \end{align} \]

So phase and group velocity are equal to $U$. Now space and time are discretized; one writes $f_l^{(m)} = f\left(l\Delta x, m\Delta t\right)$ and makes the ansatz

\[ \begin{align} f_{l}^{(m)} = e^{Cm\Delta t}e^{i\left(kl\Delta x - \omega m\Delta t\right)}\tag{25.8}\label{eq:ansatz_time_stepping_1d_analysis} \end{align} \]

with $f_{l}^{(0)} = e^{ikl\Delta x}$ as the initial state. Here $C$, $k$ and $\omega$ are real numbers. For reasons of symmetry, $k$ cannot have an imaginary part. $C$ and the dispersion $\omega\left(k\right)$ are to be determined.

25.1 Purely explicit methods

25.1.1 Explicit Euler method

First, the following discretization is examined:

\[ \begin{align} \frac{f_l^{(m + 1)} - f_l^{(m)}}{\Delta t} + U\frac{f_{l + 1}^{(m)} - f_{l - 1}^{(m)}}{2\Delta x} = 0 \end{align} \]

This is called the explicit Euler method. Insertion yields

\[ \begin{align} \frac{e^{C\Delta t}e^{-i\omega\Delta t} - 1}{\Delta t} + U \frac{e^{ik\Delta x} - e^{-ik\Delta x}}{2\Delta x} = \frac{e^{C\Delta t}\cos\left(\omega\Delta t\right) - ie^{C\Delta t}\sin\left(\omega\Delta t\right) - 1}{\Delta t} + \frac{iU\sin\left( k\Delta x\right)}{\Delta x} = 0. \end{align} \]

The real part is

\[ \begin{align} e^{C\Delta t}\cos\left(\omega\Delta t\right) = 1, \end{align} \]

from this it follows in general that $C> 0$; instability is present. This is linear instability, since Eq. (25.5) is linear in $f$. The imaginary part is

\[ \begin{align} e^{C\Delta t}\sin\left(\omega\Delta t\right) = \frac{U\Delta t\sin\left(k\Delta x\right)}{\Delta x}. \end{align} \]

From this it follows

\[ \begin{align} \tan\left(\omega\Delta t\right) = \frac{U\Delta t\sin\left(k\Delta x\right)}{\Delta x}, \end{align} \]

so the dispersion

\[ \begin{align} \omega\left(k\right) = \frac{1}{\Delta t}\arctan\left(\frac{U\Delta t\sin\left(k\Delta x\right)}{\Delta x}\right). \end{align} \]

Let $\lambda$ be the wavelength, i.e. $k = \frac{2\pi}{\lambda}$. For $\lambda\gg \Delta x$, one has $\sin\left(k\Delta x\right) \approx k\Delta x$, i.e. $\omega \approx \frac{1}{\Delta t}\arctan\left(Uk\Delta t\right)$. In the case $U\Delta t\ll \lambda$, i.e. $Uk\Delta t\ll 1$, $\omega \approx \frac{1}{\Delta t}U\Delta t k = Uk$ as in the continuous case. Long waves are resolved well with a sufficiently small time step. However, for the shortest resolved wave $\lambda = 2\Delta x$ it follows $\omega = \frac{1}{\Delta t}\arctan\left(\frac{U\Delta t\sin\pi}{\Delta x}\right) = 0$, so the wave is stationary. Fig. 25.1 shows the dispersion relation for the case $U = 10$ m/s. One sees that waves with $k\leq \frac{1}{\Delta x}$, i.e. $\lambda\geq 2\pi\Delta x$, are simulated decently.

25.1.1.1 First-order upwind method

In the case $U>0$, the discretization

\[ \begin{align} \frac{f_l^{(m + 1)} - f_l^{(m)}}{\Delta t} + U\frac{f_{l}^{(m)} - f_{l - 1}^{(m)}}{\Delta x} = 0 \end{align} \]

is called the first-order upwind method. This leads to

\[ \begin{align} \frac{e^{C\Delta t}e^{-i\omega\Delta t} - 1}{\Delta t} + U \frac{1 - e^{-ik\Delta x}}{\Delta x} = \frac{e^{C\Delta t}\cos\left(\omega\Delta t\right) - ie^{C\Delta t}\sin\left(\omega\Delta t\right) - 1}{\Delta t} + U\frac{1 - \cos\left( k\Delta x\right) + i\sin\left(k\Delta x\right)}{\Delta x} = 0. \end{align} \]

The imaginary part of this is

\[ \begin{align} &\frac{e^{C\Delta t}\sin\left(\omega\Delta t\right)}{\Delta t} = U\frac{\sin\left(k\Delta x\right)}{\Delta x}\nonumber\\ &\Leftrightarrow e^{C\Delta t}\sin\left(\omega\Delta t\right) = U\frac{\Delta t\sin\left(k\Delta x\right)}{\Delta x}\nonumber\\ &\Leftrightarrow e^{C\Delta t} = U\frac{\Delta t\sin\left(k\Delta x\right)}{\Delta x\sin\left(\omega\Delta t\right)}. \end{align} \]

The discretization is stable if and only if $C\leq 0$, i.e. for

\[ \begin{align} U\frac{\Delta t\sin\left(k\Delta x\right)}{\Delta x\sin\left(\omega\Delta t\right)}\leq 1. \end{align} \]

Because of

\[ \begin{align} c = \frac{\omega}{k}\leq\frac{\Delta x}{\Delta t}\Leftrightarrow\omega\Delta t\leq k\Delta x \end{align} \]

the following estimate further follows

This equation is called the Courant-Friedrichs-Lewy criterion (CFL criterion). One defines the Courant number or CFL number $c$ by

\[ \begin{align} c\coloneqq\frac{U\Delta t}{\Delta x}. \end{align} \]

The CFL criterion, expressed in words, reads:

25.1.2 Leapfrog method

If one uses a central temporal difference quotient, i.e. writes

\[ \begin{align} \frac{f_l^{(m + 2)} - f_l^{(m)}}{2\Delta t} + U\frac{f_{l + 1}^{(m + 1)} - f_{l - 1}^{(m + 1)}}{2\Delta x} = 0, \end{align} \]

then inserting this yields

\[ \begin{align} \frac{e^{2C\Delta t}e^{-2i\omega\Delta t} - 1}{2\Delta t} + U \frac{e^{ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t} - e^{-ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t}}{2\Delta x} &= 0 \nonumber\\ \Leftrightarrow \frac{e^{C\Delta t}e^{-i\omega\Delta t} - e^{-C\Delta t}e^{i\omega\Delta t}}{2\Delta t} + U \frac{e^{ik\Delta x} - e^{-ik\Delta x}}{2\Delta x} &= 0. \end{align} \]

The real part yields

\[ \begin{align} e^{C\Delta t}\cos\left(\omega\Delta t\right) = e^{-C\Delta t}\cos\left(\omega\Delta t\right). \end{align} \]

This means $C = 0$, so the solution is stable. The imaginary part becomes

\[ \begin{align} \frac{\sin\left(\omega\Delta t\right)}{2\Delta t}\left(e^{-C\Delta t} + e^{C\Delta t}\right) = \frac{\sin\left(\omega\Delta t\right)}{\Delta t} = U\frac{\sin\left(k\Delta x\right)}{\Delta x}. \end{align} \]

It follows

\[ \begin{align} \sin\left(\omega\Delta t\right) = U\frac{\Delta t}{\Delta x}\sin\left(k\Delta x\right)\Rightarrow\omega = \frac{1}{\Delta t}\arcsin\left[\frac{U\Delta t}{\Delta x}\sin\left(k\Delta x\right)\right], \end{align} \]

i.e. a different dispersion relation than in the explicit case. In general, the stability of the solution seems to follow from the real part, while the dispersion seems to follow from the imaginary part. Schemes are still either stable ($C \leq 0$), unstable ($C > 0$), or conditionally stable (stable for sufficiently small $\Delta t$). This is influenced by the time step method. The optimal case is $C = 0$, which results from the leapfrog method, so this method is the most suitable.

25.1.3 General multistep methods

25.1.4 Runge-Kutta method

Let $f: \mathbb{R}^l\times\mathbb{R} \to \mathbb{R}^m$ with $l, m \geq 1$ be an infinitely often continuously differentiable function that satisfies the initial value problem

\[ \begin{align} \frac{df}{dt} = F\left[f\left(t\right)\right], & {} & f_0 = f_0 \end{align} \]

So-called Runge-Kutta methods are based on computing several values for $df/dt$ and superimposing them in such a way that as many terms of the Taylor series as possible are covered.

25.1.4.1 Linear case

The Taylor expansion of $f$ is

\[ \begin{align} f\left(t\right) = \sum_{k = 0}^\infty\frac{f^{(k)}\left(0\right)}{k!}t^k. \end{align} \]

The function $F$ turns a function $f$ into its time derivative, so it follows

\[ \begin{align} f\left(t\right) = \sum_{k = 0}^\infty\frac{F^k\left(0\right)}{k!}t^k \end{align} \]

with $F^0 = f$. Assuming that $F$ is furthermore linear, one can reformulate this as

\[ \begin{align} f\left(t\right) &= \sum_{k = 0}^n\frac{F^k\left(0\right)}{k!}t^k + \sum_{k = n + 1}^\infty\frac{F^k\left(0\right)}{k!}t^k\nonumber\\ &= f_0 + tF\left\{f_0 + \frac{t}{2}F\left[f_0 + \dots + \frac{t}{n}F\left(f_0\right)\right]\right\} + O\left(t^{n + 1}\right)\tag{25.30}\label{eq:rk_lin} \end{align} \]

Thus the first $n$ derivatives of $f$ are approximated by this method, which is called the linear third-order Runge-Kutta method.

25.1.4.2 Nonlinear case

Defining

\[ \begin{align} f_1 \coloneqq \frac{t}{2}F\left[f_0 + \dots + \frac{t}{n}F\left(f_0\right)\right], \end{align} \]

one can write Eq. (25.30) in the form

\[ \begin{align} f\left(t\right) &= f_0 + tF\left(f_0 + f_1\right) + O\left(t^{n + 1}\right) \end{align} \]

If the function $F$, which determines the time evolution of the system, is nonlinear, the following no longer holds

\[ \begin{align} tF\left(f_0 + f_1\right) = tF\left(f_0\right) + tF\left(f_1\right), \end{align} \]

but rather

\[ \begin{align} tF\left(f_0 + f_1\right) = tF\left(f_0\right) + tF'\left(f_0\right)f_1 + tO\left(f_1^2\right). \end{align} \]

Thus, in this case, one has

\[ \begin{align} f\left(t\right) &= f_0 + tF\left(f_0 + f_1\right) + O\left(t^3\right). \end{align} \]

In the nonlinear case, a discretization of the form of Eq. (25.30) is therefore always of second order.

25.2 Methods with implicit parts

25.2.1 Implicit Euler method

If one uses a left-sided difference quotient for the time derivative (this is called the implicit Euler method), i.e. writes

\[ \begin{align} \frac{f_l^{(m + 1)} - f_l^{(m)}}{\Delta t} + U\frac{f_{l + 1}^{(m + 1)} - f_{l - 1}^{(m + 1)}}{2\Delta x} = 0, \end{align} \]

then inserting this yields

\[ \begin{align} & \frac{e^{C\Delta t}e^{-i\omega\Delta t} - 1}{\Delta t} + U \frac{e^{ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t} - e^{-ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t}}{2\Delta x} = 0 \nonumber\\ &\Leftrightarrow\frac{1 - e^{-C\Delta t}e^{i\omega\Delta t}}{\Delta t} + U \frac{e^{ik\Delta x} - e^{-ik\Delta x}}{2\Delta x} = 0. \end{align} \]

The real part yields

\[ \begin{align} 1 = e^{-C\Delta t}\cos\left(\omega\Delta t\right). \end{align} \]

This means $C<0$, so the solution is always stable, but it converges to zero. The imaginary part becomes

\[ \begin{align} \frac{e^{-C\Delta t}\sin\left(\omega\Delta t\right)}{\Delta t} = U\frac{\sin\left(k\Delta x\right)}{\Delta x}. \end{align} \]

It follows

\[ \begin{align} \tan\left(\omega\Delta t\right) = U\frac{\Delta t}{\Delta x}\sin\left(k\Delta x\right), \end{align} \]

i.e. the same dispersion relation as in the explicit Euler method.

25.2.2 Crank-Nicolson method

The Crank-Nicolson method is a combination of explicit and implicit Euler methods:

\[ \begin{align} \frac{f_l^{(m + 1)} - f_l^{(m)}}{\Delta t} + U\frac{f_{l + 1}^{(m)} - f_{l - 1}^{(m)} + f_{l + 1}^{(m + 1)} - f_{l - 1}^{(m + 1)}}{4\Delta x} = 0 \end{align} \]

Inserting Eq. (25.8) here, one obtains

\[ \begin{align} & \frac{e^{C\Delta t}e^{-i\omega\Delta t} - 1}{\Delta t} + U \frac{e^{ik\Delta x} - e^{-ik\Delta x} + e^{ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t} - e^{-ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t}}{4\Delta x} = 0\nonumber\\ &\Leftrightarrow\frac{e^{C\Delta t}\cos\left(\omega\Delta t\right) - ie^{C\Delta t}\sin\left(\omega\Delta t\right) - 1}{\Delta t} + U \frac{2i\sin\left(\omega\Delta t\right) + 2i\sin\left(k\Delta x\right)e^{C\Delta t}e^{-i\omega\Delta t}}{4\Delta x} = 0\nonumber\\ &\Leftrightarrow\frac{e^{C\Delta t}\cos\left(\omega\Delta t\right) - ie^{C\Delta t}\sin\left(\omega\Delta t\right) - 1}{\Delta t} + U \frac{2i\sin\left(\omega\Delta t\right) + 2ie^{C\Delta t}\sin\left(k\Delta x\right)\cos\left(\omega\Delta t\right) + 2e^{C\Delta t}\sin\left(k\Delta x\right)\sin\left(\omega\Delta t\right)}{4\Delta x} = 0. \end{align} \]

The real part yields

\[ \begin{align} &\Leftrightarrow\frac{e^{C\Delta t}\cos\left(\omega\Delta t\right) - 1}{\Delta t} + U\frac{2e^{C\Delta t}\sin\left(k\Delta x\right)\sin\left(\omega\Delta t\right)}{4\Delta x} = 0\nonumber\\ &\Leftrightarrow\frac{\cos\left(\omega\Delta t\right) - e^{-C\Delta t}}{\Delta t} + U\frac{\sin\left(k\Delta x\right)\sin\left(\omega\Delta t\right)}{2\Delta x} = 0\nonumber. \end{align} \]

The imaginary part yields

\[ \begin{align} &\Leftrightarrow-\frac{ie^{C\Delta t}\sin\left(\omega\Delta t\right)}{\Delta t} + U\frac{2i\sin\left(\omega\Delta t\right) + 2ie^{C\Delta t}\sin\left(k\Delta x\right)\cos\left(\omega\Delta t\right)}{4\Delta x} = 0\nonumber\\ &\Leftrightarrow-\frac{e^{C\Delta t}}{\Delta t} + U\frac{1 + e^{C\Delta t}\sin\left(k\Delta x\right)\cot\left(\omega\Delta t\right)}{2\Delta x} = 0\nonumber. \end{align} \]

This means $C < 0$, so the solution is stable, but it vanishes eventually. The imaginary part becomes

\[ \begin{align} \frac{e^{-C\Delta t}\sin\left(\omega\Delta t\right)}{\Delta t} = U\frac{\sin\left(k\Delta x\right)}{\Delta x}. \end{align} \]

It follows

\[ \begin{align} \tan\left(\omega\Delta t\right) = U\frac{\Delta t}{\Delta x}\sin\left(k\Delta x\right), \end{align} \]

i.e. the same dispersion relation as in the explicit Euler method.

25.3 Mixture of time step methods

One is not restricted to solving a differential equation either purely explicitly or purely implicitly. Let $\newtilde{f}$ be the discretized analogue of a continuous field $f$ and let $\mathbf{r}_i$ be the grid points. A two-step method for a differential equation of first order in time for $f$ can be written in the form

\[ \begin{align} \newtilde{f}\left(\mathbf{r}_i, \left(n + 1\right)\Delta t\right) = \newtilde{f}\left(\mathbf{r}_i, n\Delta t\right) + \Delta t\sum_{j = 1}^N\newtilde{F}_j\left[\newtilde{f}\left(\left(n + m_j\right)\Delta t\right)\right] \end{align} \]

Here $\newtilde{F}_j$ for $j = 1, \dots, N$ are the discretized forcings acting on $\newtilde{f}$. Here, the following should hold:

\[ \begin{align} m_j = \begin{cases} 0,\text{ term $j$ is treated explicitly,}\\ 1,\text{ term $j$ is treated implicitly.} \end{cases} \end{align} \]

It was tacitly assumed that $f$ is the only field appearing. More generally, $f$ can be regarded as a tuple of fields, i.e. as a state vector. Furthermore, one can allow $0 \leq m_j \leq 1$; here $m_j$ is the implicit fraction. The Crank-Nicolson method corresponds to $m_j = \frac{1}{2}$.

25.4 Aspects of partial differential equations

25.4.1 Forward-backward method

25.5 Energetic self-consistency in the ideal gas

25.5.1 Advection of momentum

25.5.1.1 Advection of kinetic energy

Since $\mathbf{v}\cdot\left[\left(\nabla\times\mathbf{v}\right)\times\mathbf{v}\right] = 0$, the generalized Coriolis term is energetically irrelevant. For a self-consistent temporal discretization of the advection of kinetic energy, one therefore uses the system of equations

\[ \begin{align} \frac{\partial\mathbf{v}}{\partial t} = -\nabla\frac{\mathbf{v}^2}{2}, & {} & \frac{\partial\rho}{\partial t} = -\nabla\cdot\left(\rho\mathbf{v}\right). \end{align} \]

In the continuous case, the following holds:

\[ \begin{align} \rho\mathbf{v}\cdot\frac{\partial\mathbf{v}}{\partial t} &= -\rho\mathbf{v}\cdot\nabla\frac{\mathbf{v}^2}{2} = -\nabla\cdot\left(\rho\mathbf{v}\frac{\mathbf{v}^2}{2}\right) + \frac{\mathbf{v}^2}{2}\nabla\cdot\left(\rho\mathbf{v}\right) = -\nabla\cdot\left(\rho\mathbf{v}\frac{\mathbf{v}^2}{2}\right) - \frac{\mathbf{v}^2}{2}\frac{\partial\rho}{\partial t}\nonumber\\ \Rightarrow \rho\mathbf{v}\cdot\frac{\partial\mathbf{v}}{\partial t} + \frac{\mathbf{v}^2}{2}\frac{\partial\rho}{\partial t} &= \frac{\partial\left(\rho\mathbf{v}^2/2\right)}{\partial t} = -\nabla\cdot\left(\rho\mathbf{v}\frac{\mathbf{v}^2}{2}\right) \end{align} \]

for the effect of the kinetic energy flux density on the local kinetic energy density. These transformation steps must be transferable to the discretization. One now defines $\delta_0, \delta_1, \delta_2, \delta_3, \delta_4, \delta_5 \in \left\{0, 1\right\}$ in order to express the temporal discretization of the kinetic energy and the mass flux density in the form

\[ \begin{align} \frac{\mathbf{v}^2}{2} = \frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}, & {} & \rho\mathbf{v} = \frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2} \end{align} \]

With this, one first obtains the system of equations

\[ \begin{align} \frac{\mathbf{v}_1 - \mathbf{v}_0}{\Delta t} = -\nabla\frac{\mathbf{v}^2}{2}, & {} & \frac{\rho_1 - \rho_0}{\Delta t} = -\nabla\cdot\left(\rho\mathbf{v}\right). \end{align} \]

From this it follows

\[ \begin{align} &\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\cdot\frac{\mathbf{v}_1 - \mathbf{v}_0}{\Delta t} = -\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\cdot\nabla\frac{\mathbf{v}^2}{2} = -\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\cdot\nabla \frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\nonumber\\ &= -\nabla\cdot\left(\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\right) + \frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\nabla\cdot\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\nonumber\\ &= -\nabla\cdot\left(\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\right) - \frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\frac{\rho_1 - \rho_0}{\Delta t}\nonumber\\ &\Leftrightarrow\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\cdot\frac{\mathbf{v}_1 - \mathbf{v}_0}{\Delta t} + \frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\frac{\rho_1 - \rho_0}{\Delta t} = -\nabla\cdot\left(\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\right)\nonumber\\ &\Leftrightarrow \frac{\rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 + \rho_1\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1} - \rho_0\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2\Delta t} = -\nabla\cdot\left(\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\right)\nonumber\\ &\Leftrightarrow \frac{\rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1} - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2\Delta t} = -\nabla\cdot\left(\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\right) \end{align} \]

The following must hold:

\[ \begin{align} \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1} - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

Inserting $\delta_0 = \delta_1 = 0$, one obtains

\[ \begin{align} \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{0}\cdot\mathbf{v}_{0} - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

With $\delta_2 = \delta_3 = 1$ it follows

\[ \begin{align} \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{0}\cdot\mathbf{v}_{0} - \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

This in turn results in $\delta_4 = 1$, $\delta_5 = 0$. The pairs $\left(\delta_2, \delta_3\right)$, $\left(\delta_4, \delta_5\right)$ are interchangeable, but this is not a new possibility since the mass flux density does not change.

Inserting $\delta_0 = \delta_1 = 1$, one obtains

\[ \begin{align} \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{1}\cdot\mathbf{v}_{1} - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{1}\cdot\mathbf{v}_{1} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

With $\delta_2 = 0$, $\delta_3 = 1$, it follows

\[ \begin{align} \rho_{0}\mathbf{v}_{1}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{1}\cdot\mathbf{v}_{1} - \rho_{0}\mathbf{v}_{1}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{1}\cdot\mathbf{v}_{1} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

This in turn results in $\delta_4 = 0$, $\delta_5 = 0$.

Inserting $\delta_0 = 0$, $\delta_1 = 1$, one obtains

\[ \begin{align} \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{0}\cdot\mathbf{v}_{1} - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{1} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

With $\delta_2 = \delta_3 = 0$ it follows

\[ \begin{align} \rho_{0}\mathbf{v}_{0}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{0}\cdot\mathbf{v}_{1} - \rho_{0}\mathbf{v}_{0}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{1} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

This in turn results in $\delta_4 = 1$, $\delta_5 = 1$.

In summary, there are the following three options:

\[ \begin{align} \frac{\mathbf{v}^2}{2} = \frac{\mathbf{v}_0\cdot\mathbf{v}_0}{2}, & \rho\mathbf{v} = \frac{\rho_1\mathbf{v}_0 + \rho_1\mathbf{v}_1}{2}\\ \frac{\mathbf{v}^2}{2} = \frac{\mathbf{v}_1\cdot\mathbf{v}_1}{2}, & \rho\mathbf{v} = \frac{\rho_0\mathbf{v}_0 + \rho_0\mathbf{v}_1}{2}\\ \frac{\mathbf{v}^2}{2} = \frac{\mathbf{v}_0\cdot\mathbf{v}_1}{2}, & \rho\mathbf{v} = \frac{\rho_0\mathbf{v}_0 + \rho_1\mathbf{v}_1}{2} \end{align} \]

25.5.1.2 Generalized Coriolis terms

For the time evolution of the kinetic energy, the following holds:

\[ \begin{align} \frac{\mathbf{v}_1^2 - \mathbf{v}_0^2}{2\Delta t} &= \frac{\left(\mathbf{v}_1 + \mathbf{v}_0\right)\cdot\left(\mathbf{v}_1 - \mathbf{v}_0\right)}{2\Delta t} = \mathbf{v}_{1/2}\cdot\frac{\mathbf{v}_1 - \mathbf{v}_0}{\Delta t} \end{align} \]

with

\[ \begin{align} \mathbf{v}_{1/2} \coloneqq \frac{\mathbf{v}_0 + \mathbf{v}_1}{2}. \end{align} \]

If the generalized Coriolis term is not to contribute to the evolution of the kinetic energy, it must therefore be computed using the velocity field $\mathbf{v}_{1/2}$. The considerations in this section (25.5.1) were first put forward by Gassmann and Herzog in [23].

25.5.1.3 Time evolution of internal energy

For the time evolution of the internal energy density $\newtilde{I}$, one can write

\[ \begin{align} \frac{\Delta\newtilde{I}}{\Delta t} &= \frac{\newtilde{I}_1 - \newtilde{I}_0}{\Delta t} = c^{(v)}\rho_0\frac{T_1 - T_0}{\Delta t} + c^{(v)}T_1\frac{\rho_1 - \rho_0}{\Delta t}\tag{25.64}\label{eq:deriv_temp_evol_model_1}\\ &= \frac{\newtilde{I}_1 - \newtilde{I}_0}{\Delta t} = c^{(v)}\rho_1\frac{T_1 - T_0}{\Delta t} + c^{(v)}T_0\frac{\rho_1 - \rho_0}{\Delta t}.\tag{25.65}\label{eq:deriv_temp_evol_model_2} \end{align} \]

First, one starts from Eq. (25.64). Using the continuity equation one obtains

\[ \begin{align} \md{\rho} = -\rho_0\nabla\cdot\mathbf{v}. \end{align} \]

From Eq. (9.14) it follows that

\[ \begin{align} \frac{T_1 - T_0}{\Delta t} = -\frac{R_d}{c^{(v)}}T_0\nabla\cdot\mathbf{v}. \end{align} \]

Both times it was assumed that the scalar variables from the old time step are used, as is the case, for example, with a forward-backward method. Inserting these two findings into Eq. (25.64), one obtains

\[ \begin{align} \frac{\newtilde{I}_1 - \newtilde{I}_0}{\Delta t} &= -c^{(v)}\rho_0\frac{R_d}{c^{(v)}}T_0\nabla\cdot\mathbf{v} - c^{(v)}T_1\rho_0\nabla\cdot\mathbf{v} = -c^{(v)}\rho_0\nabla\cdot\mathbf{v}\left(\frac{R_d}{c^{(v)}}T_0 + T_1\right)\nonumber\\ &= -c^{(p)}\rho_0\nabla\cdot\mathbf{v}\left(\frac{R_d}{c^{(p)}}T_0 + \frac{c^{(v)}}{c^{(p)}}T_1\right). \end{align} \]

From this one concludes that the temperature, in the weighting

\[ \begin{align} T^\star \coloneqq \frac{R_d}{c^{(p)}}T_0 + \frac{c^{(v)}}{c^{(p)}}T_1\tag{25.69}\label{eq:c_v_c_p} \end{align} \]

enters into the discretized evolution of the internal energy density.

However, if one starts from Eq. (25.65), this leads, with

\[ \begin{align} \md{\rho} = -\rho_1\nabla\cdot\mathbf{v} \end{align} \]

and

\[ \begin{align} \frac{T_1 - T_0}{\Delta t} = -\frac{R_d}{c^{(v)}}T_1\nabla\cdot\mathbf{v} \end{align} \]

to

\[ \begin{align} \frac{\newtilde{I}_1 - \newtilde{I}_0}{\Delta t} &= -c^{(p)}\rho_1\nabla\cdot\mathbf{v}\left(\frac{c^{(v)}}{c^{(p)}}T_0 + \frac{R_d}{c^{(p)}}T_1\right),\\ T^\star &\coloneqq \frac{c^{(v)}}{c^{(p)}}T_0 + \frac{R_d}{c^{(p)}}T_1.\tag{25.73}\label{eq:r_d_c_p} \end{align} \]

These are the two possibilities for energetically self-consistent temporal discretizations of the internal energy.