Structuring of turbulence and its impact on basic features of Ekman boundary layers

The turbulent Ekman boundary layer (EBL) has been studied in a large number of theoretical, laboratory and modeling works since F. Nansen’s observations during the Norwegian Polar Expedition 1893–1896. Nevertheless, the proposed analytical models, analysis of the EBL instabilities, and turbulence-resolving numerical simulations are not fully con istent. In particular, the role of turbulence elforganization into longitudinal roll vortices in the EBL and its dependence on the meridional component of the Coriolis force r main unclear. A new set of large-eddy simulations (LES) are presented in this study. LES were performed for eight different latitudes (from 1 ◦ N to 90 N) in the domain spanning 144 km in the meridional direction. Geo trophic winds from the west and from the east were used to drive the development of EBL turbulence. The emergence and growth of longitudinal rolls in the EBL was simulated. The simulated rolls are in good agreement with EBL stability an lysis given in Dubos et al. (2008). The destruction of rolls in the westerly flow at low latitude was observed in simulations, which agrees well with the action of secondary instability on the rolls in the EBL. This study quantifies the effect of the meridional component of the Coriolis force and the effect of rolls in the EBL on the internal EBL parameters such as friction velocity, cross-isobaric angle, param ters f the EBL depth and resistance laws. A large impact of the roll development or destruction is found. The depth of the EBL in the westerly flow is about five times less than it is in the easterly flow at low latitudes. The EBL parameters, which depend on the depth, also exhibit large difference in these two types of the EBL. Thus, this study supports the need to include the horizontal component of the Coriolis force into theoretical constructions and parameterizations of the boundary layer in models.


Introduction
The earth's rotation is an important factor affecting atmospheric and oceanic dynamics.It is known that large-scale dynamics mostly depend on f V = 2 | | sin ϕ.This Coriolis parameter, f V , is a projection of the earth's angular velocity on the direction normal to the earth's surface at latitude ϕ.The projection, f H = 2 | | cos ϕ, tangential to the surface is usually omitted in models of large-scale atmosphere and ocean dynamics.Following this tradition, f H is usually omitted in theoretical analysis (Grisogono, 1995;Tan, 2001) and turbulence parameterizations as well (e.g., Holt and Raman, 1988;Andren and Moeng, 1993;Ayotte et al., 1996;Zilitinkevich and Esau, 2005).Historically, this omission was motivated by the focus on the turbulent planetary boundary layer dynamics in high latitudes (e.g., Nansen, 1900;Ekman, 1905) where f H → 0. Later, some authors attempted to introduce f H -dependence into turbulence parameterizations (e.g., Garwood et al., 1985;McWilliams and Huckle, 1997).These modifications were not widely implemented due to algebraic difficulties created by the Coriolis terms and small amplitude of the effect estimated in the analysis of the Reynolds stress equations (Galperin et al., 1989) and in largeeddy simulations (LES) (Wang et al., 1996).
Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
Recent advances in analysis of the boundary layer instabilities (Dubos et al., 2008, hereafter referred to as D08) and more accurate LES (Zikanov et al., 2003, hereafter referred to as Z03; Salon and Armenio, 2011;Glazunov, 2010;Fricke, 2011;Esau, 2012) has reopened the problem for discussion.It was found that the primary instability caused by f V and the secondary instability caused by f H considerably restructure the turbulence in the boundary layer, resulting in development of longitudinal coherent vortices (rolls).The rolls accumulate significant turbulent kinetic energy.They modify cross-isobaric angle, surface friction velocity, turbulence anisotropy and eddy viscosity significantly more than it has been estimated in previous works.Although there are a number of papers presenting analytical and numerical studies of the rolls, their effect on the parameters of turbulent boundary layer parameterizations has not been quantified.Moreover, the majority of published studies performed simulations in too small computational domain and focused on a stationary state of the roll development, which is usually unachievable in the observed boundary layers.
This work presents LES of a non-stratified turbulent boundary layer flow over an infinite homogeneous, aerodynamically rough, flat surface in a rotated frame of reference.This idealized flow is frequently referred to as the Ekman boundary layer (EBL).The analytical solution of the stationary EBL problem with constant viscosity is known as the Ekman velocity spiral (Ekman, 1905).Qualitatively similar solutions were obtained for the stationary EBL with a more realistic height-dependent viscosity profile, ν = ν(z) (Grisogono, 1995;Tan et al., 2000).Quantitatively, this dependence, ν(z) → 0, reduces the cross-isobaric angle α (α = 45 • in the Ekman solution) between the geostrophic wind in the free flow above the turbulent EBL and the surface layer wind (Svensson and Holtslag, 2009), which agrees well with observations in the atmosphere (Zhang et al., 2003) and the ocean (Stacey et al., 1986;Wijffels et al., 1994;Chereshkin, 1995;Price and Sundermeyer, 1998), laboratory experiments (Caldwell et al., 1972;Howroyd and Slawson, 1974;Nickels and Joubert, 2000) and numerical simulations (Andren et al., 1994;Coleman, 1999;Esau, 2004a;Marlatt et al., 2012).
Here and in the following text, we refer only to a small subset of the enormous literature on the EBL, which, however, mostly neglects possible f H effects.
The stationary EBL solution is a useful reference point for further discussion.However, this solution has been found unstable with respect to turbulent perturbations (Faller, 1963;Lilly, 1966;Brown, 1970;Leibovich and Lele, 1985; see also reviews in Etling and Brown, 1993;Z03;D08).Schematically, EBL instabilities can be summarized as follows.The laminar EBL flow becomes unstable at effective Reynolds numbers, Re = Gδ ν > 50 (Faller, 1963;Lilly, 1966), where geostrophic wind speed, G; an EBL depth scale, δ; and effective eddy viscosity, ν; will be defined later.The inflection-point instability (due to the inflection point on the Ekman velocity spiral) emerges and dominates at Re > 100 (Brown, 1970).Although this instability is due to f V , its growth rate and scaling are sensitive to f H (Leibovich and Lele, 1985).The selective growth of eddies of certain length scales results in development of a secondary flow that is composed of rolls.At sufficiently large Re > 300 (D08), the rolls themselves experience secondary instability, which is also sensitive to f H and to the direction of the mean flow.The secondary instability either stabilizes or destabilizes these rolls.These complex dynamics makes EBL properties sensitive to latitude and wind direction (Z03; D08).Moreover, these dynamics and sensitivity are not represented in the Reynolds stress equations (Galperin et al., 1989) as they are highly selective with respect to the length scale and structure of the turbulent perturbations.The impact of the EBL instabilities on the integral EBL parameters remained unclear.The observed EBL is never stationary, non-stratified or barotropic.Hence, it is difficult to judge whether differences between the observed EBL and the EBL solution result from dynamical instabilities or from external factors.
Although the complete description of the planetary rotation is typically included in the LES equations, the majority of LES studies consider only a single latitude (45 • N or 90 • N) and wind direction (westerlies), as exemplified in the Andren et al. (1994) intercomparison study.The same is true for the laboratory experiments (e.g., Nickels and Joubert, 2000) and direct numerical simulations (e.g., Miyashita et al., 2006), which consider the spanwise (corresponding to ±90 • N) rotation of the domain.Much less of the literature considers the effect of f H change on the EBL.Moreover, published LES were run in small computational domains, which accommodate just a few rolls.As the spectrum of the primary EBL instability is continues (D08), one can expect emergence of rolls squeezed into a small domain, which is not proper for studies of the most unstable structures in the EBL.Analysis of the LES results, in terms of internal EBL parameters, has not been published yet.The LES runs in this study resolve these problems.They were conducted in a large domain of 144 km in the cross-flow direction, which accommodated 30 to 100 roll pairs.The analysis is done in terms of the internal EBL parameters proposed from theoretical models and parameterizations.Moreover, this analysis is presented not only for the stationary stage of the EBL development but also for the stage of the maximum difference between EBLs, which is reached by the sixth hour of the model integration.The next section formulates the problem and defines external and internal parameters to be studied.Section 3 describes numerical experiments and results of simulations.Section 4 discusses the results in the context of theoretical studies and parameters of the EBL parameterizations.Finally, Sect. 5 highlights conclusions.The brief description of the LES code LESNIC is given in the Appendix.

External control parameters
The Ekman boundary layer is controlled by three external parameters, namely the angular velocity of the frame rotation, [s −1 ]; surface roughness, z 0 [m]; and the horizontal pressure gradient, conveniently expressed through a geostrophic velocity of the free flow, is the static pressure and ϕ is the latitude of the frame of reference.The local coordinate system is given by the orthogonal axes with the unit vectors: x -directed toward the east; y -directed toward the north; and z -directed normal to the surface opposite to the gravity acceleration vector g.The geostrophic wind speed is defined as G = U g .The angle between the x axis and the geostrophic velocity is β. Figure 1 shows the coordinate system, the main vectors and the simulation domain.

Equations and internal EBL parameters
The EBL dynamics is governed by the incompressible Navier-Stokes equations (2) Here, u [m s −1 ] is the flow velocity and ν [m 2 s −1 ] is the kinematic viscosity, which in the turbulent EBL is assumed to be an eddy viscosity.In the case of the horizontally homogeneous steady-state flow, these equations admit the wellknown Ekman stationary solution (e.g., D08), This solution gives the Ekman depth scale δ = 2 √ 2ν/f V [m] as an internal length scale parameter.In the turbulent EBL above the viscous sublayer, scaling on the molecular kinematic viscosity becomes irrelevant.In order to circumvent this problem, ν should be understood as an effective viscosity of the turbulent fluid (Muschinski, 1996;Z03;Esau, 2004a;D08), which is an integral parameter of the EBL only indirectly linked to the variable eddy viscosity in the turbulent flow.
The effective viscosity can be defined through the fluxgradient relationship in the following form The effective viscosity becomes The Ekman depth scale will be expressed as and will give a simple expression for an effective Reynolds number The last equation relates the EBL bulk Reynolds number with a single internal turbulent parameter, namely with the friction velocity u * = |τ (z = 0)| The above definitions are consistent with the behavior of the turbulence in the present LES.The reader should pay attention to the obtained reciprocal dependence Re ∝ u −2 * .This dependence is known for the low Re boundary layer flows where Re is based on the constant molecular kinematic viscosity.At molecular Re > 10 3 , the boundary layer undergoes transition to well-developed turbulence and the dependence becomes the direct proportionality Re ∝ u * .Thus, the introduced effective Re characterizes the EBL as a weakly turbulent boundary layer with respect to the most energetic eddies in the flow.As LES in this study show, the simulated EBL do not undergo transition to the fully chaotic state for effective Re < 10 3 , but such a transition does occur for effective Re > 10 3 .Thus, this effective Re correctly describes the EBL flow behavior.Dynamical reasons for such a behavior remain poorly understood; specifically, the relative role of the limited grid resolution versus energy concentration in large eddies is unclear.

The logarithmic law of the wall and resistance laws
The layer just above the surface where the turbulence is the most intensive is traditionally described as a constant flux layer in theoretical constructions and parameterizations.In this layer, which typically occupies 15 to 20 % of the turbulent EBL, the non-dimensional velocity gradient is where k is the von Karman constant.The non-dimensional gradient in the non-rotated EBL is linearly scaled with height, giving = 1.The EBL cannot be considered as a non-stratified layer due to the effect of the Coriolis stratification as described in Bradshaw (1969) and subsequent publications (e.g., Esau, 2003).Galperin et al. (1989) introduced a stability correction compatible with the second order closure, which reads as where This parameterization is not fully correct as it states that f V always leads to > 1, which is incorrect for the EBL in the Southern Hemisphere and would disagree with laboratory experiments (Nickels and Joubert, 2000).However, since |C V | |C H | for zonal winds, where C V = −1.08 and C H = 3.51 cos β, and the presented LES runs were conducted in the Northern Hemisphere, this problem does not affect this study.Equation ( 9) shows that WE > 1 in the EBL driven by the flow from west to east (WE-EBL) and EW < 1 in the EBL driven by the flow from east to west (EW-EBL).
The difference | WE − EW | increases toward the equator, in proportion to cos ϕ.If the surface friction velocity is unchanged, which is not quite true as we will see later, the difference at the equator could be estimated as Here, the following values are used: At the pole the deviation from the expectation = 1 is too small to be firmly established from LES or observation data.Integration of Eq. ( 8) assuming = 1 gives the integral form of the law of the wall Then the resistance laws could be written as where Ro = G (f V z 0 ) is the roughness Rossby number (Zilitinkevich and Esau, 2005;Esau and Zilitinkevich, 2006).
The non-dimensional parameter estimates the depth of the essentially turbulent EBL (Rossby and Montgomery, 1935;Zilitinkevich et al., 2007).Empirical evidence shows that h PBL δ, which is also true for these LES.
Published studies have considered the EBL structure at the pole, as in Ekman (1905).Therefore, changes in the internal EBL parameters A, B and C R at lower latitudes were overlooked, as they were expected to be constants with the same values as at ϕ = 90 • N. As we will see here, these parameters depend on both ϕ and β, which suggests a dependence on f H .

Instability of the stationary Ekman solution
The stationary Ekman solution is known to be unstable with respect to turbulent perturbations.At Re > 100, the dominant primary instability is the inflection-point instability that results in the development of a secondary flow consisting of counter-rotating longitudinal rolls.The inflection-point instability is linked with f V , while f H modifies the stability of the rolls.This set of instabilities could affect the EBL parameters in a way unaccounted for by the Reynolds stress equations.The impact is related to selective growth of certain turbulent scales, which is not reflected in the corresponding changes of the turbulent dissipation.The cross-flow length scale of the fastest growing mode of the rolls is l roll = 4πδ ≈ π h EBL with a non-dimensional growth rate σ 1 = 0.024 at Re > 500.D08 mentioned that this growth rate should be considered as a function of the flow Reynolds number and the rotational stability parameters; that is, σ 1 = σ 1 (Re, ξ V , ξ H ). Such a theoretical analysis has been done for some combinations of parameters, but the general picture is still missing.The analysis in D08 shows that the rolls grow during the initial stage for t = 10 δ G σ 1 ∼ 6 h.This stage is followed by the saturation stage.At Re > 1000, the rolls do not saturate and, as we will show, they are destroyed by the secondary instability.The saturated roll stage was of primary interest in the earlier LES studies; for example, Z03 investigated the EBL structure after 416 h of the model run time.

Concept of numerical experiments
This study addresses the impact of the EBL structures during the growing and stationary stages of their development.Therefore, it was important to make the simulation domain sufficiently large in the cross-flow direction for undisturbed development of a statistically significant number of rolls.The longitudinal dimension along the x axis was designed to be only eight grid points (0.28125 km), whereas the cross-flow dimension along the y axis was 4096 grid points (144.0 km).The small size of the longitudinal dimension is sufficient to resolve the fastest growing mode of the expected secondary instability (D08).It is also sufficient for the dynamic-mixed turbulence closure in the model to work.The vertical dimension along z axis was 128 grid points (3.0 km).The height of the domain is sufficient to accommodate the stationary EBL at latitudes poleward of 30 • N. At lower latitudes, the EW-EBL comprised the whole domain by the sixth hour of model integration and therefore its further evolution was affected by this constraint.
The LESNIC runs are listed in Table 1.The runs were driven with the westerly (β = 0 • , therefore denoted as WE-EBL in Table 1) and easterly (β = 180 • , EW-EBL in Table 1) geostrophic winds with the speed G = 5 m s −1 .The runs were conducted in the "f plane" at eight different latitudes from 1 to 90 • N. Except for the geostrophic wind direction and the run's latitude, the other parameters were kept the same for all experiments.Surface roughness was set to z 0 = 0.1 m.The grid resolution was 35.15 m in the horizontal dimensions and 23.43 m in the vertical dimension.
Each run was integrated over 22 model hours with sampling of the horizontally averaged turbulence and mean flow statistics every 60 s and storage of the full three-dimensional fields of the model velocity at the end of each model hour.
The resulting simulation database provides information for the study of the time evolution of the EBL structures and parameters as well as for the statistical study of the EBL structures in their development.
In the following presentation of the LES results and in the discussion, we will consider the parameter change during the time evolution of the EBL as well as the EBL structure at the end of hours 3, 6 and 18 of the model integration.The results for hour six characterize the EBL structure and parameters approximately at the end of the roll growing stage in EW-EBL.

Turbulent surface drag and effective Reynolds number
We begin the discussion of the LES results with consideration of the key internal parameter -the turbulent surface drag or the surface friction velocity u * .Figure 2 reveals that in terms of Re both EW-EBL and WE-EBL reach the stationary stage at high latitudes where Re ≈ 500 (as expected for the environmental EBL) and u * ≈ 0.22 m s −1 .The WE-EBL at lower latitudes behaves differently.In these runs, Re increases for a longer period, and reaches much larger values.The Reynolds number is approximately 1000 at 5 • N.This corresponds to roughly 30 % of the turbulent stress reduction as compared to the EW-EBL at the same latitude.Smaller u * and C g indicate reduced vertical momentum transport.This difference is larger than obtained in the Reynolds stress analysis of Galperin et al. (1989), but much smaller than found in the LES of tidal flow by Salon and Armenio (2011).

The depth of the turbulent EBL
We define the depth of the essentially turbulent EBL, h PBL , as the height where the magnitude of the turbulent stress drops to 5 % of its surface value; that is, h PBL : Similarly, one can define the depth of the weakly turbulent layer as h 1 PBL : * and the depth of the strongly turbulent layer as h 10 PBL : τ (z = h 10 PBL ) = 0.10u 2 * .These three thresholds are shown in Figs.8-11 in this study.Previous studies have overlooked that h PBL and particularly the depth of the weakly turbulent layer, h 1 PBL , should be sensitive to the rolls, which are found in the upper part of the EBL. Figure 3 shows that the turbulent EW-EBL becomes deeper than the WE-EBL by 60 • N. Unfortunately, the LES domain was not sufficiently deep to accommodate the stationary EW-EBL at latitudes lower than 30 • N. Nevertheless, the faster growth of the EW-EBL depth is clearly seen in the early stages of the EBL development.
The development of the rolls in the EW-EBL has a profound impact on the EBL depth with a much deeper turbulent layer at low latitudes.This increase (by almost one order of magnitude) is not fully reflected in the surface stress and in the effective eddy viscosity, and therefore, in the Ekman scale and in the Rossby-Montgomery parameter, C R .Zilitinkevich et al. ( 2007) reviewed estimations of C R in published LES and DNS studies, as well as our own simulations with LESNIC to suggest the best fit of C R = 0.65±0.05.The new simulations, where C R was somewhat lower than in the fully three-dimensional simulations, show that this fit is valid only for latitudes poleward of 45 • N. Figure 3c shows the parameter C R1 , which characterizes the depth of the weakly turbulent layer, as it more correctly reflects the impact of the rolls on the EBL.Values of C R and C R1 decrease at low latitudes (Table 1), being just around 0.1 at 5 • N, highlighting that L V scaling becomes inadequate at these latitudes.

The cross-isobaric angle
The cross-isobaric angle defines several important phenomena (Svensson and Holtslag, 2009).Since this angle is determined by the balance between the turbulent friction and the pressure gradient, its value is sensitive to the depth of the EBL and to the vertical distribution of the eddy viscosity.In the turbulent EBL, α ≈ 19 • at the pole (Coleman, 1999;Esau, 2004a;Marlatt et al., 2012), which is defined from f V in the Ekman solution in Eq. ( 3).This angle increases in shallower layers and decreases in deeper layers (Esau and Zilitinkevich, 2006).Figure 4 shows that the deeper EW-EBL have α reduced to about 5 • at 15 • N, whereas the shallower WE-EBL still have α ≈ 20 • at this latitude.It directly follows from Eq. ( 3) that α should decrease toward the equator, which is observed in the LES data.

The resistance laws
The observed changes in the surface friction velocity, the cross-isobaric angle and the EBL depth do not compensate each other in the resistance laws.Hence, they should modify the parameters A and B of the resistance laws in Eqs. ( 15) and ( 16).When Zilitinkevich and Esau (2005) revisited the resistance laws, they focused on the effect of the stable stratification but overlooked the effects of the Coriolis stratification.This gave a formulation of the neutral limit of the resistance laws through variables where the constants C f A = C f B = 1.The proposed analytical resistance law approximations read where a = 1.4, a 0 = 1.65, b = 10 and b 0 = −2.The data scatter around these analytical expressions was quite large in the original work, and could be tolerated only when even larger systematic changes with stability, which strongly reduces C R , were considered.This scatter obscured the smaller systematic dependences of the coefficients upon latitude and wind direction.Figure 5 shows these dependences.The parameter A characterizes the resistance to the zonal, and in our case longitudinal, velocity components.It changes relatively little as it is mainly sensitive to the vertical momentum transport.The parameter B characterizes the resistance to the meridional, and in our case cross-flow, velocity components.The differences in it are sensitive to the development of the secondary flow, i.e., to the rolls.Therefore, changes in B are large, suggesting that further revision of the resistance laws is required.

Non-dimensional velocity gradient in the constant flux layer
Calculation of the non-dimensional velocity gradients, , in the LES meets certain difficulties given that the von Karman constant, k = 0.4, as prescribed at the first model level, does not necessarily have a value consistent with the simulated turbulence.We redirect the reader to Österlund et al. (2000) for a more detailed discussion on the von Karman constant values.Therefore, as the first step, the consistent value of k should be recovered.This is done using Eq. ( 14) with averaging of the obtained k over the lowest 15 % of the essentially turbulent EBL.These values of the von Karman constant are shown in Fig. 6a.The non-dimensional velocity gradients were calculated at each level using this constant k (Fig. 6b), and then averaged over the same fraction of the EBL.This procedure has proved to be robust for identification of the logarithmic regions in the velocity profiles, and is often used in analysis of laboratory experiments ( Österlund et al., 2000;Nickels and Joubert, 2000).
The non-dimensional velocity gradient ( = 1.03 in LES) corresponds well with our ad hoc estimation ( = 1.02 in Eq. 13) for the pole.Similarly, the difference | WE − EW | = 0.4-0.8 at low latitudes is in good agreement with the estimation of 0.35 from Eq. ( 12).This agreement is surprising as changes in due to f H effects and the EBL structuring must be consistent with changes in the other internal PBL parameters.The presented LES show that those changes are significant, and even large (e.g., in h PBL ), whereas Galperin et al. (1989) and subsequent publications invoking the analysis of the Reynolds stress model maintain that the effect should be small.Nevertheless, as for , the LES results and the model analysis seem to be in agreement.

The impact of rolls on the EBL parameters
A question may be asked about the impact of rolls -i.e., the large-scale, self-organized structures in the EBL -on the observed dependences of the internal EBL parameters on the wind direction and latitude.On the one hand, the Reynolds stress equations, such as in the second-order closure by Galperin et al. (1989), do not account for instability and selective growth of turbulent perturbations with a certain structure.Nevertheless, some parameters, e.g., the nondimensional velocity gradient, derived from those equations, are in good correspondence to the LES results.On the other hand, the EBL stability analysis in D08, which attributes the dependences to restructuring of the EBL, is also consistent with the LES results, including the instant flow visualizations in terms of the anomalous velocity and stream function.Esau (2004b) attempted to estimate this effect by applying constraints on the EBL depth, which would constrain the size of the largest eddies in the boundary layer.The EBL depth was constrained through imposed stable stratification.Significant effect was found but the size of the LES domain was not sufficiently large to firmly establish the effect.D08 quoted at least two possible dynamical mechanisms for the Coriolis force to impact the EBL parameters.It may act indirectly through modification of the Reynolds stress component reserving some fraction of the turbulent kinetic energy E = 0.5 u 2 + v 2 + w 2 in the cross-flow velocity fluctuations v 2 , which are otherwise much less intense.The energy in the longitudinal velocity fluctuations u 2 is quickly replenished in the sheared flow.Hence, this mechanism is important in the strongly-sheared constant-flux layer, where LES and Reynolds stress equation analysis are in good correspondence.The horizontal Coriolis force component, f H , also modifies the domain of the primary instability of the flow directly (Leibovich and Lele, 1985), but this effect should be small at Re ∼ 500, as found in LES.Instead, f H acts through the secondary instability, which stabilizes or destabilizes the rolls depending on the mean flow direction and latitude.Since the rolls are dominant structural features in the weakly turbulent layer in the upper part of the EBL, their stabilization or destabilization should have a major impact on the EBL depth and related parameters.
Figure 7 shows quasi-regular alternations of the instant integrated vertical velocity component at hour six of the simulations.The angular brackets The stream function is defined following D08 as Except for the WE-EBL at low latitudes, the closed circulation contours are clearly observed in these figures.These contours identify the rolls in the EBL with both the most intensive roll circulations (the largest amplitude | | of the stream function) and the largest spatial scales of the rolls, found in the EW-EBL at low latitudes.Centers of the roll circulation are located in the upper part of the EBL at z δ ∼ 3-4, which is consistent with theoretical analysis.
The closed stream function contours are located in the weakly turbulent layer, as comparison with the momentum flux profile confirms.Further analysis found high correlation (corr.coef.0.6) between | | y and the EBL depth.Thus, the rolls are essentially a feature of the upper EBL.The observed roll length scale is l roll = 4π δ ≈ π h PBL ≈ 3 km.We conducted four additional LES runs in a horizontally restrictive domain, L y = 0.5 h PBL (ϕ = 5 • N) = 1225 m (Table 1), to prevent development of the rolls at the most unstable mode.It was expected that without the rolls, the WE-EBL and EW-EBL at low latitudes would exhibit less different values of the internal parameters.This is, however, not the case.Figure 12 shows stream functions and velocity anomalies of four LES runs.As Fig. 12c and d show, the rolls constrained by the LES domain are organized in a twolevel structure where the size of rolls at each level fits the size of the domain.This explains the success of earlier simulations with too restrictive domains (Mason and Sykes, 1980;Z03).This result follows from the continuity of the spectrum of unstable modes in the EBL.If the most unstable mode does not fit into the LES domain, a circulation pattern consistent with the largest allowed unstable mode develops.By the other words, the instability amplifies the largest mode resolved in the domain.Thus, we are not able to resolve ambiguity around the role of the rolls in the EBL.However, these experiments suggest that the rolls have an impact on the EBL indirectly through modification of the EBL depth and subsequent change of parameters, which depend on the EBL depth.

The secondary instability of rolls
At lower latitudes, the characteristics of the WE-EBL and the EW-EBL become increasingly different, suggesting that the control of f H should be sensitive to the wind direction β.There are two mechanisms of the f H control described in the literature (D08).One mechanism (Leibovich and Lele, 1985) is supposed to be active at Re < 300.It appears due to different growth rates of the rolls in the WE-EBL (small growth rate this mechanism.The present LES demonstrated that both the EBL depth and the turbulence length scales are consistently decreasing in the WE-EBL from the pole, where this mechanism is inactive, toward the equator, where this mechanism is supposed to be responsible for the structuring of the flow.The weak but still growing rolls are not observed in the simulations. Another mechanism is secondary instability of the rolls, which becomes increasingly active at Re > 300.This instability perturbs the rolls in the longitudinal direction.But its growth rate, σ 2 , over the range of length scales permitted in the LES domain, depends on β.D08 analysis for the maximum permitted LES length scale in the zonal x direction (the minimum wave number k 2 = 2.7, latitude 45 • N and Re = 500) gives the non-dimensional growth rates equal to σ WE 2 = +0.005and σ EW 2 = −0.005.It suggests roll stabilization in the EW-EBL and destabilization, albeit only weakly, in the WE-EBL.The simulated EBL evolution suggests that there might be a run-away loop of roll destabilization in the WE-EBL.The initially weak destabilization reduces the roll circulation and decreases the momentum transport in the EBL.This raises the effective Re ∝ 1 u 2 * of the flow, making the rolls more unstable with respect to the secondary instability.These more unstable modes further destabilize the rolls, raise Re and then destabilize the rolls even more.The S-shaped colored structures in the velocity anomalies between the weak rolls are observable in Fig. 8c.These structures are similar to the velocity structures responsible for the secondary instability in D08.
The destructive secondary instability counteracts the primary instability in the EBL.The later one is, however, stronger; that is, their growth rates relate as σ 1 >σ 2 .Hence, weak and strongly perturbed rolls emerge slowly even in the WE-EBL at latitudes larger than 5 • N. Detailed examination of the LES results at 1 • N reveals a qualitatively different scenario.In close vicinity to the equator, the rolls do not emerge, leaving the WE-EBL in a chaotic state.This is probably the effect of suppression of redistribution between longitudinal and vertical components in the turbulent stress as it follows from the Reynolds stress equations.This effect has been noted in Z03 and appeared in these LES as well.

Conclusions
The new set of large-eddy simulations in this study confirmed self-organization of the EBL turbulence into a set of longitudinal vortices (rolls).The simulated rolls are in many aspects similar to those roll structures that have been previously found in analytical studies of the EBL instabilities.If our definition of the effective Re of the EBL is accepted, then not only the qualitative structure of the rolls but also their quantitative characteristics (such as the cross-flow length scale, the growth rate and sensitivity to the secondary instability), are in good correspondence with the proposed EBL theory.The rolls are found in both the WE-EBL and the EW-EBL at high latitudes.The structures of both boundary layers are nearly identical at high latitudes, but the difference between them gradually increases toward the equator.At low latitudes, the rolls in the EW-EBL quickly grow, and by hour six of simulations occupy the entire LES domain.They maintain the strong vertical momentum transport and surface stress, keeping the effective Re at about 500.The rolls in the WE-EBL evolve slowly, experiencing destabilization from the secondary instability and suppression of the momentum redistribution by the Coriolis force.The difference between the EW-EBL and the WE-EBL was found to be the largest by hour six of simulations.
The previously published analyses of the Reynolds stress equations (e.g., Galperin et al., 1989) concluded that the impact of the horizontal component of the Coriolis force, f H , on the internal EBL parameters is rather small.The theoretical studies of the EBL structures (e.g., Etling and Brown, 1983;D08) have not come to a certain conclusion on this point.These LES revealed that the impact of f H and the corresponding differences in the EBL structure is significant.The LES results are consistent with theoretical estimations in the constant flux layer.At the same time, the LES show large changes in the EBL depth and the structure of the upper EBL due to the effects of f H .At low latitudes (below 30 • N), the turbulent EW-EBL is 3-5 times deeper than the WE-EBL.Correspondingly, the internal EBL parameters, which are sensitive to the EBL depth (the Rossby-Montgomery parameter, the cross-isobaric angle, the resistance law constants), also show sensitivity to the wind direction and latitude, confirming the results in Esau (2012).
Unfortunately, this set of LES runs does not allow for unambiguous conclusions about the impact of the rolls as compared to the impact of the Coriolis terms in the Reynolds stress equations.For a large number of applications, this distinction will be irrelevant as they can rely on a simplified parameterization of the phenomenon.Some applications, e.g., the cloud-radiation processes and dispersion of air pollution, require explicit knowledge of the EBL structure and the roll length scales.A more accurate inclusion of the effects discussed here into model parameterizations remains a challenge for future studies.

Fig. 1 .
Fig. 1.The coordinate system used in this study with definitions of the main vectors and angles.The orientation of the LES computational domain is sketched by the shaded rectangle.x axis points due east, y axis points due north and z axis is normal to the planet surface at given latitude ϕ.

Fig. 2 .
Fig. 2. The geostrophic drag coefficient C g = u * /G: (a) time evolution of the effective Re = C −2 g (solid curves are for EW-EBL; dashed curves are for WE-EBL); (b) values of the coefficient at hour six of the LES integration.Shaded circles show the WE-EBL; open circles, the EW-EBL.

Fig. 3 .
Fig. 3.The depth of the essentially turbulent EBL: (a) time evolution of the non-dimensional depth h + EBL = h EBL L z , where L z = 3 km is the vertical size of the LES domain (solid curves are for EW-EBL; dashed curves are for WE-EBL); (b) values of h + EBL at hour six of the LES integration; and (c) values of the Rossby-Montgomery parameter C R1 for the weakly turbulent EBL at hour six of the LES integration.Shaded circles show the WE-EBL; open circles, the EW-EBL.

Fig. 4 .
Fig. 4. The cross-isobaric angle in the EBL: (a) time evolution (solid curves are for EW-EBL; dashed curves are for WE-EBL); (b) values of α at hour six of the LES integration.Shaded circles show the WE-EBL; open circles, the EW-EBL.

Fig. 5 .
Fig. 5.The parameters A (a) and B (b) of the resistance laws in the EBL.The solid curve is the analytical approximations with h PBL and α for the WE-EBL (dashed line) and the EW-EBL (solid line).Shaded circles show the WE-EBL; open circles, the EW-EBL.The data correspond to hour six of integration.

Fig. 6 .
Fig. 6.The non-dimensional velocity profile: (a) the von Karman constant k in LES averaged over the lowest 15 % of the EBL; (b) computed with this k and averaged over the same fraction of the EBL.Shaded circles show the WE-EBL; open circles, the EW-EBL.The data correspond to hour six of integration.

Fig. 7 .
Fig. 7.The instant cross-flow fluctuations of the vertical velocity component integrated along z and x axes at the end of the linear growth stage at the hour six of the LES run.The fluctuations are shown for latitudes 90 • N, 30 • N and 1 • N in the EW-EBL (a) and WE-EBL (b).The values are normalized by the maximum value at each latitude.The positive values correspond to upward motions.

Fig. 8 .
Fig.8.Cross-flow sections (yz planes) of the instant flow patterns in WE-EBL at 5 • N. The patterns are given at the beginning (after 3 h of LES integration at the upper panel) and the end (after 6 h at the middle panel) of the linear growth stage as well as for the quasi-stationary stage (after 18 h at the bottom panel).Contours represent the anomalies flow stream function from Eq. (30), with negative values shown as dashed contours.The color shading represents the anomalous flow velocity from Eq. (29), with blue corresponding to the slow motions.The time of the snapshots, the EBL Re at that time and the range of the stream function values are indicated on the panels.The horizontal lines show the EBL depths h 10 PBL (the lowest line), h PBL and h 1 PBL (the highest line).

Fig. 12 .
Fig. 12. Cross-flow sections (yz planes) of the instant flow patterns in four LES runs: (a) WE U5L5 and (c) WE U5L90 in the small LES domain; (b) WE U5L90 and (d) WE U5L90 in the large LES domain.The panels are composed similarly to Fig. 8.