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, native SHOC, optional EAMxx 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:

Mellor, Journal of the Atmospheric Sciences, 1973: Introduces a PBL model based on $q^2$
Mellor and Yamada, Journal of the Atmospheric Sciences, 1974: Introduces PBL Model Hierarchy (Levels 1-4)
Mellor and Yamada, Reviews of Geophysics and Space Physics, 1982: Introduces Level 2.5 Model
Nakanishi, Boundary-Layer Meteorology, 2001: Fits new model coefficients and proposes new diagnostic equation for the length scale
Nakanishi and Niino, Boundary-Layer Meteorology, 2004: Extends the MYNN PBL modeling framework for moist conditions
Nakanishi and Niino, Boundary-Layer Meteorology, 2006: Numerical stability improvements for the MYNN PBL modeling framework
Nakanishi and Niino, Journal of the Meteorological Society of Japan, 2009: Summary of MYNN model development, re-evaluation of coefficients, and additional demonstration cases
Skamarock et al., A Description of the Advanced Research WRF Model Version 4, 2021: Description of the models implemented in WRF
Olson et al., A Description of the MYNN-EDMF Scheme and the Coupling to Other Components in WRF–ARW, 2019: Description of more recent advancements upon the MYNN model
Juliano et al., Monthly Weather Review, 2022: Description of a 3D generalization Mellor-Yamada PBL models

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¶

The Simplified Higher-Order Closure (SHOC) represents unresolved turbulence, shallow convection, and subgrid-scale cloud macrophysics. It uses prognostic turbulent kinetic energy (TKE), diagnostic higher-order moments, and an assumed probability density function (PDF) for thermodynamic variability. The PDF lets the scheme represent partial cloudiness within a grid cell.

SHOC follows the approach described by Bogenschutz and Krueger (2013). The ERF-native implementation is based on the E3SM/EAMxx SHOC algorithms and source structure. It is not a bit-for-bit port of EAMxx SHOC. It is adapted to ERF data structures, AMReX loops, ERF surface-flux coupling, diagnostics, and build systems.

Selecting SHOC in a simulation¶

Native SHOC is built in tree. It needs no EAMxx, EKAT, or Kokkos setup. Select it at runtime with:

zlo.type = "surface_layer"
erf.pbl_type = NATIVE_SHOC

The surface_layer lower boundary condition supplies the surface fluxes used by SHOC. SHOC-family PBL schemes require this lower boundary type. ERF reads erf.pbl_type with the usual one-or-per-level rule. A single value applies to every level. A list can specify one value per level.

Native SHOC and EAMxx SHOC¶

ERF has two SHOC paths:

Runtime value	Meaning
`NATIVE_SHOC`	Selects the ERF-native implementation in `Source/PBL/Shoc`.
`EAMXX_SHOC`	Selects the optional EAMxx interface in `Source/PhysicsInterfaces/Shoc`.
`SHOC`	Deprecated alias for `EAMXX_SHOC`.

Native SHOC is column based. Like ERF’s other column physics, it requires each AMReX box on a SHOC-active level to span the full vertical domain. Do not use a grid decomposition that splits boxes in the vertical direction. If an input file sets vector-valued grid sizing controls, choose a vertical size at least as large as the level’s vertical cell count. With AMR, SHOC-active refined grids must also cover full vertical columns.

The implementation is AMReX-native and lives in Source/PBL/Shoc.

Use NATIVE_SHOC for the native ERF implementation. Use EAMXX_SHOC only when you build and run the optional EAMxx path.

Surface-flux and microphysics coupling¶

SHOC needs lower-boundary heat, moisture, and momentum fluxes. It does not compute those exchanges directly from a land surface model. Instead, ERF computes the lower-boundary fluxes before SHOC runs and passes the resulting flux arrays to SHOC.

Those fluxes come through ERF’s surface-layer infrastructure. They may come from Monin-Obukhov similarity theory (MOST), prescribed surface-layer inputs, or an active land or ocean surface model when that model provides fluxes. Native SHOC then consumes the ERF flux arrays. Internally it converts ERF’s host density-weighted fluxes to kinematic surface fluxes.

Set the lower boundary to surface_layer for SHOC runs:

zlo.type = "surface_layer"

SHOC diagnoses subgrid non-precipitating cloud partitioning with its assumed PDF. This includes cloud fraction and non-precipitating liquid water. To avoid double counting, ERF disables the saturation-adjustment or condensation step in the microphysics package when a SHOC-family PBL scheme is active.

Native SHOC remains a liquid-cloud macrophysics closure under this interim contract. It may use pre-existing cloud ice for phase-aware thermodynamics and buoyancy, but it does not create or repartition cloud ice.

This does not disable microphysics. Microphysics still handles precipitating processes outside SHOC’s cloud macrophysics role. Choose a moisture model that matches the case. Number-aware microphysics layouts with cloud-droplet or ice number concentrations still need an explicit number closure in their own microphysics pathways if they are coupled to SHOC. Native SHOC state_update rejects those number-aware layouts until a number closure is implemented.

Transport modes¶

Native SHOC has one scalar/cloud/TKE transport selector and one independent horizontal-momentum selector:

erf.shoc.transport_mode = state_update

or:

erf.shoc.transport_mode = host_diffusion

state_update is the default. In this mode, SHOC applies its coupled heat, moisture, non-precipitating cloud liquid, carried cloud ice, and TKE update before the dycore sees the state.

The momentum selector is:

erf.shoc.momentum_transport = host_diffusion

This is the default. In the mixed native SHOC mode, SHOC exports only the momentum diffusivity to ERF’s host diffusion path while keeping the scalar transport on the state_update path.

Other momentum choices are:

erf.shoc.momentum_transport = none
erf.shoc.momentum_transport = state_update

none disables SHOC momentum transport entirely. state_update keeps the legacy direct-velocity update available for debugging or targeted experiments, but it is not the default.

host_diffusion remains available. In this mode, SHOC exports the full vertical eddy-diffusivity block to ERF’s host diffusion path. SHOC does not apply the pre-dycore state update. This mode is currently limited to dry/no- moisture configurations and requires erf.shoc.momentum_transport = host_diffusion.

Runtime options¶

Native SHOC reads options from the erf.shoc namespace. The defaults are the recommended starting point. Most options tune the closure or enable diagnostics.

7 Native SHOC runtime options¶
Option	Default	Values	Notes
`erf.shoc.transport_mode`	`state_update`	`state_update`, `host_diffusion`	Selects the native SHOC scalar/cloud/TKE transport path. Host diffusion is currently dry/no-moisture only.
`erf.shoc.momentum_transport`	`host_diffusion`	`none`, `state_update`, `host_diffusion`	Selects the native SHOC horizontal-momentum path. `host_diffusion` is the mixed-mode default.
`erf.shoc.lambda_low`	`0.001`	Real > 0	Lower bound for the SHOC length-scale formula.
`erf.shoc.lambda_high`	`0.04`	Real >= `lambda_low`	Upper bound for the SHOC length-scale formula.
`erf.shoc.lambda_slope`	`2.65`	Real	Slope used by the SHOC length-scale formula.
`erf.shoc.lambda_thresh`	`0.02`	Real	Threshold used by the SHOC length-scale formula.
`erf.shoc.thl2tune`	`1.0`	Real	Tuning factor for liquid-water potential-temperature variance.
`erf.shoc.qw2tune`	`1.0`	Real	Tuning factor for total-water variance.
`erf.shoc.qwthl2tune`	`1.0`	Real	Tuning factor for total-water and liquid-water potential-temperature covariance.
`erf.shoc.w2tune`	`1.0`	Real	Tuning factor for vertical-velocity variance.
`erf.shoc.length_fac`	`0.5`	Real > 0	Factor used in the SHOC turbulent length scale.
`erf.shoc.c_diag_3rd_mom`	`7.0`	Real	Coefficient used by the diagnostic third-moment closure.
`erf.shoc.coeff_kh`	`0.1`	Real >= 0	Scalar diffusivity coefficient.
`erf.shoc.coeff_km`	`0.1`	Real >= 0	Momentum diffusivity coefficient.
`erf.shoc.top_taper_depth`	`0.0`	Real >= 0	Depth of the upper taper for SHOC moments.
`erf.shoc.top_taper_min_factor`	`0.0`	Real in [0, 1]	Minimum factor applied by the upper taper.
`erf.shoc.shoc_1p5tke`	`false`	Boolean	Switches to the 1.5-TKE moment formulation.
`erf.shoc.extra_shoc_diags`	`false`	Boolean	Writes extra cloud and flux diagnostics.
`erf.shoc.apply_tms`	`false`	Boolean	Applies turbulent mountain stress.
`erf.shoc.check_flux_state`	`false`	Boolean	Checks flux-state consistency in debug runs.
`erf.shoc.column_conservation_check`	`false`	Boolean	Checks column-integrated conservation in debug runs.
`erf.shoc.debug_summary`	`false`	Boolean	Prints a short runtime summary for each advance call.
`erf.shoc.allow_tendency_microphysics_overlap`	`false`	Boolean	Allows host-applied source terms to overlap microphysics coupling.
`erf.shoc.signed_tke_production`	`false`	Boolean	Keeps signed buoyancy production in the TKE budget.

ERF checks these option ranges at startup:

erf.shoc.lambda_low > 0
erf.shoc.lambda_high >= erf.shoc.lambda_low
erf.shoc.length_fac > 0
erf.shoc.coeff_kh >= 0
erf.shoc.coeff_km >= 0
erf.shoc.top_taper_depth >= 0
erf.shoc.top_taper_min_factor is in [0, 1]

Diagnostics¶

Native SHOC can write plotfile and 1D profile diagnostics when it runs in state_update mode. Request SHOC diagnostics through the normal plotfile variable lists. For example:

erf.plot_vars_1 = density rhotheta theta qv qc Kmv Khv Lturb shoc_cldfrac shoc_ql wthv_sec

Standard turbulence diagnostics that can use SHOC eddy diffusivities include:

nut
Kmv
Khv
Lturb

Native SHOC also provides these diagnostic plot variables:

shoc_cldfrac
shoc_ql
shoc_ql2
shoc_cond
wqls_sec
wthv_sec
w_sec
thl_sec
qw_sec
qwthl_sec
wthl_sec
wqw_sec
w3
brunt
isotropy
shear_prod
buoy_prod
diss_tke

These diagnostics describe the diagnosed cloud field, second and third moments, and TKE budget terms. Native-only diagnostics are meaningful only when erf.pbl_type = NATIVE_SHOC and the native driver has produced diagnostics. In other cases, these variables may be unavailable or may carry missing-value placeholders.

erf.shoc.extra_shoc_diags is a diagnostic and developer option. Do not rely on it as the only control for plotfile output. Request plotfile fields explicitly with erf.plot_vars_1 or erf.plot_vars_2.

References¶

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 boundary layer 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 implementation in ERF follows the original Hong and Pan (1996) formulation with several modern enhancements described below.

Primary Reference: Hong, S. Y., and H.-L. Pan, 1996: Nonlocal Boundary Layer Vertical Diffusion in a Medium-Range Forecast Model. Monthly Weather Review, 124, 2322-2339. https://doi.org/10.1175/1520-0493(1996)124<2322:NBLVDI>2.0.CO;2

WRF Implementation: module_bl_mrf.F https://github.com/wrf-model/WRF/blob/master/phys/module_bl_mrf.F

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 (nonlocal flux).

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}\left[\frac{\partial \theta}{\partial z} \left(\frac{\partial U}{\partial z}\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.

Note

The countergradient correction term is now optional in ERF and can be enabled via the enable_mrf_countergradient flag in the input file. By default, it is disabled to maintain backward compatibility. Set erf.enable_mrf_countergradient = true to enable this feature.

MRF Model Enhancements and Extensions¶

The MRF model has been enhanced with several optional features in ERF that extend beyond the original Hong and Pan (1996) formulation, aligning more closely with WRF’s advanced parameterizations and modern understanding of boundary layer physics. All features are disabled by default to maintain backward compatibility with existing simulations. This section documents all enhancements and their physical justification.

Important Note on Physics Fidelity: When using the original MRF scheme without enhancements, users should cite Hong and Pan (1996) and acknowledge that they are using the standard formulation. When any enhancement is enabled, proper attribution to ERF documentation and the referenced papers is required.

1. VPERT (Virtual Potential Temperature Perturbation) Correction¶

Issue Fixed: The original WRF MRF code (module_bl_mrf.F line 878) includes the following operation:

VPERT = HGAMT + EP1*THX(I,KL)*HGAMQ
VPERT = MIN(VPERT, GAMCRT)  ! <-- This is incorrect physics

Problem: After combining sensible heat (HGAMT) and latent heat (HGAMQ) contributions, limiting VPERT to GAMCRT suppresses the effect of latent heating. This violates the physics because:

HGAMT is already limited to GAMCRT (max sensible heating)
Adding positive HGAMQ can legitimately produce VPERT > GAMCRT
The latent heating from evaporation in unstable conditions is being artificially suppressed
This causes underprediction of PBL height in strongly moist convective conditions

ERF Correction: ERF can implement unbounded VPERT when enable_mrf_unbounded_vpert is enabled:

const Real VPERT = amrex::max(HGAMT + 0.61 * theta * HGAMQ, zero);

Otherwise, for strict WRF consistency, it bounds VPERT to GAMCRT (3 K):

const Real VPERT = amrex::max(amrex::min(HGAMT + 0.61 * theta * HGAMQ, GAMCRT), zero);

Enabling unbounded VPERT preserves the combined heating effect while preventing negative VPERT values (via MAX with zero). The correction produces physically more accurate PBL heights in moist environments.

References: - Original error identified during ERF development - See ERF_ComputeDiffusivityMRF.cpp source comments for detailed physics explanation

2. HGAMQ Moisture Countergradient with Proper Limiting¶

Sign Convention Clarification: The WRF convention uses: - $q_*$ is negative for upward moisture flux (evaporation from surface) - $u_*$ is positive (friction velocity magnitude) - Therefore: $-const_b \cdot u_* \cdot q_*$ gives positive HGAMQ for unstable (evaporating) conditions

Missing WRF Safeguard: WRF applies MAX limiting (HGAMQ = MAX(HGAMQ, 0.0)) but only after the MIN operation. ERF correctly applies full bounding:

HGAMQ = max(min(-const_b * u_* * q_* / w_*, GAMCRQ), zero);

This prevents unrealistic negative moisture countergradient if $q_*$ becomes positive (condensation), which would indicate upside-down countergradient flux.

Land/Water Discrimination: HGAMQ is zeroed over water surfaces because: - Evaporation over water is implicitly handled by ocean/water body parameterizations - Countergradient moisture fluxes are primarily a land phenomenon (soil moisture limitation) - This follows WRF’s approach (module_bl_mrf.F line 876)

3. Cloud-Aware Stability Function Adjustments¶

Purpose: Improves representation of cloudy boundary layers by modulating stability functions based on cloud water/ice content.

Implementation Details:

Cloud detection threshold: $q_c + q_i > 0.1 \text{ g/kg}$ (1e-4 kg/kg)
In stable layers with clouds: Reduce stability damping by 10-20% where cloud content exceeds threshold
- Physics: Clouds reduce vertical oscillations through radiative and latent effects, making the layer less suppressive to turbulence
- Improved representation of fog/stratus-topped stable layers
In unstable layers with clouds: Slightly enhance instability enhancement by 5% where cloud content exceeds threshold
- Physics: Latent heat release from condensation enhances buoyancy
- Better captures cumulus-capped boundary layers

Parameters: - Enabled via: enable_mrf_countergradient flag (default: false) - Adjustment strength: 15-20% reduction in stable, 5% boost in unstable - Can be customized via pbl_mrf_cloud_adjustment_factor parameter

Physical Justification: - Clouds modify vertical buoyancy structure through radiative cooling/warming - Latent heat release enhances convective mixing - Cloud-top entrainment zones are qualitatively different from clear-air turbulence - Conceptually similar to WRF’s IMVDIF cloud-aware parameterization (Bretherton & Park 2009)

References: - Bretherton, C. S., and S. Park, 2009: A new moist turbulence parameterization in the WRF

Advanced Research WRF (ARW) model. In Proceedings of the 9th Annual WRF Users’ Workshop.

4. Virtual Potential Temperature Treatment¶

Enhancement: Proper handling of moisture effects on buoyancy throughout the scheme.

Computes $\theta_v = \theta(1 + 0.61 \cdot q_v)$ at all levels
Used in Richardson number diagnosis: $Rig = \frac{g}{\theta_v} \frac{d\theta_v}{dz} / (shear)^2$
Critical for accurate PBL height and convective strength

This is standard meteorological practice and implicitly assumed in Hong & Pan (1996), now made explicit.

5. Free Atmosphere Mixing via YSU Stability Functions¶

Reason for Enhancement: Original MRF diffusivity formulation caused oscillations in canonical stable boundary layer tests (GABLS cases).

Solution: Use YSU scheme (Hong et al. 2006, Appendix A) Richardson number-dependent mixing above PBL:

More stable numerically
Better represents free atmosphere mixing
Standard in modern WRF configurations
Detailed equations and references provided in source code comments

References: - Hong, S. Y., Y. Noh, and J. Dudhia, 2006: A new vertical diffusion package with an explicit

treatment of entrainment processes. Monthly Weather Review, 134, 2318-2341. https://doi.org/10.1175/MWR3250.1

Older MRF Enhancements (Deprecated/Documented for Historical Completeness)¶

These features are historically documented but superseded by the above enhancements:

Iterative Thermal Excess Correction¶

Purpose: The thermal excess correction ($\theta_T$) modifies the surface virtual potential temperature to account for cumulative heating effects in the PBL. The iterative method refines this estimate through multiple passes, similar to WRF’s HGAMT/HGAMQ formulation, providing better estimates than the single-pass simple method:

\[\theta_T = -b \frac{u_* \theta_*}{w_*}\]

Parameters:

pbl.enable_mrf_iterative_thermal_excess (bool): Enable iterative refinement (default: false)
pbl.mrf_thermal_excess_iterations (int): Number of refinement iterations (default: 3)

Usage Example:

pbl.pbl_type = "mrf"
pbl.enable_mrf_iterative_thermal_excess = true
pbl.mrf_thermal_excess_iterations = 3

Physical Justification:

Simple method: Valid for most stable/neutral conditions
Iterative method: Better for strong convective conditions with multiple heating cycles
Convergence: Typically reaches solution within 3-5 iterations

Saturated Layer Handling (IMVDIF)¶

Purpose: Modulates moisture diffusivity in cloud/saturated layers using a height-dependent cloud fraction estimate. This addresses the reduced diffusivity characteristic of stratocumulus-topped boundary layers and fog/stratus evolution.

Parameters:

pbl.enable_mrf_cloudy_layers (bool): Enable cloud-aware modulation (default: false)
pbl.mrf_cloud_diffusivity_factor (Real): Diffusivity reduction factor in cloudy layers (default: 0.8)
- Range: [0, 1]
- 1.0 = no reduction (default diffusivity)
- 0.8 = 20% reduction in cloudy regions
- 0.0 = complete suppression

Implementation Details:

Uses a simple heuristic: reduces diffusivity in lower layers (z < 2000 m)
Gradual transition: full reduction near surface, tapering upward
Can be enhanced with actual cloud fraction variables from microphysics

Usage Example:

pbl.pbl_type = "mrf"
pbl.enable_mrf_cloudy_layers = true
pbl.mrf_cloud_diffusivity_factor = 0.8

Physical Justification:

Stratocumulus layers are less turbulent than well-mixed layers
Reduced diffusivity better captures cloud-top entrainment dynamics
Improves representation of fog/stratus evolution

Countergradient Term Bounding¶

Purpose: Applies bounds to countergradient (non-local) turbulent flux terms, preventing unrealistic amplification of fluxes. This follows WRF’s approach with GAMCRT and GAMCRQ parameters.

Parameters:

pbl.enable_mrf_countergradient_bounds (bool): Enable bounds (default: false)
pbl.mrf_countergradient_max_theta (Real): Maximum heat countergradient (default: 3.0)
- WRF equivalent: GAMCRT = 3
pbl.mrf_countergradient_max_q (Real): Maximum moisture countergradient (default: 0.002)
- WRF equivalent: GAMCRQ = 2E-3

Implementation Details:

The bounds are applied as:

cg_theta = min(countergradient_theta, mrf_countergradient_max_theta)
cg_q = min(countergradient_q, mrf_countergradient_max_q)

Usage Example:

pbl.pbl_type = "mrf"
pbl.enable_mrf_countergradient = true
pbl.enable_mrf_countergradient_bounds = true
pbl.mrf_countergradient_max_theta = 3.0
pbl.mrf_countergradient_max_q = 0.002

Physical Justification:

Prevents countergradient fluxes from exceeding realistic bounds
Critical for very strong convective conditions
Improves model stability in extreme parameter regimes

Iterative Thermal Excess Correction¶

Purpose: The thermal excess correction ($\theta_T$) modifies the surface virtual potential temperature to account for cumulative heating effects in the PBL. The iterative method refines this estimate through multiple passes, similar to WRF’s HGAMT/HGAMQ formulation, providing better estimates than the single-pass simple method:

\[\theta_T = -b \frac{u_* \theta_*}{w_*}\]

Parameters:

pbl.enable_mrf_iterative_thermal_excess (bool): Enable iterative refinement (default: false)
pbl.mrf_thermal_excess_iterations (int): Number of refinement iterations (default: 3)

Usage Example:

pbl.pbl_type = "mrf"
pbl.enable_mrf_iterative_thermal_excess = true
pbl.mrf_thermal_excess_iterations = 3

Physical Justification:

Simple method: Valid for most stable/neutral conditions
Iterative method: Better for strong convective conditions with multiple heating cycles
Convergence: Typically reaches solution within 3-5 iterations

Saturated Layer Handling (IMVDIF)¶

Purpose: Modulates moisture diffusivity in cloud/saturated layers using a height-dependent cloud fraction estimate. This addresses the reduced diffusivity characteristic of stratocumulus-topped boundary layers and fog/stratus evolution.

Parameters:

pbl.enable_mrf_cloudy_layers (bool): Enable cloud-aware modulation (default: false)
pbl.mrf_cloud_diffusivity_factor (Real): Diffusivity reduction factor in cloudy layers (default: 0.8)
- Range: [0, 1]
- 1.0 = no reduction (default diffusivity)
- 0.8 = 20% reduction in cloudy regions
- 0.0 = complete suppression

Implementation Details:

Uses a simple heuristic: reduces diffusivity in lower layers (z < 2000 m)
Gradual transition: full reduction near surface, tapering upward
Can be enhanced with actual cloud fraction variables from microphysics

Usage Example:

pbl.pbl_type = "mrf"
pbl.enable_mrf_cloudy_layers = true
pbl.mrf_cloud_diffusivity_factor = 0.8

Physical Justification:

Stratocumulus layers are less turbulent than well-mixed layers
Reduced diffusivity better captures cloud-top entrainment dynamics
Improves representation of fog/stratus evolution

Countergradient Term Bounding¶

Purpose: Applies bounds to countergradient (non-local) turbulent flux terms, preventing unrealistic amplification of fluxes. This follows WRF’s approach with GAMCRT and GAMCRQ parameters.

Parameters:

pbl.enable_mrf_countergradient_bounds (bool): Enable bounds (default: false)
pbl.mrf_countergradient_max_theta (Real): Maximum heat countergradient (default: 3.0)
- WRF equivalent: GAMCRT = 3
pbl.mrf_countergradient_max_q (Real): Maximum moisture countergradient (default: 0.002)
- WRF equivalent: GAMCRQ = 2E-3

Implementation Details:

The bounds are applied as:

cg_theta = min(countergradient_theta, mrf_countergradient_max_theta)
cg_q = min(countergradient_q, mrf_countergradient_max_q)

Usage Example:

pbl.pbl_type = "mrf"
pbl.enable_mrf_countergradient = true
pbl.enable_mrf_countergradient_bounds = true
pbl.mrf_countergradient_max_theta = 3.0
pbl.mrf_countergradient_max_q = 0.002

Physical Justification:

Prevents countergradient fluxes from exceeding realistic bounds
Critical for very strong convective conditions
Improves model stability in extreme parameter regimes

High-Resolution Grid-Dependent Diffusivity Bounds¶

Purpose: Applies grid-dependent (high-resolution) bounds to diffusivity coefficients, permitting stronger and more resolved mixing in fine-resolution Large Eddy Simulations (LES) or high-resolution mesoscale runs ($Delta z < 100text{m}$) compared to conservative global forecast limits.

Parameters:

pbl.pbl_mrf_highres_bounds (bool): Enable high-resolution grid-dependent bounds (default: false)

Implementation Details:

When enabled, the diffusivity bounds scale with local grid-spacing and density as formulated in Hong et al. (2006):

Conservative Bounds (Default): $K_{min} = 0.1\ \text{m}^2/\text{s}, K_{max} = 300\ \text{m}^2/\text{s}$
High-Resolution Bounds: $K_{min} = 0.001 \times \Delta z \times \rho, K_{max} = 1000\ \text{m}^2/\text{s}$

Combined Usage Examples¶

Default Configuration (All Features Off):

# Standard MRF configuration (existing parameters)
pbl.pbl_type = "mrf"
pbl.pbl_mrf_Ribcr = 0.5
pbl.pbl_mrf_const_b = 7.8
pbl.pbl_mrf_sf = 0.1

Use this configuration for:

Standard ABL simulations
Backward compatibility
Neutral/stable boundary layers

Conservative Enhancement:

pbl.pbl_type = "mrf"
pbl.enable_mrf_iterative_thermal_excess = true
pbl.mrf_thermal_excess_iterations = 3

Use for: Improving strong convection without other changes

Cloud-Aware Configuration:

pbl.pbl_type = "mrf"
pbl.enable_mrf_cloudy_layers = true
pbl.mrf_cloud_diffusivity_factor = 0.8

Use for: Stratocumulus-topped or fog/stratus simulations

Full WRF-Style Configuration:

pbl.pbl_type = "mrf"
pbl.enable_mrf_iterative_thermal_excess = true
pbl.mrf_thermal_excess_iterations = 3
pbl.enable_mrf_cloudy_layers = true
pbl.mrf_cloud_diffusivity_factor = 0.8
pbl.enable_mrf_countergradient = true
pbl.enable_mrf_countergradient_bounds = true
pbl.mrf_countergradient_max_theta = 3.0
pbl.mrf_countergradient_max_q = 0.002

Use for: Maximum realism in challenging conditions

Performance Impact¶

Iterative thermal excess: ~1-5% additional computation
Cloud-aware moisture: <1% overhead (simple height-based heuristic)
Countergradient bounds: <1% overhead (minimal arithmetic)
Total impact with all enabled: ~1-5% runtime increase

Useful References¶

Hong et al. (1996): “Nonlocal Boundary Layer Vertical Diffusion in a Medium-Range Forecast Model”
Hong et al. (2006): “A new vertical diffusion package with an explicit treatment of entrainment processes”
WRF Model Documentation: PBL Schemes
Skamarock et al., A Description of the Advanced Research WRF Model Version 4, 2021 <http://dx.doi.org/10.5065/1dfh-6p97>

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} =\]

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

Useful References¶

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

[H10] Hong, Quarterly Journal of the Royal Meteorological Society, 2010: Most up-to-date YSU model formulation as implemented in WRF, with revisions for stable boundary layers
[HND06] Hong, Noh, and Dudhia, Monthly Weather Review, 2006: Initial formulation referred to as the YSU model, adds improved entrainment formulation (relative to NCHR03) to work of TM86 and a few other modifications
[NCHR03] Noh, Cheon, Hong, and Raasch, Boundary-Layer Meteorology, 2003: Entrainment effects added to TM86
[HP96] Hong and Pan, Monthly Weather Review, 1996: Largely an implementation and evaluation of TM86
[TM86] Troen and Mahrt, Boundary-Layer Meteorology, 1986: Initial incorporation of nonlocal counter-graident term in vertical diffusion model
[WF18] Wilson and Fovell, Weather and Forecasting, 2018: Extension of YSU to handle interplay between radiation and fog, active in WRF with the ysu_topdown_pblmix = 1 option
The WRF Fortran source code for this module as of Dec. 2023. The ERF implementation supports the same physical models as this WRF implementation, with the exception of the ysu_topdown_pblmix = 1 option from WF18, i.e. the implementation in ERF largely matches the PBL scheme described in H10.