The aim of this chapter is to derive how differential and algebraic operators and boundary conditions are to be discretized. To do this, consider the atmosphere $A \subseteq \mathbb {R}^3$ as an open set and use the reversible equations of dry air with the boundary conditions as the system of equations.
\[ \begin{align} \mathbf{v}\cdot\mathbf{n} = 0 \end{align} \]
at the top and
\[ \begin{align} \mathbf{v} = \mathbf{0} \end{align} \]
at the bottom.
As prognostic variables
chosen. The Hamilton function $H$ is then written as a function of $\rho, \newtilde{s}$ and $\mathbf{v}$,
\[ \begin{align} H = H\left(\rho, \newtilde{\theta}, \mathbf{v}\right). \end{align} \]
The prognostic equations for this case were derived in Section 10.1.2, they are
\[ \begin{align} \frac{\partial\mathbf{v}}{\partial t} &= -c^{(p)}\theta\nabla\Pi - \nabla k + \mathbf{v}\times\etabi - \nabla\phi,\tag{34.4}\label{eq:prog_rev_0}\\ \frac{\partial\rho}{\partial t} &= -\nabla\cdot\left(\rho\mathbf{v}\right),\tag{34.5}\label{eq:prog_rev_1}\\ \frac{\partial\newtilde{\theta}}{\partial t} &= -\nabla\cdot\left(\newtilde{\theta}\mathbf{v}\right),\tag{34.6}\label{eq:prog_rev_2}\\ T &= f\left(\rho, \newtilde{\theta}\right).\tag{34.7}\label{eq:prog_rev_3} \end{align} \]
The function $f$ in Eq. (34.7) is the thermal equation of state of ideal gases, $\etabi\coloneqq\zetabi + \mathbf{f}$ is the absolute vorticity.
The Poisson brackets are global integrals. One therefore defines a discrete integral operator according to the substitution
\[ \begin{align} \int_A\psi\left(\mathbf{r}\right)d^3r \to \sum_{i = 1}^N\psi_iV_i, \end{align} \]
here $\psi\left(\mathbf{r}\right)$ is an integrable scalar or vector field, $\psi_i$ is the value of the discretized version of $\psi$ in grid box $i$, $V_i$ is the volume of grid box $i$ and $N$ is the number of grid boxes. The Poisson brackets are calculated according to the basic approach
\[ \begin{align} \frac{dF\left[\mathbf{u}\right]}{dt} = \frac{d}{dt}\int_Af\left(\mathbf{u}\left(t\right)\right)d^3r = \int_A\frac{\delta F}{\delta\mathbf{u}}\cdot\frac{\partial\mathbf{u}}{\partial t}d^3r \to \sum_{c, k}\left(\frac{\delta F}{\delta\mathbf{u}}\cdot\frac{\partial\mathbf{u}}{\partial t}\right)_{c, k}V_{c, k} \end{align} \]
discretized. First the Hamilton functional is needed. One discretizes Eq. (10.81) in the form
\[ \begin{align} H\left(\rho, \newtilde{s}, \mathbf{v}\right) = \sum_{c, k}\left[\frac{1}{2}\rho_{c, k}\left(\mathbf{v}\cdot\mathbf{v}\right)_{c, k} + \rho_{c, k}\phi_{c, k} + \newtilde{I}_{c, k}\right]V_{c, k}. \end{align} \]
This is divided into two parts
\[ \begin{align} H_\text{kin}\left(\mathbf{v}\right) &= \sum_{c, k}\frac{1}{2}\rho_{c, k}\left(\mathbf{v}\cdot\mathbf{v}\right)_{c, k}V_{c, k},\\ H_{I'}\left(\rho, \newtilde{s}\right) &= \sum_{c, k}\left(\rho_{c, k}\phi_{c, k} + \newtilde{I}_{c, k}\right)V_{c, k}. \end{align} \]
The proportion of kinetic energy can be divided even further:
\[ \begin{align} H_\text{kin}^{(u)}\left(u\right) &= \sum_{c, k}\frac{1}{2}\rho_{c, k}\sum_{e\in c}\frac{l_ed_e}{2A_c\Delta z_{c, k}}\Delta z_{e, k}u_{e, k}^2V_{c, k} = \sum_{c, k}\sum_{e\in c}\frac{\rho_{c, k}l_ed_e\Delta z_{e, k}}{4}u_{e, k}^2,\tag{34.13}\label{eq:h_kin_u}\\ H_\text{kin}^{(w)}\left(w\right) &= \sum_{c, k}\frac{1}{2}\rho_{c, k}\left(\frac{\Delta z_{c, k - 1/2}}{2\Delta z_{c, k}}w_u^2 + \frac{\Delta z_{c, k + 1/2}}{2\Delta z_{c, k}}w_l^2\right)V_{c, k}\nonumber\\ &= \sum_{c, k}\left(\frac{\rho_{c, k}\Delta z_{c, k - 1/2}}{4}w_{c, k + 1/2}^2 + \frac{\rho_{c, k}\Delta z_{c, k + 1/2}}{4}w_{c, k + 1/2}^2\right)A_{c}\tag{34.14}\label{eq:h_kin_w} \end{align} \]
An alternative formulation of kinetic energy is
\[ \begin{align} H_\text{kin}^\star\left(\mathbf{v}\right) &\coloneqq \sum_{c, k}\frac{1}{2}\left(\rho\mathbf{v}\cdot\mathbf{v}\right)_{c, k}V_{c, k},\tag{34.15}\label{eq:h_kin_mod}\\ H_\text{kin}^{(u, \star)}\left(u\right) &\coloneqq \sum_{c, k}\sum_{e\in c}\frac{l_ed_e\Delta z_{e, k}}{4}\newtilde{\rho_{c, k}}^{(e)}u_{e, k}^2,\tag{34.16}\label{eq:h_kin_u_mod}\\ H_\text{kin}^{(w, \star)}\left(w\right) &\coloneqq \sum_{c, k}\left(\frac{\Delta z_{c, k - 1/2}}{4}\newtilde{\rho_{c, k}}^{(k - 1/2)}w_{c, k - 1/2}^2 + \frac{\Delta z_{c, k + 1/2}}{4}\newtilde{\rho_{c, k}}^{(k + 1/2)}w_{c, k + 1/2}^2\right)A_{c}.\tag{34.17}\label{eq:h_kin_w_mod} \end{align} \]
Here $\newtilde{\rho _{c, k}}^{(e)}$ is an averaging operator on the edges and $\newtilde{\rho _{c, k}}^{(k + 1/2)}$ is an averaging operator on the points where $w$ is placed. These operators must
\[ \begin{align} \sum_{e, k}\newtilde{\rho_{c, k}}^{(e)}V_{e, k} &= \sum_{c, k}\rho_{c, k}V_{c, k},\\ \sum_{c, k}\newtilde{\rho_{c, k}}^{(k + 1/2)}V_{c, k + 1/2} &= \sum_{c, k}\rho_{c, k}V_{c, k} \end{align} \]
in order not to endanger mass conservation.
Eq. (A.122) is discretized in the form
\[ \begin{align} \lim_{\epsilon\to 0}\frac{1}{\epsilon}\sum_{c, k}\left[f\left(\mathbf{u}_{c, k} + \epsilon\mathbf{h}_{c, k}\right) - f\left(\mathbf{u}_{c, k}\right)\right]V_{c, k} \hastobe \sum_{c, k}\left(\frac{\delta F}{\delta\mathbf{u}}\right)_{c, k}\cdot\mathbf{h}_{c, k}V_{c, k}. \end{align} \]
The functional derivative is the discrete function that fulfills this for all auxiliary functions $\mathbf{h}_{c, k}$. $\mathbf{h}$ is a tuple of state variables. It can also be summed over volumes other than $V_{c, k}$.
For the functional derivatives according to the thermodynamic quantities one obtains
\[ \begin{align} \left(\frac{\delta H_\text{kin}}{\delta\rho}\right)_{c, k} = \frac{1}{2}\left(\mathbf{v}\cdot\mathbf{v}\right)_{c, k}, & {} & \left(\frac{\delta H_{I'}}{\delta\rho}\right)_{c, k} = \phi_{c, k} + \frac{\partial\newtilde{I}_{c, k}}{\partial\rho}, & {} & \left(\frac{\delta H_{I'}}{\delta q_1}\right)_{c, k} = \frac{\partial\newtilde{I}_{c, k}}{\partial q_1}. \end{align} \]
The equations (34.16) - (34.17) can be written in the form
\[ \begin{align} H_\text{kin}^{(u, \star)}\left(u\right) &\coloneqq \sum_{c, k}\sum_{e\in c}V_{e, k}\newtilde{\rho_{c, k}}^{(e)}\frac{u_{e, k}^2}{2},\\ H_\text{kin}^{(w, \star)}\left(w\right) &\coloneqq \sum_{c, k}\left(\frac{\Delta z_{c, k - 1/2}}{4}\newtilde{\rho_{c, k}}^{(k - 1/2)}w_{c, k - 1/2}^2 + \frac{\Delta z_{c, k + 1/2}}{4}\newtilde{\rho_{c, k}}^{(k + 1/2)}w_{c, k + 1/2}^2\right)A_{c}\nonumber\\ &= \sum_{c, k}\left(\frac{\Delta z_{c, k - 1/2}}{4}\newtilde{\rho_{c, k}}^{(k - 1/2)}w_{c, k - 1/2}^2 + \frac{\Delta z_{c, k + 1/2}}{4}\newtilde{\rho_{c, k}}^{(k + 1/2)}w_{c, k + 1/2}^2\right)A_{c}\nonumber\\ &= \sum_{c, k}\newtilde{\rho_{c, k}}^{(k + 1/2)}\frac{w_{c, k + 1/2}^2}{2}V_{c, k + 1/2} \end{align} \]
rewrite. There were
\[ \begin{align} V_{e, k} = \frac{l_ed_e\Delta z_{e, k}}{2} \end{align} \] and
\[ \begin{align} V_{c, k + 1/2} = A_c\Delta z_{c, k + 1/2} \end{align} \]
used. This implies
\[ \begin{align} \left(\frac{\delta H_\text{kin}^\star}{\delta\rho}\right)_{c, k} = \frac{1}{2}\left(\mathbf{v}\cdot\mathbf{v}\right)_{c, k}, \end{align} \]
since the functional derivative can also be evaluated at the boundaries of the boxes and then interpolated back to the center points using a scalar product.
From Eq. (34.13) follows
\[ \begin{align} \frac{H_\text{kin}^{(u)}\left(u + \epsilon h\right) - H_\text{kin}^{(u)}\left(u\right)}{\epsilon} &= \sum_{c, k}\sum_{e\in c}\frac{\rho_{c, k}l_ed_e\Delta z_{e, k}}{2}u_{e, k}h_{e, k} + O\left(h^2\right) = \sum_{c, k}\sum_{e\in c}\rho_{c, k}u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right).\tag{34.27}\label{eq:shallow_func_deriv_deriv_0} \end{align} \]
The sum over cells can be rewritten as a sum over edges
\[ \begin{align} \frac{H_\text{kin}^{(u)}\left(u + \epsilon h\right) - H_\text{kin}^{(u)}\left(u\right)}{\epsilon} &= \sum_{e, k}\left(\rho_{c_1(e), k} + \rho_{c_2(e), k}\right)u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right)\nonumber\\ &= 2\sum_{e, k}\frac{\rho_{c_1(e), k} + \rho_{c_2(e), k}}{2}u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right) \end{align} \]
The prefactor two stands for the two spatial directions under consideration. The two cells bordering the edge $e$ were labeled $c_0$ and $c_1$. Therefore applies
\[ \begin{align} \left(\frac{\delta H}{\delta\mathbf{v}}\right)_{e, k} = \frac{\rho_{c_1(e), k} + \rho_{c_2(e), k}}{2}u_{e, k}. \end{align} \]
Analogously it follows from Eq. (34.14)
\[ \begin{align} \frac{H_\text{kin}^{(w)}\left(w + \epsilon h\right) - H_\text{kin}^{(w)}\left(w\right)}{\epsilon} &= \sum_{c, k}\left(\frac{\rho_{c, k}\Delta z_{c, k - 1/2}}{2}w_{c, k - 1/2}h_{c, k - 1/2} + \frac{\rho_{c, k}\Delta z_{c, k + 1/2}}{2}w_{c, k + 1/2}h_{c, k + 1/2}\right)A_{c}\nonumber\\ &= \sum_{c, k}\left(\frac{\rho_{c, k}}{2}w_{c, k + 1/2}h_{c, k + 1/2} + \frac{\rho_{c, k + 1}}{2}w_{c, k + 1/2}h_{c, k + 1/2}\right)\Delta z_{c, k + 1/2}A_{c}.\tag{34.30}\label{eq:shallow_func_deriv_deriv_1} \end{align} \]
An index shift was made. With $V_{c, k + 1/2} = A_c\Delta z_{c, k + 1/2}$ you get
\[ \begin{align} \frac{H_\text{kin}^{(w)}\left(w + \epsilon h\right) - H_\text{kin}^{(w)}\left(w\right)}{\epsilon} &= \sum_{c, k}\frac{\rho_{c, k} + \rho_{c, k + 1}}{2}w_{c, k + 1/2}h_{c, k + 1/2}V_{c, k + 1/2}. \end{align} \]
Thus follows
\[ \begin{align} \left(\frac{\delta H}{\delta\mathbf{v}}\right)_{c, k + 1/2} = \frac{\rho_{c, k} + \rho_{c, k + 1}}{2}w_{c, k + 1/2}. \end{align} \]
Using the modified operator $H_\text{kin}^\star\left(\mathbf{v}\right)$ from Eq. (34.15), you get
\[ \begin{align} \left(\frac{\delta H^\star}{\delta\mathbf{v}}\right)_{e, k} &= \newtilde{\rho_{c, k}}^{(e)}u_{e, k},\tag{34.33}\label{eq:h_mod_func_deriv_0}\\ \left(\frac{\delta H^\star}{\delta\mathbf{v}}\right)_{c, k + 1/2} &= \newtilde{\rho_{c, k}}^{(k + 1/2)}w_{c, k + 1/2}.\tag{34.34}\label{eq:h_mod_func_deriv_1} \end{align} \]
Eq. (34.27) generalizes with
\[ \begin{align} V_{e, k} \coloneqq A_{e,k}\frac{d_{e, k}}{2} \end{align} \]
to the deep atmosphere as follows:
\[ \begin{align} \frac{H_\text{kin}^{(u)}\left(u + \epsilon h\right) - H_\text{kin}^{(u)}\left(u\right)}{\epsilon} &= \sum_{c, k}\sum_{e\in c}\frac{\rho_{c, k}d_eA_{e, k}}{2}u_{e, k}h_{e, k} + O\left(h^2\right) = \sum_{c, k}\sum_{e\in c}\rho_{c, k}u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right) \end{align} \]
Now write down the sum over cells again as a sum over edges:
\[ \begin{align} \frac{H_\text{kin}^{(u)}\left(u + \epsilon h\right) - H_\text{kin}^{(u)}\left(u\right)}{\epsilon} &= \sum_{e, k}\left(\rho_{c_1(e), k} + \rho_{c_2(e), k}\right)u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right)\nonumber\\ &= 2\sum_{e, k}\frac{\rho_{c_1(e), k} + \rho_{c_2(e), k}}{2}u_{e, k}h_{e, k}V_{e, k} + O\left(h^2\right) \end{align} \]
This also applies in the deep atmosphere
\[ \begin{align} \left(\frac{\delta H}{\delta\mathbf{v}}\right)_{e, k} = \frac{\rho_{c_1(e), k} + \rho_{c_2(e), k}}{2}u_{e, k}. \end{align} \]
The procedure is analogous to Eq. (34.27):
\[ \begin{align} \frac{H_\text{kin}^{(w)}\left(w + \epsilon h\right) - H_\text{kin}^{(w)}\left(w\right)}{\epsilon} &= \sum_{c, k}\left(\frac{\rho_{c, k}A_{c, k - 1/2}\Delta z_{c, k - 1/2}}{2V_{c, k}}w_{c, k - 1/2}h_{c, k - 1/2} + \frac{\rho_{c, k}A_{c, k + 1/2}\Delta z_{c, k + 1/2}}{2V_{c, k}}w_{c, k + 1/2}h_{c, k + 1/2}\right)V_{c, k}\nonumber\\ &= \sum_{c, k}\frac{\rho_{c, k}A_{c, k - 1/2}\Delta z_{c, k - 1/2}}{2}w_{c, k - 1/2}h_{c, k - 1/2} + \frac{\rho_{c, k}A_{c, k + 1/2}\Delta z_{c, k + 1/2}}{2}w_{c, k + 1/2}h_{c, k + 1/2}\nonumber\\ &= \sum_{c, k}\frac{\rho_{c, k + 1}A_{c, k + 1/2}\Delta z_{c, k + 1/2}}{2}w_{c, k + 1/2}h_{c, k + 1/2} + \frac{\rho_{c, k}A_{c, k + 1/2}\Delta z_{c, k + 1/2}}{2}w_{c, k + 1/2}h_{c, k + 1/2}\nonumber\\ &= \sum_{c, k}\left(\frac{\rho_{c, k + 1}}{2}w_{c, k + 1/2}h_{c, k + 1/2} + \frac{\rho_{c, k}}{2}w_{c, k + 1/2}h_{c, k + 1/2}\right)\Delta z_{c, k + 1/2}A_{c, k + 1/2} \end{align} \]
With $V_{c, k + 1/2} = A_c\Delta z_{c, k + 1/2}$ you get
\[ \begin{align} \frac{H_\text{kin}^{(w)}\left(w + \epsilon h\right) - H_\text{kin}^{(w)}\left(w\right)}{\epsilon} &= \sum_{c, k}\frac{\rho_{c, k} + \rho_{c, k + 1}}{2}w_{c, k + 1/2}h_{c, k + 1/2}V_{c, k + 1/2}. \end{align} \]
Thus follows
\[ \begin{align} \left(\frac{\delta H}{\delta\mathbf{v}}\right)_{c, k + 1/2} = \frac{\rho_{c, k} + \rho_{c, k + 1}}{2}w_{c, k + 1/2}. \end{align} \]
The equations (34.33) - (34.34) also apply in the deep atmosphere.
In Section 26.7 the requirement $\mathbf{v}\cdot\nabla\psi \hastobe \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v}$ was assumed in order to derive a discretization of the scalar product, the discretizations of „ Scalar field x vector field“, gradient and divergence assumed. In section 33.4.2 the scalar product was generalized to terrain-following coordinates. From this, the vertical contravariant measure could be reconstructed in Section 33.4.3. This was used in Sect. 33.4.4 to derive the shape of the divergence in terrain-following coordinates. Now the gradient must be generalized to coordinates following the terrain. For this you use
\[ \begin{align} \mathbf{v}\cdot\nabla\psi \equiv \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v} \end{align} \]
as the equation defining the gradient.
From Eq. (33.69) follow
\[ \begin{align} \psi\nabla\cdot\mathbf{v} &= \left[\sum_{e\in c}\left(l_e\frac{\Delta z_{e, k}\psi_{c, k}u_e}{A_c\Delta z_{c, k}}\right) + \frac{\psi_{c, k}w_{c, k - 1/2} - \psi_{c, k}w_{c, k + 1/2}}{\Delta z_{c, k}}\right]\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4A_c\Delta z_{c, k}}\left(\frac{\Delta z_{e, k - 1}}{\Delta z_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1} - \frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\right)\tag{34.43}\label{eq:grad_c-grid_terrain_deriv_0} \end{align} \]
and
\[ \begin{align} \nabla\cdot\left(\psi\mathbf{v}\right) &= \left[\sum_{e\in c}\left(l_e\frac{\Delta z_{e, k}\left(\psi_{o(e), k} + \psi_{c, k}\right)u_e}{2A_c\Delta z_{c, k}}\right) + \frac{\left(\psi_{c, k - 1} + \psi_{c, k}\right)w_{c, k - 1/2} - \left(\psi_{c, k + 1} + \psi_{c, k}\right)w_{c, k + 1/2}}{2\Delta z_{c, k}}\right]\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4A_c\Delta z_{c, k}}\left(\frac{\Delta z_{e, k - 1}}{\Delta z_{c, k - 1}}J_{e, k - 1}\frac{\psi_{c, k - 1} + \psi_{o(e), k - 1}}{2}u_{e, k - 1} - \frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}\frac{\psi_{c, k + 1} + \psi_{o(e), k + 1}}{2}u_{e, k + 1}\right).\tag{34.44}\label{eq:grad_c-grid_terrain_deriv_1} \end{align} \]
If you take Eq. (34.43) from Eq. (34.44), you get
\[ \begin{align} &\mathbf{v}\cdot\nabla\psi = \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v}\nonumber\\ &= \left[\sum_{e\in c}\left(l_e\frac{\Delta z_{e, k}\left(\psi_{o(e), k} - \psi_{c, k}\right)u_e}{2A_c\Delta z_{c, k}}\right) + \frac{\left(\psi_{c, k - 1} - \psi_{c, k}\right)w_{c, k - 1/2} - \left(\psi_{c, k + 1} - \psi_{c, k}\right)w_{c, k + 1/2}}{2\Delta z_{c, k}}\right]\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4A_c\Delta z_{c, k}}\left(\frac{\Delta z_{e, k - 1}}{\Delta z_{c, k - 1}}J_{e, k - 1}\frac{\psi_{c, k - 1} + \psi_{o(e), k - 1} - 2\psi_{c, k}}{2}\textcolor{red}{u_{e, k - 1}} - \frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}\frac{\psi_{c, k + 1} + \psi_{o(e), k + 1} - 2\psi_{c, k}}{2}\textcolor{red}{u_{e, k + 1}}\right)\nonumber\\ &\stackrel{\text{Glg. }\href{ch-32-kinematics.html#eq:inner_c-grid_3d_shallow_terrain}{(33.47)}}{\hastobe} \sum_{e\in c}\frac{l_ed_e}{2A_c\Delta z_{c, k}}\Delta z_{e, k}\left(\nabla\psi\right)_eu_e + \frac{\Delta z_{c, k - 1/2}}{2\Delta z_{c, k}}\left(\nabla\psi\right)_{c, k - 1/2}w_{c, k - 1/2} - \frac{\Delta z_{c, k + 1/2}}{2\Delta z_{c, k}}\left(\nabla\psi\right)_{c, k + 1/2}w_{c, k + 1/2}. \end{align} \]
The terms marked red only appear on the left side, but not on the right. $\nabla\psi$ can therefore no longer be defined in coordinates following the terrain in such a way that this equation is fulfilled. One could eliminate these terms by calculating the contravariant flux densities only via a one-sided vertical reconstruction, which would be a contradiction to the requirement of mass conservation, since in this case the flux through a surface would depend on which side the surface is viewed from. So the problem cannot be solved exactly.
Nevertheless, conservation of total energy can still be achieved. For this it is necessary that the identity
\[ \begin{align} \int_A\psi\nabla\cdot\mathbf{v}d^3r = -\int_A\mathbf{v}\cdot\nabla\psi d^3r \end{align} \]
transferred to the discretization. For this calculation, the orientations $\gamma_{e(c)}$ of the vectors must also be taken into account. It will
\[ \begin{align} \gamma_{e(c)} \coloneqq \begin{cases} 1,\text{ falls der Vektor bei }e\text{ in Bezug auf }c\text{ nach außen zeigt,}\\ -1,\text{ sonst} \end{cases} \end{align} \]
defined. Thus you get
\[ \begin{align} &\int_A\psi\nabla\cdot\mathbf{v}d^3r \to \sum_{c, k}\left(\psi\nabla\cdot\mathbf{v}\right)_{c, k}A_c\Delta z_{c, k}\nonumber\\ &= \sum_{c, k}\Bigg[\sum_{e\in c}\left(l_e\Delta z_{e, k}\psi_{c, k}\gamma_{e(c)}u_e\right) + A_c\left(\psi_{c, k}w_{c, k - 1/2} - \psi_{c, k}w_{c, k + 1/2}\right)\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4}\left(\frac{\Delta z_{e, k - 1}}{\Delta z_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1} - \frac{\Delta z_{e, k + 1}}{\Delta z_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\right)\Bigg]. \end{align} \]
With an index shift follows
\[ \begin{align} &\sum_{c, k}\left(\psi\nabla\cdot\mathbf{v}\right)_{c, k}A_c\Delta z_{c, k} = \sum_{c, k}\Bigg[\sum_{e\in c}\left(l_e\Delta z_{e, k}\psi_{c, k}\gamma_{e(c)}u_e\right) + A_c\left(\psi_{c, k}w_{c, k - 1/2} - \psi_{c, k}w_{c, k + 1/2}\right)\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4}\left(\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k + 1}u_{e, k} - \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k - 1}u_{e, k}\right)\Bigg]\nonumber\\ &= \sum_{c, k}\Bigg[\sum_{e\in c}\left(l_e\Delta z_{e, k}\psi_{c, k}\gamma_{e(c)}u_e\right) + A_c\left(\psi_{c, k}w_{c, k - 1/2} - \psi_{c, k}w_{c, k + 1/2}\right)\nonumber\\ &- \sum_{e\in c}\frac{l_ed_e}{4}\left(\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k + 1}u_{e, k} - \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} - \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k - 1}u_{e, k} + \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k}u_{e, k}\right)\Bigg]\nonumber \end{align} \] \[ \begin{align} &= \sum_{c, k}\sum_{e\in c}l_e\Delta z_{e, k}\psi_{c, k}\gamma_{e(c)}u_e + \sum_{c, k + 1/2}A_c\left(\psi_{c, k}w_{c, k - 1/2} - \psi_{c, k}w_{c, k + 1/2}\right)\nonumber\\ &- \sum_{e, k}\sum_{c\in e}\frac{l_ed_e}{4}\left(\frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k + 1}u_{e, k} - \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} - \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k - 1}u_{e, k} + \frac{\Delta z_{e, k}}{\Delta z_{c, k}}J_{e, k}\psi_{c, k}u_{e, k}\right)\nonumber\\ &= \sum_{e, k}\sum_{c\in e}l_e\Delta z_{e, k}\psi_{c, k}\gamma_{e(c)}u_e - \sum_{c, k + 1/2}A_c\Delta z_{c, k + 1/2}w_{c, k + 1/2}\frac{\psi_{c, k} - \psi_{c, k + 1}}{\Delta z_{c, k + 1/2}}\nonumber\\ &- \sum_{e, k}\sum_{c\in e}\frac{l_ed_e}{4}\left(\frac{\Delta z_{e, k}}{\Delta z_{c, k}}\Delta z_{c, k + 1/2}J_{e, k}\frac{\psi_{c, k + 1} - \psi_{c, k}}{\Delta z_{c, k + 1/2}}u_{e, k} + \frac{\Delta z_{e, k}}{\Delta z_{c, k}}\Delta z_{c, k - 1/2}J_{e, k}u_{e, k}\frac{\psi_{c, k} - \psi_{c, k - 1}}{\Delta z_{c, k - 1/2}}\right)\nonumber\\ &= \sum_{e, k}-l_ed_e\Delta z_{e, k}u_e\sum_{c\in e}-\frac{\psi_{c, k}\gamma_{e(c)}}{d_e} - \sum_{c, k + 1/2}A_c\Delta z_{c, k + 1/2}w_{c, k + 1/2}\frac{\psi_{c, k} - \psi_{c, k + 1}}{\Delta z_{c, k + 1/2}}\nonumber \end{align} \] \[ \begin{align} &- \sum_{e, k}l_ed_e\Delta z_{e, k}u_{e, k}J_{e, k}\sum_{c\in e}\frac{1}{4}\left(\frac{\Delta z_{c, k + 1/2}}{\Delta z_{c, k}}\frac{\psi_{c, k + 1} - \psi_{c, k}}{\Delta z_{c, k + 1/2}} + \frac{\Delta z_{c, k - 1/2}}{\Delta z_{c, k}}\frac{\psi_{c, k} - \psi_{c, k - 1}}{\Delta z_{c, k - 1/2}}\right)\nonumber\\ &= \sum_{e, k}-l_ed_e\Delta z_{e, k}u_e\sum_{c\in e}-\frac{\psi_{c, k}\gamma_{e(c)}}{d_e} - \sum_{c, k + 1/2}A_c\Delta z_{c, k + 1/2}w_{c, k + 1/2}\left(\nabla_z\psi\right)_{c, k + 1/2}\nonumber\\ &+ \sum_{e, k}l_ed_e\Delta z_{e, k}u_{e, k}J_{e, k}\frac{1}{2}\newoverline{\left(\frac{\Delta z_{c, k - 1/2}}{\Delta z_{c, k}}\left(\nabla_z\psi\right)_{c, k - 1/2} + \frac{\Delta z_{c, k + 1/2}}{\Delta z_{c, k}}\left(\nabla_z\psi\right)_{c, k + 1/2}\right)}^{(e)}\nonumber \end{align} \] \[ \begin{align} &= \sum_{e, k}-l_ed_e\Delta z_{e, k}u_e\sum_{c\in e}-\frac{\psi_{c, k}\gamma_{e(c)}}{d_e} - \sum_{c, k + 1/2}A_c\Delta z_{c, k + 1/2}w_{c, k + 1/2}\left(\nabla_z\psi\right)_{c, k + 1/2}\nonumber\\ &+ \sum_{e, k}l_ed_e\Delta z_{e, k}u_{e, k}J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)}\nonumber\\ &= \sum_{e, k}-l_ed_e\Delta z_{e, k}u_e\left[\sum_{c\in e}\left(-\frac{\psi_{c, k}\gamma_{e(c)}}{d_e}\right) - J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)}\right] - \sum_{c, k + 1/2}A_c\Delta z_{c, k + 1/2}w_{c, k + 1/2}\left(\nabla_z\psi\right)_{c, k + 1/2}\nonumber\\ &= -\left\{2\sum_{e, k}V_{e, k}u_e\left[\sum_{c\in e}\left(-\frac{\psi_{c, k}\gamma_{e(c)}}{d_e}\right) - J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)}\right] + \sum_{c, k + 1/2}V_{c, k + 1/2}w_{c, k + 1/2}\left(\nabla_z\psi\right)_{c, k + 1/2}\right\}. \end{align} \]
The factor 2 stands for the two spatial directions in the horizontal. This results in the discretization
\[ \begin{align} \left(\nabla\psi\right)_{e, k} = \frac{1}{d_e}\sum_{c\in e}\left(-\psi_{c, k}\gamma_{e(c)}\right) - J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)} \end{align} \]
of the gradient in terrain-following coordinates.
Even in the deep atmosphere, the discretized product rule of the Nabla operator is not locally valid in coordinates following the terrain. Here, according to Eq. (33.71)
\[ \begin{align} \int_A\psi\nabla\cdot\mathbf{v}d^3r&\to\sum_{c, k}\left(\psi\nabla\cdot\mathbf{v}\right)_{c, k}V_{c, k} = \sum_{c, k}\Bigg\lbrace\left[\sum_{e\in c}\left(A_{e, k}\psi_{c, k}u_{e, k}\right) + A_{c, k - 1/2}\psi_{c, k}w_{c, k - 1/2} - A_{c, k + 1/2}\psi_{c, k}w_{c, k + 1/2}\right]\nonumber\\ &+ \sum_{e\in c}\Bigg[-A_{c, k - 1/2}\left(\frac{A_{e, k}d_{e, k}}{4V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{A_{e, k - 1}d_{e, k - 1}}{4V_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1}\right)\nonumber\\ &+ A_{c, k + 1/2}\left(\frac{A_{e, k}d_{e, k}}{4V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{A_{e, k + 1}d_{e, k + 1}}{4V_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\right)\Bigg]\Bigg\rbrace. \end{align} \]
Except for the correction terms for the coordinates following the terrain, the product rule of the Nabla operator is also locally valid here; the derivation in Sect. 26.7.5 still applies. In addition, a global integral disappears over a divergence here too, the reason for this is the same as in the shallow atmosphere. One can therefore limit the partial integration to the terms of the terrain-following correction:
\[ \begin{align} &\sum_{c, k}\sum_{e\in c}\Bigg[-A_{c, k - 1/2}\left(\frac{A_{e, k}d_{e, k}}{4V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{A_{e, k - 1}d_{e, k - 1}}{4V_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1}\right)\nonumber\\ &+ A_{c, k + 1/2}\left(\frac{A_{e, k}d_{e, k}}{4V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{A_{e, k + 1}d_{e, k + 1}}{4V_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\right)\Bigg]\nonumber\\ &= \sum_{c, k}\sum_{e\in c}\Bigg[-A_{c, k - 1/2}\left(\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{V_{e, k - 1}}{2V_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1}\right)\nonumber\\ &+ A_{c, k + 1/2}\left(\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + \frac{V_{e, k + 1}}{2V_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\right)\Bigg]\nonumber\\ &= \sum_{c, k}\sum_{e\in c}\Bigg[-A_{c, k - 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} - A_{c, k - 1/2}\frac{V_{e, k - 1}}{2V_{c, k - 1}}J_{e, k - 1}\psi_{c, k}u_{e, k - 1}\nonumber\\ &+ A_{c, k + 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + A_{c, k + 1/2}\frac{V_{e, k + 1}}{2V_{c, k + 1}}J_{e, k + 1}\psi_{c, k}u_{e, k + 1}\Bigg]\nonumber \end{align} \] \[ \begin{align} &= \sum_{c, k}\sum_{e\in c}\Bigg(-A_{c, k - 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} - A_{c, k + 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k + 1}u_{e, k}\nonumber\\ &+ A_{c, k + 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k}u_{e, k} + A_{c, k - 1/2}\frac{V_{e, k}}{2V_{c, k}}J_{e, k}\psi_{c, k - 1}u_{e, k}\Bigg)\nonumber\\ &= \sum_{c, k}\sum_{e\in c}J_{e, k}u_{e, k}\left(-\frac{A_{c, k - 1/2}}{2V_{c, k}}\psi_{c, k} - \frac{A_{c, k + 1/2}}{2V_{c, k}}\psi_{c, k + 1} + \frac{A_{c, k + 1/2}}{2V_{c, k}}\psi_{c, k} + \frac{A_{c, k - 1/2}}{2V_{c, k}}\psi_{c, k - 1}\right)V_{e, k}\nonumber\\ &= \sum_{c, k}\sum_{e\in c}J_{e, k}u_{e, k}\left(\frac{A_{c, k - 1/2}}{2V_{c, k}}\Delta z_{c, k - 1/2}\frac{\psi_{c, k - 1} - \psi_{c, k}}{\Delta z_{c, k - 1/2}} + \frac{A_{c, k + 1/2}}{2V_{c, k}}\Delta z_{c, k + 1/2}\frac{\psi_{c, k} - \psi_{c, k + 1}}{\Delta z_{c, k + 1/2}}\right)V_{e, k}\nonumber \end{align} \] \[ \begin{align} &= \sum_{c, k}\sum_{e\in c}J_{e, k}u_{e, k}\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}V_{e, k} = -\sum_{c, k}\sum_{e\in c}-J_{e, k}u_{e, k}\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}V_{e, k}\nonumber\\ &= -\sum_{c, k}\sum_{e\in c}-J_{e, k}2u_{e, k}\frac{1}{2}\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}V_{e, k} = -\sum_{e, k}J_{e, k}u_{e, k}V_{e, k}\left(-2\sum_{c\in e}\frac{1}{2}\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}\right)\nonumber\\ &= -2\sum_{e, k}u_{e, k}V_{e, k}\left(-J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)}\right) \end{align} \]
Therefore this also applies in the deep atmosphere
\[ \begin{align} \left(\nabla\psi\right)_{e, k} = \frac{1}{d_{e, k}}\sum_{c\in e}\left(-\psi_{c, k}\gamma_{e(c)}\right) - J_{e, k}\newoverline{\newoverline{\left(\nabla_z\psi\right)_{c, k + 1/2}}^{(k)}}^{(e)}, \end{align} \]
where only this order of the averaging operators is permitted here.
\[ \begin{align} \sum_{e, k}\newoverline{\rho}_c^{(e)}v_e\frac{\partial v_{e, k}}{\partial t} &= \sum_{e, k}\newoverline{\rho}_c^{(e)}v_e\newoverline{\newoverline{\rho w\cdot\eta_t}^{(e)}}^{(k)} \end{align} \]
Since the Hamilton function does not explicitly depend on time, it is equal to the total energy. If you put in Eq. (10.66) $F = H$ is obtained
\[ \begin{align} \frac{dH}{dt} = \lbrace H, H\rbrace = \lbrace H, Z_a, H\rbrace + \lbrace H, H\rbrace_\rho + \lbrace H, H\rbrace_{\newtilde{s}}. \end{align} \]
Still apply
\[ \begin{align} \lbrace H, Z_a, H\rbrace \stackrel{\href{ch-09-application-of-the-hamilton-formalism-to.html#eq:poisson_bracket_atm_gen_deriv_5}{\text{Glg. (10.54)}}}{=} 0, & {} & \lbrace H, H\rbrace_\rho \stackrel{\href{ch-09-application-of-the-hamilton-formalism-to.html#eq:poisson_bracket_atm_gen_deriv_3}{\text{Glg. (10.50)}}}{=} 0, & {} & \lbrace H, H\rbrace_{\newtilde{s}} \stackrel{\href{ch-09-application-of-the-hamilton-formalism-to.html#eq:poisson_bracket_atm_gen_deriv_4}{\text{Glg. (10.51)}}}{=} 0. \end{align} \]
This implies conservation of energy
\[ \begin{align} H = \text{const.} \end{align} \]
In order to transfer energy conservation to spatial discretization, three requirements must be taken into account:
The first point corresponds to the fact that the generalized Coriolis force does no work, which was already discussed in Chap. 28 was detected. The second and third statements have the same mathematical requirement. Add to this the equations (10.78) - (10.79), which results in
\[ \begin{align} \int_A\underbrace{-\rho\mathbf{v}\cdot\nabla\phi}_{\text{A}}\underbrace{ - c_d^{(p)}\rho\theta\mathbf{v}\cdot\nabla\Pi}_{\text{B}}\underbrace{ - \phi\nabla\cdot\left(\rho\mathbf{v}\right)}_{\text{A}}\underbrace{ - c_d^{(p)}\Pi\nabla\cdot\left(\rho\theta\mathbf{v}\right)}_{\text{B}}d^3r = 0\tag{34.58}\label{eq:poisson_energy_deriv} \end{align} \]
leads. The integrals over the terms A and B cancel each other out. The following two facts were used:
Both facts also apply in discretization. The evidence for this has already been provided:
Another aspect should be pointed out here: one discretizes Eq. (34.58), $\theta$ only occurs at the edges of the lattice, not in the cell centers. You can freely choose the value of $\theta$ on the edge, as long as you also do this in the pressure gradient acceleration. This will be used in Section 34.4 to calculate $\theta$ on the edge in such a way that the accuracy of the flow operator $-\nabla\cdot\left(\rho\theta\mathbf{v}\right)$ is increased. The $\theta-$ conservation itself is not violated by this, since what flows out of one grid box continues to flow into the neighboring one.
Pressure gradient and gravitational acceleration as well as the generalized Coriolis force are covered here. The gradient of kinetic energy, however, cancels out in the Poisson brackets because, although it transports kinetic energy, it does not create or destroy any energy. This will be repeated here. To see this, write down the momentum equation in the form
\[ \begin{align} \frac{\partial\mathbf{v}}{\partial t} = -\nabla k + \mathbf{l}, \end{align} \]
where pressure gradient and gravitational acceleration were combined in the term $\mathbf{l}$. The Coriolis force can be ignored here. Scalar multiplication by $\mathbf{v}$ yields
\[ \begin{align} \mathbf{v}\cdot\frac{\partial\mathbf{v}}{\partial t} = -\mathbf{v}\cdot\nabla k + \mathbf{v}\cdot\mathbf{l}\nonumber & {} & \Leftrightarrow\frac{\partial k}{\partial t} = -\mathbf{v}\cdot\nabla k + \mathbf{v}\cdot\mathbf{l}. \end{align} \]
Multiplying this by the mass density $\rho$ gives
\[ \begin{align} \rho\frac{\partial k}{\partial t} &= -\rho\mathbf{v}\cdot\nabla k + \rho\mathbf{v}\cdot\mathbf{l}.\tag{34.60}\label{eq:no-split-explicit_deriv_0} \end{align} \]
The continuity equation is:
\[ \begin{align} \frac{\partial\rho}{\partial t} = -\nabla\cdot\left(\rho\mathbf{v}\right). \end{align} \]
Multiplying this by the specific kinetic energy $k$ gives
\[ \begin{align} k\frac{\partial\rho}{\partial t} = -k\nabla\cdot\left(\rho\mathbf{v}\right).\tag{34.62}\label{eq:no-split-explicit_deriv_1} \end{align} \]
If you add the equations (34.60) and (34.62), you get
\[ \begin{align} \frac{\partial\left(\rho k\right)}{\partial t} = -\nabla\cdot\left(\rho k\mathbf{v}\right) + \mathbf{v}\cdot\mathbf{l}.\tag{34.63}\label{eq:no-split-explicit_deriv_2} \end{align} \]
This means that only the global integral over $\mathbf{v}\cdot\mathbf{l}$ contributes to the evolution of the total kinetic energy of the atmosphere. The derivation only includes calculation rules for $\frac{\partial}{\partial t}$ and the rule $\nabla\cdot\left(\psi\mathbf{v}\right) = \mathbf{v}\cdot\nabla\psi + \psi\nabla\cdot\mathbf{v}$, which is already considered to be fulfilled. The acceleration $-\nabla k$ does not influence the energetic integrals even in the discretization.
Thus, the location discretization developed in this chapter is energetically self-consistent.
A hydrostatic background state as described in Sect. 13.7.5 can be used without destroying energy conservation. It must be ensured that:
\[ \begin{align} \mathbf{0} = -c^{(p)}\newoverline{\theta}\nabla\newoverline{\Pi} - \nabla\phi \end{align} \]
also applies discreetly. The corresponding terms then have each other in the discretization.
However, this is only the case without orography. However, if mountains are present, the horizontal components of the gradients do not cancel each other out exactly, resulting in a small error in the conservation of energy. However, in the interaction of spatial and time discretization, the advantages of the hydrostatic background state outweigh the disadvantages. The most important advantage here is that the horizontal orographic correction components of $-c^{(p)}\theta\nabla\Pi$ and $-\nabla\phi$ are two quantities that almost completely cancel each other out, i.e
\[ \begin{align} -c^{(p)}\newoverline{\theta_{c, k}}^{(e)}J_{e, k}\newoverline{\left(\nabla_z\Pi\right)_{c, k + 1/2}}^{(k)} - J_{e, k}\newoverline{\left(\nabla_z\phi\right)_{c, k + 1/2}}^{(k)} \approx 0. \end{align} \]
The remaining difference is largely a discretization error. This becomes smaller if you use a hydrostatic background state, since then only
\[ \begin{align} -c^{(p)}\newoverline{\theta_{c, k}}^{(e)}J_{e, k}\newoverline{\left(\nabla_z\Pi'\right)_{c, k + 1/2}}^{(k)} - c^{(p)}\newoverline{\theta'_{c, k}}^{(e)}J_{e, k}\newoverline{\left(\nabla_z\newoverline{\Pi}\right)_{c, k + 1/2}}^{(k)} \end{align} \]
must be calculated.
The kinematic boundary condition $\mathbf{v}\cdot\mathbf{n} = 0$ is used as the upper boundary condition in the context of reversible dynamics. This corresponds to the model
\[ \begin{align} w_{c, 1/2} = 0 \end{align} \]
for all cells $c$.
The lower boundary condition in the context of reversible dynamics is the no-slip boundary condition $\mathbf{v} = \mathbf{0}$. Since there are only vertical vector components there, this corresponds
\[ \begin{align} w_{c, N_L + \frac{1}{2}} = 0 \end{align} \]
for all cells $c$.
Within the framework of the formalism developed so far, there are three possible approaches to increasing the convergence order:
These remaining freedoms are not sufficient to increase the convergence order of the overall model.
The accuracy of the discretization can be increased by not using $\theta$ via the arithmetic mean
\[ \begin{align} \newtilde{\theta}^{(e)} = \frac{\theta_{c_1(e)} + \theta_{c_2(e)}}{2}, \end{align} \]
what is second order is interpolated to the edge, but in a more complex way. The vertical index $k$ is omitted here because two-dimensionality is assumed in this section and $c_1(e)$, $c_2(e)$ are, as usual, the two lines that border the edge $e$.
First of all, this is viewed one-dimensionally. To do this, one assumes the true location dependence $\theta = \theta\left(x\right)$ and assumes that the values of $\theta$ are known in the scalar grid points $i$, which are positioned at locations $x_i = i\Delta x$, i.e
\[ \begin{align} \theta_i = \theta\left(x_i\right). \end{align} \]
A suitable approximation for $\newtilde{\theta}^{(e)}$ is now sought, i.e. $\theta$ at position $x = \frac{\Delta x}{2}$. To do this, first replace the location dependence $\theta\left(x\right)$ by $\mu^{(n)}\left(x\right)$, where $n$ stands for the convergence order.
For $n = 2$ you do the approach
\[ \begin{align} \mu^{(2)}\left(0\right) &= \theta\left(0\right),\tag{34.71}\label{eq:skam_gass_deriv_1}\\ \mu^{(2)}\left(\Delta x\right) &= \theta\left(\Delta x\right).\tag{34.72}\label{eq:skam_gass_deriv_2} \end{align} \]
Since you have two linear equations for $\mu^{(2)}$, you do an approach with two coefficients, i.e
\[ \begin{align} \mu^{(2)}\left(x\right) &= a + bx. \end{align} \]
If you insert this into the equations (34.71) - (34.72), you get
\[ \begin{align} a &= \theta\left(0\right),\\ a + b\Delta x &= \theta\left(\Delta x\right). \end{align} \]
This leads to
\[ \begin{align} a &= \theta_{c_1(e)},\\ b &= \frac{\theta\left(\Delta x\right) - \theta\left(0\right)}{\Delta x} = \frac{\theta_{c_2(e)} - \theta_{c_1(e)}}{\Delta x}. \end{align} \]
In order to calculate $\newtilde{\theta}^{(e,2)}$ from this, where the superscript two stands for the order, an integration must be carried out over the time interval $\left[0, \Delta t\right]$, where $\Delta t$ is the length of a time step. You calculate:
\[ \begin{align} \newtilde{\theta}^{(e,2)} &= a + b\frac{\Delta x}{2} = \theta_{c_1(e)} + \frac{\theta\left(\Delta x\right) - \theta\left(0\right)}{2} = \frac{\theta\left(0\right) + \theta\left(\Delta x\right)}{2} \end{align} \]
To obtain a third-order approximation ($n = 3$), one includes the second derivative of $\theta$ in the grid cell $c_1\left(e\right)$ in the reconstruction of $\theta$, i.e
\[ \begin{align} \mu^{(3)}\left(x\right) &= a + bx + cx^2. \end{align} \]
Follow from this
\[ \begin{align} f_0 &= a,\\ f_1 &= a + b\Delta x + c\Delta x^2,\\ f_{-1} &= a - b\Delta x + c\Delta x^2, \end{align} \]
what is happening
\[ \begin{align} a &= f_0,\\ f_{-1} + f_1 &= 2a + 2c\Delta x^2 = 2\left(f_0 + c\Delta x^2\right) \Rightarrow \frac{f_{-1} + f_1}{2} - f_0 = c\Delta x^2 \Rightarrow c = \frac{f_{-1} - 2f_0 + f_1}{2\Delta x^2},\\ f_1 - f_{-1} &= 2b\Delta x \Rightarrow b = \frac{f_1 - f_{-1}}{2\Delta x} \end{align} \]
can be reshaped. From this it follows
\[ \begin{align} \mu^{(3)}\left(x\right) &= f_0 + \frac{f_1 - f_{-1}}{2\Delta x}x + \frac{f_{-1} - 2f_0 + f_1}{2\Delta x^2}x^2. \end{align} \]
The value of this function on the edge ($x = \frac{\Delta x}{2}$) is therefore
\[ \begin{align} \mu^{(3)}\left(\frac{\Delta x}{2}\right) &= f_0 + \frac{f_1 - f_{-1}}{4} + \frac{f_{-1} - 2f_0 + f_1}{8} = \frac{8f_0 + 2f_1 - 2f_{-1} + f_{-1} - 2f_0 + f_1}{8}\nonumber\\ &= \frac{6f_0 + 3f_1 - f_{-1}}{8} = \frac{4f_0 + 4f_1 + 2f_0 - f_1 - f_{-1}}{8} = \frac{f_0 + f_1}{2} + \frac{2f_0 - f_1 - f_{-1}}{8}\nonumber\\ &= \frac{f_0 + f_1}{2} - \frac{f_1 - 2f_0 + f_{-1}}{8} = \frac{f_0 + f_1}{2} - \Delta x^2\frac{f_1 - 2f_0 + f_{-1}}{8\Delta x^2}\nonumber\\ &= \frac{f_0 + f_1}{2} - \frac{\Delta x^2}{8}\delta_x^2f \end{align} \]
with
\[ \begin{align} \delta_x^2f \coloneqq \frac{f_1 - 2f_0 + f_{-1}}{\Delta x^2} \end{align} \]
as an approximation for the second derivative of $\theta$ in the $x-$direction at the center of the cell $c_1\left(e\right)$.
In two dimensions, the one-dimensional formulas are simply applied in every spatial direction. This is not a problem on quadrangular grids where all vectors are parallel to a global coordinate axis. On grids that do not meet this criterion, this is more complicated, which is what the next section will cover.
On the hexagonal grid, the vectors are not aligned along coordinate lines, which complicates the calculation of the second derivatives at the cell centers. Therefore, one does the approach for $\theta$
\[ \begin{align} g = g\left(x,y\right) = c_1 + c_2x + c_3y + c_4x^2 + c_5xy + c_6y^2. \end{align} \]
$x$ and $y$ are calculated on the tangent plane to the sphere in the center of the cell under consideration. The $x-$axis points in the direction of the edge on which $\newtilde{\theta}^{(e)}$ is to be calculated. Now define the coefficient vector
\[ \begin{align} \mathbf{f} \coloneqq \left(\begin{array}{c} c_1\\ c_2\\ c_3\\ c_4\\ c_5\\ c_6 \end{array}\right). \end{align} \]
Only $c_4$ is relevant for the second derivative in the $x-$direction at the center of the cell, it applies
\[ \begin{align} \frac{\partial^2 g}{\partial x^2}\left(0,0\right) = 2c_4.\tag{34.91}\label{eq:adv_hex_deriv} \end{align} \]
The goal now is to calculate the vector $\mathbf{f}$. To do this, you first define the vector
\[ \begin{align} \mathbf{s} \coloneqq \left(\begin{array}{c} \theta\left(x_1,y_1\right)\\ \vdots\\ \theta\left(x_m,y_m\right) \end{array}\right) \end{align} \]
the actual values of $\theta$ in the cell centers. Here $m \geq 6$ is the number of points that are calculated to calculate $\mathbf{f}$, these are the cell in the center of which the second derivative is to be calculated and all neighboring cells. Now define an error function
\[ \begin{align} \psi\left(\mathbf{f}\right) \coloneqq \sum_i^m\left[g\left(x_i,y_i\right) - s_i\right]^2. \end{align} \]
Now define the polynomial matrix
\[ \begin{align} P \coloneqq \left(\begin{array}{cccccc} 1 & x_1 & y_1 & x_1^2 & x_1y_1 & y_1^2\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & x_m & y_m & x_m^2 & x_my_m & y_m^2 \end{array}\right). \end{align} \]
This allows $\psi$ in the form
\[ \begin{align} \psi\left(\mathbf{f}\right) &= \left(P\mathbf{f} - \mathbf{s}\right)^2 = \left(P\mathbf{f}\right)\cdot\left(P\mathbf{f}\right) - 2\left(P\mathbf{f}\right)\cdot\mathbf{s} + \mathbf{s}^2 = \sum_{i=1}^m\left(P\mathbf{f}\right)_i^2 - 2\sum_{i=1}^m\left(P\mathbf{f}\right)_is_i + \sum_{i=1}^ms_i^2\nonumber\\ &= \sum_{i=1}^m\left(\sum_{j=1}^nP_{i, j}f_j\right)^2 - 2\sum_{i=1}^m\left(\sum_{j=1}^nP_{i, j}f_j\right)s_i + \sum_{i=1}^ms_i^2 \end{align} \]
note down. From this it follows
\[ \begin{align} \frac{\partial\psi\left(\mathbf{f}\right)}{\partial f_k} &= 2\sum_{i=1}^m\left(\sum_{j=1}^nP_{i, j}\delta_{j,k}\right)^2\left(\sum_{j=1}^nP_{i, j}f_j\right) - 2\sum_{i=1}^m\left(\sum_{j=1}^nP_{i, j}\delta_{j,k}\right)s_i = 0\nonumber\\ \Leftrightarrow\frac{\partial\psi\left(\mathbf{f}\right)}{\partial f_k} &= \sum_{i=1}^m\left(\sum_{j=1}^nP_{i,j}\delta_{j,k}\right)\left(\sum_{j=1}^nP_{i, j}f_j\right) - \sum_{i=1}^m\left(\sum_{j=1}^nP_{i, j}\delta_{j,k}\right)s_i = 0\nonumber\\ \Leftrightarrow\frac{\partial\psi\left(\mathbf{f}\right)}{\partial f_k} &= \sum_{i=1}^mP_{i,k}\left(\sum_{j=1}^nP_{i, j}f_j\right) - \sum_{i=1}^mP_{i,k}s_i = \sum_{i=1}^m\sum_{j=1}^nP_{i,k}P_{i, j}f_j - \sum_{i=1}^mP_{i, k}s_i = 0\nonumber\\ \Leftrightarrow\frac{\partial\psi\left(\mathbf{f}\right)}{\partial f_k} &= \sum_{j=1}^n\sum_{i=1}^mP^T_{k,i}P_{i,j}f_j - \sum_{i=1}^mP^T_{k,i}s_i = \sum_{j=1}^n\left(P^TP\right)_{k,j}f_j - \sum_{i=1}^mP^T_{k,i}s_i = 0\nonumber\\ \Leftrightarrow\frac{\partial\psi\left(\mathbf{f}\right)}{\partial f_k} &= \left[\left(P^TP\right)\cdot\mathbf{f}\right]_k - \left(P^T\cdot\mathbf{s}\right)_k = 0. \end{align} \]
So the equivalence applies
\[ \begin{align} \nabla_\mathbf{k}\psi = 0 \Leftrightarrow \left(P^TP\right)\cdot\mathbf{f} = P^T\cdot\mathbf{s}, \end{align} \]
what to do
\[ \begin{align} \mathbf{f} = \left[\left(P^TP\right)^{-1}P^T\right]\cdot\mathbf{s} = B\cdot\mathbf{s} \end{align} \]
with
\[ \begin{align} B \coloneqq \left[\left(P^TP\right)^{-1}P^T\right] \end{align} \]
can transform. From this we can use Eq. (34.91) calculate the second derivatives in the cell centers in the x-direction.
There is the following conflict between explicit and implicit method:
It follows that a fast implicit method would be ideal for a dynamic core. However, since the equations to be solved are nonlinear, a three-dimensional implicit method is only possible via an iterative approach. In this section, the compromise chosen is the HEVI concept,, which stands for horizontally explicit, vertically implicit.
A significant limitation with regard to the time step and thus also with regard to the efficiency of a model is the CFL criterion associated with sound waves. Therefore, the system of equations is often divided into so-called fast modes or divergent modes on the one hand and the slow modes, advective modes or non-divergent modes on the other hand. The fast modes include all terms associated with sound waves, while the slow modes include the remaining terms. With split-explicit techniques, the fast modes are now integrated with a time step that is approximately three times smaller. In the case of the equation system Equations (34.4) - (34.7), only the advection of the momentum for the larger time step would come into question. The advections of density and entropy, on the other hand, should be counted among the fast modes in order not to have to break down the divergences of the flow density into advection and compression components. However, if one places high demands on the energetic self-consistency of a time step method, splitting off momentum advection does not make sense.
In order to see this, one must first realize that the advection of kinetic energy is connected to the continuity equation and is therefore one of the fast modes. This has already been seen in Section 34.2.4. Therefore, both terms or equations must be integrated with the same time step. In order to preserve the discretized analogue of the Lamb transformation $-\left(\mathbf{v}\cdot\nabla\right)\mathbf{v} = -\nabla k - \left(\nabla\times\mathbf{v}\right)\times\mathbf{v}$, the generalized Coriolis term and thus the complete system of equations must then also be integrated with the smaller time step.
Due to the disadvantages discussed here, no split-explicit option is built into the dynamic core being developed.
Vertically propagating waves are treated implicitly for reasons of stability. To do this, you first make a note
\[ \begin{align} \theta'' &\stackrel{\theta = \newtilde{\theta}/\rho}{=} \frac{\partial\theta}{\partial\rho}\rho' + \frac{\partial\theta}{\partial\newtilde{\theta}}\newtilde{\theta}' = -\frac{\newtilde{\theta}}{\rho^2}\rho' + \frac{1}{\rho}\newtilde{\theta}',\\ \Pi'' &\stackrel{\href{ch-08-first-law-in-the-atmosphere.html#eq:exner_pressure_diag}{\text{Glg. (9.71)}}}{=} \frac{R_s}{c^{(V)}\newtilde{\theta}}\left(\frac{R_s\newtilde{\theta}}{p_0}\right)^{R_s/c^{(V)}}\newtilde{\theta}'. \end{align} \]
The two lines represent linear developments compared to the previous time step and not differences to the hydrostatic background state, which are denoted by $\theta',\Pi'$. The derivatives are abbreviated
\[ \begin{align} \alpha \coloneqq -\frac{\newtilde{\theta}}{\rho^2},& \beta \coloneqq \frac{1}{\rho},\\ \gamma \coloneqq \frac{R_s}{c^{(V)}\newtilde{\theta}}\Pi \end{align} \]
away. These are calculated in the cell centers, using the old time step (time step $n$) for the predictor step and the arithmetic mean of the steps $n$ and $n + 1$ for the corrector step. The vertical momentum equation is written down in the form
\[ \begin{align} w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left[\left(\newoverline{\newoverline{\theta}_k}^{(k + 1/2)} + \newoverline{\theta_k^{'(n + 1)^\star}}^{(k + 1/2)}\right)\frac{\Pi^{'(n + 1)}_{k} - \Pi^{'(n + 1)}_{k + 1}}{z_{k} - z_{k + 1}} + \newoverline{\theta_k^{'(n + 1)}}^{(k + 1/2)}\frac{\newoverline{\Pi}_{k} - \newoverline{\Pi}_{k + 1}}{z_{k} - z_{k + 1}}\right]\nonumber\\ &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(n + 1)}_{k} - \Pi^{'(n + 1)}_{k + 1}}{z_{k} - z_{k + 1}} + \newoverline{\theta_k^{'(n + 1)}}^{(k + 1/2)}\frac{\newoverline{\Pi}_{k} - \newoverline{\Pi}_{k + 1}}{z_{k} - z_{k + 1}}\right)\nonumber\\ &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(n + 1)}_{k} - \Pi^{'(n + 1)}_{k + 1}}{z_{k} - z_{k + 1}} + \newoverline{\theta_k^{'(n + 1)}}^{(k + 1/2)}\frac{d\newoverline{\Pi}}{dz}\vline_{k+1/2}\right).\tag{34.104}\label{eq:wave_vertical_deriv_0} \end{align} \]
Here $\epsilon$ is the implicit weight, for which $\frac{c^{(v)}}{c^{(p)}}$ is usually used based on the findings in Section 25.5.1.3. $\left(n + 1\right)^\star$ is a preliminary value for the new time step. Now write down the implicit quantities
\[ \begin{align} \theta^{'(n + 1)}_k &= \theta^{'(\text{expl.})}_k + \theta''_k = \theta^{'(\text{expl.})}_k + \alpha_k\rho'_k + \beta_k\newtilde{\theta}'_k,\\ \Pi^{'(n + 1)}_k &= \Pi^{'(\text{expl.})}_k + \Pi''_k = \Pi^{'(\text{expl.})}_k + \gamma_k\newtilde{\theta}'_{k}. \end{align} \]
The deleted sizes are the deviations from the explicit components, i.e
\[ \begin{align} \psi'_k \coloneqq \psi_k^{(n + 1)} - \psi_k^{(\text{expl.})} \end{align} \]
for $\psi = \rho, \newtilde{\theta}.$ This applies
\[ \begin{align} \rho_k' = -\Delta t\frac{M_{k - 1/2} - M_{k + 1/2}}{V_k},\\ \newtilde{\theta}_k' = -\Delta t\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}}{V_k}. \end{align} \]
Defined here
\[ \begin{align} M_{k + 1/2} \coloneqq A_{k + 1/2}\frac{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}w_{k + 1/2}^{(n + 1)} + \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}w_{k + 1/2}^{(n)}}{2}\tag{34.110}\label{eq:def_ml} \end{align} \]
for $1 \leq k \leq N_L - 1$ the vertical covariant mass flow density multiplied by the boundary surfaces $A_{k + 1/2}$ of the grid boxes, this is the unknown vector. Putting this into Eq. (34.104), you get
\[ \begin{align} w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} + \Pi''_{k} - \Pi^{'(\text{expl.})}_{k + 1} - \Pi''_{k + 1}}{z_{k} - z_{k + 1}} + \newoverline{\theta_k^{'(n + 1)}}^{(k + 1/2)}\frac{d\newoverline{\Pi}}{dz}\vline_{k+1/2}\right)\nonumber\\ \Leftrightarrow w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} + \Pi''_{k} - \Pi^{'(\text{expl.})}_{k + 1} - \Pi''_{k + 1}}{z_{k} - z_{k + 1}} + \frac{\theta_{k}^{'(n + 1)}+\theta_{k + 1}^{'(n + 1)}}{2}\frac{d\newoverline{\Pi}}{dz}\vline_{k+1/2}\right)\nonumber\\ \Leftrightarrow w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} + \gamma_{k}\newtilde{\theta}'_{k} - \Pi^{'(\text{expl.})}_{k + 1} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}}{z_{k} - z_{k + 1}} + \frac{\theta_{k}^{'(n + 1)}+\theta_{k + 1}^{'(n + 1)}}{2}\frac{d\newoverline{\Pi}}{dz}\vline_{k+1/2}\right)\nonumber\\ \Leftrightarrow w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} + \gamma_{k}\newtilde{\theta}'_{k} - \Pi^{'(\text{expl.})}_{k + 1} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}}{z_{k} - z_{k + 1}} + \frac{\theta_{k}^{'(\text{expl.})}+\theta''_{k}+\theta_{k + 1}^{'(\text{expl.})}+\theta''_{k + 1}}{2}\frac{d\newoverline{\Pi}}{dz}\vline_{k+1/2}\right)\\ \Leftrightarrow w_{k + 1/2}^{(n + 1)} &= w_{k + 1/2}^{(\text{expl.})} - \epsilon\Delta tc^{(p)}\cdot\Bigg(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} + \gamma_{k}\newtilde{\theta}'_{k} - \Pi^{'(\text{expl.})}_{k + 1} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}}{z_{k} - z_{k + 1}}\nonumber\\ &+ \frac{\theta_{k}^{'(\text{expl.})}+\alpha_{k}\rho'_{k}+\beta_{k}\newtilde{\theta}'_{k}+\theta_{k + 1}^{'(\text{expl.})}+\alpha_{k + 1}\rho'_{k + 1}+\beta_{k + 1}\newtilde{\theta}'_{k + 1}}{2}\frac{d\newoverline{\Pi}}{dz}\vline_{k+1/2}\Bigg). \end{align} \]
The deleted quantities contain the unknown vector, so put them on the left side:
\[ \begin{align} &\epsilon\Delta tc^{(p)}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\gamma_{k}\newtilde{\theta}'_{k} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}}{z_{k} - z_{k + 1}} + \frac{\alpha_{k}\rho'_{k}+\beta_{k}\newtilde{\theta}'_{k}+\alpha_{k + 1}\rho'_{k + 1}+\beta_{k + 1}\newtilde{\theta}'_{k + 1}}{2}\frac{d\newoverline{\Pi}}{dz}\vline_{k+1/2}\right)\nonumber\\ &= w_{k + 1/2}^{(\text{expl.})} - w_{k + 1/2}^{(n + 1)} - \epsilon\Delta tc^{(p)}\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{z_{k} - z_{k + 1}} - \epsilon\Delta tc^{(p)}\frac{\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}}{2}\frac{d\newoverline{\Pi}}{dz}\vline_{k+1/2}\nonumber\\ &\Leftrightarrow\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\gamma_{k}\newtilde{\theta}'_{k} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}}{z_{k} - z_{k + 1}} + \frac{\alpha_{k}\rho'_{k}+\beta_{k}\newtilde{\theta}'_{k}+\alpha_{k + 1}\rho'_{k + 1}+\beta_{k + 1}\newtilde{\theta}'_{k + 1}}{2}\frac{d\newoverline{\Pi}}{dz}\vline_{k+1/2}\nonumber\\ &= \frac{w_{k + 1/2}^{(\text{expl.})} - w_{k + 1/2}^{(n + 1)}}{\epsilon\Delta tc^{(p)}} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{z_{k} - z_{k + 1}} - \frac{\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}}{2}\frac{d\newoverline{\Pi}}{dz}\vline_{k+1/2}\nonumber\\ &\Leftrightarrow\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\left(\gamma_{k}\newtilde{\theta}'_{k} - \gamma_{k + 1}\newtilde{\theta}'_{k + 1}\right) + \frac{d\newoverline{\Pi}}{2dz}\vline_{k+1/2}\left(z_{k} - z_{k + 1}\right)\left(\alpha_{k}\rho'_{k}+\beta_{k}\newtilde{\theta}'_{k}+\alpha_{k + 1}\rho'_{k + 1}+\beta_{k + 1}\newtilde{\theta}'_{k + 1}\right)\nonumber\\ &= \frac{\left(w_{k + 1/2}^{(\text{expl.})} - w_{k + 1/2}^{(n + 1)}\right)\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta tc^{(p)}} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\left(\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}\right) - \left(z_{k} - z_{k + 1}\right)\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vline_{k+1/2}\nonumber\\ &\Leftrightarrow\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\Bigg(-\gamma_{k}\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}}{V_{k}}\nonumber\\ &+ \gamma_{k + 1}\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}}{V_{k + 1}}\Bigg)\nonumber\\ &+ \frac{d\newoverline{\Pi}}{2dz}\vline_{k+1/2}\left(z_{k} - z_{k + 1}\right)\Bigg(-\alpha_{k}\frac{M_{k - 1/2} - M_{k + 1/2}}{V_{k}}-\beta_{k}\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}}{V_{k}}\nonumber\\ &-\alpha_{k + 1}\frac{M_{k + 1/2} - M_{k + 3/2}}{V_{k + 1}}-\beta_{k + 1}\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}}{V_{k + 1}}\Bigg)\nonumber\\ &= \frac{\left(w_{k + 1/2}^{(\text{expl.})} - w_{k + 1/2}^{(n + 1)}\right)\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} - \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vline_{k+1/2}. \end{align} \]
Now take the replacements
\[ \begin{align} \alpha_k \to \frac{\alpha_k}{V_k}, \beta_k \to \frac{\beta_k}{V_k}, \gamma_k \to \frac{\gamma_k}{V_k} \end{align} \]
before. It follows
\[ \begin{align} &\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\Bigg[-\gamma_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &+ \gamma_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &+ \frac{d\newoverline{\Pi}}{2dz}\vline_{k+1/2}\left(z_{k} - z_{k + 1}\right)\Bigg[-\alpha_{k}\left(M_{k - 1/2} - M_{k + 1/2}\right)-\beta_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &-\alpha_{k + 1}\left(M_{k + 1/2} - M_{k + 3/2}\right)-\beta_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &= \frac{\left(w_{k + 1/2}^{(\text{expl.})} - w_{k + 1/2}^{(n + 1)}\right)\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} - \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vline_{k+1/2}\nonumber\\ &\Leftrightarrow\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\Bigg[\gamma_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &- \gamma_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &+ \frac{d\newoverline{\Pi}}{2dz}\vline_{k+1/2}\left(z_{k} - z_{k + 1}\right)\Bigg[\alpha_{k}\left(M_{k - 1/2} - M_{k + 1/2}\right)+\beta_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &+\alpha_{k + 1}\left(M_{k + 1/2} - M_{k + 3/2}\right)+\beta_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg] - \frac{w_{k + 1/2}^{(n + 1)}\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}}\nonumber\\ &= -\frac{w_{k + 1/2}^{(\text{expl.})}\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} + \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} + \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vline_{k + 1/2}.\tag{34.115}\label{eq:vert_solver_pre_1} \end{align} \]
$w_{k + 1/2}^{(n + 1)}$ also depends on $M_{k + 1/2}$. From Eq. (34.110) follows
\[ \begin{align} w^{(n + 1)}_{k + 1/2} &= \frac{\frac{2M_{k + 1/2}}{A_{k + 1/2}} - \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}w^{(n)}_{k + 1/2}}{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}.\tag{34.116}\label{eq:w_from_ml} \end{align} \]
Putting this into Eq. (34.115), you get
\[ \begin{align} &\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\Bigg[\gamma_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &- \gamma_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &+ \frac{d\newoverline{\Pi}}{2dz}\vline_{k+1/2}\left(z_{k} - z_{k + 1}\right)\Bigg[\alpha_{k}\left(M_{k - 1/2} - M_{k + 1/2}\right)+\beta_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &+\alpha_{k + 1}\left(M_{k + 1/2} - M_{k + 3/2}\right)+\beta_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &- \frac{\frac{2M_{k + 1/2}}{A_{k + 1/2}} - \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}w^{(n)}_{k + 1/2}}{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}\frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}}\nonumber\\ &= -\frac{w_{k + 1/2}^{(\text{expl.})}\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} + \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} + \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vline_{k + 1/2}. \end{align} \]
For the density $ \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}$ applies
\[ \begin{align} \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)} &= \frac{1}{2}\rho_{k}^{(n + 1)} + \frac{1}{2}\rho_{k + 1}^{(n + 1)}\nonumber\\ &= \frac{1}{2}\rho_{k}^{(\text{expl.})} - \Delta t\frac{M_{k - 1/2} - M_{k + 1/2}}{2V_{k}} + \frac{1}{2}\rho_{k + 1}^{(\text{expl.})} - \Delta t\frac{M_{k + 1/2} - M_{k + 3/2}}{2V_{k + 1}}\nonumber\\ &= \newoverline{\rho_{k}^{(\text{expl.})}}^{(k + 1/2)} - \Delta t\frac{M_{k - 1/2} - M_{k + 1/2}}{2V_{k}} - \Delta t\frac{M_{k + 1/2} - M_{k + 3/2}}{2V_{k + 1}}. \end{align} \]
From this it follows
\[ \begin{align} &\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\Bigg[\gamma_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &- \gamma_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &+ \frac{d\newoverline{\Pi}}{2dz}\vline_{k + 1/2}\left(z_{k} - z_{k + 1}\right)\Bigg[\alpha_{k}\left(M_{k - 1/2} - M_{k + 1/2}\right) + \beta_{k}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}\right)\nonumber\\ &+\alpha_{k + 1}\left(M_{k + 1/2} - M_{k + 3/2}\right)+\beta_{k + 1}\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}M_{k + 3/2}\right)\Bigg]\nonumber\\ &- \frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}}\frac{\frac{2M_{k + 1/2}}{A_{k + 1/2}} - w^{(n)}_{k + 1/2}\newoverline{\rho_k^{(\text{expl.})}}^{(k + 1/2)} + w^{(n)}_{k + 1/2}\Delta t\frac{M_{k - 1/2} - M_{k + 1/2}}{2V_{k}} + w^{(n)}_{k + 1/2}\Delta t\frac{M_{k + 1/2} - M_{k + 3/2}}{2V_{k + 1}}}{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}\nonumber\\ &= -\frac{w_{k + 1/2}^{(\text{expl.})}\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} + \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} + \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vline_{k + 1/2}. \end{align} \]
This is the equation to solve. This is of the form
\[ \begin{align} c_{k - 1/2}M_{k - 1/2} + d_{k + 1/2}M_{k + 1/2} + e_{k + 3/2}M_{k + 3/2} = r_{k + 1/2}. \end{align} \]
applies to the coefficients
\[ \begin{align} c_{k - 1/2} &= \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\gamma_{k}\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)} + \frac{d\newoverline{\Pi}}{2dz}\vline_{k + 1/2}\left(z_{k} - z_{k + 1}\right)\left(\alpha_{k} + \beta_{k}\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}\right)\nonumber\\ &- \frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta tc^{(p)}}\frac{w^{(n)}_{k + 1/2}}{2V_{k}\newoverline{\rho_k^{(n)}}^{(k + 1/2)}},\\ d_{k + 1/2} &= -\left(\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\right)^2\left(\gamma_{k} + \gamma_{k + 1}\right) + \frac{d\newoverline{\Pi}}{2dz}\vline_{k + 1/2}\left(z_{k} - z_{k + 1}\right)\left[\alpha_{k + 1} - \alpha_{k} + \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\left(\beta_{k + 1} - \beta_{k}\right)\right]\nonumber\\ &- \frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}\left[\frac{2}{A_{k + 1/2}} + \frac{\Delta tw^{(n)}_{k + 1/2}}{2}\left(-\frac{1}{V_{k}} + \frac{1}{V_{k + 1}}\right)\right],\\ e_{k + 3/2} &= \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\gamma_{k + 1}\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)} - \frac{d\newoverline{\Pi}}{2dz}\vline_{k + 1/2}\left(z_{k} - z_{k + 1}\right)\left(\alpha_{k + 1} + \beta_{k + 1}\newoverline{\theta_k^{(n + 1)^\star}}^{(k + 3/2)}\right)\nonumber\\ &+ \frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta tc^{(p)}}\frac{w^{(n)}_{k + 1/2}}{2V_{k + 1}\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}, \end{align} \] \[ \begin{align} r_{k + 1/2} &= -\frac{w_{k + 1/2}^{(\text{expl.})}\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}} + \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}\frac{\Pi^{'(\text{expl.})}_{k} - \Pi^{'(\text{expl.})}_{k + 1}}{\Delta t} + \frac{z_{k} - z_{k + 1}}{\Delta t}\left(\theta_{k}^{'(\text{expl.})}+\theta_{k + 1}^{'(\text{expl.})}\right)\frac{d\newoverline{\Pi}}{2dz}\vline_{k + 1/2}\nonumber\\ &- \frac{\left(z_{k} - z_{k + 1}\right)}{\epsilon\Delta t^2c^{(p)}}\frac{w^{(n)}_{k + 1/2}\newoverline{\rho_k^{(\text{expl.})}}^{(k + 1/2)}}{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}. \end{align} \]
Linear systems of equations with three-band matrices can be solved using the Thomas algorithm. From the result for $M_{k + 1/2}$ the other unknowns can be determined as follows:
\[ \begin{align} \rho^{(n + 1)}_k &= \rho^{(\text{expl.})}_k - \Delta t\frac{M_{k - 1/2} - M_{k + 1/2}}{V_k},\\ \newtilde{\theta}^{(n + 1)}_k &= \newtilde{\theta}^{(\text{expl.})}_k - \Delta t\frac{\newoverline{\theta_k^{(n + 1)^\star}}^{(k - 1/2)}M_{k - 1/2} - \newoverline{\theta_k^{(n + 1)^\star}}^{(k + 1/2)}M_{k + 1/2}}{V_k},\\ \Pi_k^{'(n+1)} &= \Pi_k^{'(\text{expl.})} + \gamma_kV_k\left(\newtilde{\theta}_k^{(n + 1)} - \newtilde{\theta}_k^{(\text{expl.})}\right),\\ w^{(n + 1)}_{k + 1/2} &\stackrel{\href{#eq:w_from_ml}{\text{Glg. (34.116)}}}{=} \frac{\frac{2M_{k + 1/2}}{A_{k + 1/2}} - \newoverline{\rho_k^{(n + 1)}}^{(k + 1/2)}w^{(n)}_{k + 1/2}}{\newoverline{\rho_k^{(n)}}^{(k + 1/2)}}. \end{align} \]
If you don't want to use a hydrostatic background state, you just have to use the replacements
\[ \begin{align} \newoverline{\theta}_k &= 0,\\ \newoverline{\Pi}_k &= 0,\\ w_{k + 1/2}^{(\text{expl.})} &\to w_{k + 1/2}^{(\text{expl.})} - \Delta t\frac{d\newoverline{\phi}}{dz}\vline_{k + 1/2} \end{align} \]
make.
The forward-backward procedure was presented in Sect. 25.4.1. It does not increase the accuracy, but it does increase the stability of the model. In the dynamic kernel to be formulated here, the horizontal momentum equation is first integrated forward, then the scalar equations are solved with the new values of the velocities.
The performance of a HEVI process with separation of horizontal fast (acute) and slow (advective) components can be generally calculated. Let $\tau$ be the acoustic time step and $T$ be the advective time step. We have $T\geq\tau$. However, when processing a complete time step, higher-order methods calculate the corresponding terms several times. $n_\tau$ is the number of calculations of the acoustic components and $n_T$ is the number of calculations of the advective components. This allows effective time steps
\[ \begin{align} \tau_\text{eff} &\coloneqq \frac{\tau}{n_\tau},\\ T_\text{eff} &\coloneqq \frac{T}{n_T} \end{align} \]
define. The larger these time steps are, the better for the performance of the model. The speed $c$ of the model, i.e. the integrated time $t$ divided by the time required for it, is calculated by
\[ \begin{align} c = \frac{t}{\frac{t}{T_\text{eff}}t_T + \frac{t}{\tau_\text{eff}}t_\tau} = \frac{1}{\frac{1}{T_\text{eff}}t_T + \frac{1}{\tau_\text{eff}}t_\tau} = \frac{\tau_\text{eff}T_\text{eff}}{\tau_\text{eff}t_T + T_\text{eff}t_\tau}. \end{align} \]
Here $t_T$ is the time required for an advective time step and $t_\tau$ is the time required for an acoustic time step. With the notation $t_T\equiv\alpha t_\tau$ this becomes
\[ \begin{align} c = \frac{\tau_\text{eff}T_\text{eff}}{\tau_\text{eff}\alpha t_\tau + T_\text{eff}t_\tau} = \frac{\tau_\text{eff}}{t_\tau}\frac{T_\text{eff}}{\alpha\tau_\text{eff} + T_\text{eff}}. \end{align} \]
The basic time-stepping procedure of the model should be explicit. The explicit Euler method is only first order and therefore not precise enough. A multi-step process produces unphysical modes, which is also disadvantageous. Therefore, a predictor-corrector method is used.
The basic architecture of the time step method is that of the predictor-corrector method. In addition, some modifications are made to stabilize and increase efficiency. The basic procedure for each of the three sub-steps is as follows:
The atmosphere is global. A regional model is still practical for conducting experiments and operational runs with increased resolution. The reason for choosing the hexagonal grid for the global model developed in this book is its numerical properties, which are similar to the quadrangular grid, combined with quasi-uniformity and orthogonality. Therefore, one could also implement a local model on a hexagonal grid. The disadvantage, however, is the effort that the grid generation process entails as well as the difficulty of parallelization and the possible additional effort that is necessary to achieve satisfactory performance. Therefore, in this section, the previous findings in this chapter will be generalized to a regional longitude-latitude grid, since the pole problem does not occur in a regional model anyway.
The term nesting describes the generation of boundary conditions for a regional model. Be
\[ \begin{align} \psi = \psi\left(\mathbf{x}, z, t\right) \end{align} \]
a model state, where $\mathbf{x}$ denotes the horizontal coordinates. A state is now generated from a model with a larger domain, which is referred to here as an edge model, for example a global model
\[ \begin{align} \phi = \phi\left(\mathbf{x}, z, t\right), \end{align} \]
which is obtained by spatial and temporal interpolation of the boundary model onto the grid of the target model. In this case, the time is usually linearized linearly between the time steps of the boundary model, i.e
\[ \begin{align} \phi\left(\mathbf{x}, z, t\right) = \frac{t_{n+1} - t}{t_{n+1} - t_n}\phi\left(\mathbf{x}, z, t_n\right) + \frac{t - t_n}{t_{n+1} - t_n}\phi\left(\mathbf{x}, z, t_{n + 1}\right). \end{align} \]
At each time step of the target model, the model state $\psi$ is moved a little towards $\phi$ at the edges of the domain, i.e
\[ \begin{align} \psi_\text{neu}\left(\mathbf{x}, z, t\right) = \left(1 - f\left(\mathbf{x}\right)\right)\psi_\text{alt}\left(\mathbf{x}, z, t\right) + f\left(\mathbf{x}\right)\phi\left(\mathbf{x}, z, t\right). \end{align} \]
Here $f = f\left(\mathbf{x}\right)$ is a function that describes the admixture of the boundary state. The simplest approach would be to set this function to one at the very edge of the domain and zero otherwise. However, this can create relatively sharp gradients, which can emit unphysical waves into the model area. For example, one approach makes more sense
\[ \begin{align} f = f\left[d\left(\mathbf{x}\right)\right] = \cos\left[\text{min}\left(\frac{d\pi}{2n\Delta x},\frac{\pi}{2}\right)\right]^2. \end{align} \]
Here, $d\left(\mathbf{x}\right)$ is the distance of the location $\mathbf{x}$ from the edge of the domain and $n$ is the number of grid cells within which the weight of the boundary conditions drops to zero. These procedures are also known as Nudging.