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:
On the molecular scale, according to Eq. (5.212) also diffusion of mass, but turbulence in a dry atmosphere does not diffuse mass. This is because turbulence diffuses only the mixing ratio, which is unity in a dry atmosphere, and thus no spatial variability to diffuse.
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 and the molecular scale. As stated in Sect. 8.2.4, momentum diffusion annihilates (dissipates). kinetic energy in favor of internal energy. The atmosphere can be conceptualized as a low pass, where wavelengths $< l_c$ are filtered, where $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 the resolution is $\Delta$
\[ \begin{align} 2\Delta\gg l_c, \end{align} \]
Such a discretization is namely a low-pass with the cutoff wavelength $2\Delta$. Accordingly, in a numerical model, the dissipation must begin earlier. According to Eq. (17.159) applies to the kinematic viscosity $\nu$ the equation
\[ \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 meter makes sense because the fluid in the model must never be tougher (more viscous) than the natural medium, otherwise no realistic simulations of the dynamics can be expected. For the vertical mixing path length, Eq. (17.75) within the Ekman layer a value of approx. 20 m is estimated, with this and with the estimates $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 that in Eq. (17.74) obtained value.
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 an approach for momentum diffusion. For this, Eq. (17.172) must now be discretized on the hexagonal grid. To do this, you first average 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, first take the term $\frac{\partial u}{\partial x}$ as an example and calculate
\[ \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 approach for diffusion is 4th order
\[ \begin{align} \frac{\partial\psi}{\partial t} = -\alpha\frac{\partial^4\psi}{\partial x^4}, \end{align} \]
here like $\alpha$ often referred to as hyperviscosity. If you do the approach $\psi = \psi_0\exp\left(-ikx-\frac{t}{\tau}\right)$, 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 the length $2\Delta x$, from which $k = \frac{2\pi}{2\Delta x} = \frac{\pi}{\Delta x}$. If you put this in, you get
\[ \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 use the volumetric TKE
\[ \begin{align} K_s \coloneqq \rho k_s \end{align} \]
introduce as an additional prognostic variable. For this size, Eq. (17.189):
\[ \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} \]
Putting this into Eq. (35.19), you get
\[ \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) from:
\[ \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 used as functions
\[ \begin{align} K = K\left[\mathbf{v}\left(\mathbf{r}\right),K_s\right] \end{align} \]
of the current wind field and $K_s$ can be expressed. About Eq. (35.23), this approach 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:
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 calculated according to Eq. (5.215) to
\[ \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 move vertically with the angular frequency $N$ according to Eq. (16.256) oscillate and have a maximum amplitude $\newhat{z}$. Your 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 applies
\[ \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 by the TKE $k_s$. According to Eq., the total specific energy $E$ of the harmonic oscillator is: (2.89) the relation
\[ \begin{align} E = \frac{1}{2}N^2\newhat{z}^2. \end{align} \]
On average, half of this is stored in kinetic energy, which is
\[ \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} \]
leads. According to the derivation, $C_l = 4$, but empirically we find that $C_l = 1$. Putting this into Eq. (35.30), you get
\[ \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)
\[ \begin{align} K_v = \frac{L_m}{3}\newoverline{v} = \frac{4}{3\pi}L_m\sqrt{k_s} \approx 0,424L_m\sqrt{k_s}.\tag{35.35}\label{eq:tke_diff_derived} \end{align} \]
Empirically, one finds that the prefactor $C \coloneqq \frac{4}{3\pi}$ is not the same size for all diffused quantities. This is because the correlations between the fluctuations of different sizes are different. At this point it will
\[ \begin{align} C_m &= 0,126,\\ C_e &= 0,34,\\ C_h &= 0,142 \end{align} \]
selected [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) can be derived using Eq. (35.25) derive a formula for the dissipation of the TKE
\[ \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$. Putting this into Eq. (35.27), you get
\[ \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:
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 proportion 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 we further assume 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. Now go 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 with the distance $r$ from the center, one therefore makes the approach
\[ \begin{align} \sigma_r = k\sqrt{r}. \end{align} \]
The standard deviation averaged over the base area of a circular cell is $\newoverline{\sigma}$
\[ \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 approach
\[ \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, in Eq. (35.42) assumed the dry adiabatic temperature gradient $\Gamma_d$. 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) the replacement
\[ \begin{align} \Gamma_d \to c_f\Gamma_h + \left(1 - c_f\right)\Gamma_d \end{align} \]
made. Here $c_f$ is the cloud-filled (i.e. supersaturated) volume fraction according to App. 36.4.