On the departure from Monin–Obukhov surface similarity and transition to the convective mixed layer

Large-eddy simulations are used to evaluate mean proﬁle similarity in the convective boundary layer (CBL). Particular care is taken regarding the grid sensitivity of the proﬁles and the mitigation of inertial oscillations in the simulation spin-up. The nondimensional gradients φ for wind speed and air temperature generally align with Monin–Obukhov similarity across cases but have a steeper slope than predicted within each proﬁle. The same trend has been noted in several other recent studies. The Businger-Dyer relations are modiﬁed here with an exponential cutoﬀ term to account for the decay in φ to ﬁrst-order approximation, yielding improved similarity from approximately 0.05 z i to above 0.3 z i , where z i is the CBL depth. The necessity for the exponential correction is attributed to an extended transition from surface scaling to zero gradient in the mixed layer, where the departure from Monin–Obukhov similarity is negligible at the surface but becomes substantial well below the conventional surface layer height of 0.1 z i .


Introduction
Within the atmospheric boundary layer (ABL), the surface layer is unsurprisingly the region directly above the Earth's surface.This layer is often described in terms of its properties -approximately constant flux, negligible rotation effects, and adherence to surface scaling -rather than formally defined (Sutton 1953;Kaimal and Finnigan 1994).One common convention is to assume the surface layer extends to the lowest 10% or so of the ABL (Stull 1988;Garratt 1994), consistent with the depth of the logarithmic (log) region in more general wallbounded flows (Pope 2000).
With respect to the surface scaling property, the scaling of mean flow statistics within the surface layer is given by Monin-Obukhov similarity theory (MOST, Monin and Obukhov 1954;Foken 2006) that extends the log law of the wall (von Kármán 1930;Prandtl 1932;Millikan 1938) to thermally stratified conditions.MOST predicts universal similarity for the nondimensional mean gradients: where the gradients are normalized using the log law definition and functions for φ must be determined empirically.Here, z is height above the surface, U (z) is the mean horizontal wind speed, θ(z) is the mean potential temperature, u * is the surface friction velocity scale, θ * is the surface temperature scale, and κ is the von Kármán constant.The velocity and temperature scales are related to the surface momentum flux as u ′ w ′ s = −u 2 * and surface heat flux as w ′ θ ′ s = −u * θ * such that Eq. ( 1) is often referred to as the flux-gradient relations.Owing to the assumed absence of other length, velocity, or temperature scales in the theory, φ m and φ h are considered functions only of ζ = z/L defined using the Obukhov (1946) length where θ s is the mean surface temperature and g is the gravitational constant.
Following the introduction of MOST, evaluations of field measurements from meteorological towers have largely corroborated the surface layer theory and universality of φ m (ζ) and φ h (ζ).For a convective ABL (CBL) with L < 0 typical of daytime conditions, several experimental campaigns and reevaluations proposed power-law relations for φ m and φ h with some variability in the fitted parameters but a consistent general form for the functions (see, e.g., Dyer and Hicks 1970;Businger et al. 1971;Carl et al. 1973;Yaglom 1977;Högström 1988;Wilson 2001;Katul et al. 2011).The most common of these empirical relations are the Businger-Dyer profiles for convective conditions (Businger et al. 1971;Dyer 1974) where the values of the constants depend modestly on the analysis.Agreement of field measurements with MOST and the empirical relations can be further improved by accounting for additional effects not considered in the idealized theory, e.g.time-dependent variability from large-scale turbulence (Salesky and Anderson 2020) and anisotropy due to complex conditions (Stiperski and Calaf 2023).More recently, direct numerical simulations and large-eddy simulations (LES) of the CBL have produced consistent trends that support MOST to first-order approximation but also reveal possible shortcomings: gradient statistics align with the Businger-Dyer profiles when comparing results across different simulated conditions (Maronga and Reuder 2017), but the decay in φ(ζ) within each individual profile is steeper than predicted by Eq. ( 3), demonstrating a lack of universality in φ(ζ) (Khanna and Brasseur 1997;Pirozzoli et al. 2017;Li et al. 2018).Accordingly, it has been proposed that surface layer gradients -particularly for velocity -may additionally depend on the boundary layer depth due to the influence of large-scale motions from the well-mixed layer that forms the bulk of the CBL (Khanna and Brasseur 1997;Johansson et al. 2001).This idea is indirectly supported by observed trends in field measurements that suggest a parameter space beyond ζ is required to account for variability in gradient statistics (Salesky and Chamecki 2012).
The connection between turbulent eddies in the mixed layer and surface layer statistics has been substantiated by several analyses of turbulent structure in the CBL.Conditional statistics show that the steep decay in φ(ζ) and deviations from MOST noted above are predominately associated with large-scale turbulent events such as downdrafts (Li et al. 2018;Fodor et al. 2019).In more general terms, these deviations are related to the modulation of near-surface turbulence by buoyancydriven eddies from aloft (Lemone 1976;Smedman et al. 2007;Gao et al. 2016;Salesky and Anderson 2018;Liu et al. 2019;Dupont and Patton 2022).A signature of the boundary layer depth also appears in the velocity and temperature spectra within and above the surface layer (McNaughton et al. 2007;Chowdhuri et al. 2019).
There have been relatively few attempts to model the modulation of surface layer gradients.Santoso and Stull (1998) modeled mean wind profiles with a power law and exponential decay to achieve a smooth transition from the surface to the uniform mixed layer.Gryning et al. (2007) combined surface scaling with a constant mixed layer length scale, but the foremost goal was to extend similarity to higher positions above the surface layer.Salesky and Anderson (2020) corrected Eq. ( 3) for local-in-time deviations due to large-scale fluctuations.Li et al. (2021) quantified the deviation as a nonlocal transport through the framework of eddy diffusivity models.Cheng et al. (2021) and Liu et al. (2022) both introduced a correction to φ prescribed as a function of z i /L, where z i is the base height of the stable capping inversion and is typically used to define the CBL depth.
One cautionary note regarding the previous findings and models is that several of the studies used simulations confined to relatively low Reynolds number compared to the ABL (Pirozzoli et al. 2017;Li et al. 2018;Fodor et al. 2019;Cheng et al. 2021).Considering the log law only emerges for high Reynolds numbers (Marusic et al. 2013;Sillero et al. 2013;Lee and Moser 2015), the results may reflect a combination of buoyancy effects and finite Reynolds number corrections to the log law.Additionally, for wall-modeled LES studies the grid convergence of surface layer statistics is often not scrutinized.The presence of these limitations precludes a careful quantitative comparison of deviations from MOST observed across the literature.
In the present work, recurring trends in φ(ζ) observed from simulations are further evaluated using new LES of the idealized dry CBL.In consideration of the effects noted above, the LES is for the inviscid limit and includes a detailed test of grid sensitivity.To account for the observed trends, the dependencies of φ m and φ h in Eq. ( 1) are expanded to include z i following the suggestion of Khanna and Brasseur (1997).The approach also considers recent evidence of outer layer stratification influencing surface similarity under stable conditions (Heisel and Chamecki 2023), except in this case the similarity is influenced by free convection in the mixed layer.The goal is to empirically explain the behavior of the simulated mean profile statistics in the broader context of the transition from the surface layer to free convection in the mixed layer, while also reconciling the simulation results with the widespread support for MOST from field experiments discussed above.The proposed explanation -an approximately exponential decay in φ as a function of z/z i -has a limited effect on statistics very close to the surface where many field measurements are acquired, is qualitatively consistent with the profile shapes seen in recent simulation studies, and accounts for the profile transition between the surface and mixed layers with reasonable accuracy.The remainder of the article is organized as follows: the new LES cases are described in Sect.2; mean profile similarity is assessed in Sect.3; implications for resistance laws in the mixed layer and for the definition of the surface layer are discussed in Sect.4; finally, a summary is given in Sect. 5.

Large-eddy simulations
The present simulations were conducted using standard practices for representing an idealized dry convective ABL (see, e.g., Deardorff 1972;Moeng and Sullivan 1994;Sullivan et al. 1994;Noh et al. 2003;Salesky et al. 2017): a range of unstable conditions was achieved by imposing different combinations of fixed surface heat flux Q * = u ′ w ′ s and geostrophic wind speed U g , and the boundary layer was confined by a stable capping inversion.The inversion was introduced through the initial temperature profile following Sullivan and Patton (2011), which includes a uniform temperature in the boundary layer, a strong lapse rate Γ = 0.08 K m −1 in the range z = 1000 to 1100 meters, and a weaker lapse rate Γ = 0.003 K m −1 to form the top-most capping inversion.
Additional imposed parameters include the aerodynamic roughness length z o = 0.1 m and Coriolis frequency f = 1 × 10 −4 s −1 .Six primary cases with varying Q * and U g and the resulting scaling parameters are summarized in Table 1.A seventh case (E1) is the same as case E, but with the initial capping inversion positioned 200 m lower as seen in the z i values that were determined from the height of the minimum heat flux.The range of simulated conditions span from relatively weak (−z i /L = 2.5) to moderately strong (−z i /L = 39) convection.All cases employed a numerical grid with 1024×1024×512 points and corresponding domain dimensions of 12×12×2 km.
Dimensional profiles of the mean horizontal wind speed and potential temperature are respectively shown in Figs.1a and 1b for the seven simulated cases.Every case exhibits an extensive mixed layer with approximately uniform wind speed and air temperature, as well as an entrainment layer centered around z i and an overlying temperature inversion.The temperature in the mixed layer increases with Q * and is highest for case E1 with a shallower CBL that can be heated more quickly.Hereafter, each case is referred to using its alphabetical label A-F indicated in the figure legends alongside the z i /L values.Unless otherwise noted, "velocity" refers to the magnitude of the horizontal components U x and U y , and "temperature" refers to the potential temperature.Further details on the numerical code and procedure used to generate the Fig. 1 profiles are given below.x + U 2 y ) 1/2 (a,c) and potential temperature θ (b,d).Rows correspond to dimensional height (a,b) and relative height z/z i with logarithmic spacing (c,d) to show the near-surface trends.In this and later figures, the legend corresponds to the z i /L stability parameter and alphabetical label for the cases in Table 1.The temperature profiles in cases C, D, E, and F overlap in the outer layer due to the matching Q * and initial z i .
The LES flow solver and numerics are a modified version of the code by Albertson and co-authors (Albertson 1996;Albertson and Parlange 1999;Porté-Agel et al. 2000).The solver utilizes a staggered vertical grid for the vertical velocity component to optimize the vertical differentiation with second-order-accurate finite differencing.Horizontal derivatives are computed spectrally with full de-aliasing using zero padding.The flow is evolved through time integration using the second-order-accurate Adams-Bashforth method.The boundary conditions are periodic in the horizontal directions and zero penetration along the bottom and top of the domain, with the top additionally employing stress-free conditions.Fluctuations in the top 25% of the domain are damped to achieve the upper boundary conditions and mitigate gravity waves (Nieuwstadt et al. 1993).Subgrid-scale (SGS) stresses are estimated using a Lagrangian averaged scale-dependent model to represent the SGS eddy viscosity (Bou-Zeid et al. 2005).The SGS heat flux is estimated using the same value as the local SGS eddy viscosity along with a constant Prandtl number Pr sgs = 0.4.A more detailed description is given elsewhere (Kumar et al. 2006;Salesky et al. 2017).
The surface stress in the wall-modeled LES is estimated locally in space and time using MOST and Eq.(3) for momentum (Moeng 1984).The wall model for u * is subject to so-called overshoot and log-law mismatch (Mason and Thomson 1992;Brasseur and Wei 2010;Larsson et al. 2015).To mitigate this effect, the wall model was evaluated using velocity values at 0.05z i rather than the first grid point (Kawai and Larsson 2012), where z i in this case is based on the initial temperature profile and for simplicity does not change with time.Based on a comparative test, it was found that using values farther from the surface reduces the wall model overshoot but does not noticeably improve the convergence discussed below.
The final methodological considerations are the procedure to spin up the simulations from initial conditions and the sensitivity of the results to grid resolution.Both of these aspects are crucial to the repeatability and validity of the study and an extended discussion is given for each in the two subsections below.Statistical uncertainty is discussed later in Sect. 3 in the context of the results.

Simulation spin-up and inertial oscillations
A common approach for simulating the CBL is to impose an initially geostrophic velocity field with no boundary layer, i.e. with U x (z) = U g and U y (z) = 0.The simulations are then spun up for a duration in the range of tw * /z i ≈ 5 to 30 turnover times until statistics appear stationary (see, e.g., Moeng and Sullivan 1994;Noh et al. 2003;Salesky et al. 2017;Maronga and Reuder 2017;Li et al. 2018;Liu et al. 2023, among others), where z i and the convective velocity w * (Deardorff 1970) represent properties of the large-scale eddies that govern dynamics in the bulk of the CBL.While the reasoning in this approach is sound, it neglects the presence of inertial oscillations resulting from the initial velocity field (e.g., Shibuya et al. 2014;Momen and Bou-Zeid 2017), noting the oscillations are not relevant to free convection cases with no geostrophic wind.
When the U x and U y fields differ significantly from the quasi-equilibrium condition, an oscillatory response arises from the Coriolis force and the oscillations dampen in time due to the surface stress (Schröter et al. 2013).An example is shown in Fig. 2, where the trivial initial profile in 2a (spin-up 1, yellow) results in an oscillating time evolution of average surface shear stress τ s (2b), mixed layer wind speed U m (2c), and stability z i /L (2d).The same time evolution occurs during spin-up of the conventionally neutral ABL (Pedersen et al. 2014;Liu et al. 2021).
For the present example, case A in Table 1 is employed on a coarser 200×200×100 numerical grid with f = 1.4 × 10 −4 s −1 increased to a typical polar value and the initial inversion height reduced by 200 m in the same manner as case E1.The latter two changes were necessary to prevent the CBL from reaching the damping layer before completing an inertial period 2π/f .Parameters with overbars ( • ) indicate the long-term average value over the full period.
Ending the simulations near tw * /z i ≈ 20 when the statistics are momentarily stationary would lead to shear stress and wind speed values that are significantly out of equilibrium with the conditions given by Q * , z o , U g , f , and the initial z i .Of particular importance to the present analysis is the Obukhov length L which varies by approximately 50% within the first half oscillation for the example in Fig. 2d.
The amplitude of inertial oscillations can be reduced by using initial U x and U y velocity fields that are closer to the quasi-equilibrium.To this end, an ad hoc procedure using multiple spin-up trials was developed to determine appropriate initial conditions for the final simulations.For the initial velocity profiles of the second spin-up, the mean conditions at tf = π in the first spin-up were rescaled along z to match the original z i as shown in the inset of Fig. 2a.Additionally, the original initial temperature profile (Sullivan and Patton 2011) was used rather than a rescaled temperature profile in order to restore the two initial inversion layers.The profile re-scaling and new spin-ups can be repeated as necessary until the initial condition is close to the equilibrium.In Fig. 2, the oscillations are significantly reduced for spin-up 2, and there are limited changes in the initial profiles and time evolution between spin-ups 3 and 4, suggesting the proper initial velocity field has been reached.It is also important to note that the wind speed and scaling parameters slowly increase over time in Figs.2b-d in response to the growth of the CBL through entrainment.
One consequence of the long spin-up trials is that a low-level jet can develop within the entrainment layer directly above the heat flux minimum, particularly for weaker convection.This feature is not apparent from U x and U y for the profiles at tf = π in the inset of Fig. 2a (yellow lines), and is more evident from the horizontal wind magnitude U .The jet is excluded from the rescaled profiles by assuming a linear trend in the top 100 m of U x and U y as seen in Fig. 2a.The jet is far from the region of interest below the mixed layer and understanding its emergence is outside the scope of the present study.
While Fig. 2 is predominately for demonstration purposes, the same spin-up procedure was employed for the cases in Table 1 to approximately determine the quasi-equilibrium wind profiles.For each simulated condition, four spin-up trials were completed on a grid with 200×200×100 points: the first with initial conditions U x = U g and U y = 0, and the subsequent trials using the rescaled profiles as described above.The rescaled profiles resulting from the fourth trial were then used as initial conditions for the final simulations.It may be possible to expedite the determination of quasi-equilibrium conditions using existing methods not employed here.These options include initializing U x and U y based on an approximation of the cross-isobar angle and generating profiles using a single-column model (Ghan et al. 2000).Inertial oscillations in the latter option can be reduced by separating the equations based on the geostrophic "background" flow and the boundary layer deviations (Gutman 1973).
The final simulations were spun-up on a larger grid with 400×400×200 points for approximately tw * /z i ≈ 15 turnover times.This duration ranges from 140 to 180 physical minutes for the different cases and is long enough for the flux profiles to develop, but short enough to avoid interference of the damping layer on the growing CBL.Because the initial U x and U y profiles are close to the equilibrium in the final spin-up, the inertial oscillations are minimized and it is not necessary to simulate a full inertial period.
After the spin-up, the velocity and temperature fields were interpolated onto the final grid with 1024×1024×512 points, and the simulations were continued for an additional 70 physical minutes.A short period is required for small-scale turbulence to develop within the finer grid, such that the first 10 minutes are excluded from the analysis.The time-averaged statistics for each case were computed across the last 60 minutes of the simulations.The time-averaged statistics are featured in Fig. 1 and all later results.

Grid sensitivity of near-surface statistics
For wall-modeled LES, flow statistics near the surface can be significantly biased by the wall model and grid resolution.To avoid inclusion of such biases in the analysis, the present section assesses the specific range of heights near the surface where the mean and gradient statistics are approximately converged with respect Table 2 Numerical grids used to test the resolution convergence of flow statistics for cases A, C, and E in table 1. Grid properties include the number of nodes N j along each direction j = x, y, and z, the domain size L j , corresponding grid resolution ∆ j , and effective resolution ∆ = (∆x∆y∆z) 1/3 relative to the inversion height for case A. Results reported in the later results are based on the finest resolution.  1 were repeated for the series of grid sizes summarized in Table 2.For a given case, all resolutions used the same initial velocity profiles and spin-up as determined from Sect.2.1, with the different resolution introduced for the final 70 minutes of simulation.The sensitivity analysis only directly uses the two finest grids, but the additional coarser grids are useful for identify general trends discussed later in the context of the results.
Mean profiles of wind speed and air temperature are shown in Fig. 3 for all tested grid sizes, with a vertical discplacement between different cases for visualization.The profiles are plotted as the diabatic term (Panofsky 1963) that quantifies the difference between the local mean value and the log law, where ψ increases with convection.The value used here for the von Kármán constant is κ = 0.39 (Marusic et al. 2013).
To more directly compare the profiles across resolution, the normalization in Fig. 3 uses fixed scaling parameters u * and θ * from the finest grid size rather than the parameter values resulting from each individual resolution.While there is an order of magnitude increase in resolution from the coarsest to finest grids in Table 2, the corresponding increase in u * is only 6.0% for case A and 2.6% for case E. The difference in u * between the two largest grids is 0.14% for case A and 0.15% for case E.
The criterion used to determine the converged region of the mean profiles is the percent difference in ψ between the grids with 800×800×400 points and 1024×1024×512.Specifically, the difference in ψ must be less than 1% for this 28% increase in resolution.The vertical lines in Fig. 3 indicate the lowest height where this criterion is met, which varies from 0.04z i (case E) to 0.06z i (case A).
The converged regions for cases B and D were inferred through interpolation, and the lower height in meters for case E was directly used for cases E1 and F. Later figures either exclude statistics in the near-surface region where convergence is not observed or clearly differentiate the unconverged results.Finally, while it is not apparent from Fig. 3, the mean wind speed and air temperature in the mixed layer are converged for the 600×600×300 grid in case A, and at coarser resolutions for stronger convection.
A similar assessment is made for the nondimensional gradient profiles φ(z) in Fig. 4. The general concave shape of the profiles at the surface matches closely with previous observations (e.g., Bou-Zeid et al. 2005;Maronga and Reuder 2017) and demonstrates that the influence of the wall model, SGS model, and near-surface resolution extends well beyond the first grid points.The criterion used here for convergence of the gradient statistics is that the difference in φ must be less than 0.01 between the grids with 800×800×400 points and 1024×1024×512.Using the same 1% threshold as above was found to be overly conservative given the very small magnitude of the derivatives far from the surface and would result in the defined converged regions beginning at moderately higher positions.The trends beginning near 0.8z i in Figs. 3 and 4 correspond to the start of the entrainment layer that is outside the scope of the present similarity analysis.Under weak convection in case A, the gradient statistics do not converge within the traditional surface layer below 0.1z i .There is also some uncertainty in the convergence across the lower half of the CBL, as the slope in φ appears to decrease with increasing resolution across the region with positive gradient.In moderate convection, both the upper portion of the surface layer and the outer layer exhibit grid convergence for the relatively fine resolution employed here.Figure 4 demonstrates the challenge in using wall-modeled LES to critically evaluate near-surface behavior with a high degree of certainty.Accordingly, the conclusions drawn in the present work are confined to robust trends that extend beyond the surface layer and across all cases.The general findings are not contingent on the weakly convective cases that may not have reached mesh independence above the surface layer.With the threshold of 0.01, the converged regions of φ begin approximately 20 to 30 points away from the surface.The exact number of points is expected to depend on numerous variables including the LES wall and SGS models, grid size, grid aspect ratio, and convergence criterion, such that the quantitative outcomes of Figs. 3 and 4 are considered to be specific to this study.

Mean velocity and temperature similarity
The profiles in Fig. 1, generated using the spin-up procedure outlined in Sect.2.1, exhibit convergence with respect to grid resolution in the upper portion of the surface layer as discussed in Sect.2.2.Similarity in the region with converged statistics, including in the convective matching layer above the surface layer, is now evaluated in further detail.The upper half of the CBL is excluded from the results due to weak top-down effects from the entrainment layer that yield negative φ values as seen in Fig. 4.

Comparison with Businger-Dyer relations
The diabatic mean wind speed and air temperature profiles predicted by the Businger-Dyer relations in Eq. ( 3) result directly from the integration in Eq. ( 4) (Paulson 1970): Here, x = φ −1 m = (1 − 16ζ) 0.25 and the contribution ψ(−z o /L) from the lower limit of the integral is neglected under the condition −z o /L ≪ 1.In Figs.5a and  5b, the relations in Eq. ( 5) are compared with ψ computed directly from the LES profiles using Eq. ( 4).The range of heights included for each profile in Figs.5a and 5b spans from the bottom of the converged region identified from Fig. 3 up to 0.3z i .The upper limit extends beyond the traditional definition for the surface layer and is included for consistency with the later analysis.The height 0.1z i is indicated in all figures where necessary to facilitate the distinction of trends within and above the traditional surface layer.While Eq. ( 5) appears to conform to the general shape of each LES profile, there is a distinct stability trend in which the profiles become increasingly offset from the prediction with weakening convection.
The displacement between cases is amplified in the gradient φ profiles shown in Figs.5c and 5d.Consistent with previous observations, the general trend across cases appears to follow a curve similar to the Businger-Dyer profiles, but the decay in φ along each individual profile is significantly steeper (Maronga and Reuder 2017;Pirozzoli et al. 2017;Li et al. 2018).Unlike the referenced studies, however, the entirety of each φ curve is below the Businger-Dyer profiles.This difference is likely due to a combination of excluding the near-surface points from the present cases and the significant effect of inertial oscillations on the relevant normalization parameters u * and L as seen in Fig. 2. The shaded regions in Figs.5c and 5d represent 95% confidence intervals for the mean gradient statistics.The intervals are estimated from the formula t 95 σ φ / √ N , where σ φ is the standard deviation of observed gradients over the 60-minute simulation period, N is the number of independent observations, and t 95 is the Student's t value that depends on the confidence level and N .The independent observations N vary between 10 (moderate convection) and 20 (weak convection) across cases based on integral time scales computed from turbulent spectra.The steep decay in φ noted above and the separation of the curves in Fig. 5 both exceed the statistical uncertainty.
Figure 5 demonstrates that MOST and ζ alone cannot fully account for the LES profile trends.Scaling adjustments to the definitions of ζ and/or φ are required for the profiles within the surface layer to collapse along a common curve within the uncertainty bounds.At the same time, alignment with Businger-Dyer relations in field experiments and across different cases in simulation studies suggest that Eq.
(3) provides a reasonably accurate foundation for similarity in the mean profiles.In the following section, the trends in φ are considered in the context of the extended profile to identify a possible similarity framework that is compatible with these past and present findings.

Gradient profile trends
The nondimensional φ profiles in Figs.5c and 5d are replotted in Fig. 6 as a function of relative position z/z i within the CBL.The axes are shown with both logarithmic (top) and log-linear (bottom) scaling to facilitate interpretation of the profile shapes in different regions.The transparent extension of each curve is the region that did not meet the resolution convergence criterion detailed in Sect.2.2 and Fig. 4. The velocity gradients approaching the mixed layer exhibit incomplete statistical convergence, leading to modest random error along the profiles.Recalling that the observed values φ O(0.1) are the product of the gradient and the height z ∼ O(100 m) in Eq. ( 1), the dimensional gradient is exceedingly small in this region.Additional computing resources are not currently available to continue the simulations for a longer duration.Conclusions made from the gradient statistics are limited to trends that exceed the observed variability within each profile.
For the region near and below 0.1z i in Figs.6a and 6b, there is noticeable curvature in the profile for weak convection in case A, but φ(z) increasingly resembles a power law with increasing convection.Further, the slope of the profiles, i.e. the exponent of the power law, does not vary across the tested cases.These trends are consistent with the form of the Businger-Dyer profiles in Eq. ( 3), where the contribution of 1 vanishes with increasing convection and the power law exponent is assumed constant.With respect to grid resolution, the approximate power law in φ(z) under moderate convection is only apparent for the two largest grids tested in Table 2.
At higher positions, the mean gradients in Figs.6a and 6b appear to be governed by an sharp decay in φ that results in the cutoff φ(z=0.4zi ) ≈ 0. The position of the cutoff is approximately constant as a fraction of z i , whereas the start of the cutoff seen in Figs.5c and 5d unambiguously varies with L. This has important implications for the approach to free convection and −ζ → ∞: rather than the mixed layer starting at lower positions and free convection reaching closer to the surface with decreasing L, Fig. 6 suggests the near-surface region maintains a fixed height and φ decreases with increasing convection until it vanishes in the free convection limit.
The shape of φ approaching the mixed layer in Figs.6c and 6d is approximately linear, particularly for the temperature.The linear trend implies a logarithmic decay in φ and a dimensional gradient that resembles − log(z)/z.The height where the logarithmic decay gains leading order importance over the approximate power law appears to depend on stability and emerges at lower heights for the weaker convection cases.Importantly, the decay observed here is specific to a CBL with a well-defined mixed layer as seen in Fig. 1.For the tested grid sizes in Table 2, the logarithmic trend begins to appear with an overestimated slope for the grid with 400×400×200 points, and the slope is approximately converged for the grid with 600×600×300 points.
The trends in Fig. 6, i.e. the consistency with Eq. ( 3) closer to the surface and the logarithmic decay approaching the mixed layer height, show the potential to augment the existing Businger-Dyer profiles with an additional term that enforces the decay in φ with increasing z/z i .In this sense, the term is a correction for the transition between the surface and mixed layers, where previous findings (e.g., Salesky and Chamecki 2012;Pirozzoli et al. 2017;Li et al. 2018) and Fig. 5 indicate the correction is necessary even within the traditional surface layer below 0.1z i .The existing model corrections for φ discussed in the introduction do not account for the specific trends observed in Fig. 6.For instance, the corrections based on z i /L displace φ by a constant value for a given case and do not consider the shape of the cutoff (Cheng et al. 2021;Liu et al. 2022).The cutoff is also not well described by an inverse summation of length scales (Gryning et al. 2007).A preliminary effort to empirically model the gradient cutoff is given in the following section.

Preliminary model for extended similarity
Deeper within the ABL in the roughness sublayer (RSL), the mean profiles are influenced by complex turbulent drag and mixing interactions associated with the local surface roughness.In this sublayer, it is standard to correct for the mean similarity in a multiplicative manner as φ m (ζ)ϕ m (z/z * ), where φ m accounts for atmospheric stability and ϕ m is a correction based on the relative position within the RSL depth z * (e.g., Garratt 1980;Cellier and Brunet 1992;Mölder et al. 1999).The most common functional form for ϕ is an exponential that models the decay in the gradient within the RSL as the surface is approached (Garratt 1980;Harman and Finnigan 2007;Mo et al. 2022).The same form is adopted here, except the purpose of the exponential in this case is to ensure the gradient de-creases towards zero approaching the mixed layer rather than within the RSL.Santoso andStull (1998, 2001) also used an exponential cutoff correction to the mean profiles to enforce a smooth transition to uniform mixed layer values.The revised similarity relations are given as where a, b, and c are fitted constants and the leading constant a h for temperature accounts for the turbulent Prandtl number.The exponent values in Eq. ( 3) are adopted here, noting that testing a wide range of exponents resulted in values within ±0.05 of 0.25 and 0.5 and only a nominal increase in the coefficient of determination R 2 for the fit.The exponential is included in the definition of φ in Eq. ( 6) rather than as a separate ϕ function because it is part of the same stability correction.The constant c determines the height relative to z i where the exponential becomes small.A value c > 2 is expected such that the cutoff function decreases to small values within the lower half of the CBL.
To evaluate the applicability of the revised similarity relations, the expressions for φ in Eq. ( 6) and their integral ψ defined in Eq. ( 4) were fitted to the LES profiles.The integral for ψ was computed numerically in the absence of a simple analytical solution.The reason for including ψ in the fitting procedure is to assess the extension of φ down to z = z o .While the near-surface region cannot be fitted directly, Eq. ( 6) must have the correct cumulative magnitude below the fitted region in order to align with the LES profiles for ψ.
The cost function for the nonlinear fitting algorithm was the total residual between the predicted and observed φ and ψ values compiled for all cases simultaneously.The fit result therefore represents the range of convective conditions rather than any individual case.The fitting procedure was conducted separately for the velocity and temperature statistics.Due to the complexity of the equations and the use of numerical integration, the algorithm was unable to converge when multiple parameters were undefined.Accordingly, the fit was designed to optimize c with a and b as prescribed inputs, and was repeated for a range of a and b values.The values presented here are those with the highest resulting R 2 , with a m = 1 assumed for velocity.As noted above, the power law exponents in Eq. ( 6) were also varied before the traditional values were selected for simplicity.Finally, heights up to 0.3z i were included in the fits under the assumption that Eq. ( 6) approximates the transition across an extended range up to the convective mixed layer.
The preliminary values resulting from the fit to the velocity profiles are b m ≈ 22, c m ≈ 3.7, and R 2 = 0.974.The values for the temperature profiles are a h ≈ 0.93, b h ≈ 14, c h ≈ 2.9, and R 2 = 0.992.The higher R 2 for temperature is likely due in part to the additional fitted parameter a h and the better statistical convergence of φ in Fig. 6 relative to velocity.Owing to the lack of near-surface points, the present values for a h and b are not suggested as replacements to existing values.Further, the difference between c m ≈ 3.7 and c h ≈ 2.9 is not given a physical interpretation here.It may result simply from the fact that the smaller Businger-Dyer exponent −0.25 for momentum requires a larger cutoff correction to have similar φ near the mixed layer.The main outcome of the fit is to demonstrate the systematic improvement of the profile prediction with the inclusion of a z/z i cutoff.
Figure 7 compares the mean profiles predicted from the integral of Eq. (6) (dashed lines) with the LES profiles.The dotted lines in Figs.7a and 7b are Eq.( 6) with the fitted values, but excluding the exponential cutoff.When plotted as ψ, these dotted lines collapse along the solid black lines in Figs.7c and 7d defined by the integral of the equations given in the legend.
The expression with the exponential cutoff leads to significant improvements in the alignment with the LES profiles in Figs.7a and 7b across all heights within the converged region, particularly for heights above 0.1z i .While there are some discrepancies, most notably for cases A and F, the predictions resulting from Eq. ( 6) provide a reasonable approximation of the mean profiles for an extended range from below 0.1z i to above 0.3z i and near the start of the convective mixed layer.
The diabatic term ψ in Figs.7c and 7d demonstrates the departure of the mean profiles from surface layer similarity as a result of the decay in the gradients along z/z i .The dashed lines representing Eq. ( 6) all begin at ψ(z=z o ) = 0 and become increasingly dissimilar from Monin-Obukhov scaling similarity with increasing z/z i , i.e. the ψ values spread farther apart.This dissimilarity is well predicted by the exponential cutoff correction.Figure 8 evaluates the gradient profiles predicted from Eq. (6) (dashed lines) in the same manner as the mean values in Fig. 7.As before, the dotted lines in Figs.8a and 8b exclude the exponential cutoff and are equivalent to the solid black lines in Figs.8c-f defined by the equations in the legends.
Equation ( 6) matches closely with the φ curves in Figs.8a and 8b compared to the Businger-Dyer profiles without a cutoff correction.However, the exponential cutoff does not fully account for the entirety of the gradient decay, as seen in the deviations from the LES profiles that emerge between 0.2-0.3zi .As discussed previously, the dimensional value of the gradients in this range of heights is very small such that the discrepancy may be of limited practical importance.For instance, the discrepancy in φ approaching the mixed layer is not readily seen in the Figs.7a and 7b mean profiles.
The plots of φ(ζ) in Figs.8c and 8d show the departure from Monin-Obukhov similarity in the LES profiles compared with the prediction from the exponential cutoff.The exponential cutoff closely approximates the decay in the gradients within and above the surface layer for the fitted cases.To collapse the nondimensional gradients along a single curve, it is necessary to group the gradient and cutoff correction as φ exp (cz/z i ).The product, shown in Figs.8e and 8f, now aligns reasonably well with the Businger-Dyer profiles.Most of the residual differences occur near 0.3z i and are due to the discrepancies seen in Figs.8a and 8b and discussed above.
Figure 8 provides promising evidence for expanding the similarity parameter space to include z/z i .The correction in Eq. ( 6) provides profile predictions for an extended range approaching the mixed layer and may account for the departures from MOST observed in previous simulation studies (Khanna and Brasseur 1997;Pirozzoli et al. 2017;Li et al. 2018).However, the generality of the results remain unproven at this point and further comparison with other simulations and measurements is required.In particular, the exponential form of the correction may be specific to the use of the Businger-Dyer profiles to functionally represent MOST.Alternate empirical relations with a steeper decay for strong convection with large ζ would require a different correction.However, the necessity to expand the parameter space to achieve universality of profiles in Fig. 8 is independent of the empirical MOST function used.
If Eq. ( 6) is applicable beyond the present cases, one interesting note is that the cutoff correction is independent of stability.The LES cases in Table 1 span the transition from relatively weak convection with thermal rolls to moderately strong convection with cells (Etling and Brown 1993;Atkinson and Zhang 1996;Khanna and Brasseur 1998;Salesky et al. 2017).For the current results, the correction at a given height z/z i is the same regardless of the convective regime, indicating that the different large-scale structures (i.e.rolls or cells) impinging on the surface layer reduce the average gradient by the same fraction.

Implications for mixed layer resistance dependencies
While there is extensive theory for resistance laws in the geostrophic drag and heat transfer across the entire ABL (see, e.g.Monin 1970;Yamada 1976;Arya 1977), there are relatively fewer studies relating mean wind speed U m and temperature θ m in the convective mixed layer to surface properties.The predicted scaling for U m /u * and (θ m − θ s )/θ * depends on underlying assumptions, in particular regarding the bottom height z m of the mixed layer.The dependencies of the mean values can be inferred by evaluating the mean profiles at the base of the mixed layer z = z m : Equation ( 7) is a direct outcome of the definition for the diabatic term in Eq. ( 4), noting that traditional approaches do not include the z/z i correction for ψ.If z m ∝ z i is assumed in Eq. ( 7), the resulting mixed layer values depend directly on both log (z i /z o ) and ψ(z i /L) (Garratt et al. 1982).Alternatively, if the z/z i correction is excluded and z m ∝ −L is assumed based on arguments of local free convection (Wyngaard et al. 1971), the result in Eq. ( 7) depends solely on log (−L/z o ) (Zilitinkevich et al. 1992;Liu et al. 2023).This same dependency has been derived for a mixed layer velocity scale using matched asymptotic expansions with three scaling layers spanning the surface layer and no assumptions for z m (Tong and Ding 2020).
The results in Sect.3.2 and Fig. 6 yield z m ≈ 0.4z i for all LES cases based on the logarithmic decay of the gradients.Considering Eq. ( 6) aligns with the mean profiles up to the mixed layer in Figs.7a and 7b, the present findings indicate that the mixed layer values in Eq. ( 7) should depend on both z i /z o and z i /L in a complex manner.To test this, the mixed layer mean values at z m = 0.4z i are shown in Fig. 9.The values in Fig. 9 are insensitive to the choice of z m /z i within the range 0.4−0.8due to the uniformity of the mixed layer seen in Fig. 1.Because U m and θ m are mean values, the purpose of Eq. ( 7) and Fig. 9 is to assess the dependencies discussed above and not to define new velocity and temperature scales for the mixed layer.Included in Fig. 9 for comparison are the predicted values from numerical integration of Eq. ( 6) up to z m .It may be possible to evaluate further the integral definition for ψ(z m /L, z m /z i ) (see, e.g., Physick and Garratt 1995;De Ridder 2010), but for the present discussion the primary concern is the log (z m /z o ) term in Eq. ( 7) that is already uncoupled from ψ.
The velocity statistics in Fig. 9 are supplemented with results from two recent LES studies (Tong and Ding 2020;Liu et al. 2023).The Tong and Ding (2020) points only represent their LES cases with the Kosović (1997) SGS model.Further, the U m values reported in Tong and Ding (2020) are a velocity scale for the mixed layer that differs in value from the mean velocity (see, e.g., their Fig.3), such that the U m values used in Fig. 9 were inferred from their CBL profiles.The differing log-linear slope in Fig. 9a for the present LES and the cases from Liu et al. (2023) indicate that log (−L/z o ) is not the sole determinant of U m /u * .In particular, increasing z i /z o leads to a vertical shift in the resistance value.
The predicted values (closed symbols) for the external references use the same parameter values fitted to the present LES.The combination of the corrected similarity expression in Eq. ( 6) and the fixed height z m = 0.4z i lead to an accurate prediction of the differing trends noted above for Fig. 9a.The alignment with the reference studies is promising regarding the potential generality of the exponential cutoff correction.
The roughness dependence leading to a vertical shift in the mixed layer resistance value can be offset by plotting the diabatic term ψ m as in Fig. 9c.In this format, the data are approximately aligned along a common curve, noting that the residual differences may be due in part to the spin-up procedure that can significantly affect L (Fig. 2) and differences between SGS models as observed in Tong and Ding (2020).Importantly, ψ will also vary with roughness as z o becomes large, but z o ≪ z m for these data and its contribution to ψ is negligible.Fig. 9 supports z m ∝ z i and the dependence of U m /u * on both z i /z o and z i /L.However, the combination of Eqs. ( 6) and ( 7) do not form a drag law and further work is required to develop the result into a velocity scale.Further, no conclusion can be made for temperature due to lack of independent data across a range of z i /z o .The mixed layer temperature resistance is included for the present LES in Figs.9b and 9d for completeness.

Implications for the surface layer height
The extended logarithmic decay of φ as the gradients vanish to zero in Fig. 6 provide a consistent criterion for defining the bottom of the mixed layer where free convection occurs, but defining the top of the surface layer z SL is more ambiguous.The z/z i cutoff correction in Eq. 6 is non-negligible within the full range of heights analyzed here, from approximately 0.05z i up to z m , where the correction is required to account for simulation trends observed in Fig. 8 and in previous studies (e.g., Khanna and Brasseur 1997;Pirozzoli et al. 2017;Li et al. 2018).Upon further considering the nonlocal contribution of large-scale eddies to the decay in the gradients at these heights (Li et al. 2018;Fodor et al. 2019), it is possible that the extended range from z SL (not yet defined) to z m is a transition region resulting from the coexistence of local and nonlocal eddies governed by different mechanisms and scales.Here, the exponential cutoff approximates the transition in the mean profiles from Monin-Obukhov similarity in the surface layer to the zero gradient condition characterizing the mixed layer.This transition is depicted in Fig. 10.
Assuming the exponential cutoff and fitted constant c in Fig. 10a represents the appropriate correction to surface similarity, the correction is less than 10% only for heights in the lowest few percent of the CBL depth.For z i ∼ O(1 km), these heights correspond to the lowest 30 m or so of the atmosphere where many field measurements are sampled.On one hand, this indicates that the correction is not significant for many previous field campaigns and explains why the original Businger-Dyer profiles align well with experimental data to within the uncertainty and scatter of the results.On the other hand, a 10% correction at 0.03z i is nonnegligible and emphasizes the need to reconsider the conventional height for z SL .
Derivations for the generalized log law U/u * = κ −1 [log (z/z o ) − ψ m (ζ)] in the inertial sublayer of the ABL often rely on asymptotic matching between the surface conditions and the outer layer where the velocity defect follows Rossby-number similarity (see, e.g., Blackadar and Tennekes 1968;Hess 1973).These derivations use the assumption z/z i → 0 (i.e.z ≪ z i ) to define z SL ∼ 0.1z i and do not consider the extensive uniformly mixed layer in convective conditions.Owing to the influence of the mixed layer on near-surface statistics, a new assumption z/z m → 0 for the surface layer depth is suggested here.While it may be possible to formalize a revised asymptotic matching between the surface conditions and the mixed layer, the matching is complicated by the dependence of the mixed layer resistance on surface properties as seen in Eq. ( 7).
In Fig. 10, the approximation z SL ≈ 0.1z m ≈ 0.04z i is applied under the assumption that z SL ≪ z m is required for the contribution of nonlocal eddies to the gradient -and the resulting correction factor -to be small.However, the correction is greater than 10% at this height such that a more stringent definition may be warranted.A formal framework is required to address this issue of surface layer height in a more rigorous manner.Regardless of the exact definition for z SL , the traditional estimate z SL = 0.1z i is insufficient for convective conditions based on the growing body of evidence discussed in the introduction.
The question of defining z SL is complemented by recent evidence that mean profile similarity in the surface layer of the stable ABL also depends on the boundary layer depth (Heisel and Chamecki 2023).While the flow structure for stable conditions is considerably different, the same general reasoning applies: z-less stratification above the surface layer indirectly influences the gradients below 0.1z i .Here, free convection turbulence above the surface layer influences the gradients below 0.1z i in a more direct manner.
The transition region depicted in Fig. 10b coincides with the regime previously discussed as a local free convection layer (Tennekes 1970;Kader and Yaglom 1990) or convection matching layer (Panofsky 1978;Kaimal and Finnigan 1994).However, formulations for those layers do not transition smoothly into the uniformly mixed layer (Panofsky 1978), unlike the present work and similar studies of the so-called radix layer (Santoso andStull 1998, 2001).There is also a distinction in the observed scaling, at least for the first-order statistics.In the traditional local convection layer with −L ≪ z ≪ z i , the z scaling is still relevant but the velocity and temperature variables no longer depend on u * (Wyngaard et al. 1971).In the present study, the exponential cutoff reduces the gradient based on the relative position within the transition region, but does not alter the scaling for φ or ζ.As an analogy, the fluxes −u ′ w ′ /u 2 * decay with z/z i but scale with u * throughout the CBL depth.In this sense, the gradients maintain their surface scaling up to z m despite the incomplete similarity due to the z/z i correction demonstrated in Fig.

8.
This z/z i correction corresponds to the presence of large-scale eddies such as downdrafts (Li et al. 2018) and updrafts (Fodor et al. 2019) whose governing scale w * is oriented along z (Deardorff 1970).Within the mixed layer, the ensemble of these predominately vertical motions results in zero mean gradient.If the free convection events that extend into the transition region also have a collective mean gradient close to 0, the events would contribute to a decrease in the overall φ without incurring a statistically meaningful transition in scaling from u * to w * .The decrease in φ would then directly depend on the relative probability of free convection events that increases with height in a manner consistent with the exponential in Fig. 10a.Importantly, the same argument cannot be extended to the higher-order variance statistics and an analysis of the variances is outside the scope of the present work.
There is uncertainty in Fig. 10 with respect to changes in the behavior with increasing instability.The present analysis suggests the exponential cutoff and parameter c do not vary with z i /L.As noted previously, this suggests the correction and the relative gradient contributions are independent of changes in the eddy topology from roll structures to vertical cells.Further, with z m ∝ z i the surfacescaling region in Fig. 10b must increasingly resemble the mixed layer turbulent structure in order for the surface layer to vanish in free convection.These impli-cations should be evaluated across a wider range of z i /L and for a greater number of cases before conclusions are drawn.

Summary
The present work uses a series of seven LES cases to study mean profile similarity in the lower half of the convective boundary layer.The cases represent a dry, barotropic idealized CBL under weak to moderately strong convection with midlatitude Coriolis frequency, a stable capping inversion, and a well-defined mixed layer.An ad hoc spin-up procedure is used to mitigate inertial oscillations in the final simulations (Fig. 2), and the grid converge of near-surface profile statistics is closely examined (Figs. 3 and 4).
The new LES cases reveal the same qualitative trends in the nondimensional gradients φ seen in other recent simulation-based studies (Fig. 5): the results generally align with Monin-Obukhov similarity across the different cases, but the individual profiles each exhibit a steeper slope that precludes universal similarity in φ(ζ) (Pirozzoli et al. 2017;Maronga and Reuder 2017;Li et al. 2018).In other words, MOST captures variability in φ across a range of L and fixed z, but does not fully account for the variability across z for fixed L. The behavior of the φ profiles above the surface layer indicates that the steeper slope is associated with a broader trend in z/z i that reduces the gradient towards 0 at the height of the mixed layer (Fig. 6).
To account for this trend, the well-known Businger-Dyer profiles are revised in Eq. ( 6) with an exponential cutoff similar to corrections for similarity in the roughness sublayer (Garratt 1980;Harman and Finnigan 2007).The revised expressions, with fitted parameters c m ≈ 3.7 and c h ≈ 2.9 for the exponential term exp (−cz/z i ), result in significantly improved similarity for the mean (Fig. 7) and gradient (Fig. 8) profiles from approximately 0.05z i up to the mixed layer near 0.4z i .The correction is expected to be small close to the surface where most point measurements are acquired in field experiments, which may explain why the consistent z/z i trend seen in simulations is not readily apparent from field observations.Further, the parameter space probed by field measurements often spans a wide range of L at a limited series of fixed heights, which as noted above can yield curves that closely follow MOST relations.
In addition to the improved similarity, there are three important implications arising from the revised relations in Eq. ( 6).First, the mean values in the mixed layer depend on both log (z i /z o ) and ζ (Fig. 9).Second, the exponential correction accounts for an extended transition region in the mean profiles between the surface layer and the mixed layer (Fig. 10), where this region is strongly influenced by large-scale buoyancy-driven eddies (Li et al. 2018;Fodor et al. 2019).Owing to the effect of these eddies, the correction only becomes small for z/z i ∼ O(0.01), such that the common assumption of 0.1z i for the surface layer height is too large for idealized convective conditions.Third, the start of the mixed layer at a fixed fraction of z i (Fig. 6) suggests the surface layer turbulent structure changes with increasing convection in order to match the mixed layer under free convection, but extending the analysis to stronger convection is required to validate the last point.
While the fitted results in Figs.7 and 8 are promising, the present analysis includes a limited number of cases and lacks reliable statistics in the bottom half of the surface layer.In particular, modest grid dependence is observed across the surface layer for weak convection in case A (Fig. 4).Equation ( 6) is thus considered to be a preliminary effort to model the trends in φ observed in Figs. 5 and 6.Several other functional forms were evaluated to more accurately account for the logarithmic decay in φ, but the collective profile data were found to be prone to overfitting, where the functional dependencies were borne by the fitted parameters.Accordingly, the proposed correction is limited to an extension of the widely-tested Businger-Dyer profiles, and the correction has only one nondimensional parameter (c m or c h ) that has a clear physical interpretation corresponding to the height where the cutoff reaches a given magnitude.Importantly, the present quantitative correction depends on the empirical profiles being corrected, i.e.Businger-Dyer, O'KEYPS, or another alternative.However, the lack of universality in φ(ζ) in Fig. 5 and the necessity to include z/z i in mean profile similarity are general findings that do not depend on the specific empirical relations.Further analysis with additional datasets may support a more sophisticated model that better matches the full gradient cutoff up to 0.4z i in Figs.8a and 8b, but for the available data in the present study only a simpler model is warranted.

Fig. 2
Fig.2Demonstration of the employed spin-up procedure and the impact of initial conditions on inertial oscillations.For the initial velocity profiles in (a), the resulting average surface shear stress τs (b), mixed layer wind speed Um (c), and stability z i /L (d) are plotted over time for one inertial period.During a given spin-up, the velocity profiles at tf = π are rescaled to the original z i and used as the initial profiles for the subsequent spin-up, as shown in the inset of (a).Parameters with overbars ( • ) indicate the long-term average value over the full period.

Fig. 3
Fig. 3 Nondimensional mean wind speed (a) and air temperature (b) with varying grid resolution for cases A (top), C (middle) and E (bottom).Means are expressed as the diabatic term ψ defined in Eq. (4).Short vertical lines indicate the start of the converged region where the change in ψ is less than 1% between the finest two grids.The means are normalized using the von Kármán constant κ = 0.39.

Fig. 4
Fig. 4 Nondimensional gradients of velocity (a) and air temperature (b) with varying grid resolution for cases A (top), C (middle) and E (bottom).Gradients are expressed as φ defined in Eq. (1).Short vertical lines indicate the start of the converged region where the change in φ is less than 0.01 between the finest two grids.

Fig. 5
Fig. 5 Nondimensional profiles as a function of ζ = z/L, compared to Businger-Dyer relations (thin lines) from Eqs. (3) and (5).Rows correspond to the mean ψ (a,b) computed from Eq. (4) and gradient φ (c,d) computed from Eq. (1), and columns correspond to momentum (a,c) and heat (b,d).Shaded regions in (c,d) represent 95% confidence intervals for statistical uncertainty in the mean gradients.Each profile includes heights from the bottom of the converged region determined in Sect.2.2 up to 0.3z i , and short lines indicate the height z = 0.1z i for reference.

Fig. 6
Fig. 6 Nondimensional gradient profiles as a function of z/z i .Rows correspond to logarithmic (a,b) and log-linear (c,d) axis scaling, and columns correspond to momentum φm (a,c) and heat φ h (b,d).Transparent regions of each curve are considered unconverged for the present resolution based on the assessment in Fig. 4.

Fig. 7
Fig.7Comparison of LES mean profiles (thick lines) with the best fit of the Eq.(6) integral (dashed lines) and the same equation without the exponential cutoff (dotted lines).Rows correspond to the mean value (a,b) and the diabatic term ψ (c,d), and columns correspond to momentum (a,c) and heat (b,d).Transparent regions of the curves in (a,b) are considered unconverged for the present resolution based on the assessment in Fig.3.Each profile in (c,d) includes heights from the bottom of the converged region up to 0.3z i , and short lines indicate the height z = 0.1z i for reference.

Fig. 8
Fig. 8 Comparison of LES φ profiles (thick lines) with the best fit of Eq. (6) (dashed lines) and the same equation without the exponential cutoff (dotted lines).Rows correspond to the profiles plotted versus z/z i (a,b), plotted versus ζ (c,d), and after compensating for the exponential cutoff (e,f).Columns correspond to momentum (a,c,e) and heat (b,d,f).Transparent regions of the curves in (a,b) are considered unconverged for the present resolution based on the assessment in Fig. 4.Each profile in (c-f) includes heights from the bottom of the converged region up to 0.3z i , and short lines indicate the height z = 0.1z i for reference.

Fig. 9
Fig. 9 Dependence of the mean velocity Um and temperature θm in the convective mixed layer.Rows correspond to the mean values versus L/zo (a,b) and their diabatic terms versus z i /L (c,d), and columns correspond to velocity (a,c) and temperature (b,d).Results are included for the present LES and two reference LES studies indicated by open symbols, where each are compared to predicted values based on Eq. (6) indicated by closed transparent symbols.Observed and predicted values are based on the mixed layer bottom height zm = 0.4z i .Symbol color corresponds to z i /zo.
Depiction of the transition region between the top of the surface layer z SL and the bottom of the mixed layer zm.(a) Value of the exponential cutoff in the gradient momentum and heat.Contribution of both surface scaling and mixed layer free convection events to the mean gradients within the transition region, where the relative contribution of each varies with height in conjunction with the cutoff function in (a).

Table 1
Key scaling parameters for large-eddy simulations (LES) of the convective atmospheric boundary layer (CBL) on a 1024×1024×512 numerical grid: imposed geostrophic wind speed Ug, imposed surface heat flux Q * , friction velocity u * , surface temperature scaling θ * , Obukhov length L, boundary layer depth based on the capping inversion height z i , and bulk instability parameter z