Obukhov

class abcmodel.atmos.surface_layer.obukhov.ObukhovState(ustar, z0m, z0h, drag_m=<factory>, drag_s=<factory>, uw=<factory>, vw=<factory>, temp_2m=<factory>, q2m=<factory>, u2m=<factory>, v2m=<factory>, e2m=<factory>, esat2m=<factory>, thetasurf=<factory>, thetavsurf=<factory>, qsurf=<factory>, obukhov_length=<factory>, rib_number=<factory>, ra=<factory>)

Bases: AbstractSurfaceLayerState

Standard surface layer model state.

ustar: Array

Surface friction velocity [m/s].

z0m: Array

Roughness length for momentum [m].

z0h: Array

Roughness length for scalars [m].

drag_m: Array

Drag coefficient for momentum [-].

drag_s: Array

Drag coefficient for scalars [-].

uw: Array

Surface momentum flux u [m2 s-2].

vw: Array

Surface momentum flux v [m2 s-2].

temp_2m: Array

2m temperature [K].

q2m: Array

2m specific humidity [kg kg-1].

u2m: Array

2m u-wind [m s-1].

v2m: Array

2m v-wind [m s-1].

e2m: Array

2m vapor pressure [Pa].

esat2m: Array

2m saturated vapor pressure [Pa].

thetasurf: Array

Surface potential temperature [K].

thetavsurf: Array

Surface virtual potential temperature [K].

qsurf: Array

Surface specific humidity [kg kg-1].

obukhov_length: Array

Obukhov length [m].

rib_number: Array

Bulk Richardson number [-].

ra: Array

Aerodynamic resistance [s/m].

class abcmodel.atmos.surface_layer.obukhov.ObukhovModel(*args: Any, **kwargs: Any)

Bases: AbstractSurfaceLayerModel

Standard surface layer model with atmospheric stability corrections.

Calculates surface-atmosphere exchange using Monin-Obukhov similarity theory with stability functions and iterative solution for Obukhov length.

init_state(ustar=0.3, z0m=0.02, z0h=0.002)

Initialize the model state.

Parameters:

• ustar (float) – Friction velocity [m s-1]. Default is 0.3.

• z0m (float) – Surface roughness length for momentum [m]. Default is 0.02.

• z0h (float) – Surface roughness length for heat [m]. Default is 0.002.

Returns:

The initial surface layer state.

run(state)

Run the model.

Parameters:

state (AbstractCoupledState[TypeVar(RadT, bound= AbstractRadiationState), TypeVar(LandT, bound= AbstractLandState), DayOnlyAtmosphereState[ObukhovState, TypeVar(MixedT, bound= AbstractMixedLayerState), TypeVar(CloudT, bound= AbstractCloudState)]])

Returns:

The updated state.

compute_ra(u, v, wstar, drag_s)

Calculate aerodynamic resistance from wind speed and drag coefficient.

Parameters:

• u (Array) – zonal wind speed [m s-1].

• v (Array) – meridional wind speed [m s-1].

• wstar (Array) – convective velocity scale [m s-1].

• drag_s (Array) – drag coefficient for scalars [-].

Returns:

Aerodynamic resistance [s m-1].

Notes

The aerodynamic resistance is given by

\[r_a = \frac{1}{C_s u_{\text{eff}}}\]

where \(C_s\) is the drag coefficient for scalars and \(u_{\text{eff}}\) is the effective wind speed.

compute_effective_wind_speed(u, v, wstar)

Compute effective wind speed ueff.

Parameters:

• u (Array) – zonal wind speed \(u\).

• v (Array) – meridional wind speed \(v\).

• wstar (Array) – convective velocity scale \(w_*\).

Returns:

Effective wind speed \(u_{\text{eff}}\).

Notes

The effective wind speed is given by

\[u_{\text{eff}} = \sqrt{u^2 + v^2 + w_*^2}.\]

A minimum value of 0.01 m/s is enforced to avoid division by zero afterwards.

compute_zsl(h_abl)

Compute surface layer height.

Parameters:

h_abl (Array) – Atmospheric boundary layer height [m].

Returns:

Surface layer height [m], defined as 10% of ABL height.

Notes

The surface layer height is conventionally taken as 10% of the atmospheric boundary layer height.

compute_thetasurf(theta, wtheta, drag_s, ueff)

Compute surface potential temperature.

Parameters:

• theta (Array) – mixed layer potential temperature \(\theta\).

• wtheta (Array) – surface kinematic heat flux \(w'\theta'\).

• drag_s (Array) – surface drag coefficient \(C_s\).

• ueff (Array) – effective wind speed \(u_{\text{eff}}\).

Returns:

Surface potential temperature \(\theta_s\).

Notes

The surface potential temperature is given by

\[\theta_s = \theta + \frac{w'\theta'}{C_s u_{\text{eff}}}.\]

compute_qsurf(q, thetasurf, surf_pressure, rs, drag_s, ueff)

Compute surface specific humidity.

Parameters:

• q (Array) – mixed layer specific humidity \(q\).

• thetasurf (Array) – surface potential temperature \(\theta_s\).

• surf_pressure (Array) – surface pressure \(p\).

• rs (Array) – surface resistance \(r_s\).

• drag_s (Array) – surface drag coefficient \(C_s\).

• ueff (Array) – effective wind speed \(u_{\text{eff}}\).

Returns:

Surface specific humidity \(q_s\).

Notes

The surface specific humidity is a weighted average between the air and the saturated value at the surface

\[q_s = (1 - c_q) q + c_q q_{sat}(\theta_s, p),\]

where \(c_q = [1 + C_s u_{\text{eff}} r_s]^{-1}\) and \(q_{sat}\) is the saturation specific humidity.

compute_thetavsurf(thetasurf, qsurf)

Compute surface virtual potential temperature.

Parameters:

• thetasurf (Array) – surface potential temperature \(\theta_s\).

• qsurf (Array) – surface specific humidity \(q_s\).

Returns:

Surface virtual potential temperature \(\theta_{v,s}\).

Notes

The surface virtual potential temperature is

\[\theta_{v,s} = \theta_s (1 + 0.61 q_s)\]

compute_richardson_number(ueff, zsl, g, thetav, thetavsurf)

Compute bulk Richardson number.

Parameters:

• ueff (Array) – effective wind speed \(u_{\text{eff}}\).

• zsl (Array) – surface layer height \(z_{sl}\).

• g (float) – gravity \(g\).

• thetav (Array) – virtual potential temperature at reference height \(\theta_v\).

• thetavsurf (Array) – Surface virtual potential temperature \(\theta_{v,s}\).

Notes

The bulk Richardson number is given by

\[Ri_b = \frac{g}{\theta_v} \frac{z_{sl} (\theta_v - \theta_{v,s})}{u_{\text{eff}}^2}.\]

The value is capped at 0.2 for numerical stability.

compute_rib_function(zsl, oblen, rib_number, z0h, z0m)

Compute the Richardson number function for iterative solution of Obukhov length.

Notes

This function computes the difference between the bulk Richardson number and its Monin-Obukhov similarity theory estimate, used in the Newton-Raphson iteration for finding the Obukhov length.

The function is:

\[f(L) = Ri_b - \frac{z_{sl}}{L} \frac{\psi_h(z_{sl}/L) - \psi_h(z_{0h}/L) + \ln(z_{sl}/z_{0h})}{[\psi_m(z_{sl}/L) - \psi_m(z_{0m}/L) + \ln(z_{sl}/z_{0m})]^2}\]

where \(Ri_b\) is the bulk Richardson number, \(z_{sl}\) is the surface layer height, \(L\) is the Obukhov length, \(z_{0h}\) and \(z_{0m}\) are roughness lengths for scalars and momentum, and \(\psi_h\), \(\psi_m\) are stability correction functions.

ribtol(zsl, rib_number, z0h, z0m)

Iteratively solve for the Obukhov length given the Richardson number.

Notes

Uses a Newton-Raphson method to find the Obukhov length \(L\) such that the Monin-Obukhov similarity theory estimate matches the bulk Richardson number.

The iteration continues until the change in \(L\) is below a threshold or a maximum value is reached.

compute_drag_m(zsl, k, obukhov_length, z0m)

Compute drag coefficient for momentum with stability corrections.

Parameters:

• zsl (Array) – surface layer height \(z_{sl}\).

• k (float) – Von Kármán constant \(k\).

• obukhov_length (Array) – Obukhov length \(L\).

• z0m (Array) – roughness length for momentum \(z_{0m}\).

Returns:

Drag coefficient for momentum \(C_m\).

Notes

The drag coefficient for momentum is given by

\[C_m = \frac{k^2}{[\psi_m(z_{sl}/L) - \psi_m(z_{0m}/L) + \ln(z_{sl}/z_{0m})]^2}\]

where \(\psi_m\) (see compute_momentum_correction_term()) is the stability correction function for momentum.

compute_drag_s(zsl, k, obukhov_length, z0h, z0m)

Compute drag coefficient for scalars with stability corrections.

Parameters:

• zsl (Array) – surface layer height \(z_{sl}\).

• k (float) – Von Kármán constant \(k\).

• obukhov_length (Array) – Obukhov length \(L\).

• z0h (Array) – roughness length for scalars \(z_{0h}\).

• z0m (Array) – roughness length for momentum \(z_{0m}\).

Returns:

Drag coefficient for scalars \(C_s\).

Notes

The drag coefficient for scalars is given by

\[C_s = \frac{k^2}{[\psi_m(z_{sl}/L) - \psi_m(z_{0m}/L) + \ln(z_{sl}/z_{0m})] [\psi_h(z_{sl}/L) - \psi_h(z_{0h}/L) + \ln(z_{sl}/z_{0h})]}\]

where \(\psi_m\) (see compute_momentum_correction_term()) and \(\psi_h\) (see compute_scalar_correction_term()) are stability correction functions for momentum and scalars.

compute_ustar(ueff, drag_m)

Compute surface friction velocity.

Parameters:

• ueff (Array) – effective wind speed \(u_{\text{eff}}\).

• drag_m (Array) – drag coefficient for momentum \(C_m\).

Returns:

Friction velocity \(u_*\).

Notes

The friction velocity \(u_*\) is given by

\[u_* = \sqrt{C_m} u_{\text{eff}}.\]

compute_uw(ueff, u, drag_m)

Compute zonal momentum flux.

Parameters:

• ueff (Array) – effective wind speed \(u_{\text{eff}}\).

• u (Array) – zonal wind speed \(u\).

• drag_m (Array) – drag coefficient for momentum \(C_m\).

Returns:

Zonal momentum flux \(\overline{u'w'}\).

Notes

The zonal momentum flux is given by

\[\overline{u'w'} = -C_m u_{\text{eff}} u.\]

compute_vw(ueff, v, drag_m)

Compute meridional momentum flux.

Parameters:

• ueff (Array) – effective wind speed \(u_{\text{eff}}\).

• v (Array) – meridional wind speed \(v\).

• drag_m (Array) – drag coefficient for momentum \(C_m\).

Returns:

Meridional momentum flux \(\overline{v'w'}\).

Notes

The meridional momentum flux is given by

\[\overline{v'w'} = -C_m u_{\text{eff}} v.\]

compute_temp_2m(thetasurf, wtheta, ustar, k, z0h, obukhov_length)

Compute 2m temperature diagnostic.

Parameters:

• thetasurf (Array) – surface potential temperature \(\theta_s\).

• wtheta (Array) – surface kinematic heat flux \(w'\theta'\).

• ustar (Array) – friction velocity \(u_*\).

• k (float) – Von Kármán constant \(k\).

• z0h (Array) – roughness length for scalars \(z_{0h}\).

• obukhov_length (Array) – Obukhov length \(L\).

Returns:

Temperature at 2 meters above surface [K].

Notes

Uses Monin-Obukhov similarity theory with stability corrections to extrapolate surface temperature to 2m height.

compute_q2m(qsurf, wq, ustar, k, z0h, obukhov_length)

Compute 2m specific humidity diagnostic.

Parameters:

• qsurf (Array) – surface specific humidity \(q_s\).

• wq (Array) – surface kinematic moisture flux \(w'q'\).

• ustar (Array) – friction velocity \(u_*\).

• k (float) – Von Kármán constant \(k\).

• z0h (Array) – roughness length for scalars \(z_{0h}\).

• obukhov_length (Array) – Obukhov length \(L\).

Returns:

Specific humidity at 2 meters above surface [kg kg-1].

Notes

Uses Monin-Obukhov similarity theory with stability corrections to extrapolate surface specific humidity to 2m height.

compute_u2m(uw, ustar, k, z0m, obukhov_length)

Compute 2m zonal wind diagnostic.

Parameters:

• uw (Array) – zonal momentum flux \(\overline{u'w'}\).

• ustar (Array) – friction velocity \(u_*\).

• k (float) – Von Kármán constant \(k\).

• z0m (Array) – roughness length for momentum \(z_{0m}\).

• obukhov_length (Array) – Obukhov length \(L\).

Returns:

Zonal wind at 2 meters above surface [m s-1].

Notes

Uses Monin-Obukhov similarity theory with stability corrections to extrapolate surface momentum flux to 2m wind speed.

compute_v2m(vw, ustar, k, z0m, obukhov_length)

Compute 2m meridional wind diagnostic.

Parameters:

• vw (Array) – meridional momentum flux \(\overline{v'w'}\).

• ustar (Array) – friction velocity \(u_*\).

• k (float) – Von Kármán constant \(k\).

• z0m (Array) – roughness length for momentum \(z_{0m}\).

• obukhov_length (Array) – Obukhov length \(L\).

Returns:

Meridional wind at 2 meters above surface [m s-1].

Notes

Uses Monin-Obukhov similarity theory with stability corrections to extrapolate surface momentum flux to 2m wind speed.

compute_e2m(q2m, surf_pressure)

Compute 2m vapor pressure.

Parameters:

• q2m (Array) – specific humidity at 2m \(q_{2m}\).

• surf_pressure (Array) – surface pressure \(p\).

Returns:

Vapor pressure at 2 meters above surface [Pa].

Notes

Converts specific humidity to vapor pressure using:

\[e_{2m} = \frac{q_{2m} \cdot p}{0.622}\]

compute_esat2m(temp_2m)

Compute 2m saturated vapor pressure.

Parameters:

temp_2m (Array) – temperature at 2m [K].

Returns:

Saturated vapor pressure at 2 meters above surface [Pa].

Notes

Uses the Tetens formula:

\[e_{sat,2m} = 611 \exp\left(\frac{17.2694(T_{2m} - 273.16)}{T_{2m} - 35.86}\right)\]

compute_scalar_correction_term(z, oblen, z0h)

Compute scalar stability correction term.

Parameters:

• z (Array | float) – height above ground level \(z\).

• oblen (Array) – Obukhov length \(L\).

• z0h (Array) – roughness length for heat \(z_{0h}\).

Returns:

The scalar stability correction.

Notes

This term is used in Monin-Obukhov similarity theory for scalars, and is given by

\[\ln\left(\frac{z}{z_{0h}}\right) - \psi_h\left(\frac{z}{L}\right) + \psi_h\left(\frac{z_{0h}}{L}\right)\]

where \(\psi_h\) is the stability correction function for scalars (see compute_psih()).

compute_momentum_correction_term(z, oblen, z0m)

Compute momentum stability correction term.

Parameters:

• z (Array | float) – height above ground level \(z\).

• oblen (Array) – Obukhov length \(L\).

• z0m (Array) – roughness length for momentum \(z_{0m}\).

Returns:

The momentum stability correction.

Notes

This term is used in Monin-Obukhov similarity theory for momentum, and is given by

\[\ln\left(\frac{z}{z_{0m}}\right) - \psi_m\left(\frac{z}{L}\right) + \psi_m\left(\frac{z_{0m}}{L}\right)\]

where \(\psi_m\) is the stability correction function for momentum (see compute_psim()).

compute_psim(zeta)

Compute momentum stability function from Monin-Obukhov similarity theory.

Parameters:

zeta (Array) – stability parameter z/L \(\zeta\).

Returns:

Momentum stability correction [-].

Notes

This function calculates the integrated stability correction function for momentum \(\Psi_m\), which is used to adjust wind profiles based on atmospheric stability.

The function is piecewise, depending on the stability parameter \(\zeta = z/L\).

1. Unstable conditions (ζ ≤ 0):

Based on Businger-Dyer relations, an intermediate variable

\[x = (1 - 16\zeta)^{1/4}\]

is used to write the stability function as

\[\Psi_m(\zeta) = \ln\left( \frac{(1+x)^2 (1+x^2)}{8} \right) - 2 \arctan(x) + \frac{\pi}{2}.\]

2. Stable conditions (ζ > 0):

This uses an empirical formula (e.g., Holtslag and De Bruin, 1988) with constants:

• \(\alpha = 0.35\),

• \(\beta = 5.0 / \alpha\),

• \(\gamma = (10.0 / 3.0) / \alpha\).

The stability function is then given by

\[\Psi_m(\zeta) = -\frac{2}{3}(\zeta - \beta)e^{-\alpha \zeta} - \zeta - \gamma.\]

compute_psih(zeta)

Compute scalar stability function from Monin-Obukhov similarity theory.

Parameters:

zeta (Array) – stability parameter z/L \(\zeta\).

Returns:

The scalar stability correction.

Notes

This function calculates the integrated stability correction function for scalars \(\Psi_h\), like heat and humidity, which is used to adjust temperature and humidity profiles based on atmospheric stability.

The function is piecewise, depending on the stability parameter \(\zeta = z/L\).

1. Unstable conditions (ζ ≤ 0):

Based on Businger-Dyer relations, an intermediate variable (same as above)

\[x = (1 - 16\zeta)^{1/4}\]

is used to write the integrated stability function

\[\Psi_h(\zeta) = 2 \ln\left( \frac{1+x^2}{2} \right).\]

2. Stable conditions (ζ > 0):

This uses a corresponding empirical formula with the same constants (\(\alpha\), \(\beta\), \(\gamma\)) as above to write

\[\Psi_h(\zeta) = -\frac{2}{3}(\zeta - \beta)e^{-\alpha \zeta} - \left(1 + \frac{2}{3}\zeta\right)^{3/2} - \gamma + 1.\]