Journal topic
Wind Energ. Sci., 5, 355–374, 2020
https://doi.org/10.5194/wes-5-355-2020
Wind Energ. Sci., 5, 355–374, 2020
https://doi.org/10.5194/wes-5-355-2020

Research article 26 Mar 2020

Research article | 26 Mar 2020

# Rossby number similarity of an atmospheric RANS model using limited-length-scale turbulence closures extended to unstable stratification

Rossby number similarity of an atmospheric RANS model using limited-length-scale turbulence closures extended to unstable stratification
Maarten Paul van der Laan, Mark Kelly, Rogier Floors, and Alfredo Peña Maarten Paul van der Laan et al.
• DTU Wind Energy, Technical University of Denmark, Risø Campus, Frederiksborgvej 399, 4000 Roskilde, Denmark

Correspondence: Maarten Paul van der Laan (plaa@dtu.dk)

Abstract

The design of wind turbines and wind farms can be improved by increasing the accuracy of the inflow models representing the atmospheric boundary layer. In this work we employ one-dimensional Reynolds-averaged Navier–Stokes (RANS) simulations of the idealized atmospheric boundary layer (ABL), using turbulence closures with a length-scale limiter. These models can represent the mean effects of surface roughness, Coriolis force, limited ABL depth, and neutral and stable atmospheric conditions using four input parameters: the roughness length, the Coriolis parameter, a maximum turbulence length, and the geostrophic wind speed. We find a new model-based Rossby similarity, which reduces the four input parameters to two Rossby numbers with different length scales. In addition, we extend the limited-length-scale turbulence models to treat the mean effect of unstable stratification in steady-state simulations. The original and extended turbulence models are compared with historical measurements of meteorological quantities and profiles of the atmospheric boundary layer for different atmospheric stabilities.

1 Introduction

Wind turbines operate in the turbulent atmospheric boundary layer (ABL) but are designed with simplified inflow conditions that represent analytic wind profiles of the atmospheric surface layer (ASL). The ASL corresponds to roughly the first 10 % of the ABL, typically less than 100 m, while the tip heights of modern wind turbines are now sometimes beyond 200 m. Hence, there is a need for inflow models that represent the entire ABL in order to improve the design of wind turbines and wind farms. Such a model should be simple enough to efficiently improve the chain of design tools used by the wind energy industry.

The ABL is complex and changes continuously over time. Idealized, steady-state models can represent long-term-averaged velocity and turbulence profiles of the real ABL, including the effects of Coriolis, atmospheric stability, capping inversion, homogeneous surface roughness and flat terrain; here we exclude the effects of flow inhomogeneity and nonstationarity, which are typically considered by mesoscale and three-dimensional time-varying models. In this work, we investigate idealized ABL models that are based on one-dimensional Reynolds-averaged Navier–Stokes (RANS) equations, where the only spatial dimension is the height above ground. The output of the model can be used as inflow conditions for three-dimensional RANS simulations of complex terrain and wind farms . Turbulence is modeled here by two limited-length-scale turbulence closures, the mixing-length model of and the two-equation kε model of . These turbulence models can simulate one-dimensional stable and neutral ABLs without the necessity of a temperature equation and a momentum source term of buoyancy. In other words, all temperature effects are represented by the turbulence model. The limited-length-scale turbulence models depend on four parameters: the roughness length, the Coriolis parameter, the geostrophic wind speed, and a chosen maximum turbulence length scale that is related to the ABL depth. We show that the normalized profiles of wind speed, wind direction and turbulence quantities are only dependent on two dimensionless parameters that represent the ratio of the inertial to the Coriolis force, based on two different length scales: the roughness length and the maximum turbulence length scale. These dimensionless parameters are Rossby numbers. The Rossby number based on the roughness length is known as the surface Rossby number, as introduced by , while the Rossby number based on the maximum turbulence length is a new dimensionless parameter. The obtained model-based Rossby number similarity is used to validate a range of simulations with historical measurements of the geostrophic drag coefficient and cross-isobar angle. In addition, we show that both RANS models' solutions are bounded by two analytic solutions of the idealized ABL.

The limited-length-scale turbulence closures of and can model the effect of stable and neutral stability but cannot model the unstable atmosphere. We propose simple extensions to solve this issue and validate the results of the extended kε model with measurements of wind speed and wind direction profiles. The model extensions lead to a third Rossby number, where the length scale is based on the Obukhov length. The limited mixing-length model is not considered in the comparison with measurements because we are mainly interested in the limited-length-scale kε model. The kε model is more applicable to wind energy applications because it can also provide an estimate of the turbulence intensity, which is not available from the limited mixing-length model of . The limited mixing-length model is applied here to show that the same model-based Rossby number similarity as obtained for the kε model is recovered.

The article is structured as follows. The background and theory of the idealized ABL are discussed in Sect. 2. Extensions to unstable surface layer stratification are presented in Sect. 3. Section 4 presents the methodology of the one-dimensional RANS simulations. The model-based Rossby similarity is illustrated in Sect. 5. The simulation results of the limited-length-scale kε model including the extension to unstable conditions are compared with measurements in Sect. 6.

2 Background and theory – idealized models of the ABL

We model the mean steady-state flow in an idealized ABL. Here idealized refers to flow over homogeneous and flat terrain under barotropic conditions such that the geostrophic wind does not vary with height. This flow can be described by the incompressible RANS equations for momentum, where the contribution from the molecular viscosity is neglected due to the high Reynolds number:

$\begin{array}{}\text{(1)}& \begin{array}{rl}& \frac{\mathrm{D}U}{\mathrm{D}t}={f}_{\mathrm{c}}\left(V-{V}_{\mathrm{G}}\right)+\frac{\mathrm{d}}{\mathrm{d}z}\left({\mathit{\nu }}_{T}\frac{\mathrm{d}U}{\mathrm{d}z}\right)=\mathrm{0},\\ & \frac{\mathrm{D}V}{\mathrm{D}t}=-{f}_{\mathrm{c}}\left(U-{U}_{\mathrm{G}}\right)+\frac{\mathrm{d}}{\mathrm{d}z}\left({\mathit{\nu }}_{T}\frac{\mathrm{d}V}{\mathrm{d}z}\right)=\mathrm{0},\end{array}\end{array}$

where U and V are the mean horizontal velocity components, UG and VG are the corresponding mean geostrophic velocities, fc=2Ωsin (λ) is the Coriolis parameter with Ω as Earth's angular velocity and λ as the latitude, and z is the height above ground. In addition, the Reynolds stresses $\stackrel{\mathrm{‾}}{{u}^{\prime }{w}^{\prime }}$ and $\stackrel{\mathrm{‾}}{{v}^{\prime }{w}^{\prime }}$ are modeled by the linear stress–strain relationship of : $\stackrel{\mathrm{‾}}{{u}^{\prime }{w}^{\prime }}=-{\mathit{\nu }}_{T}\mathrm{d}U/\mathrm{d}z$ and $\stackrel{\mathrm{‾}}{{v}^{\prime }{w}^{\prime }}=-{\mathit{\nu }}_{T}\mathrm{d}V/\mathrm{d}z$, where νT is the eddy viscosity. The boundary conditions for U and V are $U=V=\mathrm{0}$ at z=z0 and U=UG and V=VG at z→∞, where z0 is the roughness length. Note that it is possible to write the two momentum equations as a single ordinary differential equation:

$\begin{array}{}\text{(2)}& \frac{\mathrm{d}}{\mathrm{d}z}\left({\mathit{\nu }}_{T}\frac{\mathrm{d}W}{\mathrm{d}z}\right)=i{f}_{\mathrm{c}}W,\end{array}$

where $W\equiv \left(U-{U}_{\mathrm{G}}\right)+i\left(V-{V}_{\mathrm{G}}\right)$ is a complex variable and ${i}^{\mathrm{2}}=-\mathrm{1}$.

The eddy viscosity, νT, needs to be modeled in order to close the system of equations. The eddy viscosity can be written as ${\mathit{\nu }}_{T}={u}_{*}\mathrm{\ell }$, where u* and  represent turbulence velocity and turbulence length scales. For a constant eddy viscosity, the equations can be solved analytically, and the solution is known as the Ekman spiral (Ekman1905), which includes the wind direction change with height due to Coriolis effects. The Ekman spiral can also be considered a laminar solution, since one can neglect the turbulence in the momentum equations and set the molecular viscosity to determine the rate of mixing. For an eddy viscosity that increases linearly with height, the equations can also be solved analytically, as introduced by and discussed by and . The two analytic solutions are provided in Appendix A. One can relate the analytic solution of to the (neutral) ASL (zzi), while the Ekman spiral is more valid for altitudes around the ABL depth zi. Neither of the two analytic solutions result in a realistic representation of the entire (idealized) ABL. A combination of both a linear eddy viscosity for zzi and a constant eddy viscosity for zzi should provide a more realistic solution. For example, the eddy viscosity could have the form ${\mathit{\nu }}_{T}=\mathit{\kappa }{u}_{*\mathrm{0}}z\mathrm{exp}\left(-z/h\right)$, where νT increases linearly with height for zh as expected in the surface layer; then it reaches a maximum value at z=h, and decreases to zero for z>h. Note that u*0 is the friction velocity near the surface. derived a number of solutions for a variable eddy viscosity, although an explicit solution for the entire idealized ABL with a realistic eddy viscosity (in the previously mentioned form) has not been found yet. Hence, numerical methods are still necessary, and one of the simplest numerical models for the idealized ABL is given by using Prandtl's mixing-length model:

$\begin{array}{}\text{(3)}& {\mathit{\nu }}_{T}={\mathrm{\ell }}^{\mathrm{2}}\mathcal{S},\end{array}$

where $\mathcal{S}=\sqrt{\left(\mathrm{d}U/\mathrm{d}z{\right)}^{\mathrm{2}}+\left(\mathrm{d}V/\mathrm{d}z{\right)}^{\mathrm{2}}}=|\mathrm{d}W/\mathrm{d}z|$ is the magnitude of the strain-rate tensor, and   is prescribed as a turbulence length scale,

$\begin{array}{}\text{(4)}& \mathrm{\ell }=\frac{\mathit{\kappa }z}{\mathrm{1}+\frac{\mathit{\kappa }z}{{\mathrm{\ell }}_{\mathrm{max}}}},\end{array}$

where κz is the turbulence length scale in the neutral surface layer with κ as the von Kármán constant, and max is a maximum turbulence length scale. It is also possible to model the eddy viscosity with a two-equation turbulence closure, e.g., the kε model:

$\begin{array}{}\text{(5)}& {\mathit{\nu }}_{T}={C}_{\mathit{\mu }}\frac{{k}^{\mathrm{2}}}{\mathit{\epsilon }},\end{array}$

with Cμ as a model parameter, k as the turbulent kinetic energy, and ε as the dissipation of k. Both k  and ε are modeled by a transport equation:

$\begin{array}{}\text{(6)}& & \frac{\mathrm{D}k}{\mathrm{D}t}=\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{{\mathit{\nu }}_{T}}{{\mathit{\sigma }}_{k}}\frac{\mathrm{d}k}{\mathrm{d}z}\right)+\mathcal{P}-\mathit{\epsilon },\text{(7)}& & \frac{\mathrm{D}\mathit{\epsilon }}{\mathrm{D}t}=\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{{\mathit{\nu }}_{T}}{{\mathit{\sigma }}_{\mathit{\epsilon }}}\frac{\mathrm{d}\mathit{\epsilon }}{\mathrm{d}z}\right)+\left({C}_{\mathit{\epsilon },\mathrm{1}}\mathcal{P}-{C}_{\mathit{\epsilon },\mathrm{2}}\mathit{\epsilon }\right)\frac{\mathit{\epsilon }}{k},\end{array}$

where 𝒫 is the turbulence production, and σk, σε, Cε,1, and Cε,2 are model constants that should follow the relationship ${\mathit{\kappa }}^{\mathrm{2}}={\mathit{\sigma }}_{\mathit{\epsilon }}\sqrt{{C}_{\mathit{\mu }}}\left({C}_{\mathit{\epsilon },\mathrm{2}}-{C}_{\mathit{\epsilon },\mathrm{1}}\right)$. When using the standard kε model calibrated for atmospheric flows , the turbulence length scale or eddy viscosity will keep increasing until a boundary layer depth is formed and the analytic solution of is approximated. proposed modifying the transport equation of ε, such that a maximum turbulence length scale is enforced by replacing the constant Cε,1 with a variable parameter ${C}_{\mathit{\epsilon },\mathrm{1}}^{*}$:

$\begin{array}{}\text{(8)}& {C}_{\mathit{\epsilon },\mathrm{1}}^{*}={C}_{\mathit{\epsilon },\mathrm{1}}+\left({C}_{\mathit{\epsilon },\mathrm{2}}-{C}_{\mathit{\epsilon },\mathrm{1}}\right)\frac{\mathrm{\ell }}{{\mathrm{\ell }}_{\mathrm{max}}},\end{array}$

where the turbulence length scale is modeled as $\mathrm{\ell }={C}_{\mathit{\mu }}^{\mathrm{3}/\mathrm{4}}{k}^{\mathrm{3}/\mathrm{2}}/\mathit{\epsilon }$. This limited-length-scale kε model behaves very similarly to the mixing-length model of (Eqs. 3 and 4). For ℓ≪ℓmax, the surface layer solution is obtained, while for ℓ∼ℓmax, the source terms in the transport equation of ε cancel each other out (${C}_{\mathit{\epsilon },\mathrm{1}}^{*}\mathcal{P}\sim {C}_{\mathit{\epsilon },\mathrm{1}}\mathit{\epsilon }$), and the turbulence length scale is limited. For a given z0, G, and fc, the ABL depth can be controlled by max. This means that max is related to zi; noted that ${\mathrm{\ell }}_{\mathrm{max}}\sim {z}_{i}/\mathrm{3}$ applies to typical neutral ABLs. However, the simulated boundary layer depth using the kε model of has an approximate dependence of ${z}_{i}\propto \left(G/|{f}_{\mathrm{c}}|{\right)}^{\mathrm{1}-a}{\mathrm{\ell }}_{\mathrm{max}}^{a}$ with a≈0.6, which we will further discuss in Sect. 5. A summary of the discussed eddy viscosity closures is listed in Table 1. Figure 1 compares the analytic solutions of and with the numerical solutions of the limited mixing-length model of and the limited-length-scale kε of in terms of wind speed, wind direction, $\mathit{\theta }=\mathrm{arctan}\left(V/U\right)$, and eddy viscosity. The Ekman spiral is depicted with two constant eddy viscosities, which only translates the solution vertically. In addition, we have chosen ${f}_{\mathrm{c}}={\mathrm{10}}^{-\mathrm{4}}$ s−1, G=10 m s−1, and ${z}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{2}}$ m. The numerical solutions are shown for a range of max values. It is clear that the ABL depth decreases for lower values of max, for both numerical models, and their solutions behave similarly. A lower max also results in a higher shear and wind veer and a lower eddy viscosity, which are characteristics of a stable ABL. Hence, the limited-length-scale turbulence closures can model the effects of stable stratification by solely limiting the turbulence length scale, without the need of a temperature equation or buoyancy source terms. When max→0 m (note that there is minimal limit of max in order to obtain numerically stable results), the solution approaches the Ekman spiral because the eddy viscosity in the ABL can be approximated by a constant eddy viscosity. Hence, the maximum change in wind direction simulated by the kε model of is that of the Ekman spiral: 45. For large max values, the numerical solution approximates the analytic solution of but does not match it because their eddy viscosities are different for zzi.

Figure 1Comparison of analytic and numerical solutions of existing idealized ABL models using ${f}_{\mathrm{c}}={\mathrm{10}}^{-\mathrm{4}}$ s−1, G=10 m s−1, and ${z}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{2}}$ m for different model parameters. (a) Wind speed. (b) Wind direction. (c) Eddy viscosity.

Table 1Eddy viscosity closures for the idealized ABL.

3 Extension to unstable surface layer stratification

The two limited-length-scale turbulence closures discussed in Sect. 2 can be used to model neutral and stable ABLs without the need of a temperature equation and buoyancy forces. However, it is not possible to model the unstable ABL because the turbulence length scale is only limited, not enhanced, i.e., ℓ≤κz. In order to model unstable conditions, we need to extend the models such that the turbulence length scale is enhanced in the surface layer, ℓ>κz, which we present in the following sections for each turbulence closure.

## 3.1 Limited mixing-length model

One can generically parameterize the turbulence length scale  as a “parallel” combination of ASL and ABL scales,

$\begin{array}{}\text{(9)}& \frac{\mathrm{1}}{\mathrm{\ell }}=\frac{\mathrm{1}}{{\mathrm{\ell }}_{\mathrm{ASL}}}+\frac{\mathrm{1}}{{\mathrm{\ell }}_{\mathrm{ABL}}}.\end{array}$

chose ASL=κz and ABL=ℓmax to arrive at Eq. (4). If we choose to set

$\begin{array}{}\text{(10)}& {\mathrm{\ell }}_{\mathrm{ASL}}=\frac{\mathit{\kappa }z}{{\mathit{\varphi }}_{\mathrm{m}}}\end{array}$

following the turbulence length scale that is a result of Monin–Obukhov Similarity Theory (MOST, ) – where

$\begin{array}{}\text{(11)}& {\mathit{\varphi }}_{\mathrm{m}}={\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}z/L\right)}^{-\mathrm{1}/\mathrm{4}}\end{array}$

is the dimensionless velocity gradient for unstable conditions, with γ1≈16 as shown by Dyer (1974), and L is the Obukhov length – then it is possible to extend the limited mixing-length model of to unstable surface layer stratification, as

$\begin{array}{}\text{(12)}& \mathrm{\ell }=\frac{\mathit{\kappa }z}{{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}z/L\right)}^{-\mathrm{1}/\mathrm{4}}+\mathit{\kappa }z/{\mathrm{\ell }}_{\mathrm{max}}}.\end{array}$

Approaching neutral conditions, ${L}^{-\mathrm{1}}\to \mathrm{0}$, the original length-scale model of is obtained. Note that in stable conditions, ${\mathit{\varphi }}_{\mathrm{m}}=\mathrm{1}+\mathit{\beta }z/L$, so the resulting turbulence length can also be rewritten in the form of Eqs. (4) and (9), using an effective maximum turbulence length scale of

$\begin{array}{}\text{(13)}& {\mathrm{\ell }}_{\mathrm{ABL},\mathrm{stable}}^{-\mathrm{1}}={\mathrm{\ell }}_{\mathrm{max},\mathrm{eff}}^{-\mathrm{1}}\equiv {\mathrm{\ell }}_{\mathrm{max}}^{-\mathrm{1}}+\mathit{\beta }/\left(\mathit{\kappa }L\right).\end{array}$

Thus we can simply use the original length-scale model of for stable and neutral conditions; the stable ϕm function simply informs the selection of max,eff, following Eq. (13).

## 3.2 Limited-length-scale k–ε model

argued that, in stable conditions, the limited-length-scale kε model of overpredicts  in the surface layer compared to MOST, where ${\mathrm{\ell }}_{\mathrm{max}}=L\mathit{\kappa }/\mathit{\beta }$ and β≈5. They proposed a more complicated expression for ${C}_{\mathit{\epsilon },\mathrm{1}}^{*}$ in the transport equation of ε compared to the original model of . alternatively prescribed a coefficient in the buoyant term of the ε equation, depending on ℓ∕ℓmax and being similar to the production-related term that gives results consistent (at least asymptotically) with MOST. We find that the correction of provides a better match of the turbulence length scale within the surface layer compared to MOST with respect to the original kε model of . However, we also find that a larger overshoot of the turbulence length scale around the ABL depth is found when Coriolis is included. Alternatively, one could improve the surface layer solution of the original model of by simply reducing max by roughly 20 %. Therefore, we choose to use the model of as our starting point.

In order to account for the increase in turbulence length scale in the surface layer under unstable conditions, we add a buoyancy source term  in the kε transport equations:

$\begin{array}{}\text{(14)}& & \frac{\mathrm{D}k}{\mathrm{D}t}=\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{{\mathit{\nu }}_{T}}{{\mathit{\sigma }}_{k}}\frac{\mathrm{d}k}{\mathrm{d}z}\right)+\mathcal{P}-\mathit{\epsilon }+\mathcal{B},\text{(15)}& & \frac{\mathrm{D}\mathit{\epsilon }}{\mathrm{D}t}=\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{{\mathit{\nu }}_{T}}{{\mathit{\sigma }}_{\mathit{\epsilon }}}\frac{\mathrm{d}\mathit{\epsilon }}{\mathrm{d}z}\right)+\left({C}_{\mathit{\epsilon },\mathrm{1}}^{*}\mathcal{P}-{C}_{\mathit{\epsilon },\mathrm{2}}\mathit{\epsilon }+{C}_{\mathit{\epsilon },\mathrm{3}}^{*}\mathcal{B}\right)\frac{\mathit{\epsilon }}{k}.\end{array}$

Here  is modeled as

$\begin{array}{}\text{(16)}& \mathcal{B}=-{\mathit{\nu }}_{T}\left[{\left(\frac{\mathrm{d}U}{\mathrm{d}z}\right)}^{\mathrm{2}}+{\left(\frac{\mathrm{d}V}{\mathrm{d}z}\right)}^{\mathrm{2}}\right]\frac{z}{L},=-{\mathit{\nu }}_{T}{\mathcal{S}}^{\mathrm{2}}\frac{z}{L}\end{array}$

following MOST, using the similarity functions of Dyer (1974) as discussed in . We use the flow-dependent parameter ${C}_{\mathit{\epsilon },\mathrm{3}}^{*}\equiv \mathrm{1}+{\mathit{\alpha }}_{B}\left({C}_{\mathit{\epsilon },\mathrm{1}}-{C}_{\mathit{\epsilon },\mathrm{2}}\right)$ of , which for unstable conditions includes the prescription

$\begin{array}{}\text{(17)}& {\mathit{\alpha }}_{B}=\mathrm{1}-\left[\mathrm{1}+\frac{\left({C}_{\mathit{\epsilon },\mathrm{2}}-\mathrm{1}\right)}{\left({C}_{\mathit{\epsilon },\mathrm{2}}-{C}_{\mathit{\epsilon },\mathrm{1}}\right)}\right]\frac{\mathrm{\ell }}{{\mathrm{\ell }}_{\mathrm{max}}},\end{array}$

amenable to the free-convection limit: $\mathit{\epsilon }/\mathcal{B}\to \mathrm{1}$ for $\mathcal{P}/\mathcal{B}\to \mathrm{0}$. Further, αB→1 as ℓ→0, matching neutral conditions since zL also vanishes then. The prescription (Eq. 17) results in

$\begin{array}{}\text{(18)}& {C}_{\mathit{\epsilon },\mathrm{3}}^{*}=\mathrm{1}+{C}_{\mathit{\epsilon },\mathrm{1}}-{C}_{\mathit{\epsilon },\mathrm{2}}+\left(\mathrm{2}{C}_{\mathit{\epsilon },\mathrm{2}}-{C}_{\mathit{\epsilon },\mathrm{1}}-\mathrm{1}\right)\frac{\mathrm{\ell }}{{\mathrm{\ell }}_{\mathrm{max}}},\end{array}$

which also means that $C{*}_{\mathit{\epsilon },\mathrm{3}}\to {C}_{\mathit{\epsilon },\mathrm{2}}$ approaching the effective ABL top (ℓ→ℓmax), so sources and sinks of ε balance in Eq. (15), i.e., $\mathcal{P}-\mathit{\epsilon }+\mathcal{B}$, all have the same coefficient Cε,2.

4 Methodology of numerical simulations

The one-dimensional numerical simulations are performed with EllipSys1D , which is a simplified one-dimensional version of EllipSys3D, initially developed by and . EllipSys1D is a finite volume solver for incompressible flow, with collocated storage of flow variables. It is assumed that the vertical velocity is zero and the pressure gradients are constant, which is valid in an idealized ABL, as discussed in Sect. 2. As a consequence, it is not necessary to solve the pressure correction equation that is normally used to ensure mass conservation.

## 4.1 Ambient turbulence in the limited-length-scale k–ε turbulence model

The limited-length-scale kε model typically simulates an eddy viscosity that decays to zero for z→∞, which can lead to numerical instability. While, for example,  chose to set upper limits for k and ε to prevent numerical instabilities, we prefer a more physical method, including ambient source terms Sk,amb and Sε,amb in the k and ε transport equations, respectively. Following , we set

$\begin{array}{}\text{(19)}& {S}_{k,\mathrm{amb}}={\mathit{\epsilon }}_{\mathrm{amb}},\phantom{\rule{0.25em}{0ex}}{S}_{\mathit{\epsilon },\mathrm{amb}}={C}_{\mathit{\epsilon },\mathrm{2}}\frac{{\mathit{\epsilon }}_{\mathrm{amb}}^{\mathrm{2}}}{{k}_{\mathrm{amb}}}.\end{array}$

When all sources of turbulence are zero ($\mathcal{P}=B=\mathrm{0}$) and the diffusion terms are zero ($\mathrm{d}k/\mathrm{d}z=\mathrm{d}\mathit{\epsilon }/\mathrm{d}z=\mathrm{0}$), then k=kamb and ε=εamb. To be consistent with the equations solved, we define the ambient turbulence quantities in terms of the driving parameters, G and max:

$\begin{array}{}\text{(20)}& \begin{array}{rl}& {\mathrm{\ell }}_{\mathrm{amb}}={C}_{\mathrm{amb}}{\mathrm{\ell }}_{\mathrm{max}},\phantom{\rule{0.25em}{0ex}}{k}_{\mathrm{amb}}=\frac{\mathrm{3}}{\mathrm{2}}{I}_{\mathrm{amb}}^{\mathrm{2}}{G}^{\mathrm{2}},\\ & {\mathit{\epsilon }}_{\mathrm{amb}}={C}_{\mathit{\mu }}^{\mathrm{3}/\mathrm{4}}\frac{{k}_{\mathrm{amb}}^{\frac{\mathrm{3}}{\mathrm{2}}}}{{\mathrm{\ell }}_{\mathrm{amb}}}={C}_{\mathit{\mu }}^{\mathrm{3}/\mathrm{4}}\frac{\mathrm{3}}{\mathrm{2}}\sqrt{\frac{\mathrm{3}}{\mathrm{2}}}\frac{{I}_{\mathrm{amb}}^{\mathrm{3}}}{{C}_{\mathrm{amb}}}\frac{{G}^{\mathrm{3}}}{{\mathrm{\ell }}_{\mathrm{max}}}.\end{array}\end{array}$

Here Iamb is the total turbulence intensity1 above the (simulated) ABL, and Camb is the ratio of the turbulence length scale above the ABL (amb) to maximum turbulence length scale (max). We choose small values for ${I}_{\mathrm{amb}}={\mathrm{10}}^{-\mathrm{6}}$ and ${C}_{\mathrm{amb}}={\mathrm{10}}^{-\mathrm{6}}$, such that the ambient turbulence does not affect the solution for U and V, while the numerical stability is maintained. It should be noted that the overshoot in ℓ∕ℓmax that can occur near the ABL depth is still affected by the ambient values. and chose to use a limiter on ε to avoid an overshoot in , but we choose not to use it. In general, we prefer to avoid limiters because they can break the Rossby number similarity that is presented in Sect. 5.

## 4.2 Numerical setup

The flow is driven by a constant pressure gradient using a prescribed constant geostrophic wind speed. The initial wind speed is set to the geostrophic wind speed at all heights. During the solving procedure, the ABL depth grows from the ground until convergence is achieved, which occurs when the growth rate of the ABL depth is negligible because a balance between the prescribed pressure gradient, the Coriolis forces, and the turbulence stresses is obtained. The flow that we are solving is relatively stiff, and we choose to include the transient terms using a second-order three-level implicit method with a large time step that is set as $\mathrm{1}/|{f}_{\mathrm{c}}|$ s. All spatial gradients are discretized by a second-order central-difference scheme. Convergence is typically achieved after 105 iterations, which takes about 10 s on a single 2.7 GHz CPU. The domain height is set to 105 m to ensure that the ABL depth is significantly smaller than the domain height for all flow cases considered. The numerical grid represents a line, where the first cell height is set to 10−2 m. The cells are stretched for increasing heights using an expansion ratio of about 1.2. The grid consists of 384 cells, which is based on a grid refinement study presented in Sect. 4.3. A rough-wall boundary condition is set at the ground, as discussed by . For the length-scale-limited kε model, this means that we set ε at the first cell, use a Neumann condition for k, and the shear stress at the wall is defined by the neutral surface layer. The first cell is placed on top of the roughness length, which allows us to choose the first cell height independent of the roughness length. This means that we add the roughness length to all relations that include z, i.e., z+z0. For the limited mixing-length model, we simply set the eddy viscosity from the neutral surface layer at the first cell. Neumann conditions are set for all flow variables at the top boundary.

The turbulence model constants of the kε model are set as (Cμ, σk, σε, Cε,1, Cε,2, κ)=(0.03, 1.0, 1.3, 1.21, 1.92, 0.4). The chosen Cμ value is based on neutral ASL measurements, as discussed by , and Cε,1 is used to maintain the neutral ASL solution of the kε model.

## 4.3 Grid refinement study

A grid refinement study of the numerical setup is performed for the limited-length-scale kε model of , using 48, 96, 192, 384, and 768 cells. We choose ${f}_{\mathrm{c}}={\mathrm{10}}^{-\mathrm{4}}$ s−1, ${z}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$ m, and G=10 m s−1 for max=100 and max=1 m. The results in terms of the wind speed of each grid are depicted in Fig. 2 for both values of max. For max=100 m, the largest difference with respect to the finest grid is 0.5 %, 0.2 %, 0.09 %, and 0.03 % for 48, 96, 192, and 384 cells, respectively, located at the first cell near the wall boundary. When using max=1 m, a small ABL depth of 100 m is simulated with a sharp low-level jet. In the enlarged plot of Fig. 2b, one can see how the grid size affects the low-level jet, where the largest difference with respect to the finest grid is 1 %, 0.2 %, 0.04 %, and 0.01 %, for 48, 96, 192, and 384 cells, respectively. We find similar results for the limited mixing-length model of . In addition, the turbulence model extensions to unstable surface layer stratification typically show smaller differences between the grids due to the enhanced mixing and the use of a high max value that represents a convective ABL. Hence, our choice of using 384 cells is conservative.

Figure 2Grid refinement study of the one-dimensional RANS simulation using the limited-length-scale kε model, for 48, 96, 192, 384, and 768 cells. (a) max=100 m. (b) max=1 m.

5 Rossby number similarity in numerical and analytical solutions

The numerical solution of the original limited-length-scale turbulence closures of and depend on four parameters: fc (s−1), G (m s−1), max (m), and z0 (m). Applying the Buckingham π theorem, it is clear that there should exist two dimensionless numbers that define the entire solution, since the four dimensional parameters only have two dimensions (meters and seconds). This can be shown by writing a nondimensional momentum equation in complex form (Eq. 2) using the nondimensional variables ${W}^{\prime }\equiv W/\mathcal{U}$, ${\mathit{\nu }}_{{T}_{\mathrm{0}}}\equiv {\mathit{\nu }}_{T}/\left(\mathcal{U}\mathcal{L}\right)$ and ${z}_{\mathrm{0}}\equiv z/\mathcal{L}$, where 𝒰 and  are characteristic velocity and length scales, respectively:

$\begin{array}{}\text{(21)}& Ro\frac{\mathrm{d}}{\mathrm{d}{z}^{\prime }}\left({\mathit{\nu }}_{T}^{\prime }\frac{\mathrm{d}{W}^{\prime }}{\mathrm{d}{z}^{\prime }}\right)=i{W}^{\prime }.\end{array}$

Here, Ro is the Rossby number, $Ro=\mathcal{U}/\left(|{f}_{\mathrm{c}}|\mathcal{L}\right)$, which describes the ratio of the inertial (advective) tendency to the Coriolis force. If we apply the original mixing-length model of for ${\mathit{\nu }}_{T}^{\prime }$ using Eqs. (3) and (4), then Eq. (21) can be written as

$\begin{array}{}\text{(22)}& Ro\frac{\mathrm{d}}{\mathrm{d}{z}^{\prime }}\left({\left[\frac{\mathit{\kappa }{z}^{\prime }}{\mathrm{1}+\mathit{\kappa }{z}^{\prime }\mathcal{L}/{\mathrm{\ell }}_{\mathrm{max}}}\right]}^{\mathrm{2}}\left|\frac{\mathrm{d}{W}^{\prime }}{\mathrm{d}{z}^{\prime }}\right|\frac{\mathrm{d}{W}^{\prime }}{\mathrm{d}{z}^{\prime }}\right)=i{W}^{\prime },\end{array}$

where $-\mathcal{L}/{\mathrm{\ell }}_{\mathrm{max}}$, is a second dimensionless number. If we choose 𝒰=G and ℒ=z0, we may define two Rossby-like numbers, with characteristic length scales based on z0 and max, respectively:

$\begin{array}{}\text{(23)}& R{o}_{\mathrm{0}}\equiv \frac{G}{|{f}_{\mathrm{c}}|{z}_{\mathrm{0}}},\phantom{\rule{0.25em}{0ex}}R{o}_{\mathrm{\ell }}\equiv \frac{G}{|{f}_{\mathrm{c}}|{\mathrm{\ell }}_{\mathrm{max}}}.\end{array}$

Here, we have obtained Ro by rewriting the second dimensionless number ℒ∕ℓmax as the ratio of the two Rossby numbers: $\mathcal{L}/{\mathrm{\ell }}_{\mathrm{max}}=\left[\mathcal{U}/\left(|{f}_{\mathrm{c}}|{\mathrm{\ell }}_{\mathrm{max}}\right)\right]/\left[\mathcal{U}/\left(|{f}_{\mathrm{c}}|\mathcal{L}\right)\right]$. Ro0 is known as the surface Rossby number, first introduced by ; it also resembles a ratio of (inertial) boundary layer depth to z0. Analogously, Ro is like the ratio of two boundary layer depths, ${f}_{\mathrm{c}}/{u}_{*\mathrm{0}}$ and zi (e.g., Arya and Wyngaard1975); here max is a proxy for zi, acting as a “lid” for the ABL. Considering Eq. (23), we have reduced the number of dependent parameters from four to two: f(fc, G, max, z0)→f(Ro0, Ro). For a fixed surface roughness z0, the ratio of the two Rossby numbers is then the only dependent parameter:

$\begin{array}{}\text{(24)}& {\mathrm{\ell }}_{\mathrm{max}}=\frac{R{o}_{\mathrm{0}}}{R{o}_{\mathrm{\ell }}}{z}_{\mathrm{0}};\end{array}$

i.e., the ratio of simulated ABL depth to z0 is the lone parameter. found a characteristic maximum ABL turbulence length scale of 0.00027 $G/|{f}_{\mathrm{c}}|$ for the Leipzig wind profile (Lettau1950), which equating with max corresponds to Ro≃3700.

Figure 3 depicts the Rossby number similarity of our one-dimensional RANS simulations using the original limited-length-scale turbulence closures of , Fig. 3a–c, and , Fig. 3d–g. Four combinations of Ro0 (106 and 109) and Ro (103 and 105) are used, each simulated with four combinations of G (10 and 20 m s−1) and fc ($\mathrm{5}×{\mathrm{10}}^{-\mathrm{5}}$ and 10−4 s−1). The roughness length and maximum turbulence length scale follow from Eq. (23) and cover a wide range of z0 from 10−4 to 0.4 m and max of 100–400 m. Figure 3 shows that normalized wind speed, wind direction, and turbulence quantities for both turbulence closures are only dependent on Ro0 and Ro. Both turbulence closures produce similar results in terms of wind speed, wind direction, and eddy viscosity. The limited-length-scale kε model of also predicts a total turbulence intensity I (Fig. 3g) and a turbulence length scale  (not shown in Fig. 3), which are only dependent on the two Rossby numbers. In addition, the total turbulence intensity close to the surface only depends on Ro0, while further away, it is mainly influenced by Ro with a weaker dependence on Ro0.

Figure 3Rossby number similarity of the original turbulence closures. (a–c) Limited mixing-length model. (d–g) Limited-length-scale kε model.

Considering the non-neutral ABL with Coriolis effects but ignoring the strength of capping inversion (entrainment), in the micrometeorological literature the parameter ${u}_{*\mathrm{0}}/\left(|{f}_{\mathrm{c}}|L\right)$ is typically invoked (e.g., Arya1975; Zilitinkevich1989). This can also be considered as a third Rossby number, which in our context of using G instead of u*0 is

$\begin{array}{}\text{(25)}& R{o}_{{L}_{-}}\equiv \frac{-G}{|{f}_{\mathrm{c}}|L};\end{array}$

here the subscript (${}_{{L}_{-}}$) denotes that Eq. (25) is defined for unstable conditions, i.e., L≤0. For the convective boundary layer, ${u}_{*\mathrm{0}}/\left(-|{f}_{\mathrm{c}}|L\right)$ is generally replaced by the dimensionless inversion height $-{z}_{i}/L$, because the convective ABL depth does not have a significant dependence on ${u}_{*\mathrm{0}}/{f}_{\mathrm{c}}$ (Arya1975). However, we note that $R{o}_{{L}_{-}}$ functions as a “bottom-up” parameter in the non-neutral RANS equation set, with the Obukhov length L in Eq. (16) specified as a surface layer quantity; in effect $R{o}_{{L}_{-}}$ dictates the relative increase in mixing length (i.e., in the dimensionless coordinate $z|{f}_{\mathrm{c}}|/G$). Our length-scale-limited turbulence closures extended to unstable surface layer stratification, as presented in Sect. 3, are dependent on $R{o}_{{L}_{-}}$. This becomes clear when we substitute the mixing-length model extended to unstable surface layer stratification from Eq. (12) into the nondimensional momentum equation from Eq. (21):

$\begin{array}{}\text{(26)}& Ro\frac{\mathrm{d}}{\mathrm{d}{z}^{\prime }}\left({\left[\frac{\mathit{\kappa }{z}^{\prime }}{{\left(\mathrm{1}-{\mathit{\gamma }}_{\mathrm{1}}{z}^{\prime }\mathcal{L}/L\right)}^{-\mathrm{1}/\mathrm{4}}+\mathit{\kappa }{z}^{\prime }\mathcal{L}/{\mathrm{\ell }}_{\mathrm{max}}}\right]}^{\mathrm{2}}\left|\frac{\mathrm{d}{W}^{\prime }}{\mathrm{d}{z}^{\prime }}\right|\frac{\mathrm{d}{W}^{\prime }}{\mathrm{d}{z}^{\prime }}\right)=i{W}^{\prime },\end{array}$

where ℒ∕L is a third nondimensional number, which can also be written as the ratio of two Rossby numbers: $R{o}_{{L}_{-}}/R{o}_{\mathrm{0}}$. For $R{o}_{{L}_{-}}=\mathrm{0}$, the extended models return to the original models. Figure 4 depicts the Rossby number similarity of the extended turbulence closures using six combinations of the three Rossby numbers, which are each simulated with four combinations of G and fc. We use two values of Ro0 (106 and 109) and three values of $R{o}_{{L}_{-}}$ (0, 5×102, and 2×103) for Ro=103. For these Rossby number combinations, $R{o}_{{L}_{-}}=\mathrm{5}×{\mathrm{10}}^{\mathrm{2}}$ and $R{o}_{{L}_{-}}=\mathrm{2}×{\mathrm{10}}^{\mathrm{3}}$ correspond to near-unstable conditions ($-\mathrm{1}/L=\mathrm{0.00125}$–0.005 m−1) and unstable to very unstable conditions ($-\mathrm{1}/L=\mathrm{0.005}$–0.02 m−1), respectively. Figure 4 shows that both extended turbulence closures only depend on Ro0 and $R{o}_{{L}_{-}}$ for a given Ro. Although not shown in Fig. 4, changing Ro would not break the Rossby number similarity. Note that it does not make sense to include combinations of nonzero values of $R{o}_{{L}_{-}}$ that correspond to unstable conditions and large values of Ro that correspond to stable conditions.

Figure 4Rossby number similarity of the turbulence closures extended to unstable surface layer conditions. (a–c) Limited mixing-length model. (d–g) Limited-length-scale kε model.

The extended limited-length-scale mixing-length model (Fig. 4a–c) is less sensitive to $R{o}_{{L}_{-}}$ compared to the extended limited-length-scale kε model (Fig. 4d–g) because of the buoyancy production in the transport equations of k and ε, which is not present in the extended mixing-length model. Both models predict a deeper ABL (larger zi) that is more mixed for stronger unstable surface layer stratification (increasing $R{o}_{{L}_{-}}$). The wind veer is also reduced in stronger unstable conditions for the extended kε model (Fig. 4e), but it does not always decrease with increasing unstable conditions for the extended mixing-length model (Fig. 4b).

One could choose to use the friction velocity at the surface, u*0, as a velocity scale in the Rossby numbers instead of the geostrophic wind speed. However, the friction velocity depends on height z and is a result of the model, not an input. In other words, the height at which the friction velocity needs to be extracted to obtain a collapse is also dependent on the ABL profiles, since the height scales with friction velocity. Hence it is more sensible to use geostrophic wind speed as a velocity scale in the model-based Rossby number similarity – consistent also with classic Ekman theory (which relates the wind speed in terms of G). Nevertheless, it is possible to obtain a Rossby similarity using u*0 as the velocity scale, which is presented in Appendix B.

The Rossby number similarity can be employed to generate a library of ABL profiles for a range of Ro0, Ro, and $R{o}_{{L}_{-}}$. The library contains all possible model solutions for the range of chosen Rossby numbers, and it can be used to determine inflow profiles for three-dimensional RANS simulations, without the need for running one-dimensional precursor simulations.

The obtained Rossby number similarity can only be achieved for a grid-independent numerical setup, as we have shown in Sect. 4.3. In addition, the ambient source terms should also be scaled by the relevant input parameters (G and max), as discussed in Sect. 4.1.

The ABL depth zi predicted by the original limited-length-scale turbulence closures is mainly dependent on the maximum turbulence length scale max. The normalized ABL depth ($\left({z}_{i}+{z}_{\mathrm{0}}\right)|{f}_{\mathrm{c}}|/G$) is mainly dependent on Ro, which is depicted in Fig. 5, where results of the limited-length-scale kε model extended to unstable surface layer stratification are shown for $\mathrm{3}×\mathrm{6}×\mathrm{3}$ combinations of the three Rossby numbers Ro0, Ro, and $R{o}_{{L}_{-}}$. We have chosen G=10 m s−1 and ${f}_{\mathrm{c}}={\mathrm{10}}^{-\mathrm{4}}$ s−1, but the results are independent of G and fc due to the Rossby number similarity. The normalized ABL depth is defined as the height at which the wind direction (relative to the geostrophic wind direction) becomes zero for the second time, i.e., above the mean jet and associated turning as in Ekman theory. For the Ekman solution (Sect. A1), this definition results in an ABL depth equal to ${z}_{i}=\mathrm{2}\mathit{\pi }\sqrt{\mathrm{2}{\mathit{\nu }}_{T}/|{f}_{\mathrm{c}}|}$. The normalized ABL depth in the RANS model increases for stronger unstable surface layer conditions (larger $R{o}_{{L}_{-}}$), i.e., for larger values of the surface heat flux. For neutral and stable conditions ($R{o}_{{L}_{-}}=\mathrm{0}$) and moderate to shallow ABL depths, i.e., $\mathrm{3}×{\mathrm{10}}^{\mathrm{3}}\le R{o}_{\mathrm{\ell }}\le \mathrm{3}×{\mathrm{10}}^{\mathrm{4}}$ – corresponding to ${z}_{i}<\sim \mathrm{2000}$ m as seen in Fig. 5 – we find that ${\mathrm{log}}_{\mathrm{10}}\left(\left[{z}_{i}+{z}_{\mathrm{0}}\right]|{f}_{\mathrm{c}}|/G\right)\propto -a{\mathrm{log}}_{\mathrm{10}}\left(R{o}_{\mathrm{\ell }}\right)$, with a=0.57–0.62 for Ro0 over the range of 109–105. Hence for moderate to shallow ABLs the effective depth modeled in neutral and stable conditions is roughly ${z}_{i}\propto {\mathrm{\ell }}_{\mathrm{max}}^{a}\left(G/|{f}_{\mathrm{c}}|{\right)}^{\mathrm{1}-a},$ with a≈0.6. As seen by the solid lines in Fig. 5, under neutral conditions and with large ABL depths, the zi dependence on max softens ($a<\mathrm{2}/\mathrm{3}$) and deviates from a power law, while for unstable conditions a is similar to the previously found value of 0.6.

Figure 5Normalized boundary layer depth zi predicted by limited-length-scale kε model extended to unstable surface layer stratification, as a function of the three Rossby numbers.

6 Validation and model limits

We employ the Rossby similarity from Sect. 5 to validate a range of results simulated by the original limited-length-scale kε model of including our proposed extension to unstable surface layer stratification. Historical measurements of the geostrophic drag coefficient ${u}_{*\mathrm{0}}/G$ and the cross-isobar angle (the angle between the surface wind direction and the geostrophic wind direction), as summarized by , and measured profiles of the ASL and ABL for different atmospheric stabilities from , respectively, are used as validation metrics. The limited mixing-length model of and its extension are not considered in the comparison with measurements, since we are mainly interested in the kε model.

## 6.1 Geostrophic drag coefficient

The geostrophic drag law (GDL) is a widely used relation in boundary layer meteorology and wind resource assessment (after Troen and Petersen1989), which connects the surface layer properties as z0 and u*0 with the driving forces on top of the ABL proportional to $|{f}_{\mathrm{c}}|G$:

$\begin{array}{}\text{(27)}& G=\frac{{u}_{*\mathrm{0}}}{\mathit{\kappa }}\sqrt{{\left[\mathrm{ln}\left(\frac{{u}_{*\mathrm{0}}}{|{f}_{\mathrm{c}}|{z}_{\mathrm{0}}}\right)-A\right]}^{\mathrm{2}}+{B}^{\mathrm{2}}},\end{array}$

where A and B are empirical constants. The GDL can be derived from Eq. (1), where the Reynolds stresses do not need to be modeled explicitly (as in, for example, Zilitinkevich1989) and can be expressed as an implicit relation for the geostrophic drag coefficient ${u}_{*\mathrm{0}}/G$ and Ro0:

$\begin{array}{}\text{(28)}& \frac{{u}_{*\mathrm{0}}}{G}=\frac{\mathit{\kappa }}{\sqrt{{\left[\mathrm{ln}\left(R{o}_{\mathrm{0}}\right)+\mathrm{ln}\left(\frac{{u}_{*\mathrm{0}}}{G}\right)-A\right]}^{\mathrm{2}}+{B}^{\mathrm{2}}}}.\end{array}$

Figure 6Reproduced from . Geostrophic drag coefficient simulated by the limited-length-scale kε model extended to unstable surface layer stratification, taken at a normalized height of $\left(z+{z}_{\mathrm{0}}\right)|{f}_{\mathrm{c}}|/G=\mathrm{5}×{\mathrm{10}}^{-\mathrm{5}}$ (i.e., in the surface layer), for different Ro0, Ro, and $R{o}_{{L}_{-}}$. Black markers represent measurements from . $R{o}_{\mathrm{\ell }}=\mathrm{3.7}×{\mathrm{10}}^{\mathrm{3}}$ represents max from . Analytic results of and are summarized in Appendix A.

Figure 6 is a reproduction from , where the geostrophic drag coefficient is depicted as a function of surface Rossby number Ro0. The black markers are measurements summarized by , where the dots are near-neutral and near-barotropic conditions, the triangles and squares reflect less idealized atmospheric conditions, and the open circles are measurements with a relative high uncertainty. Results of the limited-length-scale kε model including the extension to unstable surface layer stratification are shown as colored markers, where the colors represent a range of Ro. For the two smallest values of Ro, two additional results are plotted for $R{o}_{{L}_{-}}=\mathrm{5}×{\mathrm{10}}^{\mathrm{2}}$ and $R{o}_{{L}_{-}}=\mathrm{2}×{\mathrm{10}}^{\mathrm{3}}$, representing unstable ($-\mathrm{1}/L=\mathrm{0.005}$ m−1) and very unstable conditions ($-\mathrm{1}/L=\mathrm{0.02}$ m−1) for the chosen values of G=10 m s−1 and ${f}_{\mathrm{c}}={\mathrm{10}}^{-\mathrm{4}}$ s−1. The colored lines are fitted A and B constants from the GDL as defined in Eq. (28). The analytic solutions from and , as summarized in Appendix A, are shown as black and gray lines, respectively. For $R{o}_{{L}_{-}}=\mathrm{0}$, the geostrophic drag coefficient predicted by the limited-length-scale kε model is bounded by the analytic solutions. For Ro→0, the geostrophic drag coefficient of is approximated. For increasing Ro or decreasing ABL depths, the $\mathit{\left\{}{u}_{*}\mathrm{0}/G$, log (Ro0)} relationship becomes more linear. In addition, for $R{o}_{\mathrm{\ell }}=\mathrm{3.7}×{\mathrm{10}}^{\mathrm{3}}$, as used by , and $R{o}_{{L}_{-}}=\mathrm{0}$, most of the near-neutral and near-barotropic measurements are captured quite well. used the measurements of the geostrophic drag coefficient to validate a number of models, which often have only one result for each Ro0. The geostrophic drag coefficients predicted by the limited-length-scale kε model can cover all measurements by varying Ro. In addition, the extension to unstable surface layer conditions can also explain the trend of the more uncertain measurements (black dots). Since Ro and $R{o}_{{L}_{-}}$ influence the ABL depth, as previously shown in Fig. 5, the model suggests that the measurements were conducted for a range of ABL depths that could reflect a range of atmospheric stabilities, although the geostrophic wind shear can play a role here as shown by .

The fitted A and B parameters in Fig. 6 are dependent on Ro and $R{o}_{{L}_{-}}$, which both influence the ABL depth. This is not a surprising result, since many authors have shown that A and B are dependent on atmospheric stability (see, for example, Arya1975; Zilitinkevich1989; Landberg1994). For moderate roughness lengths over land, the measured values tabulated by generally fall between the blue and yellow lines for neutral conditions, which are consistent with the typically used values in wind energy, i.e., A=1.8 and B=4.5 (e.g., Troen and Petersen1989). Assuming max is a measure of the ABL depth, then in the actual atmosphere over land we have $R{o}_{\mathrm{0}}/R{o}_{\mathrm{\ell }}\sim {\mathrm{10}}^{\mathrm{3}}$–105, while over sea the ratio is roughly 106–107. Thus one can see that the typical wind energy values of A and B are a compromise for applicability over both land and sea. The real-world limits mean that the result for Ro=102 (red line) can extend only from Ro0∼105–107, while the over-sea regime (large Ro0) tends to involve a smaller range of Ro. We note that the GDL from Eq. (27) limits how large B can be; generally ${u}_{*\mathrm{0}}/G<\mathit{\kappa }/B$, so values of B greater than those shown are not physical. The model results in Fig. 6 do not violate this limit.

Figure 7Reproduced from . Cross-isobar angle simulated by the limited-length-scale kε model extended to unstable surface layer stratification, taken at a normalized height of $\left(z+{z}_{\mathrm{0}}\right)|{f}_{\mathrm{c}}|/G=\mathrm{5}×{\mathrm{10}}^{-\mathrm{5}}$ for different Ro0, Ro, and $R{o}_{{L}_{-}}$. Black markers represent measurements from . Ro=3.7 represents max from . Analytic results of and are summarized in Appendix A.

## 6.2 Cross-isobar angle

Figure 7 is a reproduction of , where the angle between surface wind direction and the geostrophic wind direction is plotted as a function of the surface Rossby number. This angle is known as the cross-isobar angle, θ0. The black markers, analytic solutions, and model results follow the same definition as used in Fig. 6, where additional black diamond markers are added that correspond to climatological measurements, as discussed by . For $R{o}_{{L}_{-}}=\mathrm{0}$, the model results of the cross-isobar angle are bounded by the analytic solutions, as also found for the geostrophic drag coefficient in Fig. 6. All measurements summarized by can be simulated by the limited-length-scale kε model by varying the ABL depth using Ro. Most of the measurements are well predicted for $R{o}_{{L}_{-}}=\mathrm{0}$ and Ro=103–104, which is the range used by ($R{o}_{\mathrm{\ell }}=\mathrm{3.7}×{\mathrm{10}}^{\mathrm{3}}$). For $R{o}_{{L}_{-}}\ne \mathrm{0}$, smaller values of the cross-isobar angle can be simulated compared with the analytic solution of due to the enhanced rate of mixing. The model cannot predict larger values of the cross-isobar angle compared to the analytic solution of  (45).

## 6.3 Atmospheric surface layer profiles

used measurements of the wind speed components from 10 to 160 m, from The National Test Station for Wind Turbines at Høvsøre, a coastal site in Denmark, characterized as flat grassland. The Coriolis parameter for the test location is $\mathrm{1.21}×{\mathrm{10}}^{-\mathrm{4}}$ s−1. The measurements were taken from sonic anemometers over 1 year, and a wind direction sector was selected to avoid the influence of the coastline and wind turbine wakes. also calculated a “mixing” (turbulence) length scale $\stackrel{\mathrm{^}}{\mathrm{\ell }}$ using a local friction velocity u* and the wind speed gradient:

$\begin{array}{}\text{(29)}& \stackrel{\mathrm{^}}{\mathrm{\ell }}=\frac{{u}_{*}}{\mathrm{d}U/\mathrm{d}z}.\end{array}$

Seven cases were defined based on the atmospheric stability, and these are listed in Table 2 in terms of the Obukhov length, roughness length, and friction velocity. In order to apply the limited-length scale kε, we need to set the geostrophic wind speed and the maximum turbulence length scale, which are both unknown. We choose to use G and max as free parameters, which we fit to a reference wind speed and a turbulence length scale, at a reference height of 60 m. The wind speed gradient is obtained from a central-difference scheme taking the wind speed at 40, 60, and 80 m. The fitted parameters are obtained by running the numerical simulations with a gradients-based optimizer, and the results are listed in Table 2. The maximum max is set to 103 m, which corresponds to an ABL depth on the order of 5 km, as depicted in Fig. 5. The unstable cases are also simulated with the extended limited-length-scale kε model using the measured L and refitted G and max, which are listed in Table 2 as values in parentheses.

Table 2ASL validation cases. Fitted G, fitted max, and modeled u*0 in parentheses represent values for the extended model.

Figure 8ASL measurements of compared to simulation results of the original limited-length-scale kε model of . (a) Wind speed. (b) Turbulence length scale from Eq. (29). Unstable cases are also simulated with our extension to unstable surface layer stratification with L from Table 2.

Figure 8 depicts the wind speed and turbulence length scale of the measurements and numerical simulations using the original and extended limited-length-scale kε models. The turbulence length scale from the numerical simulation is calculated by Eq. (29), instead of the usual definition $\mathrm{\ell }={C}_{\mathit{\mu }}^{\mathrm{3}/\mathrm{4}}{k}^{\mathrm{3}/\mathrm{2}}/\mathit{\epsilon }$. The original limited-length-scale kε model of can capture the wind speed and turbulence length scale for the stable and neutral cases. Note that for the very stable case, the shear is underestimated and the model predicts an ABL depth of about 100 m, which results in a spike in $\stackrel{\mathrm{^}}{\mathrm{\ell }}$, since dU∕dz is zero around the ABL depth. As expected, the original limited-length-scale kε model cannot predict a lower shear and a larger turbulence length scale compared to neutral atmospheric conditions (where $\mathrm{d}U/\mathrm{d}z={u}_{*}/\mathrm{\ell }$ and ℓ=κz), and the optimizer used to fit G and max sets max to our chosen maximum value of 103 m. Note that therefore the lines corresponding to unstable conditions of the original kε model largely overlap in Fig. 8. Higher values of max would not improve the results. The limited-length-scale kε model extended to unstable surface layer stratification is able to predict turbulence length scales larger than ℓ=κz and shows improved results for both the shear and the turbulence length scale.

Table 2 also shows the measured and simulated friction velocity at a height of 10 m. The simulated friction velocity is calculated as ${u}_{*}=\left({\stackrel{\mathrm{‾}}{{u}^{\prime }{w}^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{{v}^{\prime }{w}^{\prime }}}^{\mathrm{2}}{\right)}^{\mathrm{1}/\mathrm{4}}$. For the unstable cases, it is clear that the extended model predicts friction velocities that are closer to the measurements compared to the original limited-length-scale k-ε model due to the enhanced mixing.

It should be noted that the validation presented in Fig. 8 could be considered as a best-possible simulation-to-measurement comparison because we have allowed ourselves to tune both G and max. When G is provided by the measurements, it is more difficult to obtain a good match, as shown in Sect. 6.4.

Table 3ABL validation cases based on .

## 6.4 Atmospheric-boundary-layer profiles

performed lidar measurements of the horizontal wind speed components from 10 to 1200 m at the same test site as discussed in Sect. 6.3. selected 10 cases that differ in geostrophic forcing and atmospheric stability. The cases were selected to challenge the validation of numerical models. Since our numerical setup can only handle a constant geostrophic wind speed, we select the barotropic cases from : cases 4, 5, and 9 and the corresponding values of the Obukhov length, geostrophic wind, roughness length, friction velocity, Ro0, and $R{o}_{{L}_{-}}$ are listed in Table 3. For convenience, we keep the case names as introduced by . Cases 4 and 5 represent a stable and a neutral ABL with high forcing, respectively, where Ro0=107. Case 9 is characterized by a low forcing and very unstable stratification, where $R{o}_{\mathrm{0}}=\mathrm{2.8}×{\mathrm{10}}^{\mathrm{6}}$.

In Case 6 from it is observed that the lidar measurements do not approach the geostrophic wind speed at large heights above the surface. This is because the geostrophic wind speed in is derived from outputs of the Weather Research and Forecasting (WRF) model over a large area, potentially leading to bias. Therefore, we use a slightly different approach to estimate the geostrophic wind; because the wind speed above the ABL is nearly always in geostrophic balance we can just assume the wind speed measured by the wind lidar above the boundary layer depth to be equal to the geostrophic wind speed, thereby avoiding possible prediction errors in wind speed from the WRF model. Instead, only the ABL depth is estimated from the WRF model outputs. The ABL depth is available as a diagnostic variable predicted by the Yonsei University ABL scheme in the WRF model. To be sure that the lidar wind speed is close to the geostrophic wind speed, we always estimate it from the level that is higher than the modeled ABL depth during all 30 min means, which constitute the three cases.

Since G is known, we can use Rossby similarity for model validation. While one could try to find an max to get the best comparison with the measurements, we find that it is difficult to define a good metric. For example, we could attempt to find an max that results in an equivalent ABL depth equal to that of the measurement cases; however, the ABL depth was not directly measured and only estimated from a model. Instead of finding a single max value, we choose to simulate a range of max values. We note that part of this difficulty is due to the limited extent of the model. There is no “top-down” information; i.e., we lack entrainment effects and the impact of the strength of the capping inversion. An extra length scale could be introduced to account for such effects; examples are the nonlocal static stability scale found in , the “mid-ABL” scale of (generalized by for matching G), and the “top-down” scale of .

Figure 9 depicts the measured wind speed and wind direction, for each validation case. Since cases 4 and 5 have the same G (within 5 %) and thus same surface Rossby number Ro0≃107, we can plot them together because the normalized model results are the same for both cases. The error bars represent the standard error of the mean. The original limited-length-scale kε model of is employed with a range of Ro. The unstable ABL case (Case 9) is also simulated with the model extension to unstable surface layer stratification using $R{o}_{{L}_{-}}$ from Table 3 and the two smallest values of Ro. Case 4 has a strong wind shear and a wind veer that leads to a cross-isobar angle of 50. The limited-length-scale kε model can predict a maximum cross-isobar angle of 45 for extremely shallow ABL depths, as shown in Sect. 6.2. Hence, the measured ABL from Case 4 is not a possible numerical solution. The measured ABL from Case 5 can be predicted by the original limited-length-scale kε model, while this is not the case for the wind speed close to the ground of Case 9 due to the strong unstable stratification. When the limited-length-scale kε model including the extension for unstable surface layer conditions is employed, the prediction of the wind speed near the ground is improved, although it is difficult to correctly predict both wind speed and wind direction. It should be noted that the extended (unstable) model only improves the wind speed in the surface layer (at 10 m), noting the dotted and solid lines crossing in Fig. 9c.

Figure 9ABL measurements from compared to simulation results of the original limited-length-scale kε model of , for a range of Ro. (a, c) Wind speed. (b, d) Wind direction. Unstable Case 9 (c, d) is also simulated with our extension to unstable surface layer stratification (dashed lines), with $R{o}_{{L}_{-}}$ (i.e., 1∕L) from Table 3.

From the measurements during Case 9 it was observed that the WRF-modeled ABL depth grew from 300 m to nearly 1200 m, which indicates that the conditions were largely transient; such nonstationary conditions are difficult for a RANS model. More unstable cases are necessary to further validate the extended model, including measurements of turbulence quantities such as the (total) turbulence intensity. It is possible to use validation cases based on turbulence-resolving methods, such as large-eddy simulations, in future work.

7 Conclusions

The idealized ABL was simulated with a one-dimensional RANS solver, using two different turbulence closures: a limited mixing-length model and a limited-length-scale kε model. While these models require four input parameters, we have shown that the simulated ABL profiles collapse to a dependence upon two Rossby numbers, which are defined by the roughness length and the maximum turbulence length scale, respectively. The Rossby number based on the maximum turbulence length scale is a new dimensionless number and is related to the ABL depth. The model-based Rossby number similarity obtained herein is valid for both turbulence models. We have employed Rossby number similarity to compare the range of model solutions with historical measurements of relevant associated meteorological quantities, such as the geostrophic drag coefficient and cross-isobar angle. The measured variation in these measurements can be explained by dependence upon the new Rossby number (i.e., ABL depth). In addition, we have shown how two classic analytic solutions of the idealized ABL (Ekman1905; Ellison1956) act as bounds on the results obtainable by the limited-length-scale kε model.

The limited-length-scale turbulence closures can represent the effects of stable and neutral stratification but cannot model unstable conditions. We have proposed simple extensions to overcome this issue, without adding a temperature equation . The extended models require an additional input, the Obukhov length, which can be used to define a third Rossby number. We have shown that the extension of the kε model compares well with measurements of seven ASL profiles, representing a range of atmospheric stabilities, including three unstable cases. The kε model further offers turbulence intensity, whose profile is also found to collapse according to the developed similarity theories. A model validation of the ABL for a stable, a neutral, and an unstable case is performed, with less success for the non-neutral cases. In the very stable case, the measured wind veer of 50 was larger than the maximum wind veer of 45 that the kε model can simulate. In addition, the very unstable case was characterized by nonstationary conditions, which are difficult to capture with a RANS model. More validation cases based on the convective ABL are necessary to quantify the performance of the turbulence model extension to unstable conditions beyond the surface layer.

The application of the one-dimensional RANS simulations to generate inflow profiles for three-dimensional RANS simulations is not performed here and it should be investigated in future work. Ongoing and future work also includes the incorporation of the effect of the capping-inversion strength to accommodate entrainment at the ABL top (softening the ABL lid, one might say); this can be considered as an introduction of an additional length scale. In addition, the effects of length-scale limitation and neglecting the buoyancy force in the momentum equation need to be quantified for three-dimensional RANS simulations of complex terrain and wind farms.

Appendix A: Analytic solutions for the idealized ABL

## A1 Constant eddy viscosity – Ekman spiral

The analytic solution of , known as the Ekman spiral, can be expressed as a function of a single variable, the normalized height $\mathit{\xi }\equiv z\sqrt{|{f}_{\mathrm{c}}|/\left(\mathrm{2}{\mathit{\nu }}_{T}\right)}$. The wind speed $S=\sqrt{{U}^{\mathrm{2}}+{V}^{\mathrm{2}}}$ and the wind direction θ can be written as

$\begin{array}{}\text{(A1)}& \begin{array}{rl}& S\left(\mathit{\xi }\right)=G\sqrt{\mathrm{1}-\mathrm{2}\mathrm{cos}\left(\mathit{\xi }\right)\mathrm{exp}\left(-\mathit{\xi }\right)+\mathrm{exp}\left(-\mathrm{2}\mathit{\xi }\right)},\\ & \mathit{\theta }\left(\mathit{\xi }\right)=\mathrm{arctan}\left(\frac{\mathrm{sin}\left(\mathit{\xi }\right)}{\mathrm{exp}\left(\mathit{\xi }\right)-\mathrm{cos}\left(\mathit{\xi }\right)}\right).\end{array}\end{array}$

The cross-isobar angle θ(0) is found to be 45, and the geostrophic drag coefficient is zero (since there is no roughness or u* within Ekman theory).

## A2 Linear eddy viscosity – Ellison

The analytic solution of for the U and V velocity components can be written in terms of the Kelvin functions ker and kei, as discussed by :

$\begin{array}{}\text{(A2)}& U=cG\text{ker}\left(x\right)+{U}_{\mathrm{G}},\phantom{\rule{0.25em}{0ex}}V=cG\text{kei}\left(x\right)+{V}_{\mathrm{G}},\end{array}$

where x is a normalized height $x\equiv \mathrm{2}\sqrt{z|{f}_{\mathrm{c}}|/\left(\mathit{\kappa }{u}_{*\mathrm{0}}\right)}$ and c is a constant. For zz0 (and assuming ${z}_{\mathrm{0}}\ll {u}_{*\mathrm{0}}/|{f}_{c}|$), the Kelvin functions can be expanded, and the solution can be written as

$\begin{array}{}\text{(A3)}& \begin{array}{rl}& U\approx -cG\left[\frac{\mathrm{1}}{\mathrm{2}}\mathrm{ln}\left(\frac{{z}_{\mathrm{0}}|{f}_{\mathrm{c}}|}{\mathit{\kappa }{u}_{*\mathrm{0}}}\right)+{\mathit{\gamma }}_{\mathrm{e}}\right]+{U}_{\mathrm{G}}=\mathrm{0},\\ & V\approx -cG\frac{\mathit{\pi }}{\mathrm{4}}+{V}_{\mathrm{G}}=\mathrm{0},\end{array}\end{array}$

where γe≈0.57721 is the Euler–Mascheroni constant. We can set the geostrophic wind G through the constant c:

$\begin{array}{}\text{(A4)}& \begin{array}{rl}& c=-{\left({\left[\frac{\mathrm{1}}{\mathrm{2}}\mathrm{ln}\left(\frac{{z}_{\mathrm{0}}|{f}_{\mathrm{c}}|}{\mathit{\kappa }{u}_{*\mathrm{0}}}\right)+{\mathit{\gamma }}_{\mathrm{e}}\right]}^{\mathrm{2}}+\frac{{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{16}}\right)}^{-\mathrm{1}/\mathrm{2}},\\ & {U}_{\mathrm{G}}=cG\left[\frac{\mathrm{1}}{\mathrm{2}}\mathrm{ln}\left(\frac{{z}_{\mathrm{0}}|{f}_{\mathrm{c}}|}{\mathit{\kappa }{u}_{*\mathrm{0}}}\right)+{\mathit{\gamma }}_{\mathrm{e}}\right],\phantom{\rule{0.25em}{0ex}}{V}_{\mathrm{G}}=cG\frac{\mathit{\pi }}{\mathrm{4}}.\end{array}\end{array}$

Note that chose $cG=-\mathrm{2}{u}_{*\mathrm{0}}/\mathit{\kappa }$ (so his c is 5 times the geostrophic drag coefficient ${u}_{*\mathrm{0}}/G$ for κ=0.4), which follows from the Neumann condition,

$\begin{array}{}\text{(A5)}& \frac{\mathrm{d}U}{\mathrm{d}z}=\frac{{u}_{*\mathrm{0}}}{\mathit{\kappa }z}=-\frac{cG}{\mathrm{2}z},\end{array}$

by taking dU∕dz in Eq. (A2) for zz0. As as consequence, the geostrophic wind becomes a dependent parameter. We prefer to keep the geostrophic wind as an independent parameter by using c as defined in Eq. (A4). Then, the effective u*0 is calculated as ${u}_{*\mathrm{0},\mathrm{eff}}=cG\mathit{\kappa }/\mathrm{2}$.

A GDL can be derived in the form of Eq. (28) by using the Neumann conditions of Eq. (A5) and the constant c from Eq. (A4), where $A=-\mathrm{ln}\left(\mathit{\kappa }\right)+\mathrm{2}{\mathit{\gamma }}_{\mathrm{e}}\approx \mathrm{2.07}$ and $B=\mathit{\pi }/\mathrm{2}\approx \mathrm{1.57}$, as also shown by . The friction velocity in Eq. (A4) can now be calculated by solving the GDL for ${u}_{*\mathrm{0}}/G$. Hence, the analytic solution of is only dependent on Ro0.

The cross-isobar angle (angle between the geostrophic wind direction and surface wind direction) can be written as a function of the geostrophic drag coefficient ${u}_{*\mathrm{0}}/G$ and the Rossby number Ro0 using Eq. (A4):

$\begin{array}{}\text{(A6)}& {\mathit{\theta }}_{\mathrm{0}}=\mathrm{arctan}\left(\frac{{V}_{\mathrm{G}}}{{U}_{\mathrm{G}}}\right)=\mathrm{arctan}\left(\frac{\mathit{\pi }/\mathrm{2}}{\mathrm{2}{\mathit{\gamma }}_{\mathrm{e}}-\mathrm{ln}\left(R{o}_{\mathrm{0}}\right)-\mathrm{ln}\left(\mathit{\kappa }{u}_{*\mathrm{0}}/G\right)}\right),\end{array}$

where the GDL can be used to solve for ${u}_{*\mathrm{0}}/G$.

Appendix B: Rossby number similarity based on the friction velocity

In this article, we have shown a Rossby similarity of two limited-length-scale turbulence closures using the geostrophic wind speed, G, as a velocity scale, instead of the friction velocity near the ground, u0*. It is more convenient to use G because it is a constant and a model input, while u0* is a model result that depends on height. However, it is possible to obtain a Rossby similarity based on u0* using the geostrophic drag coefficient ${u}_{*\mathrm{0}}/G$ from Fig. 6, since we can write the following:

$\begin{array}{}\text{(B1)}& \begin{array}{rl}& R{o}_{\mathrm{\ell }}^{*}\equiv \frac{{u}_{*\mathrm{0}}}{|{f}_{\mathrm{c}}|{\mathrm{\ell }}_{\mathrm{max}}}=\frac{{u}_{*\mathrm{0}}}{G}R{o}_{\mathrm{\ell }},\\ & R{o}_{\mathrm{0}}^{*}\equiv \frac{{u}_{*\mathrm{0}}}{|{f}_{\mathrm{c}}|{z}_{\mathrm{0}}}=\frac{{u}_{*\mathrm{0}}}{G}R{o}_{\mathrm{0}},\\ & R{o}_{{L}_{-}}^{*}\equiv \frac{-{u}_{*\mathrm{0}}}{|{f}_{\mathrm{c}}|L}=\frac{{u}_{*\mathrm{0}}}{G}R{o}_{{L}_{-}}.\end{array}\end{array}$

Figure 6 can be transformed into an explicit relation of ${u}_{*\mathrm{0}}/G$ as a function of $R{o}_{\mathrm{\ell }}^{*}$, $R{o}_{\mathrm{0}}^{*}$, and $R{o}_{{L}_{-}}$; the result is depicted in Fig. B1 for $R{o}_{{L}_{-}}^{*}=\mathrm{0}$.

The Rossby similarity based on u0* is illustrated in Fig. B2 for four combinations of $R{o}_{\mathrm{0}}^{*}$ (105 and 108) and $R{o}_{\mathrm{\ell }}^{*}$ (102 and 104) for $R{o}_{{L}_{-}}=\mathrm{0}$, using four combinations of u*0 (0.2 and 0.4 m s−1) and fc ($\mathrm{5}×{\mathrm{10}}^{-\mathrm{5}}$ and 10−4 s−1). Only results of the limited-length-scale kε model are shown for brevity, although the Rossby similarity based on u0* also applies to the limited mixing-length model and for the unstable extension (where $R{o}_{{L}_{-}}\ne \mathrm{0}$).

It should be noted that u0* in Fig. 6 was extracted at a normalized height of $\left(z+{z}_{\mathrm{0}}\right)|{f}_{\mathrm{c}}|/G=\mathrm{5}×{\mathrm{10}}^{-\mathrm{5}}$, which represents the surface layer. If a perfect Rossby similarity based on u0* is desired, one would need to extract u0* at a constant normalized height equal to $\left(z+{z}_{\mathrm{0}}\right)|{f}_{\mathrm{c}}|/{u}_{\mathrm{0}*}$, which requires an iterative process of finding a geostrophic wind speed that results in a RANS simulation with a desired u0*, at a constant normalized height. This is beyond the scope of the present work.

Figure B1Geostrophic drag coefficient simulated by the limited-length-scale k-ε model as a function of $R{o}_{\mathrm{\ell }}^{*}$ and $R{o}_{\mathrm{0}}^{*}$ with $R{o}_{{L}_{-}}^{*}=\mathrm{0}$.

Figure B2Rossby number similarity of the limited-length-scale kε model using the friction velocity as the velocity scale, for $R{o}_{{L}_{-}}=\mathrm{0}$.

Code and data availability
Code and data availability.

The numerical results are generated with proprietary software, although the data presented can be made available by contacting the corresponding author.

Author contributions
Author contributions.

MPVDL performed the simulations, obtained the model-based Rossby number similarity for the kε model, produced all figures, and drafted the article. MPVDL and MK proposed the extension to unstable conditions. MK added connections and relations to meteorological theory and interpretations. AP proposed the Rossby number similarity of the mixing-length model. AP and RF contributed to the model validation. All authors contributed to the methodology and finalization of the paper.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Review statement
Review statement.

This paper was edited by Johan Meyers and reviewed by Javier Sanz Rodrigo and one anonymous referee.

References

Apsley, D. D. and Castro, I. P.: A limited-length-scale kε model for the neutral and stably-stratified atmospheric boundary layer, Bound.-Lay. Meteorol., 83, 75–98, https://doi.org/10.1023/A:1000252210512, 1997. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v

Arya, S. P. S.: Geostrophic drag and heat transfer relations for the atmospheric boundary layer, Q. J. Roy. Meteorol. Soc., 101, 147–161, 1975. a, b, c

Arya, S. P. S. and Wyngaard, J. C.: Effect of baroclinicity on wind profiles and the geostrophic drag law for the convective boundary layer, J. Atmos. Sci., 32, 767–778, 1975. a

Blackadar, A. K.: The vertical distribution of wind and turbulent exchange in a neutral atmosphere, J. Geophys. Res., 67, 3095–3102, 1962. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u

Boussinesq, M. J.: Théorie de l'écoulement tourbillonnant et tumultueux des liquides, Gauthier-Villars et fils, Paris, France, 1897. a

Constantin, A. and Johnson, R. S.: Atmospheric Ekman Flows with Variable Eddy Viscosity, Bound.-Lay. Meteorol., 170, 395–414, https://doi.org/10.1007/s10546-018-0404-0, 2019. a, b

Dyer, A. J.: A review of flux-profile relationships, Bound.-Lay. Meteorol., 7, 363–372, https://doi.org/10.1007/BF00240838, 1974. a, b

Ekman, V. W.: On the influence of the earth's rotation on ocean-currents, Arkiv Mat. Astron. Fysik, 2, 11, 1905. a, b, c, d, e, f, g, h, i

Ellison, T. H.: Atmospheric Turbulence in Surveys of mechanics, Cambridge University Press, Cambridge, UK, 1956. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Floors, R., Peña, A., and Gryning, S.-E.: The effect of baroclinicity on the wind in the planetary boundary layer, Q. J. Roy. Meteorol. Soc., 141, 619–630, https://doi.org/10.1002/qj.2386, 2015. a

Gryning, S.-E., Batchvarova, E., Brümmer, B., Jørgensen, H., and Larsen, S.: On the extension of the wind profile over homogeneous terrain beyond the surface boundary layer, Bound.-Lay. Meteorol., 124, 371–379, https://doi.org/10.1007/s10546-007-9166-9, 2007. a

Hess, G. D. and Garratt, J. R.: Evaluating Models Of The Neutral, Barotropic Planetary Boundary Layer Using Integral Measures: Part II. Modelling Observed Conditions, Bound.-Lay. Meteorol., 104, 359–369, https://doi.org/10.1023/A:1016525332683, 2002. a, b, c, d, e, f, g, h, i, j, k, l

Hong, S.-Y., Noh, Y., and Dudhia, J.: A New Vertical Diffusion Package with an Explicit Treatment of Entrainment Processes, Mon. Weather. Rev., 134, 2318–2341, https://doi.org/10.1175/MWR3199.1, 2006. a

Kazanskii, A. and Monin, A.: Dynamic interaction between atmosphere and surface of earth, Akademiya Nauk Sssr Izvestiya Seriya Geofizicheskaya, 24, 786–788, 1961. a

Kelly, M. and Troen, I.: Probabilistic stability and “tall” wind profiles: theory and method for use in wind resource assessment, Wind Energy, 19, 227–241, 2016. a

Kelly, M. C., Cersosimo, R. A., and Berg, J.: A universal wind profile for the inversion-capped neutral atmospheric boundary layer, Q. J. Roy. Meteorol. Soc., 145, 982–992, https://doi.org/10.1002/qj.3472, 2019. a

Koblitz, T., Bechmann, A., Sogachev, A., Sørensen, N., and Réthoré, P.-E.: Computational Fluid Dynamics model of stratified atmospheric boundary-layer flow, Wind Energy, 18, 75–89, https://doi.org/10.1002/we.1684, 2015. a, b, c

Krishna, K.: The planetary-boundary-layer model of Ellison (1956) – A retrospect, Bound.-Lay. Meteorol., 19, 293–301, https://doi.org/10.1007/BF00120593, 1980. a, b, c, d

Landberg, L.: Short-term prediction of local wind conditions, Tech. rep., Risø National Laboratory, Roskilde, Denmark, 1994. a

Lettau, H.: A Re‐examination of the “Leipzig Wind Profile” Considering some Relations between Wind and Turbulence in the Frictional Layer, Tellus, 2, 125–129, https://doi.org/10.1111/j.2153-3490.1950.tb00321.x, 1950. a

Lettau, H. H.: Wind Profile, Surface Stress and Geostrophic Drag Coefficients in the Atmospheric Surface Layer, Adv. Geophys., 6, 241–257, 1959. a, b

Michelsen, J. A.: Basis3D – a platform for development of multiblock PDE solvers, Tech. Rep. AFM 92-05, Technical University of Denmark, Lyngby, Denmark, 1992. a

Monin, A. S. and Obukhov, A. M.: Basic laws of turbulent mixing in the surface layer of the atmosphere, Tr. Akad. Nauk. SSSR Geophiz. Inst., 24, 163–187, 1954. a

Peña, A., Gryning, S.-E., and Mann, J.: On the length-scale of the wind profile, Q. J. Royal Meteorol. Soc., 136, 2119–2131, https://doi.org/10.1002/qj.714, 2010. a, b

Peña, A., Floors, R., and Gryning, S.-E.: The Høvsøre Tall Wind-Profile Experiment: A Description of Wind Profile Observations in the Atmospheric Boundary Layer, Bound.-Lay. Meteorol., 150, 69–89, https://doi.org/10.1007/s10546-013-9856-4, 2014. a, b, c, d, e, f, g, h, i, j, k

Richards, P. J. and Hoxey, R. P.: Appropriate boundary conditions for computational wind engineering models using the kε turbulence model, J. Wind Eng. Indust. Aerodynam., 46, 145–153, 1993. a, b

Sogachev, A., Kelly, M., and Leclerc, M. Y.: Consistent Two-Equation Closure Modelling for Atmospheric Research: Buoyancy and Vegetation Implementations, Bound.-Lay. Meteorol., 145, 307–327, https://doi.org/10.1007/s10546-012-9726-5, 2012. a, b, c

Sørensen, N. N.: General purpose flow solver applied to flow over hills, Ph.D. thesis, Risø National Laboratory, Roskilde, Denmark, 1994. a

Sørensen, N. N., Bechmann, A., Johansen, J., Myllerup, L., Botha, P., Vinther, S., and Nielsen, B. S.: Identification of severe wind conditions using a Reynolds Averaged Navier–Stokes solver, J. Phys.: Conf. Ser., 75, 1–13, https://doi.org/10.1088/1742-6596/75/1/012053, 2007. a

Spalart, P. and Rumsey, C.: Effective inflow conditions for turbulence models in aerodynamic calculations, AIAA J., 45, 2544–2553, 2007. a

Sumner, J. and Masson, C.: The Apsley and Castro Limited-Length-Scale kε Model Revisited for Improved Performance in the Atmospheric Surface Layer, Bound.-Lay. Meteorol., 144, 199–215, https://doi.org/10.1007/s10546-012-9724-7, 2012. a, b

Troen, I. and Petersen, E. L.: European Wind Atlas, Risø National Laboratory, Roskilde, Denmark, 1989. a, b

van der Laan, M. P. and Sørensen, N. N.: A 1D version of EllipSys, Tech. Rep. DTU Wind Energy E-0141, Technical University of Denmark, Roskilde, Denmark, 2017a.  a

van der Laan, M. P. and Sørensen, N. N.: Why the Coriolis force turns a wind farm wake clockwise in the Northern Hemisphere, Wind Energ. Sci., 2, 285–294, https://doi.org/10.5194/wes-2-285-2017, 2017b. a

van der Laan, M. P., Kelly, M. C., and Sørensen, N. N.: A new kϵ model consistent with Monin–Obukhov similarity theory, Wind Energy, 20, 479–489, https://doi.org/10.1002/we.2017, 2017. a, b

Zilitinkevich, S. S.: Velocity profiles, the resistance law and the dissipation rate of mean flow kinetic energy in a neutrally and stably stratified planetary boundary layer, Bound.-Lay. Meteorol., 46, 367–387, 1989. a, b, c

Zilitinkevich, S. S. and Esau, I. N.: Resistance and heat-transfer laws for stable and neutral planetary boundary layers: old theory advanced and re-evaluated, Q. J. Roy. Meteorol. Soc., 131, 1863–1892, 2005. a

From the two-equation kε model (which is isotropic), the total turbulence intensity is calculated by $I=\sqrt{\mathrm{2}/\mathrm{3}k}/\sqrt{{U}^{\mathrm{2}}+{V}^{\mathrm{2}}}$.