The role of turbulent mixing in the formation of the structure of stratocumulus clouds is investigated using a Lagrangian–Eulerian parcel cloud model. The model contains approximately 2000–5000 adjacent parcels with the linear size of 25–40 m, moving with a turbulent-like velocity field with observed energetic and statistical properties. The process of turbulent mixing of Lagrangian parcels is parameterized using the k-epsilon theory extended to the case of mixing nonconservative values. The model includes the interaction of cloud and the overlying inversion layer. The stratocumulus clouds observed during flight RF01 of the Second Dynamics and Chemistry of the Marine Stratocumulus field study (DYCOMS II) are simulated.
Effects of turbulent mixing are analyzed by comparing simulations with and without mixing. When mixing between parcels is included, the thermodynamical and microphysical structure of the measured stratocumulus clouds is properly reproduced. Mixing leads to a more uniform cloud structure with well-defined borders. Good agreement is found between Paluch diagrams calculated in the model and those reproduced from measurements. The radius of correlation of liquid water content and other variables calculated in the model is on the order of several hundred meters and agrees well with observations. When mixing is not included, the radius of correlation is on the scale of a single parcel and the cloud layer contains dry entrained parcels, making the microphysical structure unrealistic. It is also shown that turbulent mixing leads to an increase in the effective radius and facilitates and accelerates drizzle formation. The time in which a 40-m air parcel preserves its identification is estimated from the results and is found to be on the order of 25 min.
The effect of turbulent mixing between cloudy volumes as well as between cloudy volumes and droplet-free environment air on the structure and microphysical cloud properties has long been a matter of investigation. Within this field much attention has been directed to the formation of raindrops in convective clouds and drizzle in stratocumulus clouds. Formation of raindrops and drizzle is usually related to the existence of large droplets with radii of about 22 μm. Several mechanisms have been proposed to explain the formation of such large droplets—for example, enhanced collisions between cloud droplets (Pinsky and Khain 2002), fluctuations of turbulent velocity accompanied by in-cloud nucleation (Erlick et al. 2005), supersaturation fluctuations in turbulent flows (Khvorostyanov and Curry 1999), and turbulent mixing between drops in updrafts and downdrafts (Korolev et al. 2013). Baker and Latham (1979, 1982), Jensen and Baker (1989), Baker et al. (1980), and Blyth et al. (1980) pioneered studies in which the formation of large droplets was related to the process of entrainment and mixing of cloudy and cloud-free air. They introduced a concept of homogeneous and inhomogeneous mixing, which is typically referred to in studies dedicated to the effect of mixing on the droplet size distributions (DSD) (Devenish et al. 2012). The type of mixing depends on the relationship between the characteristic scales of mixing and the phase relaxation time (Lehmann et al. 2009; Kumar et al. 2012; Andrejczuk et al. 2009). In case of extremely homogeneous mixing, homogenization due to turbulent mixing is much faster than adaptation of supersaturation, so all droplets experience the same supersaturation. Simple estimations show that extremely homogeneous mixing takes place at spatial scales lower than 0.5 m. At larger scales mixing is always inhomogeneous. In the case of extremely inhomogeneous mixing, supersaturation in different parts of the mixing volumes is different and some fraction of the droplets evaporates in dry air.
There are contradictory opinions concerning the role played by extremely inhomogeneous mixing in the formation of large droplets and first raindrops in cumulus clouds. Numerical results of Cooper et al. (2013) seem to lend some support for the idea that inhomogeneous mixing leads to an increase in drop size. In contrast, a numerical study by Schlüter (2006) has shown that extremely inhomogeneous mixing does not result in formation of larger droplets, but leads to spectrum broadening toward smaller drop sizes. These results were reached even though quite different methods were used; Schlüter (2006) used the method known as “explicit mixing parcel model” (Krueger et al. 1997), while Kumar et al. (2012) used a three-dimensional, direct numerical simulations (DNS) method in an Eulerian–Lagrangian framework.
Hill et al. (2009) simulated the evolution of a stratocumulus cloud using a large-eddy simulation (LES) model, in which mixing was described in two ways: either through the traditional use of 1.5-order equation closure or with additional corrections of droplet concentration required by the inhomogeneous mixing hypothesis. The results in both cases were indistinguishable from one another. The insensitivity of evolution of stratocumulus clouds to the subgrid mixing assumption was explained by the difference in the rates of condensation and evaporation caused by resolved dynamics, in contrast to the rates caused by subgrid processes, which are smaller by two orders of magnitude.
It is difficult to distinguish between the types of mixing in clouds merely through observations. The uncertainty is well formulated in the study of Gerber et al. (2008), who examined the microphysical structure of small cumulus clouds observed in the Rain in Cumulus over the Ocean (RICO) field program and concluded the analysis by stating that “Entrainment causes primary dilution of the drops, without significant size changes, thus either extreme inhomogeneous mixing or more likely homogeneous mixing resulting from mixing with cool and humid entrained air takes place” (p. 88). Their in situ measurements are supported by the radar observations of Knight and Miller (1998), as well as by Bar-Or et al. (2012), who reported the existence of a 100–150-m layer of high relative humidity (RH) around cumulus clouds. At high RH the effects of homogeneous and inhomogeneous mixing become indistinguishable (Burnet and Brenguier 2007).
Observational studies and numerical simulations of deep convective clouds (Paluch 1986; Prabha et al. 2011; Benmoshe et al. 2012; Korolev et al. 2013) have shown that undiluted (i.e., adiabatic) and slightly diluted air volumes have maximum liquid water content (LWC) and contain the highest amount of the largest drops. Furthermore, the formation of first raindrops takes place in the undiluted and weakly diluted parcels.
It is worth noting that the effect of turbulent mixing on cloud microphysics goes far beyond the problem of whether the mixing is homogeneous or inhomogeneous. Several alternative mechanisms have been proposed by means of which mixing can accelerate the formation of large droplets (e.g., Korolev and Isaac 2000; Korolev et al. 2013). As will be shown below, the role of mixing involves much more than just the possible acceleration of large droplet formation: mixing is necessary in order to realistically describe microphysical and statistical cloud properties.
Significant efforts have been made in recent years to improve the representation of microphysical processes in numerical models, especially by decreasing numerically induced DSD broadening, which accelerates formation of drizzle and raindrops (Khain et al. 2000). A variety of methods have been proposed to decrease this artificial DSD broadening, the most effective of which is to avoid the use of a regular mass grid and to allow the droplets to be of all possible sizes, in accordance with the equation of diffusion growth. There are two types of models that allow such an approach. Although they are both either “Lagrangian” or “hybrid Eulerian–Lagrangian” models, they are based on quite different approaches.
The first approach simulates the motion and growth of individual droplets (Lagrangian droplets) in a given velocity field. These models result from the development of DNS (Devenish et al. 2012) to include calculation of supersaturation and droplet growth. For instance, this type of model was used by Kumar et al. (2012) to simulate DSD evolution in a turbulent flow.
To simulate droplet motion and growth on real cloud scales, some models simulate droplet movement within a flow field generated in LES. Such models, sometimes referred to as Lagrangian cloud models (LCM), have also been described in studies by Andrejczuk et al. (2009, 2008, 2010), Shima et al. (2009), and Riechelmann et al. (2012). Use of these models requires handling an extremely large number of droplets. To solve this problem, the concept of superdroplets was introduced (Shima et al. 2009). Each superdroplet represents a large number of real droplets of equal size and location. The process of collisions is reduced to collisions between superdroplets. These potentially powerful models are currently at the development stage and further efforts are required to properly include processes related to collisions, droplet nucleation, and the formation of rain drops. Mixing in such models is treated explicitly, without using any parameterization method, and the accuracy of the representation of this process depends on the ability of the corresponding LES model to adequately reproduce the turbulent structure, and on the LES grid spacing. The grid spacing, often set to a few tens of meters, introduces the minimal scale of turbulent vortices involved in the mixing process.
The second approach of Lagrangian models is widely used, either in single-parcel models (Bower and Choularton 1993; Pinsky and Khain 2002) or in trajectory ensemble models (TEM) (Stevens et al. 1996; Feingold et al. 1998; Harrington et al. 2000; Cooper et al. 2011; Erlick et al. 2005). In contrast to LCM, the basic element in these models is a Lagrangian air parcel, often treated as a moving adiabatic cloud volume. The advantage of these models, from a microphysical point of view, is that they avoid utilization of simplified procedures for droplet nucleation. The use of movable mass grids, in which the mass of bins changes according to the equation for diffusion growth, eliminates artificial DSD broadening. Usually, TEMs comprise two submodels: the LES model and the Lagrangian parcel model. The DSDs are calculated in several hundreds of Lagrangian cloud parcels moving within the velocity field produced by an LES model. This approach allows the simulation of DSD formation in parcels with different histories and to investigate the reasons for DSD variability. Note that in TEM cloud parcels are usually separated by significant distances and sedimentation is not taken into account. In addition, the averaging of separate cloud parcels located at the same height is not equal to the averaging over the horizontal levels within the cloud layer, which makes the comparison of TEM results with observations difficult.
Khain et al. (2008) and Pinsky et al. (2008) have developed a new generation TEM, in which the Lagrangian parcels are adjacent to one another and cover the entire computational area. The time-dependent turbulent-like velocity field is generated by a model that reproduces observed turbulent intensity and correlation properties. In these studies diffusion growth and collisions of droplets, as well as drop sedimentation from one parcel to another, were calculated on mass grids containing 500 mass bins. Since drop sedimentation is calculated within the Eulerian framework, the model is referred to as a hybrid Lagrangian–Eulerian model of stratocumulus clouds.
Using the model, Magaritz et al. (2009) first simulated the formation of drizzle drops in stratocumulus clouds and their fall to the surface within a “nonmixing” limit; that is, turbulent mixing between the parcels was not included. The model simulated observed microphysical parameters and drizzle size quite well. Moreover, the correlation properties of LWC and droplet concentration were in good agreement with observations. These results could create the false impression that mixing is insignificant in the formation of a dynamical and microphysical cloud structure. It is to be noted that in the model used by Magaritz et al. (2009) no inversion layer and no entrainment of dry air from above were allowed. The lack of entrainment from the upper boundary artificially decreased the role of turbulent mixing between clouds (see discussion below).
Here we present the new model version, in which several processes have been added, with turbulent mixing between the Lagrangian parcels being the most important one. This paper focuses on the effects of turbulent mixing and entrainment of dry parcels on macroscopic properties of stratocumulus clouds, such as cloud geometry, vertical profiles of average characteristics, and correlation properties of major microphysical quantities. This model will also be used in future work, where the effect of the mixing process on the DSD and drizzle formation will be analyzed, in particular the presence of “lucky parcels,” in which the first drizzle drops form.
2. New model configuration and experimental design
For this study a spectral bin microphysics hybrid Lagrangian–Eulerian model of a stratocumulus topped boundary layer (STBL) is used. The initial version of the model is described in detail in Pinsky et al. (2008) and Magaritz et al. (2009). The model is constructed of about 2000 adjacent parcels with a characteristic (average) linear size of 40 m. Supplemental simulations contained over 5000 parcels with a linear size of 25 m. The parcels cover the entire 2D computational area of 2500 × 1250 m2 and describe all parts of a STBL, from the ocean surface where surface fluxes are calculated, to the top of an approximately 300 m deep warm and dry inversion layer. The parcels are advected by a turbulent-like velocity field.
The velocity field is represented as the sum of a large number of harmonics with random time-dependent amplitudes. It is statistically uniform in the horizontal direction and obeys the Kolmogorov law, in agreement with available turbulent measurements in maritime STBL. The time-dependent turbulent-like velocity field is assumed quasi stationary during the entire simulation, with observed energetic and statistical properties. To reproduce the observed statistical and energetic properties, the amplitudes of the harmonics have been calculated using two measured quantities: the vertical profile of rms of velocity fluctuations, (where w′ are the fluctuations of vertical wind velocity and brackets indicate horizontal averaging) and the longitude structure function (Magaritz et al. 2009; Pinsky et al. 2008). Below the base of the inversion layer the dissipation rate is set at 10 cm2 s−3. This value is typical of stratocumulus clouds under consideration (Siebert et al. 2006). A sharp decrease of the dissipation rate to zero was assumed within the inversion layer above cloud top (Katzwinkel et al. 2011).
All parcels at t = 0 min are subsaturated and contain only wetted aerosols (haze particles). As the parcels move in the BL some may cross the lifting condensation level, and droplets will form within them; these parcels now contain cloud droplets as well as wetted aerosols. The droplets may continue to grow or evaporate and release the aerosol back into the parcel. In each parcel all microphysical processes are calculated. Nucleation, diffusion growth of haze particles and droplets, evaporation, and denucleation are described using the diffusion growth equation (Pinsky et al. 2008; Magaritz et al. 2010). To accurately calculate diffusion growth of wetted aerosols a small time step of 0.01 s is used. In each parcel the process of growth by collisions is described using the stochastic equation for collisions and high-resolution tables for collision efficiencies presented by Pinsky et al. (2001). Both aerosol and drop distributions are calculated using a single 500-bin mass grid with a 0.01–1000-μm-radius range. The single grid of high resolution allows explicit separation between haze particles and cloud drops without parameterization of droplet nucleation. During all of these processes the aerosol and drop distributions in each individual parcel are calculated separately and change during the simulation.
The model has been substantially modified as compared to previous studies. Most prominently, the new model version contains two types of interactions between the parcels: sedimentation and turbulent mixing. Sedimentation allows larger droplets that form in the cloud parcels to act as drop collectors during their fall and reach the surface as drizzle. Since the effects of turbulent mixing are the focus of this study, we describe the algorithm in more detail.
In the Eulerian LES models based on the system of averaged equations, the effect of turbulent mixing is represented by terms such as resulting from Reynolds averaging (a is a generic scalar value). Usually these terms are parameterized by means of the semi-empiric K theory, according to which turbulent fluxes are expressed as gradients of mean values: . The turbulent coefficient K is calculated using the 1.5-order closure approach, based on solving the equation for turbulent kinetic energy (Benmoshe et al. 2012). We apply the mixing algorithm between adjacent Lagrangian parcels proposed in the study by Pinsky et al. (2010). This algorithm represents an expansion of K theory for the case of mixing of both conservative and nonconservative values (such as DSD) given on a nonregular spatial grid formed by parcel centers.
Let us consider two air volumes (parcels) located at different heights, where their centers are separated by the distance l. Fig. 1 schematically illustrates the reciprocal penetration of turbulence-induced filaments belonging to parcel 1 (P1) into parcel 2 (P2). Turbulent fluxes between parcels are calculated based on the analogy with K theory, where the quantities in the parcel centers are regarded as the mean values while the quantities in filaments are considered as turbulent fluctuations in vortices moving from one level to another.
Let a be a nonconservative value to be mixed. Denote initial values of a in the centers of P1 and P2 as and , where index denotes the number of the parcel and arguments mark the level at which these values are calculated. As mass ascends–descends in the filaments, the value of a changes because of the nonconservative properties of the parameter. During an ascent of mass in a filament from the center of P1 to the level of P2, the value changed to . Similarly, during the descent of mass in a filament belonging to P2, the value changes to . Changes over time of quantity a due to turbulent fluxes between two volumes are represented as (Pinsky et al. 2010)
where is the turbulent exchange coefficient. The terms and in Eqs. (1) represent the changes in the value a along the mixing length l and as such the non-conservativeness of the quantity a. When conservative values are mixed, these terms are equal to zero. We assume that mixing takes place in isotropic turbulence within the inertial subrange, where (Richardson’s law), is the turbulent kinetic energy dissipation rate, and C = 0.2 (Monin and Yglom 1975). The distance between the parcel centers l is regarded as the mixing length. In an Eulerian framework, the right-hand sides of Eqs. (1) are reduced to the standard mixing operator . Solutions of Eqs. (1) take the form (Pinsky et al. 2010)
where and are the initial values (before mixing), the terms and express the effect of non conservativeness of the mixing quantities, and
is the characteristic mixing time corresponding to the mixing length l.
Temperature and humidity change as a result of the mixing of conservative values: total water mixing ratio and liquid potential temperature , where , , and are water vapor mixing ratio, liquid water mixing ratio, and potential temperature, respectively; also is the latent heat of evaporation and is the specific heat capacity.
The calculation of changes of DSD during the mixing is reduced to an estimate of the drop radii changes caused by supersaturation fluctuations, which in turn are caused by turbulent velocity fluctuations. The calculations are performed using approximate solutions of the diffusion growth–evaporation equation and the equation for supersaturation, under the assumption that turbulent velocity fluctuations obey turbulent laws in the inertial subrange with a given dissipation rate.
The mixing algorithm has the following steps: 1) mixing of and using Eqs. (2) taken in their conservative form (i.e., and ), 2) calculation of changes in the drop radii and water liquid mixing ratio in filaments during their displacement by the mixing length l for each parcel pair, 3) mixing of DSD between parcels and filaments using Eqs. (2) and according to the scheme in Fig. 1 and 4) determining after mixing in each parcel by integrating droplet mass over all the bins, and 5) calculation of the water vapor mixing ratio and the potential temperature in each parcel after mixing as and .
The maximum distance between the centers of the parcels, at which mixing is still performed, is chosen as 50 and 30 m for cases with the characteristic parcel size of 40 and 25 m, respectively.
To investigate the possible effects of the mixing procedure, two different scenarios have been compared. The first scenario takes into account the change in drop size and corresponding latent heat release during mixing, as described above. In the second scenario no change in droplet size is assumed during the mixing procedure (i.e., DSDs were mixed as conservative variables). In the latter case there is no latent heat release during mixing and it is referred to as passive mixing. Note that in the mixing algorithm used in the study superstaturation is the same for all droplets in a given parcel so the mixing obeys the properties of homogeneous mixing. In the mixing scenario in which DSDs are considered as nonconservative values, the DSDs in filaments and in parcels obey some properties of inhomogeneous mixing. As mentioned in the introduction, the strictly homogeneous mixing takes place at scales lower than 0.5 m. However, mixing between parcels within the cloud, where relative humidity is high, makes the mixing types indistinguishable one from the other.
Other modifications of the model include implementation of sensible and latent surface fluxes as well as radiative longwave cooling. Background wind at 10 m was set equal to 10 m s−1. The sensible and latent heat surface fluxes were calculated using bulk-aerodynamic formulas with a Dalton number of CE = 0.002 (Smith 1988). It is assumed that the model computational area is perpendicular to the background wind, so no background wind component appears. The effect of background wind speed manifests itself only through the effect on the surface fluxes. Implementation of the surface fluxes in the model allows us to begin the simulations when all Lagrangian parcels are undersaturated and to simulate formation of a cloud with continuously increasing water vapor mixing ratio within the computational area.
Longwave radiative cooling was included using the two-stream approximation following Khvorostyanov (1995) and Khvorostyanov et al. (2003). This parameterization of radiation was successfully used for calculation of radiative cooling in stratocumulus clouds of different types.
In this study, stratocumulus clouds were simulated under conditions close to those observed during the research flight RF01 of the Second Dynamics and Chemistry of the Marine Stratocumulus field study (DYCOMS II) (Stevens et al. 2003). During this night flight the stratocumulus cloud measured was ~300 m thick, capped by a strong inversion at 850 m. No drizzle was measured at the surface.
To generate the turbulent-like velocity field the vertical profile of was taken from measurements (Stevens et al. 2005) and implemented in the model. The quantity was assumed zero above z = 800 m, as shown in Fig. 2, to account for the inversion layer. The magnitude of at the boundary between the cloud layer and inversion determines the intensity of droplet-free parcels penetrating the cloud layer and cloudy parcels penetrating the inversion layer (i.e., affects the intensity of entrainment). In contrast, small-scale turbulent mixing between parcels is determined by Eqs. (1)–(2).
The initial aerosol distribution was derived from observations (total concentration 200 cm−3, radius range 0.01–1.3 μm) and for most simulations was assumed to be the same for all parcels at t = 0 (Magaritz et al. 2009). Initial humidity and temperature were selected to allow the simulation to begin under subsaturated conditions. Cloud development in the model initiates as some of the parcels become saturated. During the simulation the BL becomes well mixed, because of the parcel motion, and surface fluxes increase humidity to the observed values.
The initial temperature and humidity values are assumed to be horizontally uniform at t = 0. Relative humidity is set to approximately 90% below the inversion level and rapidly decreases with height within the inversion layer.
Several simulations were performed in order to reveal the role of turbulent mixing between Lagrangian parcels. In three simulations the linear parcel size was set at 40 m: mixing (MI), no mixing (NoMI), and passive mixing (MP), in which DSDs were mixed like conservative values and no latent heat was released during the mixing. In contrast to these simulations, an additional simulation (MIA) was performed, in which the aerosol concentration was only constant below the inversion layer; above it the concentration decreased linearly to zero at the top of the computational area of 1250 m. This simulation was performed in order to test the hypothesis that the constant profile can lead to penetration of new aerosols from the inversion layer and their activation may increase droplet concentration. Finally, two simulations were performed to examine the effect of parcel size on the structure of simulated clouds (MI25 and NoMI25). In these simulations the linear size of the parcels was set to be equal to 25 m. The list of simulations and their conditions are presented in Table 1. The model simulation starts at t = 0 with a clear BL and the simulation of clouds continued for a period of approximately 8 h. Note that surface fluxes are included during the entire simulation, leading to a graduate increase in the STBL humidity and an increase in the depth of the cloud layer. For the analysis we use a time period when all of these changes are very gradual and the process of drizzle formation has not started.
Measurements from flight RF01 of the DYCOMS II field experiment were used for validation of the model.
a. Effects of mixing on the cloud geometric structure
All simulations initially start in a subsaturated state; because of parcel motion, cloudy parcels form both in the MI and NoMI simulations. However, the cloud structure in these two cases is quite different. After about 150 min, a uniform cloud is established in the MI case, and a quasi-stationary state remains until the end of the simulation. The first impression of the extent of the effect of turbulent mixing on the cloud layer structure can be obtained by analyzing the LWC and droplet concentration fields in the MI and NoMI simulations presented in Fig. 3. It is evident that turbulent mixing has a profound effect on the geometrical and microphysical structure of the cloud layer. In the MI simulation a uniform cloud deck forms with a coherent LWC maximum near cloud top and a clear sharp decrease of LWC at cloud boundaries. The cloud base and cloud top are well defined. In the NoMI case the cloud layer contains many droplet-free air volumes and the cloud structure is highly nonuniform. In this case the cloud boundaries are not easy to recognize. In the MI case the droplet concentration in the cloud layer is nearly constant horizontally and vertically except at the cloud boundaries. This agrees with measurements taken during flight RF01 and as a general feature of stratocumulus clouds (Pawlowska et al. 2000). This feature is not seen in the NoMI simulation where many parcels with very high or low concentrations can be found. Also presented in the figure are the profiles of standard deviation for these fields. As expected, the standard deviation profiles indicate that the variability in LWC in both cases reaches its maximum near cloud top in the interface area between very different parcels from the cloud and the inversion layers. Below this maximum the standard deviation is substantially lower in the MI case. Figure 3 is a clear indication of the dramatic role of mixing in case of entrainment of inversion air into the cloud. As will be shown below, the areas of reduced LWC and droplet concentration are largely formed by dry and warm air parcels entrained into the cloud through the cloud top; in the MI simulation these parcels mix with cloudy parcels and become an integral part of the cloud while in the NoMI simulation they remain as “holes” in the cloud layer.
b. The effect of mixing on the mean vertical profiles
Figure 4 compares the vertical profiles of horizontally averaged main microphysical and thermodynamical parameters in some of the simulations; all the profiles are averaged during the period of 170–310 min of the simulations. There is good agreement between the calculated profiles and the measurements of flight RF01 (also presented in Fig. 4). The LWC increases with height nearly linearly until a maximum of 0.6 g m−3 at the upper cloud boundary. The horizontally averaged concentration is nearly constant throughout the cloud. The vertical profiles of the temperature and the liquid water potential temperature are in good agreement with the measurements as well. The total water mixing ratio is constant throughout the boundary layer and is in agreement with the measured profile.
Even though the fields of LWC and droplet concentration in the MI and NoMI presented in Fig. 3 appear quite different, the average profiles calculated from these simulations do not show a significant difference. A strong effect of mixing on the averaged profiles is seen in the inversion layer, where the penetration of cloudy volumes into this dry layer and subsequent mixing lead to droplet evaporation and to an increase in air humidity.
Figure 4 shows that mixing leads to a small decrease in LWC and humidity and to an increase in the temperature within the cloud layer. These differences can be attributed to the effects of heating and drying caused by entrainment and mixing with warm and dry parcels that penetrated from the inversion layer. Small differences within the cloud layer can be seen in the LWC profiles in the MI and NoMI runs, possibly indicating that mixing within the cloud layer is largely of “mechanical” nature, and that the turbulent mixing itself is not accompanied by substantial changes in temperature or humidity, which is in agreement with observations by Gerber et al. (2008). In the MI simulation values in each parcel, such as humidity or LWC, represent mixtures of cloudy and penetrated air. The mixture has the same proportion as in the case of simple horizontal averaging of noninteracting cloudy and drop-free volumes (as in NoMI). Thus, horizontal averaging leads to nearly similar results, even if particular values in the NoMI simulations are unrealistic. To justify this explanation, a supplemental simulation was performed, in which the latent heat release during the mixing process was excluded (MP). All vertical profiles were very similar to the ones plotted in Fig. 4 (not presented). This result indicates that the main heating and cooling effects are controlled by the rates of condensation/evaporation caused by resolved dynamics during updrafts and downdrafts of the parcels. Similar conclusions were reached by Hill et al. (2009).
The LWC profile in both simulations is close to adiabatic (Fig. 4). In the NoMI simulation each parcel is adiabatic so the averaged profile should also be close to adiabatic. In the MI simulation the dry parcels entrained from the inversion layer are mixed with cloudy parcels within the cloud layer and decrease the averaged LWC. Most cloudy parcels, even with different rates of dilution, remain supersaturated or saturated. Moreover, cloudy parcels contain droplets that increase humidity to nearly saturated values in cloud downdrafts. Accordingly, the LWC profile in the MI simulation is also close to adiabatic, but with a lower value of LWC, which corresponds to a slightly higher location of the lifting condensation level (cloud base).
Mixing increases the droplet concentration in the simulated clouds (Fig. 4). Results from the MIA case, in which the CCN concentration decreased in the inversion layer, indicate that this increase is a consequence of activation of aerosols entering the cloud layer through entrainment at cloud top. The droplet concentration obtained in the MIA run is lower than in the MI run and substantially closer to observations. Despite the changes in the concentration, the LWC profiles in the MI and MIA cases are very similar (not shown). As in the comparison between the MI and NoMI, one might suggest that LWC, unlike droplet concentration, is mostly determined by parcels ascending from cloud base and not by parcels entering through cloud top.
c. Mixing diagrams (Paluch diagram)
For analysis of turbulent mixing, conservative variable mixing diagrams (Paluch diagram) are often used (Paluch 1979; Burnet and Brenguier 2007). These diagrams express a relationship between conservative variables, which do not change in the adiabatic process, such as the total water content (the sum of water vapor content and LWC) and liquid water potential temperature . Since mixing is a nonadiabatic process, a parcel that mixes with its surroundings will be seen on the diagram as having a specific “trajectory” determined by the properties of the two air masses that mix. Larger shifts of the parcel on the Paluch diagram indicate stronger mixing or a greater difference between the two air volumes. When turbulent mixing is not included in the model, all parcels are nearly adiabatic. Deviations from adiabatic are due to droplet sedimentation and surface fluxes. On a Paluch diagram, a single point ( versus ) in the MI case representing one air volume will not change its location during cloud evolution.
Parcel paths on the Paluch diagram in the MI and NoMI cases are illustrated in Fig. 5. A single parcel that initially starts from the inversion layer and then enters the cloud is selected. The parcel has the same trajectory in both MI and NoMI simulations (left panel). On the Paluch diagram the path of the parcel for both cases is very different (right panel). Each point on the diagram denotes the parcel location every 5 min from 170 to 445 min of simulation; the color (in the MI case) indicates the height of the parcel in the model domain. At t = 0 min the parcel is dry and warm and located at the bottom of the Paluch diagram. In the MI case, as the parcel approaches the cloud layer it becomes colder and moister and its location on the diagram shifts. This parcel eventually becomes part of the cloud layer reaching high values of . In the NoMI case, the parcel locations at the same time steps are presented. However, all the points appear at the same location on the Paluch diagram since the parcel remains warm and dry throughout the simulation.
In Fig. 6 the Paluch diagrams for the MI simulation (Fig. 6a) and for measurements during RF01 (Fig. 6b) are presented. Each dot represents a single air parcel in a 300-m-thick layer at the top of the cloud and bottom of the inversion layer during 100 min of simulation. Dark round markers represent cloudy parcels, in which the droplet concentration is higher than 3 cm−3 and gray triangles denote noncloudy parcels. There are two areas that can be identified in the diagram. The first is the cloud layer with cool and moist parcels; the second is the inversion layer with warmer and dryer parcels. In between these two sections a mixing area is seen connecting the extreme values. The mixing area consists of parcels that either ascend from the cloud into the inversion layer, where they become dryer and warmer, or parcels that are entrained into the cloud where they become moister and colder. By the end of the mixing process a parcel that either penetrates the cloud or the inversion layer becomes inseparable from its environment. The slope of the mixing line is determined by the parameters of the cloud and the inversion layer. One can see a reasonable agreement of the Paluch diagram obtained in the MI simulation with that plotted using the observed data (Fig. 6b).
Figure 7 shows the mixing diagram for the NoMI simulation and again demonstrates the unrealistic structure of the cloud in this case. This figure is dramatically different from that presented in Fig. 6. In the NoMI case a lot of scatter can be seen in the entire diagram and the clear mixing line disappears. As stated before, in the case of adiabatic processes an air parcel will keep its conservative quantities and essentially remain at the same point on the diagram. In this simulation parcels are nearly adiabatic. Changes to the conservative values can arise from two processes that are also represented in the model: surface fluxes and radiative cooling. Parcels with high qt in the diagram obtained humidity because of surface moisture flux when they were located near the surface. In the mixing simulation these parcels are responsible for the moistening of the entire BL. On the diagram, parcels with low values of θl demonstrate the effect of radiative cooling. These parcels are located at the upper levels and contain significant LWC. A comparison of this figure with Fig. 6 is a clear indication that simulations without mixing cannot properly describe the thermodynamics of the cloud-topped boundary layer.
d. The lifetime of air volumes
When discussing the mixing between Lagrangian parcels one must inquire as to the characteristic lifetime of a single parcel, during which the parcel keeps its identification and can be traced among the surrounding air volumes. To evaluate the time frame of this process several parcels that ascend from the ocean surface or penetrate the cloud from above were chosen. Initially, temperature and humidity in these parcels are substantially different from those in the parcels in their immediate surroundings. During the mixing process, the parcel becomes similar to its environment and the difference reduces to approximately zero. The time period from when a parcel has a maximum difference from adjacent parcels until it blends with its environment can be considered as the lifetime of a particular air volume. Figure 8 presents six different examples of time dependencies of differences in between selected parcels and their environment. Each line represents a single parcel. Parcels with additional humidity from surface fluxes will initially have a positive difference, while parcels that penetrate the cloud through cloud top will initially have a negative difference. Using the conservative eliminates changes caused by parcel motion so only the mixing process is represented. Figure 8 shows that the characteristic time during which the air in the parcel preserves its identity is about 25–30 min. The characteristic time can also be defined as the time during which the initial difference in the parcel variables (say or ) decreases e-folding time. In this case the characteristic lifetime can be estimated as ~20 min. As could be expected, this time is close to the mixing time used in the mixing procedure and determined by Eq. (3).
e. Correlation functions and other statistical characteristics
Figure 9 presents the normalized horizontal correlation functions of LWC, droplet concentration, and temperature in the MI and NoMI simulations. The horizontal correlation function is defined as , where is the length of the series and and are LWC (or other parameter) at the locations and . The normalized horizontal correlation functions are defined as the ratio . All the correlation functions are calculated at the height of 700 m during the time period between 170 and 310 min. The correlation function of the vertical velocity, reproduced from observations, is also presented in the figure. For all parameters the correlation lengths in the MI case are considerably larger than the same functions in the NoMI case. While in the MI case the correlation length is on the scale of a few hundred meters, close to the correlation length of vertical velocity, in the NoMI run the correlation length is on the order of the linear size of the parcel. This result is expected since in the NoMI case there are no interactions between the different parcels, so entrainment of dry parcels from the inversion lead to a “broken” cloud structure. The difference in the correlation lengths reflects the difference in the cloud structures simulated with and without mixing as seen in Fig. 3.
In the MI case the correlation length of LWC is larger than the correlation length of droplet concentration. This difference arises from the nature of the measured parameter. LWC is an integral parameter determined by the total effect of supersaturation on droplet size. For a single parcel the LWC is determined mostly by the parcel’s history: the lifting condensation level and its residence time in the cloud. As such, LWC responses are relatively slow to fluctuations of supersaturation caused by turbulent velocity changes. Indeed, the correlation length of LWC turns out to be larger than that of vertical velocity (Fig. 9). This reasoning applies also to temperature. Temperature is also a slowly changing parameter that is mostly affected by the height of the parcel; it has a large correlation length. The droplet concentration is the most rapidly changing value and has a shorter correlation length (~300 m). Droplet concentration is largely determined by the concentration of the smallest droplets, which in turn is sensitive to comparatively small and high-frequency fluctuations of supersaturation. Small droplets may form as a result of nucleation in entrained parcels; the concentration of small droplets can change because of in-cloud nucleation. These processes affecting the concentration of the smallest droplets do not affect LWC and it will remain nearly constant.
Figure 10 compares the correlation functions of LWC, temperature, and droplet concentration derived from the model simulation with those derived from in situ measurements. The correlation functions from measurements were calculated using 25-Hz-resolution data. Three flight lags were chosen within the height range 620–750 m (the average height of all chosen lags was 670 m). In each lag stationary segments were chosen and the correlation functions were calculated. The correlation functions presented in the figure were calculated by averaging the correlation functions at each one of the flight lags. Figure 10 also shows maximal and minimal values of the correlation functions in these lags (marked by a grayscale). In agreement with the results of simulations, the correlation length for the temperature is the longest, and the correlation function of droplet concentration is the shortest. The agreement between the correlation lengths calculated in the model with those derived from observations indicates that the mixing is simulated realistically in the model.
To support the conclusion that the correlation length is determined by the properties of the microphysical quantities, supplemental simulations were performed with a characteristic linear parcel size of 25 m (simulations MI25 and NoMI25). The correlation functions for these simulations are presented in Fig. 11. One can see that correlation lengths in the NoMI25 case are again on the order of the parcel size, while correlation lengths in MI25 are comparable with those shown in Fig. 9 (i.e., a few hundred meters).
In this study we wish to examine to what extent different microphysical variables of stratocumulus clouds are affected by turbulent mixing, both within the cloud and at the boundaries. Figure 12 presents histograms of humidity and temperature in the MI and NoMI cases. In these histograms only parcels at cloud top and the inversion layer, between 700 and 950 m, are presented. The histograms quantitatively show the effects of mixing mentioned above in the analysis of the Paluch diagrams. In the histograms, we separate between entrained parcels that were initially in the inversion layer (white bars) and other cloudy parcels (dark bars). This separation allows us to better understand the transformation of cloud thermodynamical and microphysical structure due to entrainment and mixing. In the NoMI case both temperature and humidity histograms indicate the existence of two peaks clearly separating parcels located initially within the inversion layer and all other parcels. Extremely dry parcels in the NoMI case, which were initially located at upper levels, remain as gaps in the cloud (Fig. 3). Histograms in the MI case indicate a single narrower distribution. In this case the entrained parcels become part of the cloud and will eventually become indistinguishable from the rest of the BL parcels.
It is interesting to focus on the changes in the cloud layer itself (dark bars, Fig. 12). As parcels from the cloud top penetrate the cloud they will warm and dry the parcels in their environment and will lead to a shift in the histograms seen in the figure. For instance, in the NoMI case cloudy parcels may have temperatures lower than 270 K; these parcels do not exist in the MI simulation, in which the minimum temperature of parcels is 280 K. This effect is due to the fact that the entire cloud mixes with the warm entrainment air. Another example for the importance of the mixing process can be seen when examining the humidity histogram. In the NoMI case the humidity histogram spreads into high values reaching ~15 g m−3. These high values are associated with parcels that passed near the surface during the simulation; they ascend into the cloud layer but retain their high humidity. In the MI simulation the same parcels share the humidity with neighboring parcels in the process of turbulent mixing that leads to disappearing parcels with such a high humidity and to a shift in the location of the maximum of the histogram. Including the mixing process leads to a formation of realistic distributions of humidity and temperature, decreasing the number of extremely humid parcels as well as extremely dry ones. The histograms plotted in Fig. 12 indicate again that no-mixing clouds have an unrealistic thermodynamical structure and that the turbulent mixing is crucial for the proper representation of a stratocumulus topped boundary layer with an overlaying inversion layer.
Figure 13 shows histograms of main microphysical parameters: LWC and droplet concentration for MI and NoMI simulations. Here only parcels with a droplet concentration higher than 3 cm−3 are presented. By using this filter we remove all dry parcels from the histogram and show only the cloud layer. The histograms of LWC (top panels) are very similar, and values in both cases span from zero to the adiabatic value with all the variability expected in the cloud layer. A different picture arises when looking at the histograms for the droplet concentration (bottom panels). When mixing is not included the histogram depicts extremely high variability of drop concentration. Indeed, an isolated adiabatic parcel can be either droplet free or have different values of drop concentration depending on vertical velocity during CCN activation. In the NoMI case different parcels have different lifting condensation levels, which further leads to high variability of concentration. In the MI case the concentration histogram is narrow, demonstrating the homogenization of concentration by the mixing process. Analysis of observation data shows that the standard deviation of the concentration in this area is close to 20%–30%, which agrees with the width of the presented histogram. The exception is parcels with very low droplet concentration in the MI case. Analysis shows that these parcels are located at the cloud boundaries or around areas of downdraft with reduced LWC (refer to Fig. 3). It is interesting to note that there is low dispersion of concentration in spite of the small radius of correlation (i.e., of high spatial variability). Such homogenization of drop concentration might be attributed to the formation of a clearer cloud base with horizontally uniform temperature and humidity; as a result, most parcels have similar lifting condensation level. Parcels ascending from the same height will have nearly the same droplet concentration, which will change only slightly with height (say, via in-cloud nucleation). Such behavior differs from that of LWC, which increases with height from zero to maximum near the upper cloud boundary. Of course, mixing between parcels also tends to homogenize droplet concentration. As was discussed above, mixing within clouds is close to mechanical mixing.
f. Comparison of DSD properties and drizzle formation
In previous sections we looked at the larger scales—the averaged cloud properties—but one must also examine the smaller-scale changes that arise because of turbulent mixing. Much of the discussion about turbulent mixing revolves around the effect that this process has on the DSD. This topic is of great importance and will be investigated in future studies in more detail. Some initial results of the effects of mixing are presented in this section.
Figure 14 presents the effective radius distribution in the MI and the NoMI simulations. Here again only parcels with droplet concentrations higher than 3 cm−3 are presented to isolate the cloud layer itself. The white bars represent parcels that at t = 0 were located in the inversion layer but in the presented time steps have a droplet concentration higher than the applied filter. These parcels are entrained parcels. Two time instances are plotted: the first in the early cloud formation stage (top panels) and the second when the cloud is mature (bottom panels). In the early stages of the cloud the histograms in both simulations are similar in width with a slightly smaller mean peak value in the MI case. At the late stages of the cloud the differences between the two simulations becomes clearer. The distribution in the NoMI case does not change, indicating formation of a stationary state at which the histogram is centered at 10 μm. The distribution in the MI simulation is much narrower and the peak of the distribution is closer to 12 μm. From Fig. 14 it appears that parcels penetrated from the inversion layer into the cloud (lighter-colored bars) reach all values of effective radius and essentially become part of the cloud and indistinguishable from other parcels. No preferential formation of large droplets in the parcels penetrated from the inversion into the cloud layer is seen. Mixing leads to narrowing of the range of the effective radius variation. Low dispersion of the effective radius in stratocumulus, small cumulus, and deep convective clouds were reported earlier from observations (Gerber 2000; Gerber et al. 2008; Paluch and Knight 1984; Paluch 1986; Prabha et al. 2011; Freud et al. 2008; Freud and Rosenfeld 2012) as well as in numerical studies using spectral bin microphysics cloud models (Benmoshe et al. 2012; Korolev et al. 2013). Figure 14 demonstrates that this narrowing in the distribution of the effective radius is caused by in-cloud mixing. The processes leading to the low dispersion are unclear and require further research.
Figure 14 also shows that mixing leads to an increase (in average) in the effective radius in the cloud. The cloud measured during RF01 and simulated here was a nondrizzling cloud, and as such the mean values of the effective radius remained low and did not increase beyond the drizzling threshold of 12–13 μm. One possible mechanism to explain this result has been suggested by Korolev et al. (2013). They proposed that turbulent mixing between ascending and descending branches of oscillating parcels may lead to an increase in the effective radius. In zones of downdrafts, cloud droplets are larger than in updrafts because of several reasons, such as asymmetry in the diffusion growth–evaporation (Pinsky et al. 2013) or because of collisions between cloud droplets (Magaritz et al. 2010). According to the proposed mechanism, turbulent mixing can cause some drops from downdraft areas to enter updrafts’ areas (drop recycling). In case this penetration takes place in the lower part of the cloud where supersaturation is high, the drops will continue growing, fostering the formation of larger drops. Analysis of parcel trajectories in the model indicates that about 10% of parcels spend significant time within the cloud layer vertically oscillating from the lower to the upper cloud levels.
Figure 15 shows the cumulative precipitation at the bottom 200 m in three simulations: MI, MP, and NoMI. In all simulations the amount of drizzle that forms is very small and if measured will be negligible at the ground. This agrees with observations according to which clouds measured during RF01 were not drizzling ones. Nevertheless, there is a clear tendency of the formation of rain flux (large droplets) in simulations with mixing included. No drizzle reached the lower levels in the NoMI run.
One possible reason that no drizzle forms in the NoMI simulation is that the largest droplets cross droplet-free dry volumes during their fall, where they partially evaporate and lose their ability to serve as drop collectors. In the MI and MP simulations all parcels within the cloud are saturated and contain small droplets that can be collected by falling droplets.
An interesting note concerning turbulent mixing in the cloud layer arises when comparing the two mixing simulations: MP and MI. In the MP simulation all quantities are considered conservative and do not change during the mixing process (DSD in the filaments remains the same as it was in the parent parcel). Even though the comparison of the averaged vertical profiles did not indicate any significant difference, less drizzle reaches the ground in the MP case. Formation of drizzle is a fine process, where even small differences being accumulated may be of importance. With all other parameters equal between the simulations, when change in drop size during mixing is turned off (MP) the formation of drizzle is delayed by approximately 1.5 h, as compared to the case (MI) when droplet size changes during the mixing process.
The possible reason for the earlier formation of large drops in a case when the DSD is considered as nonconservative was discussed by Pinsky et al. (2010). In that study the evolution of DSDs in adjacent parcels located at different heights was simulated in the cases of nonpassive and passive mixing. It was shown that in the nonpassive mixing there was a broadening of the DSD toward larger sizes in the higher parcel, while in the lower parcel, mixing led to a shift of the DSD toward smaller sizes. In a case of passive mixing the DSD did not show a shift toward larger or smaller drop sizes. The results shown in Fig. 15 indicate that a continuous shift toward the larger sizes in the DSD in parcels located at the upper levels may result in an earlier formation of drizzle droplets.
4. Discussion and summary
This study represents the first part of a research dedicated to the investigation of thermodynamic and microphysical structures of stratocumulus clouds and to the formation of drizzle. In this paper we describe a new version of a novel hybrid Lagrangian–Eulerian model of stratocumulus clouds. The model represents the next important step in the development of trajectory ensemble models of stratocumulus clouds, described earlier in studies by Pinsky et al. (2008) and Magaritz et al. (2009). These steps were the development and implementation of the process of turbulent mixing between Lagrangian parcels and the inclusion of an inversion layer above the cloud. These two additions allowed simulating nonadiabatic interactions between cloudy parcels, as well as simulating the effects of entrainment and mixing on cloud thermodynamics and microphysics. To improve the energetic and moisture budget description in the model, heat and moisture surface fluxes and longwave radiation cooling were included.
There were two goals for this paper: first, to demonstrate that the model in its new configuration realistically reproduces the thermodynamical and microphysical structure of stratocumulus clouds, and second, to reveal the role of turbulent mixing in the formation of major thermodynamic and microphysical cloud properties. The model was tested using in situ measured data collected during flight RF01 of the DYCOMS II field experiment.
It was demonstrated that turbulent mixing between the parcels has a profound effect on cloud structure. Mixing forms a more uniform cloud base and cloud top and determines the spatial variability (correlation lengths) of many main variables such as LWC, droplet concentration, temperature, and humidity. Mixing leads to formation of nearly constant droplet concentration with height.
The model realistically reproduces the mixing diagram (the Paluch diagram), which characterizes the interaction between the cloud and inversion layer, indicating that entrainment and mixing processes are adequately represented in the model. At the same time, excluding the mixing between parcels makes the thermodynamic and microphysical cloud structure, as well as the mixing diagram, unrealistic. It is shown that when an inversion layer and entrainment are included in the model, mixing is the main factor affecting cloud structure, both at a macroscale and at a microscale.
A radius of correlation of a few hundred meters found in simulations with mixing is on the same scale as the radius calculated using the observed data. In the nonmixing simulations the radius of correlation is independent of the physical processes and is determined by the size of the parcel.
In the study presented by Magaritz et al. (2009), cloud structure (including the correlation lengths of all variables) and drizzle formation were realistically simulated in a nonmixing limit. The question arises: How were these results obtained if no-mixing simulations seen in this paper lead to unrealistic cloud structure? The answer lies in the simplifying assumption that the upper boundary of the computational area coincides with the upper cloud boundary. As a result, the entrainment of dry warm parcels from above was prohibited. Analysis of parcel tracks shows that most droplet-free parcels within clouds in the NoMI case (Fig. 3) arise because of penetration of dry parcels through the upper boundary. These penetrating parcels were excluded in the study by Magaritz et al. (2009), which led to an increase in the correlation lengths, decrease in variability of thermodynamic and microphysical parameters, and created more realistic results. From results presented in the current study we may conclude that the process of turbulent mixing is crucially important if one wishes to simulate stratocumulus clouds with an overlying inversion layer.
We would like to stress that the Paluch diagram as well as correlation lengths of different variables were simulated by the model using parcels (volumes) with characteristic linear scales of 25–40 m. The model results agree well with observations performed with a 25-Hz frequency. The mixing algorithm assumes that supersaturation in parcels is the same for all droplets, implying that results were obtained by applying an algorithm with homogeneous mixing properties. At the scales of 25–40 m, mixing might strictly be considered as inhomogeneous. However, high humidity within clouds and a comparatively rapid increase in humidity of entrained parcels makes these two types of mixing indistinguishable. In our opinion, this makes the utilization of such parcel sizes justified. In addition, the mixing algorithm is an extension of the well-known K theory (1.5 closure) that parameterizes effects of turbulent mixing via gradients of values taken at larger scales. This approach is appropriate for the consideration of parcels of this size. In addition, note that the assumption that supersaturation over all droplets is the same in a single-parcel does not exclude the ability of the model to simulate effects typically attributed to inhomogeneous mixing. For instance, first drops mixed into entrained parcels are fully or partially evaporated. The sequential ascent of such parcels, could, in principle, lead to formation of larger drops. This study indicates that turbulent mixing tends to increase the effective radius and make its distribution narrower. It is also shown that turbulent mixing fosters the formation of large droplets and creates facilitating conditions for drizzle to fall to the surface. We did not find any pronounced differences in the values of the effective radius between parcels entering from the inversion layer and mixed into the cloud and other cloud parcels.
It was also shown that while simulations with passive or nonpassive mixing show similar averaged values, finer-scale processes affect the rate of drizzle formation. Disregarding the nonconservative properties of DSD during mixing inhibits the formation of drizzle in the cloud. It is important to note that in the present study the nondrizzling clouds observed in RF01 were simulated. This case is not very suitable for investigation of processes of drizzle formation. We plan to test the proposed hypotheses concerning the effects of mixing on drizzle formation in simulations of drizzling clouds—for instance, those observed in RF07 of DYCOMS II. We will reexamine the concept of lucky parcels introduced by Magaritz et al. (2009), according to which drizzle forms first in parcels that ascend from the ocean surface, and spend significant time near cloud top, where the LWC is maximum, and see how turbulent mixing and entrainment may affect these results.
The study is supported by the Office of Science (BER), U.S. Department of Energy Award DE-SC0006788 and the Binational U.S.–Israel Science Foundation (Grant 2010446). DYCOMS II data used in this study were obtained from the NCAR EOL database.