(theory_abl_stratified)= # The thermally-stratified atmospheric boundary layer The {ref}`turbulent-inflow chapter ` treats the atmospheric boundary layer (ABL) as mechanically neutral: the mean wind is the logarithmic profile set by friction and roughness, and the turbulence is whatever the {ref}`SEM ` or a {ref}`precursor ` supplies. That picture holds only when the ground and the air are at the same temperature. On a sunny day the ground is warmer than the air above it and buoyancy *adds* turbulence; on a clear night the ground is colder and buoyancy *suppresses* it. This page restores that thermodynamic degree of freedom to the ABL, through the thermal route already established in the {ref}`thermodynamics chapter `: a buoyancy-carrying **active scalar** coupled back to momentum by a single Boussinesq body force. It underpins the stratified-ABL and urban pollutant-dispersion validation cases. ## 1. Motivation and stratification regimes Wind loads, pedestrian comfort and, above all, pollutant dispersion in cities are set by the turbulence in the lowest tens of metres of the atmosphere, and that turbulence is not governed by shear alone. Whenever the ground and the near-surface air differ in temperature, a vertical density gradient develops and gravity does work on vertical motions. This buoyant production of turbulent kinetic energy (TKE) adds to, or subtracts from, the shear production, and its sign fixes the regime {footcite:t}`turner1994workbook`: - **Unstable (convective)** stratification. The surface is warmer than the air, so a parcel displaced upward is lighter than its surroundings and keeps rising: buoyancy *produces* TKE, thickens the mixed layer and intensifies vertical exchange. This is the daytime regime, and the one both benchmarks below target {footcite:t}`jiang2018large,zhou2021large`. - **Neutral** stratification. Surface and air temperatures match, buoyancy neither produces nor destroys TKE, and the ABL reduces to the purely mechanical logarithmic layer of the {ref}`inflow chapter `. - **Stable** stratification. The surface is colder than the air, so a displaced parcel is heavier than its surroundings and sinks back: buoyancy *consumes* TKE, damps vertical motion and can laminarise the near-surface flow. This is the night-time regime, the hardest to reproduce and the one most associated with pollutant trapping. ```{figure} /_static/img/theory/stratified_abl/stratification_regimes.svg --- width: 80% align: center --- The sign of near-surface buoyancy selects the stratification regime: a warm floor (unstable) produces TKE and mixes momentum into a fuller, less-sheared wind profile, equal temperatures (neutral) leave the plain logarithmic profile, and a cold floor (stable) consumes TKE and leaves a more strongly sheared profile, matching the sign of the buoyant production $P_b$ and of the gradient Richardson number. ``` ```{admonition} Two numbers name the regime --- class: note --- The **gradient Richardson number** compares the buoyant destruction of TKE to its shear production, using local gradients, $$ \mathrm{Ri} = \frac{\dfrac{g}{\theta_0}\,\dfrac{\partial \bar\theta}{\partial z}} {\left(\dfrac{\partial \bar u}{\partial z}\right)^{2}}, $$ (abl_gradient_ri) with $\theta_0$ a reference (potential) temperature and $g$ the gravitational acceleration. $\mathrm{Ri} < 0$ is unstable, $\mathrm{Ri} = 0$ neutral, $\mathrm{Ri} > 0$ stable. Because $\partial\bar u/\partial z$ and $\partial\bar\theta/\partial z$ are awkward to evaluate at a wall, wind-tunnel and CFD studies more often report a **bulk Richardson number** built from finite differences across a reference height. The two benchmarks below each define their own, so bulk-Richardson values are not directly comparable between studies. ``` The **Monin-Obukhov** framing of Section 2 makes the same statement with a length rather than a ratio: the Obukhov length $L$ is the height at which buoyant and shear production become comparable, so $z/L$ plays the role of $\mathrm{Ri}$ through the surface layer. The classical **Pasquill-Gifford** stability classes A-F used in regulatory Gaussian-plume dispersion {footcite:t}`turner1994workbook` are a coarse binning of exactly this axis, from strongly convective (A) through neutral (D) to strongly stable (F). ## 2. Monin-Obukhov similarity theory Monin-Obukhov similarity theory (MOST) extends the neutral law of the wall ({ref}`derived here `) to a stratified surface layer. The neutral argument had a single velocity scale $u^{*}$ and a single length $z$; MOST adds a second length, the **Obukhov length** $$ L = -\,\frac{u^{*3}\,\theta_0}{\kappa\, g\, \overline{w'\theta'}_0}, $$ (abl_obukhov_length) where $\overline{w'\theta'}_0$ is the surface kinematic heat flux, $\kappa \approx 0.41$ the von Karman constant and $\theta_0$ a reference potential temperature. An upward (warming) surface flux $\overline{w'\theta'}_0 > 0$ makes the length negative, $L < 0$, which is the unstable regime; $L > 0$ is stable and $|L| \to \infty$ is neutral. The dimensionless stability parameter is $\zeta = z/L$. This sign convention is used consistently throughout the page. Nassu carries the buoyancy through the Boussinesq scalar of Section 3, where the buoyant forcing per unit scalar is $\beta\,g$ rather than $g/\theta_0$. The wall model therefore uses the buoyancy strength $B \equiv \beta\,|\mathbf g|$ read from the coupling scalar's Boussinesq block in place of $g/\theta_0$, so {eq}`abl_obukhov_length` becomes $L = -u^{*3}/(\kappa\,B\,\overline{w'\theta'}_0)$. Sharing $\beta g$ with the fluid body force guarantees the wall model and the resolved buoyancy cannot disagree on the stratification strength. Dimensional analysis then fixes the **stability-corrected** gradients. The mean-wind and mean-temperature gradients become $$ \frac{\partial \bar u}{\partial z} = \frac{u^{*}}{\kappa z}\,\phi_m\!\left(\zeta\right), \qquad \frac{\partial \bar\theta}{\partial z} = \frac{\theta^{*}}{\kappa z}\,\phi_h\!\left(\zeta\right), \qquad \theta^{*} \equiv -\,\frac{\overline{w'\theta'}_0}{u^{*}}, $$ (abl_most_gradients) with $\phi_m, \phi_h$ the universal stability functions and $\theta^{*}$ the temperature scale. Integrating from the roughness length $z_0$ upward gives the **stability-corrected log law** that generalises the neutral profile, $$ \bar u(z) = \frac{u^{*}}{\kappa}\left[\ln\!\left(\frac{z}{z_0}\right) - \psi_m\!\left(\frac{z}{L}\right)\right], \qquad \bar\theta(z) - \theta_s = \frac{\theta^{*}}{\kappa}\left[\ln\!\left(\frac{z}{z_{0h}}\right) - \psi_h\!\left(\frac{z}{L}\right)\right], $$ (abl_most_loglaw) where $\theta_s$ is the aerodynamic surface temperature, $z_{0h}$ the roughness length for heat, and $\psi_m, \psi_h$ the integrated stability corrections, $\psi(\zeta) = \int_0^{\zeta}[1-\phi(\zeta')]\,\mathrm{d}\zeta'/\zeta'$. The heat roughness $z_{0h}$ is generally smaller than the momentum $z_0$; their ratio is an empirical input that a validation case must fix from a cited source rather than assume equal. Setting $\psi_m = \psi_h = 0$ recovers the neutral profile of the {ref}`inflow chapter ` exactly, so the stratified law is a strict extension of the neutral one. ```{figure} /_static/img/theory/stratified_abl/mo_loglaw_profiles.svg --- width: 65% align: center --- Stability-corrected log-law mean-wind profiles for one $u_*$ and $z_0$: the unstable (convective) profile is less sheared and well mixed, the stable (night-time) profile more sheared, and neutral recovers the plain log law when $\psi_m = 0$. This is the profile ordering Nassu's GPU differential test reproduces. ``` The universal functions are empirical. The widely used **Businger-Dyer / Paulson** forms are, for the **unstable** branch ($\zeta < 0$), with $x = (1 - \gamma\,\zeta)^{1/4}$, $$ \psi_m(\zeta) = 2\ln\!\frac{1+x}{2} + \ln\!\frac{1+x^{2}}{2} - 2\arctan x + \frac{\pi}{2}, \qquad \psi_h(\zeta) = 2\ln\!\frac{1+x^{2}}{2}, $$ (abl_psi_unstable) and for the **stable** branch ($\zeta > 0$) the linear form $\psi_m(\zeta) = \psi_h(\zeta) = -\beta_s\,\zeta$. The coefficients $\gamma$ and $\beta_s$ are empirical and source-dependent (commonly $\gamma \approx 16$, $\beta_s \approx 5$, but $\gamma = 15$ and $\beta_s = 4.7$ are also in use), so each validation case must cite the specific values it adopts. ```{figure} /_static/img/theory/stratified_abl/mo_stability_functions.svg --- width: 65% align: center --- The integrated stability corrections $\psi_m$ (momentum) and $\psi_h$ (heat): both vanish at neutral ($\zeta = 0$), rise positive on the unstable side where they subtract from $\ln(z/z_0)$ to fill out the mean profile, and fall along the common linear branch $-\beta_s\,\zeta$ on the stable side. Curves are evaluated from the same `nassu.IBM.mo_stability` code the wall model uses. ``` The surface-layer relation between the local $\mathrm{Ri}$ of {eq}`abl_gradient_ri` and $\zeta = z/L$ follows from {eq}`abl_most_gradients`, $\mathrm{Ri} = \zeta\,\phi_h(\zeta)/\phi_m^{2}(\zeta)$, so $\mathrm{Ri} \to \zeta$ in the neutral limit. This is the tie between the two framings of Section 1. ```{admonition} Why this grounds the IBM wall-model correction --- class: important --- Equation {eq}`abl_most_loglaw` is the stratified generalisation of the log law the {ref}`IBM wall model ` imposes at building and ground surfaces. Under stratification the friction velocity recovered from the matching-point tangential velocity acquires the $\psi_m(z/L)$ correction. Nassu drives this **flux-driven**: the surface kinematic heat flux $\overline{w'\theta'}_0$ is prescribed directly and forms $L$ through {eq}`abl_obukhov_length`, so no surface-temperature or temperature-scale inversion is needed and the momentum wall model never samples the scalar field. Because $u^{*}$ and $L$ are coupled through {eq}`abl_obukhov_length` and {eq}`abl_most_loglaw`, they are converged together by a short fixed-point (Picard) iteration seeded from the neutral $u^{*}$. The argument $\zeta = z/L$ is clamped to the surface-layer validity range: besides respecting MOST validity this breaks the positive feedback of the stable-branch $u^{*}$ iteration, which would otherwise drive $u^{*}\to 0$, $L\to 0$ and $\zeta\to\infty$. ``` ```{admonition} Implementation status --- class: note --- The Monin-Obukhov correction of this section and the height-varying (stratified) scalar inlet are **implemented** in epic #1034 P4 (the IBM `EqLog` wall model applies the $\psi_m(z/L)$-corrected log law, and equation-form scalar inlet / wall values are baked per node into the scalar BC kernels). The correction is **flux-driven** (fix #1043): the surface kinematic heat flux is prescribed directly through `surface_heat_flux` and $u^{*}$ and $L$ are converged by a short Picard iteration. The stability-function coefficients are hard-coded to $\gamma = 16$, $\beta_s = 5$; as noted above, a validation case must cite the coefficients it adopts. The wall model has been **GPU-validated** by a differential neutral / unstable / stable sign check (epic #1034 P5). ``` ## 3. Boussinesq scalar buoyancy in the LBM Nassu carries stratification through the **constant-density Boussinesq** route of the {ref}`modeling hierarchy `. A single **active scalar** $b$ - a buoyancy variable proportional to the potential-temperature perturbation, $b \equiv \bar\theta - \theta_{\text{ref}}$ - is transported on the existing {doc}`double-distribution scalar path `, obeying the plain advection-diffusion equation $$ \partial_t b + u_\alpha \partial_\alpha b = D\,\nabla^2 b + S_b, $$ (abl_scalar_ade) with $D$ the diffusivity (molecular plus the {doc}`subgrid ` contribution when LES is active) and $S_b$ any surface or volumetric source. The scalar is *active*, not passive, because it feeds back to the momentum equation through one term only, the **Boussinesq body force** $$ F^{\text{buoy}}_\alpha = \rho_0\,\beta\,b\,g\,\hat e_z, $$ (abl_buoyancy_force) with $\beta$ the thermal expansion coefficient, $\rho_0$ the constant reference density and $\hat e_z$ the vertical unit vector, so a warm parcel ($b>0$) is pushed upward. The force {eq}`abl_buoyancy_force` enters the fluid **not** by modifying the equilibrium but through the **Guo forcing scheme** as the body-force density $g_\alpha$. It contributes the half-step velocity correction $\rho u_\alpha = \sum_i f_i c_{i,\alpha} + \tfrac{\Delta t}{2} F^{\text{buoy}}_\alpha$ and the Guo source populations documented on the {ref}`lattice-unit page `. The fluid collision operator itself is untouched: the coupling is one additive force term, evaluated from the current scalar field each step. ```{figure} /_static/img/theory/stratified_abl/boussinesq_scalar_coupling.svg --- width: 65% align: center --- The fluid velocity advects both scalars on one shared double-distribution path, but only the active buoyancy scalar $b$ feeds back to momentum, through the single Boussinesq body force; the pollutant is passive with no back-coupling. ``` ```{admonition} Why the scalar-Boussinesq route, not an energy equation or the EOS --- class: important --- Three thermal routes exist in the {ref}`hierarchy `; the ABL uses the weakest deliberately. - **The perturbations are small.** A daytime near-surface air-ground difference of a few to a few tens of kelvin gives $\beta\,\Delta T \sim 10^{-2}$ to $10^{-1}$, well inside the Boussinesq small-$\Delta T$ envelope $\beta\,\Delta T \ll 1$. Reproducing the *density* field to fire-scale accuracy is unnecessary; reproducing the buoyant *forcing* is the whole physics, and that is exactly what {eq}`abl_buoyancy_force` supplies. - **It is more numerically stable.** The Boussinesq route keeps the fluid zeroth moment as the density and touches the collision only through an additive force. The variable-density low-Mach closure instead reinterprets the zeroth moment as reduced *pressure* and slaves density to temperature through the equation of state, which admits the checkerboard pressure mode analysed on the {ref}`stability page `. The ABL does not need that machinery. - **It reuses one code path.** $b$ rides the same scalar transport, multiblock {doc}`flux rescaling ` and {doc}`LES subgrid closure ` ($D_{\text{SGS}} = \nu_{\text{SGS}}/Sc_t$) as any tracer. There is no energy solver and no equation of state to maintain. ``` ```{admonition} Surface boundaries: IBM for momentum, voxelization for the scalar --- class: note --- The building and ground momentum boundaries are imposed by the {ref}`immersed boundary method ` (optionally with the stratified wall-model correction of Section 2), while the **surface heat/scalar condition** on $b$ - a prescribed wall value or a prescribed wall-normal flux - is imposed by {doc}`voxelization ` on the near-surface band rather than by an equilibrium wall-model BC. A prescribed surface heat flux is the natural driver here, because it sets $\overline{w'\theta'}_0$ and therefore the Obukhov length {eq}`abl_obukhov_length` directly. ``` The mean part of the inflow is the stratified profile {eq}`abl_most_loglaw`, and the fluctuations come from the {ref}`SEM ` or a {ref}`precursor ` as in the neutral case; the scalar carries a matching mean temperature profile and, where available, its own fluctuation statistics. As the {ref}`inflow chapter ` stresses, the inlet and the ground are one boundary condition: under stratification they must agree not only on $z_0$ but on the surface heat flux, or the boundary layer re-equilibrates over the fetch. ## 4. Dispersion in the stratified ABL Pollutant transport rides the same scalar machinery: a **passive** tracer $c$ obeys {eq}`abl_scalar_ade` with its own source, decoupled from momentum (no buoyancy force). Field studies and wind-tunnel benchmarks report the **normalised concentration** $$ K \equiv \frac{\langle c\rangle\, U_H\, H^{2}}{q} = \frac{\langle c\rangle}{C_0}, \qquad C_0 = \frac{q}{U_H H^{2}}, $$ (abl_normalized_conc) with $q$ the volumetric emission rate, $U_H$ the reference velocity at building height $H$, and $C_0$ the reference concentration; this is the normalisation both benchmarks use {footcite:t}`jiang2018large,zhou2021large`. The mechanism by which stratification changes dispersion is read most cleanly from the **flux decomposition** of the mean scalar transport across a surface, $$ Q_\alpha = \underbrace{\langle u_\alpha\rangle\langle c\rangle}_{\text{advection } Q_{a}} + \underbrace{\langle u_\alpha' c'\rangle}_{\text{turbulent } Q_{t}} - \underbrace{D\,\partial_\alpha\langle c\rangle}_{\text{molecular } Q_{m}}, $$ (abl_flux_decomposition) as used by {footcite:t}`jiang2018large`, {footcite:t}`zhou2021large` and {footcite:t}`bazdiditehrani2019large`. Three robust, benchmark-supported effects of increasing instability follow. - **Enhanced vertical mixing lowers near-ground concentration.** Buoyant TKE production raises the vertical turbulent flux $Q_t$, venting the tracer out of street canyons and wakes and reducing $K$ near the source. {footcite:t}`jiang2018large` show for a 3-D street canyon that the roof-level exchange is turbulence-dominated. - **Wakes and recirculation shrink.** {footcite:t}`zhou2021large` report that for an isolated rectangular building the streamwise length of the leeward recirculation region shrinks by $32\%$ as $\mathrm{Ri}$ goes from $0$ to $-1.5$, and that inside that region the *advection* flux $Q_a$ reaches roughly four times the turbulent-diffusion flux $Q_t$ at $\mathrm{Ri} = -1.5$. The mean flow, not turbulent diffusion, therefore controls transport in the near wake even as instability grows. - **Turbulent transport can run counter to the mean gradient.** {footcite:t}`bazdiditehrani2019large`, decomposing the convective and turbulent diffusion fluxes around a high-rise building across the stable, neutral and unstable regimes, show that LES reproduces the *counter-gradient* mechanism - a turbulent flux directed up the mean concentration gradient - which the gradient-diffusion closures of RANS cannot represent. This is a further reason the buoyancy-driven dispersion cases here are run with LES rather than RANS. ```{admonition} Where this sits in the CFD-dispersion landscape --- class: note --- Urban pollutant dispersion is now routinely tackled with CFD, but the systematic review of {footcite:t}`pantusheva2022air` finds that many published studies omit model-setup details needed to reproduce them and that few carry out formal verification and validation. LES outperforms RANS for the unsteady, buoyancy-driven transport of these flows, but it is far more sensitive to the **inflow turbulence**: the tracer field responds directly to the resolved eddies, so a poorly conditioned inlet corrupts the concentration field even when the mean wind is correct {footcite:t}`ghanetehrani2025comparative` - an effect first shown for a non-isothermal high-rise building by {footcite:t}`yoshie2011cfd`, whose LES matched the wind-tunnel gas dispersion only once inflow turbulence was supplied. That sensitivity, together with the need for a *divergence-free* stratified inlet, is why the ABL-dispersion cases lean on the well-conditioned synthetic and spectral inflow generators - digital-filter {footcite:t}`klein2003digital`, divergence-free synthetic-eddy {footcite:t}`poletto2013divergence,kim2013divergence`, spectral {footcite:t}`melaku2021divergence` and optimised-ABL {footcite:t}`lamberti2018optimizing` methods - reviewed in the {ref}`inflow chapter `, and validate against the stratified wind-tunnel benchmarks of {footcite:t}`jiang2018large` and {footcite:t}`zhou2021large` (see also the accidental-release review {footcite:t}`dou2022atmospheric`). ``` ```{eval-rst} .. footbibliography:: ```