Estimating the Diffusive Heat Flux across a Stable Interface Forced by Convective Motions

Entrainment at the top of the convectively-driven boundary layer (CBL) is revisited using data from a high-resolution large-eddy simulation (LES). In the range of values of the bulk Richardson number Ri B studied here (about 15–25), the entrainment process is mainly driven by the scouring of the interfacial layer (IL) by convective cells. We estimate the length and time scales associated with these convective cells by computing one-dimensional wavenum-ber and frequency kinetic energy spectra. Using a Taylor assumption, based upon transport by the convective cells, we show that the frequency and wavenumber spectra follow the Kolmogorov law in the inertial range, with the multiplica-tive constant being in good agreement with previous measurements in the atmosphere. We next focus on the heat flux at the top of the CBL, F i , which is parameterized in classical closure models for the entrainment rate w e at the interface. We show that F i can be computed exactly using the method proposed by Winters et al. (1995), from which the values of a turbulent diffusivity K across the IL can be inferred. These values are recovered by tracking particles within the IL using a Lagrangian stochastic model coupled with the LES. The relative difference between the Eulerian and Lagrangian values of K is found to be lower than 10%. A simple expression of w e as a function of K is also proposed. Our results are finally used to assess the validity of the classical " first-order " model for w e. We find that, when Ri B is varied, the values for w e derived from the " first-order " model with the exact computation of F i agree to better than 10% with those computed directly from the LES (using its definition). The simple expression we propose appears to provide a reliable estimate of w e for the largest values of Ri B only.


Introduction
An interfacial layer (IL) divides the clear convective atmospheric boundary layer (CBL) and the stably-stratified free atmosphere (FA) above.The IL is forced by turbulent motions, which are primarily triggered by ground surface heating.Indeed, the main mechanism of turbulence production within the well-mixed part of the CBL (referred to as the mixed layer) is buoyant convection, with a possible windshear contribution.The penetration of rising thermals into the FA is associated with an entrainment of air down into the mixed layer (e.g.Sorbjan, 1996).As a result, the CBL deepens or equivalently the IL raises.A mixed layer with similar structure and dynamics also forms in the upper ocean when cooling occurs at the surface.As pointed out for instance by Stevens and Lenschow (2001), the modeling of the entrainment process is an essential issue in any attempt to parameterize the CBL in large-scale models, for both atmospheric and oceanic applications.Indeed, only a few parameterization schemes of boundary-layer flow within meso-scale models represent explicitly the entrainment process (e.g.Hong et al., 2006).The representation of the entrainment process is also an issue for air quality prediction.Mean vertical gradients of concentrations of atmospheric constituents are close to zero within the mixed layer.Hence, the rising rate of the mixed layer into the FA determines partly the concentrations at the ground surface (e.g.Cai and Luhar, 2002).
The entrainment process across a buoyancy interface (such as the IL) due to turbulent motions has been studied extensively in laboratory experiments.Hopfinger (1987) and Fernando (1991) gave a thorough review for an IL that is forced by grid turbulence.In grid turbulence experiments, the entrainment process is discussed classically as a function of a bulk Richardson number at the interface, R i = bl 0 /(σ u ) 2 C. Chemel et al.: Estimating the diffusive heat flux across a stable interface forced by convective motions is the buoyancy jump across the interface, and l 0 and (σ u ) 2 0 are the integral length scale and variance of the turbulence in the absence of the interface.Entrainment and resulting mixing of entrained fluid occur roughly for 1<R i <50.Hannoun and List (1988) showed that for this range of R i values, mixing results from internal gravity wave breaking at the interface.Below this range, mixing occurs as if the fluid were homogeneous, and above that range, mixing occurs through pure molecular diffusion.As stressed for instance by Sullivan et al. (1998), the extension of the results from grid turbulence experiments is actually debatable since the CBL contains large-scale organized structures, which are not present in such experiments.
The thermally-driven convection tank experiment of Deardorff et al. (1980) was designed to mimic the CBL dynamics.The stability of the buoyancy interface was also characterized by a bulk Richardson number at the interface, defined as where g is the gravitational acceleration, β the coefficient of thermal expansion, the potential temperature1 jump across the interface, z i the mixing depth, and w * the convective velocity (defined below).The quantities β, , z i and w * refer to horizontally averaged quantities.The convective velocity within the mixed layer is expressed as where F s = w s is the horizontally averaged heat flux just above the surface.Typical vertical profiles of potential temperature and heat flux F = w are depicted in Fig. 1a (and are further discussed in Sect.4.4).For a broad range of Ri B values (about 2-85), which are usually observed in the atmosphere, the entrainment rate, w e = d t z i (with d t ≡ d/dt), was found to vary as where the dimensionless parameter A is the entrainment ratio and is close to 0.25.There is actually a wide spread in the values of A reported in the literature, though it is often found to be around 0.2 in the regime of equilibrium entrainment (i.e. when the mixed layer dynamics has reached a quasisteady state).The parameterization of A, and more generally of w e , is at the heart of the debate on entrainment.
The most common parameterizations of w e are the socalled "zero-order" and "first-order" jump models, which were proposed by Lilly (1968) and Betts (1974), respectively.In the "zero-order" model, the thickness of the IL is assumed infinitesimal, while the potential temperature profile exhibits a jump across that interface.In "first-order" models, the finite thickness of the IL is taken into account (see for instance Fig. 1b).In both models, the altitude of the IL is the mixing depth z i and is defined as the level where the heat flux F i is minimum (being negative).The main issue in these models is to derive a closure for this flux.Fedorovich et al. (2004) presents more general formulations for the entrainment law and reviewed methods for determining entrainment parameters from large-eddy simulation (LES) outputs.Sullivan et al. (1998) used LESs to investigate the convective entrainment process and the structure of the IL over a wide range of Ri B values (about 15-45).The authors showed that the finite thickness of the IL needs to be considered in an entrainment law formulation derived from a jump model.In other terms, the "zero-order" jump model was found insufficient, especially at low Ri B .Conversely, the "first-order" jump model was found to work well.Fedorovich et al. (2004) also found that the "zero-order" parameterization is insufficient outside the regime of equilibrium entrainment.
In this study, we present results from a high-resolution LES of the convectively-driven boundary layer initialized by a commonly used sounding of Day 33 of the Wangara experiment (Clarke et al., 1971).Our main purpose is to show that the heat flux at the interface (i) can be computed exactly, using the method proposed by Winters et al. (1995), and (ii) can be expressed in terms of a vertical turbulent diffusivity K, which we also obtain from Lagrangian particle tracking.This allows us to assess the validity of the commonly used "first-order" model to parameterize w e and to provide a simple expression of w e in terms of K.
Several LES studies have been conducted to investigate the entrainment process in the CBL (see for instance Stevens and Lenschow, 2001).Sorbjan (1996) carried out LES experiments to analyse the effects, on entrainment, of the vertical potential temperature gradient in the FA, hereafter denoted by FA (see Fig. 1).The entrainment rate was found to depend on FA but the entrainment ratio varied only slightly in the range 0.2-0.3 for values of FA from 1 to 10 K km −1 (which are usually observed in the atmosphere).In addition, the statistical moments in the lower 90% portion of the mixed layer were found almost independent of FA .Lewellen and Lewellen (1998) also examined the convective entrainment process and stressed that the entrainment rate is controlled by the turbulent transport at the scale of the boundary layer and is relatively insensitive to the smaller scales of mixing near the IL.This confirms earlier findings by Linden (1975) (for grid-generated turbulence), Manins and Turner (1978), and Schmidt and Schumann (1989).Otte and Wyngaard (2001) focused on the properties of the IL and found that turbulence there behaves as in stably-stratified flows, consistent with the work of Hannoun and List (1988).
The outline of the paper is as follows.A description of the LES and Lagrangian stochastic models is given in Sect. 2. In Sect.3, we focus on characteristics of the mixed layer turbulence which forces the IL.We compute one-dimensional wavenumber as well as frequency spectra and discuss the length and time scales involved in the entrainment process.
2 Model description and setup

The LES model
The numerical experiments presented in this paper were conducted with the Advanced Regional Prediction System (ARPS), a non-hydrostatic, compressible LES code devoted to meso-scale and small-scale atmospheric flows.Xue et al. (2000;2001) gave an extensive description of the model formulation and applications.
The basic idea of physical LES is the 'filtering approach' to separate the small scales from the large scales (see Lesieur and Métais, 1996, for a review).In this approach, a lowpass spatial filter (denoted by a tilde hereafter) is applied to the turbulent fields.In the present study the characteristic width of the filter ∆ is equal to the geometric average of the grid size in the three spatial directions.The application of this filter to the mass-and momentum-conservation equations, assuming that the filtering operation commutes with differentiation, results in where u, p, and ρ are the velocity, the pressure and the density fields, respectively, µ is the dynamic viscosity, and the subscripts (i, j) ∈ {1, 2, 3} refer to the geometrical coordinates.For convenience, we will also adopt the following notation: (x 1 , x 2 , x 3 ) ≡ (x, y, z) and (u 1 , u 2 , u 3 ) ≡ (u, v, w).The terms τ ij = u i u j − u i u j and −2Ω × u, with Ω being the Earth's angular velocity, represent the subgridscale (SGS) turbulent stress and the Coriolis acceleration, respectively.The SGS term must be parameterized as a function of the filtered variables.For this purpose, an eddyviscosity model is used, namely where δ ij is the Kronecker delta symbol and S ij = (∂ j u i + ∂ i u j ) /2 is the filtered strain-rate tensor.ν t is the SGS turbulent viscosity, which is expressed as a function of the filtered variables through a mixing length formulation.This formulation yields ν t = 0.1 ℓ e 1/2 , where e = τ kk /2 is the turbulent kinetic energy of the subgrid scales and ℓ is a typical subgrid length scale, which accounts for the effects of stratification (Deardorff, 1980).For a grid size with an aspect ratio in the order of unity, ℓ is equal to ∆ for unstable or neutral cases and min( ∆, 0.76 √ e N −1 ) for stable case, where is the buoyancy frequency.Note that for a larger aspect ratio, we need to set the vertical length scale apart from the horizontal one.The prognostic equation for e is In Sect.4, we show that the heat flux at the interface can be computed exactly and that a turbulent diffusivity across the interface can be inferred, whose values are recovered from the tracking of fluid particles within the interface.These results are finally used to assess the validity of the "first-order" model.Conclusions are given in Sect. 5.
2 Model description and setup

The LES model
The numerical experiments presented in this paper were conducted with the Advanced Regional Prediction System (ARPS), a non-hydrostatic, compressible LES code devoted to meso-scale and small-scale atmospheric flows.Xue et al. (2000;2001) gave an extensive description of the model formulation and applications.
The basic idea of physical LES is the "filtering approach" to separate the small scales from the large scales (see Lesieur and Métais, 1996, for a review).In this approach, a lowpass spatial filter (denoted by a tilde hereafter) is applied to the turbulent fields.In the present study the characteristic width of the filter is equal to the geometric average of the grid size in the three spatial directions.The application of this filter to the mass-and momentum-conservation equations, assuming that the filtering operation commutes with differentiation, results in where u, p, and ρ are the velocity, the pressure and the density fields, respectively, µ is the dynamic viscosity, and the subscripts (i,j ) ∈ {1,2,3} refer to the geometrical coordinates.For convenience, we will also adopt the following notation: (x 1 , x 2 , x 3 ) ≡ (x, y, z) and (u 1 , u 2 , u 3 ) ≡ (u, v, w).
The terms τ ij = u i u j − u i u j and −2 × u, with being the Earth's angular velocity, represent the subgrid-scale (SGS) turbulent stress and the Coriolis acceleration, respectively.
The SGS term must be parameterized as a function of the filtered variables.For this purpose, an eddy-viscosity model is used, namely where δ ij is the Kronecker delta symbol and S ij = ∂ j u i + ∂ i u j /2 is the filtered strain-rate tensor.ν t is the SGS turbulent viscosity, which is expressed as a function of the filtered variables through a mixing length formulation.This formulation yields ν t = 0.1 e 1/2 , where e = τ kk /2 is the turbulent kinetic energy of the subgrid scales and is a typical subgrid length scale, which accounts for the effects of stratification (Deardorff, 1980).For a grid size with an aspect ratio in the order of unity, is equal to for unstable or neutral cases and min ,0.76 √ eN −1 for stable case, where The prognostic equation for e is where ε = C ε e 3/2 / is the dissipation rate of turbulent kinetic energy and the coefficient C ε has the value 3.9 at the lowest vertical level and 0.93 otherwise.The turbulent Prandtl number is parameterized as Pr t =1/ 1 + 2 / , from where the SGS turbulent thermal diffusivity κ t =ν t /Pr t can be inferred.
The energy-conservation equation for is written as where λ is the thermal conductivity, c p is the specific heat at constant pressure, and ϕ j = u j − u j is the SGS turbulent heat flux, which is expressed as a function of the filtered potential temperature gradient ∂ j , namely

Model setup
The model is initialized using vertical profiles of potential temperature, horizontal wind, and vapor mixing ratio taken at 09:00 EST during Day 33 of the Wangara experiment held in Hay, Australia (Clarke et al., 1971).The wind profile is almost shear-free up to the top of the domain and the vertical potential temperature gradient in the FA, FA , is about 10 K km −1 .The ground surface is heated through the absorption of solar radiation.This results in a diurnal variation in the ground surface temperature and turbulent heat fluxes, which trigger convective motions.
A good representation of land surface characteristics was found necessary to reproduce realistically the atmospheric boundary-layer structure and its evolution.The land-surface energy budget was calculated by a simplified soil-vegetation model (Noilhan and Planton, 1989;Pleim and Xiu, 1995).The soil type was loam and the vegetation type was desert.The roughness length was 0.24 m, the leaf area index was 0.1 and the fractional vegetation coverage was 5%.The ground surface temperature was initialized to its observed value at 09:00 EST (278.7 K).There was no direct measurement of the deep surface temperature, so that its initial value was evaluated from the ground surface temperature and soil heat flux.This flux was about zero at 08:00 EST (Clarke et al., 1971) indicating that the deep soil temperature was approximately the same as the ground surface temperature at that time (274.0K).Assuming that the deep soil temperature did not vary from 08:00 EST to 09:00 EST, it was initialized to 274.0 K.Both ground surface and deep soil moisture were set to the wilting point as suggested by Clarke et al. (1971) since it had not rained for many days.
In the numerical code, periodic lateral boundary conditions are prescribed and a rigid wall condition is applied at the bottom and top of the domain (with a Rayleigh sponge close to the top boundary).The computations are performed on a 5.12 km×5.12 km×4.535km domain with 256 grid points in each direction.The vertical resolution is 20 m over the bulk of the boundary layer, 5 m within the IL and 50 m far above.A gradually-varying mesh size is employed near the transition zones.Such a rather fine grid has been selected to let turbulence develop with minimal bias due to the aspect ratio of the grid size and to have a fair representation of the IL.Indeed, earlier LES investigations show that only high-resolution LESs would provide reliable estimates of the entrainment rate (see for instance Bretherton et al., 1999;Stevens and Lenschow, 2001).

The Lagrangian stochastic model
A Lagrangian particle dispersion model has been implemented in the ARPS code to track a large number of particles, following Weil et al. (2004) and Vinkovic et al. (2006).Let x p0 be the particle position at initial time and x p x p0 ,t its position at time t.The trajectory of the fluid particles is computed by integrating the equation where v is the Lagrangian velocity of the particles.This velocity is decomposed into (Lamb, 1978) It involves an Eulerian filtered part u x p ,t and a fluctuating SGS contribution v x p ,t , which is modeled by a modified three-dimensional Langevin model.The ith component of the Lagrangian velocity v is given by the stochastic differential equation (Thomson, 1987) where γ i x p ,v,t + α ij x p ,t v j − u j is a deterministic forcing function composed of a filtered contribution γ i x p ,v,t and a fluctuating SGS contribution α ij x p ,t v j − u j .The last term in Eq. ( 11) is a random forcing with dη j (t) being an isotropic Gaussian white noise with variance dt (namely dη i t dη j t = δ ij δ t − t dt where stands for time correlation).In the present study, the SGS turbulence is assumed homogeneous and isotropic, so that these terms are given by (see Weil et al., 2004;Vinkovic et al., 2006, for details) Nonlin.Processes Geophys., 17, 187-200, 2010 www.nonlin-processes-geophys.net/17/187/2010/ where C 0 is the Lagrangian constant (Thomson, 1987).It follows that the Lagrangian velocity is obtained by integrating the equation The filtered velocity u at the position of the particle was obtained from the gridded computed Eulerian velocity by a cubic spline interpolation procedure.The time integration was performed using a fourth-order Runge-Kutta scheme.At the boundaries, particles that moved out of the domain were forgotten.

Boundary-layer structures
The development of the clear and shear-free CBL in a stablystratified atmosphere has been studied in several papers (see Fedorovich et al., 2004, for a review).In this situation the warm underlying surface is the unique source that triggers convection.Therefore convective cells do not oscillate as Rayleigh-Benard cells do (Matthews and Cox, 2000).These cells initiate from the heated ground surface, grow and decay after a finite lifetime.Though the spatial distribution of the cells is determined by interactions with the surrounding cells, their location is unpredictable.However, all the key length scales (e.g.horizontal extension and distance between cells) must be related to the height of the CBL, since it is the only length scale in this problem.
An overview of the boundary-layer structure is shown in Fig. 2, where contours of the field in the range 285.7< <287.7 K are displayed in a (x,z) plane located near an updraft at 15:00 EST.The mixed layer as well as the IL are strongly turbulent, implying that the thickness of the IL has a high variability.Fast-rising localized updrafts impact the IL (which leads to a folding of the interface) or erode the interface by a scouring mechanism.Similar downward plumes transport heat from the cooled sea surface toward the bottom of the oceanic mixed layer (D'Asaro et al., 2002).These entrainment events are localized and turbulent motions mix the entrained air downwards.The typical horizontal size of the convective cells at z/z i =0.25 is found to be in the range 1500-2000 m, which is in good agreement with that found in previous studies (e.g.Schmidt and Schumann, 1989), and is in the order of the mixing depth z i .
To identify the instantaneous structure of these cells, we use the Q-criterion (Hunt et al., 1988).This criterion is derived from the second invariant of the fluctuating velocity gradient tensor ∇ u, denoted Q, which is expressed as where R ij = (1/2) ∂ j u i − ∂ i u j and S ij are the antisymmetric and symmetric parts of ∇ u, respectively.The Q-criterion may be regarded as the competition between the rotation rate R 2 = R ij R ij and the strain rate S 2 = S ij S ij .Thus, positive Q isosurfaces highlight areas where the rotation rate overcomes the strain rate, which are therefore eligible as vortex envelopes.
An isosurface of Q of positive value is displayed for an isolated convective cell pattern at 15:00 EST in Fig. 3a.The highlighted structures are traces of the fast rising updrafts, which are characterized by strong vorticity components.A horizontal cross-section of these updrafts is visible in Fig. 3b, where the vertical velocity scaled by w * is plotted at the same time.It is then possible to discuss the degree of organization of the convective cells.They consist of well organized updrafts, which vanish and diffuse at the interface creating broad ring-shaped patterns.The air mass is gradually mixed downwards in the center of the pattern.Downdrafts are not organized compared to updrafts, and lead to downward mixing associated with small-scale turbulent structures.The fast rising updrafts occupy a smaller fraction (about 40%) of the CBL horizontal cross-sectional area than the slowly broader downdrafts (see Fig. 4), due to the vanishing of the vertical mass flux averaged over a horizontal surface.The values of this fraction are consistent with observational data (Lenschow, 1998).

Mixed-layer statistics
The statistical properties of the mixed layer have been studied in several papers (e.g.Moeng and Wyngaard, 1988; www.nonlin-processes-geophys.net/17/187/2010/ Nonlin.Processes Geophys., 17, 187-200, 2010  3a has a horseshoe shape, which is clearly visible in Fig. 3b for approximately 0.2≤y≤0.6 and 0.2≤(−x/L + 1)≤0.6.Peltier et al., 1996;Kelly and Wyngaard, 2006), leading to the conclusions that the kinetic energy density and temperature variance spectra obey Kolmogoroff and Corrsin-Oboukhov laws, respectively, the universal constants in these spectra being also recovered.Our interest in this section is to show that the frequency spectrum for the horizontal velocity matches precisely the spatial one-dimensional longitudinal spectrum when a Taylor assumption, based upon transport of fluctuations by the convective cells, is made.The computation of the spatial and frequency spectra also provides the typical length and time scales, respectively, of the mixed layer.The one-dimensional longitudinal spectra of kinetic energy density for u and v, denoted by E xx (k x ,t) and E yy k y ,t , respectively, are displayed in Fig. 5a.For homogeneous and isotropic turbulence, the spectra should behave as where k n is the wavenumber, the subscript n denoting x or y, and C 1 =(18/55)C K ≈ 0.49 for C K = 1.5 (e.g.Champagne et al., 1977;Moeng and Wyngaard, 1988).The computed constant C 1 averaged for the u and v compensated spectra is displayed versus k n in Fig. 5b.The value of C 1 thus obtained agrees well with the theoretical prediction of 0.49 for the smallest scales of the inertial range.We also note that these spectra are nearly the same, as expected from local isotropy.Two wavenumbers are indicated in Fig. 5a and b, denoted by k i and k v .The former, k i , is the wavenumber at which the two-dimensional spectra of ũ and ṽ peak (not shown).Therefore i = 2π/k i is the integral scale of turbulent motions.The latter, k v , is defined as 2π/ v , where v is the effective dissipative scale, namely v = ν t 3 / 1/4 , with and ν t being inferred from the SGS model.The computed integral scale i is close to 1900 m, which corresponds to the typical size of the convective cells.In the present LES, the dissipative scale v is in the order of 5 m.Since the large-scale flow within the mixed layer consists of convective cells, the scales contributing to the inertial range may be assumed to be advected by those cells.In other terms, the Taylor's (1938) frozen turbulence hypothesis may be assumed to hold.This reasoning also requires the magnitude of the velocity fluctuations to be much smaller than the convective velocity (see Peltier et al., 1996, p. 55, for a discussion of this point).Under these assumptions, a f −5/3 power law, with f being the frequency of the motions, is expected in the inertial range for the velocity components.The frequency spectrum of u, denoted by S u (f ), computed from 12:00 EST to 15:00 EST at z=500 m in the center of the (x,y) plane, is displayed versus f in Fig. 5c.A f −5/3 power law is obtained over almost a decade in the inertial range.Beyond f =0.03 s −1 , turbulence is significantly damped by the SGS turbulent viscosity.This frequency corresponds to a characteristic time scale of 30 s, which is approximately the time for disturbances to travel across a grid element in the mesh.
The eddy-turnover time τ i associated with the integral scale i may be estimated by τ i = i / u rms , where u rms is the root-mean-square of u.This yields a time of 15 min, which is the typical time for air to circulate between the ground surface and the top of the mixed layer, namely roughly z i /w * , as we checked it.The corresponding frequency f i = 2π/τ i is indicated in Fig. 5c.
Using the Taylor's hypothesis, frequency spectra can be converted to one-dimensional wavenumber spectra by substituting the frequency f for k x | u|.The one-dimensional wavenumber spectrum thus obtained is superimposed upon E nn in Fig. 5a.Both spectra remarkably coincide over the inertial range.This demonstrates the reliability of the Taylor's hypothesis within the mixed layer.
Thus, the turbulence within the mixed layer may be assumed to be locally homogeneous and isotropic over a broad range of scales in the inertial range.We checked that the IL is forced by the largest scales of the mixed layer by investigating the vertical evolution of the two-dimensional heat flux spectrum, as done by Schmidt and Schumann (1989) from LES results and by Kaiser and Fedorovich (1998) from wind tunnel measurements.In agreement with these authors, we found that the heat flux spectrum becomes negative at the largest scales as the IL is approached from below (not shown).This implies that heat is transferred down from the IL and that the largest scales of the mixed layer are involved in this process.

Entrainment rate formulation
In this section, the focus is directed onto the IL where entrainment events take place.As recalled in the introduction, the parameterization of the entrainment rate at the top of the CBL, w e , involves the (unknown) heat flux F i at the interface and hence a closure for this flux.In this section, we show that F i can be computed exactly from the method of Winters et al. (1995).Then we introduce a turbulent thermal diffusivity from F i , which we also compute by tracking Lagrangian fluid particles within the interface.This analysis is finally applied to the "first-order" model discussed in the introduction.
We first compute the characteristics of the IL from our LES, which are needed in the analysis of entrainment.    1 for the 256 3 resolution run, at successive times during the mixed layer growth.The mixing depth z i is defined as the level where the heat flux is minimum as in the standard flux method (e.g.Fedorovich and Mironov, 1995;Sullivan et al., 1998).The values of z i obtained in this way were compared with those computed from the gradient method, described for instance by Sullivan et al. (1998), for which z i corresponds to the height above ground level where ∂ 3 is maximum.Relative differences were found to be lower than 10%.The lower and upper limits of the IL were more difficult to determine.First, we have used computed values of the second derivative ∂ 2 3 .Indeed, ∂ 2 3 is expected to reach a maximum and a minimum at the lower and upper limits of the IL, respectively.Since large ∂ 2 3 values often occur close to the ground surface, it was computed from z i upward and downward to search the first minimum and maximum values, respectively.Nonetheless, this method was found to be not so accurate because of non representative local extrema of ∂ 2 3 .Thus, the thickness of the IL, δ 1 , was computed as z 2 −z 1 , where z 1 coincides with the zero-crossing height of the heat flux profile and z 2 is the vertical position where the heat flux first goes to zero above z i , as illustrated in Fig. 1a.Note that, consistent with the convection tank measurements of Deardorff et al. (1980), δ 1 /z i is close to 0.2 for strong enough stratification of the interface (see Table 1).The potential temperature jump was calculated as (z 2 ) − (z 1 ).The entrainment velocity was computed from the time derivative of z i using a centered difference scheme.

Computing the diffusive heat flux at the interface
Mixing results from a diffusive heat flux.Indeed, a purely advective heat flux displaces the surfaces without modifying their value.The diffusive heat flux occurs across, and normal to, the constant surfaces (since there cannot be any diffusive flux along those surfaces).
One way to compute the diffusive heat flux is to average the actual advective heat flux in space or time.The idea in doing so is that the oscillations due to reversible (wave) motions are filtered out by the averaging process and the residual non zero value gives the diffusive flux.However this method is not very precise because the residual value is usually very small relative to the maximum advective flux.An alternative method to access directly this residual diffusive contribution is provided by Winters et al. (1995).The principle of this method is to compute the hydrostatic equilibrium temperature profile associated with the minimum potential energy of the fluid at a given time.Conceptually this equilibrium state is reached by moving adiabatically and instantaneously the fluid particles towards their hydrostatic equilibrium position.Let s (z,t) be the temperature profile of this virtual equilibrium state, which is stable by construction.s evolves in time because of diffusive processes only and satisfies an equation of the form: ∂ t s = −∂ 3 ϕ d .The flux ϕ d is responsible for the temporal variation in s and is therefore the diffusive heat flux responsible for mixing.Hence, in the present context of interfacial mixing by convective motions, In practice, the stable temperature profile s (z,t) at a given time is computed by a simple adiabatic sorting of the temperature field at that time.More precisely the s profile is retrieved from each instantaneous field as follows.Let us consider a volume V, fixed with time, extending on both sides of the interface, from a level above the ground surface (at z/z i = 0.5) up to a level far above the upper boundary of the IL (at z/z i = 1.5).The instantaneous profiles are "sorted" over V, so that the fluid elements are moved adiabatically according to their value of , the lowest element Nonlin.Processes Geophys., 17,[187][188][189][190][191][192][193][194][195][196][197][198][199][200]2010 www.nonlin-processes-geophys.net/17/187/2010/ being the coldest 2 .Then, the change in time of the resulting sorted profile gives access to the diffusive heat flux.As shown by Winters et al. (1995), ϕ d has the following theoretical expression where I denotes an average along a -surface and κ t is the SGS turbulent diffusivity.Note that ϕ d is negative as expected.If pure laminar diffusion occurs, Eq. ( 17) reduces to the standard flux-gradient relation ϕ d,lam = −κ ∂ 3 s , where κ is the molecular diffusivity.Equation ( 17) therefore provides an expression to compute the turbulent diffusive flux ϕ d at any time.
To examine whether the horizontally averaged heat flux of the resolved scales F = w is a good approximation of the diffusive heat flux, we compare its value with that given by Eq. ( 17).For consistency with the definition of the interfacial heat flux, we should compare F and ϕ d at the altitudes where they each reach a minimum value, namely at z = z * i , say, for ϕ d and at z = z i for F. If F is a good approximation for ϕ d , these minimum values as well as z * i and z i should be very close.
The vertical profiles of F and ϕ d are compared in Fig. 6 at 12:00 EST and 13:30 EST.The flux ϕ d is negative, by definition, and is slightly smaller than F: the minimum value of ϕ d is 5% smaller than that of F at 12:00 EST and 16% smaller at 13:30 EST.The altitude z * i where ϕ d reaches its absolute minimum is 20 m lower than z i , the relative difference in altitude being less than 2%.This shows that w (z i ) is a very good approximation for the diffusive heat flux at the interface.In the following we take ϕ d z * i as the reference value for this flux.In other words, we define F i by ϕ d z * i .We now compare the values of w (z i ) and ϕ d z * i scaled by the surface heat flux for the times displayed in Table 1 (see Fig. 7).The constant value -0.2 is also indicated since a commonly used closure for F i is that it is proportional to F s with an empirical -0.2 coefficient.(The heat flux based upon the Lagrangian turbulent diffusivity, discussed in the next section, is also displayed in Fig. 7.) Figure 7 shows that the good agreement found between the two fluxes in Fig. 6 holds at all times, regardless of the value of the Richardson number, the relative difference ranging between 3% and 2 Incompressibility is assumed in the sorting method and we checked that this assumption is verified here.Indeed, the vertical displacements of fluid particles in the sorting process are at most equal to the thickness of the interfacial layer, that is 250 m or so.Since the sorting process is adiabatic, the change in the volume V of the fluid particles before (state 1) and after (state 2) sorting can be estimated by writing that p 1 .V γ 1 = p 2 .V γ 2 , where γ = 1.4 is the heat capacity ratio.If one assumes that the pressure is dominated by its hydrostatic component, one finds that the change in volume of the fluid particles during the sorting process is at most 3%.20%.Note that ϕ d z * i /F s varies by at most 9% during the 4 h of simulation reported here while w (z i )/F s changes twice more.Hence, not surprisingly, the diffusive heat flux is much less sensitive to large scale fluctuations than the advective heat flux.Figure 7 also shows that the simple closure F s = −0.2F s is an acceptable lower bound of the diffusive heat flux at the IL during this period of time.

Computation of the turbulent diffusivity from ϕ d
A turbulent diffusivity K ϕ can be inferred from the turbulent diffusive heat flux ϕ d , namely Note that, using Eq. ( 17) for ϕ d , K ϕ can also be expressed directly in terms of the temperature field.If the scale of the vertical gradient of s is much larger than the turbulent overturning scale, relation ( 18) is linear i.e.K ϕ is (nearly)  18) at z = z * i , using the turbulent diffusivity K L averaged over the interfacial layer (×); by −0.2 F s (dashed line).
uniform in z (e.g.Gregg, 1987).Under this condition, the turbulence may be assumed locally homogeneous.
The turbulent overturning scale within the IL is usually quantified by the buoyancy length scale b , defined by the ratio of the rms fluctuating vertical velocity σ w and the buoyancy frequency (e.g.Hopfinger, 1987).b is the largest vertical distance a fluid particle can move in a stably-stratified fluid against the potential temperature gradient.We computed b for the times indicated in Table 1, using the rms vertical velocity with the mean referring to an average over the IL (see Table 1) and using the buoyancy frequency defined as N B = (g β /δ 1 ) 1/2 .We found that b varies between 25 and 58 m, with a mean value of 41 m.This is consistent with the analysis of the IL dynamics by Otte and Wyngaard (2001), which yields b ≈20 m for conditions close to our LES (see their cases 19 to 22, in which is twice stronger than in the present case, all parameters being otherwise comparable).The values of b that we found have to be compared with the length scale associated with the mean vertical gradient of , which is δ 1 .Since δ 1 ≈225 m in our LES (see Table 1), local homogeneity may be assumed.
We computed K ϕ from our LES during the regime of equilibrium entrainment, which lasts from 12:00 EST to 13:30 EST as the ground surface heat flux is nearly constant during this period (see the values of F s in Table 1).During this period, the values of K ϕ obtained from Eq. (18) (and averaged over the IL) are between 3.52 and 4.17 m 2 s −1 depending upon time, and average 3.8 m 2 s −1 .In terms of SGS turbulent diffusivity κ t , we found that K ϕ is in the range 10-25 κ t , implying that the IL is turbulent.

Computation of the turbulent diffusivity from the dispersion of particles
An alternative method, based upon the dispersion of fluid particles within the IL, can be used to retrieve the turbulent diffusivity.This diffusivity will be denoted K L hereafter to make it distinct from that computed from ϕ d (though we expect K L K ϕ ).By "fluid particles", as usual, we mean non buoyant particles, which are advected by the velocity field (see Eq. 9).Let (δz) ms (t) be the mean square vertical displacement of fluid particles at time t for a given release of a particle cloud.(δz) ms (t) is defined by where N p is the number of particles of the release, z n (t) is the vertical position of the particle n and z G (t) the vertical position of the center of gravity of the particle cloud at time t.If the turbulence is locally homogeneous and stationary, and for t ≥ 2T L , with T L being the Lagrangian time scale of the turbulence, K L can be inferred from the growth rate of (δz) ms (see Taylor, 1921;Hunt, 1985, for a review), namely Since the IL is continuously forced by the quasi-stationary convective cells, the turbulence within this layer may be assumed stationary.In this case, the Lagrangian time scale T L is in the same order of magnitude as the Eulerian time scale T E (e.g.Hanna, 1981;Yeung, 2002;Dosio et al., 2005).This result is valid also in the presence of a stable stratification (Hunt, 1985).Let us assume that T E = 2π/N B .Hence, T E ≈40 s implying that 2T L is in the order of 1 min.Particles were released for z−z i =±100 m, that is, within the bulk of the IL.Note that some of the particles were released below and above the IL since its thickness varies over a wide range within the computational domain.The releases were made at 4 equally-spaced times from 11:55 EST to 13:25 EST over 10-min periods and resulted in a total of 57 500 particles per release.As an example, the time evolution of (δz) ms for the release carried out around 12:00 EST is displayed in Fig. 8.A quasi-linear growth occurs after about 1 min, whose growth rate is 2 K L according to Eq. ( 20).Values for K L between 3.24 and 3.83 m 2 s −1 were obtained depending upon the time of the release and average 3.6 m 2 s −1 .This range of values for K L is in very good agreement with that computed for K ϕ from the diffusive heat flux, the relative differences being lower than 10% on average.This is attested in Fig. 7 where −K L ∂ 3 s scaled by F s is displayed at altitude z * i for the times reported in SGS turbulence is likely to contribute to dispersion within the IL where small-scale turbulence dominates and T L is rather small.Thus, we conducted a simulation with a single release at 11:55 EST and switched off the SGS contribution in the Lagrangian stochastic model.In this case, (δz) ms was found to increase less rapidly during 1−2T L but reached a quasi-linear regime with the same slope (not shown), giving the same value of K L .Therefore the SGS contribution appears to play a negligible role in dispersing the particles.This result is consistent with the findings of Gopalakrishnan and Avissar (2000), and Cai et al. (2006) for a passive tracer.

Evaluation of the "first-order" model
Within the framework of the "first-order" jump model proposed by Betts (1974), the entrainment heat flux at the interface F i is related to the entrainment rate by where δ i = z 2 − z i and * = (z i ) + (z i + δ) /2.In the limit of infinitely small thickness of the IL, i.e. δ i →0, Eq. ( 21) reduces to the "zero-order" approximation for the interfacial heat flux derived by Lilly (1968), namely −F i = w e .Our purpose in this subsection is to evaluate the first-order model for w e by comparing the LES computation of w e from its definition (namely d t z i ) with its prediction by Eq. ( 21) using F i = ϕ d z * i .In order to compare also with the experimental data of Deardorff et al. (1980) we rather consider the parameter A = (w e /w * )Ri B instead of w e .This parameter is plotted in Fig. 9 versus time (see Fig. 9a) and versus Ri −1 B (see Fig. 9b) for the times reported in Table 1.The convection tank measurements of Deardorff et al. (1980) are included in Fig. 9b.The two quantities w * and Ri B are computed from the LES.A very good agreement is obtained between the LES values of w e and its prediction by the first-order model, the relative difference for A ranging between 2% and 18% with a mean value of 8%.(The relative difference for A with −F i computed as w (z i ) averages 11%.)The parameter A is in the range 0.21-0.26and averages 0.24, which is in good agreement with values reported in previous studies.This also shows that the standard parameterization w e /w * = 0.2 Ri −1 B for the entrainment rate is consistent with the present analysis of mixing.The coefficient A was shown to be equal to the efficiency of the mixing process by Chemel and Staquet (2007).
As pointed out by Fedorovich et al. (2004), the computation of the different terms in a given model should be consistent with the model order.The "first-order" model relies upon the finite thickness of the IL.Replacing F i by ϕ d , whose computation involves the depth of the CBL through the sorting process, may not fulfill this consistency condition.However, the temperature profile being (quasi-) uniform in the mixed layer and stable above the IL, this sorting process involves actually only the thickness of the mixed layer (see Fig. 6).Hence, estimating F i by ϕ d z * i in Eq. ( 21) is consistent with a "first-order" model.

Expression of the "first-order" model in terms of the turbulent diffusivity
With F i = ϕ d z * i , Eq. ( 18) becomes Approximating ∂ 3 s z * i by /δ 1 , the "first-order" model (21) becomes where K = K ϕ or K L .The approximation ∂ 3 s z * i /δ 1 should be discussed.The relative differences between ∂ 3 s z * i and ∂ 3 (z i ) were found to be less that 5% in our simulation.The relative differences between ∂ 3 (z i ) and /δ 1 range from 6% to 27% while those between ∂ 3 (z i ) and /δ i range from 14% to more than 100%.Hence, /δ 1 is a better approximation of ∂ 3 (z i ) than is /δ i .It is worth noting that by introducing δ 1 in Eq. ( 23), we extend the "first-order" model beyond that proposed by Betts (1974).Equation ( 23) is actually a mixture between what Sun and Wang (2008) have called the models "FOM1" and "FOM2", which differ only in the definition of the IL thickness (equal to δ i and to δ 1 , respectively).With this expression for w e and K = K L , the parameter A = (w e /w * )Ri B is plotted in Fig. 10   definition (w e = d t z i ).A very good agreement is found, the relative difference being lower than 10% on average.Also plotted in Fig. 10 are the results from a simple expression for w e , namely It is well-known (Sullivan et al., 1998) that the term δ i ∂ t * / in the "first-order" model is not negligible compared to −F i / in the range of Ri B values considered in our LES, its contribution here being up to 40% for the lowest Ri B values.However, Fig. 10 shows that, when the stratification is strong enough (Ri B approximately larger than 15 according to our data), the simple expression (24) accounts for the actual value of the entrainment rate to better than 25%.  1 as a function of Ri −1 B .w e is computed by three methods: from the LES (using its definition d t z i ) ( ); from Eq. ( 23), with K = K L averaged over the interfacial layer (×); by the simple model (24) ( ).Convection tank measurements of Deardorff et al. (1980) are included for comparison (•).The filled area represents A in the range 0.2-0.3.

Concluding remarks
In the present paper, the entrainment at the top of the convectively-driven boundary layer is reexamined using data from a high-resolution LES initialized by a commonly used sounding of Day 33 of the Wangara experiment and the analysis of mixing proposed by Winters et al. (1995).Note than an analysis along the same lines was conducted by D 'Asaro et al. (2002) for the oceanic convective mixed layer.
The mixed layer turbulence which forces the IL is first analysed in the present case of a "realistic" initialization.We found that the turbulence follows precisely the Kolmogorov spectral law for the velocity field over almost a decade in the inertial range.The multiplicative constant in this law is found to be in good agreement with previous measurements in the atmosphere.The Kolmogorov spectral law also holds for the frequency spectrum, when the Taylor's frozen turbulence hypothesis is used.To our knowledge, this is the first time that this hypothesis is verified properly in the context of the atmospheric boundary layer.Hence, the turbulence within the mixed layer may be assumed to be locally homogeneous and isotropic over a broad range of scales in the inertial range.
This turbulence forces and mixes the IL at the top of the convective layer.The parameterization of the heat flux at the IL, which is responsible for mixing, and of the resulting entrainment rate has been the subject of intensive research since Lilly (1968).We showed that the heat flux at the IL can be computed exactly from the analysis of Winters et al. (1995).The exact expression of this flux is denoted ϕ d .We defined the heat flux at the interface, usually referred to as F i , by the minimum value of ϕ d (consistent with entrainment models Nonlin.Processes Geophys., 17,[187][188][189][190][191][192][193][194][195][196][197][198][199][200]2010 www.nonlin-processes-geophys.net/17/187/2010/ in the atmospheric context) and we denoted z * i the altitude at which this minimum value is reached.This allowed us to show that the standard closure for F i , namely the minimum value of the horizontally averaged advective heat flux, agrees well with ϕ d z * i , to about 10%.The exact computation of F i along with a properly defined temperature profile within the interface (namely through a "sorting" process, following Winters et al., 1995) naturally yields a turbulent thermal diffusivity.The values of this turbulent diffusivity were recovered from the dispersion of fluid particles within the IL, which were tracked by a Lagrangian stochastic model coupled with the LES.The values thus derived agree indeed to better than 10% on average with those computed from ϕ d (whether SGS turbulence is included or not in the Lagrangian stochastic model).
These different estimates for the interfacial heat flux F i were next applied to the parameterization of the entrainment rate w e within the framework of the "first-order" model.This model basically relies on the thickness of the IL (as opposed to the "zero-order" model, where this thickness is infinitely small) and provides an expression for w e involving both F i and the thickness of the IL.We examined different predictions for w e from this model, depending upon the estimate for F i , which we compared with the LES value of w e .Overall, whether F i is computed from its exact expression ϕ d z * i , from its approximation using the horizontally-averaged advective heat flux or when the Lagrangian turbulent diffusivity is introduced, the prediction of w e by the "first-order" model agrees to about 10% with that computed from the LES using its definition (the best agreement being found for F i = ϕ d z * i ).A simple expression was also proposed for the entrainment rate, for which w e is equal to the Lagrangian turbulent diffusivity divided by the IL thickness.We showed that the values thus obtained differ from the LES values by 25% for strong enough stratification only (this relative difference being larger otherwise).However, for this expression to be of any use, one needs to access both the IL thickness and the turbulent diffusivity.Remote sensing techniques can provide values of the IL thickness (e.g.Steyn et al., 1999;Cohn and Angevine, 2000).Measurement of the turbulent diffusivity is likely to be more difficult but, as discussed by Winters and D'Asaro (1996), its values can still be retrieved from finescale-resolving vertical temperature profiles.

Fig. 1 .
Fig. 1.(a) Typical vertical profiles of potential temperature Θ and heat flux F for the convectively-driven boundary layer.(b) Same as (a) for the 'first-order' model proposed byBetts (1974) (referred to  as 'FOM1' in  § 4.4).All parameters are defined in the text.

Fig. 1 .
Fig. 1.(a) Typical vertical profiles of potential temperature and heat flux F for the convectively-driven boundary layer.(b) Same as (a) for the "first-order" model proposed byBetts (1974) (referred to as "FOM1" in Sect.4.4).All parameters are defined in the text.

Fig. 2 .
Fig. 2. Visualization of the structure of the boundary layer using potential temperaturecontours in a (x,z) plane located in the vicinity of an updraft at 15:00 EST.The distances along x and z are scaled by the domain length L and the mixed layer depth z i , respectively.The grayscale color table indicates variations at the interface (lower and higher appear white).The profiles, measured during the Wangara experiment (•) and computed from the LES results as a horizontally-averaged profile over the computational domain (-), are also included for comparison.

Fig. 3 .
Fig. 3. (a) Iso-surface Q=0.0015 s −2 for an isolated convective cell at 15:00 EST.The displayed view encompasses a horizontal domain of about 1.5 km×1.5 km, the height ranging from z=0 to z i .(b) Contour plot of the dimensionless vertical velocity w/w * at the same time in the horizontal plane z=z i /2.At the displayed time, w * =1.70 m s −1 and z i =1300 m.The contour lines correspond to w/w * =±1, ±2, ... Solid and dashed lines represent positive and negative contour values, respectively, with darkest zones being associated with updrafts and lightest with downdrafts.The distances along x and y are normalized by the domain length L. The structure visible in Fig. 3a has a horseshoe shape, which is clearly visible in Fig. 3b for approximately 0.2≤y≤0.6 and 0.2≤(−x/L + 1)≤0.6.

Fig. 4 .
Fig. 4. Relative spatial coverage of updrafts as a function of z/z i , as computed from the surface occupied by positive values of the vertical velocity.

Fig. 5 .
Fig. 5. (a)One-dimensional longitudinal velocity spectra E ii (k i ) of u and v computed for the 256 3 resolution run (-) at 15:00 EST and averaged over the range 0.4<z/z i <0.6.The spectra computed for a 128 3 resolution run (---) are superimposed for comparison.The dotted line (• • •) represents the spectrum deduced from the frequency velocity spectrum S u (f ) of u, displayed in plot (c) and computed for the 256 3 resolution run, from 12:00 EST to 15:00 EST at z=500 m in the center of the (x,y) plane.(b) Constant C 1 in Eq. (15) computed for the 256 3 resolution run and averaged for the u and v spectra as a function of k i .

Fig. 6 .
Fig. 6.Vertical profile of the heat fluxes w (-) and ϕ d ( ---) scaled by the surface heat flux F s at 12:00 EST (a) and at 13:30 EST (b), as a function of the vertical coordinate z scaled by z i .The filled area represents the interfacial layer (IL).

Fig. 9 .
Fig. 9. (a) Dimensionless parameter A = (w e /w * )Ri B for the different times displayed in Table1, with w e computed by two methods: from the LES (using its definition d t z i ) ( ); from the firstorder model with F i = ϕ d z * i (•).(b) same as (a) except that (w e /w * ) Ri B is now plotted as a function of Ri −1 B .Convection tank measurements ofDeardorff et al. (1980) are included for comparison (•).The filled area represents A in the range 0.2-0.3.

Fig. 10 .
Fig. 10.Dimensionless parameter A = (w e /w * )Ri B for the different times displayed in Table1as a function of Ri −1 B .w e is computed by three methods: from the LES (using its definition d t z i ) ( ); from Eq. (23), with K = K L averaged over the interfacial layer (×); by the simple model (24) ( ).Convection tank measurements ofDeardorff et al. (1980) are included for comparison (•).The filled area represents A in the range 0.2-0.3.

Table 1 .
Characteristics of the convective boundary layer for the 256 3 resolution run.