ICON-LEM parameterizations

The ICON-LEM configuration is based on the ICON-NWP branch and offers different parameterizations for SGS processes. To get some insight on the Turbulence Closures and Cloud Cover Schemes, this page will give some explanations on what the schemes do, how to steer them, and where you can find the regarding code. The newly implemented Deardorff turbulence (and in the future also the cloud cover) scheme can be found in the ICON-AES repository.

Due to filtering of the governing equations for momentum, heat and moisture covariance terms appear. Since these terms contain unknown information through the coupling to the subgrid-scale,they can not be solved explicitly and therefore have to be parameterized. In the literature, this problem is associated with different terms like the closure problem or as subgrid-scale modeling. The most common approach for Large-eddy Modeling (LEM) of the planetary boundary layer (PBL) is to link the energy transport by SGS motions, e.g. the flux of momentum $\overline{u_{i}'u_{j}'}$ or flux of some scalar $\overline{u_i'\psi'}$, to the mean flow gradient $$\begin{align} \overline{u_{i}'u_{j}'} &= -\,K_\mathrm{m}\,\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \text{ ,} \label{eq:momentumflux} \tag{1} \\ \overline{u_i'\psi'} &= -\,K_\mathrm{h}\,\frac{\partial \psi}{\partial x_i} \text{ .} \label{eq:scalarflux} \tag{2} \end{align} $$ In this, $K_\mathrm{m}$ and $K_\mathrm{h}$ are the (local) SGS eddy diffusivities for momentum and heat, $\psi$ can be the temperature, moisture, liquid water content or any passive scalar. The diffusivities have a direct relationship to the SGS-TKE ($e = \frac{1}{2}\,\overline{u_i'u_i'}$) and therefore can be derived from it. In the following two subsections the already implemented Smagorinsky model (modified version by Lilly [1962]) will be explained as well as the newly added 1.5-order turbulence closure following Deardorff (1980). Both schemes base on the prognostic equation for the SGS-TKE $e$ that reads $$\begin{equation} \frac{\partial e}{\partial t} = \frac{1}{\rho}\frac{\partial(\rho\,\overline{u}_i\, e)}{\partial x_j} - \overline{(u_i' u_j')}\frac{\partial \overline{u}_i}{\partial x_j} + \frac{g}{\theta_\mathrm{v,0}} \overline{u_3'\theta'_\mathrm{v}} - \frac{1}{\rho} \frac{\partial}{\partial x_i} \left(\overline{\rho\,u'_i(e'+\frac{p'}{\rho})}\right) - \epsilon \text{ .} \label{eq:fullTKE} \tag{3} \end{equation}$$ From left to right the terms are the time-rate of change, the advection of TKE, the production of TKE by shear, by buoyancy, the diffusive transport of TKE and the diffusion $\epsilon$. How the different closure schemes are linked to equation \ref{eq:fullTKE} will be further described in the following subsections.

The Smagorinsky approach is based on a local equilibrium constraint that implies a balance between the production of TKE (through shear and buoyancy) and the dissipation of it. With this assumption, the prognostic equation \ref{eq:fullTKE} simplifies to $$\begin{equation} \overline{(u_i' u_j')}\frac{\partial \overline{u}_i}{\partial x_j} + \frac{g}{\theta_\mathrm{v,0}}\overline{u_3'\theta'_\mathrm{v}} = \epsilon \text{ .} \label{eq:smago} \tag{4} \end{equation} $$ One way to solve equation \ref{eq:smago} is the substitution of the covariance terms $$\begin{align} \overline{(u_i' u_j')}\,\frac{\partial \overline{u}_i}{\partial x_j} &= -\,K_\mathrm{m}\,S^2 \text{ ,} \label{eq:shearterm} \tag{10} \\ \frac{g}{\theta_\mathrm{v,0}}\,\overline{u_3'\theta'_\mathrm{v}} &= -\,K_\mathrm{h}\,N^2 \text{ ,} \label{eq:buoyterm} \tag{11} \end{align} $$ with the deformation $S^2$ and the $N^2$ squared Brunt-Vaisala frequency. Additional information about these terms can be found in the technical implementation in chapter 3. Using the parameterized subgrid-scale stress tensor following Lilly (1962) the subgrid diffusivities for momentum and heat can be written as $$\begin{align} K_\mathrm{m} &= 2\,\lambda^2\rho\,||S||\,\left(1 - \frac{\mathrm{Ri}}{\mathrm{Pr}_\mathrm{t}} \right)^{1/2} \label{eq:smagodiffmom} \tag{5} \qquad \text{for} \qquad \frac{\mathrm{Ri}}{\mathrm{Pr}_\mathrm{t}} > 0 \text{ ,} \\ K_\mathrm{h} &= \frac{K_\mathrm{m}}{\mathrm{Pr}_\mathrm{t}} \text{ .} \label{eq:smagodiffheat} \tag{6} \end{align} $$ Here, $\lambda$ is the subgrid length scale for the Smagorinksy approach, which is proportional to the grid volume amongst other factors. $\rho$ describes the density, $||S|| = \sqrt{S^2} = (S_\mathrm{ij}S_\mathrm{ij})^{1/2}$ is the deformation through shear, $\mathrm{Ri}$ the Richardson and $\mathrm{Pr}_\mathrm{t}$ the turbulent Prandtl number. For more information about the Smagorinsky model, the reader is referred to the work of Lilly [1962] or for the implementation into the ICON-LEM to Dipankar et al. (2015).

To capture all SGS processes the approach following Deardorff (1980) takes the full prognostic equation for the SGS-TKE into consideration. This allows tying the SGS eddy diffusivities directly to the SGS-TKE

$$\begin{align} K_\mathrm{m} &= c_\mathrm{m}\, l \, \sqrt{e} \text{ ,} \label{eq:deardiffmom} \tag{7} \\ K_\mathrm{h} &= \left( 1 + 2\,\frac{l}{\Delta} \right)\, K_\mathrm{m} \text{ ,} \label{eq:deardiffheat} \tag{8} \end{align} $$ with a constant $c_\mathrm{m} = 0.1$ . The average grid spacing $\Delta = (\Delta x\Delta y\Delta z)^{1/3}$ is the composite of the grid spacing of the three spatial directions $x$, $y$ and $z$ and $l$ the SGS mixing length. The latter depends on the height above the ground limiting the maximum diameter of eddies near the surface, on the stratification and the average grid spacing $\Delta$. In this work, the modified version of the prognostic equation by Moeng et al. (2000) is being used. It is formulated as

$$\begin{align} l = \begin{cases} \mathrm{min} \left(1.8\,z\,,\,\Delta\,, 0.76\sqrt{e}\,\left( \frac{g}{\theta_\mathrm{v}}\,\frac{\partial\,\theta\mathrm{v}}{\partial \,z}\right)^{-1/2} \right)&, \text{for } \frac{\partial \,\theta\mathrm{v}}{\partial \,z} > 0 \\ \mathrm{min} \,\bigg( 1.8\,z\,,\,\Delta \bigg)&, \text{for } \frac{\partial \,\theta\mathrm{v}}{\partial\, z} \leq 0 \\ \label{eq:mixinglength} \tag{9} \end{cases} \end{align} $$ Equation \ref{eq:fullTKE} uses the same substitution for the production terms as in the Smagorinsky model (Eq. \ref{eq:shearterm} and \ref{eq:buoyterm}). The pressure term is parameterized as $$\begin{align} \frac{1}{\rho}\frac{\partial}{\partial x_i} \left(\overline{\rho\,u_i'(e'+\frac{p'}{\rho})} \right ) &= - \frac{1}{\rho}\frac{\partial}{\partial x_i} \left( 2\,\rho\,K_\mathrm{m}\frac{\partial e}{\partial x_i} \right ) \label{eq:differm} \tag{12} \end{align}$$ and the dissipation $\epsilon$ is given by $$\begin{align} \epsilon = \left( 0.19 + 0.74\,\frac{l}{\Delta} \right)\,\frac{e^{3/2}}{l} \text{ .} \label{eq:dissipationterm} \tag{13} \end{align} $$

Therefore all terms in equation \ref{eq:fullTKE} are defined and can be used as turbulence closure.

The turbulence schemes can be picked in the namelist of your run file in the part &nwp_phy_nml under

inwp_turb = 4 or 5 .

The default LES turbulence scheme is 5 which is the Smagorinsky model. If you want to use the Deardorff scheme, you also have to activate the advection of TKE as a passive tracer under &transport_nml with

iadv_tke = 2

to allow both horizontal and vertical advection. That's it already!

The code for the turbulence closure schemes can be found in your ICON directory in the files <YOUR-ICON-REPO>/src/atm_phy_les/mo_sgs_turbulence.f90 and <YOUR-ICON-REPO>/src/atm_phy_les/mo_sgs_turbmetric.f90. As of today's status, two versions exist: one with and one without metric terms. While the code with metric terms yields the same results as the one without in a no topography case, the development in the framework without metric terms can be easier and therefore this file coexists. In the future, both files could be merged to reduce redundancy and simplify maintenance.