Rossby number similarity of an atmospheric RANS model using limited-length-scale turbulence closures extended to unstable stratification
- DTU Wind Energy, Technical University of Denmark, Risø Campus, Frederiksborgvej 399, 4000 Roskilde, Denmark
Correspondence: Maarten Paul van der Laan (email@example.com)
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.
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 (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 two-equation 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 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 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 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.
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, 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 and are modeled by the linear stress–strain relationship of Boussinesq (1897): and , where νT is the eddy viscosity. The boundary conditions for U and V are 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:
where is a complex variable and .
The eddy viscosity, νT, needs to be modeled in order to close the system of equations. The eddy viscosity can be written as , 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≪zi), 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 z≪zi and a constant eddy viscosity for z∼zi should provide a more realistic solution. For example, the eddy viscosity could have the form , 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 necessary, and one of the simplest numerical models for the idealized ABL is given by Blackadar (1962) using Prandtl's mixing-length model:
where 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 𝒫 is the turbulence production, and σk, σε, Cε,1, and Cε,2 are model constants that should follow the relationship . 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 :
where the turbulence length scale is modeled as . 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 (), 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; Apsley and Castro (1997) noted that 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 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 limited-length-scale k–ε of Apsley and Castro (1997) in terms of wind speed, wind direction, , and eddy viscosity. The Ekman spiral is depicted with two constant eddy viscosities, which only translates the solution vertically. In addition, we have chosen s−1, G=10 m s−1, and 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 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≥zi.Ekman (1905)Ellison (1956)Blackadar (1962)Apsley and Castro (1997)
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,
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, , the original length-scale model of Blackadar (1962) is obtained. Note that in stable conditions, , 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
3.2 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 and β≈5. They proposed a more complicated expression for 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). 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 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 ℬ in the k–ε transport equations:
Here ℬ is modeled as
following MOST, using the similarity functions of Dyer (1974) as discussed in van der Laan et al. (2017). We use the flow-dependent parameter of Sogachev et al. (2012), which for unstable conditions includes the prescription
amenable to the free-convection limit: for . Further, αB→1 as ℓ→0, matching neutral conditions since z∕L also vanishes then. The prescription (Eq. 17) results in
which also means that approaching the effective ABL top (ℓ→ℓmax), so sources and sinks of ε balance in Eq. (15), i.e., , all have the same coefficient Cε,2.
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.
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, 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 Sk,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 () and the diffusion terms are zero (), 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:
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 and , 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.
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 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 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+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 Richards and Hoxey (1993), 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 Apsley and Castro (1997), using 48, 96, 192, 384, and 768 cells. We choose s−1, 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 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.
The numerical solution of the original limited-length-scale turbulence closures of Blackadar (1962) and Apsley and Castro (1997) 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 , and , where 𝒰 and ℒ are characteristic velocity and length scales, respectively:
Here, Ro is the Rossby number, , 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 using Eqs. (3) and (4), then Eq. (21) can be written as
where , 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:
Here, we have obtained Roℓ by rewriting the second dimensionless number ℒ∕ℓmax as the ratio of the two Rossby numbers: . Ro0 is known as the surface Rossby number, first introduced by Lettau (1959); it also resembles a ratio of (inertial) boundary layer depth to z0. Analogously, Roℓ is like the ratio of two boundary layer depths, and zi (e.g., Arya and Wyngaard, 1975); 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:
i.e., the ratio of simulated ABL depth to z0 is the lone parameter. Blackadar (1962) found a characteristic maximum ABL turbulence length scale of 0.00027 for the Leipzig 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. 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 ( 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 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 Ro0, while further away, it is mainly influenced by Roℓ with a weaker dependence on Ro0.
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 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 () denotes that Eq. (25) is defined for unstable conditions, i.e., L≤0. For the convective boundary layer, is generally replaced by the dimensionless inversion height , because the convective ABL depth does not have a significant dependence on (Arya, 1975). However, we note that 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 dictates the relative increase in mixing length (i.e., in the dimensionless coordinate ). Our length-scale-limited turbulence closures extended to unstable surface layer stratification, as presented in Sect. 3, are dependent on . 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 is a third nondimensional number, which can also be written as the ratio of two Rossby numbers: . For , 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 (0, 5×102, and 2×103) for Roℓ=103. For these Rossby number combinations, and correspond to near-unstable conditions (–0.005 m−1) and unstable to very unstable conditions (–0.02 m−1), respectively. Figure 4 shows that both extended turbulence closures only depend on Ro0 and 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 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 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 ). 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 . 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 () 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 combinations of the three Rossby numbers Ro0, Roℓ, and . We have chosen G=10 m s−1 and 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 . The normalized ABL depth in the RANS model increases for stronger unstable surface layer conditions (larger ), i.e., for larger values of the surface heat flux. For neutral and stable conditions () and moderate to shallow ABL depths, i.e., – corresponding to m as seen in Fig. 5 – we find that , 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 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 () and deviates from a power law, while for unstable conditions a is similar to the previously found value of 0.6.
We employ the Rossby similarity from Sect. 5 to validate a range of results simulated by the original limited-length-scale k–ε model of Apsley and Castro (1997) including our proposed extension to unstable surface layer stratification. Historical measurements of the geostrophic drag coefficient 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. (2010, 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.
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 Petersen, 1989), which connects the surface layer properties as z0 and u*0 with the driving forces on top of the ABL proportional to :
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 and Ro0:
Figure 6 is a reproduction from Hess and Garratt (2002), where the geostrophic drag coefficient is depicted as a function of surface Rossby number Ro0. 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-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 and , representing unstable ( m−1) and very unstable conditions ( m−1) for the chosen values of G=10 m s−1 and 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 , 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 , log (Ro0)} relationship becomes more linear. In addition, for , as used by Blackadar (1962), and , 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 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 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 , 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 –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 , so values of B greater than those shown are not physical. The model results in Fig. 6 do not violate this limit.
6.2 Cross-isobar angle
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 , 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 and Roℓ=103–104, which is the range used by Blackadar (1962) (). For , 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 values of the cross-isobar angle compared to the analytic solution of Ekman (1905) (45∘).
6.3 Atmospheric surface layer profiles
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 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:
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.
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 . 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 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 . 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.
6.4 Atmospheric-boundary-layer profiles
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 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, Ro0, and 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 Ro0=107. Case 9 is characterized by a low forcing and very unstable stratification, where .
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 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 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 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 turbulence-resolving methods, such as large-eddy simulations, in future work.
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 (Ekman, 1905; Ellison, 1956) 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 (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 three-dimensional RANS simulations of complex terrain and wind farms.
A1 Constant eddy viscosity – Ekman spiral
The analytic solution of Ekman (1905), known as the Ekman spiral, can be expressed as a function of a single variable, the normalized height . The wind speed and the wind direction θ can be written as
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
where x is a normalized height and c is a constant. For z→z0 (and assuming ), the Kelvin functions can be expanded, and the solution can be written as
where γe≈0.57721 is the Euler–Mascheroni constant. We can set the geostrophic wind G through the constant c:
by taking dU∕dz in Eq. (A2) for z→z0. 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 .
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 and , as also shown by Krishna (1980). The friction velocity in Eq. (A4) can now be calculated by solving the GDL for . Hence, the analytic solution of Ellison (1956) 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 and the Rossby number Ro0 using Eq. (A4):
where the GDL can be used to solve for .
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 from Fig. 6, since we can write the following:
The Rossby similarity based on u0* is illustrated in Fig. B2 for four combinations of (105 and 108) and (102 and 104) for , using four combinations of u*0 (0.2 and 0.4 m s−1) and fc ( 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 ).
It should be noted that u0* in Fig. 6 was extracted at a normalized height of , 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 , 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.
The numerical results are generated with proprietary software, although the data presented can be made available by contacting the corresponding author.
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.
The authors declare that they have no conflict of interest.
This paper was edited by Johan Meyers and reviewed by Javier Sanz Rodrigo and one anonymous referee.
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. 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
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
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
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., 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
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
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
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 .