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

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 Reynoldsaveraged 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.


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-termaveraged 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 (Koblitz et al., 2015) and wind farms (van der Laan and Sørensen, 2017b). Turbulence is modeled here by two limited-length-scale turbulence closures, the mixing-length model of Blackadar (1962) and the twoequation k-ε model of Apsley and Castro (1997). 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 maxi-mum 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 Lettau (1959), 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 Blackadar (1962) and Apsley and Castro (1997) 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 Blackadar (1962). 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 onedimensional 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.

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: where U and V are the mean horizontal velocity components, U G and V G are the corresponding mean geostrophic velocities, f c = 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 u w and v w are modeled by the linear stress-strain relationship of Boussinesq (1897): u w = −ν T dU/dz and v w = −ν T dV /dz, where ν T is the eddy viscosity. The boundary conditions for U and V are U = V = 0 at z = z 0 and U = U G and V = V G at z → ∞, where z 0 is the roughness length. Note that it is possible to write the two momentum equations as a single ordinary differential equation: where The eddy viscosity, ν T , needs to be modeled in order to close the system of equations. The eddy viscosity can be written as ν T = u * , 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 (Ekman, 1905), 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 Ellison (1956) and discussed by Krishna (1980) and Constantin and Johnson (2019). The two analytic solutions are provided in Appendix A. One can relate the analytic solution of Ellison (1956) to the (neutral) ASL (z z i ), while the Ekman spiral is more valid for altitudes around the ABL depth z i . 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 z z i and a constant eddy viscosity for z ∼ z i should provide a more realistic solution. For example, the eddy viscosity could have the form ν T = κu * 0 z exp(−z/ h), where ν T increases linearly with height for z h 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. Constantin and Johnson (2019) 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 Wind Energ. Sci., 5, 355-374, 2020 www.wind-energ-sci.net/5/355/2020/ necessary, and one of the simplest numerical models for the idealized ABL is given by Blackadar (1962) using Prandtl's mixing-length model: where S = (dU/dz) 2 + (dV /dz) 2 = |dW/dz| is the magnitude of the strain-rate tensor, and is prescribed as a turbulence length scale, 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: 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: where P is the turbulence production, and σ k , σ ε , C ε,1 , and C ε,2 are model constants that should follow the relationship κ 2 = σ ε C µ (C ε,2 − C ε,1 ). When using the standard k-ε model calibrated for atmospheric flows (Richards and Hoxey, 1993), the turbulence length scale or eddy viscosity will keep increasing until a boundary layer depth is formed and the analytic solution of Ellison (1956) is approximated. Apsley and Castro (1997) 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 * ε,1 : where the turbulence length scale is modeled as = C 3/4 µ k 3/2 /ε. This limited-length-scale k-ε model behaves very similarly to the mixing-length model of Blackadar (1962) (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 * ε,1 P ∼ C ε,1 ε), and the turbulence length scale is limited. For a given z 0 , G, and f c , the ABL depth can be controlled by max . This means that max is related to z i ; Apsley and Castro (1997) noted that max ∼ z i /3 applies to typical neutral ABLs. However, the simulated boundary layer depth using the k-ε model of Apsley and Castro (1997) has an approximate dependence of z i ∝ (G/|f c |) 1−a a max 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 Ekman (1905) and Ellison (1956) with the numerical solutions of the limited mixing-length model of Blackadar (1962) and the limitedlength-scale k-ε of Apsley and Castro (1997) in terms of wind speed, wind direction, θ = arctan(V /U ), 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 c = 10 −4 s −1 , G = 10 m s −1 , and z 0 = 10 −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-lengthscale 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 Apsley and Castro (1997) is that of the Ekman spiral: 45 • . For large max values, the numerical solution approximates the analytic solution of Ellison (1956) but does not match it because their eddy viscosities are different for z ≥ z i .

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.

Limited mixing-length model
One can generically parameterize the turbulence length scale as a "parallel" combination of ASL and ABL scales, Blackadar (1962) chose ASL = κz and ABL = max to arrive at Eq. (4). If we choose to set www.wind-energ-sci.net/5/355/2020/ Wind Energ. Sci., 5, 355-374, 2020   Blackadar (1962) Limited-length-scale k-ε model ν T = C µ k 2 /ε = C 3/4 µ k 3/2 /ε Numerical Apsley and Castro (1997) following the turbulence length scale that is a result of Monin-Obukhov Similarity Theory (MOST, Monin and Obukhov, 1954) 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 Blackadar (1962) to unstable surface layer stratification, as Approaching neutral conditions, L −1 → 0, the original length-scale model of Blackadar (1962) is obtained. Note that in stable conditions, φ m = 1+β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 −1 Thus we can simply use the original length-scale model of Blackadar (1962) for stable and neutral conditions; the stable φ m function simply informs the selection of max,eff , following Eq. (13).

Limited-length-scale k-ε model
Sumner and Masson (2012) argued that, in stable conditions, the limited-length-scale k-ε model of Apsley and Castro (1997) overpredicts in the surface layer compared to MOST, where max = Lκ/β and β ≈ 5. They proposed a more complicated expression for C * ε,1 in the transport equation of ε compared to the original model of Apsley and Castro (1997). Sogachev et al. (2012) 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 Sumner and Masson (2012) provides a better match of the turbulence length scale within the surface layer compared to MOST with respect to the original k-ε model of Apsley and Castro (1997). How-Wind Energ. Sci., 5, 355-374, 2020 www.wind-energ-sci.net/5/355/2020/ ever, 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 Apsley and Castro (1997) by simply reducing max by roughly 20 %. Therefore, we choose to use the model of Apsley and Castro (1997) 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 B in the k-ε transport equations: Here B is modeled as following MOST, using the similarity functions of Dyer (1974) as discussed in van der Laan et al. (2017). We use the gachev et al. (2012), which for unstable conditions includes the prescription amenable to the free-convection limit: ε/B → 1 for P/B → 0. Further, α B → 1 as → 0, matching neutral conditions since z/L also vanishes then. The prescription (Eq. 17) results in which also means that C * ε,3 → C ε,2 approaching the effective ABL top ( → max ), so sources and sinks of ε balance in Eq. (15), i.e., P −ε +B, all have the same coefficient C ε,2 .

Methodology of numerical simulations
The one-dimensional numerical simulations are performed with EllipSys1D (van der Laan and Sørensen, 2017a), which is a simplified one-dimensional version of EllipSys3D, initially developed by Sørensen (1994) and Michelsen (1992). 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.

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, Koblitz et al. (2015) chose to set upper limits for k and ε to prevent numerical instabilities, we prefer a more physical method, including ambient source terms S k,amb and S ε,amb in the k and ε transport equations, respectively. Following Spalart and Rumsey (2007), we set When all sources of turbulence are zero (P = B = 0) and the diffusion terms are zero (dk/dz = dε/dz = 0), then k = k amb and ε = ε amb . To be consistent with the equations solved, we define the ambient turbulence quantities in terms of the driving parameters, G and max : Here I amb is the total turbulence intensity 1 above the (simulated) ABL, and C amb is the ratio of the turbulence length scale above the ABL ( amb ) to maximum turbulence length scale ( max ). We choose small values for I amb = 10 −6 and C amb = 10 −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. Sogachev et al. (2012) and Koblitz et al. (2015) 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.

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 1/|f c | s. All spatial gradients are discretized by a second-order central-difference scheme. Convergence is typically achieved after 10 5 iterations, which takes about 10 s on a single 2.7 GHz CPU. The domain height is set to 10 5 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 Sørensen et al. (2007). 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 + z 0 . 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.

Grid refinement study
A grid refinement study of the numerical setup is performed for the limited-length-scale k-ε model of Apsley and Castro (1997), using 48, 96, 192, 384, and 768 cells. We choose f c = 10 −4 s −1 , z 0 = 10 −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 lowlevel 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 Blackadar (1962). 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.

Rossby number similarity in numerical and analytical solutions
The numerical solution of the original limited-length-scale turbulence closures of Blackadar (1962) and Apsley and Castro (1997) depend on four parameters: f c (s −1 ), G (m s −1 ), max (m), and z 0 (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 ≡ W/U, ν T 0 ≡ ν T /(UL) and z 0 ≡ z/L, where U and L are characteristic velocity and length scales, respectively: Here, Ro is the Rossby number, Ro = U/(|f c |L), which describes the ratio of the inertial (advective) tendency to the Coriolis force. If we apply the original mixing-length model of Blackadar (1962) for ν T using Eqs. (3) and (4), then Eq. (21) can be written as where −L/ max , is a second dimensionless number. If we choose U = G and L = z 0 , we may define two Rossbylike numbers, with characteristic length scales based on z 0 and max , respectively: Here, we have obtained Ro by rewriting the second dimensionless number L/ max as the ratio of the two Rossby numbers: . Ro 0 is known as the surface Rossby number, first introduced by Lettau (1959); it also resembles a ratio of (inertial) boundary layer depth to z 0 . Analogously, Ro is like the ratio of two boundary layer depths, f c /u * 0 and z i (e.g., Arya and Wyngaard, 1975); here max is a proxy for z i , acting as a "lid" for the ABL. Considering Eq. (23), we have reduced the number of dependent parameters from four to two: f (f c , G, max , z 0 ) → f (Ro 0 , Ro ). For a fixed surface roughness z 0 , the ratio of the two Rossby numbers is then the only dependent parameter: i.e., the ratio of simulated ABL depth to z 0 is the lone parameter. Blackadar (1962) found a characteristic maximum ABL turbulence length scale of 0.00027 G/|f c | for the Leipzig Wind Energ. Sci., 5, 355-374, 2020 www.wind-energ-sci.net/5/355/2020/ wind profile (Lettau, 1950), 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 Blackadar (1962), Fig. 3a-c, and Apsley and Castro (1997), Fig. 3dg. Four combinations of Ro 0 (10 6 and 10 9 ) and Ro (10 3 and 10 5 ) are used, each simulated with four combinations of G (10 and 20 m s −1 ) and f c (5 × 10 −5 and 10 −4 s −1 ). The roughness length and maximum turbulence length scale follow from Eq. (23) and cover a wide range of z 0 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 Ro 0 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 Apsley and Castro (1997) 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 Ro 0 , while further away, it is mainly influenced by Ro with a weaker dependence on Ro 0 .
Considering the non-neutral ABL with Coriolis effects but ignoring the strength of capping inversion (entrainment), in the micrometeorological literature the Kazanskii and Monin (1961) parameter u * 0 /(|f c |L) is typically invoked (e.g., Arya, 1975;Zilitinkevich, 1989). This can also be considered as a third Rossby number, which in our context of using G instead of u * 0 is here the subscript ( L − ) denotes that Eq. (25) is defined for unstable conditions, i.e., L ≤ 0. For the convective boundary layer, u * 0 /(−|f c |L) is generally replaced by the dimensionless inversion height −z i /L, because the convective ABL depth does not have a significant dependence on u * 0 /f c (Arya, 1975). However, we note that Ro 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 Ro L − dictates the relative increase in mixing length (i.e., in the dimensionless coordinate z|f c |/G). Our length-scale-limited turbulence closures extended to unstable surface layer stratification, as presented in Sect. 3, are dependent on Ro 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): where L/L is a third nondimensional number, which can also be written as the ratio of two Rossby numbers: Ro L − /Ro 0 . For Ro L − = 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 f c . We use two values of Ro 0 (10 6 and 10 9 ) and three values of Ro L − (0, 5×10 2 , and 2×10 3 ) for Ro = 10 3 . For these Rossby number combinations, Ro L − = 5×10 2 and Ro L − = 2×10 3 correspond to near-unstable conditions (−1/L = 0.00125-0.005 m −1 ) and unstable to very unstable conditions (−1/L = 0.005-0.02 m −1 ), respectively. Figure 4 shows that both extended turbulence closures only www.wind-energ-sci.net/5/355/2020/ Wind Energ. Sci., 5, 355-374, 2020 depend on Ro 0 and Ro 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 Ro L − that correspond to unstable conditions and large values of Ro that correspond to stable conditions. The extended limited-length-scale mixing-length model (Fig. 4a-c) is less sensitive to Ro 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 z i ) that is more mixed for stronger unstable surface layer stratification (increasing Ro 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.
Wind Energ. Sci., 5, 355-374, 2020 www.wind-energ-sci.net/5/355/2020/ The Rossby number similarity can be employed to generate a library of ABL profiles for a range of Ro 0 , Ro , and Ro 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 z i predicted by the original limited-lengthscale turbulence closures is mainly dependent on the maximum turbulence length scale max . The normalized ABL depth ((z i + z 0 )|f c |/G) is mainly dependent on Ro , which is depicted in Fig. 5, where results of the limited-lengthscale k-ε model extended to unstable surface layer stratification are shown for 3 × 6 × 3 combinations of the three Rossby numbers Ro 0 , Ro , and Ro L − . We have chosen G = 10 m s −1 and f c = 10 −4 s −1 , but the results are independent of G and f c 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 2ν T /|f c |. The normalized ABL depth in the RANS model increases for stronger unstable surface layer conditions (larger Ro L − ), i.e., for larger values of the surface heat flux. For neutral and stable conditions (Ro L − = 0) and moderate to shallow ABL depths, i.e., 3 × 10 3 ≤ Ro ≤ 3 × 10 4 -corresponding to z i <∼ 2000 m as seen in Fig. 5 -we find that log 10 ([z i +z 0 ]|f c |/G) ∝ −alog 10 (Ro ), with a = 0.57-0.62 for Ro 0 over the range of 10 9 -10 5 . Hence for moderate to shallow ABLs the effective depth modeled in neutral and stable conditions is roughly z i ∝ a max (G/|f c |) 1−a , with a ≈ 0.6. As seen by the solid lines in Fig. 5, under neutral conditions and with large ABL depths, the z i dependence on max softens (a < 2/3) and deviates from a power law, while for unstable conditions a is similar to the previously found value of 0.6.

Validation and model limits
We employ the Rossby similarity from Sect. 5 to validate a range of results simulated by the original limited-lengthscale k-ε model of Apsley and Castro (1997) including our proposed extension to unstable surface layer stratification. Historical measurements of the geostrophic drag coefficient u * 0 /G and the cross-isobar angle (the angle between the surface wind direction and the geostrophic wind direction), as summarized by Hess and Garratt (2002), and measured profiles of the ASL and ABL for different atmospheric stabilities from Peña et al. (2010Peña et al. ( , 2014, respectively, are used as validation metrics. The limited mixing-length model of Blackadar (1962) and its extension are not considered in the comparison with measurements, since we are mainly interested in the k-ε model.

Geostrophic drag coefficient
The geostrophic drag law (GDL) is a widely used relation in boundary layer meteorology and wind resource assessment (after Troen and Petersen, 1989), which connects the surface layer properties as z 0 and u * 0 with the driving forces on top of the ABL proportional to |f c |G: 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, Zilitinkevich, 1989) and can be expressed as an implicit relation for the geostrophic drag coefficient u * 0 /G and Ro 0 : (28) Figure 6 is a reproduction from Hess and Garratt (2002), where the geostrophic drag coefficient is depicted as a function of surface Rossby number Ro 0 . The black markers are measurements summarized by Hess and Garratt (2002), 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-lengthscale 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 Ro L − = 5 × 10 2 and Ro L − = 2 × 10 3 , representing unstable (−1/L = 0.005 m −1 ) and very unstable conditions (−1/L = 0.02 m −1 ) for the chosen values of G = 10 m s −1 and f c = 10 −4 s −1 . The colored lines are fitted A and B constants from the GDL as defined in Eq. (28). The analytic solutions from Ekman (1905) and Ellison (1956), as summarized in Appendix A, are shown as black and gray lines, respectively. For Ro L − = 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 Ellison (1956) is approximated. For increasing Ro or decreasing ABL depths, the {u * 0/G, log(Ro 0 )} relationship becomes more linear. In addition, for Ro = 3.7 × 10 3 , as used by Blackadar (1962), and Ro L − = 0, most of the near-neutral and near-barotropic measurements are captured quite well. Hess and Garratt (2002) used the measurements of the geostrophic drag coefficient to validate a number of models, which often have only one result for each Ro 0 . The geostrophic drag coefficients predicted by the limited-lengthscale 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 Ro 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 Floors et al. (2015). The fitted A and B parameters in Fig. 6 are dependent on Ro and Ro 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, Arya, 1975;Zilitinkevich, 1989;Landberg, 1994). For moderate roughness lengths over land, the measured values tabulated by Hess and Garratt (2002) 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 Petersen, 1989). Assuming max is a measure of the ABL depth, then in the actual atmosphere over land we have Ro 0 /Ro ∼ 10 3 -10 5 , while over sea the ratio is roughly 10 6 -10 7 . 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 = 10 2 (red line) can extend only from Ro 0 ∼ 10 5 -10 7 , while the over-sea regime (large Ro 0 ) tends to involve a smaller range of Ro . We note that the GDL from Eq. (27) limits how large B can be; generally u * 0 /G < κ/B, so values of B greater than those shown are not physical. The model results in Fig. 6 do not violate this limit. Figure 7 is a reproduction of Hess and Garratt (2002), 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 Hess and Garratt (2002). For Ro L − = 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 Hess and Garratt (2002) 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 Ro L − = 0 and Ro = 10 3 -10 4 , which is the range used by Blackadar (1962) (Ro = 3.7 × 10 3 ). For Ro L − = 0, smaller values of the cross-isobar angle can be simulated compared with the analytic solution of Ellison (1956) due to the enhanced rate of mixing. The model cannot predict larger  Hess and Garratt (2002). Cross-isobar angle simulated by the limited-length-scale k-ε model extended to unstable surface layer stratification, taken at a normalized height of (z + z 0 )|f c |/G = 5 × 10 −5 for different Ro 0 , Ro , and Ro L − . Black markers represent measurements from Hess and Garratt (2002). Ro = 3.7 represents max from Blackadar (1962). Analytic results of Ekman (1905) and Ellison (1956) are summarized in Appendix A. values of the cross-isobar angle compared to the analytic solution of Ekman (1905) (45 • ). Peña et al. (2014) 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 1.21 × 10 −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. Peña et al. (2014) also calculated a "mixing" (turbulence) length scaleˆ using a local friction velocity u * and the wind speed gradient:

Atmospheric surface layer profiles
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 10 3 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. 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 = C 3/4 µ k 3/2 /ε. The original limited-length-scale kε model of Apsley and Castro (1997) 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ˆ , 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 dU/dz = u * / and = κz), and the optimizer used to fit G and max sets max to our chosen maximum value of 10 3 m. Note that therefore the lines corresponding Wind Energ. Sci., 5, 355-374, 2020 www.wind-energ-sci.net/5/355/2020/  . Unstable cases are also simulated with our extension to unstable surface layer stratification with L from Table 2. 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 * = (u w 2 +v w 2 ) 1/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 simulationto-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. Peña et al. (2014) 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. Peña et al. (2014) selected 10 cases that differ in geostrophic forcing and atmospheric Table 3. ABL validation cases based on Peña et al. (2014).

Atmospheric-boundary-layer profiles
Stable, strongly forced 4.5 × 10 −3 20.5 1.6 × 10 −2 0.45 1.0 × 10 7 -5 Neutral −5 × 10 −4 19.5 1.6 × 10 −2 0.70 1.0 × 10 7 -9 Very unstable, weak forcing −4.0 × 10 −2 5.02 1.6 × 10 −2 0.26 2.8 × 10 6 1.7 × 10 3 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 Peña et al. (2014): cases 4, 5, and 9 and the corresponding values of the Obukhov length, geostrophic wind, roughness length, friction velocity, Ro 0 , and Ro L − are listed in Table 3. For convenience, we keep the case names as introduced by Peña et al. (2014). Cases 4 and 5 represent a stable and a neutral ABL with high forcing, respectively, where Ro 0 = 10 7 . Case 9 is characterized by a low forcing and very unstable stratification, where Ro 0 = 2.8 × 10 6 . In Case 6 from Peña et al. (2014) 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 Peña et al. (2014) 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 (Hong et al., 2006) 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 Zilitinkevich and Esau (2005), the "mid-ABL" scale of Gryning et al. (2007) (generalized by Kelly and Troen (2016) for matching G), and the "top-down" scale of Kelly et al. (2019). 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 Ro 0 10 7 , 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 Apsley and Castro (1997) 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 Ro 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.
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 turbulenceresolving methods, such as large-eddy simulations, in future work.

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 mea-surements 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 (Ekman, 1905;Ellison, 1956) act as bounds on the results obtainable by the limitedlength-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 (van der Laan et al., 2017). 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 threedimensional RANS simulations of complex terrain and wind farms.
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.