PBL Schemes

Planetary Boundary Layer (PBL) schemes are used to model unresolved transport in the vertical direction within the planetary boundary layer when mesh resolutions are too coarse to resolve any of the turbulent eddies responsible for this transport (~1 km grid resolution or larger). The PBL scheme is used to provide closure for vertical turbulent fluxes (i.e., \(\widetilde{w'\phi'} = \widetilde{w\phi} - \widetilde{w}\widetilde{\phi}\), for any quantity \(\phi\)). PBL schemes may be used in conjunction with an LES model that specifies horizontal turbulent transport, in which case the vertical component of the LES model is ignored.

Right now, ERF supports several PBL schemes: MYNN Level 2.5, MYJ, SHOC, MRF, and YSU.

The MYNN Level 2.5 model is the Mellor-Yamada-Nakanishi-Niino Level 2.5 model, largely matching the original forumulation proposed by Nakanishi and Niino in a series of papers from 2001 to 2009.

MYNN Level 2.5 PBL Model

In this model, the vertical turbulent diffusivities are computed in a local manner based on a gradient diffusion approach with coefficients computed based on a transported turbulent kinetic energy value. The implementation and description here largely follows Nakanishi and Niino, Journal of Meteorological Society of Japan, 2009, but has also been influenced by the full series of papers that led to the development of this model and by a few documents published since then, as listed in the Useful References section below.

The prognostic equation for \(q^2 = \widetilde{u_i u_i} - \widetilde{u}_i\widetilde{u}_i\) is

\[\frac{\partial \bar{\rho} q^2}{\partial t} + \left[ \frac{\partial \bar{\rho} \widetilde{u}_i q^2}{\partial x_i} \right] = \frac{\partial}{\partial z} \left(K_{q,v} \frac{\partial q^2}{\partial z} \right) + 2\bar{\rho} \left(-\widetilde{u'w'} \frac{\partial \widetilde{u}}{\partial z} - \widetilde{v'w'}\frac{\partial \widetilde{v}}{\partial z} + \beta g \widetilde{w'\theta'} - \frac{q^3}{B_1 l} \right)\]

where \(B_1\) is a model parameter, \(\beta\) is the thermal expansion coefficient and l is a lengthscale. The vertical turbulent transport coefficients are then computed:

\[K_{m,v} = l q S_m, K_{q,v} = 3 l q S_m, K_{\theta, v} = l q S_\theta\]

where \(S_m\) and \(S_\theta\) are stability parameters thaat account for buoyancy effects. These coefficients are then applied in evaluating the vertical component of turbulent fluxes in a similar manner as is described for the Smagorinsky LES model. Computation of the stability parameters and lengthscale depend on the Obukhov length and surface heat flux, which are obtained from the sec:MOST. Further detail on these computations can be found in the cited works. Several model coefficients are required, with default values in ERF taken from the work of Nakanishi and Niino.

Useful References

The following references have informed the implementation of the MYNN PBL model in ERF:

Discussions with Branko Kosovic (NCAR) and Joseph B. Olson (NOAA) have also played a major role in informing the implementation of MYNN PBL models in ERF.

MYNN-EDMF Level 2.5 PBL Model

Warning

Implementation is in progress with basic support.

More recent advancements that add significant complexity to the MYNN scheme have been incorporated into WRF, as described in Olson et al. 2019. These advancements are not included in ERF, but may be in the future.

MYJ PBL Model

Warning

Implementation is in progress with basic support.

The Mellor-Yamada-Janjic (MYJ) scheme is a 1.5-order turbulence closure that solves a prognostic equation for turbulent kinetic energy (TKE). It uses a local closure approach with no counter-gradient terms, making it particularly effective for stable and neutral boundary layers.

The turbulent fluxes are computed using gradient diffusion:

\[\overline{w'\phi'} = -K_\phi \frac{\partial \phi}{\partial z}\]

The vertical turbulent transport coefficients are computed from TKE and a master length scale:

\[K_{m,v} = \rho L q S_m, \quad K_{\theta,v} = \rho L q S_h, \quad K_{q,v} = \rho L q S_h\]

where \(q = \sqrt{2\cdot\text{TKE}}\), \(L\) is the master length scale, and \(S_m\), \(S_h\) are stability functions that account for buoyancy effects and depend on the gradient Richardson number. The master length scale \(L\) is diagnosed based on the PBL height, von Kármán’s constant, and height above the surface within the PBL, transitioning to a local mixing length in the free atmosphere.

The prognostic TKE equation includes production by shear and buoyancy, and dissipation:

\[\frac{\partial \text{TKE}}{\partial t} + \nabla \cdot (\mathbf{u} \text{TKE}) = P_s + P_b - \epsilon + \nabla \cdot (K_q \nabla \text{TKE})\]

where \(P_s\) is shear production, \(P_b\) is buoyancy production, and \(\epsilon\) is dissipation.

Closure coefficients are taken from Janjić (2002) NCEP Office Note 437. The implementation in ERF follows Janjić (1994, 2002) and uses the Mellor-Yamada (1982) length scale formulation.

References

  • Janjić, Z. I. (1994): “The Step-Mountain Eta Coordinate Model: Further developments of the convection, viscous sublayer, and turbulence closure schemes”, Monthly Weather Review, 122(5), 927-945.

  • Janjić, Z. I. (2002): “Nonsingular implementation of the Mellor-Yamada Level 2.5 Scheme in the NCEP Meso model”, NCEP Office Note No. 437.

  • Mellor, G. L., & Yamada, T. (1982): “Development of a turbulence closure model for geophysical fluid problems”, Reviews of Geophysics, 20(4), 851-875.

SHOC PBL Model

Warning

Implementation is in progress with basic support.

The Simplified Higher-Order Closure (SHOC) is a unified parameterization that represents both turbulent mixing and shallow convection in a single framework. Originally developed for the Community Atmosphere Model (CAM) and now used in E3SM, SHOC uses prognostic TKE with diagnostic second and third-order moments and assumed probability density functions (PDFs) to represent subgrid-scale variability.

SHOC computes vertical turbulent fluxes for momentum, heat, and moisture, along with subgrid-scale cloud fraction and liquid water content. The assumed PDFs allow the scheme to predict partial cloudiness and transitions between clear and cloudy conditions. The implementation uses higher-order closure equations to diagnose eddy diffusivities and turbulent fluxes, with special treatment for cloud-top entrainment.

References

  • Golaz, J.-C., et al. (2002): “A PDF-based model for boundary layer clouds. Part I: Method and model description”, Journal of the Atmospheric Sciences, 59(24), 3540-3551.

  • Bogenschutz, P. A., & Krueger, S. K. (2013): “A simplified PDF parameterization of subgrid-scale clouds and turbulence for cloud-resolving models”, Journal of Advances in Modeling Earth Systems, 5(2), 195-211.

  • E3SM SHOC Documentation: https://github.com/E3SM-Project/E3SM/tree/master/components/eamxx/src/physics/shoc

MRF PBL Model

Warning

Implementation is in progress with basic support. Need to be tuned in future for real flows.

The Medium Range Forecast (MRF) PBL model is a nonlocal PBL scheme that was originally developed for the MRF model, which was used in the NCEP global forecast system. It is a nonlocal scheme that uses a countergradient diffusion approach to model vertical turbulent transport within the PBL.

The turbulent diffusion for prognostic variables (\(C= u, v, \theta, q_k\)), where \(q_k\) includes all moisture variables is given by

\[\frac{\partial C}{\partial t} = \frac{\partial}{\partial z} \left[ K_c \left( \frac{\partial C}{\partial z} - \gamma_c \right) \right]\]

Here \(K_c\) is the turbulent diffusion coefficient, and \(\gamma_c\) is the countergradient correction term.

The turbulent diffusion coefficient in the mixed layer is given by:

\[K_m = \kappa w_s z \left( 1- \frac{z}{h} \right)^2\]
\[w_s = \frac{u_*}{\phi_m}\]

where \(\kappa\) is the von Karman constant, \(w_s\) is a representative velocity scale in the mixed layer, and \(h\) is the PBL height. The stability function \(\phi_m\) is computed to be consistent with the surface layer bottom. For unstable regime (\(u_*\theta_* < 0\)), it is calculated as follows:

\[\phi_m = \left(1 - 8 sf \frac{h}{L}\right)^{-1/3}\]
\[\phi_{t,q} = \left(1 - 16 sf \frac{h}{L}\right)^{-1/2}\]

and for stable regime (\(u_*\theta_* > 0\)), it is calculated as:

\[\phi_{m,t,q} = \left(1 + 5 sf \frac{h}{L}\right)\]

where \(sf\) is a fraction of the surface layer and atmospheric boundary layer height and \(L\) is the Monin-Obukhov length, which is computed from the surface heat fluxes. The turbulent coefficient for temperature and moisture is given by:

\[K_t = K_q = \frac{K_m}{Pr}\]
\[Pr = \left(\frac{\phi_t}{\phi_m}+ b \kappa sf\right)\]

where \(K_t\) is the turbulent diffusion coefficient for temperature, \(K_q\) is the turbulent diffusion coefficient for moisture and \(Pr\) is the Prandtl number.

The turbulent diffusion coefficient in the free atmosphere is computed from the YSU model as the MRF expressions showed oscillations in the canonical stable boundary layer tests.

\[K_{m,t} = l^2 f_{m,t}(Rig)\left|\frac{\partial U}{\partial z}\right|\]
\[l = \frac{\kappa z \lambda}{\kappa z + \lambda}\]

where \(l\) is the length scale, \(f_{m,t}\) is a stability function for momentum and temperature (or moisture), \(Rig\) is the gradient Richardson number, and \(U\) is the horizontal wind speed. The gradient Richardson number is computed as:

\[Rig = \frac{g}{\theta_v}\left[\frac{\partial \theta_v}{\partial z} \left(\frac{\partial z}{\partial U}\right)^2\right]\]

A different expression is used for the stability function \(f_{m,t}\) for stable and unstable regimes. For stable regime we have,

\[f_t = f_m (1+2.1 Rig) = \frac{1}{\left(1 + 5 Rig\right)^2}\]

For the unstable regime, we have:

\[f_t = 1 - \frac{8 Rig}{1+1.286\sqrt{-Rig}}\]
\[f_m = 1 - \frac{8 Rig}{1+1.746\sqrt{-Rig}}\]

The countergradient correction term is given by:

\[\gamma_c = b \frac{ u_* \theta_*}{w_s}\]

where \(b=7.8\) is a constant, \(u_*\) is the surface frictional velocity scale, \(\theta_*\) is the surface potential temperature scale.

YSU PBL Model

Warning

Implementation is in progress, this option is not yet supported

The Yonsei University (YSU) PBL model is another commonly use scheme in WRF. It includes nonlocal mixing with contergradient diffusion within the PBL, and a local mixing treatment outside the PBL.

Turbulent diffusion for prognostic variables (\(C, u, v, \theta, q_k\)), where \(q_k\) includes all moisture variables and \(C\) any additional scalars (other terms in the equations omitted for brevity):

\[\frac{\partial C}{\partial t} = \frac{\partial}{\partial z} \left[ K_c \left( \frac{\partial C}{\partial z} - \gamma_c \right) - \overline{\left(w'c' \right)_h} \left( \frac{z}{h} \right)^3 \right]\]

Note

Not applied for vertical velocity?

Where for each variable the turbulent diffusion coefficient \(K_c\), countergradient correction \(\gamma_c\), and entrainment flux at the PBL top \(\overline{\left(w'c' \right)_h}\) must be diagnosed for each variable. The main controlling parameter is the PBL height \(h\). Notably, a nonlocal model for turbulent diffusion is used for \(z \leq h\), but a local model is used for \(z \ge h\).

The first step is to determine the PBL height \(h\). This is defined as the smallest value of \(z\) where the bulk Richardson number equals the critical value, which is taken to be 0:

\[{\rm Rib}(z) = \frac{ g \left[ \theta_m(z) - \theta_s\right] }{\theta_{ma} U(z)^2}z\]
\[{\rm Rib}(h) = {\rm Rib_{cr}} = 0\]

where

  • \(\theta_m\) is the moist potential temperature,

  • \(\theta_{ma}\) is the value at the lowest vertical cell in a column,

  • \(U = \sqrt{u^2 + v^2}\) is the horizontal wind speed,

  • \(\theta_s = \theta_{ma} + \theta_T\) is the virtual temperature near the surface,

  • \(\theta_T = a\frac{\overline{\left(w'\theta_m' \right)_0}}{w_{s0}}\) is the excess virtual temperature near the surface,

  • \(a\) is a constant taken to be 6.8 per HND06 (matching the \(b\) constant that appears elsewhere in the YSU model)

  • \(\overline{\left(w'\theta_m' \right)_0}\) is the surface virtual heat flux (determined from the MOST surface layer model),

  • \(w_{s}(z) = \left(u_*^3 + 8 k w_{*b}^3z/h \right)^{1/3}\) is a representative velocity scale in the mixed layer, with \(w_{s0} = w_s(h/2)\) (note this equation matches the WRF implementation and description in H10, but differs from HND06, where \(\phi_m\) appears in place of the constant 8),

  • \(u_*\) is the surface frictional velocity scale determined from the MOST surface layer model,

  • \(k = 0.4\) is the von Karman constant

  • \(w_{*b} = \left[ g/\theta_{ma} \overline{\left(w'\theta_m' \right)_0} h \right]^{1/3}\) for \(\overline{\left(w'\theta_m' \right)_0} > 0\), \(w_{*b} = 0\) otherwise, is a convective velocity scale for moist air

In practice, an approximate value of \(h\) is determined through a two-step process. First, \(\theta_T\) is set to be zero and a provisional value of \(h\) is estimated. Then this provisional value of \(h\) is used to compute \(\theta_T\), which is in turn used to provide an improved estimate of \(h\), which is the value used in subsequent calculations.

Note

This two step-process matches the WRF implementation, but this could be extended iteratively to reach convergence.

Countergradient corrections are computed as follows:

\[\gamma_\theta =\]
\[\gamma_u =\]
\[\gamma_v =\]
\[\gamma_{q_k} = \gamma_C = 0\]

Entrainment fluxes are computed:

\[\overline{\left(w'c' \right)_h} =\]
\[\overline{\left(w'c' \right)_h} =\]

Within the PBL (\(z \leq h\)),

Useful References

The following references have informed the implementation of the MRF and YSU model in ERF: