On the Non-monotonic Variation of the Entrainment Buoyancy Flux with Wind Shear

The magnitude of the entrainment buoyancy flux, and hence the growth rate of the convective boundary layer, does not increase monotonically with wind shear. Explanations for this have previously been based on wind-shear effects on the turbulence kinetic energy. By distinguishing between turbulent and non-turbulent regions, we provide an alternative explanation based on two competing wind-shear effects: the initial decrease in the correlation between buoyancy and vertical velocity fluctuations, and the increase in the turbulent area fraction. The former is determined by the change in the dominant forcing; without wind shear, buoyancy fluctuations drive vertical velocity fluctuations and the two are thus highly correlated; with wind shear, vertical velocity fluctuations are partly determined by horizontal velocity fluctuations via the transfer of kinetic energy through the pressure–strain correlation, thus reducing their correlation with the buoyancy field. The increasing turbulent area fraction, on the other hand, is determined by the increasing shear production of turbulence kinetic energy inside the entrainment zone. We also show that the dependence of these conditional statistics on the boundary-layer depth and on the magnitude of the wind shear can be captured by a single non-dimensional variable, which can be interpreted as an entrainment-zone Froude number.


Introduction
The atmospheric boundary layer constitutes an important component of many weather forecasting and climate models. Significant progress has been made in our understanding of the B Juan Pedro Mellado juan.pedro.mellado@uni-hamburg.de atmospheric boundary layer (LeMone et al. 2019), particularly the convective boundary layer (CBL), yet models still struggle to appropriately represent the entrainment process, by which air from the free atmosphere is incorporated into the CBL (Edwards et al. 2020).
Of particular importance for models, due to its relevance to CBL growth, is the entrainment buoyancy flux. The entrainment buoyancy flux refers to the negative buoyancy flux at the CBL top as a consequence of entrainment and the stable stratification in the free atmosphere, and it is often characterized by its minimum value, b w z i, f . This value behaves differently depending on the magnitude of the wind shear and the boundary-layer depth, as illustrated in Fig. 1a in terms of the Froude number, Fr 0 , and the normalized CBL depth, z enc /L 0 (exact definitions are provided in Sect. 2). For weak shear conditions, the magnitude of the entrainment buoyancy flux, and concomitantly the CBL growth rate, has been noted to remain the same or even decrease by approximately 10−20% compared to shear-free conditions. This has been observed in laboratory studies ) and numerical simulations (Conzemius and Fedorovich 2006;Pino and de Arrellano 2008;. Figure 1 shows a 10% decrease in the magnitude of the entrainment buoyancy flux from the case Fr 0 = 0 to the case Fr 0 = 5. Although the observed decrease is small and could arguably be ignored in models, understanding why it occurs may shed light on the entrainment-zone dynamics and help to parametrize other entrainment-zone properties. When wind shear is moderate, there is little change from shear-free conditions, whilst under strong shear conditions, the magnitude of the entrainment buoyancy flux increases (Pino et al. 2003;Fedorovich and Conzemius 2008).
Various explanations have been given for these individual behaviours for weak, moderate, and strong wind shear. The decrease in the magnitude of the entrainment flux under weak shear conditions has been associated with a shear-sheltering effect (Hunt and Durbin 1999;Fedorovich and Thäter 2001) and with the redistribution of kinetic energy from vertical velocity to horizontal velocity fluctuations ). The increase under strong shear conditions is typically attributed to the additional shear production of turbulence kinetic energy (TKE), which creates larger amplitude undulations of the CBL top and stronger buoyancy and vertical velocity fluctuations (Kim et al. 2003). In this paper, we extend this view by conditioning the analysis of the entrainment zone into turbulent and non-turbulent regions.
To this aim, we use recent results of Fodor and Mellado (2020), who found that, although turbulent regions contribute the most to the entrainment buoyancy flux, wind shear does Fig. 1 The buoyancy flux as a function of, a height at z enc /L 0 = 20, b time or CBL depth at z = zi, f, and c entrainment-zone Froude number, Fr e (Eq. 6). Data in panel c are normalized by the shear-free value. Here and in subsequent figures, lines correspond to a running average over an interval z enc /L 0 = 2 and shadow regions indicate an interval of two standard deviations around that average not much change the magnitude of the flux itself within those regions. Although buoyancy fluctuations do become stronger with increasing wind shear, the effect is cancelled out by a reduction in the correlation between buoyancy and vertical velocity fluctuations. The primary factor determining the shear enhancement of entrainment is the increase in the turbulent area fraction. That study, however, considered only two regimes, namely, a shear-free case and a strong shear case. Here we extend it to also provide an explanation for the decrease in the entrainment flux under weak shear conditions, covering in this way a complete range of wind-shear effects on the barotropic CBLs in the quasi-steady regime.
A second contribution concerns how to parametrize these wind-shear effects as a simple function of the environmental conditions. Conventional statistics can be expressed in terms of a single non-dimensional variable that can be interpreted as a local Froude number in the entrainment zone . This is illustrated in Fig. 1c for the case of the entrainment buoyancy flux. Here we assess the ability of this entrainment-zone Froude number to also characterize conditional statistics.

Simulations
We perform direct numerical simulations of a cloud-free barotropic CBL that develops over a flat, aerodynamically smooth surface and evolves into a linearly stratified free atmosphere with constant buoyancy gradient, N 2 0 . Forcing is provided by a constant and homogeneous surface buoyancy flux, B 0 , and a constant wind speed in the free atmosphere, U 0 . Full details of the governing equations, dimensional analysis, and simulation set-up are provided in .
An overview of the simulation properties is given in Table 1 . All simulations are conducted at a reference buoyancy Reynolds number of Re 0 = 25, where and ν is the kinematic viscosity. The magnitude of wind shear is characterized by the reference Froude number where L 0 is a reference Ozmidov length defined as This length characterizes the entrainment-zone thickness and thereby other entrainment-zone properties in shear-free conditions (Garcia and Mellado 2014;. Due to statistical homogeneity in the horizontal directions, statistical properties depend only on height, z, and time, t. The upper limit of the mixed layer is well characterized by the encroachment length where z ∞ is located far enough into the free atmosphere for the integral to be approximately independent of z ∞ , angled brackets denote averaging in the horizontal planes, and b is the buoyancy field. Integrating vertically the evolution equation for the mean buoyancy, one can obtain the following one-to-one relationship between the encroachment length and time where t 0 is a constant of integration. Hence, the independent variables {z, t} are expressed in the non-dimensional form {z/z enc , z enc /L 0 } throughout. For the sake of generality, we present the results in a non-dimensional form. This allows us to formally embed the dependence on several parameters and variables into a reduced set, which is convenient for the proposed parametrizations. It also facilitates the comparison of data from simulations and measurements under different environmental conditions. Considering N 0 ≈ 0.006 − 0.018 s −1 , B 0 ≈ 0.001 − 0.01 m 2 s −3 , U 0 ≈ 0 − 20 m s −1 , and z enc ≈ 500 − 2000 m as typical midday conditions over land, one finds L 0 ≈ 20 − 200 m, z enc /L 0 ≈ 5 − 50, and Fr 0 ≈ 0 − 85. We are interested in the transition from weak to strong wind-shear effects and, as shown below by the results, this is well covered by the range of Froude numbers Fr 0 = 0 − 25 considered herein. The independent variable z enc /L 0 characterizes the state of the CBL development. The CBL reaches a quasi-steady state from z enc /L 0 ≈ 15, but we also consider earlier states of development where the magnitude of the entrainment buoyancy flux has been noted to decrease under weak shear conditions. To better characterize such conditions, we use an approximate fourfold larger domain size for the cases Fr 0 = 0 and Fr 0 = 5 to increase statistical convergence, namely 420 L 0 × 420 L 0 instead of 215 L 0 × 215 L 0 (see Table 1).
The only non-dimensional parameter whose atmospheric value cannot be matched in the simulations is the Reynolds number. Turbulence-resolving simulations, both direct numerical simulations and large-eddy simulations, are restricted to moderate Reynolds numbers because of limited computational resources. However, simulations can reach Reynolds numbers that are high enough for relevant properties to show some degree of Reynolds number similarity, which allows certain extrapolation of results to atmospheric conditions (Mellado et al. 2018;Mellado 2019). Here, we choose direct numerical simulations to reduce the uncertainty associated with turbulence models and numerical artefacts. This choice proves particularly convenient for entrainment-zone properties, where small scales become more relevant than in mixed-layer properties.
The entrainment flux ratio is defined as the negative minimum buoyancy flux normalized by the surface buoyancy flux (Fig. 1b). Since we are primarily investigating how wind shear affects this quantity, much of our analysis will take place at the height of minimum buoyancy flux, z i, f . A further reference height we consider is the height of maximum mean buoyancy gradient, z i,g , which lies near the top of the entrainment zone.
One last variable that we use is the entrainment-zone Froude number where measures the velocity difference across the entrainment zone between the free atmosphere and the mixed layer of the CBL, u is the streamwise component of the velocity field, and is the entrainment-zone thickness in free convection.  showed that wind-shear effects on entrainment-zone properties can be well represented by this Froude number (see also Fig. 1c), and this result was used to derive bulk parametrizations that smoothly vary with wind speed ). This previous work, however, only considered conventional statistics. Here we assess how well this entrainment-zone Froude number characterizes entrainment-zone properties conditioned into turbulent and non-turbulent regions.

Conditional Analysis
To partition the flow into regions of turbulent and non-turbulent fluid, we place a threshold on the enstrophy, below which the fluid may be considered approximately irrotational (da Silva et al. 2014). We choose this threshold based on where the probability density function (p.d.f.) of enstrophy has a saddle point (Borrell and Jiménez 2016). This saddle point has been shown to correspond to where the turbulent volume fraction becomes approximately insensitive to the chosen threshold (Watanabe et al. 2018). The p.d.f.s are shown in Fig. 2, with the enstrophy normalized by which is a reference value for the enstrophy in the mixed layer (Fodor and Mellado 2020).
To locate the saddle point, first, we find the median of the p.d.f. at each height (black profiles in Fig. 2). The saddle point is then approximated as the location of the minimum p.d.f. value along those median profiles between z enc and 1.5 z enc . We verified that results do not depend on order-of-one changes of these two limits. With these two limits, we exclude the near-surface region, where there is another local minimum (Fig. 2), and we also exclude the region next to the top boundary of the computational domain where the free-slip boundary condition introduces another local minimum. The development of the height and enstrophy value of the saddle point is shown in Fig. 3 . When Fr 0 is less than Fr 0 ≈ 10, wind shear has little effect on the height of the saddle point, z i,ω , and one finds that the average value of z i,ω is comparable to z i, f (see also Table 1). Under stronger shear conditions, z i,ω rapidly moves up beyond z i,g . The enstrophy value of the saddle point is correspondingly smaller at higher Fr 0 , as the saddle point is closer to the non-turbulent region (Fig. 3b). The actual threshold that we use for the partitioning in each case is the mean value of the respective curve in Fig. 3b. As a sensitivity analysis, we verified that using the same threshold for Fr 0 = 5 and Fr 0 = 10 as in the shear-free  Fig. 3 Evolution of, a the height, and b the enstrophy value of the saddle point in the enstrophy p.d.f.s (Fig. 2) case only affects the magnitudes of the quantities we look at by 10% at most and does not alter our conclusions. A detailed sensitivity analysis to the chosen threshold was provided for the Fr 0 = 0 and Fr 0 = 20 cases in Fodor and Mellado (2020), where it is shown that a comparable threshold in both cases does not alter any of the conclusions. The same outcome may be expected for all cases studied here, given that the thresholds at Fr 0 = 5 and 10 are similar to that at Fr 0 = 0, and the thresholds at Fr 0 = 15 and 25 are similar to that at Fr 0 = 20.
For stratified flows in general, potential vorticity may serve as a more suitable turbulence indicator than enstrophy, as in the absence of diffusion, gravity waves possess zero potential vorticity and there is thus a greater scale separation between the turbulent and non-turbulent regions (Riley and Lelong 2000;Watanabe et al. 2016). For the CBL in particular, however, Fodor and Mellado (2020) showed that at the Reynolds numbers we are able to achieve in simulations, the gain in scale separation when considering the potential vorticity is marginal. In any case, in the real atmosphere the scale separation between enstrophy values in the mixed layer and in the free atmosphere is on the order of 10 6 , making it just as good a turbulence indicator as potential vorticity. Moreover, Fodor and Mellado (2020) showed that the potential vorticity leads to an unrealistic turbulent area fraction profile due to the small buoyancy gradient in the mixed layer. For these reasons, we use the enstrophy and not the potential vorticity as the conditioning variable.

Results
Once the flow has been partitioned into turbulent and non-turbulent regions, we can express the buoyancy flux as follows Here, w represents the vertical velocity field, primes denote fluctuations from the horizontal mean, a T and a NT = 1 − a T are the turbulent and non-turbulent area fractions, respectively, · T and · NT denote the mean inside turbulent and non-turbulent regions, respectively, and b w T = bw T − b T w T and b w NT = bw NT − b NT w NT denote the buoyancy flux within turbulent and non-turbulent regions, respectively. The first term on the right-hand side of Eq. 10 is the turbulent contribution to the total buoyancy flux, the second term is the non-turbulent contribution, and the third term is the contribution resulting from the difference in mean properties between turbulent and non-turbulent regions. Fodor and Mellado (2020) showed that at both Fr 0 = 0 and Fr 0 = 20, turbulent regions contribute the most to the entrainment buoyancy flux, i.e., the term a T b w T dominates in Eq. 10 and hence b w ≈ a T b w T .
As observed in Fig. 4 , the turbulent contribution represents more than 80 − 90% of the total flux for all Froude numbers considered in this study except for Fr 0 = 10, where it is 60 − 70%. Hence, we focus on the term a T b w T and study the two contributions a T and b w T separately.
The changes in b w with Fr 0 observed in Fig. 1 cannot be attributed to the flux itself within turbulent regions, b w T , which is seen in Fig. 5 to behave rather differently with increasing wind shear. In particular, we observe in Fig. 5b a large decrease of the magnitude of b w T with Fr 0 up to Fr 0 ≈ 10 − 15 for all values of the CBL depth. This decrease is more clearly illustrated when plotted in terms of the entrainment-zone Froude number in Fig. 5c, where the flux for different combinations of Froude number and CBL depths approximately collapse in one single curve. This decrease is caused by the decreasing magnitude of the correlation coefficient between buoyancy and vertical velocity fluctuations that comes with increasing wind shear (Fig. 6). We note that the magnitude of the buoyancy flux within turbulent regions increases again after Fr e ≈ 1 mainly due to the increase of the intensity of the buoyancy fluctuations, the vertical velocity fluctuations playing a minor role. This is shown in Figs. 7 and 8, where the root-mean-square (r.m.s.) of the buoyancy is observed to increase 30 − 70% beyond Fr e ≈ 1, whereas the corresponding increase in the vertical velocity r.m.s. is only 10 − 20%.
The variation of the turbulent area fraction with increasing wind shear is shown in Fig. 9. In contrast to the magnitude of the turbulent buoyancy flux, a T increases with increasing Fr 0 until almost the entire area is turbulent. This variation is illustrated in Fig. 10 by means of horizontal cross-sections of the enstrophy field. The increase in turbulent area fraction is more significant than the increase of the buoyancy r.m.s., particularly at Fr 0 ≈ 10 − 15 when a T approximately doubles, whereas b rms increases by 30%. We also note that the increase in a T occurs quite rapidly around Fr 0 ≈ 10 − 15, which is consistent with the rapid increase of the entrainment-zone enstrophy with increasing Fr 0 observed in the p.d.f.s in Fig. 2.
In summary, we have two competing wind-shear effects on the entrainment buoyancy flux: the decrease in the magnitude of the correlation between buoyancy and vertical velocity fluctuations (Fig. 6), and the increase in the turbulent area fraction (Fig. 9). These two competing effects help us understand the behaviour of the entrainment buoyancy flux observed in Fig. 1 and previously reported in the literature Conzemius and Fedorovich 2006;Pino and de Arrellano 2008). The decrease in the magnitude of the entrainment buoyancy flux from Fr 0 = 0 to Fr 0 = 5 up to z enc /L 0 ≈ 20 is due to the turbulent area fraction remaining comparable to that in the shear-free case, but the weak wind shear already reducing the correlation between buoyancy and vertical velocity fluctuations. This is seen more clearly in Figs. 6c and 9c, where the correlation coefficient and turbulent area fraction are compared to the shear-free values. (We note that time increases from right to left in those panels, since ( z i ) c grows faster in time than u and the entrainment-zone Froude number, Fr e , decreases in time.) At early times, for CBL depths z enc /L 0 ≈ 10, the magnitude of the correlation coefficient at Fr 0 = 5 decreases by up to 20% of the shear-free value, whilst at later times, for CBL depths z enc /L 0 ≈ 30, the magnitude of the correlation coefficient becomes approximately equal between the two cases. This convergence of the Fr 0 = 5 case towards the Fr 0 = 0 case is expected, as wind shear is limited by the fixed velocity in the free atmosphere in the barotropic configuration considered herein, whereas the turbulence intensity caused by convection grows with the CBL depth.
The reason for the reduction of the magnitude of the correlation between buoyancy and vertical velocity fluctuations with wind speed can be explained by studying how TKE is generated and transferred between the vertical and the streamwise components. Under shearfree conditions, the kinetic energy is predominantly generated in the vertical direction and transferred to the streamwise direction, as indicated by the positive sign of the horizontal component of the pressure-strain correlation in Fig. 11 b, and the negative sign of the vertical component in Fig. 11c. In contrast, under sheared conditions, the buoyancy is no longer the only forcing that causes vertical velocity fluctuations. Shear production (Fig. 11a) generates kinetic energy in the streamwise direction, which, if strong enough, can change the signs of the pressure-strain correlations and thus invert the direction in which kinetic energy is being transferred between the vertical and the streamwise components. Figure 11b, c indicate that this change occurs at Fr 0 ≈ 5 − 10; for larger values of the Froude number, kinetic energy in the entrainment zone is transferred from the streamwise to the vertical direction, as indicated by the negative sign in the horizontal component of the pressure-strain correlation (Fig. 11b), and the positive sign in the vertical component (Fig. 11c). The transfer of energy from the streamwise to the vertical direction implies that vertical velocity fluctuations start to co-vary with streamwise velocity fluctuations, and thus explains the reduction in the magnitude of the correlation between the buoyancy and vertical velocity.
The decrease in the magnitude of the correlation between buoyancy and vertical velocity fluctuations may also be associated with a change in organization that takes place within the entrainment zone under weak shear conditions. In particular,  discuss how wind shear stretches thermals in the streamwise direction and that this stretching restricts their vertical motion to a shallow layer, resulting in reduced vertical turbulent exchange. However, those authors also concluded that "the kinetic energy of thermals [...] is primarily redistributed into horizontal velocity fluctuations" and our results suggest otherwise: TKE is produced locally by shear production, which is then transferred to the vertical direction by means of the pressure-strain correlation. In addition, we note that under the shear-sheltering hypothesis of Hunt and Durbin (1999), buoyancy and vertical velocity fluctuations should stay the same or possibly decrease in magnitude compared to the shear-free case. Our data lend some support to this hypothesis; at Fr 0 = 5, the buoyancy r.m.s. remains similar to the shear-free case (Fig. 7), whilst the vertical velocity r.m.s. decreases slightly (Fig. 8 ).
Returning to Fig. 1, we observe practically no wind effect on the total buoyancy flux at Fr 0 = 10. For this Froude number, the magnitude of the correlation coefficient between buoyancy and vertical velocity decreases by 30 − 50% of the shear-free value (Fig. 6b, c). The increase in the turbulent area fraction (Fig. 9b, c) and buoyancy r.m.s. (Fig. 7b, c) partly compensates this decrease, although the non-turbulent contributions in Eq. 10 are nonnegligible at Fr 0 = 10 (Fig. 4). For Fr 0 > 10, the magnitude of the correlation coefficient decreases by up to 60% of the shear-free value, whilst the turbulent area fraction more than doubles. The turbulent buoyancy r.m.s. also increases, although not as much as the turbulent area fraction and more gradually with Fr e (compare Figs. 9c, 7c). The increase of the r.m.s. of the vertical velocity component with wind shear is less relevant (Fig. 8c). Thus, the increase in the turbulent area fraction and, to a lesser extent, in the turbulent buoyancy r.m.s. become the main drivers of the strengthening of the entrainment buoyancy flux for the cases considered herein.
The turbulent area fraction increases significantly beyond Fr 0 = 10 because this is when shear production of turbulence becomes the dominant source term in the TKE budget, as shown by . More concretely, for conditions where the entrainment-zone Froude number is larger than Fr e ≈ 1, their study shows that the shear production of TKE in the entrainment zone is larger than the turbulent transport of TKE from the mixed layer into the entrainment zone. Since the non-dimensional variable Fr e was derived based on an integral analysis of the TKE budget equation restricted to the entrainment zone, it is not surprising that it successfully characterizes the behaviour of the turbulent area fraction with increasing wind shear (Fig. 9c). Indeed, panels c in Figs. 5-9 demonstrate that Fr e is capable of characterizing a variety of entrainment-zone properties conditioned into turbulent regions under different wind-shear conditions.
One may also ask what happens for even larger reference Froude numbers than are considered here. Figure 9b suggests that the turbulent area fraction at z i, f reaches a maximum between 0.9 and 1, but Fig. 6b indicates that the magnitude of the correlation coefficient may actually increase, as suggested by its initially larger value at Fr 0 = 25 compared to 10 ≤ Fr 0 ≤ 20. One may expect that, under very strong wind conditions when the dynamical effect of the buoyancy fluctuations introduced by the surface flux and the stratification in the free troposphere becomes negligibly small, the boundary layer approaches neutral conditions and turbulence fluctuations are increasingly dominated by shear instabilities only. In this case, the buoyancy behaves like a passive scalar. Both the buoyancy and vertical velocity then become determined by the shear forcing and follow the same flow organization imposed by the shear instabilities. Subsequently, their correlation increases.

Conclusions
Explanations for the behaviour of the entrainment buoyancy flux under different wind-shear conditions in the CBL have, in the past, been dominated by the viewpoint that changes are caused by the effect of wind shear on the TKE. Here we have extended that view by a conditional analysis of turbulent and non-turbulent regions and by observing the competing effects of wind shear on the correlation between buoyancy and vertical velocity, on the intensity of their fluctuations, and on the turbulent area fraction.
The correlation between buoyancy and vertical velocity in the turbulent regions varies non-monotonically and is determined by the dominant forcing inside the entrainment zone. In the absence of flow, buoyancy is the only external forcing and the organization is dominated by convective plumes. Kinetic energy is transferred to the horizontal direction by means of the pressure-strain correlation and the horizontal velocity is determined by the forcing in the vertical direction. However, when wind is introduced as an additional external forcing, shear production generates kinetic energy in the streamwise direction, which starts to compete with the energy received from the buoyancy forcing. Vertical velocity fluctuations become partly determined by the redistribution of energy from the streamwise direction to the vertical direction, thus reducing the magnitude of the correlation between the buoyancy and vertical velocity. Under very strong wind conditions, the buoyancy virtually becomes a passive scalar. Both the buoyancy and vertical velocity become determined by the shear forcing and follow the same flow organization imposed by the shear instabilities. This then increases the magnitude of their correlation again.
The turbulent area fraction, on the other hand, increases with wind shear until almost the entire area at the height of minimum buoyancy flux is turbulent. This is determined by the TKE budget. In shear-free conditions, the buoyancy flux term is a sink of TKE in the entrainment zone and turbulence is not generated there, it is only transported from the mixed layer. This causes turbulence to rapidly decay with height within the entrainment zone. As the wind speed is increased, shear production of TKE increases, generating turbulence inside the entrainment zone and thus increasing the turbulent area fraction. The strength of buoyancy fluctuations in turbulent regions also increases substantially with wind shear, but less so and more gradually than the turbulent area fraction. The magnitude of the vertical velocity fluctuations increases only moderately.
These findings may also have implications for the parametrization of entrainment fluxes other than the buoyancy flux. In particular, they suggest obtaining parametrizations in terms of the turbulent area fraction and the flux inside the turbulent regions. Here we have shown that the entrainment-zone Froude number Fr e ≡ u/[N 0 ( z i ) c ] captures the dependence of these properties on both buoyancy and wind external forcings, thus enabling entrainmentzone properties to be expressed as a function of easily measured variables: the velocity jump between the mixed layer and the free atmosphere, u, the buoyancy frequency in the free atmosphere, N 0 , and the mixed-layer depth, z enc , through the shear-free entrainment-zone thickness ( z i ) c ≡ 0.25 z enc .