1. Introduction
Wall-bounded turbulent flows have been studied extensively for more than a century, and continue to be an active area of research (Kovasznay Reference Kovasznay1970; Smits, McKeon & Marusic Reference Smits, McKeon and Marusic2011; Jiménez Reference Jiménez2012, Reference Jiménez2018; Graham & Floryan Reference Graham and Floryan2021). A fundamental challenge in wall-bounded turbulent flow is to accurately predict mean flow properties such as velocity and stress distributions across the layer (Smits Reference Smits2022). For example, Prandtl (Reference Prandtl1925) derived that the mean velocity profile in the inertial sublayer follows the logarithmic law of the wall using his mixing length theory combined with the assumption of constant stress. Using dimensional analysis, von Kármán (Reference von Kármán1930) also obtained the logarithmic law of the wall. Coles (Reference Coles1956) showed using measurement data that the velocity profile in the turbulent boundary layer is described as the sum of the logarithmic law and the law of the wake. By now, the logarithmic law view is the generally accepted view, and its universality has been supported by experimental measurements (Marusic et al. Reference Marusic, Monty, Hultmark and Smits2013) and numerical simulations (Stevens, Wilczek & Meneveau Reference Stevens, Wilczek and Meneveau2014).
The atmospheric boundary layer (ABL) is the lowest part of the atmosphere that is directly influenced by the Earth’s surface (Monin Reference Monin1970; Panofsky Reference Panofsky1974; Wyngaard Reference Wyngaard1992; Mahrt Reference Mahrt2014). Its vertical structure plays a crucial role in the exchange of heat, mass and momentum between the surface and the free atmosphere. Therefore, it holds significant fundamental interest and is highly relevant for meteorological applications (Holtslag et al. Reference Holtslag2013; Stevens & Meneveau Reference Stevens and Meneveau2017; LeMone et al. Reference LeMone2019; Veers et al. Reference Veers2019). In the daytime, the ABL is usually convective due to the surface heating by solar radiation, whereas at night-time, it is usually stable due to radiative cooling of the ground surface (Katul, Konings & Porporato Reference Katul, Konings and Porporato2011). In convective ABLs, the flow is nearly uniform due to the intense vertical mixing, and its turbulent structures have been described very well by the mixed-layer scaling (Kaimal et al. Reference Kaimal, Wyngaard, Haugen, Coté, Izumi, Caughey and Readings1976; Tong & Ding Reference Tong and Ding2020; Liu, Gadde & Stevens Reference Liu, Gadde and Stevens2023). In contrast, the stable ABL remains one of the least understood of all types of ABLs because it may range from fully turbulent to intermittently turbulent or even non-turbulent at a variety of heights, temporal scales and spatial locations (Fernando & Weil Reference Fernando and Weil2010; Mahrt Reference Mahrt2014; Kit et al. Reference Kit, Hocut, Liberzon and Fernando2017).
The mean flow properties of the ABL are much more difficult to describe than other turbulent boundary layers due to the coexistence of the shear stress, buoyancy force and Coriolis force. In general, the ABL is vertically partitioned into the surface (Prandtl) layer and the Ekman layer. In the surface layer, which is approximately the lowest 10 % of the ABL, the thermal stratification is important but the Earth’s rotation is usually neglected. Note that this does not necessarily mean that the rotational effect is unimportant. For instance, Kostelecky & Ansorge (Reference Kostelecky and Ansorge2024) showed numerically that significant wind veer can be present in the surface layer. To account for the thermal stratification effect on the mean wind speed profile, Monin & Obukhov (Reference Monin and Obukhov1954) introduced a stability correction function
$\psi$
to the logarithmic law of the wall (Prandtl Reference Prandtl1925), and based on dimensional analysis proposed that
$\psi$
is only a function of the atmospheric stability parameter
$z/L_s$
. Here,
$z$
is the vertical height above the ground, and
$L_s$
is the Obukhov length based on the surface potential temperature flux (Obukhov Reference Obukhov1946). This is the well-known Monin–Obukhov similarity theory (MOST), which forms the basis for parametrisations of turbulent exchange in virtually all numerical models of atmospheric flows (Foken Reference Foken2006; Stiperski & Calaf Reference Stiperski and Calaf2023).
The stability correction function
$\psi$
in MOST cannot generally be inferred from dimensional considerations alone and is thus usually determined from empirical data (Businger et al. Reference Businger, Wyngaard, Izumi and Bradley1971; Dyer Reference Dyer1974). To predict the shape of
$\psi$
, Katul et al. (Reference Katul, Konings and Porporato2011) generalised the theoretical linkage between the mean velocity profile and the spectrum of turbulence by including the effects of thermal stratification on the turbulent kinetic energy dissipation rate and eddy-size anisotropy, where the length scale of the latter is selected as the distance from the wall. Their theory succeeded in predicting the mean velocity profile for unstable stratification, but failed in stable conditions. Li, Salesky & Banerjee (Reference Li, Salesky and Banerjee2016) hypothesised that this failure is due to the change in the characteristic size of turbulent eddies, where the Ozimidov scale becomes a more stringent limitation for the size of momentum transporting eddies than the distance from the wall under mildly to moderately stable conditions. More recently, Salesky & Anderson (Reference Salesky and Anderson2020) improved the performance of MOST under unstable conditions by introducing an additional dimensionless parameter that accounts for the effects of large coherent structures, and Stiperski & Calaf (Reference Stiperski and Calaf2023) generalised MOST to complex atmospheric turbulence, based on the inclusion of turbulence anisotropy as an additional non-dimensional term. Additionally, Heisel & Chamecki (Reference Heisel and Chamecki2023) demonstrated that a mixed scaling parameter incorporating both the Obukhov length and the stable ABL depth improves similarity for the mean profiles of wind speed and air temperature compared to MOST.
In the Ekman layer lying above the surface layer, the influence of the Coriolis force becomes apparent. A fundamental feature of the flow therein is that the direction of the velocity varies with height, which is known as the Ekman spiral. In addition, the flow in the Ekman layer usually displays a low-level jet with a velocity higher than the free-stream velocity (Deusebio et al. Reference Deusebio, Brethouwer, Schlatter and Lindborg2014), which is also different from other wall-bounded turbulent flows. As there are relatively few reliable observations, the Ekman layer is more difficult to understand (Panofsky Reference Panofsky1974). For truly neutral ABLs, Ekman (Reference Ekman1905) first provided an analytical solution for the velocity profiles by assuming a constant eddy viscosity. Parmhed, Kos & Grisogono (Reference Parmhed, Kos and Grisogono2005) extended the Ekman model analytically for any smoothly varying eddy viscosity using the Wentzel–Kramers–Brillouin method. Constantin & Johnson (Reference Constantin and Johnson2019) revisited the Ekman-type solutions for arbitrary height-dependent eddy viscosities, and showed that the mean geostrophic wind deficit profile always decays in magnitude and turns clockwise with increasing height in the Northern Hemisphere.
For conventionally neutral ABLs, the wind speed profile is usually described by the logarithmic law of the wall with an additional linear or quadratic term with height to account for the effect of free-atmosphere stratification (Zilitinkevich & Esau Reference Zilitinkevich and Esau2005; Gryning et al. Reference Gryning, Batchvarova, Brümmer, Jørgensen and Larsen2007; Kelly, Cersosimo & Berg Reference Kelly, Cersosimo and Berg2019). Based on the similarity theory, Liu, Gadde & Stevens (Reference Liu, Gadde and Stevens2021) derived an analytical expression for the mean wind speed profile that accounts for the effect of local thermal stratification. By solving the Ekman equations, Liu & Stevens (Reference Liu and Stevens2022) further derived the analytical expressions for the mean turbulent momentum fluxes and geostrophic wind deficits, which agree well with simulation results and atmospheric measurements. For convective ABLs, Tong & Ding (Reference Tong and Ding2020) derived the convective logarithmic friction law from first principles, which connects the mixed-layer mean velocity scale and the friction velocity in the convective-roll-dominated regime. They also analytically derived the mean velocity profile using the method of matched asymptotic expansions. However, this model is valid only in approximately the lower half of the boundary layer, as it does not account for the effects of the entrainment zone. Recently, Liu et al. (Reference Liu, Gadde and Stevens2023) developed an analytical model of the mean velocity profile that connects the surface layer, the mixed layer and the entrainment zone, using a perturbation method approach in conjuncture with the convective logarithmic friction law (Tong & Ding Reference Tong and Ding2020) and MOST (Monin & Obukhov Reference Monin and Obukhov1954). For stable ABLs, the characteristics of the Ekman layer are not as well understood now as they are in other types of ABL. Therefore, they require further study.
Stable ABLs prevail at nocturnal conditions under relatively clear skies, but may also survive for several hours after sunrise over continental lands, and persist continuously for days in polar regions (Mahrt Reference Mahrt2014; Steeneveld Reference Steeneveld2014; Jiang & Wang Reference Jiang and Wang2021). They can be influenced by multiple processes such as turbulence, radiation, gravity waves, katabatic winds, fog and dew formation. The stable ABL is observed to exist in at least two regimes, commonly referred to as the very stable ABL (vSBL) and the weakly stable ABL (wSBL) (Mahrt et al. Reference Mahrt, Sun, Blumen, Delany and Oncley1998; Wyngaard Reference Wyngaard2010; Mahrt Reference Mahrt2014). The wSBL is characterised by weak stratification, strong winds and sustained turbulence, with turbulence decreasing with height, and near-surface profiles described well by MOST. In contrast, the vSBL is marked by strong stratification, weak winds and intermittent turbulence, where turbulence can decouple from the surface and MOST breaks down. Transitions between these two regimes are frequently observed (Vercauteren, Mahrt & Klein Reference Vercauteren, Mahrt and Klein2016; van de Wiel et al. Reference van de Wiel, Vignon, Baas, van Hooijdonk, van der Linden, Antoon van Hooft, Bosveld, de Roode, Moene and Genthon2017; Abraham, Holdsworth & Monahan Reference Abraham, Holdsworth and Monahan2019; van der Linden et al. Reference van der Linden, Edwards, van Heerwaarden, Vignon, Genthon, Petenko, Baas, Jonker and van de Wiel2019; Abraham & Monahan Reference Abraham and Monahan2020). In this study, we use relatively small constant cooling rates and relatively strong geostrophic winds. Therefore, our results are expected to be valid only for wSBL, not for vSBL.
The transition from wSBL to vSBL is typically driven by radiative cooling at the surface over land, while over water it is caused by the advection of warm air aloft. These processes are relatively well understood (van de Wiel et al. Reference van de Wiel, Vignon, Baas, van Hooijdonk, van der Linden, Antoon van Hooft, Bosveld, de Roode, Moene and Genthon2017; Abraham et al. Reference Abraham, Holdsworth and Monahan2019; Kaiser, Vercauteren & Krumscheid Reference Kaiser, Vercauteren and Krumscheid2024). However, the vSBL to wSBL transition is less well understood as various mechanisms are associated with such events, including shear instabilities, cloud cover and intermittent turbulence (Abraham et al. Reference Abraham, Holdsworth and Monahan2019; Abraham & Monahan Reference Abraham and Monahan2020). To understand regime transitions, van de Wiel et al. (Reference van de Wiel, Vignon, Baas, van Hooijdonk, van der Linden, Antoon van Hooft, Bosveld, de Roode, Moene and Genthon2017) introduced a conceptual atmospheric-surface model that combines a surface energy budget with a bulk parametrisation for turbulent transport. Based on theoretical analysis, the authors demonstrated that the stable ABL exhibits bistability within a certain range of parameters; this means that under identical wind forcing conditions, the ABL can exist in either the wSBL or vSBL state, depending on factors such as initial conditions, turbulence dynamics and surface properties. Bistability has significant implications for the predictability of the stable ABL as small perturbations or uncertainties in initial and forcing conditions can result in vastly different mean states of the stable ABL (Kaiser et al. Reference Kaiser, Vercauteren and Krumscheid2024). Nevertheless, van der Linden et al. (Reference van der Linden, Edwards, van Heerwaarden, Vignon, Genthon, Petenko, Baas, Jonker and van de Wiel2019) observed at the Dome C site that although the trajectories of both transitions do not overlap for intermediate wind speeds, they do appear to overlap for both low and high wind speeds. This suggests the presence of well-defined robust steady states that determine whether the system settles into a vSBL or wSBL state, which is the focus of our analysis. Finally, we highlight that analysing these steady states is beneficial for enhancing theoretical predictions and gaining deeper insights into boundary layer dynamics.
The stable ABL can also be classified as nocturnal or long-lived, depending on whether the flow immediately aloft is neutrally or stably stratified (Zilitinkevich Reference Zilitinkevich2002; Zilitinkevich, Esau & Baklanov Reference Zilitinkevich, Esau and Baklanov2007). The long-lived stable ABL situates immediately below the stably stratified free atmosphere, and hence, the surface-layer turbulence is affected by the free-atmosphere stability (Zilitinkevich & Calanca Reference Zilitinkevich and Calanca2000). In contrast, the nocturnal stable atmospheric boundary layer (NSBL) is decoupled from the free atmosphere by a thick neutral or weakly stable residual layer (Kim & Mahrt Reference Kim and Mahrt1992; Shen et al. Reference Shen, Liu, Lu and Stevens2024). The residual layer prevents internal-wave interactions between the stable boundary layer and the free atmosphere. Within the surface layer, the flow is well described by MOST (Obukhov Reference Obukhov1946; Monin & Obukhov Reference Monin and Obukhov1954; Foken Reference Foken2006). However, within the Ekman layer, which is situated above the surface layer, the Coriolis effect becomes significant, thus MOST no longer applies (Optis et al. Reference Optis, Monahan and Bosveld2014, Reference Optis, Monahan and Bosveld2016; Liu et al. Reference Liu, Gadde and Stevens2021). We emphasise that although the NSBLs prevail at nocturnal conditions, which would typically only occur and develop for less than 12 h (i.e. from sunset to sunrise), they can also survive for long periods in middle latitudes in the presence of a large flat glacier (Vercauteren et al. Reference Vercauteren, Mahrt and Klein2016), the entire day in valleys and basins with low winter sun angles (Mahrt Reference Mahrt2014), and for many days in polar regions in winter in the absence of the diurnal cycle (Bon, Cal & Meyers Reference Bon, Cal and Meyers2024).
Hitherto, many studies have addressed the prediction of mean wind profiles within the Ekman layer of stable ABLs. These profiles are typically obtained by solving the Ekman equations using an assumed eddy viscosity (Shapiro & Fedorovich Reference Shapiro and Fedorovich2010; Van de Wiel et al. Reference van de Wiel, Moene, Steeneveld, Baas, Bosveld and Holtslag2010; Kadantsev, Mortikov & Zilitinkevich Reference Kadantsev, Mortikov and Zilitinkevich2021) or turbulent shear stress (Gryning et al. Reference Gryning, Batchvarova, Brümmer, Jørgensen and Larsen2007; Basu & Holtslag Reference Basu and Holtslag2023; Narasimhan, Gayme & Meneveau Reference Narasimhan, Gayme and Meneveau2024). In particular, Nieuwstadt (Reference Nieuwstadt1984b
) introduced the local scaling hypothesis, suggesting that the dimensionless turbulence characteristics can be expressed as functions of the single parameter
$z/L$
, where
$L$
is the local Obukhov length. This hypothesis could be regarded as an extension of MOST to the entire stable ABL, which describes the measurement (Nieuwstadt Reference Nieuwstadt1984b
; Sorbjan Reference Sorbjan1986; Basu et al. Reference Basu, Porté-Agel, Foufoula-Georgiou, Vinuesa and Pahlow2006) and simulations (Lacser & Arya Reference Lacser and Arya1986; Derbyshire Reference Derbyshire1990; Coleman, Ferziger & Spalart Reference Coleman, Ferziger and Spalart1992; van Dop & Axelsen Reference van Dop and Axelsen2007) well. Furthermore, Nieuwstadt (Reference Nieuwstadt1984b
) derived an analytical model of stable ABLs by assuming the flux Richardson number
$R_{\!f}$
, the ratio of the buoyancy flux and turbulent production by mean shear, to be constant throughout the boundary layer. This leads to a power law of the total turbulent shear stress
$\tau$
with exponent
$\alpha =3/2$
. However, later studies (Coleman et al. Reference Coleman, Ferziger and Spalart1992; Kosović & Curry Reference Kosović and Curry2000; Saiki, Moeng & Sullivan Reference Saiki, Moeng and Sullivan2000; Basu & Porté-Agel Reference Basu, Porté-Agel, Foufoula-Georgiou, Vinuesa and Pahlow2006; Huang & Bou-Zeid Reference Huang and Bou-Zeid2013) showed that
$R_{\!f}$
is not constant but continuously increases all the way up to the top of the boundary layer, implying that
$\tau$
may damp more rapidly than previously assumed. Hence different values of
$\alpha$
have been reported in the literature (Caughey, Wyngaard & Kaimal Reference Caughey, Wyngaard and Kaimal1979; Yokoyama, Gamo & Yamamoto Reference Yokoyama, Gamo and Yamamoto1979; Nieuwstadt Reference Nieuwstadt1984b
; Lacser & Arya Reference Lacser and Arya1986; Sorbjan Reference Sorbjan1986, Reference Sorbjan1988; Lenschow et al. Reference Lenschow, Li, Zhu and Stankov1988; Wittich Reference Wittich1991; Grant Reference Grant1997). For instance, Gryning et al. (Reference Gryning, Batchvarova, Brümmer, Jørgensen and Larsen2007) assumed
$\alpha =2$
in their formulation of the wind profile beyond the surface layer. In addition, we note that the velocity deficit profiles in steady state can be obtained by taking the vertical derivatives of shear stresses (Garratt Reference Garratt1992). However, these profiles have seldom been discussed due to a lack of high-fidelity data. This shows that the vertical structure of stable ABLs is still not well understood, which limits the effective estimation of wind characteristics at high altitudes for wind energy applications (Archer & Caldeira Reference Archer and Caldeira2009; Bechtle et al. Reference Bechtle, Schelbergen, Schmehl, Zillmann and Watson2019). Therefore, predicting analytically and understanding deeply the mean turbulent momentum fluxes and geostrophic wind deficits in the Ekman layers of the NSBLs are the main concern of this study.
Improvements in computational power and techniques now enable highly accurate simulations. As the Reynolds numbers in the atmospheric flows are very high, large-eddy simulations (LES) are commonly employed to model these flows (Kim & Leonard Reference Kim and Leonard2024). However, accurately simulating stable ABLs is still not an easy task. For example, the grid convergence of stable ABL simulations has been investigated and was found to depend on various factors, such as the wall modelling (Maronga, Knigge & Raasch Reference Maronga, Knigge and Raasch2020), the roughness length (Couvreux et al. Reference Couvreux2020), the sub-grid scale model (Dai et al. Reference Dai, Basu, Maronga and de Roode2021), and the simulation duration (Maronga & Li Reference Maronga and Li2022). Although this issue remains partially unresolved, it is widely accepted that the lowest grid level should exceed the lower limit of the inertial sublayer (Garratt Reference Garratt1983; Basu & Lacser Reference Basu and Lacser2017) and fall below the Ozmidov length (Dougherty Reference Dougherty1961; Ozmidov Reference Ozmidov1965; Huang & Bou-Zeid Reference Huang and Bou-Zeid2013; Sullivan et al. Reference Sullivan, Weil, Patton, Jonker and Mironov2016).
Taking these lessons into account, Shen et al. (Reference Shen, Liu, Lu and Stevens2024) performed a unique set of LES of quasi-steady, barotropic and horizontally homogeneous NSBLs covering a wider range of the dimensionless parameters
$\mu =u_*/(f L_s) \in [16.7, 193.3]$
and
$Ro = u_*/(fz_0) \in [1.68\times 10^4\!, 2.42 \times 10^6]$
than previously considered. Here,
$\mu$
is the Kazanski–Monin parameter (Kazanski & Monin Reference Kazanski and Monin1961),
$\textit{Ro}$
is the friction Rossby number (Rossby & Montgomery Reference Rossby and Montgomery1935; Tennekes & Lumley Reference Tennekes and Lumley1972),
$u_*$
is the friction velocity,
$f$
is the Coriolis parameter, and
$z_0$
is the roughness height. Note that in the quasi-steady situation, the mean profiles of velocity and turbulence statistics remain statistically steady, whereas only the potential temperature profile keeps changing (Brost & Wyngaard Reference Brost and Wyngaard1978; Bon et al. Reference Bon, Cal and Meyers2024). This is the case since a constant cooling rate is set as the bottom boundary condition of the potential temperature. Utilising this set of simulations, Shen et al. (Reference Shen, Liu, Lu and Stevens2024) investigated the global properties such as the boundary layer height and the geostrophic drag coefficient of NSBLs, and validated their newly derived analytical expressions for geostrophic drag law coefficients. In this study, we go beyond by focusing on the mean flow properties in the Ekman layers of the quasi-steady, barotropic and horizontally homogeneous NSBLs. In particular, we determine analytically the mean turbulent momentum fluxes and wind deficits profiles in NSBLs. This high-fidelity database thus enables us to evaluate the performance of these analytical profiles and deepen our understanding of the mean flow properties of NSBLs.
The organisation of the paper is as follows. In § 2, we present an analytical model for the mean turbulent momentum fluxes and geostrophic wind deficits in NSBLs. In § 3, the simulation database of Shen et al. (Reference Shen, Liu, Lu and Stevens2024) is introduced briefly. In § 4, the model predictions are compared with the high-fidelity simulations and field measurements. We conclude with a summary of the main findings in § 5.
2. Analytical model
This study focuses on the quasi-steady, barotropic and horizontally homogeneous NSBLs over flat surfaces. For simplicity, we focus on the Northern Hemisphere with latitude
$\phi \gt 0$
. From dimensional analysis, these flows are governed by two dimensionless parameters, i.e. the friction Rossby number
$Ro=u_*/(fz_0)$
(Rossby & Montgomery Reference Rossby and Montgomery1935) and the Kazanski–Monin parameter
$\mu = u_*/(f L_s)$
(Kazanski & Monin Reference Kazanski and Monin1961). Here,
$u_*$
is the friction velocity,
$f=2\Omega \sin \phi$
is the Coriolis parameter,
$\Omega$
is the Earth’s rotation speed,
$z_0$
is the roughness length,
$L_s = - u_*^3 /(\beta q_s)$
is the surface Obukhov length (Obukhov Reference Obukhov1946),
$\beta$
is the buoyancy parameter, and
$q_s$
is the potential temperature flux at the surface. For
$\mu =O(10^3)$
, the flow is intermittent and may no longer fully satisfy the usual conditions for the definition of turbulence (Mahrt Reference Mahrt2014); while for
$|\mu |\ll 10$
, the flow is traditionally considered neutral (Zilitinkevich & Esau Reference Zilitinkevich and Esau2002). As will be shown below, for nearly neutral conditions
$q_s$
approaches to zero therein, thus it can no longer be used as a scaling to derive our analytical model. Therefore, in this study we consider
$\mu = 30{-}200$
.
Under these conditions, the NSBLs are governed by the Ekman equations, which describe the force balance between the pressure gradient force, Coriolis force and turbulent drag (Ekman Reference Ekman1905; Narasimhan et al. Reference Narasimhan, Gayme and Meneveau2024):
\begin{equation} \frac {\textrm {d} \tau _x}{\textrm {d} z} = -f(V-V_g), \quad \frac {\textrm {d} \tau _y}{\textrm {d} z} = f(U-U_g), \end{equation}
where
$\tau _x, \tau _y$
are the streamwise and spanwise turbulent momentum fluxes,
$U, V$
are the streamwise and spanwise velocities, and
$U_g, V_g$
are the corresponding geostrophic wind components. The streamwise axis is parallel to the wind direction at the surface, the vertical axis is pointing upwards, and the streamwise (
$x$
), spanwise (
$y$
) and vertical (
$z$
) axes form a right-handed coordinate system (see figure 1). In the literature, the boundary layer height is commonly defined as the height where the total turbulent momentum flux
$\tau = (\tau _x^2+\tau _y^2)^{1/2}$
attains 5 % of its surface value, which we denote as
$h_{0.05}$
. Here, we define the boundary layer height as the height where
$\tau$
first reduces to zero, and denote it as
$h$
. The relation between
$h$
and
$h_{0.05}$
will be explained later. Under this definition and the coordinate system sketched in figure 1, the boundary conditions of (2.1) are
\begin{equation} \tau _x(0)=u_*^2, \quad \tau _x(h)=\tau _y(0)=\tau _y(h)=0. \end{equation}
From (2.1) and (2.2), the turbulent momentum fluxes are usually approximated as (Nieuwstadt Reference Nieuwstadt1984b ; Derbyshire Reference Derbyshire1990; Wyngaard Reference Wyngaard2010)
\begin{equation} \frac {\tau _x + \textrm{i} \tau _y}{u_*^2} = \left ( 1- \frac {z}{h} \right )^{\alpha + \textrm{i} \gamma }, \end{equation}
where
$\textrm{i}$
is the imaginary unit, and
$(\alpha , \gamma )$
are empirical constants. Equation (2.3) indicates that the total turbulent momentum flux
$\tau$
satisfies the power law
$\tau /u_*^2 = (1-z/h)^\alpha$
. Thus the height
$h$
can be determined from the height
$h_{0.05}$
:
\begin{equation} h = \frac { h_{0.05} }{ 1-0.05^{1/\alpha } }. \end{equation}
Substituting (2.3) into (2.1), we can obtain the geostrophic wind deficits:
\begin{equation} \frac {h f}{u_*^2} [ (U-U_g) + \textrm{i} (V-V_g)] = \textrm{i} (\alpha + \textrm{i} \gamma ) \left (1- \frac {z}{h} \right )^{\alpha + \textrm{i} \gamma -1 }. \end{equation}
Equation (2.5) indicates that the wind deficits scale as
$u_*^2/(hf)$
, which is consistent with our recent finding that the wind gradients scale as
$u_*^2/(h^2 f)$
(Shen et al. Reference Shen, Liu, Lu and Stevens2024). Note that (2.3) and (2.5) are not closed since the values of
$\alpha$
and
$\gamma$
are unknown yet.
Schematic representation of NSBLs (not to scale). It shows the velocity and potential temperature profiles as functions of height, and the typical Ekman spiral in the boundary layer. Here,
$f$
is the Coriolis parameter,
$h$
is the boundary layer height, and
$(U,V)$
are the velocity components along the
$(x,y)$
axes, respectively. The
$x$
axis is parallel to the wind on the surface, the
$z$
axis points upwards, and
$(x,y,z)$
form a right-handed coordinate system. Adapted from Liu & Stevens (Reference Liu and Stevens2022), who adapted and modified the figure from Allaerts & Meyers (Reference Allaerts and Meyers2015).
To determine the value of
$\gamma$
, we further assume the eddy viscosity approach to be valid in the NSBLs such that
\begin{equation} \tau _x = -K_m \frac { \textrm {d} U}{ \textrm {d} z}, \quad \tau _y = -K_m \frac { \textrm {d} V}{ \textrm {d} z}, \end{equation}
where
$K_m$
is the eddy viscosity, which is a function of height. Thus from (2.3), (2.5) and (2.6), the eddy viscosity
$K_m$
can be determined as
\begin{equation} K_m = \frac {h^2 f}{ (2\alpha -1) \gamma - \textrm{i}\left[\alpha (\alpha -1) - \gamma ^2\right] } \left (1- \frac {z}{h} \right )^{2}. \end{equation}
As
$K_m$
should be a real number, one obtains
\begin{equation} \gamma = \sqrt {\alpha (\alpha -1)}. \end{equation}
We note that (2.7) is valid only in the Ekman layer, and does not represent surface layer conditions, which has also been reported by Basu & Holtslag (Reference Basu and Holtslag2023). In the surface layer, the eddy viscosity
$K_m$
follows MOST (Monin & Obukhov Reference Monin and Obukhov1954; Brost & Wyngaard Reference Brost and Wyngaard1978):
\begin{equation} K_m = \frac {\kappa u_* z}{1 + c_\psi \kappa z/L_s}, \end{equation}
where
$\kappa =0.4$
is the von Kármán constant, and
$c_\psi$
is an empirical constant to be determined later (see (2.16) below).
To determine the value of
$\alpha$
, we introduce the flux Richardson number
\begin{equation} R_f \equiv \frac {- \beta q }{ \displaystyle \tau _x \frac {\textrm {d} U}{\textrm {d} z} + \tau _y \frac {\textrm {d} V}{\textrm {d} z} } = \frac {c_s^2} { (2\alpha -1) \gamma } \left ( 1- \frac {z}{h} \right )^{ 2-2\alpha } \frac {q}{q_s}, \end{equation}
which represents the ratio of the buoyancy flux and turbulent production by mean shear in the turbulent kinetic energy equation (Williams et al. Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017). Here,
\begin{equation} c_s = \sqrt {(h^2 f) / (L_s u_*)} \end{equation}
has been shown to be a constant for a large value of
$\mu$
, and is called the Zilitinkevich constant (Derbyshire Reference Derbyshire1990; Shen et al. Reference Shen, Liu, Lu and Stevens2024), and
$q$
is the potential temperature flux, which follows a linear profile (Nieuwstadt Reference Nieuwstadt1984b
),
\begin{equation} \frac {q}{q_s} = 1- \frac {z}{h}. \end{equation}
Nieuwstadt (Reference Nieuwstadt1984b
) assumed
$R_{\!f}$
to be constant throughout the stable boundary layer, thus from (2.10) and (2.12), we have
$\alpha =3/2$
. However, later studies showed that
$R_{\!f}$
is not constant but continuously increases all the way up to the top of the boundary layer (Coleman et al. Reference Coleman, Ferziger and Spalart1992; Kosović & Curry Reference Kosović and Curry2000; Saiki et al. Reference Saiki, Moeng and Sullivan2000; Basu & Porté-Agel Reference Basu, Porté-Agel, Foufoula-Georgiou, Vinuesa and Pahlow2006; Huang & Bou-Zeid Reference Huang and Bou-Zeid2013). This is because the potential temperature flux
$q$
vanishes at a slightly higher elevation than the turbulent momentum fluxes
$(\tau _x, \tau _y)$
; see e.g. Basu & Porté-Agel (Reference Basu, Porté-Agel, Foufoula-Georgiou, Vinuesa and Pahlow2006), Shen et al. (Reference Shen, Liu, Lu and Stevens2024) and figures 2–5 below. Therefore, to provide a more accurate parametrisation of
$R_{\!f}$
, especially in the upper part of the boundary layer, we propose using a reciprocal curve with the origin located at the boundary layer top,
$(1-z/h)^{-1}$
, to approximate the flux Richardson number
$R_{\!f}$
. This ensures that
$R_{\!f}$
approaches infinity at the boundary layer height
$z/h=1$
. Under this assumption,
$R_{\!f}$
can be parametrised as
\begin{equation} R_f = \frac {c_s^2}{(2\alpha -1) \gamma } \left ( 1- \frac {z}{h} \right )^{ -1 }. \end{equation}
Finally, by combining (2.10)–(2.13), we find
\begin{equation} \alpha =2. \end{equation}
We remark that
$R_f=0.2$
proposed by Nieuwstadt (Reference Nieuwstadt1984b
) can be regarded as a constant approximation, and (2.13) proposed by the present study can be regarded as a linear approximation of the inverse flux Richardson number at moderate thermal stratification with
$\mu = 30{-}200$
. In addition, (2.14) rationalises the choice of Gryning et al. (Reference Gryning, Batchvarova, Brümmer, Jørgensen and Larsen2007), who took
$\alpha =2$
for the purpose of simplicity. We further remark that for nearly neutral conditions with
$\mu \ll 10$
, (2.13) can no longer be used to close (2.3) and (2.5) as
$q_s$
approaches zero. In addition, within the inner layer,
$R_{\!f}$
scales as (Monin & Obukhov Reference Monin and Obukhov1954; Wyngaard Reference Wyngaard2010, p. 282)
\begin{equation} R_f = \frac {\kappa z/L_s}{1 + c_\psi \kappa z/L_s}. \end{equation}
By asymptotic matching (2.13) as
$z/h \rightarrow 0$
with (2.15) as
$z/L_s \rightarrow \infty$
, we obtain
\begin{equation} c_\psi = \frac {(2\alpha -1) \gamma }{c_s^2}. \end{equation}
Finally, we remark that the wind speed profile in high altitudes can be predicted analytically from the combination of (2.5), the boundary layer parametrisation (Zilitinkevich Reference Zilitinkevich1972; Basu & Holtslag Reference Basu and Holtslag2023)
\begin{equation} \frac {hf}{u_*} = \frac {c_s}{ \sqrt {\mu } }, \end{equation}
and the geostrophic drag law for NSBLs that relates the friction velocity and the geostrophic wind (Zilitinkevich & Esau Reference Zilitinkevich and Esau2005; Shen et al. Reference Shen, Liu, Lu and Stevens2024),
\begin{align}& A = \ln Ro - \frac {\kappa U_g}{u_*} = \ln \mu - a_1 \sqrt {\mu } - a_2, \end{align}
\begin{align}&\qquad B = - \frac {\kappa V_g}{u_*} = b_1 \sqrt {\mu }+ \frac { b_2 }{ \sqrt {\mu } }. \end{align}
Here,
$A, B$
are the geostrophic drag law coefficients for NSBLs, and
$a_1, a_2, b_1, b_2$
are the empirical constants. Shen et al. (Reference Shen, Liu, Lu and Stevens2024) calculated the boundary layer height
$h$
using (2.4) with the exponent
$\alpha =3/2$
, and determined the Zilitinkevich constant
$c_s=0.62$
using a least squares fitting procedure. In this study, the boundary layer height
$h$
is calculated using (2.4) with the exponent
$\alpha =2$
. Since the same dataset is used for both studies, the model-independent height
$h_{0.05}$
should be unchanged. Thus in this study, the Zilitinkevich constant is
$c_s= (1-0.05^{3/2})/(1-0.05^{1/2}) \times 0.62= 0.69$
. We note that the value of
$c_s$
has no effect on the geostrophic drag law parametrisation, therefore the values of the four empirical constants
$(a_1, a_2, b_1, b_2)=(0.72, 0.04, 1.02, 4.41)$
are exactly the same as those in Shen et al. (Reference Shen, Liu, Lu and Stevens2024). Finally, we note that the wind speed profile within the surface layer follows MOST (Monin & Obukhov Reference Monin and Obukhov1954):
\begin{equation} \frac {\sqrt {U^2 +V^2}}{u_*} = \frac {1}{\kappa } \ln \left ( \frac {z}{z_0} \right ) + c_\psi \frac {z}{L_s}. \end{equation}
From (2.8), (2.14) and (2.16), the empirical constant is
$c_\psi =8.9$
, which happens to be very close to the value (
$c_\psi =9$
) obtained by Shah & Bou-Zeid (Reference Shah and Bou-Zeid2014) using direct numerical simulations at lower Reynolds numbers.
3. Simulation database
The simulation database used for this study has been described extensively in Shen et al. (Reference Shen, Liu, Lu and Stevens2024), therefore it is only briefly outlined here. The database was generated using LES, which have the benefit that all quantities can be measured under perfectly known conditions, allowing different effects to be disentangled. The simulations were performed with an in-house code (Liu & Stevens Reference Liu and Stevens2021), which employs a pseudo-spectral discretisation in the horizontal directions, and a second-order finite difference method in the vertical direction. At the top boundary, the vertical velocity and the vertical potential temperature gradient were enforced to zero. At the bottom boundary, a wall model based on MOST for the velocity field and a constant cooling rate for the potential temperature field were employed. The advanced anisotropic minimum dissipation model was used to parametrise the sub-grid scale momentum flux and potential temperature flux (Abkar & Moin Reference Abkar and Moin2017). The computational domain size was
$800\;\rm m \times\, 800\;\rm m \times \,800\;\rm m $
in the streamwise, spanwise and vertical directions, with grid resolution
$192\times 192\times 201$
. It has been shown that the mean wind speed and total turbulent momentum flux profiles change only slightly with grid resolution. Additionally, the variation with grid resolution becomes less important when the results are rescaled and height is expressed relative to the boundary layer depth. Consequently, the chosen grid resolution has minimal impact on the comparison between model and simulation results.
The flows were initialised with a constant geostrophic wind speed and a constant potential temperature, with random noise added to these profiles below 50 m to trigger the turbulence. To ensure that the quasi-steady states were reached, and to minimise the effects of inertial oscillations on the mean wind profiles (Momen & Bou-Zeid Reference Momen and Bou-Zeid2017; Maronga & Li Reference Maronga and Li2022), the simulations were performed for 2.5 inertial periods
$(5\unicode{x03C0} /f)$
, and the statistics were collected over the final inertial period. Because the statistics are collected over one inertial period, NSBLs are implicitly assumed to persist for at least 24 h at middle latitudes, and 12 h at high latitudes. This suggests that the simulated mean profiles are more likely to occur in the real world at high latitudes, where the night usually lasts more than 12 h in winter. Nevertheless, the simulation database covers the values of
$\textit{Ro}$
and
$\mu$
found in typical NSBLs at mid to high latitudes (Caughey et al. Reference Caughey, Wyngaard and Kaimal1979; Nieuwstadt Reference Nieuwstadt1984a
; Lenschow et al. Reference Lenschow, Li, Zhu and Stankov1988). For our purposes, we consider only the cases with
$\mu =O(10^2)$
, i.e.
$\mu \in [34.0, 193.3]$
and
$Ro \in [1.68\times 10^4, 2.42 \times 10^6]$
. A summary of simulations considered in this study is presented in table 1, which is taken from Shen et al. (Reference Shen, Liu, Lu and Stevens2024) and rearranged such that the value of
$\mu$
increases monotonically.
Summary of LES performed by Shen et al. (Reference Shen, Liu, Lu and Stevens2024). Here,
$\phi$
is the latitude,
$C_r$
is the cooling rate,
$z_0$
is the roughness height,
$G$
is the geostrophic wind speed,
$h_{0.05}$
is the boundary-layer height where the total turbulent momentum flux attains 5 % of its surface value,
$L_s = - u_*^3 /(\beta q_s)$
is the surface Obukhov length,
$\mu$
is the Kazanski–Monin parameter, and
$\textit{Ro}$
is the friction Rossby number.
The mean vertical profiles of the normalised potential temperature flux
$q/q_s$
. Filled circles indicate simulation data from table 1; open diamonds indicate field data of Nieuwstadt (Reference Nieuwstadt1984b
); dashed line indicates theoretical prediction by Nieuwstadt (Reference Nieuwstadt1984b
); solid line indicates theoretical prediction by present model of (2.12) with
$\alpha =2$
.
4. Comparisons to field and simulation data
4.1. Mean turbulent momentum fluxes
Figure 2 shows the temporally and horizontally averaged vertical profile of the normalised potential temperature flux
$q/q_s$
, where the height
$z$
is normalised by the boundary layer height
$h$
. Note that the height
$h$
is calculated using (2.4) with
$\alpha =2$
. For fair comparisons, the height
$h$
of Nieuwstadt’s model with
$\alpha =3/2$
has been adjusted by the factor
$(1-0.05^{1/2})/(1-0.05^{2/3})$
to ensure that the model-independent height
$h_{0.05}$
remains unchanged. The circles are the simulation data, the diamonds are the measurement data of Nieuwstadt (Reference Nieuwstadt1984b
), the dashed line is the theoretical prediction by Nieuwstadt (Reference Nieuwstadt1984b
), and the solid line is the theoretical prediction of (2.12) with
$\alpha =2$
. Clearly, the theoretical prediction of Nieuwstadt’s model with
$\alpha =3/2$
(dashed line) slightly underestimates the simulations, while the theoretical prediction of (2.12) with
$\alpha =2$
(solid line) agrees well with all simulations and field data. Note that the profile of (2.12) only depends on
$\alpha$
through the determination of
$h$
via (2.4). The error near the top of the boundary layer therefore is mainly due to the
$q$
profiles not following the linear relationship in (2.12) at that point. This implies that there will always be a mismatch for any parametrisation of
$h$
. In other words, the mismatch near the boundary layer top cannot be fully eliminated by optimising the value of
$\alpha$
. While
$\alpha =2$
appears to be more accurate, as indicated by figure 2, further analysis of other flow quantities, such as the turbulent momentum fluxes, is needed to confirm this finding.
Figure 3(a) shows the mean vertical profile of the flux Richardson number
$R_{\!f}$
in the Ekman layer. The figure reveals that the assumption
$R_f=0.2$
(dashed line) made by Nieuwstadt (Reference Nieuwstadt1984b
) is only approximately valid in the lower part of the boundary layer. Near the top of the boundary layer,
$R_{\!f}$
becomes very large as the potential temperature flux vanishes at a slightly higher elevation than the turbulent momentum fluxes (see figure 2, and figure 5 below). In contrast, the prediction of (2.13) (solid line) agrees well with the simulation data (circles) and is nearly the same as assumed by Nieuwstadt (Reference Nieuwstadt1984b
) within the surface layer (dashed line). This result confirms the correctness of (2.13) and hence (2.14). Note that in (2.13),
$c_s=0.69$
, while in Shen et al. (Reference Shen, Liu, Lu and Stevens2024)
$c_s=0.62$
. This difference is due to the fact that Shen et al. (Reference Shen, Liu, Lu and Stevens2024) calculated the boundary layer height
$h$
using (2.4) with the exponent
$\alpha =3/2$
, while in this study the height
$h$
is calculated using (2.4) with the exponent
$\alpha =2$
; see also the remark at the end of § 2. Figure 3(b) shows the mean vertical profile of the flux Richardson number
$R_{\!f}$
in the surface layer. According to Shen et al. (Reference Shen, Liu, Lu and Stevens2024), the height of the surface layer is approximately
$L_s/2.25$
, therefore the range of
$z/L_s$
shown in the figure spans from 0 to 0.5. Clearly, all data (filled circles) collapse to a single curve of the dimensionless parameter
$z/L_s$
described by (2.15) (solid line). This confirms that the flux Richardson number
$R_{\!f}$
profile in the surface layer follows MOST.
The mean vertical profiles of the flux Richardson number
$R_{\!f}$
in (a) the Ekman layer and (b) the surface layer. Filled circles indicate simulation data from table 1; open diamonds indicate field data of Nieuwstadt (Reference Nieuwstadt1984b
); dashed line indicates theoretical prediction by Nieuwstadt (Reference Nieuwstadt1984b
); solid line indicates theoretical prediction by (a) present model of (2.13) and (b) MOST of (2.15), with
$\alpha =2$
,
$c_s=0.69$
.
Figure 4 shows the mean vertical profile of the total turbulent momentum flux
$\tau$
normalised by the scaling
$u_*^2$
. All simulation data in table 1 are shown in the figure, which shows that the profiles of the dimensionless total turbulent momentum flux
$\tau /u_*^2$
in NSBLs collapse to a single curve that is nearly independent of the dimensionless parameters
$\textit{Ro}$
and
$\mu$
. In addition, figure 4 shows that the theoretical prediction of
$\tau /u_*^2$
by both the Nieuwstadt’s model (dashed line) and the present model given by (2.3), (2.8) and (2.14) (solid line) agree very well with the simulation data (filled circles), as well as the field data from Nieuwstadt (Reference Nieuwstadt1984b
) (open diamonds) and Wittich (Reference Wittich1991) (open triangles). This fact reveals that investigating the total turbulent momentum flux profile alone may not be able to determine faithfully the value of the exponent
$\alpha$
. At the same time, it also explains the reason why Shen et al. (Reference Shen, Liu, Lu and Stevens2024) still used the power law of the total turbulent momentum flux with exponent
$\alpha =3/2$
as a reference to validate their simulation database. Therefore, it is necessary to further investigate the streamwise and spanwise turbulent momentum fluxes.
The mean vertical profile of the normalised total turbulent momentum flux
$\tau /u_*^2$
. Filled circles indicate simulation data from table 1; open diamonds indicate field data of Nieuwstadt (Reference Nieuwstadt1984b
); open triangles indicate field data of Wittich (Reference Wittich1991); dashed line indicates theoretical prediction by Nieuwstadt (Reference Nieuwstadt1984b
); solid line indicates theoretical prediction by present model of (2.3), (2.8) and (2.14).
The mean vertical profiles of the normalised (a) streamwise turbulent momentum flux
$\tau _x/u_*^2$
and (b) spanwise turbulent momentum flux
$\tau _y/u_*^2$
. Filled circles indicate simulation data from table 1; open triangles indicate field data of Caughey et al. (Reference Caughey, Wyngaard and Kaimal1979); dashed line indicates theoretical prediction by Nieuwstadt (Reference Nieuwstadt1984b
); solid line indicates: theoretical prediction by present model of (2.3), (2.8) and (2.14).
Figure 5 shows the mean vertical profiles of the streamwise and spanwise turbulent momentum fluxes (
$\tau _x, \tau _y$
) normalised by the scaling
$u_*^2$
. All simulation data in table 1 are shown in the figure, which shows that the profiles of the dimensionless turbulent momentum fluxes
$\tau _x/u_*^2$
and
$\tau _y/u_*^2$
in NSBLs collapse to single curves that are nearly independent of the dimensionless parameters
$\textit{Ro}$
and
$\mu$
. In particular, figure 5(a) shows that the theoretical predictions of
$\tau _x/u_*^2$
by both Nieuwstadt’s model (dashed line) and the present model given by (2.3), (2.8) and (2.14) (solid line) agree very well with the simulation data (filled circles), and they also agree reasonably well with the field data from Caughey et al. (Reference Caughey, Wyngaard and Kaimal1979) (open triangles). This is in our expectation since the streamwise turbulent momentum flux
$\tau _x$
contributes most to the total turbulent momentum flux
$\tau$
. On the other hand, figure 5(b) shows that the theoretical prediction of
$\tau _y/u_*^2$
by Nieuwstadt’s model (dashed line) significantly underestimates the magnitude of
$\tau _y/u_*^2$
obtained by simulations (filled circles). In contrast, the theoretical prediction of
$\tau _x/u_*^2$
by the present model given by (2.3), (2.8) and (2.14) (solid line) agrees very well with the simulation data (filled circles). This highlights the benefit of analysing simultaneously the various turbulent momentum fluxes
$\tau$
,
$\tau _x$
and
$\tau _y$
, which enables a more accurate determination of the value of the exponent
$\alpha$
.
The mean vertical profiles of the normalised (a) streamwise velocity deficit and (b) spanwise velocity deficit. Filled circles indicate simulation data from table 1; dashed line indicates theoretical prediction by Nieuwstadt (Reference Nieuwstadt1984b ); solid line indicates theoretical prediction by the present model of (2.5), (2.8) and (2.14).
4.2. Mean wind deficits
To further evaluate the performance of our proposed model, figure 6 shows the mean vertical profiles of the streamwise and spanwise velocity deficits
$U-U_g, V-V_g$
normalised by the scaling
$u_*^2/(hf)$
. The filled circles are the simulation data of table 1, and the dashed line is the theoretical prediction by Nieuwstadt (Reference Nieuwstadt1984b
). The simulation data show good collapse throughout most of the boundary layer, except very close to the surface, indicating that the geostrophic wind deficits
$U-U_g, V-V_g$
within the Ekman layer scale with
$u_*^2/(hf)$
. Note that this result directly leads to the conclusion that the wind gradients
$\textrm {d} U/\textrm {d} z, \textrm {d} V/\textrm {d} z$
scale with
$u_*^2/(h^2f)$
, which is an essential assumption made by Shen et al. (Reference Shen, Liu, Lu and Stevens2024) to derive their analytical expressions of the geostrophic drag law coefficients of NSBLs. Recall that Nieuwstadt’s model assumes the exponent
$\alpha =3/2$
, while the boundary layer height
$h$
in figure 6 is calculated using (2.4) with the exponent
$\alpha =2$
. As a consequence, Nieuwstadt’s model underestimates the boundary layer height
$h$
, and has significant deviations from the simulated velocity deficit profiles in both the surface layer and Ekman layer (dashed line). Figure 6 also shows that the present model of (2.5), (2.8) and (2.14) accurately captures the wind deficit profiles in the Ekman layer and the upper part of the surface layer (solid line), which is due to a more accurate parametrisation of the flux Richardson number
$R_{\!f}$
(see figure 2
b). Therefore, the present model can accurately predict all components of the turbulent momentum flux profiles and the velocity deficit profiles within the Ekman layer of NSBLs.
We note that the streamwise velocity has a log-linear form in the surface layer, while our model and that of Nieuwstadt (Reference Nieuwstadt1984b
) reduce to a linear profile due to the assumption of nearly constant
$R_{\!f}$
in that region (Garratt Reference Garratt1992). Additionally, we note that the term
$(1-z/h)^{i \gamma }= \cos \, [\gamma \ln (1-z/h)] + \textrm{i} \sin \, [\gamma \ln (1-z/h)]$
in (2.5) indicates rapid oscillations near the top of the boundary layer. Although most of these oscillations are suppressed by the term
$(1-z/h)^{\alpha -1}$
in (2.5), one oscillation is still visible in the theoretical curves (figure 6
a,b). Despite these small differences, the overall performance of our model in the whole Ekman layer is very good and helps us to better understand the flow physics of quasi-steady, barotropic and horizontally homogeneous NSBLs.
The MAE of the normalised (a) streamwise velocity deficit and (b) spanwise velocity deficit. Dashed line indicates theoretical prediction by Nieuwstadt (Reference Nieuwstadt1984b ); solid line indicates theoretical prediction by the present model of (2.5), (2.8) and (2.14).
To quantify the improved performance of our proposed model compared to Nieuwstadt’s model, figure 7 shows the mean absolute error (MAE) of the wind deficit profiles in figure 6. Here, the MAE is calculated as
\begin{equation} \text{MAE} (X) = \left|X_{{model}} - \overline X \right|, \quad \overline X = \frac {1}{n} \sum _{i=1}^n X_{{LES}}^i, \end{equation}
where
$X = hf (U-U_g)/u_*^2$
and
$hf (V-V_g)/u_*^2$
,
$\overline X$
is the sample mean (i.e. the average of all cases) profile of
$X$
(see figure 6), the subscripts ‘LES’ and ‘model’ denote the results obtained from numerical simulation and model prediction, respectively,
$i$
is the case number, and
$n=16$
is the total number of all simulated cases (see table 1). Figure 7 shows clearly that the MAE obtained by the present model is generally much smaller than in Nieuwstadt’s model. In particular, above the surface layer (
$z/h \gtrsim 0.1$
), the maximum MAEs of normalised streamwise and spanwise velocity deficits obtained by the present model are approximately 0.10 and 0.12, respectively, whereas those obtained by Nieuwstadt’s model are approximately 0.45 and 0.51. This quantitatively determines the error of the model prediction and the improved performance of our proposed model compared to Nieuwstadt’s model.
The universal Ekman layer structure has important potential impact on meteorological applications. For instance, wind turbines are becoming larger, and the average turbine hub height is now approximately 100 m (see table 1 for the typical NSBL height, and GWEC (2024) for the typical turbine hub height). Consequently, modern turbines will increasingly operate in less-explored atmospheric regions, emphasising the need for a more thorough understanding of atmospheric physics beyond the surface layer. Thus to explore the wind characteristics in high latitudes, figure 8(a) shows the mean vertical profiles of the wind speed in the Ekman layer obtained from simulations (open circles) and the combination of (2.5), (2.8) and (2.14)–(2.19) (solid line). Clearly, the theoretically predicted wind speed profiles agree well with the simulation results in the Ekman layer and the upper part of the surface layer. Therefore, the wind speed profile in high altitude of quasi-steady, barotropic and horizontally homogeneous NSBLs can be estimated theoretically and accurately from
$u_*, L_s, z_0$
that can be conveniently measured near the surface. In addition, figure 8(b) shows the mean wind speed profile in the surface layer obtained from simulations (filled circles) and MOST of (2.20) with
$c_\psi =4.8$
(dashed line), a value suggested by Beare et al. (Reference Beare2006) based on atmospheric observations and adopted in the wall-modelled LES by Shen et al. (Reference Shen, Liu, Lu and Stevens2024), as well as MOST of (2.20) with
$c_\psi =8.9$
(solid line), a value predicted by the present theory (see (2.16)). The simulation results at the first grid point show good agreement with the dashed line, as expected, since this point reflects the imposed boundary condition. At higher elevations, deviations become apparent due to the log-layer mismatch commonly observed in wall-modelled LES (Mason & Thomson Reference Mason and Thomson1992; Cabot & Moin Reference Cabot and Moin1999; Gadde, Stieren & Stevens Reference Gadde, Stieren and Stevens2021). The elevated value of
$c_\psi$
in the present model results from fitting the model parameters to LES data affected by this mismatch. Overall, the comparison shows that the model captures the key characteristics of the MOST profile in the surface layer well.
The normalised mean wind speed profiles of NSBLs in (a) the Ekman layer and (b) the surface layer. Open and filled circles indicate simulation data from table 1; solid line indicates theoretical prediction by (a) the present model of (2.5), (2.8), and (2.14)–(2.19), and (b) MOST of (2.20) with
$c_\psi =8.9$
; dashed line indicates MOST with
$c_\psi =4.8$
.
4.3. Eddy viscosity profile
Before concluding the paper, we remark again that the derivation present in § 2 results in the same expression for the eddy viscosity
$K_m$
in the Ekman layer as that in Basu & Holtslag (Reference Basu and Holtslag2023). With
$h^2 \sim L_f L_s = L_s u_*/f$
(Shen et al. Reference Shen, Liu, Lu and Stevens2024), we find that the eddy viscosity
$K_m$
scales with
$u_* L_s$
, as proposed before by e.g. Zilitinkevich et al. (Reference Zilitinkevich, Esau and Baklanov2007). Given the prevalence of formulations for
$K_m$
and the fact that such expressions are often used to obtain mean wind profiles (see § 1), it may be very relevant to evaluate the performance of (2.7). Therefore, the mean vertical profile of the eddy viscosity
$K_m$
normalised by
$h^2 f$
is shown in figure 9(a), where the filled circles are the simulation data from table 1, the dashed line is the theoretical prediction of (2.7) with
$\alpha =3/2$
, and the solid line is the theoretical prediction of (2.7) with
$\alpha =2$
. Note that the eddy viscosity
$K_m$
from the simulation data is calculated from (2.6), i.e.
$K_m = \{(\tau _x^2+\tau _y^2)/[(\textrm {d} U/\textrm {d} z)^2+(\textrm {d} V/\textrm {d} z)^2]\}^{1/2}$
, and the boundary layer height
$h$
involved in the theoretical prediction with
$\alpha =3/2$
(dashed line) has been adjusted by the factor
$(1-0.05^{1/2})/(1-0.05^{2/3})$
for fair comparisons (see also the discussion on figure 2). The figure shows clearly that the prediction of (2.7) with
$\alpha =2$
(solid line) matches well with the
$K_m$
profiles from the simulation data above the surface layer
$z/h \gtrsim 0.1$
, while that with
$\alpha =3/2$
(dashed line) overestimates the results in the region
$z/h \lesssim 0.6$
. This result also indicates that the exponent
$\alpha =2$
should be more accurate. Furthermore, figure 9(b) shows the mean vertical profile of
$K_m/(\kappa u_* L_s)$
in the surface layer. It shows clearly that all simulation data (filled circles) collapse to a single curve of the dimensionless parameter
$z/L_s$
described by (2.9) (solid line). This confirms that the normalised eddy viscosity
$K_m/(\kappa u_* L_s)$
profile in the surface layer follows MOST.
The mean vertical profile of the normalised eddy viscosity (a)
$K_m/(h^2 f)$
in the Ekman layer, and (b)
$K_m/(\kappa u_* L_s)$
in the surface layer. Filled circles indicate simulation data from table 1; dashed line indicates theoretical prediction by (2.7) with
$\alpha =3/2$
; solid line indicates theoretical prediction by (a) (2.7) with
$\alpha =2$
, and (b) MOST of (2.9).
5. Conclusions
In this study, we investigate the mean flow properties of the quasi-steady, barotropic and horizontally homogeneous NSBLs through theoretical analysis. By solving the Ekman equations with the eddy viscosity approach and a novel parametrisation of the flux Richardson number
$R_{\!f}$
, we derive an analytical model for the vertical profiles of the turbulent momentum fluxes and geostrophic wind deficits. The model reveals that the mean vertical profiles of turbulent momentum fluxes
$\tau , \tau _x, \tau _y$
and wind deficits
$U-U_g, V-V_g$
reduce to a set of universal curves when normalised with the scales
$u_*^2$
and
$u_*^2/(hf)$
, respectively. To validate and evaluate the model performance, we use the high-fidelity database of the quasi-steady, barotropic and horizontally homogeneous NSBLs generated by Shen et al. (Reference Shen, Liu, Lu and Stevens2024) using advanced LES, and focus on the cases under moderate thermal stratification, i.e. with the Kazanski–Monin parameter
$\mu \in [34.0, 193.3]$
. The theoretical predictions agree well with simulation results and atmospheric measurement data, confirming the validity of the proposed model. This model, when combined with the geostrophic drag law (Zilitinkevich & Esau Reference Zilitinkevich and Esau2005; Shen et al. Reference Shen, Liu, Lu and Stevens2024) and the boundary layer parametrisation (Zilitinkevich Reference Zilitinkevich1972; Basu & Holtslag Reference Basu and Holtslag2023), enables us to predict theoretically and accurately the wind speed in high altitudes or the Ekman layer of NSBLs, which has also been confirmed by the simulation data. Our study advances the fundamental understanding of NSBLs with direct relevance for meteorological applications.
The present model accurately reproduces the wind speed profile in the Ekman layer, but does not capture the log-linear structure characteristic of the surface layer. This limitation arises because the flux Richardson number
$R_{\!f}$
of (2.13) and the eddy viscosity
$K_m$
of (2.7) are valid only within the Ekman layer, whereas the surface layer is governed by Monin–Obukhov similarity theory (MOST; Monin & Obukhov Reference Monin and Obukhov1954; Brost & Wyngaard Reference Brost and Wyngaard1978; Wyngaard Reference Wyngaard2010). Asymptotic matching between the surface and Ekman layers is crucial not only for deriving geostrophic drag laws (e.g. Shen et al. Reference Shen, Liu, Lu and Stevens2024) but also for obtaining accurate velocity profiles (e.g. Kadantsev et al. Reference Kadantsev, Mortikov and Zilitinkevich2021; Narasimhan et al. Reference Narasimhan, Gayme and Meneveau2024). Therefore, by combining with MOST, the presented model can also effectively reproduce the essential features of the flux Richardson number
$R_{\!f}$
, the eddy viscosity
$K_m$
, and the wind speed profile within the surface layer.
Finally, we note that using a small cooling rate at the surface with strong geostrophic forcing means that we are considering a weakly stable ABL regime. The fact that the quasi-steady state is reached only after a long time suggests that the simulated mean profiles may occur more frequently in polar regions, where subsidence could have a significant impact (Bon et al. Reference Bon, Cal and Meyers2024). Further research is required to determine the applicability of the approach to the very stable ABL regime (Chinita, Matheou & Miranda Reference Chinita, Matheou and Miranda2022) and scenerios with subsidence (Bon et al. Reference Bon, Cal and Meyers2024).
Acknowledgements
We acknowledge the anonymous referees whose insightful comments and suggestions have substantially enhanced the paper’s overall quality. The numerical calculations in this paper have been done on the supercomputing system in the Supercomputing Center of the University of Science and Technology of China and the Hefei Advanced Computing Center. We also acknowledge the Dutch national e-infrastructure of SURFsara, a subsidiary of the SURF cooperation, the collaborative ICT organisation for Dutch education and research.
Funding
This work was supported by the National Natural Science Foundation of China (grant nos 12388101, 12272371) and the National Natural Science Fund for Excellent Young Scientists Fund Programme (Overseas), and has received funding from the European Research Council under the European Union’s Horizon Europe Programme (grant no. 101124815).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
All study data are included in this paper.
