A large-eddy simulation model is adopted to investigate the evolution of scalars transported by atmospheric cloud-free convective boundary layer flows. Temperature fluctuations due to the ground release of sensible heat and concentration fluctuations of a trace gas emitted at the homogeneous surface are mixed by turbulence within the unstable boundary layer. On the top, the entrainment zone is varied to obtain two distinct situations: (i) the temperature inversion is strong and the trace gas increment across the entrainment region is small, yielding to a small top flux with respect to the surface emission; (ii) the temperature inversion at the top of the convective boundary layer is weak, and the scalar increment large enough to achieve a concentration flux toward the free atmosphere that overwhelms the surface flux. In both cases, an estimation of the entrainment flux is obtained within a simple model, and it is tested against numerical data. The evolution of the scalar profiles is discussed in terms of the different entrainment–surface flux ratios.
Results show that, when entrainment at the top of the boundary layer is weak, temperature and trace gas scalar fields are strongly correlated, particularly in the lower part of the boundary layer. This means that they exhibit similar behavior from the largest down to the smallest spatial scales. However, when entrainment is strong, as moving from the surface, differences in the transport of the two scalars arise.
Finally, it is shown that, independently of the scalar regime, the temperature field exhibits more intermittent fluctuations than the trace gas.
Surface–atmosphere interactions are key issues for environmental and climate modeling, since surface fluxes and scalar transport in the lower atmospheric boundary layer provide boundary conditions for numerical models of weather prediction and general circulation. This has motivated an extensive experimental, theoretical, and numerical research activity about scalar fluctuations and related fluxes in atmospheric boundary layers, at varying thermal stability and surface conditions (see, e.g., Wyngaard 1990; Mahrt 1991; Jonker et al. 1999; Van Dop et al. 2005; Huang et al. 2009). More generally, turbulent transport and mixing of scalars is of fundamental relevance for a broad variety of systems (Dimotakis 2005).
Active scalars affect the flow dynamics; this is the case for temperature fluctuations in convective or stably stratified boundary layers. Passive scalars, on the contrary, do not modify the flow field; this is the situation for transport of carbon monoxide or, for example, that of the humidity field under neutral conditions. In the atmospheric boundary layer, in addition to the physics of turbulent transport, both surface and entrainment physics are crucial to scalar mixing. Scalar fields are commonly emitted close to the ground and are dispersed by turbulence within the mixed layer. The capping inversion acts as a lid to the rise of thermals and to diffusion of scalars in the free atmosphere. On the other hand, entrainment is the mechanism through which stably stratified air is mixed from above into the unstable boundary layer. It can represent a source or a sink of scalar fluctuations depending on the scalar increment across the temperature inversion region. The relative weight of the entrainment to surface emission may lead to different regimes, thus the importance of a reliable estimation of entrainment fluxes from the free atmosphere (see, e.g., Stull 1976; Deardorff 1979; Deardorff et al. 1980; Sullivan et al. 1998; Huang et al. 2011; Fedorovich et al. 2004). The strong dependence of scalars on boundary conditions is also responsible for the failure of simple parameterizations based on bulk quantities of the mixed layer scales only (Deardorff 1970).
In several works, the role of boundary constraints was accounted for by considering scalar mixing as a superposition of diffusion originated from the surface “bottom up” and from the entrainment zone “top down” processes (Wyngaard and Brost 1984; Wyngaard 1984; Moeng and Wyngaard 1984; Wyngaard 1990; Sorbjan 1999; Moene et al. 2006). The current understanding is that the statistics of the two processes evidently depart from each other, and are regulated by scalar fluxes or gradients evaluated at the (top or lower) boundaries. A comprehensive investigation of the upper portion of the mixed layer during both free-convection and forced-convection regimes [see Sorbjan (2005) and Sorbjan (2006), respectively] suggests that statistical moments at the top of the mixed layer can be parameterized in terms of interfacial scales, defined through peak values of the temperature and passive scalar gradients in the entrainment region (Sorbjan 2004) and on the velocity gradients at the interface. In addition, the strengths of the capping inversion and of the geostrophic wind are shown to have a deep impact on the passive scalar profile and on its trend in the entire boundary layer.
From a statistical point of view, it is important to characterize typical fluctuations at the large scale of the system. Such fluctuations do enter turbulent fluxes and similarity functions. It is also important, however, to characterize bursty fluctuations taking place on short time lags or small spatial scales. Understanding and modeling the statistics of intermittent fluctuations is of obvious importance if one considers, for example, the probability of having a local spatial fluctuation of pollutant gas concentration exceeding some tolerable level as it is emitted from the source. Equally relevant is the role of concentration fluctuations much larger than mean ones in controlling the rate of slow chemical reactions—for example, in the process of atmospheric ozone destruction, as was suggested in Edouard et al. (1996).
Near the boundaries, the organization of the turbulent fields is strongly nonuniversal: momentum transfer is dominated by shear at the bottom and by entrainment at the top. This has consequences on the scalar transport also. Indeed, turbulent structures can be responsible for the flux imbalance of heat and a bottom-up tracer, while the same relations do not hold for top-down diffusion (J. Huang et al. 2008). Far from the boundaries, however, turbulent scalar transport should possess more universal features. For example, ramp-cliff patterns in the scalar field have been observed in a variety of conditions (Tong and Warhaft 1994; Shraiman and Siggia 2000; Warhaft 2000). Their presence is not the footprint of the coherent structures of the vorticity field or of the shear, but is an inherent feature of the turbulent nature of the scalar field and scalar mixing (Shraiman and Siggia 2000). However, their spatial organization can vary depending on the turbulent velocity field and on the presence of anisotropy and inhomogeneities in the system [see the discussion in Gao et al. (1989), Paw et al. (1992), and references therein]. Ramp-cliff patterns in passive fields or sharp edges in the turbulent thermal plumes are associated with small-scale frontal structures of the scalar field, generally present in any turbulent flow; in a nutshell, large spatial homogeneous regions for the scalar field are separated by thin regions where very intense scalar gradients are located. These localized violent fluctuations are responsible for the large deviations from mean values, and leave also a statistical fingerprint in the strongly non-Gaussian behaviors of scalar fields (Frisch 1995; Pope 2000; Shraiman and Siggia 2000; Mazzitelli and Lanotte 2011; Tong and Warhaft 1994; Warhaft 2000).
In this paper, by means of large-eddy simulations, we focus on similarities and differences in the statistics of a passive trace gas field, and the temperature active scalar field, in cloud-free convective atmospheric boundary layers (CBLs), by varying the structure of the inversion layer. We detail typical fluctuations, in terms of the vertical profiles of mean potential temperature and mean concentration, and profiles of turbulent fluxes. For this, it is relevant to account for the role of the entrainment, for which we suggest a simple parameterization. Entrainment and mixing with clear warm air is however also responsible for differences in the scalars’ turbulent fluctuations at smaller scales and for front structures. These are quantified by studying the statistical moments of scalar fluctuations of high order. While most works focused on scalar variances and covariances, there are few results about the statistical signature of small-scale fronts in boundary layers (see, e.g., Antonelli et al. 2003). The interest is twofold. First, atmospheric dispersion is one case study of scalar turbulence at high Reynolds and Péclet numbers (Warhaft 2000). Second, it is also evident the relevant practical importance of being able to quantify the probability of extremely large spatial fluctuations, sampled by high-order statistical moments of scalar increments.
Laboratory experiments (Zhou and Xia 2008) investigated active and passive scalar fields in turbulent Rayleigh–Bénard convection at high Rayleigh numbers. They found quite different statistics for the local concentration of a fluorescent dye and the local temperature. In particular, the latter was found to be much more intermittent and anisotropic than the former. Mazzitelli and Lanotte (2011) presented an analysis of the intermittency of both temperature and tracer fields in CBLs. Results confirm what is found in laboratory experiments—namely, that the temperature is more intermittent than the trace gas field. In section 2, the numerical method is applied and the simulation’s parameters are discussed. A complete description of the dynamic model adopted to describe subgrid-scale fluctuations is given in the appendix. In section 3, the CBL characteristics, mean profiles, and fluxes are presented. The entrainment process and the different regimes that may develop are investigated. We also discuss a simple parameterization of the entrainment flux, and compare it with the numerical data. For two selected case studies, in section 4, we contrast the behavior of temperature and trace gas fluctuations, and in particular, we discuss the role of small-scale scalar fronts. We quantify their statistical weight and highlight differences occurring between temperature and concentration fields. The last section contains concluding remarks.
2. Large-eddy simulations
Large-eddy simulations have become a standard tool for numerically studying atmospheric boundary layer dynamics (Sagaut 2006). Particularly, they have the capability of correctly describing the atmospheric boundary layer physics in a convective regime (Deardorff 1972, 1974; Moeng 1984; Mason 1989) where the energy transfer is dominated by the almost self-similar turbulent energy cascade (Pope 2000). In addition, according to Kolmogorov theory of isotropy restoration at small scales, turbulent eddies are expected to become more and more homogeneous and isotropic when going to smaller and smaller scales (Frisch 1995; Biferale and Procaccia 2005; Antonelli et al. 2007), where they are ultimately dissipated by molecular viscosity or diffusivity in the case of scalar fluctuations. In a nutshell, the application of the large-eddy simulation approach is based on self-similarity and on the restoration of local symmetries [see, e.g., Meneveau and Katz (2000) for a review on the subject]. The extent of validity of the local isotropy restoration has been thoroughly investigated in turbulent flows (Biferale and Procaccia 2005) and in the specific context of the convective boundary layer also (see Antonelli et al. 2007; Chen and Tong 2006).
The incompressible Navier–Stokes equations, with Boussinesq approximation, for the velocity field u(x, t), and the advection–diffusion equations for the potential temperature θ(x, t) and trace gas concentration c(x, t) scalars are integrated. The large-eddy simulation (LES) model equations, resulting from low pass filtering of the previous governing equations, are
Here the indexes i, j are running over x, y, z (ux = u, streamwise; uy = υ, spanwise; uz = w, wall-normal component), and repeated indexes are retained summed. In the above equations, g is the acceleration due to gravity directed along z, θ′ is fluctuation with respect to the mean of the potential temperature, θ0 is a reference potential temperature, fc is the Coriolis parameter, and Ugj is the j component of the geostrophic wind; is the distortional part of the residual [or subgrid scale (SGS)] strain tensor τij, which is defined according to
The isotropic component of the strain (5) is included in the pressure term
with p the fluctuating component of the physical pressure and ρ0 the density of air. An advection–diffusion equation describes the time evolution of the resolved potential temperature and of the trace gas field ; however, while concentration is passive, temperature is actively being coupled to the velocity field via the buoyancy term .
The SGS stress for the scalar θ (or c) is defined as
Molecular viscosity and diffusivities have been neglected, assuming that, at high Reynolds (Peclét) numbers, and with the filter cutoff well inside the inertial range, any molecular effect is negligible with respect to turbulent dissipation. Closures of Eqs. (1)–(4) are obtained by modeling the SGS terms through the resolved fields. Most SGS models are based on the following linear eddy viscosity assumption:
with km the eddy viscosity to be modeled, and the resolved rate-of-strain tensor . To avoid the use of empirically determined values of the model coefficients, we adopted the dynamic subgrid-scale model (Germano et al. 1991); this requires the introduction of another filter scale, known as test filter scale. In the appendix, we detail the formulation of the dynamic model implemented in the LES model originally developed by Moeng (Moeng 1984; Moeng and Sullivan 1994).
The physical domain has horizontal dimensions Lx = 5000 m and Ly = 5000 m, while the vertical one is Lz = 2000 m. Simulations were carried out with two different numerical resolutions, with a mesh spacing of about 78 m × 78 m × 31 m along the x, y, and z, directions respectively, and with a mesh spacing of 20 m × 20 m × 8 m. The coarse resolution is used to study various instances for the time evolution of trace gas field, depending on the initial and boundary conditions. The high resolution with N3 = 2563 grid points is applied to compute high-order statistics in two selected cases; see below. In all the runs to be described below, a constant in time, homogeneous in space, surface heat flux H0 = 0.24 m s−1 K−1 is imposed, so that a quasi-steady state is reached for the potential temperature and velocity fields. Similarly, a constant in time, homogeneous in space, surface concentration flux Q0 = 0.2 ppm m s−1 guarantees a quasi-steady state for the concentration of the trace gas. A geostrophic wind Ug = 10 m s−1 is imposed in all runs. The initial parameters and numerical features common to all runs are summarized in Table 1.
3. Scalar evolution with different entrainment fluxes
Entrainment is the process of mixing of a weakly turbulent flow into a turbulent one. In the atmosphere, the engulfment of nonturbulent air parcels from the stable free atmosphere within the unstable mixed layer results in the deepening of the CBL. Such air parcels are efficiently mixed and transported by turbulent structures and, upon certain conditions, can penetrate down to the lower boundary. Recently, convective boundary layer entrainment has been reviewed in connection to Doppler lidar measurements (Träumner et al. 2011) and in LES studies (see, e.g., Fedorovich et al. 2004; Huang et al. 2011). While the distinction of different entrainment mechanisms can be useful for a physical understanding of the ongoing processes, parameterizations mostly rely on dimensional arguments and scaling relations of general validity. Since a first-principle theory is still lacking, the agreement between experimental measurements and parameterizations is only empirically determined and is often situation dependent.
Comparison between in situ airborne observations and large-eddy simulations results about entrainment has been discussed in Canut et al. (2012). Two types of parameterizations are evaluated: namely, zero-order jump models (ZOM) where entrainment zone is described as a sharp discontinuity, and first-order jump models (FOM) where the discontinuity is replaced by a finite-depth region. Results show that the FOM approach works generally better that ZOM in reproducing observed values for the entrainment velocity in the Sahelian boundary layers characterized by large scalar jumps and deep inversion layers.
In this section, we discuss a zero-order jump model for the intensity of the turbulent concentration flux across the entrainment zone and the influence of the modeled flux on the boundary layer temporal evolution. The reason for adopting a ZOM approach is that it allows a priori estimations. Moreover, provided that we have a good estimation of the scalar jump across the top layer, ZOM is expected to give a good approximation of the observed entrainment fluxes (see, e.g., Huang et al. 2011).
The rate of deepening of the mixed layer is quantified by the entrainment velocity we. The influence of large-scale advection translated into a synoptic-scale vertical velocity is explicitly neglected in our work. In the case of forced convection—that is, a convective boundary layer with imposed geostrophic wind—an estimation of we results from the energetics theory (Stull 1976):
where d1 is the difference between the boundary layer height zi and the height of zero heat flux, (ΔEZθ) is the average (over the horizontal direction) temperature difference over the entrainment zone, (ΔEZU) is the difference for the magnitude of the wind, w* = (g ziH0/θ0)1/3 is the convective velocity, and u* the surface friction velocity. In the above formula, empirical nondimensional constants are c1 = 0.0167, c2 = 0.5, and c3 = 0.0006 (Stull 1988). Assuming a linear profile for the heat flux H(z), the distance d1 can be evaluated as
The heat flux at z = zi is
under the assumption that the interfacial layer has infinitesimal depth (zero-order model; Lilly 1968). Equations (9) and (10) provide a closed system that we solved for the unknown we and d1, once temperature and velocity increments, and the boundary layer driving velocity scales, are specified. In the numerical experiments that we performed, such parameterization satisfactorily agrees with the numerically obtained value for ∂〈zi〉/∂t.
In the free-convection regime, the turbulence at the top of the boundary layer is directly related to the surface forcing and the entrainment heat flux is often assumed to be a fraction of the surface flux (i.e., Hi = −βH0), with the parameter β assumed in the range 0.1 < β < 0.3 (see, e.g., Tennekes 1973). Indeed, the application of Eq. (9), in the case of free convection (u* = 0, ΔEZU = 0), and the use of Eq. (10) or Eq. (11), yield respectively
For β ≃ 0.2 these coincide, and also agree with, the result obtained by Deardorff et al. (1980) in convection tank experiments, where the entrainment rate was connected to a bulk Richardson number Ri*, according to
On similar dimensional grounds, we propose the following parameterization of the entrainment flux of the trace gas at the top of the mixed layer:
with ΔEZc the scalar increment over the entrainment region. The angle brackets indicate the average over the homogeneous horizontal planes. With respect to Lilly (1968)’s jump relation, we explicitly include a proportionality constant χ, which may be predicted to be equal or larger than 1. Indeed, it is reasonable to expect that the concentration entrainment occurs at a faster rate than temperature entrainment [see Eq. (11)], as the latter one is resisted by buoyancy effects. The constant χ needs to be estimated empirically.
The intensity of Qi with respect to the surface concentration flux Q0 has a crucial relevance upon the boundary layer regime. From experimental and numerical studies on moist boundary layers (see, e.g., Mahrt 1991; Sorbjan 2005, 2006), it is well known that, depending on the degree of instability and on the strength of the surface moisture flux, an entrainment-moistening or -drying boundary layer will develop. It is the ratio Qi/Q0, between the top and the surface flux, that governs the boundary layer evolution in time. When Qi/Q0 is larger than one, the amount of dry air entrained overwhelms the concentration flux from the surface and the boundary layer dries in time. On the other hand, for Qi/Q0 smaller than one, the mixed layer tendency is toward higher moist mixing ratios. Since Eq. (15) requires the knowledge of the interfacial increments only, it allows us to estimate a priori the concentration trend in the boundary layer by making use of quantities of easy access for large-eddy simulation models, as well as for experiments. Other parameterization schemes for the entrainment variables based for instance on the gradients at the interfacial layer (e.g., Sorbjan 2005, 2006) are able to provide an analytic expression for the fluxes at all heights within the mixed layer, but make use of inner parameters of the CBLs, which are more difficult to estimate in advance. In the following, the parameterization proposed by Sorbjan 2006 for the passive scalar entrainment flux, under a forced-convection regime, will be tested against our numerical data and compared with Eq. (15).
We studied different regimes of scalar dilution or increased mixing ratio by restricting the parameter space, which is in principle multidimensional, to a two-dimensional one: we varied the entrainment velocity we and the scalar increment at the capping interface ΔEZc only. These two parameters account for all relevant mechanisms; by changing the entrainment velocity, different instability and shear effects are included, while by changing the scalar increment, we explore the effects of different entrainment zones.
Five large-eddy simulations with spatial resolution 78 m × 78 m × 31 m are performed for the analysis of the scalar field evolution and the scalar flux. In these runs, initial conditions are a constant profile for the trace gas concentration in the mixed layer and linear decrease in the entrainment zone. The scalar mixing ratio is set to zero in the free atmosphere. The simulations differ for the initial value of ΔEZc and for the strength of the capping inversion. The parameters for these runs are listed in Table 2. From the knowledge of we and of the nondimensional constant χ, it is possible to obtain an estimation of the threshold concentration jump between the mixed layer and the free atmosphere that discriminates between the two entrainment regimes: |ΔEZc|thr = Q0/(χwe). When |ΔEZc| > |ΔEZc|thr the entrainment flux Qi exceeds the surface flux Q0, and the boundary layer will dilute in time. On the other hand, when |ΔEZc| < |ΔEZc|thr, the opposite trend is likely to occur.
In Figs. 1 and 2, we report the evolution of concentration profiles in strong and weak inversion runs. The SGS fluxes and the total concentration fluxes are also plotted in the insets. They present the correct qualitative vertical dependence for entrainment-diluting and -increasing case, respectively. The plots confirm that it is the combined effect of entrainment rate and concentration jump that determines the regime. Consider, for example, Fig. 1: it shows that, even with a very strong capping inversion, large values of |ΔEZc| may lead to a regime of decreasing concentration by entrainment. Note also that in such case the total scalar flux exhibits a sharp maximum near the top of the boundary layer: this is the combined effect of the capping inversion and of the scalar jump ΔEZc.
It is instructive to consider a further numerical experiment with a strong inversion and initial trace gas concentration |ΔEZc| = 10 ppm. The outcome is presented in Fig. 3. As it reads from Table 2, the entrainment velocity increases in time owing to the more intense convection. Hence, whereas the concentration initially increases, later on the top scalar flux Qi overpasses Q0 and a boundary layer–diluting regime settles in. This run permits the estimation of the constant χ, introduced in Eq. (15): the numerically obtained value belongs to the range 1.5 < χ < 2.2. We put χ = 1.7 ±0.3, in order to get good agreement between the numerical results and the theoretical estimation, during the first large-eddy turnover times of the run (see Fig. 4). By using such value to estimate the threshold concentration jump |ΔEZc|thr, the results obtained from all simulations are consistent with the predicable tendency of concentration in the mixed layer. Moreover, this value for χ is robust, since the runs at higher resolutions, presented in next paragraphs, are consistent with the same value. This suggests that it is possible to distinguish a priori between different flow regimes by means of a simple parameterization for the entrainment flux of the passive scalar, like Eq. (15).
We can now compare the estimated value for the top trace gas flux by means of Eq. (15), with the value measured in the simulations. In Fig. 4, it is plotted the time evolution of the ratio Qi/Q0, evaluated in the run with strong inversion, and with initial scalar concentration |ΔEZc| = 10 ppm. The results are also compared to the parameterization proposed by Sorbjan 2006 [Eq. (9b) in their paper], which is calculated starting from our results. It appears that the two models correctly capture the trend, but overestimate the increase of the flux ratio Qi/Q0, with comparable accuracy. Only for the short time evolution (i.e., for t/τ* ≤ 16), predictions of both models are rather good. However, we remark that taking into account the uncertainty associated to the value χ = 1.7 ±0.3, the predictions of model (15) agree with observations within the error bars.
We have performed different runs at changing the resolution (see below) and at varying the physical parameters: we generally observe that at short time there is a good agreement between the model (15) and the observations, which, however, deteriorates at later stages. Hence, we conclude that a simple model like the one proposed is not fully able to quantitatively capture the entrainment dynamics when nonlinear or saturation effects become important.
4. Temperature versus concentration statistics in two case studies
A detailed comparison of temperature and concentration fluctuations, from large down to smaller, inertial range scales, is carried out by selecting two case studies. These correspond to entrainment-diluting and -increasing concentration CBLs, respectively, and are realized by varying the thermal inversion and the concentration jump at the top of the boundary. For this purpose, we performed high-resolution large-eddy simulations. The initial parameters of the runs are collected in Table 1. The internal parameters characterizing the simulations are listed in Table 3. For both runs, the quasi-stationary state for the CBL is approached within six large-eddy turnover times, when the resolved turbulent kinetic energy stops to oscillate (Nieuwstadt et al. 1991). At this point the passive concentration field is injected in the system, via a constant in time, homogeneous in space surface flux equal to Q0 = 0.2 ppm m s−1. Quasi-steady equilibrium for the passive scalar settles in after six large-eddy turnover times since the injection. This is tested by computing, on horizontal planes, the probability density function (PDF) of the resolved scalar concentration and checking the time independence of its functional shape. Statistics are then accumulated for 20τ* by saving 80 realizations of the model fields, equally spaced in time. Results presented in the following are obtained by averaging over all realizations.
In the strong inversion run labeled S256, we impose an initial concentration increment |ΔEZc| = 0.4 ppm across the entrainment zone, and we then observe a mean concentration profile increasing with time. Under these conditions, the entrainment velocity we is low, owing to the strong stability at the inversion, and the initial scalar top flux Qi ≃ 0.006 ppm m s−1—estimated by Eq. (15) with χ = 1.7—is extremely small with respect to the surface concentration flux Q0. Differently, in the run labeled W256, the entrainment-diluting boundary layer is obtained starting with a weak temperature inversion and a large concentration jump: |ΔEZc| = 10 ppm. The top scalar flux estimation is, in such a case, Qi ≃ 0.3 ppm m s−1, implying that entrainment of clean air exceeds the surface scalar flux (i.e., Qi > Q0). The time evolution of the concentration profiles is presented in Fig. 5; in the inset, the scalar vertical fluxes are also plotted.
It is interesting to remark the opposite effect of entrainment for active and passive scalars in the convective boundary layer. Typically, scalar concentration approaches zero above the mixed layer, whereas the temperature field presents a positive free atmosphere lapse rate. Therefore, in convective boundary layers, entrainment of passive scalar tends to dilute the internal concentration, whereas entrainment of active scalar introduces positive temperature fluctuations. Accordingly, the concentration flux Qi at the top of the mixed layer is positive, while the heat flux Hi is negative. In other words, the sign of the scalar vertical gradients is the signature of the physics underneath. Close to the surface, negative gradients indicate the decrease of scalar fluctuations with height, thus a scalar source from below. On the contrary, at the top of the mixed layer, gradients are negative or positive depending on whether scalar fluctuations decrease or increase with height (i.e., on whether the top boundary behaves as a scalar sink or source, respectively). From the inspection of the derivative skewness plots (not shown), it appears that similar spatial structures of temperature and concentration are observed in the lower part of the boundary layer, signaling a strong correlation at small scales between the two scalar fields (Mazzitelli and Lanotte 2011). In this respect, in Fig. 6, we plot the skewness of the passive scalar fluctuations:
averaged over the horizontal planes (denoted by the angle brackets). A skewness larger than zero indicates that positive fluctuations are more intense than negative fluctuations. This is indeed the case for narrow rising hot plumes in the mixed layer. In our setup, the passive scalar is emitted at the surface: similarly to what has been observed for the temperature in convective boundary layers, it tends to collect in the low velocity regions and rises inside the convective plumes. This phenomenology leads to a positive skewness for the concentration field. On the other hand, entrainment effects cause negative concentration fluctuations at the top of the mixed layer, which are well captured by the skewness negative sign (see, e.g., Mahrt 1991). In the competition between surface emission and entrainment, in run S256, the first effect dominates upon the second, and concentration skewness is positive in the lower part of the boundary layer. On the contrary, in run W256, the negative fluctuations originated at the top are capable of reaching the surface layer. Skewness is negative at all heights, and an entrainment-diluting boundary layer rapidly develops.
a. Scalars’ correlation
Because they are turbulent, temperature and concentration display fluctuations over a wide range of scales from the large ones associated to the injection mechanisms to the small ones typically associated to small-scale frontal spatial structures. In Fig. 7, simultaneous horizontal snapshots of temperature and concentration fluctuations for runs S256 and W256 are presented. These are evaluated in the bottom part of the mixed layer, at z/zi ≃ 0.17. It is evident that, at this distance from the surface, similar large-scale structures of the two scalar fields show up in run S256, while the same is not true for run W256.
A priori, it is not obvious if two fields, correlated on the large scales, will continue to be correlated at intermediate and smaller scales, thus exhibiting similar statistical properties. Neither it is clear if differential entrainment can be responsible of breaking correlations among the two scalars. To assess to which extent the two fields are statistically different, we start by measuring the linear correlation coefficient between temperature and concentration resolved fields,
where the average is performed over the horizontal planes. The coefficient profiles, as a function of the normalized vertical coordinate z/zi, are plotted in the top panel of Fig. 8. In run S256, rθc is very close to 1 in the bottom part of the boundary layer, thus signaling the strong bound between the two fields. Anticorrelation between temperature and concentration occurs only in the upper part of the mixed layer. In simulation W256, the correlation coefficient profile indicates a lower degree of correlation between temperature and concentration at all heights in the boundary layer. In the top panel of Fig. 8, the same analysis is carried out for the scalar increments over the smallest resolved scale Δx ≃ 20 m along the x direction. Again, results for run S256 show high correlation in the lower part of the CBL. Similar plots (not shown) were obtained when considering directions y and z. We can conclude that both runs exhibit on average similarities between temperature and concentration, with spatial correlations from the large down to the smallest turbulent resolved scales, and whose importance is modulated by entrainment effects. Such result can be useful in parameterization of turbulent fluctuations.
b. Intense fluctuations of the scalar fields: Definition and statistical relevance
Turbulent scalar fields organize themselves in spatial structures—for example, in the form of large-scale patches. Within these, scalar fields are homogeneous and do not significantly vary. Such structures are separated by thin edges where sharp variations take place: these are the so-called small-scale fronts or cliffs. Such turbulent spatial structures at small scales can be investigated by looking at the statistics of the scalar increments measured over a distance r: . In an isotropic situation, the separation vector r can be taken along any direction. In the CBL, considering the symmetries of the system, we choose the separation vector r = (0; ry; 0) lying on a horizontal plane at distance z from the surface.
Temperature structures in the convective surface layer are plumes, characterized by a sharp trailing edge and a more diffuse leading edge (see, e.g., Antonia et al. 1979; Stull 1988). The presence of a plume, moving along the mean wind direction, can be identified by detection of a front with sharp temperature increase followed by a ramp with slow temperature decrease. In passive scalar turbulence, ramp-cliff structures are displayed by the scalar field, independently of the statistics of the advecting flow (Warhaft 2000). To study the behavior of intense scalar fluctuations in the boundary layer region, a useful observable is the probability of detecting a spatial increment larger than the standard deviation as function of the distance from the surface. This is reported in Fig. 9 for both runs S256 and W256. The cumulated probability curves show that the number of intense fluctuations is generally larger for temperature than for concentration, at all heights within the boundary layer, and more importantly that the lower part of the boundary layer is the region where the strongest fluctuations take place. In the run with the weak inversion, W256, the cumulated probability of observing very large concentration increments is smaller; in such a case, the dilution caused by the entrainment of clear air has smoothed down large scalar increments near the surface, thanks to an efficient mixing. The main result is that the region near the surface is mostly affected by these very intense and intermittent fluctuations; indeed, the scalar fields experience, from one point to another, excursions much larger than the standard deviation.
It is also interesting to examine whether the intense fluctuations of the two scalars are detected at the same spatial location or not. In Fig. 10, we plot the conditional probability of observing a large trace gas increment at the same location where we have observed a large temperature increment. It appears that there are regions where the majority of large concentration fluctuations are connected to large temperature fluctuations; this is the case for run S256, particularly in the bottom part of the CBL. This implies that here the active scalar is driving the most intense fluctuations of the passive scalar.
A much smaller value of the conditioned probability is observed in the case of run W256, where trace gas fluctuations are not carried within thermal plumes.
Such intense fluctuations of the scalar increments are connected to quasi singularities of the field, also dubbed cliffs. A cliff is observed wherever there is a large jump of the scalar field occurring across a small-length scale r (Celani et al. 2000, 2001):
For the present system, a visualization of the spatial location of the cliffs is given in Fig. 11, where we plot the union of the locations—in a horizontal plane—where resolved temperature increments at scale r are larger than :
with λ > 1, and similarly for the trace gas. The separation vector r is taken along direction y [i.e., (0, r, 0) = (0, 2Δy, 0)] to avoid the residual asymmetry due to the surface-layer plumes; the amplification factor is here chosen λ = 3.5 (see Antonelli et al. 2003). It appears that cliffs are aligned with the mean wind and break the symmetry between the x and y direction within the surface layer. The snapshots are taken at z/zi ≃ 0.17 so to test the physics of a region where concentration and temperature are strongly correlated.
It is clear that, in general, temperature displays a density of intense fluctuations higher than the trace gas. The plot referring to run S256 (left panel) shows that concentration cliffs, in the lower boundary layer region, are frequently associated to temperature jumps. Differently, in the case of run W256 (right panel), the trace gas cliffs are rare and their locations are uncorrelated from those of the temperature.
A prominent characteristic of passive scalar statistics is the so-called saturation of intermittency, which has been attributed to the presence of cliffs in the scalar field spatial organization (Celani et al. 2000, 2001). From a physical point of view, saturation means that the very large fluctuations of scalar increments are dominated by the occurrence of frontlike discontinuities. The statistical signature of saturation manifests in the tails of the scalar increment probability density functions (Celani et al. 2001; Antonelli et al. 2003), which show an algebraic behavior with a single exponent ζ∞. In Mazzitelli and Lanotte (2011), it is shown that both the most intense fluctuations of concentration and temperature show saturation of intermittency in the mixed layer of CBLs; however, the saturation scaling exponents are different and the temperature field appears to be more intermittent than the concentration one. The scaling exponent characterizing the tails of scalar increment PDFs also has an interesting geometrical meaning. It is related to the dimension of the fractal set DF where cliffs live (Frisch 1995). In particular, the following relation applies ζ∞ = 2 − DF; that is, the saturation scaling exponent is related to the fractal dimension of the set where cliffs stay (we recall that we measured cliffs only in the horizontal planes). We verified the validity of the above relation by measuring the fractal dimension of the fronts, and we obtained ζ∞ = 1.2 and ζ∞ = 0.6 for concentration and temperature, respectively, in agreement with what we measured from the algebraic scaling of the tails of the scalar increment probability density function (not shown). Fractal dimensions have been measured not too close to the surface to avoid wall effects and high spatial correlation (see Fig. 8). This behavior is observed in both runs S256 and W256, indicating that the different features displayed by two scalars statistics are robust with respect to changes in the capping inversion and surface fluxes.
The spatial statistics of temperature and trace gas concentration fields was compared in a numerical study of cloud-free convective boundary layers, driven by different surface and entrainment fluxes. To investigate different CBL regimes, we focused on the dependence of the entrainment fluxes on the entrainment velocity and the scalar jump across the entrainment zone. For the trace gas concentration, entrainment strength is due both to the intensity of the capping inversion and to the mean trace gas concentration in the CBL. We showed that a simple parameterization based on a zero-order model can be a useful tool to predict scalar fluxes, particularly when experiments do not allow for high-resolution spatial measurements. A reason why a zero-order model, such as the one proposed, works satisfactorily well is that in the CBLs here considered the entrainment layer thickness is small. This is in contrast with the boundary layers observed, for example, in Canut et al. (2012) during the African Monsoon Multidisciplinary Analysis experiment, where a thick entrainment layer arises owing to the presence of a southwesterly moist, cold monsoon flow topped by a northeasterly dry, warm flow. We remark that in our study, we did not consider the role of wind shear, although it is certainly important (see, e.g., Sorbjan 2006; Fedorovich and Conzemius 2008). We leave the investigation of the wind shear dependency for future work. Additionally, we remark that a ZOM such as the one discussed correctly captures the qualitative behavior but tends to overestimate the observed entrainment flux. This is because, at long times, nonlinear and saturation effects start to become important in the entrainment dynamics and as a result, entrainment rate slows down. A linear model cannot account for these features. This indicates that a complete understanding of the entrainment dynamics is still an open issue, since more complex models do not substantially improve the estimation of the entrainment fluxes.
We have shown that entrainment influences the spatial organization of the scalar fields in the boundary layer, from the large down to the smallest resolved scales. If sufficiently strong, the effect of entrainment is to weaken the spatial correlation between temperature and trace gas concentration. In particular, the trace gas is no longer confined into thermal plumes, and tongues of clear air can travel from the top to the surface of the CBL, as it appears from the observed negative skewness profile. Moreover, the entrained clear air that gets to the surface is capable of reducing the number of cliffs in the trace gas distribution. Since the trace gas is no longer transported inside the plumes, the spatial location of strong temperature fluctuations (exceeding the variance) is different from the location of the similarly strong trace gas fluctuations. Such effects are reduced for a smaller entrainment flux.
Independently of the entrainment case study, the temperature fluctuations appear to be much more intermittent than the trace gas ones; this seems to be a quite general result, since it is observed also in laboratory experiments of Rayleigh–Bénard convection (Zhou and Xia 2008). This suggests a universal behavior of scalar fields in the mixing regions of convection flows, which is worthy of deeper investigation.
We note that in many convective boundary layer studies (see, e.g., Sorbjan 2006), humidity is described as a passive scalar field. A direct consequence is that our findings can be applied to the entrainment-drying or entrainment-moistening cases as well. In particular, humidity variance and spatial distribution in the lower part of the CBL is regulated by entrainment—evidence that should be taken into account when applying similarity laws. Finally, we have evidenced that in CBLs the probability to observe temperature or trace gas increments much larger than the standard deviation (i.e., cliffs) is maximal in the bottom part of the convective boundary layer, hence the importance of a correct description of turbulent fluctuations in the vicinity of the surface.
We thank Anna Maria Sempreviva for valuable comments and suggestions. We also thank Mariana Cassol and Umberto Rizza for useful discussions. This work was developed with the support of CNR Grant “Ricerca Spontanea a Tema Libero 2007.” Numerical simulations were performed at the Supercomputing Center CINECA (Italy), under the project “Large Eddy Simulation for Planetary Boundary Layer.” This work was partially funded by the project “SSD PESCA.”
Numerical Details and Formulation of the Dynamic Subgrid-Scale Model (DSM)
We recall some details of the numerical code first. In the integration of Eq. (4), pseudospectral integration is adopted in the homogeneous directions x and y, with the periodic boundary conditions. Finite-centered differences are adopted along the nonhomogeneous vertical direction z. Dealiasing is performed on horizontal planes applying the ⅔ rule to the nonlinear terms in the equations of motion and to the SGS model terms. Time integration is carried on using a third-order Runge–Kutta algorithm. A two-dimensional sharp spectral cutoff kernel is applied for both the grid and the test filters in the homogeneous kx and ky directions. There is no explicit filtering in the kz direction: the use of central finite differencing along z direction is equivalent to a top-hat filter applied in physical space. The grid filter width is , where, taking into account the dealiasing procedure, , i = x, y, and . The test filter width is , i = x, y. No explicit test filtering is applied along z.
The traditional Smagorinsky SGS model uses standard empirical values of the so-called Smagorinsky coefficient CS. In the DSM, two filter operators are defined: the grid filter, whose length is , and a test filter, with different width . The latter is generally taken to be twice the former, which corresponds to α = 2. The same functional form is assumed here for the two operators. Test filtering is denoted by a tilde; that is,
We choose to make the test filtering, as well as the grid filtering, in Fourier space with a cutoff filter. This has the advantage that the application to any function of the grid filter followed by the test filter has the same result as the direct convolution with the test filter. Therefore and the effective double filter width is .
We emphasize that the Dynamic Model was separately applied to the momentum, temperature, and concentration equations; this implies that three independent model coefficients were computed. We start discussing the application of the DSM to the momentum equation. By filtering the Navier–Stokes equations, with the grid filter first, and then the test filter, the new subgrid-scale stress tensor results in
The tensor Lij is entirely known in terms of the resolved velocity , and it can be viewed as the contribution to the SGS tensor (A2) from the largest unresolved scales. Within the Smagorinsky parameterization, one obtains
Note that the model coefficient C(x, t) is a real function of space and time and, as it will be clear from the sequel, it is not strictly positive. The substitution of terms (A4) and (A5) into the distortional part of Eq. (A3) gives
The final step is to obtain model coefficient C(x, t): this is done by minimizing the mean-square error between and its approximation as given by the Smagorinsky model (A6). This yields to the following expression for the model coefficient (Lilly 1992):
A direct calculation of the model coefficient with the previous formula would give numerical instabilities because of the strong fluctuations of both the numerator and the denominator. Averaging procedures, both in space and time, are thus generally applied (Germano et al. 1991; Ghosal et al. 1995; Meneveau et al. 1996). Since we are considering a homogeneous surface at the bottom of the CBL, the coefficient C(x, t) is computed by space averaging over the homogeneous horizontal directions (here the angle brackets indicate horizontal mean):
As a result, C depends on the z coordinate and on time t only. In addition, negative values of C(z, t) are set to zero. This procedure, known as clipping (Geurts 2003), is seldom used to limit instabilities. Backscattering of energy from the unresolved to the resolved scales of motion is ruled out only in an average sense. Local negative values of the numerator LijMij do contribute to the reduction of the plane average (Piomelli 1993).
In the past, the dynamic model has been applied under different stability conditions. The convective boundary layer was shown to be not very sensitive to the SGS parameterization (Nieuwstadt et al. 1991). Within it, the surface layer and the viscous sublayer near the bottom of the domain are shallow, and the application of the base (scale invariant) dynamic model is justified (see, e.g., H.-Y. Huang et al. 2008). On the other hand, a scale-dependent version of the model has been formulated, where the coefficients C appearing in Eqs. (A4) and (A5) are not the same but vary with the filter sizes. Such reformulation has proved to be very successful in neutral and stable boundary layers (see, e.g., Porté-Agel et al. 2000; Bou-Zeid et al. 2005; Stoll and Porté-Agel 2008).
For the scalar equations, the same procedure is repeated, giving the expression for the dynamic coefficient for the scalar SGS tensors closure. In the case of potential temperature, the effect of the unresolved scales on the filtered variable is described by the SGS temperature flux at the grid scale:
A classical parameterization gives
with the eddy diffusivity equal to . In many LES applications, the eddy diffusivity model coefficient is chosen to be Cθ = C/Prsgs, which is proportional via the subgrid-scale Prandtl number Prsgs to the model coefficient for the momentum (A8). In the present application, however, it is determined independently. In particular, the temperature flux is evaluated at the test scale also:
Repeating the analogous of the Germano identity (A3) for scalars, it is obtained
The eddy diffusivity model coefficient is computed as
(here the angle brackets indicate horizontal mean). A final point to be discussed is the behavior of the eddy viscosity and diffusivities in the upper part of the atmospheric boundary layer: there, stratification effects are important, particularly when going through the entrainment zone. From the point of view of numerical modeling, it is desirable to have eddy viscosity and diffusivities equal to zero in the stably stratified free atmosphere. Since there is no mechanical or thermal production of turbulence, we do not expect any turbulent dissipation there.
Stratification effects can be included by multiplying the eddy viscosity for a correction factor CB (Lilly 1962):
The Richardson number,
is evaluated runtime during the simulation, while the Prandtl number is kept constant Pr = 0.4.
Note that an instantaneous Prandtl number can be computed through the ratio of the dynamic eddy viscosity to the dynamic eddy diffusivity: it would therefore be a function of the coordinate z and time t. Here, local fluctuations of Pr were disregarded by choosing its value consistently with previous boundary layer studies (e.g., Moeng 1984). We mention that a posteriori analysis of our results indicates that, for instance in run S256, the time-averaged Prandtl number is Prsgs ≃ 0.33 in the mixed layer, while it increases up to 0.95 in the entrainment region.
Taking into account the buoyancy correction, the eddy viscosity for momentum equation used in the present CBL simulation is
with C evaluated through Eq. (A9). The eddy diffusivity for the heat equation is
with Cθ evaluated through Eq. (A14). Accordingly, the eddy diffusivity for the concentration equation is
with Cc evaluated by the expression obtained from Eq. (A14) after replacement of θ with c.
Simulation results for the eddy coefficients (again referring to case S256) are plotted in Figs. A1 and A2. The dynamic procedure to compute the closure coefficients for the velocity and the scalar equations is repeated once every five time steps.