In this chapter, the reversible discretization developed in the previous chapter is extended to include irreversible effects. It is still assumed that the air is dry. The irreversible effects in a dry atmosphere are the diffusions. There are two types of diffusion:
Temperature diffusion
Momentum diffusion
On the molecular scale, according to Eq. (5.212), there is also diffusion of mass; however, turbulence in a dry atmosphere does not diffuse any mass. This is because turbulence diffuses only the mixing ratio, which is unity in a dry atmosphere, so there is no spatial variability that could be diffused.
Accordingly, the system of equations Equations (34.4) - (34.7) must be modified as follows:
\[ \begin{align} \frac{\partial\mathbf{v}}{\partial t} &= -c^{(p)}\theta\nabla\Pi - \nabla k + \mathbf{v}\times\etabi - \nabla\phi + \textcolor{red}{\mathbf{f}_R}\tag{35.1}\label{eq:prog_irrev_0}\\ \frac{\partial\rho}{\partial t} &= -\nabla\cdot\left(\rho\mathbf{v}\right)\tag{35.2}\label{eq:prog_irrev_1}\\ \frac{\partial\newtilde{\theta}}{\partial t} &= -\nabla\cdot\left(\newtilde{\theta}\mathbf{v}\right) + \textcolor{red}{\frac{q^{(V)}}{c^{(p)}\Pi}}\tag{35.3}\label{eq:prog_irrev_2}\\ \Pi &= f\left(\newtilde{\theta}\right)\tag{35.4}\label{eq:prog_irrev_3} \end{align} \]
The mass diffusion is included in the continuity equation in the form of the convergence of a flux density $\nabla\cdot\left(\kappa_\rho\nabla\rho\right)$, where $\kappa_\rho$ is the associated molecular diffusion coefficient. Two effects are included in the thermal power density $q^{(V)}$:
\[ \begin{align} q^{(V)} = q^{(V)}_\text{diss} + q^{(V)}_\text{Tdiff} = \rho\epsilon + \nabla\cdot\left(c_d^{(v)}\rho\kappa_T\nabla T\right) \end{align} \]
Here $\epsilon$ is the dissipation. The term $\nabla\cdot\left(c_d^{(v)}\rho\kappa_T\nabla T\right)$ describes the convergence of the heat flux density resulting from the temperature diffusion.
The temperature equation Eq. (35.4) is the same as Eq. (34.7), the diffusive effects also enter into the temperature via $\Delta\rho$ and $\Delta\newtilde{s}$. In analogy to the previous chapter, the horizontal and vertical discretizations are also derived separately in this chapter, even if there is no longer any distinction of the vertical spatial direction in a non-hydrostatic model.
All terms developed in this chapter are only calculated in the Predictor step and kept constant in the Corrector step.
The diffusive (irreversible) terms represent a coupling of the continuous world of the Navier-Stokes equations to the molecular scale. As stated in Sect. 8.2.4, momentum diffusion destroys (dissipates) kinetic energy in favor of internal energy. The atmosphere can be conceptualized as a low-pass filter, where wavelengths $< l_c$ are filtered out; here $l_c$ is the characteristic length according to Eq. (17.159).
So far only actual physical equations have been discretized, this is called direct numerical simulation (DNS). However, this approach is problematic if, for the resolution $\Delta$,
\[ \begin{align} 2\Delta\gg l_c, \end{align} \]
holds, since such a discretization is a low-pass filter with cutoff wavelength $2\Delta$. Accordingly, in a numerical model, the dissipation must begin earlier. According to Eq. (17.159), the following equation holds for the kinematic viscosity $\nu$:
\[ \begin{align} \nu = \left(l_c^4\epsilon\right)^{1/3}. \end{align} \]
Accordingly, one expects that
\[ \begin{align} \frac{\nu_\text{numerisch}}{\nu_\text{molekular}} \sim \left(\frac{\text{min}\left(l_m, 2\Delta r\right)}{l_c}\right)^{4/3} \end{align} \]
is a useful estimate for the numerical viscosity; here $l_m$ is the mixing path length. The restriction to the mixing path length in the numerator makes sense because the fluid in the model must never be more viscous than the natural medium, otherwise no realistic simulations of the dynamics can be expected. For the vertical mixing path length, a value of approx. 20 m within the Ekman layer was estimated using Eq. (17.75); with this and with the estimate $l_c = 5$ mm, one obtains
\[ \begin{align} \text{log}_{10}\left(\frac{\nu_\text{numerisch}}{\nu_\text{molekular}}\right) = 4,8. \end{align} \]
This roughly corresponds to the value obtained in Eq. (17.74).
Outside the planetary boundary layer, a value for the vertical mixing path length of approx. 100 m can be assumed, as larger eddies exist there. The horizontal mixing path length can be estimated to be approximately 10 km outside the planetary boundary layer (two orders of magnitude greater than the vertical mixing path length). Simulations whose resolution is less than the mixing path length are referred to as large eddy simulations (LES). As the resolution is further increased, the numerical and molecular diffusion coefficients converge.
The classical Smagorinsky model is used as the ansatz for momentum diffusion. For this, Eq. (17.172) must now be discretized on the hexagonal grid. To do this, one first averages the discretized terms:
\[ \begin{align} K_\Delta = \rho c_S^2\Delta^2\sqrt{\newoverline{\left(\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y}\right)}^2 + \newoverline{\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)}^2} = \rho c_S^2\Delta^2\sqrt{\left(\newoverline{\frac{\partial u}{\partial x}} - \newoverline{\frac{\partial v}{\partial y}}\right)^2 + \left(\newoverline{\frac{\partial u}{\partial y}} + \newoverline{\frac{\partial v}{\partial x}}\right)^2} \end{align} \]
To determine the averaged terms, one first picks out, as an example, the term $\frac{\partial u}{\partial x}$ and computes
\[ \begin{align} \newoverline{\frac{\partial u}{\partial x}} = \frac{1}{A_c}\int_{A_c}\frac{\partial u}{\partial x}d^2r. \end{align} \]
This integral is now evaluated in Cartesian coordinates on the tangential plane:
\[ \begin{align} \newoverline{\frac{\partial u}{\partial x}} = \frac{1}{A_c}\int_{y_1}^{y_2}\int_{x_1\left(y\right)}^{x_2\left(y\right)}\frac{\partial u}{\partial x}dxdy = \frac{1}{A_c}\int_{y_1}^{y_2}u\left(x_2,y\right) - u\left(x_1,y\right)dy. \end{align} \]
This can be reformulated as a line integral:
\[ \begin{align} \newoverline{\frac{\partial u}{\partial x}} = \frac{1}{A_c}\int_{\partial A_c}u\left(x,y\right)\mathbf{e}_y\cdot d\mathbf{r} = \frac{1}{A_c}\int_{\partial A_c}u\left(x,y\right)\cos\left(\phi\right)ds, \end{align} \]
here $\phi$ is the angle between $d\mathbf{r}$ and $\mathbf{e}_y$.
From physical derivations one usually obtains diffusive terms proportional to the second derivative of the field to be diffused. However, for numerical reasons, so-called Hyperdiffusion is often used in models, which is proportional to a higher power of the Nabla operator. In the one-dimensional case, the ansatz for fourth-order diffusion reads
\[ \begin{align} \frac{\partial\psi}{\partial t} = -\alpha\frac{\partial^4\psi}{\partial x^4}, \end{align} \]
here $\alpha$ is often referred to as the hyperviscosity. Making the ansatz $\psi = \psi_0\exp\left(-ikx-\frac{t}{\tau}\right)$ here, it follows
\[ \begin{align} -\frac{1}{\tau}\psi &= -\alpha k^4\psi\nonumber\\ \Leftrightarrow \frac{1}{\tau} &= \alpha k^4. \end{align} \]
The shortest resolved wave has length $2\Delta x$, from which $k = \frac{2\pi}{2\Delta x} = \frac{\pi}{\Delta x}$ follows. Substituting this, one obtains
\[ \begin{align} \frac{1}{\tau} &= \alpha\frac{\pi^4}{\Delta x^4}\nonumber\\ \Leftrightarrow \alpha &= \frac{A^2}{\pi^4\tau} \end{align} \]
for the relationship between the diffusion coefficient and the decay time of the shortest resolved wave, where $A = \Delta x^2$ is the area of the grid cells. Expressing $\tau$ as a multiple of the model time step $\Delta t$, one obtains
\[ \begin{align} \alpha &= \frac{A^2}{c\Delta t\pi^4}. \end{align} \]
The subscale viscosity describes the influence of the TKE on the resolved scales. One can introduce the volumetric TKE
\[ \begin{align} K_s \coloneqq \rho k_s \end{align} \]
as an additional prognostic quantity. For this quantity, Eq. (17.189) holds:
\[ \begin{align} \frac{\partial K_s}{\partial t} = -\nabla\cdot\left(K_s\mathbf{v}\right) - \rho\newoverline{u'w'}\frac{\partial\newoverline{u}}{\partial z} - \rho\newoverline{v'w'}\frac{\partial\newoverline{v}}{\partial z} + \rho g\newoverline{\frac{\theta'}{\newoverline{\theta}}w'} - \rho\newoverline{\epsilon}\tag{35.19}\label{eq:tke_num} \end{align} \]
To calculate the covariance terms of the velocity components, the K-theory approach is used:
\[ \begin{align} \newoverline{u'w'} &= -K\frac{\partial\newoverline{u}}{\partial z},\\ \newoverline{v'w'} &= -K\frac{\partial\newoverline{v}}{\partial z},\\ g\newoverline{\frac{\theta'}{\newoverline{\theta}}w'} &= -K\frac{g}{\newoverline{\theta}}\frac{\partial\newoverline{\theta}}{\partial z} = -KN^2. \end{align} \]
Substituting this into Eq. (35.19), one obtains
\[ \begin{align} \frac{\partial K_s}{\partial t} = -\nabla\cdot\left(K_s\mathbf{v}\right) + \rho K\left[\left(\frac{\partial\newoverline{u}}{\partial z}\right)^2 + \left(\frac{\partial\newoverline{v}}{\partial z}\right)^2\right] - \rho KN^2 - \rho\newoverline{\epsilon}.\tag{35.23}\label{eq:tke_num_mod} \end{align} \]
To determine $\newoverline{\epsilon}$, one starts from Eq. (17.159):
\[ \begin{align} l_c = \left(\frac{\nu^3}{\newoverline{\epsilon}}\right)^{1/4} \Leftrightarrow l_c^4 = \frac{\nu^3}{\newoverline{\epsilon}} \Leftrightarrow \newoverline{\epsilon} = \frac{\nu^3}{l_c^4}. \end{align} \]
With the substitutions $\nu\to K$ and $l_c\to L$ one obtains
\[ \begin{align} \newoverline{\epsilon} = \frac{K^3}{L^4}.\tag{35.25}\label{eq:diff_delta} \end{align} \]
The viscosity coefficients can be expressed as functions
\[ \begin{align} K = K\left[\mathbf{v}\left(\mathbf{r}\right),K_s\right] \end{align} \]
of the current wind field and of $K_s$. Through Eq. (35.23), this ansatz also contains information about the past wind field.
From Eq. (35.23), neglecting advection, one obtains an approximate diagnostic relation for the dissipation of the TKE:
\[ \begin{align} \newoverline{\epsilon} = K\left[\left(\frac{\partial\newoverline{u}}{\partial z}\right)^2 + \left(\frac{\partial\newoverline{v}}{\partial z}\right)^2 - N^2\right]\tag{35.27}\label{eq:tke_diag} \end{align} \]
The TKE approach is ultimately about calculating the vertical diffusion coefficient. There are basically three approaches to this:
In first-order closure, the vertical diffusion coefficient is calculated using the vertical wind shear multiplied by a length scale.
In the second-order closure, prognostic variables are solved for both the TKE and the dissipation of the TKE. The vertical diffusion coefficient can be derived from this.
With the 1.5th order closure a prognostic equation is solved for the TKE, but not for the dissipation. The vertical diffusion coefficient is derived algebraically from the TKE.
Only the 1.5th order closure is pursued here. To determine the vertical diffusion coefficient, the kinetic gas model according to Sect. 5.4.0.2 is used as a model. In this case, particles fly a mean free path between two collisions, which is here called mixing path length $L_m$ , their average speed is $\newoverline{v}$. The resulting diffusion coefficient is given, according to Eq. (5.215), by
\[ \begin{align} K_v = L_m\newoverline{v}.\tag{35.28}\label{eq:k_eff_deriv_1} \end{align} \]
The factor $\frac{1}{3}$ from Eq. (5.215) does not occur here because only vertical movements are considered here. It is now further assumed that the fluid particles oscillate vertically with the angular frequency $N$ according to Eq. (16.256) and thereby have a maximum amplitude $\newhat{z}$. Their trajectory is then
\[ \begin{align} z\left(t\right) = \newhat{z}\exp\left(iNt\right). \end{align} \]
From this it follows
\[ \begin{align} \newoverline{v} = \frac{4\newhat{z}}{T} = \frac{4\newhat{z}}{2\pi}\frac{2\pi}{T} = \frac{2\newhat{z}N}{\pi},\tag{35.30}\label{eq:k_v_deriv_1} \end{align} \]
since the particle covers the distance $\newhat{z}$ four times during each period. For the maximum speed, one has
\[ \begin{align} v_\text{max} = N\newhat{z} = \frac{\pi}{2}\newoverline{v}. \end{align} \]
The mean free path $L_m$ here is $2\newhat{z}$. Now $\newhat{z}$ must be expressed in terms of the TKE $k_s$. For the specific total energy $E$ of the harmonic oscillator, the relation
\[ \begin{align} E = \frac{1}{2}N^2\newhat{z}^2. \end{align} \]
holds according to Eq. (2.89). On average, half of this is stored as kinetic energy, which leads to
\[ \begin{align} k_s = \frac{1}{4}N^2\newhat{z}^2 \Rightarrow \newhat{z} = \frac{2}{N}\sqrt{k_s} \Rightarrow L_m = \frac{4}{N}\sqrt{k_s} = \frac{C_l}{N}\sqrt{k_s} \end{align} \]
According to the derivation, $C_l = 4$, but empirically one rather finds $C_l = 1$. Substituting this into Eq. (35.30), one obtains
\[ \begin{align} \newoverline{v} = \frac{2\newhat{z}N}{\pi} = \frac{4}{\pi}\sqrt{k_s}.\tag{35.34}\label{eq:tke_vbar} \end{align} \]
Assuming that a third of the particles fly along the direction of one of the coordinate axes, it follows for the diffusion coefficient with Eq. (35.28)
Empirically, one finds that the prefactor $C \coloneqq \frac{4}{3\pi}$ is not the same for all diffused quantities. This is because the correlations between the fluctuations of different quantities differ. At this point,
\[ \begin{align} C_m &= 0,126,\\ C_e &= 0,34,\\ C_h &= 0,142 \end{align} \]
is chosen [39]. Here, $C_m$ is the prefactor for the vertical diffusion of momentum, $C_e$ is the prefactor for the vertical diffusion of TKE and $C_h$ is the prefactor for the vertical diffusion of humidity and potential temperature (passive tracers).
From Eq. (35.35), using Eq. (35.25), a formula for the dissipation of the TKE can be derived:
\[ \begin{align} \newoverline{\epsilon} = \frac{\left(C_m L_m\sqrt{k_s}\right)^3}{L_m^4} = C_m^3\frac{k_s^{3/2}}{L_m} = C_\epsilon\frac{k_s^{3/2}}{L_m}. \end{align} \]
with $C_\epsilon\coloneqq C_m^3$. Substituting this into Eq. (35.27), one obtains
\[ \begin{align} C_\epsilon\frac{k_s^{3/2}}{L_m} &= K\left[\left(\frac{\partial\newoverline{u}}{\partial z}\right)^2 + \left(\frac{\partial\newoverline{v}}{\partial z}\right)^2 - N^2\right]\nonumber\\ \Leftrightarrow C_\epsilon\frac{k_s^{3/2}}{L_m} &= L_mC_m\sqrt{k_s}\left[\left(\frac{\partial\newoverline{u}}{\partial z}\right)^2 + \left(\frac{\partial\newoverline{v}}{\partial z}\right)^2 - N^2\right]\nonumber\\ \Leftrightarrow k_s &= \frac{C_mL_m^2}{C_\epsilon}\left[\left(\frac{\partial\newoverline{u}}{\partial z}\right)^2 + \left(\frac{\partial\newoverline{v}}{\partial z}\right)^2 - N^2\right]. \end{align} \]
The basic problem of convection parameterizations has already been outlined in Sect. 17.9.4. The approach to bringing these considerations into the model is that the convection has to be represented in two parts:
The resolved fraction is the one that is „seen“ by a non-hydrostatic model in the prognostic variables. This fraction does not need to be treated further, but is simulated explicitly.
The unresolved portion consists of unstable sub-regions (area fraction $0\leq c\leq 1$), which are smaller than the grid cells. This proportion is not taken into account by the model without further action.
Convection is vertical mixing. An attempt is therefore made to represent the unresolved portion by modifying the mixing path length:
\[ \begin{align} L_m \to cL_\text{conv} + \left(c - 1\right)L_m, \end{align} \]
here $L_\text{conv}$ is a convective mixing path length. $c$ is the fraction of the base area of the grid cell in which thermal instability prevails:
\[ \begin{align} c = \frac{1}{A_c}\int_{A_c}\Theta\left(\Gamma-\Gamma_d\right)dA\tag{35.42}\label{eq:ansatz_c_conv} \end{align} \]
Let $\newoverline{\Gamma}$ be the mean negative vertical temperature gradient, i.e. the one calculated from the averaged (i.e. resolved) variables. If one further assumes that $\Gamma$ is normally distributed, it follows from Eq. (35.42)
\[ \begin{align} c &= \int_{\Gamma_d}^\infty f_\sigma\left(\Gamma - \newoverline{\Gamma}\right)d\Gamma = \int_{\Gamma_d - \newoverline{\Gamma}}^\infty f_\sigma\left(\Gamma\right)d\Gamma = \int_{\frac{\Gamma_d - \newoverline{\Gamma}}{\sigma}}^\infty f_1\left(\Gamma\right)d\Gamma = \int_{\frac{\Gamma_d - \newoverline{\Gamma}}{\sigma}}^\infty\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{\Gamma^2}{2}\right)d\Gamma\nonumber\\ &= \frac{1}{\sqrt{\pi}}\int_{\frac{\Gamma_d - \newoverline{\Gamma}}{\sqrt{2}\sigma}}^\infty\exp\left(-\Gamma^2\right)d\Gamma = \frac{1}{\sqrt{\pi}}\left(\int_{0}^\infty\exp\left(-\Gamma^2\right)d\Gamma - \int_{0}^{\frac{\Gamma_d - \newoverline{\Gamma}}{\sqrt{2}\sigma}}\exp\left(-\Gamma^2\right)d\Gamma\right)\nonumber\\ &= \frac{1}{\sqrt{\pi}}\left(\frac{\sqrt{\pi}}{2} - \int_{0}^{\frac{\Gamma_d - \newoverline{\Gamma}}{\sqrt{2}\sigma}}\exp\left(-\Gamma^2\right)d\Gamma\right) = \frac{1}{2} - \frac{1}{\sqrt{\pi}}\int_{0}^{\frac{\Gamma_d - \newoverline{\Gamma}}{\sqrt{2}\sigma}}\exp\left(-\Gamma^2\right)d\Gamma = \frac{1}{2}\left[1 - \text{erf}\left(\frac{\Gamma_d - \newoverline{\Gamma}}{\sqrt{2}\sigma}\right)\right]. \end{align} \]
The standard deviation of the stratification $\sigma$ now needs to be determined. In a random walk, the expected value of the distance from the origin is proportional to the square root of the number of steps. One now proceeds radially from the center of a grid cell towards the edge in steps of $\Delta r$. For the standard deviation $\sigma_r$ of the probability distribution of the stratification $\Gamma_r$ at a point at distance $r$ from the center, one therefore makes the ansatz
\[ \begin{align} \sigma_r = k\sqrt{r}. \end{align} \]
The standard deviation $\newoverline{\sigma}$ averaged over the base area of a circular cell thus satisfies
\[ \begin{align} \newoverline{\sigma} = \frac{\int_A\sigma dA}{\int_A\sigma dA} = \frac{\int_0^R2\pi r\sigma_rdr}{\pi R^2} = \frac{\int_0^R2\pi rk\sqrt{r}dr}{\pi R^2} = \frac{2\pi k\int_0^R r^{3/2}dr}{\pi R^2} = 2\pi k\frac{2R^{5/2}}{5\pi R^2} = 2\pi k\frac{2R^{1/2}}{5\pi} = \frac{4k}{5\pi^{1/4}}A_c^{1/4}. \end{align} \]
For $k$, one makes a heuristic ansatz
\[ \begin{align} k = \frac{\delta\Gamma}{\sqrt{d}}. \end{align} \]
Here $\delta\Gamma$ is a typical fluctuation range, for which $\delta\Gamma = \frac{\Gamma_d}{3}$ is assumed. $d$ is the typical length scale associated with this variation, which is assumed to be 1000 km.
So far, the dry-adiabatic temperature gradient $\Gamma_d$ has been assumed in Eq. (35.42). In reality, however, a subset of a grid box is saturated; the moist-adiabatic temperature gradient $\Gamma_h$ must be used there. In order to take this into account approximately, from Eq. (35.42) onward, the replacement
\[ \begin{align} \Gamma_d \to c_f\Gamma_h + \left(1 - c_f\right)\Gamma_d \end{align} \]
is made. Here $c_f$ is the cloud-filled (i.e. supersaturated) volume fraction according to Sect. 36.4.