## Abstract

The vertical turbulent velocity variance normalized by the convective velocity squared as a function of the boundary layer depth–normalized height [i.e., ] in the convective boundary layer (CBL) over a homogeneous surface exhibits a near-universal profile, as demonstrated by field observations, laboratory experiments, and numerical simulations. The profile holds over a wide CBL stability range set by the friction velocity, CBL depth, and surface heating. In contrast, the normalized horizontal turbulent velocity variance increases monotonically with decreasing stability. This study investigates the independence of the profile to changes in CBL stability, or more precisely, wind shear. Large-eddy simulations of several convective and neutral cases are performed by varying surface heating and geostrophic winds. Analysis of the turbulent kinetic energy budgets reveals that the conversion term between and depends almost entirely on buoyancy. This explains why does not vary with shear, which is a source to only. Further analysis through rotational and divergent decomposition suggests that the near-universal profile of is fundamentally related to the dynamics and interactions of local and nonlocal CBL turbulence. Specifically, the preferential interactions between local wavenumbers and the downscale energy cascade of CBL turbulence offer plausible explanations to the universal profile of .

## 1. Introduction

Lenschow et al. (1980, hereafter L80) presented the profile of vertical velocity variance in the convective boundary layer (CBL) based on empirical fitting to aircraft observations made during the Air Mass Transformation Experiment (AMTEX) over the East China Sea (Lenschow and Agee 1976),

where *w* (m s^{−1}) is the perturbation vertical velocity, (m s^{−1}) is the convective velocity scale, (m) is the boundary layer depth, *g* (m s^{−2}) is the gravitational acceleration constant, (K) is the reference potential temperature, and (K m s^{−1}) is the surface sensible heat flux. The overline is an ensemble-averaging operator.

Since then, Eq. (1) has been tested and verified in independent field campaigns. For example, profiles of taken from the Phoenix 78 experiment, as well as four earlier field experiments in the 1970s, were well described by Eq. (1) (Young 1988). Lenschow et al. (2012, hereafter L12) analyzed profiles collected during the Lidars in Flat Terrain (LIFT) experiment (Cohn et al. 1998), and showed that they also observe the empirical scaling of L80. Berg et al. (2017) presented yearlong vertical velocity statistics derived from Doppler lidar data for the cloud-free CBL over the U.S. Southern Great Plains. The composite profile based on several thousand individual 1-h periods again has a similar shape as Eq. (1), although its magnitude is slightly smaller. Lareau et al. (2018) extended Berg et al.’s work to include shallow cumulus days at the same site. They showed that shallow cumulus CBL variance profiles also agree with L80’s scaling almost irrespective of the cloud fraction.^{1} A brief summary of all the abovementioned field experiments is presented in Table 1, and the list is by no means exhaustive.

Besides field experiments, L80’s profile has also been verified in laboratory experiments (Willis and Deardorff 1974; Kumar and Adrian 1986; Adrian et al. 1986; Hibberd and Sawford 1994; Fedorovich et al. 1996), and numerical studies using large-eddy simulation (LES) and direct numerical simulation of dry CBL cases (e.g., Deardorff 1974; Moeng 1984; Schmidt and Schumann 1989; Sullivan and Patton 2011; Garcia and Mellado 2014; Salesky et al. 2017), as well as shallow cumulus and stratocumulus cases (e.g., Deardorff 1980; Siebesma et al. 2003; Zhu and Albrecht 2003; Heinze et al. 2015). With abundant evidence from the field, laboratory, and numerical simulations, it is extremely likely that under rather general circumstances, of the CBL is a universal function of .

The robustness of over a broad CBL stability range is a useful property that has been exploited in the physical parameterizations of the planetary boundary layer (PBL). The most frequent use of Eq. (1) is found in the formulation of the countergradient term in PBL schemes, in which a diagnostic profile of is required (Deardorff 1972; Holtslag and Moeng 1991; Zilitinkevich et al. 1999). It has also been used to parameterize the mass-flux term in the eddy diffusivity mass flux scheme of Siebesma et al. (2007). In addition, second- and higher-order statistics of *w* are also key to quite a few cumulus schemes used in regional and global models (see Berg et al. 2017, and references therein).

Judging from observations, an Eq. (1) type of scaling in the CBL requires the following external conditions: first, a relatively flat and horizontally homogeneous surface. Locations of all field experiments listed in Table 1 including ocean, farmland, and plain satisfy this criterion. Both topography and heterogeneous surface can influence the dynamics of the boundary layer, affecting vertical velocity variances in complex ways (Kaimal and Finnigan 1994, chapters 4 and 5). Second, quasi-steady-state conditions, as almost all field experiments selected data from the early noon to afternoon period when the CBL is well developed, such that second-order turbulence statistics reach steady state despite the slowly evolving mean flow. Berg et al. (2017) showed that acquired during morning and evening transition periods, when the CBL evolves rapidly, do not conform to the classical L80 profile. Third, negligible mesoscale circulations, as the mesoscale environment also affects the profile of . As noted by L12, in all three “outlier” cases of the LIFT experiment, either mesoscale roll circulation or mesoscale cellular convection (see Atkinson and Wu Zhang 1996) were observed. The additional energy at scales larger than the CBL increases , pushing it away from the classical scaling profile.

In the absence of mesoscale disturbances and over horizontally homogeneous land surfaces, do internal CBL dynamics affect ? Turbulence in the well-developed CBL is driven by both buoyancy and shear. The relative strength of the mechanical (production by vertical wind shear) to thermal (production by buoyancy) forcings determines the organization of CBL convection (see Salesky et al. 2017, and references therein). A measure of the relative importance of thermal to mechanical forcings is the bulk stability parameter^{2 }, where (m) is the Obukhov length, (m s^{−1}) is the surface friction velocity, and *κ* is the von Kármán constant. The range of *ζ* extends from zero to infinity indicating neutral stability and free convection, respectively.

Among the field experiments listed in Table 1, the observed *ζ* varies from near-neutral to highly convective values . A universal profile of over this broad stability range is intriguing. This is because the normalization factor is only a measure of the buoyancy forcing of the CBL, and does not account for shear, yet collapses over a wide range of wind shear. In other words, to a large extent, buoyancy production alone determines the vertical component of the turbulent kinetic energy (TKE), irrespective of the other TKE source term, which is shear production.

Unlike the vertical, the horizontal component of TKE in the CBL does not exhibit such universal scaling when normalized by . L80 stated that “the height dependence of and is somewhat obscured by the scatter,” although this could be because horizontal turbulent velocity is susceptible to fluctuations at the low-frequency (i.e., mesoscale) end of the spectrum (Caughey and Palmer 1979; L12). Free of mesoscale inhomogeneity, and from LES studies of idealized CBL do exhibit sensitivity to bulk CBL stability. For example, Gibbs and Fedorovich (2014) performed LES of a free convective and a sheared CBL. The latter case had a *ζ* of approximately 10. While the profiles of the two cases almost overlap, the and profiles exhibit significant differences throughout the depth of the CBL. Salesky et al. (2017) performed a series of LESs with *ζ* ranging from 0 to 1041, and showed monotonic decrease of with increasing *ζ*.

Motivated by the universal profile of , this study sets out to investigate the insensitivity of to CBL stability, or more precisely, wind shear. Data analyzed in this study are from LESs of several idealized convective and neutral cases as presented in section 2. In section 3, the simulated vertical and horizontal variances and their budgets are analyzed, followed by a rotational and divergent decomposition of the flow field to further examine the responses of the horizontal and vertical components of TKE to shear and buoyancy. An explanation that relates the universal profile to the dynamics and interactions of local and nonlocal eddy structures in the CBL is given in section 4.

## 2. Case description and numerical setup

This study is based on LESs of seven idealized CBL cases and a neutral boundary layer (NBL) case. A list of model parameters is presented in Table 2. The model setup follows Shin and Hong (2013) with prescribed surface heating and large-scale pressure gradient forcings. The “B” cases are driven by a constant heat flux of 0.20 K m s^{−1}. Case B1 represents an almost free-convective boundary layer with a weak 1 m s^{−1} geostrophic wind. The B5 and B10 cases have 5 and 10 m s^{−1} geostrophic winds, respectively. Case B10 has a bulk stability parameter *ζ* around 20, almost the same as that of the “moderately convective” case of L12. The prescribed surface heat flux for the four “S” cases is 0.05 K m s^{−1}. Case S10 with a geostrophic wind of 10 m s^{−1} has , similar to the “least convective” case of L12. While the first six cases in Table 2 feature barotropic conditions, the SG10 is a baroclinic case with a prescribed linear geostrophic wind. The neutral case N10 is forced by a constant 10 m s^{−1} geostrophic wind. All runs are initialized with the potential temperature profile from Shin and Hong (2013). The initial has a sharp inversion between 925 and 1075 m, which effectively caps the growth of the boundary layer in all simulations. The initial wind profile is geostrophic. The Coriolis parameter *f* is 10^{−4} s^{−1}, and the surface roughness length is 0.1 m.

All CBL cases are performed on a doubly periodic domain of 5.04 km × 5.04 km × 2 km, with Rayleigh damping applied to the top 500 m. Uniform 10-m horizontal and 4-m vertical grid spacings are adopted below 1.3 km, so that the CBL is well resolved with a negligible contribution from the subgrid scale (SGS), except for the first few grid points above the ground. Above 1.3 km, the vertical grid is stretched from 4 to 50 m over 24 grid points. The NBL case is performed with a finer 8-m isotropic grid spacing since the energy-containing turbulence length scale under neutral stratification is smaller than that under convective conditions. The domain size for the NBL is 9.216 km × 9.216 km × 1.5 km. The 1.5-order TKE closure of Moeng (1984) is adopted to parameterize SGS turbulence. The CBL reaches a statistical quasi-steady state after 6*τ*, where is the large-eddy turnover time, approximately 600 and 800 s for the B and S cases (Shin and Hong 2013). Therefore, the B cases are run for 4 h while the S cases are run for 6 h. The NBL case is run for 10 h to reach steady state. Data are output at every 10 min, and a total of 13 snapshots from the last 2 h of the simulations are analyzed.

The Advanced Regional Prediction System (ARPS) is used for the simulations. ARPS was developed at the Center for Analysis and Prediction of Storms at the University of Oklahoma. It is a nonhydrostatic mesoscale and convective-scale finite-difference numerical weather prediction model that is also suitable for LES (Chow et al. 2005). More details on ARPS are documented in Xue et al. (2000, 2001).

## 3. Results and discussion

The objective of this study is to understand the universal profile of the buoyancy-normalized vertical velocity variance over a broad range of CBL stability. As such, the sole premise of this work is that profiles of collapse to *some universal curve* in a well-developed CBL. Properties of that particular curve including its shape, magnitude, and empirical fit, however, are not the focus of this study. In fact, curve fits other than Eq. (1) have been proposed in the literature. For example, Eq. (15a) in Holtslag and Moeng (1991) is an alternative formulation that also takes the friction velocity into account, although its profile (see their Fig. 4) is similar to that of L80.

Traditionally, Eq. (1) has been regarded as a mixed-layer scaling (Stull 1988, chapter 9.6.3). The mixed layer is located between the unstable surface layer (occupying roughly the bottom 10% of the CBL) and the stably stratified entrainment zone (usually in the top 15%–20% of the CBL), and is characterized by intense turbulent mixing. Nevertheless, in many field observations, Eq. (1) often fits well over the depth of the CBL. Among those, the most convincing evidences are from the yearlong data of Berg et al. (2017, their Fig. 4a) and Lareau et al. (2018, their Fig. 6a), where the observed profiles exhibit universality between 0.05*z*_{i} and 1.2*z*_{i}.

Next, we define the parameter range of *ζ* over which the universal scaling applies. As discussed in section 1, observational evidence suggests that L80 scaling applies for . The exact lower bound of *ζ* is unclear, and is most likely related to the organization of turbulence as the boundary layer transitions from convective to neutral conditions. Salesky et al.’s (2017) LES results suggest that when *ζ* falls below 7.2, begins to deviate from the universal profile, exhibiting a second peak in the surface layer. However, L12 found universal scaling for their weakly unstable case with *ζ* = 5.9, and our results show that this lower bound falls somewhere between 5.8 (S10) and 2.4 (S15). In this study, we adopt *ζ* ~ *O*(1) as the qualitative lower bound for the universal profile.

Furthermore, we clarify the extent of “universality” of . Across the stability range of *ζ* > *O*(1), small differences are indeed found among the profiles. In L12, the most convective case has a slightly larger than those of their moderately and least convective cases. Berg et al. (2017) found small (~10% or less judging from their Fig. 9a) but statistically significant differences between the very unstable and moderately unstable cases based on a threshold *ζ* value of 30, at elevations between 0.05–0.45*z*_{i} and 0.65–0.75*z*_{i}. As shown in Fig. 1, such differences are small especially when compared to the changes in . However, recognizing their statistical significance, the term “near universal” is used hereafter.

### a. Mean profiles

Horizontally and time-averaged resolved variance profiles of and as functions of the normalized height are presented in Fig. 1. Angle brackets represent horizontal and temporal averaging. In this study, horizontal and time averaging are performed instead of ensemble (Reynolds) averaging. The two are equivalent when the size of the dominant eddies, which is on the order of the boundary layer depth in the CBL, is much less than the horizontal size of the domain (Wyngaard 2004). Lowercase *u*, *υ*, and *w* represent the difference between the grid-resolved velocity and its horizontal mean; is diagnosed as the elevation of the minimum sensible heat flux. The SGS variances (i.e., twice the SGS TKE) are much smaller than their resolved counterparts, and are not shown here. The L80 profile is also presented in Fig. 1a as reference.

The bulk stability parameter *ζ* of the seven CBL cases ranges from 2.4 for the weakly unstable case S15 to 1.5 × 10^{3} for the nearly free-convective case B1. General agreement in is obtained among all cases in Fig. 1a. The difference between LES and the L80 profiles can be largely attributed to the elevated location of peak variances in the LES profiles, except for the most unstable case B1, which overlaps the L80 profile below 0.4*z*_{i}. Case S15 deviates from the rest by exhibiting a bulge in the surface layer. The same transitional behavior of is also documented in Salesky et al. (2017). Overall, differences of the profiles among the CBL cases are small. (This will become more evident in Figs. 4 and 5 where the differences are directly computed and presented.) In contrast, differences of among cases are quite substantial compared to that of (notice the different range of *x* axes in Figs. 1a and 1b). Normalized horizontal variances in Fig. 1b increase monotonically with decreasing *ζ* at all elevations, also in agreement with the LES results of Salesky et al. (2017).

The vertical velocity spectra illustrate the partition of variances among scales. Figure 2 presents the two-dimensional time-averaged spectra of *w* obtained through a horizontal Fourier transformation following Sullivan and Patton (2011). The wavenumber-weighted spectra *kE* are presented to reveal the energy-containing scales (Kaimal and Finnigan 1994, chapter 2.3). Spectra are plotted at two selected elevations, one slightly above the surface layer at 0.15*z*_{i} and the other well inside the mixed layer at 0.5*z*_{i}. The normalized spectra exhibit −5/3 inertial subrange scaling at both elevations. The spectra of all seven cases agree well with each other except at the high-wavenumber end where the spectra drop off sharply. Such differences approaching the grid-cutoff resolution could be an artifact of the finite-difference numerical scheme. The B1 spectrum is somewhat different than the other six cases with a higher spectral density at wavenumbers on the order of the CBL depth (i.e., ). This is most evident at 0.15*z*_{i}, and less so at 0.5*z*_{i}. In the higher-wavenumber range ( between 5 and 30), the S15 spectrum is associated with larger spectral density than the rest. Examination of the spectra reveals that the causes for the larger values in the B1 and S15 profiles in the CBL surface layer shown in Fig. 1a are different. In the buoyancy-dominated case, the presence of more energetic boundary layer–scale eddies are responsible for the increase of . While in the shear-controlled case, the excess is associated with more vigorous inertial subrange-scale motions.

Figures 1a and 2 together reveal that not only the vertical velocity variances collapse to a near-universal profile, so do their spectral densities, although some moderate differences are found at the high and low ends of *ζ* (i.e., cases B1 and S15). Next, responses of turbulent velocity variances to shear and buoyancy forcings are quantified. In our idealized case setup (see section 2), shear and buoyancy are determined by the prescribed geostrophic winds and surface heat fluxes. The normalized profiles of wind shear are presented in Fig. 3. In the CBL, shear is concentrated in the surface layer and the entrainment zone. In the mixed layer in between, the magnitude of shear is small due to vigorous convective mixing. This is most evident by comparing the case-averaged CBL shear profile with that of the NBL. The heat flux profiles of all CBL cases show the classic linear shape and overlap each other (results not shown).

Next, we move on to investigate the response of turbulent velocity variances to the bulk stability parameter *ζ*. Since *ζ* is a measure of the interplay between buoyancy and shear forcings in the CBL, the effects of buoyancy and shear are first isolated and then examined separately. First, the influence of buoyancy is isolated by taking the difference of the turbulent velocity variances from three pairs of cases with the same geostrophic winds but different surface heating. Figure 4 presents the differences of the vertical and horizontal velocity variances between (B5, S5), (B10, S10), and (S10, N10). The former member of each pair is driven with a larger prescribed surface heat flux than the latter (see Table 1). With fixed geostrophic winds, turning up the surface heat flux increases both the vertical and horizontal turbulent velocity variances throughout the depth of the CBL. For all three cases, shows large increases near the surface and toward the top of the CBL, while the most significant incremental is found close to 0.4*z*_{i} near the center of the mixed layer. It is also interesting to note that the incremental change for the (B5, S5) pair is nearly identical to that of the (B10, S10) pair, even though the geostrophic winds differ by 5 m s^{−1} between the two pairs.

Next, the effect of wind shear is examined by comparing cases with the same surface heating but different geostrophic winds. Our case setup offers quite a few possible combinations, and we present five such pairs^{3} in Fig. 5. In each pair, the former member is forced with a stronger geostrophic wind. The S10 and SG10 pair is differentiated by the presence of geostrophic wind shear. Compared to Fig. 4, responses of turbulent velocity variances to shear are rather different from those to buoyancy. With stronger shear, profiles increase almost uniformly at all elevations, except for the (S10, SG10) pair where the incremental becomes negative at about 0.8*z*_{i} upward. This is likely due to the enhanced production by shear in the stably stratified entrainment zone in the SG10 case, where the geostrophic wind shear extends all the way to the domain top. This also implies that the strength of the capping inversion, which affects the entrainment rate, is another potential factor that may influence the horizontal variance profiles. In general, obtains the largest values near the surface, decreases upward until 0.4*z*_{i}, and remains more or less constant up until 0.9*z*_{i}. At the CBL top, has another smaller local maximum. The general shape of is related to the profiles of shear production shown in the following subsection. In comparison, profiles of remain nearly fixed as shear varies, as indicated by the near-zero profiles for all five pairs. Among them, the (B5, B1) pair has a noticeable negative between the surface and ~0.4*z*_{i}. But even so, is still small compared to for this particular pair.

To summarize the findings from Figs. 4 and 5, actively responds to both thermal and mechanical forcings. Its magnitude increases when either one of the forcings is strengthened. In contrast, is most strongly affected by changes in the thermal forcing, and largely disregards changes in wind shear. It is exactly such behavior that gives rise to the near-universal profile of over a broad stability range. Stated in another way, although thermal and mechanical productions are both TKE sources, the latter is only able to affect the horizontal component of TKE, while the former influences both the horizontal and vertical TKE. Therefore, the key to understanding the near-universal profile of lies in the reason why shear production does not share its energy input with the vertical TKE.

Here we clarify an important point. Given the characteristic shape of the CBL mean shear profiles in Fig. 3, it is tempting to think that shear does not affect the vertical velocity variance profiles simply because shear production is small in the mixed layer. This would be a valid argument if the nature of turbulence in the CBL were local. In other words, forcings (in this case shear) at a certain elevation were only able to affect the turbulence intensity in the vicinity of that particular elevation. However, the nonlocal nature of CBL turbulence undermines the localness explanation. In the CBL, organized convective eddies span the depth of the boundary layer providing vigorous turbulent mixing (Hunt et al. 1988; Lothon et al. 2006). Such nonlocal characteristics are evidenced in the incremental variance profiles in Fig. 5; in the mixed layer (0.2*z*_{i}–0.8*z*_{i}) where shear is close to zero, is certainly positive and much larger than .

### b. Budget analysis

To examine the effects of buoyancy and shear productions in the CBL, we turn to the Reynolds-averaged turbulent velocity variance equations. Under horizontally homogenous conditions and adopting the Boussinesq approximation for the boundary layer,

where the overline represents Reynolds-averaging operator and upper- and lowercase letters represent the mean and turbulent variables, respectively. We note that the vertical and horizontal turbulent velocity variances are the same as twice the vertical and horizontal components of the TKE . Therefore, “turbulent velocity variance” and “TKE” are used interchangeably. In Eq. (2), the right-hand-side (rhs) terms represent buoyancy production, turbulent transport, pressure correlation, return to isotropy,^{4} and dissipation, respectively. On the rhs of Eq. (3), the first two terms together represent shear production, the third and fourth are turbulent diffusion and return-to-isotropy terms, and the last term is dissipation (Stull 1988, chapter 4.3). The TKE equation [Eq. (4)] is obtained by adding Eqs. (2) and (3).

According to Eq. (4), buoyancy and shear production are both sources to the TKE, although the former contributes to only [Eq. (2)], and the latter to only [Eq. (3)]. Here, and are coupled through the return-to-isotropy terms, which sum up to zero under the Boussinesq approximation,

Among the remaining forcings, the turbulent transport and pressure correlation terms redistribute variances within CBL, and the dissipation terms are sinks of the TKE. Only resolved forcings are computed in this study. Dissipation is computed as the residual of its respective budget equation (see Heinze et al. 2015, and references therein). The residual method works well as indicated by obtaining a ratio of 2 for , in agreement with the generally accepted theory of local isotropy at small scales (see the appendix for details). Under quasi-steady-state conditions, the time tendency terms on the left-hand sides of Eqs. (2)–(4) are close to zero.

Based on the budget equations, we reinterpret the results of Figs. 4 and 5. When surface heat flux is increased, according to Eq. (2), the additional buoyancy production goes straight to increasing (Fig. 4b). However, part of the incremental thermal forcing must be channeled to the horizontal TKE through the return-to-isotropy terms, because is increased along with (Fig. 4a). Likewise, when the geostrophic winds are stronger, enhancement of follows directly from shear production [Eq. (3) and Fig. 5a]. However, in this case, the return-to-isotropy terms must remain nearly unchanged; otherwise would have been affected (Fig. 5b).

To evaluate the above interpretation, the return-to-isotropy term in Eq. (2) for all cases is presented in Fig. 6, can be obtained directly from according to Eq. (5) and therefore is not presented. The normalized profiles overlap quite well among all CBL cases except S15. This is consistent with our speculation that the return-to-isotropy terms are largely determined by the thermal forcing. Most CBL cases have negative indicating net transfer of energy from to . On the contrary, acts as a source term to for case N10, indicating conversion of shear production to in the absence of heating (results not shown). Despite the overall similarities, differences in profiles emerge as *ζ* decreases. A small increase in was first observed in case S10 close to the surface, where shear production is the strongest. As *ζ* drops further in S15, the vertical extent of the increase in expands to almost the top of the CBL. Below ~0.1*z*_{i}, even becomes positive, indicating conversion of to . Some moderate differences are also found in the entrainment layer, which is likely due to the increased shear production at that level.

Responses of the rest of the budget terms to changes in heating and shear are presented through three selected cases (B5, S5, S10) overlaid in Fig. 7. All terms are normalized by the surface buoyancy forcing . S5 serves as the control case whose differences with S10 and B5 are due to changes in shear and buoyancy forcings respectively. From S5 to S10 in Fig. 7a, shear production exhibits a marked increase from the surface up to about 0.4*z*_{i} as a result of the increase in geostrophic winds. The dissipation rate is also increased accordingly to balance the increase in shear production. Turbulent transport of is apparently insensitive to the changes of shear as shown in Fig. 7a. In the budget in Fig. 7c, buoyancy production, turbulent transport, and pressure redistribution terms are almost identical for S5 and S10, given the same normalization factor for the two cases. The only differences are observed in the return-to-isotropy and dissipation terms, consistent with the changes in the budget.

In comparison, changes in surface heating between S5 and B5 have almost no effect on the normalized budget terms for both and , especially the latter (see Figs. 7b and 7d). Small differences are found in the normalized shear production term below ~0.2*z*_{i} in Fig. 7b. This is mainly caused by normalization with . The unnormalized shear production profiles of the two cases are nearly the same (results not shown). The difference in normalized shear production between the two cases is largely balanced by changes in the normalized dissipation term. Overall, the results suggest that when shear is fixed, buoyancy only changes the magnitude of budget terms without qualitatively affecting the balance among them.

By analyzing the responses of TKE forcing terms to shear and buoyancy, it is found that the return-to-isotropy terms, which facilitate the conversion between and , are largely set by buoyancy alone. Changes in shear do not cause appreciable changes to the return-to-isotropy terms until *ζ* is approximately *O*(1). This explains why does not respond (through the return-to-isotropy terms) to changes in shear. On the contrary, when buoyancy production is changed, all budget terms including the return-to-isotropy terms adjust their magnitudes accordingly such that their buoyancy-normalized profiles remain unchanged. Therefore, both and are affected by changes in buoyancy production.

Why does the return-to-isotropy term stay invariant when changes with shear production? One plausible hypothesis is that the incremental turbulent velocity vector (*δu*,*δυ*) as a result of the changes in shear is nondivergent (i.e., ), such that the return-to-isotropy term remains unchanged,

To evaluate this hypothesis, we perform a rotational and divergent (R&D) decomposition of the horizontal turbulent velocity, and analyze their respective energy budgets in the next subsection.

### c. A perspective from rotational and divergent winds

The idea of decomposing the horizontal wind vector into rotational and divergent parts,

dates back to the early of work of Wipperman (1957), Wiin-Nielsen (1968), and Chen and Wiin-Nielsen (1976). The method has been applied to the analyses of global circulation (Chen and Wiin-Nielsen 1976), synoptic-scale monsoon circulation (Krishnamurti and Ramanathan 1982), mesoscale convection (Buechler and Fuelberg 1986), midlatitude cyclones (Chen et al. 1978), and tropical cyclones (Cao et al. 2014). Previous studies have found that the kinetic energy of the divergent flow is much smaller than that of the rotational flow . The percentage divergent energy is estimated to be about 1% globally (Chen and Wiin-Nielsen 1976), and about 2% in mesoscale convective systems (Buechler and Fuelberg 1986). Despite its trivial contribution to the total kinetic energy, the divergent flow plays a key role in the energy cycle, because all energy converted from available potential energy to must go through . Therefore, the role of the divergent flow has been described as a “catalytic” one that regulates the dominant rotational flow. To the best of our knowledge, R&D decomposition has not been applied to PBL turbulence before; neither have the characteristics and interactions of the R&D winds in the PBL been systematically documented.

The Helmholtz decomposition is applied to obtain the rotational/nondivergent (*u*_{r}, *υ*_{r}) and divergent/irrotational (*u*_{d}, *υ*_{d}) winds.^{5} The horizontal TKE is decomposed into R&D parts where

In a horizontally periodic domain, the planar average of is zero, so that

Vertical profiles of normalized R&D variances are presented in Fig. 8. The normalized divergent variances agree well among all CBL cases. The general profile of obtains its maximum value at the surface, decreases upward until a local minimum is reached at ~0.5*z*_{i}, and increases again toward the CBL top. Above ~0.9*z*_{i}, the divergent variances reach an elevated local maximum and drop quickly upward to zero at ~1.1*z*_{i}. Given the near-universal profile of , Eq. (7) suggests that the normalized rotational variances must be responsible for the patterns of the horizontal variances presented in Fig. 1b. Indeed, in Fig. 8b increases with decreasing *ζ*. The largest increase is found near the surface where vertical wind shear is the strongest.

The general patterns of the R&D variances presented in Fig. 8 support the hypothesis proposed at the end of section 3b that changes in shear only induce perturbations in the nondivergent (i.e., rotational) part of the horizontal turbulent velocity. As such, the return-to-isotropy term remains unchanged, and as a result, also stays the same despite changes in , or more specifically, its rotational component . While is affected by shear, is largely determined by the buoyancy forcing alone as shown in Fig. 8.

#### 1) Characteristics of the rotational and divergent winds in the CBL

Before further investigation, it is worthwhile to briefly document the characteristics of and in the CBL, with a focus on their relation to shear and buoyancy. A snapshot of the horizontal flow field of case B10 at 0.15*z*_{i} at the end of 4 h is presented in Fig. 9. Other cases show qualitatively similar features and are not presented here. The vertical velocity *w* in Fig. 9a exhibits a mixture of streaky and cellular patterns, because the convective organization of B10 is in transition from cells to roll vortices (*ζ* smaller than about 15) (Shin and Hong 2013). The turbulent spanwise velocity *υ* and its divergent and rotational components are presented in Figs. 9b–d. Components and show qualitatively similar features and are not presented here. Comparing Fig. 9c to Fig. 9d reveals striking differences between and . The former exhibits the same organized patterns found in *υ*, although much smoother. The small-scale fluctuations in *υ* are found in , which lacks apparent coherence.

The spatial scales of the R&D winds are characterized in Fig. 10 through the energy spectra. At 0.15*z*_{i}, the divergent and rotational energy spectra contribute to the horizontal energy spectrum at different scales. The divergent part is dominant at the small-wavenumber end (i.e., the boundary layer scale), which explains the presence of large coherent structures in the divergent wind in Fig. 9c. As the wavenumber increases, drops off quickly. Over the inertial subrange, is about an order of magnitude smaller than . The low energy density at the intermediate and small scales gives the divergent wind its smooth appearance in Fig. 9c. The rotational part , on the other hand, is dominant over the inertial subrange, but is small at boundary layer scale. The crossover of and at 0.15*z*_{i} occurs at about , where is the wavelength. Higher up at 0.5*z*_{i} where shear is small, is approximately equal to across the entire spectrum, while remains small at all scales.

To facilitate the upcoming discussion on momentum fluxes, the vertical energy spectra of B10 from Fig. 2 is presented again in Fig. 10. At 0.15*z*_{i}, is small at the boundary layer scale since the vertical scale of *w* is restricted by the distance to the ground. The crossover of and also occurs at about . Over the inertial subrange, is much larger than . Higher up at 0.5*z*_{i}, is comparable to and .

Height–wavenumber contours of the energy spectra of case B10 are presented in Fig. 11. As in Fig. 10, the boundary layer–scale energy peaks of in the surface layer and the entrainment zone (Fig. 11a) reside almost entirely in the divergent wind (Fig. 11c). These two concentrated patches of intense spectral density are the only places where is of an appreciable magnitude; in Fig. 11d comprises the rest of . Its most energetic scale increases with height from the surface upward, reaches an asymptotic value of about at 0.4*z*_{i}, and stays nearly constant up to ~0.8*z*_{i}. Such behavior is similar to the classic Blackadar length-scale formulation for the eddy viscosity in the purely shear-driven NBL (Blackadar 1962). The vertical energy spectra are presented in Fig. 11b. The most intense spectral density is found in the mixed layer at length scales of around 1.0*z*_{i}, where vertical motions are most vigorous.

The horizontal energy spectra of case N10 is also examined to allow a general qualitative comparison of the R&D energy partition in the NBL and the CBL. Comparing Fig. 12 to Fig. 10, the most distinct difference is that is now much smaller than at both elevations. In fact, this is true over the entire depth of the NBL in the absence of heating (results not shown).

The streamwise momentum fluxes associated with the R&D winds are presented in Fig. 13. The spanwise momentum fluxes are close to zero and are not presented. Momentum fluxes for B1 and B5 are not presented since they are very small, and show large wiggles when normalized with . Minimizing such uncertainties requires a wider domain and/or longer integration time that are not available in the present dataset. The vertical profiles of of the five CBL cases agree well with each other in Fig. 13c, except for case S5, whose profile overlaps with that of N10 from the surface to ~0.6*z*_{i}. The exceptional case of S5 is also likely subjected to uncertainties associated with normalization of a small momentum flux. In general, the normalized momentum fluxes of the CBL are only slightly larger (in terms of magnitude) than that of the NBL over the depth of the boundary layer.

The divergent momentum fluxes in Fig. 13a exhibit two interesting features. First, unlike , the sign of is positive indicating *upward* transfer of momentum against the positive mean wind shear . Similar upward-directing divergent flux is also found in the linear Rayleigh–Bénard convection (RBC) with constant background shear (Emanuel 1994, chapter 3.4). This is likely one of the many similarities between the extremely high Rayleigh number nonlinear boundary layer convection and the linear RBC, which is an active area of research (see Pandey et al. 2018, and references therein), but out of the scope of this study. As a consequence, by depleting the divergent TKE, the countergradient feeds TKE back into the mean flow through turbulent production, or in this case “destruction.” Unfortunately, we do not have a theory to explain the upward-directed and have to leave this for future research. To compensate for the positive , the rotational momentum fluxes must transfer extra amounts of momentum downward, hence their larger magnitudes than as shown in Fig. 13b.

The second notable feature of Fig. 13a is that is substantially smaller than and . The normalized decreases with *ζ* and approaches zero in the neutral limit. Its small magnitude can be understood in the light of scale interactions between and *w* in wavenumber space. As shown in Figs. 10 and 11, the spectra have sharp peaks at boundary layer scales in the surface layer and the entrainment zone, but drop off rapidly at smaller scales. In these two regions of the CBL, the *w* spectra have little energy density at boundary layer scales due to the constraints on vertical motions imposed by the physical ground and the overcapping inversion (Kaimal and Finnigan 1994, chapter 2.6). Scale separation of the energetic spectral modes of and *w* is most evident in Fig. 10. Although still subjected to debate, scale interaction of spectral modes in turbulent flow has been shown to be predominately local (i.e., between wavenumbers and where ) (Pope 2000; Cardesa et al. 2017), hence the small correlation between and *w* due to their scale separation. In the mixed layer in between, is also small because the spectral density of the divergent flow is much smaller than that of *w* (i.e., ).

We conclude this subsection by briefly summarizing the characteristics of the R&D winds in the CBL. The divergent wind contains the large coherent CBL structures that are essentially buoyancy driven. Without heating, is suppressed at all elevations in the NBL. Similar to , the vertical profiles of are determined by alone, which is a measure of the CBL buoyancy. The divergent momentum fluxes are small compared to the total momentum fluxes, and are directed upward against the mean shear. In comparison, the rotational wind is most energetic over the inertial subrange. With heating fixed, its associated TKE increases with shear. The rotational momentum fluxes are largely (solely) responsible for the turbulent mixing of momentum in the CBL (NBL). It is worth noting that the R&D spectra in Fig. 10 are very similar to those presented in Hellsten and Zilitinkevich (2013, their Fig. 3), who separated the semiorganized convective structures from turbulence with several empirical averaging and filtering approaches.

#### 2) Budget analysis of the rotational and divergent TKE

With the characteristics of the R&D winds in mind, we continue to investigate the reason why changes in shear are only able to affect horizontally nondivergent wind. We first present the R&D TKE budgets (the derivations largely follow Buechler and Fuelberg (1986) and are not repeated here for brevity):

where (s^{−1}) is the vertical vorticity. The first terms in square brackets of Eqs. (8) and (9) have the same form with opposite signs. They represent conversion between and . The second terms in parenthesis represent shear production/destruction. The third and fourth terms of Eq. (9) are the same turbulent diffusion and return-to-isotropy terms in Eq. (3). The last terms of both equations are dissipation rates, which are obtained based on the residual method.

The and budgets averaged over all CBL cases except for S15 are plotted in Fig. 14 to present their general characteristics; S15 is excluded due to its transitional nature as shown in Fig. 6. In Fig. 14a, the budget is effectively characterized as a balance of conversion, turbulent diffusion, and return to isotropy. The conversion term acts as a sink to and a source to except close to the surface where the direction of energy transfer is reversed. The return-to-isotropy term converts TKE from to , which is then passed along to through the conversion term. Turbulent diffusion tends to homogenize over the depth of the CBL. For the two remaining terms, shear production by the divergent momentum flux is slightly negative close to the surface and zero above, primarily due to the small magnitude of as shown in Fig. 13a. Dissipation of is also close to zero, which is expected given the spectral characteristics of the divergent wind as presented in Figs. 10 and 12, specifically the dominance of boundary layer–scale organized convection and the lack of inertial subrange eddies. In agreement with the theory of turbulence cascades, direct dissipation of the boundary layer–scale large eddies is negligible.

The budget of in Fig. 14b consists of only three terms, with conversion and shear production as sources and dissipation as the sink. Note that in the absence of turbulent diffusion in the budget, almost all shear production of is dissipated locally, although a small portion close to the surface is converted to . This rejects the speculation made at the end of section 3a on the vertical redistribution of shear production by nonlocal turbulent mixing. In the mixed layer, the energy source of comes mostly from through conversion, rather than the suspected upward turbulent mixing of surface layer shear production.

Budget analyses of Figs. 7 and 14 together suggest the following schematic of energy transfer in the CBL as illustrated in Fig. 15. Thermal energy enters the system by generating through buoyancy production. It is then converted to through the return-to-isotropy forcing. After being somewhat homogenized vertically by turbulent diffusion, the energy is further passed on to , mostly in the mixed layer where it finally gets dissipated. For mechanical energy, the transfer route is much simpler. It goes almost exclusively into below 0.4*z*_{i} and in the entrainment zone, and gets dissipated directly. In the CBL, the role of in the overall energy cycle is to channel buoyant energy into , much like the previous findings at the global-, synoptic-, and mesoscales. Note that in the NBL without heating, processes III and IV in Fig. 15 are reversed; receives part of the mechanical energy injected to through the conversion term and passes it on to through return to isotropy (results not shown).

Going back to the question posed at beginning of this subsection, we examine the possible routes through which shear could have affected . As shown in Fig. 15, mechanical energy could influence through either *direct* shear production or *indirect* conversion from . Figure 14a has shown that direct shear production plays a minimal role in the budget. As explained earlier, this is likely due to the spectral gap between the energetic modes of *u*_{d} and *w* leading to the very limited shear production by divergent momentum fluxes.

The indirect route involves shear production of and subsequent conversion to . As shown in Fig. 14, net conversion of TKE occurs in the -to- direction except close to the surface. The exact cause of the preferential -to- conversion is unclear to us, but given the spectral characteristics of the R&D winds presented in Figs. 10 and 11, transfer of energy from to would be in the upscale direction from small to large eddies. This seems highly unlikely in the presence of vigorous turbulent mixing in the CBL. Inverse cascade occurs in quasi-2D turbulence and is the dominant process in large-scale circulation (Gage and Nastrom 1986), but for boundary layer–scale eddies, downscale cascade should dominate. As such, downscale energy transfer from to is much more plausible (furthermore, Fig. A3 in the appendix reveals that there is little spread among the vertical profiles of the conversion term of all cases except S15, irrespective of the changes in shear). Together, they suggest that TKE conversion occurs primarily in the -to- direction, while the conversion rates are largely determined by buoyancy alone. The absence of direct production and the restricted indirect effects together explain why shear has a very limited influence on the horizontally divergent winds.

## 4. Summary

The classic L80 profile of in the CBL has been observed repeatedly in field observations, laboratory experiments, and numerical simulations. Over a horizontally homogeneous land surface and in the absence of mesoscale variability, the state of the CBL is determined by surface heating and vertical wind shear, whose combined effects can be grouped into a bulk stability parameter *ζ*. The reason why the profile in the CBL is considered near universal is that it is applicable to a broad spectrum of *ζ* ranging from free-convective to near-neutral conditions. This study aims to understand the resilience of to changes in CBL stability, or more precisely, wind shear, since effects of buoyancy are already absorbed into the normalization factor .

Large-eddy simulations of seven CBL and one NBL cases are performed with different thermal and mechanical forcings. By constructing the horizontal and vertical TKE budgets, it is shown that buoyancy and shear production are respective sources to and , which are coupled through the return-to-isotropy terms. Shear production only affects indirectly through -to- conversions, and likewise for buoyancy to act on . Analysis of the return-to-isotropy terms reveals that the net energy conversion is from to throughout the depth of the CBL for . Furthermore, the magnitude of the conversion rate is largely set by buoyancy alone except close to the surface where some influences of shear is found as *ζ* approaches *O*(1). This explains why shear (production) does not share (its energy input with ).

To understand the null response of the return-to-isotropy terms to shear, it is hypothesized that the incremental horizontal turbulent velocity as a result of the changes in shear is divergence free. Separation of the horizontal turbulent winds into rotational **u**_{r} and divergent **u**_{d} components and subsequent analysis confirms this hypothesis. It is found that similar to , the normalized divergent variance profiles collapse to a near-universal curve, while stability dependence is found in the normalized rotational (divergence free) variance.

Further analysis of the spectral characteristics of the divergent and rotational winds relates to the former with boundary layer–scale organized convection and the latter with inertial subrange eddies. This suggests a connection of the reason why changes in shear do not induce horizontally divergent velocity perturbations to the dynamics and interactions of eddy structures. With the physical representation of **u**_{r} and **u**_{d} in mind, examination of the rotational and divergent TKE budgets reveals the following:

Direct shear production of is almost negligible because of the insignificant divergent momentum flux . A plausible explanation is the separation of energetic scales between and

*w*, hence their poor correlation. In both the CBL and NBL, direct shear production produces almost exclusively horizontally divergence-free TKE.Indirect effects of shear production on through -to- conversion is limited, as TKE is converted preferentially from to . Physically, energy is converted from the boundary layer–scale convective motions to the inertial subrange-scale eddies, in accordance with the generally expected downscale energy transfer. The conversion rate is largely determined by buoyancy alone except in the surface layer.

With both limited direct and indirect influences of shear production on , the resulting horizontal velocity perturbations in response to changes in shear are largely divergence free. As conversion between and requires horizontally divergent winds, remains unaffected by changes in shear, giving rise to the near-universal profile of . Dynamics and interactions of local and nonlocal eddies, specifically the preferential interactions between local wavenumbers and the downscale energy cascade under convective conditions, offer plausible explanations for such behavior at the fundamental level. However, neither characteristic of turbulence has been rigorously and fully established, especially the local interaction theory. Therefore, our explanation remains a plausible speculation at this point. Hopefully, continued study of the fundamental characteristics of turbulence will offer more support to our theory in the near future.

## Acknowledgments

We thank Hao Fu for suggesting the R&D decomposition and his MATLAB code, and Prof. Ming Xue for helpful discussions. We thank Dr. Don Lenschow and two other reviewers for their constructive comments. This work was supported by the National Key R&D Program of China (Grant 2018YFC1506802), and the National Natural Science Foundation of China (Grant 41875011). Simulations are performed at the High Performance Computing Center of Nanjing University.

### APPENDIX

#### Variance Budgets

Vertical budget profiles of , , and are presented in Fig. A1. The dissipation rate *ε* is obtained from the 1.5-order TKE closure, and *ε*/3 and 2*ε*/3 are assigned to the dissipation of and , respectively. The residual term is obtained by summing up all forcings. The forcings are case averaged except for S15, since it deviates significantly from the rest of the CBL cases. The general balance among budget terms agree with classic textbook findings (Wyngaard 2010, chapter 11), and is not discussed further. The residual terms are close to zero except for the first few grid points above the surface and around the CBL top, where resolved turbulence are suppressed. This suggests that the TKE budgets are quite well closed with resolved forcings for the bulk of the CBL.

In section 3b, the residual method is used to compute the dissipation rate, which effectively absorbs the residual term into *ε*. Ratios of the resulting and are presented in Fig. A2 for all CBL cases, along with the case-averaged profile. The ratio is close to 1 in the CBL, in agreement with the generally accepted theory of local isotropy at small scales, except near the surface and in the entrainment zone where subgrid-scale forcings could be important.

Vertical profiles of the normalized conversion term for all CBL cases are presented Fig. A3. Case S15 aside, the normalized conversion terms exhibit near universality among the CBL cases above the surface layer. The positive signs of the conversion terms suggest that TKE flows from to . Systematic influences of shear are observed mainly below 0.2*z*_{i}. As the bulk stability *ζ* decreases, the conversion term becomes smaller in the surface layer, and is eventual negative when *ζ* reaches 10 (i.e., case SG10). In S15 where *ζ* is of *O*(1), the vertical extent of -to- conversion further increases to 0.2*z*_{i}, and the magnitude of the conversion term is also much larger.

## REFERENCES

*Atmospheric Convection*. Oxford University Press, 580 pp.

*Atmospheric Boundary Layer Flows: Their Structure and Measurement*. Oxford University Press, 304 pp.

*Turbulent Flows*. Cambridge University Press, 771 pp.

*An Introduction to Boundary Layer Meteorology*. Kluwer Academic, 670 pp.

*Turbulence in the Atmosphere*. Cambridge University Press, 406 pp.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

Lareau et al. (2018) noted the exception that shallow cumulus of medium cloud cover (i.e., cloud fraction between 0.3 and 0.5) has a somewhat larger than the classical scaling profile.

^{2}

The bulk stability parameter is sometimes also expressed as , where .

^{3}

Other combinations can be inferred from these five pairs, for example, the (S15, S5) pair is obtained by adding (S10, S5) and (S15, S10).

^{4}

The return-to-isotropy term is sometimes also referred to as the pressure redistribution term (Stull 1988).

^{5}

Strictly, there should be a third harmonic component (*u*_{h}, *υ*_{h}) from the Helmholtz–Hodge decomposition. However, the harmonic component is identically zero for periodic domains (Cao et al. 2014).