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:

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 boundary-layer turbulence, shallow convection, and subgrid-scale cloud macrophysics. SHOC uses prognostic turbulent kinetic energy (TKE), diagnosed turbulent moments, and an assumed probability density function (PDF) for thermodynamic variability. The PDF allows a grid cell to contain a diagnosed fractional cloud field rather than requiring the entire cell to be either clear or cloudy.

SHOC follows the approach described by Bogenschutz and Krueger (2013). ERF provides both an ERF-native SHOC implementation and an optional interface to EAMxx SHOC. The native implementation follows the E3SM/EAMxx SHOC algorithmic structure but is adapted to ERF state variables, AMReX data structures, ERF surface-flux coupling, and ERF time integration. It should not be expected to reproduce EAMxx SHOC bit-for-bit.

Selecting a SHOC implementation

ERF keeps both paths because they serve different integration goals. Native SHOC is adapted to ERF’s AMReX data structures, surface coupling, and time integration and therefore is not expected to reproduce EAMxx SHOC bit-for-bit; the EAMxx path is retained for users who need the E3SM/EAMxx reference implementation path.

Native SHOC is part of the ERF source tree and does not require EAMxx, EKAT, or Kokkos setup.

ERF recognizes the following SHOC-family PBL selections:

Runtime value

Meaning

NATIVE_SHOC

ERF-native SHOC in Source/PBL/Shoc.

EAMXX_SHOC

Optional EAMxx SHOC interface. The executable must be built with ERF_ENABLE_EAMXX_SHOC=ON.

SHOC

Deprecated alias for EAMXX_SHOC. It does not select native SHOC.

Native SHOC quick start

A minimal Native-SHOC configuration is

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

The Native transport defaults are already

erf.shoc.transport_mode = state_update
erf.shoc.momentum_transport = state_update

so those two lines do not need to be specified unless an input file benefits from making the choice explicit.

Before running Native SHOC, note three important restrictions:

  • the lower boundary must use zlo.type = "surface_layer";

  • every AMReX box on a Native-SHOC-active level must span the full vertical domain;

  • moisture layouts carrying cloud-liquid or cloud-ice number concentrations are currently rejected because Native SHOC does not yet provide the required number closure.

The remaining sections describe these requirements and the available tuning and diagnostic controls in more detail.

SHOC-family PBL schemes require zlo.type = "surface_layer". The surface layer supplies the lower-boundary heat, moisture, and momentum fluxes consumed by SHOC.

Vertical grid decomposition

Native SHOC is a column scheme. Every AMReX box on a native-SHOC-active level must span that level’s full vertical domain. A decomposition that splits a box in the vertical direction is rejected at runtime.

When setting vector-valued AMR grid-size controls, choose a vertical grid size large enough to keep each SHOC box full height. With mesh refinement, SHOC-active refined grids must also consist of full vertical columns; setting amr.refine_whole_domain_dir = 2 makes every refined level cover the full depth of the domain regardless of where cells are tagged. See Mesh Refinement for the general ERF refinement rules and Refining the Whole Domain in One Direction for that parameter.

Surface fluxes and moisture coupling

Native SHOC consumes the lower-boundary flux arrays produced by ERF’s surface-layer infrastructure. Those fluxes may be based on Monin-Obukhov similarity theory, prescribed surface-layer forcing, or an active surface model that supplies the ERF surface fluxes. ERF stores the host fluxes in density-weighted form; native SHOC converts them internally to the kinematic fluxes used by the SHOC closure.

The native SHOC PDF diagnoses nonprecipitating cloud liquid and cloud fraction. For SHOC-family configurations, ERF suppresses the separate microphysics condensation/saturation-adjustment step that would otherwise duplicate SHOC’s cloud macrophysics role. This does not disable the rest of the microphysics scheme; precipitating processes remain the responsibility of microphysics.

Native SHOC uses pre-existing cloud ice in its phase-aware thermodynamics. The current PDF reconstruction preserves that existing ice but does not create or repartition cloud ice. The PDF repartitions the nonprecipitating liquid component.

For moist native SHOC, use the native state_update transport path described below. Native state_update currently rejects moisture layouts containing cloud-water or cloud-ice number concentration components because a compatible number closure has not yet been implemented.

Native SHOC transport modes

Native SHOC separates thermodynamic/moisture/TKE transport from the choice used for horizontal momentum.

The default native configuration is:

erf.shoc.transport_mode = state_update
erf.shoc.momentum_transport = state_update

This is the normal starting point for native SHOC.

erf.shoc.transport_mode = state_update

Native SHOC applies its coupled thermodynamic, moisture/cloud, and TKE column update to the ERF state before the dycore advances that state. This is the required native transport mode for moist SHOC configurations.

erf.shoc.transport_mode accepts only state_update. The former host_diffusion value is removed and is rejected at startup with migration guidance.

The former value erf.shoc.transport_mode = tendencies has been removed and is rejected at startup.

Horizontal momentum is controlled independently:

erf.shoc.momentum_transport = state_update

Apply the native SHOC horizontal-velocity column increment directly to the ERF face velocities. This is the default production path.

erf.shoc.momentum_transport = none

Disable the Native-SHOC horizontal-momentum transport contribution. ERF retains generic host ownership of vertical momentum diffusion in this mode. This does not by itself create a nonzero host momentum diffusivity; it only leaves the generic ERF momentum-diffusion path eligible to act when another configured host closure supplies one.

The former host_diffusion momentum value is removed. Use state_update to retain SHOC momentum transport or none to disable it.

PBL height and turbulent structure

Native SHOC diagnoses planetary boundary layer height from a bulk-Richardson criterion, with an additional convective correction for unstable surface forcing. Native SHOC pblh is reported in metres above local ground (AGL). It is not an absolute geometric height.

The diagnosed PBL depth also enters the native surface second-moment scaling under positive surface heat flux, so the convective velocity scale uses the current diagnosed boundary-layer depth rather than a fixed depth.

Native SHOC diagnoses a turbulent mixing length, Brunt-Vaisala frequency squared, momentum and heat diffusivities, an isotropy timescale, second and third moments, and the PDF cloud quantities used by the cloud macrophysics closure.

TKE behavior and reduced 1.5-order mode

Native SHOC uses ERF’s prognostic TKE state. At model initialization ERF sets TKE from erf.tke_min (default 1.e-6 m^2/s^2) before any problem-specific or surface-layer-based TKE initialization is applied.

erf.tke_min is an initialization control, not the internal Native-SHOC TKE floor. During each Native-SHOC update the closure currently bounds TKE to

\[4\times10^{-4} \le e \le 50 \qquad \mathrm{m^2\,s^{-2}},\]

where \(e\) denotes TKE. Consequently a run initialized below \(4\times10^{-4}\ \mathrm{m^2\,s^{-2}}\) can move to the Native-SHOC internal floor on the first SHOC update.

In the default higher-order form (erf.shoc.shoc_1p5tke = false), the Native-SHOC buoyancy-production term is

\[P_b = \frac{g}{300\ \mathrm{K}}\,\overline{w'\theta_v'}.\]

Setting

erf.shoc.shoc_1p5tke = true

selects the reduced 1.5-order form. In that mode the buoyancy term is

\[P_b = -K_h N^2,\]

using the carried heat diffusivity and the diagnosed squared Brunt–Väisälä frequency \(N^2\). The higher-order thermodynamic variance/covariance structure and vertical-velocity third moment are suppressed in this mode, so shoc_1p5tke changes the closure rather than merely changing diagnostics.

The treatment of the total production term is controlled independently by erf.shoc.signed_tke_production. With the default value false,

\[P = \max(0,\,P_s + P_b),\]

where \(P_s\) is shear production. With

erf.shoc.signed_tke_production = true

the clipping is removed:

\[P = P_s + P_b.\]

The latter permits negative buoyancy production in stable conditions to reduce TKE directly. The default remains false; use a non-default value only when the intended experiment calls for the signed-production formulation.

Native stability and diffusivity tuning

Native SHOC diagnoses the squared Brunt–Väisälä frequency from virtual potential temperature,

\[N^2 = \frac{g}{\theta_v} \frac{\partial\theta_v}{\partial z}.\]

Positive \(N^2\) denotes stable stratification in this diagnostic.

On the general Native diffusivity path, the kinematic momentum and heat diffusivities are

\[K_m = C_{K_m}\,\tau\,e, \qquad K_h = C_{K_h}\,\tau\,e,\]

where \(e\) is TKE, \(\tau\) is the diagnosed isotropy timescale, erf.shoc.coeff_km supplies \(C_{K_m}\), and erf.shoc.coeff_kh supplies \(C_{K_h}\). Native SHOC also contains a special stable-mixing branch with fixed coefficients, so these two runtime coefficients do not control every possible diffusivity calculation.

The full input table is given in Inputs. For normal use, begin with the documented defaults and change closure coefficients only for a deliberate sensitivity study.

  • erf.shoc.lambda_low, lambda_high, lambda_slope, and lambda_thresh control the stability-dependent coefficient in the isotropy-timescale calculation.

  • erf.shoc.length_fac scales the diagnosed mixing length inversely: larger values reduce the otherwise diagnosed length.

  • The Native mixing-length expression is limited to the range 20 m–20 km before an additional horizontal-grid cap of \(\sqrt{\Delta x\,\Delta y}\) is applied. On sufficiently fine horizontal grids that final grid cap can therefore be smaller than 20 m.

  • erf.shoc.thl2tune, qw2tune, qwthl2tune, and w2tune tune the diagnosed second moments.

  • erf.shoc.c_diag_3rd_mom is the damping coefficient used in the diagnostic vertical-velocity third-moment closure.

Upper-boundary taper

The native upper taper is disabled by default:

erf.shoc.top_taper_depth = 0.0

For a positive taper depth \(d_{\mathrm{taper}}\), define

\[x = \operatorname{clip} \left( \frac{z_{\mathrm{top}}-z}{d_{\mathrm{taper}}}, 0, 1 \right).\]

Native SHOC uses the smooth multiplier

\[f = f_{\min} + (1-f_{\min})x^2(3-2x),\]

where \(f_{\min}\) is erf.shoc.top_taper_min_factor. Thus the multiplier is top_taper_min_factor at the model top and smoothly reaches one at the bottom of the taper layer.

The taper acts on native TKE evolution and TKE-budget quantities, isotropy and eddy diffusivities, and diagnosed higher-order moments. It is therefore a numerical/physical upper-boundary control, not only a plot-diagnostic filter.

Diagnostics

Native SHOC diagnostics are produced after the native driver runs in the supported state_update transport path. The diagnosed diffusivities remain available as diagnostics even though they are no longer exported for generic host reapplication.

Request 3-D fields through the ordinary plotfile variable list. For example:

erf.plot_vars_1 = density theta KE qv qc Kmv Khv Lturb pblh shoc_cldfrac shoc_ql wthv_sec brunt buoy_prod diss_tke

The standard turbulence fields Kmv and Khv are the density-weighted vertical momentum and heat diffusivities. For native SHOC these are \(\rho K_m\) and \(\rho K_h\), respectively. Lturb is the native SHOC mixing length, and nut is the kinematic momentum diffusivity obtained from Kmv / density.

Native SHOC also exposes the following 3-D diagnostics:

  • pblh

  • 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

See 3D plotfile variable reference for the exact definitions and units.

erf.shoc.extra_shoc_diags controls additional SHOC diagnostic calculation, not plot-variable selection. In the native implementation, shoc_cond is zero when this option is false. When the option is true, shoc_cond is diagnosed from positive PDF condensate change divided by the SHOC time step. A field must still be requested with erf.plot_vars_1 or erf.plot_vars_2 to appear in a plotfile.

When native SHOC diagnostics are not available, the native-only 3-D SHOC fields are written with the standard -999 missing value.

The 2-D diagnostic catalog also provides native SHOC surface/PBL diagnostics including pblh, shoc_u_star, shoc_Olen, and shoc_wthv_sfc. See 2D plotfile diagnostic reference for selection and missing-value behavior.

Advanced debugging

erf.shoc.debug_summary = true prints a synchronized min/max summary after each native SHOC advance. This can add synchronization and output overhead and is intended for debugging rather than routine production runs.

Only options documented in the native SHOC input table are part of the supported native user-facing runtime contract. Some historical/development erf.shoc keys remain parsed internally but do not currently alter native SHOC and should not be treated as native model features.

EAMxx SHOC notes

EAMXX_SHOC is a separate optional implementation path. It uses the shared SHOC closure-tuning input names documented in Inputs, but its data structures and coupling are provided through the EAMxx interface.

The EAMxx interface additionally supports:

  • erf.shoc.apply_tms to apply its turbulent mountain stress contribution.

  • erf.shoc.check_flux_state to run its flux/state consistency check.

These are EAMxx-interface controls. They are not active native-SHOC controls.

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 (a column whose bulk Richardson number never reaches the critical value takes the top of the scan, see the YSU PBL height scan). 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.

This mixed-layer profile is used below the PBL top in every stability regime, as in WRF (module_bl_mrf.F has no regime test below KPBL); stability enters only through \(w_s\). A Richardson-number closure applied inside the PBL instead makes the diffusivity follow the local shear, which excites a grid-scale, sign-flipping wind mode under the implicit vertical solve that the A-stable integrators (MidPoint, RK3) do not damp.

Note

Using the profile in every regime also inherits WRF MRF’s over-mixing of stable boundary layers: on GABLS1 the peak \(K_m\) is about twice that of YSU. This is the behaviour that Hong et al. (2006) addressed in YSU, so stable cases that need the weaker mixing should set erf.pbl_type = YSUNew.

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:

  1. HGAMT is already limited to GAMCRT (max sensible heating)

  2. Adding positive HGAMQ can legitimately produce VPERT > GAMCRT

  3. The latent heating from evaporation in unstable conditions is being artificially suppressed

  4. 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

6. Scale-Aware PBL-LES Blending

Motivation: In the grey-zone grid spacing (100–1000 m), mesoscale-calibrated PBL schemes like MRF overestimate vertical diffusivity, leading to excessively strong mixing. Two corrections are applied together to smoothly reduce K_h at fine resolution while maintaining full PBL behaviour at coarse scales.

Blending Factor (Boutle et al. 2014):

The Boutle blending factor smoothly interpolates from full PBL diffusivity at mesoscale to minimal PBL contribution at LES scales:

\[f_{blend} = \frac{\left(\frac{dx}{L_{blend}}\right)^2}{1 + \left(\frac{dx}{L_{blend}}\right)^2}\]
where:

  • \(dx\) is horizontal grid spacing [m]

  • \(L_{blend}\) is the blending length scale [m]. Typical: 750 m.

  • When \(dx >> L_{blend}\) (mesoscale): \(f_{blend} \to 1\) (full PBL)

  • When \(dx << L_{blend}\) (LES): \(f_{blend} \to 0\) (no PBL contribution)

\(f_{blend}\) is the weight on the parameterized part of the turbulence, so it increases with \(dx\): the coarser the grid, the less of the boundary layer overturning is resolved and the more the PBL scheme has to supply.

K_h Ceiling (Smagorinsky or Power-Law Fallback):

A second constraint prevents K_h from exceeding the subgrid-scale dissipation capacity:

  1. Smagorinsky ceiling (when strain rate available):

    \[K_{h,max} = (C_s \cdot dx)^2 \sqrt{2 \, S_{mn}S_{mn}}\]

    where \(C_s = 0.17\) (default) and \(S_{mn}S_{mn}\) is the strain rate squared [s⁻²].

  2. Power-law fallback (when strain rate unavailable, used in MRF):

    \[K_{h,max} = c_{max} \cdot dx^{4/3}\]

    where \(c_{max} = 0.1\) [m^(2/3) s⁻¹] (default).

Combined Effect:

\[K_{h,eff} = \min\left(K_h \cdot f_{blend}, K_{h,max}\right)\]

Applied only to vertical diffusivity for heat (\(K_{h,Theta}\)) and moisture (\(K_{h,Q}\)). Momentum diffusivity (\(K_{h,Mom}\)) is not modified, maintaining stratification-dependent mixing ratios.

Parameters:

  • pbl_blend_length (Real): Boutle blending length [m]. Default: 0.0 (disabled). Set > 0 to enable blending. Typical: 750 m for grey-zone applications.

  • pbl_blend_cs (Real): Smagorinsky coefficient for K_h ceiling. Default: 0.17.

  • pbl_blend_c_max (Real): Power-law ceiling coefficient [m^(2/3) s⁻¹]. Default: 0.1.

  • pbl_blend_use_smag (bool): Use Smagorinsky ceiling when strain rate available. Default: true. Set false to always use power-law.

Gating: All blending is a strict no-op when pbl_blend_length <= 0.0 (default). This ensures backward compatibility and allows selective enabling per simulation.

Proof-of-Concept Implementation: Currently applied to MRF scheme. The blending functions are defined in header Source/PBL/ERF_PBLScaleAwareBlending.H as GPU-inline free functions with no PBL-specific dependencies, enabling reuse across YSU, MYNN, or future schemes without code duplication.

Regression Tests:

Two ABL tests in Exec/CanonicalTests/ABL/MRF_Enhancements/ verify the implementation:

  1. blending_disabled: Confirms pbl_blend_length=0 is a strict no-op (identical to standard MRF).

  2. blending_active: Verifies blending reduces K_h by expected factor (0.2 at dx=375 m, L=750 m).

Physical References:

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.

On the grid, YSUNew finds \(h\) by scanning the bulk Richardson number of each column upward from its first fluid cell \(k_s\) to the top cell \(k_{hi}\), as WRF does (bl_ysu.F90):

  • If \({\rm Rib}\) first reaches \({\rm Rib_{cr}}\) at a level \(k > k_s\), \(h\) is interpolated linearly in \({\rm Rib}\) between \(z_{k-1}\) and \(z_k\). This includes a crossing at the top cell \(k_{hi}\).

  • The first cell’s \({\rm Rib}\) only seeds the lower end of the first interval. If it and the next level are both at or above \({\rm Rib_{cr}}\), \(h\) is the floor \(h_{min}\) (WRF’s brint = 0): the height of the first cell’s centre on a flat mesh, or of its top face on a terrain-fitted mesh, and at least 10 m.

  • If \({\rm Rib}\) stays below \({\rm Rib_{cr}}\) up to \(k_{hi}\), \(h\) is the top of the scan, \(z_{k_{hi}}\) (WRF’s brint = 1).

The result is then kept in \([h_{min},\ h_{max}]\) with \(h_{max} = 0.9\, z_{k_{hi}}\). The PBL top index is \(k_{hi}\) in the no-crossing case. The K profile then reaches \(h\), and the cells between \(h\) and \(k_{hi}\) keep the minimum diffusivity. No entrainment flux is added at the top cell of such a column.

MRF uses the same treatment for a column with no crossing, in its first guess, its corrector and its zero-Ri extent (WRF module_bl_mrf.F). Its PBL top index is then \(k_{hi}\) as well, so the K profile and, with erf.enable_mrf_countergradient, the countergradient terms reach \(h\). When the first two levels are both at or above \({\rm Rib_{cr}}\), MRF takes the first level’s height, as WRF does, rather than extrapolating above it.

Before erf-model/ERF#4458, both schemes gave a column with no crossing \(h_{min}\), and YSUNew also sent a crossing at \(k_{hi}\) there. MRF also set its index to \(k_s + 1\), so only the first cell used the K profile.

In YSUNew the sign of the surface virtual heat flux selects the form of the diffusivity inside the PBL. A column counts as stable when the flux is negative and as unstable or neutral otherwise, as WRF does with its surface bulk Richardson number (sfcflg). Over a neutral surface the flux is round-off. With a plane-averaged surface layer (the default erf.most.average_policy) it comes from an averaged \(\theta\) that differs from the surface value only in the last bits, and how depends on the box and tile decomposition. A flux within \(64\,\epsilon\, u_* \theta_1\) of zero therefore counts as neutral, where \(\epsilon\) is the machine epsilon of amrex::Real and \(\theta_1\) the dry potential temperature of the lowest cell, so that the choice does not depend on the decomposition. With \(u_* = 0.3\) m/s and \(\theta_1 = 300\) K the band is about \(1.3\times 10^{-12}\) K m/s in double precision and \(7\times 10^{-4}\) K m/s (about 0.8 W/m2) in single precision.

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:

YSUNew Cloudy PBL Extensions

The YSUNew implementation includes extensions for handling cloud-topped boundary layers without requiring radiation coupling. These extensions are controlled by the enable_ysu_topdown flag and the use_moisture condition.

Liquid Potential Temperature PBLH Extension

When enabled, an additional upward scan of the PBL height is performed using liquid potential temperature (\(\theta_{li}\)) as the stability criterion. This extension allows detection of deeper cloud-topped boundary layers where the traditional bulk Richardson number criterion may underestimate the mixed layer depth. The scan uses an unstable threshold (zero Richardson number) to extend the PBL height upward through layers where cloud liquid water provides buoyancy support.

Configuration: enable_ysu_topdown = true, use_moisture = true

Cloudy Entrainment Correction

When liquid water content plus ice content at the layer below PBL top exceeds a threshold (0.01 g/kg), an adjusted entrainment coefficient is computed using cloud buoyancy considerations. This correction replaces the clear-sky entrainment velocity with a value derived from the liquid-theta buoyancy jump, accounting for the reduced stability at cloud top. The entrainment efficiency is computed from cloud liquid water content and limited to 0.4. This path does not require radiation coupling; the calculation uses only surface buoyancy flux.

Configuration: enable_ysu_topdown = true, use_moisture = true