2.1. Governing equations

We use LES to simulate the NSBL flow over an infinite flat surface with homogeneous roughness. We integrate the spatially filtered continuity equation, momentum equation and potential temperature equation (Albertson Reference Albertson1996; Albertson & Parlange Reference Albertson and Parlange1999; Liu et al. Reference Liu, Gadde and Stevens2021a; Liu, Gadde & Stevens Reference Liu, Gadde and Stevens2021b)

(2.1)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}} = 0, \end{gather}$$

(2.2)$$\begin{gather}\partial_t \tilde{\boldsymbol{u}} + \tilde{\boldsymbol{\omega}} \times \tilde{\boldsymbol{u}} = f \boldsymbol{e}_z \times (\boldsymbol{G} - \tilde{\boldsymbol{u}}) + \beta (\tilde{\theta} - \langle \tilde{\theta} \rangle) \boldsymbol{e}_z - \boldsymbol{\nabla} \tilde{p} - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau}, \end{gather}$$

(2.3)$$\begin{gather}\partial_t \tilde{\theta} + \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\theta} ={-} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{q}}, \end{gather}$$

where the tilde denotes spatial filtering, $\langle \boldsymbol {\cdot } \rangle$ represents horizontal averaging, $\tilde {\boldsymbol {u}}$ is the velocity, $\tilde {\boldsymbol {\omega }} = \boldsymbol {\nabla } \times \tilde {\boldsymbol {u}}$ is the vorticity, $\tilde {p}$ is the modified pressure departure from equilibrium, $\tilde {\theta }$ is the potential temperature, $\beta$ is the buoyancy parameter, $\boldsymbol {G}=(U_g, V_g, 0)$ is the geostrophic wind velocity, $\boldsymbol {\tau }$ denotes the deviatoric part of the SGS shear stress and $\boldsymbol {q}$ represents the SGS potential temperature flux. Molecular viscosity is neglected as the Reynolds number of atmospheric flows is very high.

For closure of (2.1)–(2.3), we parameterize the SGS shear stress and potential temperature flux as

(2.4a–c)\begin{equation} \boldsymbol{\tau} ={-} 2 \nu_t \tilde{\boldsymbol{S}}, \quad \boldsymbol{q} ={-} \nu_\theta \boldsymbol{\nabla} \tilde{\theta}, \quad \tilde{\boldsymbol{S}} \equiv \frac{1}{2} [ \boldsymbol{\nabla} \tilde{\boldsymbol{u}} + (\boldsymbol{\nabla} \tilde{\boldsymbol{u}})^T ], \end{equation}

where $\nu _t$ and $\nu _\theta$ are the eddy viscosity and eddy diffusivity, respectively, and $\tilde {\boldsymbol {S}}$ is the resolved strain-rate tensor with the superscript $T$ denoting a matrix transpose. We use the anisotropic minimum dissipation (AMD) model (Abkar & Moin Reference Abkar and Moin2017), which has been extensively validated (Gadde, Stieren & Stevens Reference Gadde, Stieren and Stevens2021), to perform the simulations. In the AMD model, the eddy viscosity and eddy diffusivity are determined so that the energy of the sub-filter scales of the LES solution does not increase with time.

2.2. Numerical method

Our code is an updated version of the one used by Albertson & Parlange (Reference Albertson and Parlange1999). The grid points are uniformly distributed, and the computational planes for horizontal and vertical velocities are staggered in the vertical direction. The first vertical velocity grid plane is located at the ground, and the first streamwise and spanwise velocities and potential temperature grid planes are located at half a grid distance away from the ground (Albertson Reference Albertson1996). In the vertical direction, we use a second-order finite difference method, while we use a pseudo-spectral discretization with periodic boundary conditions in the horizontal directions. Time integration is performed using the second-order Adams–Bashforth method. The projection method is used to ensure the divergence-free condition of the velocity field.

At the top boundary the vertical velocity, the SGS shear stress and the SGS potential temperature flux are enforced to zero, while the potential temperature is imposed by a constant value. In the top 25 % of the domain a Rayleigh damping layer with the damping frequency $0.001\,{\rm s}^{-1}$ determined by trial and error is applied every time step to reduce the effects of gravity waves (Klemp & Lilly Reference Klemp and Lilly1978; Liu & Stevens Reference Liu and Stevens2021). At the bottom boundary, we employ a wall model based on the Monin–Obukhov similarity theory (MOST) for both the velocity field and the potential temperature field (Monin & Obukhov Reference Monin and Obukhov1954; Moeng Reference Moeng1984; Stoll & Porté-Agel Reference Stoll and Porté-Agel2008; Liu & Stevens Reference Liu and Stevens2021)

(2.5a,b)\begin{equation} u_{*} = \frac{ \kappa (\bar{\tilde{u}}^2 + \bar{\tilde{v}}^2)^{1/2} }{ \ln ( z / z_0 ) - \psi_m }, \quad q_{s} = \frac{ \kappa u_* ( \theta_s - \bar{\tilde{\theta}} ) }{ \ln ( z / z_0 ) - \psi_h }, \end{equation}

where the overline denotes filtering at a scale of $2\Delta$ with $\Delta = (\Delta x \Delta y \Delta z)^{1/3}$ the filter scale of (2.1)–(2.3) and $\Delta x, \Delta y, \Delta z$ the grid spacings in the streamwise, spanwise and vertical directions, respectively, and where $\theta _s$ is the potential temperature at the surface, and $\psi _m$ and $\psi _h$ are, respectively, the stability corrections for the momentum and potential temperature fluxes, which for NSBLs can be parameterized as (Businger et al. Reference Businger, Wyngaard, Izumi and Bradley1971; Stull Reference Stull1988; Huang & Bou-Zeid Reference Huang and Bou-Zeid2013; Sullivan et al. Reference Sullivan, Weil, Patton, Jonker and Mironov2016; Abkar & Moin Reference Abkar and Moin2017)

(2.6a–c)\begin{equation} \psi_m ={-}4.8 \kappa z/L_s, \quad \psi_h ={-}7.8 \kappa z/L_s, \quad z/L_s\ge 0. \end{equation}

Note that the definition of $L_s$ in this study is the same as Zilitinkevich & Esau (Reference Zilitinkevich and Esau2005) but different from Liu & Stevens (Reference Liu and Stevens2021), and thus the constant $\kappa$ appears in (2.6a–c). In addition, we note that the MOST, which is the basis of the wall model (2.5a,b), has been confirmed by various direct numerical simulations (e.g. Shah & Bou-Zeid Reference Shah and Bou-Zeid2014) and field experiments (e.g. Businger et al. Reference Businger, Wyngaard, Izumi and Bradley1971). Since the MOST is generally valid in the surface layer of the ABLs in weakly stable regime (Grachev et al. Reference Grachev, Fairall, Persson, Andreas and Guest2005; Katul, Konings & Porporato Reference Katul, Konings and Porporato2011; Stoll et al. Reference Stoll, Gibbs, Salesky, Anderson and Calaf2020), the wall model based on the MOST is used in this study.

2.3. Computational set-up

To obtain simulation results independent of the computational domain, the domain size in each direction has to be several times the size of the largest eddies or the boundary-layer height (Couvreux et al. Reference Couvreux2020). The selected computational domain size is $800\,{\rm m}\times 800\,{\rm m}\times 800\,{\rm m}$ in the streamwise, spanwise and vertical directions, respectively. As shown in table 1, the deepest boundary-layer height is approximately 266 m, which is less than one third of the computational domain size. Since the MOST-based wall model is employed, the lowest grid level $z_1$ should be above the lower limit of the inertial sublayer $z_*$ (Basu & Lacser Reference Basu and Lacser2017). The parameter $z_*$ can be related to the roughness length $z_0$, and their ratio depends on the atmospheric stability. In particular, for stably stratified conditions, the ratio $z_*/z_0$ for the wind profile was found to be significantly lower, approximately in the range of $11 \sim 55$ (Garratt Reference Garratt1983; Basu & Lacser Reference Basu and Lacser2017). Meanwhile, $z_1$ should be below the upper limit of the inertial sublayer, i.e. the Ozmidov length $L_O=(\varepsilon /N^3)^{1/2}$, where $\varepsilon$ is the mean viscous dissipation and $N$ is the Brunt–Väisälä frequency (Dougherty Reference Dougherty1961; Ozmidov Reference Ozmidov1965; Huang & Bou-Zeid Reference Huang and Bou-Zeid2013; Sullivan et al. Reference Sullivan, Weil, Patton, Jonker and Mironov2016). We select the grid resolution as $192\times 192\times 201$ such that the grid is nearly isotropic (Gadde et al. Reference Gadde, Stieren and Stevens2021), and, more importantly, $z_1/z_0\ge 20$ since $z_0 \le 0.1$ m (see table 1). At the same time, this grid also satisfies the grid resolution requirement $z_1/L_O<1$ (see figure 3 below).

Table 1.

Summary of simulations performed on the domain size $800\,{\rm m}\times 800\,{\rm m}\times 800\,{\rm m}$ and $192\times 192\times 201$ grid points, where $h_{0.05}$ is defined at the height where the total shear stress reduces to 5 % of its surface value and $\mu =L_f/L_s$ with the Obukhov length $L_s=-u_*^3/(\beta q_s)$.

The simulations are initialized with a constant potential temperature $\theta =265$ K and a constant geostrophic wind speed $G=8\sim 12\,{\rm m}\,{\rm s}^{-1}$. Turbulence is triggered by adding random perturbations. A random noise term with a magnitude of $3\,\%$ of the geostrophic wind speed is added to velocities below $50$ m, and for the potential temperature a noise term with an amplitude of $0.1$ K is added. The reference temperature is set to $\theta _0=263.5$ K. The latitude $\phi$ varies between $30^\circ$ and $73^\circ$, and the surface cooling rate is set to $C_r=0.05\sim 0.35\,{\rm K}\,{\rm h}^{-1}$. The roughness length is $z_0=0.001\sim 0.1$ m. To remove the effects of slowly decaying inertial oscillations on the top of the domain (Van de Wiel et al. Reference Van de Wiel, Moene, Steeneveld, Baas, Bosveld and Holtslag2010; Maronga & Li Reference Maronga and Li2022), the initial transient period lasts for 1.5 inertial time periods and the statistics are collected over the following one inertial period, i.e. $ft \in [3{\rm \pi}, 5{\rm \pi} ]$ (see figure 2). This ensures the flow has reached its quasi-stationary state and the changes of $u_*$ and $\alpha _0$ are less than 5 % of their mean values (see figure 2). We note that many time steps ($\sim O(10^7)$) are required to resolve the large range of relevant time scales in the problem, spanning from small-scale turbulence to large-scale inertial oscillations. A summary of these simulations is given in table 1, where the cases are arranged such that the value of $\phi$ increases monotonically from upper to lower rows.

Figure 2.

The time history of the domain-averaged (a) friction velocity $u_*$ and (b) cross-isobaric angle $\alpha _0$. The solid lines are the simulation cases of table 1, and the dashed line denotes the beginning of the time-averaging window.

We remark that the simulation input parameters of case 20 are the same as the standard Global Energy and Water Cycle Experiment Atmospheric Boundary Layer Study case (Beare et al. Reference Beare2006; Liu & Stevens Reference Liu and Stevens2021), with the only exception being that the free-atmosphere lapse rate has been set as zero to ensure the flow belongs to the NSBL regime (Maronga et al. Reference Maronga, Knigge and Raasch2020). Besides, we remark that the simulations of Zilitinkevich & Esau (Reference Zilitinkevich and Esau2005) that will be used for comparison were performed for $ft=2.1{\rm \pi}$, with time averaging over the period $ft\in [1.3{\rm \pi}, 2.1{\rm \pi} ]$. As shown in figure 2, a time averaging should start around $ft \approx 3{\rm \pi}$ to ensure reliable determination of the friction velocity $u_*$ and the cross-isobaric angle $\alpha _0$. This is the main reason for the large scatter of the simulation data of Zilitinkevich & Esau (Reference Zilitinkevich and Esau2005). In contrast, the simulation data obtained from our LES with much longer simulation time show very limited scatter (see below).

It is well known that sufficient grid resolution is important in simulating stable ABLs with low turbulence levels (Huang & Bou-Zeid Reference Huang and Bou-Zeid2013; Maronga & Li Reference Maronga and Li2022). Therefore, a quantitative investigation into the appropriateness of the chosen grid resolution is performed. Figure 3 shows that for all simulated cases the ratio $z_1/L_O <1$, which justifies the use of the MOST-based wall model. Similar to Sullivan et al. (Reference Sullivan, Weil, Patton, Jonker and Mironov2016), figure 3 shows that the Ozmidov length $L_O$ decreases as the cooling rate $C_r$ increases. Thus, if a much larger cooling rate, e.g. $C_r=2.5\,{\rm K}\,{\rm h}^{-1}$ (Huang & Bou-Zeid Reference Huang and Bou-Zeid2013), is considered, a higher grid resolution may be needed. However, for the cooling rate range $C_r=0.05\sim 0.35\,{\rm K}\,{\rm h}^{-1}$ considered here, the chosen grid resolution meets the required standards.

Figure 3.

The ratio $z_1/L_O$ vs the Kazanski–Monin parameter $\mu$. The simulation data are from table 1 and the Ozmidov length $L_O$ is calculated at $z=z_1$ (the lowest grid level).

To further illustrate the effect of grid resolution on simulation results, figure 4 shows the mean vertical profiles of the wind speed, potential temperature, total shear stress and potential temperature flux for case 2 in table 1, which has the largest value of $\mu =193.3$ and hence may have the strongest dependence on the grid resolution. The figure indicates that the magnitudes of the total shear stress and potential temperature flux in the boundary layer, as well as the boundary-layer height, decrease slightly as the grid resolution increases. This lack of strictly grid convergence in LES of the stable boundary layer is well known and remains partially unresolved (Maronga & Li Reference Maronga and Li2022). Nevertheless, figure 4 shows that the mean vertical profiles with medium and dense grid resolutions do not change significantly. This finding confirms the appropriateness of the adopted grid resolution of $192\times 192\times 201$ in this study.

Figure 4.

The grid dependence of the mean vertical profiles of the (a) wind speed $\sqrt {U^2+V^2}$, (b) potential temperature $\varTheta$, (c) total shear stress $\tau$ and (d) potential temperature flux $q$. The simulated case is case 2 with $\mu =193.3$ in table 1.

3.1. Boundary-layer height

In stable ABL flows the profile of the magnitude of total turbulent shear stress $\tau$ is often expressed as follows (Nieuwstadt Reference Nieuwstadt1984; Sorbjan Reference Sorbjan1986; Basu & Holtslag Reference Basu and Holtslag2023):

(3.1)\begin{equation} \frac{\tau}{u_*^2} = \left( 1-\frac{z}{h} \right)^\alpha. \end{equation}

We remark that (3.1) was proposed under quasi-steady, barotropic and horizontally homogeneous conditions (Nieuwstadt Reference Nieuwstadt1984). Thus, it might not be valid for unsteady (Mahrt Reference Mahrt2014; Momen & Bou-Zeid Reference Momen and Bou-Zeid2017), baroclinic (Arya & Wyngaard Reference Arya and Wyngaard1975; Floors, Peña & Gryning Reference Floors, Peña and Gryning2015; Momen et al. Reference Momen, Bou-Zeid, Parlange and Giometto2018) or surface heterogeneous (Schmid et al. Reference Schmid, Lawrence, Parlange and Giometto2019; Cooke, Jerolmack & Park Reference Cooke, Jerolmack and Park2024) flows. In addition, we remark that many different definitions of the stable ABL height $h$ are available in literature, such as (a) the layer through which surface-based turbulence extends, (b) the top of the layer with downward potential temperature flux, (c) the height of the low-level jet and (d) the top of the potential temperature inversion layer (Vickers & Mahrt Reference Vickers and Mahrt2004).

In this study, we adopt definition (a). Specifically, the boundary-layer height $h$ is defined as the height at which the total shear stress first vanishes. In particular, $h$ is proportional to the height $h_{0.05}$, where the total shear stress reduces to 5 % of its surface value (van Dop & Axelsen Reference van Dop and Axelsen2007; Liu et al. Reference Liu, Gadde and Stevens2021b)

(3.2)\begin{equation} \frac{h_{0.05}}{h} = 1 - 0.05^{{1}/{\alpha}}. \end{equation}

By utilizing his local scaling hypothesis and the constant assumption of the flux Richardson number, Nieuwstadt (Reference Nieuwstadt1984) determined analytically that the exponent $\alpha = 3/2$. Note that the value of $\alpha$ may not be a universal constant and other values of $\alpha$ have also been reported in the literature (Caughey, Wyngaard & Kaimal Reference Caughey, Wyngaard and Kaimal1979; Yokoyama, Gamo & Yamamoto Reference Yokoyama, Gamo and Yamamoto1979; Lacser & Arya Reference Lacser and Arya1986; Sorbjan Reference Sorbjan1986; Lenschow et al. Reference Lenschow, Li, Zhu and Stankov1988; Byun Reference Byun1991; Delage Reference Delage1997).

Figure 5 shows the profile of normalized total turbulent shear stress in NSBLs. The open circles are the simulation data from table 1, the open triangles are the field data taken from Nieuwstadt (Reference Nieuwstadt1984), the open inverted triangles are the field data of Wittich (Reference Wittich1991) and the solid line is the prediction of (3.1) with $\alpha =3/2$. The good agreement of simulation data with field measurements validates our simulation results. Furthermore, our simulations indicate that the 3/2-law of total turbulent shear stress proposed by Nieuwstadt (Reference Nieuwstadt1984) describes the data across the entire parameter range considered here.

Figure 5.

The profile of normalized total turbulent shear stress. Open circles: simulation data from table 1; open triangles: field data of Nieuwstadt (Reference Nieuwstadt1984); open inverted triangles: field data of Wittich (Reference Wittich1991); solid line: prediction of (3.1) with $\alpha =3/2$.

For NSBLs with large values of $\mu$, the parameterization of the boundary-layer height $h$ given by (1.4) can be approximated as

(3.3)\begin{equation} \frac{h}{L_f} = \frac{c_s}{ \sqrt{\mu} }, \end{equation}

where $c_s$ is the Zilitinkevich constant (Derbyshire Reference Derbyshire1990), which depends on the definition of the boundary-layer height $h$. Equation (3.3) was originally proposed by Zilitinkevich (Reference Zilitinkevich1972) based on boundary-layer scaling arguments and implies that the boundary-layer height $h$ scales as $\sqrt {L_f L_s}$. Based on this definition, Zilitinkevich (Reference Zilitinkevich1989b) summarized $c_s$ from several studies and found it to vary between 0.55 and 1.58. This large variation affects the values of other empirical constants in the GDL formulation of Zilitinkevich & Esau (Reference Zilitinkevich and Esau2005). As shown below, this is not the case for our formulation, where only the scaling of (3.3) matters while the value of $c_s$ is irrelevant. This is also the case for the value of $\alpha$ since it can be absorbed into $c_s$, see (3.2) and (3.3). Recently, Basu & Holtslag (Reference Basu and Holtslag2023) derived two different eddy viscosity profile formulations, both leading to the boundary-layer height parameterization of Zilitinkevich (Reference Zilitinkevich1972). The two derivations give $c_s=0.922$ and $c_s=0.658$, respectively, and they argued that the Zilitinkevich constant $c_s=0.658$ seems to be a better option. More recently, Narasimhan et al. (Reference Narasimhan, Gayme and Meneveau2024) also reported $c_s=0.78$ based on their LES data.

Figure 6 compares the dimensionless boundary-layer height $h/L_f$ obtained from the present simulations of table 1 (filled circles), the simulation data of Zilitinkevich et al. (Reference Zilitinkevich, Esau and Baklanov2007) (open circles) and the prediction of (3.3) with $c_s=0.62$, which is determined using a least-squares fitting procedure based on the present LES database. All symbols collapse to a single curve, indicating that the dimensionless boundary-layer height $h/L_f$ indeed only depends on the Kazanski–Monin parameter $\mu$. The good agreement between the simulation data and theoretical prediction confirms the validity of the boundary-layer height parameterization of (3.3) over a much wider range of $\mu$ than Zilitinkevich et al. (Reference Zilitinkevich, Esau and Baklanov2007) considered. Note that the simulation data of Zilitinkevich et al. (Reference Zilitinkevich, Esau and Baklanov2007), which have a smaller range of $\mu$ than the present study, are very close to our theoretical model and simulation results.

Figure 6.

The dimensionless boundary-layer height $h/L_f$ vs the Kazanski–Monin parameter $\mu$. Filled circles: simulation data of table 1; open circles: simulation data of Zilitinkevich et al. (Reference Zilitinkevich, Esau and Baklanov2007); solid line: prediction of (3.3) with $c_s=0.62$ obtained by a least-squares fit.

3.2. Wind gradient profiles

Recall that we focus on the quasi-steady and barotropic NSBLs in weakly stable regime, whose dynamics is governed by the Ekman equations

(3.4a,b)\begin{equation} \frac{{\rm d} \tau_x}{{\rm d} z} ={-}f (V-V_g), \quad \frac{{\rm d} \tau_y}{{\rm d} z} = f (U-U_g), \end{equation}

where $\tau _x$ and $\tau _y$ are the horizontally averaged streamwise and spanwise turbulent shear stresses. The vertical derivative of (3.4a,b) is

(3.5a,b)\begin{equation} \frac{{\rm d^2}}{{\rm d} \xi^2} \left( \frac{ \tau_x}{u_*^2} \right) ={-} \frac{h^2 f}{u_*^2} \frac{{\rm d} V}{{\rm d} z}, \quad \frac{{\rm d^2}}{{\rm d} \xi^2} \left( \frac{ \tau_y}{u_*^2} \right) = \frac{h^2 f}{u_*^2} \frac{{\rm d} U}{{\rm d} z}, \end{equation}

where $\xi =z/h$ is the normalized coordinate. Since the total turbulent shear stresses nearly collapse onto a single curve as a function of $z/h$, see (3.1) and figure 5, it is intuitive to conjecture that the wind gradients $\textrm {d} U/\textrm {d} z$ and $\textrm {d} V/\textrm {d} z$ normalized with $u_*^2/(h^2 f)$ should also approximately collapse. To confirm this conjecture, figure 7 shows the streamwise and spanwise velocity gradients normalized by $u_*^2/(h^2 f)$. Both tend to collapse to a single curve and are $O(1)$ above the surface layer. This indicates that $u_*^2/(h^2 f)$ is indeed the appropriate scaling for the wind gradients above the surface layer (Lin & Segel Reference Lin and Segel1988). Additionally, neither $\textrm {d} U/\textrm {d} z$ nor $\textrm {d} V/\textrm {d} z$ shows a uniform region, i.e. a region independent of $z$. This indicates that a $z$-less stratification layer is absent in NSBLs in the studied parameter range of $\mu$ and $Ro$.

Figure 7.

The vertical profiles of the normalized (a) streamwise and (b) spanwise velocity gradients. The solid line is the corresponding sample mean profile.

To quantify the effectiveness of the scaling $u_*^2/(h^2 f)$, figure 8 shows the mean absolute error (MAE) of the wind gradient profiles in figure 7. Here, the MAE is calculated as

(3.6)\begin{equation} \textrm{MAE} (X) = \frac{1}{n} \sum_{i=1}^n |X_i - \bar{X}|, \end{equation}

where $X=(h^2 f/u_*^2)(\textrm {d} U/\textrm {d} z)$ and $(h^2 f/u_*^2)(\textrm {d} V/\textrm {d} z)$, $\bar {X}$ is the sample mean (i.e. the average of all cases) profile of $X$ (see figure 7), $i$ is the case number and $n=20$ is the total number of all simulated cases (see table 1). Figure 8 shows clearly that the MAEs of wind gradients normalized by the scaling $u_*^2/(h^2 f)$ are $O(0.1)$ above the surface layer. This quantitatively confirms the effective collapse of the wind gradients by the scaling $u_*^2/(h^2 f)$ in the studied parameter range.

Figure 8.

The MAE of (3.6) of the normalized wind gradients. The dashed line denotes the mean location of the surface-layer top.

We remark that Zilitinkevich & Esau (Reference Zilitinkevich and Esau2005) assumed that the wind gradients $\textrm {d} U/\textrm {d} z$ and $\textrm {d} V/\textrm {d} z$ above the surface layer scale as $u_*/L_s$ and $h^2 f/L_s^2$, respectively. For NSBLs with large values of $\mu$, we have shown that the boundary-layer height $h$ scales as $\sqrt {L_f L_s}$, see (3.3) and figure 6. Under this condition, the scalings $u_*/L_s$, $h^2 f/L_s^2$ and $u_*^2/(h^2 f)$ are equivalent since

(3.7a,b)\begin{equation} \frac{u_*/L_s}{u_*^2/(h^2 f)} = \frac{h^2}{L_f L_s} = c_s^2, \quad \frac{h^2 f/L_s^2}{u_*^2/(h^2 f)} = \frac{h^4}{L_f^2 L_s^2} = c_s^4. \end{equation}

However, choosing the scaling $u_*^2/(h^2 f)$ may still be more convenient in practice, as it is directly derived from the vertical derivatives of the Ekman equations (3.5a,b) and hence ensures the normalized wind gradients above the surface layer are $O(1)$. In addition, for smaller values of $\mu$, the difference between these scalings may exist since the boundary-layer height $h$ ceases to scale as $\sqrt {L_f L_s}$ but scales as $L_f$.

4.1. Analytical expression of A

We first determine the analytical expression of the GDL coefficient $A$. In the surface layer, of which the height is typically approximately 10 % of the boundary-layer height $h$ (Basu & Lacser Reference Basu and Lacser2017), the mean streamwise velocity $U$ is assumed to be described by the MOST (Monin & Obukhov Reference Monin and Obukhov1954; Businger et al. Reference Businger, Wyngaard, Izumi and Bradley1971)

(4.1)\begin{equation} \frac{\kappa U}{u_*} = \ln \left( \frac{z}{z_0} \right) + c_\psi \frac{z}{L_s}, \end{equation}

where $c_\psi$ is an empirical constant and, according to (2.6a–c), $c_\psi = 4.8 \kappa$ in this study. Note that (4.1) is also the basis of the wall model employed in the current LES.

Above the surface layer or in the Ekman layer, figure 7(a) shows that the streamwise velocity gradient $\textrm {d} U/\textrm {d} z$ scales as $u_*^2/(h^2 f)$. Thus, we can assume

(4.2)\begin{equation} \kappa \frac{\textrm{d} U}{\textrm{d} z} = \frac{u_*^2}{ h^2 f} f_u(\xi), \quad \xi=\frac{z}{h}, \end{equation}

where $f_u$ is independent of $Ro$ and $\mu$. Integrating (4.2) from below to the top of the boundary layer and noting that $L_f=u_*/f$, we find

(4.3)\begin{equation} \frac{\kappa }{u_*} (U_g - U) = \frac{ L_f }{h} \int_\xi^1 f_{u} (\xi') \,{\rm d} \xi'. \end{equation}

We further assume the mean streamwise velocity $U$ given by (4.1) and (4.3) matches at the interface $\xi =L_s/(c_1 h)$ or $z=L_s/c_1$, where $c_1$ is an empirical constant. Thus, by substituting (4.1) into (4.3), we find

(4.4)\begin{equation} \frac{\kappa U_g}{u_*} = \ln \left( \frac{L_s}{c_1 z_0} \right) + \frac{c_\psi}{c_1} + \frac{ L_f }{h} \int_{{L_s}/{c_1 h} }^1 f_{u} (\xi') \,{\rm d} \xi'. \end{equation}

Finally, substituting (1.1a) and (3.3) into (4.4) and noting that $\mu =L_f/L_s$, we obtain

(4.5)\begin{equation} A = \ln \mu - a_1 \sqrt{\mu} + a_2, \end{equation}

where

(4.6)\begin{equation} a_1 = \frac{1}{c_s} \int_{{1}/{c_1 c_s \sqrt{\mu}}}^1 f_{u} (\xi')\, \textrm{d} \xi', \quad a_2= \ln c_1 - \frac{c_\psi}{c_1}. \end{equation}

Without knowing the form of $f_u$, the analytical expression of $a_1$ and thus of $A$ cannot be obtained. However, figures 6, 7(a) and 9 indicate that $c_s = O(1)$, $c_1=O(1)$ and $f_u=O(1)$ above the surface layer. Thus, $\textrm {d} a_1/\textrm {d} \mu = f_u / (2 c_1 c_s^2 \mu ^{3/2}) = O(\mu ^{-3/2})= O(10^{-3})$ when $\mu =O(10^2)$, indicating that $a_1$ depends only very weakly on $\mu$ and hence can be assumed as a constant for simplicity. In this sense, the analytical expression of $A$ is given by (4.5), which involves only two empirical constants $(a_1, c_1)$ or equivalently $(a_1, a_2)$. Note that (4.5) has the same structure as equation (24) of Zilitinkevich & Esau (Reference Zilitinkevich and Esau2005). However, in our derivation we divide the boundary layer only into the surface layer and the Ekman layer, and requiring a $z$-less stratification layer is unnecessary. This makes the physics underlying the GDL more clear: it is a direct result of the asymptotic matching of the velocity profile between the surface layer and the Ekman layer. In addition, by making a comparison between (1.2a) and (4.5), we see the term $\ln (a_0 + m_A)$ in (1.2a) should be $a_0 + \ln m_A$. In this way, (1.2a) in the limit $\mu \gg 10$ approaches $A = \ln \mu - a c_s \sqrt {\mu } + a_0$, which is the same as (4.5) if $a c_s =a_1$ and $a_0 =a_2$.

Figure 9.

The GDL coefficient $A$ vs the Kazanski–Monin parameter $\mu$. Filled circles: simulation data of table 1; open circles: simulation data of Zilitinkevich & Esau (Reference Zilitinkevich and Esau2005); dashed line: prediction of (1.2a) proposed by Zilitinkevich & Esau (Reference Zilitinkevich and Esau2005); solid line: prediction of (4.5) with $a_1=0.72$, $c_1=2.25$ obtained by a least-squares fitting procedure.

4.2. Analytical expression of B

To determine the analytical expression of $B$, we recall that

(4.7)\begin{equation} \frac{{\rm d} \tau_y}{{\rm d} z} = f (U-U_g), \end{equation}

where $\tau _y$ is the spanwise component of the total shear stress tensor. We focus on the surface layer. Then, by substituting (1.1a) and (4.1) into (4.7) we obtain

(4.8)\begin{equation} \frac{\kappa}{f u_*} \frac{{\rm d} \tau_y}{{\rm d} z} = A + \ln \left( \frac{z}{L_f} \right) + c_\psi \frac{z}{L_s}. \end{equation}

The bottom boundary condition of (4.8) is $\tau _y(0)=0$. Integrating (4.8) from 0 to $z$, one obtains

(4.9)\begin{equation} \frac{\kappa \tau_y}{f u_*} = (A-1) z + z \ln \left( \frac{z}{L_f} \right) + c_\psi \frac{z^2}{2L_s}. \end{equation}

In the surface layer the eddy viscosity approach is assumed to be valid, such that (Zilitinkevich & Esau Reference Zilitinkevich and Esau2005)

(4.10)\begin{equation} \tau_y = K_m \frac{{\rm d} V}{{\rm d} z}, \quad K_m = \kappa u_* l, \quad \frac{1}{l} = \frac{1}{z} + \frac{c_\psi}{L_s}, \end{equation}

where $V$ is the mean spanwise velocity. We remark that there is an extensive discussion in the literature about the eddy viscosity approach used here (see, e.g. Constantin & Johnson Reference Constantin and Johnson2019), where the eddy viscosity $K_m$ can be a function of height (Grisogono Reference Grisogono1995), stability (Nieuwstadt Reference Nieuwstadt1983) and baroclinicity (Berger & Grisogono Reference Berger and Grisogono1998). In these studies, the velocity profiles were derived based on assumed $K_m$ profiles. However, correct $K_m$ profiles in the whole boundary layer are very difficult to obtain and sometimes the eddy viscosity approach itself might become invalid in some regions of the boundary layer (Liu & Stevens Reference Liu and Stevens2022; Liu et al. Reference Liu, Gadde and Stevens2023). Because the eddy viscosity approach is generally valid in the surface layer, it is adopted here. In particular, in the surface layer $K_m$ in (4.10) is determined from (4.1) and the approximation $K_m {\rm d} U/{\rm d}z \approx u_*^2$.

By combining (4.9) and (4.10), we find

(4.11)\begin{equation} \frac{\kappa^2 }{f } \frac{{\rm d} V}{{\rm d} z} = \left[ A-1 + \ln \left( \frac{z}{L_f} \right) \right] + c_\psi \frac{z}{L_s} \left[ \left(A-\frac{1}{2}\right) + \ln \left( \frac{z}{L_f} \right) + \frac{c_\psi}{2} \frac{z}{L_s} \right]. \end{equation}

The bottom boundary condition of (4.11) is $V(0)=0$. By integrating (4.11) from 0 to $z$ one can determine the mean spanwise velocity $V$ in the surface layer as

(4.12)\begin{equation} \frac{\kappa^2 V}{f} = z \left\{ \left[ A-2 + \ln \left( \frac{z}{L_f} \right) \right] + \frac{c_\psi }{2 } \frac{ z}{ L_s} \left[ A-1 + \ln \left( \frac{z}{L_f} \right) + \frac{c_\psi}{3} \frac{z}{L_s} \right] \right\}. \end{equation}

Above the surface layer, figure 7(b) shows that the spanwise velocity gradient $\textrm {d} V/\textrm {d} z$ scales as $u_*^2/(h^2 f)$ for NSBLs in the studied parameter range of $Ro$ and $\mu$. Thus, similar to the derivation of $A$, we can assume

(4.13)\begin{equation} \kappa \frac{\textrm{d} V}{\textrm{d} z} ={-}\frac{u_*^2}{ h^2 f} f_v(\xi), \end{equation}

where $f_v$ is independent of $Ro$ and $\mu$. Integrating (4.13) from below to the top of the boundary layer and noting that $L_f=u_*/f$, we obtain

(4.14)\begin{equation} \frac{\kappa}{u_*} (V_g - V) ={-} \frac{L_f}{h} \int_\xi^1 f_v (\xi') \,{\rm d} \xi'. \end{equation}

We further assume the mean spanwise velocity $V$ given by (4.12) and (4.14) matches at the interface $\xi =L_s/(c_2 h)$ or $z=L_s/c_2$, where $c_2$ is an empirical constant. Thus, by substituting (4.12) into (4.14) and noting that $\mu =L_f/L_s$, we obtain

(4.15)$$\begin{gather} \frac{\kappa V_g}{u_*} ={-} \frac{ L_f }{h} \int_{{L_s}/{c_2 h}}^1 f_{v} \,\textrm{d} \xi' + \frac{1}{\kappa c_2 \mu} \left[ A-2 - \ln \left( c_2 \mu \right) \right]\nonumber\\ + \frac{1}{\kappa c_2 \mu} \frac{c_\psi}{2 c_2} \left[ A-1 - \ln \left( c_2 \mu \right) + \frac{c_\psi}{3 c_2} \right]. \end{gather}$$

Substituting (1.1b), (3.3) and (4.5) into (4.15), we find

(4.16)\begin{equation} B = b_1 \sqrt{\mu}+ \frac{ b_2 }{ \sqrt{\mu} } + \frac{b_3}{\mu} , \end{equation}

where

(4.17)\begin{equation} \left. \begin{gathered} b_1 = \frac{1}{c_s} \int_{{1}/{c_2 c_s \sqrt{\mu}} }^1 f_{v} (\xi') \,\textrm{d} \xi', \quad b_2 = \frac{ a_1 }{\kappa c_2} \left( 1 + \frac{c_\psi}{2 c_2} \right), \\ b_3 = \frac{1}{\kappa c_2} \left[ \left( 1 + \ln c_2 - a_2 \right) \left( 1 + \frac{c_\psi}{2 c_2} \right) + 1 - \frac{c_\psi^2}{6 c_2^2} \right]. \end{gathered} \right\} \end{equation}

Finally, since the final term on the right-hand side of (4.16) is much smaller than the other two terms for $\mu \gg 10$, (4.16) can be written as

(4.18)\begin{equation} B = b_1 \sqrt{\mu}+ \frac{ b_2 }{ \sqrt{\mu} }. \end{equation}

Similar to $a_1$ we can assume $b_1$ to be constant. Therefore, the analytical expression of $B$ is given by (4.16), which also involves two empirical constants $(b_1, c_2)$ since $(a_1, c_1)$ are already given by (4.5). Note that (4.18) is different from equation (37) of Zilitinkevich & Esau (Reference Zilitinkevich and Esau2005). The reason is as follows: in their derivation Zilitinkevich & Esau (Reference Zilitinkevich and Esau2005) divided the surface layer into a logarithmic layer and a $z$-less stratification layer, and determined the spanwise velocity in the $z$-less stratification layer by only keeping the term proportional to $z^2$. This leads to a discontinuity of the spanwise velocity across the interface of the logarithmic layer and the $z$-less stratification layer. As a result, their obtained equation (37) only included a single term similar to $b_1 \sqrt {\mu }$, which underestimates the values of $B$ for small and moderate $\mu$ and thus a correction constant $b_0$ is introduced to their general expression of $B$ in (1.2b). Actually, (1.2b) reduces to $B = b c_s^3 \sqrt {\mu } + b_0 c_s / \sqrt {\mu }$ in the limit $\mu \gg 10$, which is the same as (4.18) when $b c_s^3 = b_1$ and $b_0 c_s=b_2$. Because we treat the surface layer as a whole, the continuity of the spanwise velocity is satisfied automatically. As a result, our derivation gives a higher-order term, and hence, there is no need to introduce an additional correction constant. The effectiveness of this treatment also reveals the essential physics of the GDL more clearly.

4.3. Further remarks

We remark that when deriving the analytical expressions of $A$ and $B$, i.e. (4.5) and (4.16) we have made two ansatzes: first, the streamwise and spanwise velocity profile between the inner and outer layers matches at a dimensionless height $1/(c_1c_s\sqrt {\mu })$ and $1/(c_2c_s\sqrt {\mu })$, respectively; and second, $a_1$ and $b_1$ can be approximated as constants that are independent of $\mu$. These ansatzes enable us to predict well $(A, B)$ without knowing the forms of $(f_u, f_v)$ and with reasonable values of $(c_1, c_2)$, see § 5 below. This fact also verifies the ansatzes themselves. However, if one assumes that the streamwise and spanwise velocity profiles between the inner and outer layers match at dimensionless heights independent of $\mu$, then one would find these heights are too close to the top of boundary layer, and therefore, this assumption should be abandoned. On the other hand, if one assumes that $a_1$ and $b_1$ can be expanded in power series of $1/\sqrt {\mu }$, then the values of $(c_1, c_2)$ cannot be determined without knowing the forms of $(f_u, f_v)$ and thus this assumption should also be discarded. We further remark that the two velocities $U_g$ and $V_g$ can be obtained by substituting (4.5) and (4.16) into (1.1a,b). At first sight, these two velocities seem to be expressed to different orders of $\mu$. However, we emphasize that the high-order term in $B$ or $V_g$ is obtained directly by ensuring the continuity of the spanwise velocity at the matching height, which differs from the formal asymptotic expansion. We also emphasize that, since (4.5) and (4.16) predict $(A,B)$ accurately, they also predict well $(U_g, V_g)$ and hence the geostrophic drag coefficient and the cross-isobaric angle.