Long-resident droplets at the stratocumulus top

Turbulence models predict low droplet-collision rates in stratocumulus clouds, which should imply a narrow droplet size distribution and little rain. Contrary to this expectation, rain is often observed in stratocumuli. In this paper, we explore the hypothesis that some droplets can grow well above the average because small-scale turbulence allows them to reside at cloud top for a time longer than the convective-eddy time t. Long-resident droplets can grow larger because condensation due to longwave radiative cooling, and collisions have more time to enhance droplet growth. We investigate the trajectories of 1 billion Lagrangian droplets in direct numerical simulations of a cloudy mixed-layer configuration that is based on observations from the flight 11 from the VERDI campaign. High resolution is employed to represent a well-developed turbulent state at cloud top. Only one-way coupling is considered. We observe that 70 % of the droplets spend less than 0.6t at cloud top before leaving the cloud, while 15 % of the droplets remain at least 0.9t at cloud top. In addition, 0.2 % of the droplets spend more than 2.5t at cloud top and decouple from the large-scale convective eddies that brought them to the top, with the result that they become memoryless. Modeling collisions like a Poisson process leads to the conclusion that most rain droplets originate from those memoryless droplets. Furthermore, most long-resident droplets accumulate at the downdraft regions of the flow, which could be related to the closed-cell stratocumulus pattern. Finally, we see that condensation due to longwave radiative cooling considerably broadens the cloud-top droplet size distribution: 6.5 % of the droplets double their mass due to radiation in their time at cloud top. This simulated droplet size distribution matches the flight measurements, confirming that condensation due to longwave radiation can be an important mechanism for broadening the droplet size distribution in radiatively driven stratocumuli.


Introduction
Stratocumulus clouds are characterized by small droplets (radius, R ∼ 10 µm) and by low turbulence intensities (dissipation rates, ∼ 20 cm 2 s −3 ), which implies low collision rates between cloud droplets (Grabowski and Wang, 2013;Onishi et al., 2015).For microphysical models, low collision rates imply slow autoconversion (where autoconversion is the process in which cloud droplets collide forming rain droplets), narrow droplet size distributions (DSD), and sporadic rain (Rogers and Yau, 1989).Contrary to this expectation, stratocumulus rain is actually quiet commonly observed (Leon et al., 2008).This apparent paradox is known as the size gap problem, and it is relevant for modeling rain in most low clouds (Grabowski and Wang, 2013).There are several possible explanations for the size gap paradox, such as the existence of giant aerosols or fluctuations in the condensation process, but the debate is still open (Wang, 2013).The insufficient knowledge of what causes stratocumulus rain explains why current parameterizations of precipitation rates produce quite different results (Wood, 2012).
In this paper, we explore the hypothesis of whether stratocumulus rain formation can be enhanced if a small fraction of the droplets spend long residence times at the cloud top, as proposed, for example, by Stevens and Feingold (1996).Stratocumulus layers can exist for several days, and there is the possibility that some droplets remain close to the cloud top for long times.Given that droplet growth due to condensation by longwave radiative cooling and due to collisions reaches its peak at the cloud top, long-resident droplets could grow well above the average size and lead to the formation of rain droplets (Roach, 1976;Barkstrom, 1978).In order to stay close to the cloud top, long-resident droplets need to escape the large-scale convective eddies that drive them out of the cloud.This is possible due to the chaotic nature of tur-Published by Copernicus Publications on behalf of the European Geosciences Union.
A. de Lozar and L. Muessle: Long-resident droplets at the stratocumulus top bulence, which introduces some apparent randomness in the droplet motion.Similar hypotheses based on the chaotic nature of turbulence are the eddy-hopping mechanism (Cooper, 1989;Grabowski and Wang, 2013), which predicts that rain formation is enhanced because turbulence mixes droplets with different condensational growth histories, and the recycling mechanisms (Kogan, 2006;Naumann, 2015;Naumann and Seifert, 2016), which predict that rain droplets are mostly nucleated in lucky parcels/droplets that recirculate through the cloud and thus have more time for collisions.Also similar to these hypotheses, Yang et al. (2013Yang et al. ( , 2015) ) have recently proposed that large ice crystals in stratiform mixed clouds originate from ice particles that either recirculate through the cloud or are trapped in the updraft regions.
The enhancement of rain in parcels with different histories has been studied with trajectory ensemble models (TEM).In those models, the evolution of the droplet size distribution (DSD) is calculated in several parcels that are driven by an external model, typically a large-eddy simulation (LES).The parcels vary their position following the flow dynamics, so that the TEM can be considered as a Lagrangian model.Harrington (2000) and Hartman and Harrington (2005) used a TEM driven by a LES to show that stratocumulus drizzle is sped up by radiative cooling at parcels that stay close to the top.Similarly, Magaritz et al. (2009) used a TEM driven by two-dimensional synthetic turbulence and found that drizzle is initiated in lucky parcels that stay for a long time close to the cloud top.All these TEM show that rain is mostly initiated at cloud top, but they have three main shortcomings for a proper quantification of the effect.First, parcel mixing is not accounted for.Real parcels mix with the environment and lose their identity after a time comparable to the convective-eddy time, as most droplets that were initially in the parcels have already left.Neglecting mixing causes an overestimation of the fluctuations in the liquid field and a tooearly rain formation (Magaritz-Ronen et al., 2014).Second, most cases are based on a small number of parcels (∼ 1000), which do not allow the investigation of extreme behaviors that might be important for rain formation.Third, turbulence is usually not accurately modeled because it is computationally too expensive to solve in detail turbulence and the microphysical processes that determine the evolution of the DSD at the same time.This lack of accuracy can be an important issue for the study of long-resident parcels/droplets, giving that they rely on turbulence to stay at the cloud top.
In this paper, we use direct numerical simulations (DNS) to quantify the time that droplets reside at the top of a cloudy mixed layer driven by radiative and evaporative cooling.Our purpose is to complement past Lagrangian studies (Stevens and Feingold, 1996;Kogan, 2006) by focusing on the turbulent cloud-top region, where rain droplets are expected to be formed.As in these past studies, we restrict to small cloud droplets and neglect sedimentation, which could have a relevant effect for larger droplets.In our DNS, we explicitly solve part of the complex turbulent dynamics at the cloud-top boundary, which we think is necessary for a proper quantification of the cloud-top resident times.This investigation is partly motivated by recent studies in the convective boundary layer (Mellado et al., 2015), which have shown that eddies of multiple sizes are relevant for the flow dynamics close to the boundary where convection is initiated.In addition, our Lagrangian scheme follows 1 billion individual droplets and therefore provides better statistics for extreme cases than in previous studies.The plausibility of our results is tested by comparing the DSD from our simulations with observations from the VERDI campaign.
2 Formulation in a mixed layer Today's computational resources do not allow us yet to simulate at the same time all physical processes relevant for the stratocumulus dynamics plus the turbulent flow length scales that are relevant for the droplets' movement.Since the focus of this paper is on turbulence, we decided to use a relatively high resolution to simulate a simple configuration that mimics the stratocumulus dynamics: a cloudy mixed layer.This configuration consists of a cloud layer (moist and cool) that lays below a layer of dry and warm air that represents the free atmosphere.Both layers are characterized by their total water content q t and temperature T , and are at the same thermodynamic pressure p.The cloud layer is bounded at the bottom by a strong stratification that defines the cloud base (this is the main variation from the cloud-top mixing layer configuration that we used in previous entrainment studies; see de Lozar and Mellado, 2015a).Cloud-top mixing is driven by radiative and evaporative cooling.Other driving mechanisms like mean cloud-top shear or cloud-base fluxes are neglected.The longwave radiation is characterized by the divergence of the radiative flux at cloud top F 0 , and by the radiative-extinction length, λ, that defines the region cooled by radiation (Larson et al., 2007).
The flow dynamics are calculated by solving the evolution equations for momentum, total water, and enthalpy on a fixed grid.We use the name Eulerian to refer to this first part of the calculations.We also track the evolution of 1 billion droplets that follow the tendencies dictated by the Eulerian calculations.We use the name Lagrangian to refer to this second part of the calculations.We consider only one-way coupling, which means that the Eulerian calculations are independent from the Lagrangian ones.This means that the Lagrangian calculations can be considered as complex diagnostics statistics from the Eulerian simulation.In this section, we present the Eulerian and Lagrangian formulations employed in this study.

Eulerian formulation
We use the Eulerian formulation described in de Lozar and Mellado (2015a) and summarized in Appendix A. This for-mulation uses the linearized forms of the buoyancy and saturated-vapor functions, infinitely fast thermodynamics, the Boussinesq approximation, and Prandtl and Schmidt numbers equal to 1.
We define the normalized liquid water as = q l /q c l , where q l is the liquid water content and q c l is the initial liquid water content in the cloud layer.A diagnostic exact equation for can be obtained from the equations introduced in de Lozar and Mellado (2015a): where ν is the kinematic viscosity of air, s cond is a source term due to the condensation of liquid due to longwave radiation, and s evap is a sink of liquid due to evaporation.The exact forms for the source and sink terms are given in Appendix A1.
The first term on the right-hand side of Eq. ( 1) represents the mixing of liquid content due to diffusivity.In the original exact formulation the diffusivity of liquid droplets originates from the thermal fluctuations in the droplets' movement, and it is negligible when compared to the other terms of the equation (Mellado et al., 2010;de Lozar and Mellado, 2014).However, this term is not negligible in our formulation due to the assumption that the liquid phase diffuses like heat or momentum.This assumption is partly justified when using a much higher viscosity than in air because the diffusive mixing then captures the mixing by the unresolved turbulence, similarly as the subgrid diffusion in LES.This approximation is well justified when the statistics of interest are independent of ν, which is an indication that mixing is independent of the small scales, as it is often observed in turbulent flows (Dimotakis, 2005) (Reynolds number independence).

Lagrangian formulation
Parallel to the Eulerian calculations, we track the trajectories and mass of 1 billion individual droplets.We call these droplets Lagrangian because the equation of motion is solved independently for each droplet.The Lagrangian droplets represent only a small fraction of the number of cloud droplets in a real cloud, and are introduced to obtain extra information from the simulations.
We impose that the Lagrangian droplet dynamics follow the mean local tendencies given by the velocity and source terms from the Eulerian formulation.This assumption neglects all fluctuations that arise from length scales between the simulated Kolmogorov length scale (η sim ∼ 40 cm) and the atmospheric Kolmogorov length scale (η atm ∼ 1 mm).
The assumption is consistent with DNS, where the smallest resolved length scale is comparable to η sim .In Appendix B, we investigate the effect of neglecting fluctuations below the smallest simulated length scale, by comparing simulations with lower viscosities in which the Kolmogorov length scale is reduced to η sim = 14 cm.
We initially place each Lagrangian droplet at a random position inside the cloud (defined by > 0).These cloud droplets move as tracers, following the Eulerian-flow velocity: where X i denotes the position of the Lagrangian droplet (i = 1, 2, . .., 10 9 ), and v(X i ) is the Eulerian velocity at the position X i .In general, we use capital letters for Lagrangian properties and lowercase letters for the Eulerian properties.Equation ( 2) is consistent with the assumption that each droplet moves with the mass-averaged mean velocity of a volume given by the grid size, and neglects all unresolved turbulent fluctuations.This assumption also neglects the effect of inertia and gravitational settling on the droplets' movement.
In order to study the evolution of the droplet size distribution, we assign a mass M i to each Lagrangian droplet.We define a normalized mass for each droplet as L i = M i /M r , where M r is the mean initial mass of a droplet in the cloud layer.The mass of each droplet at the initial condition is extrapolated from the Eulerian field.Analogous as with the velocity, we impose that the evolution of the mass of each Lagrangian droplet follows the mean condensation and evaporation imposed by the Eulerian calculations: where the evaporative and condensation tendencies are again calculated at the particle's position.The change of droplet's mass due to collisions cannot be extrapolated from the Eulerian tendencies, and has not yet been included in this formulation.Again Eq. ( 3) neglects all fluctuations in the liquid field that arise from the unresolved scales.Furthermore, we set the mass of each droplet to 0 if the normalized Eulerian liquid field at the droplet position (X i ) is lower than 10 −5 .This clipping eliminates Lagrangian droplets that exit the cloud without fully evaporating as a consequence of the different tendencies of the Eulerian and Lagrangian fields.Droplets grow as they ascend through the clouds and shrink when they descend, due to the change of temperature in the atmospheric boundary layer.Equations ( 1) and (3) neglect the condensation/evaporation of droplets with changing height because this is identically zero in a mixed layer with constant pressure.We decided not to include this process because its effect on the cloud-top DSD of a well-mixed stratocumulus is also negligible when assuming that all droplets follow the mean condensation/evaporation.The reason is that this assumption implies that condensation/evaporation affects equally all droplets independently of their size, and all of them grow by the same amount when they reach the cloud top.In reality, larger droplets condensate/evaporate faster than smaller ones leading to some broadening of the DSD (Cooper, 1989;Lanotte et al., 2009), but this mechanism is not considered in this study.
The Lagrangian equation for the liquid mass (Eq. 3) is equivalent to the Eulerian one (Eq.1) except for the diffusion term, which is neglected in the Lagrangian formulation.While in the Eulerian formulation this term represents the mixing of mean liquid water content at the smallest scales, there is not an equivalent physical process when the Lagrangian particles represent single droplets.In order to explain this concept better, let us consider two saturated parcels with same volume, same number of droplets, same temperature, but different liquid water content.We also consider that both parcels are monodisperse, so that the parcel with higher liquid water content contains larger droplets.If the number of droplets in each parcel is high enough, small-scale turbulence will mix these two parcels until the liquid water content in both parcels is the same.In the Eulerian formulation, this process is modeled by the diffusion term, which homogenizes the liquid water content.As both parcels are initially at the same temperature, there is no condensation/evaporation and individual droplets do not change their size in this process.The composition of each parcel is now polydisperse, containing small and large droplets.This information about the DSD is not captured by the Eulerian bulk formulation, but it is retained by the Lagrangian formulation without the diffusion term.

Caveats of the Lagrangian formulation
Our simulations show too little evaporation of the Lagrangian droplets, when compared with the Eulerian mean values (reduction up to 90 %).This strong difference is caused by neglecting the diffusion term in Eq. (3).Extra calculations in which the diffusion term was included in the Lagrangian tendencies show little differences from the Eulerian results, proving that numerical errors are much smaller than the differences introduced by the diffusion term.
In order to explain the differences introduced by diffusion, we sketch in Fig. 1 a mixing process between two parcels with saturated and unsaturated air.The left and right parts of Fig. 1 represent parcels before and after the mixing process, respectively.In the top-right diagram, we represent the mixing process according to the Eulerian formulation (with diffusion).First, vapor, heat, and liquid water diffuse until both parcels are unsaturated.Next, cloud droplets evaporate in both parcels.In the bottom-right diagram, we represent the mixing process according to the Lagrangian formulation.Here, droplets do not diffuse with heat and vapor and remain in the originally saturated parcel.Next, the Lagrangian droplets evaporate following the Eulerian tendencies in the originally saturated parcel, but they do not evaporate in the second parcel because there are no Lagrangian droplets there.Consequently, the average evaporation of Lagrangian droplets is always weaker than the mean evapora- tion provided by the Eulerian tendencies.This mismatch introduces an uncertainty in the DSD, specially at the small droplets which arise from the evaporation process.These small droplets are not the main goal of this investigation, but nevertheless this problem should be improved in further Lagrangian studies that use the formulation introduced in this paper.

Reference case
Our reference case is based on the diurnal measurements from the reference flight 11 in the VERDI campaign, taken on 15 May 2012 in the Beaufort Sea area (Klingebiel et al., 2015).The measurements show a solid stratocumulus deck with very low ice particle number concentration N ice = 0.75 L −1 .The measured profile of liquid potential temperature is not perfectly well mixed, and shows some indications of decoupling at cloud base (θ top l − θ base l 3 K).The DSD was measured with a cloud droplet probe and a cloud imaging probe, which cover the size range of about 2-960 µm.The droplet number concentration is approximately constant during the flight across the stratocumulus N = 70 cm −3 .
The cloud and dry layer in the simulations mimic the cloud-top thermodynamical properties measured by this flight: cloud layer at temperature of T = −5 • C, with total water content q t = 3.15 g kg −1 , and liquid water content q C l = 0.25 g kg −1 ; and free atmosphere at T = −0.5 • C, with total water content q t = 2.4 g kg −1 .We do not consider the ice phase in the simulations.The radius of cloud-top droplets with average mass M r is R = (3ρ a q C l ) 1/3 (4π N ρ l ) −1/3 = 10 µm, where ρ a and ρ l are the densities of air and water, respectively.The atmospheric pressure at cloud top is 905 hPa.The cloud-top longwave net radiative flux F 0 = 60 W m −2 matches the measured data of a second airplane that flew on the top of the reference flight 11.This radiative flux defines a reference buoyancy flux B 0 = (F 0 g)/(ρC P T ) = 17 cm 2 s −3 that serves to scale the entrainment velocity and turbulence dissipation rate in the cloud layer (de Lozar and Mellado, 2015b).The longwave radiative extinction length is chosen as λ = 15 m.
The reference case is characterized by a relatively wet free atmosphere that limits the impact of evaporation on turbulence ( de Lozar and Mellado, 2015b).A weak evaporativecooling forcing relaxes the resolution requirements to solve the turbulence drivers (see Appendix B for details), which allows us to explore larger configurations and longer simulation times than in our previous entrainment studies ( de Lozar and Mellado, 2015a).The radiative parameters and stratification are quite typical for stratocumuli.For example, the cloud-top net radiative flux in most test cases reviewed by Stevens (2002) are in the interval 50 W m −2 < F 0 < 74 W m −2 .For this reason, we expect that the investigation of the reference case can serve to understand droplet dynamics in other stratocumuli that are mostly radiatively driven.

Flow description
The main simulation was performed in a horizontally periodic mixed layer 750 × 750 m 2 wide.The cloud base is placed z * = 300 m below cloud top, where we impose a stratification that does not allow for a deep penetration of the flow (see Appendix A2 for details).The flow is thus confined to the vertical extent of the cloud, like in a strongly decoupled stratocumulus-topped boundary layer.Consequently, the largest eddies, here called convective eddies, have a fixed size that is comparable to the cloud depth z * .The convective-eddy velocity is characterized by the integral velocity w * = (2.5 w b dz) 1/3 = 0.67 m s −1 introduced by Deardorff (1970).This defines the convective-eddy time t * = z * /w * = 447 s.The simulation runs for 83 min, which correspond to 11.1t * .The Lagrangian droplets are initiated at the beginning of the simulation at t = 0.The kinematic viscosity in the DNS is such that the reference Reynolds number is Re 0 = B 1/3 0 λ 4/3 ν −1 = 200.The resulting Kolmogorov length scale (η sim ∼ 40 cm) is considerably smaller than the length scales that characterize radiative cooling and entrainment, both of the order of ∼ 15 m ( de Lozar and Mellado, 2013;Gerber et al., 2013).Finally, the integral Reynolds number Re * = w * z * /ν = 8770 and the Reynolds number based on the Taylor microscale Re λ 200 are large enough to expect Reynolds number independence in many turbulence statistics.

Numerics
The numerical algorithm is based on high-order, spectrallike compact finite differences (Lele, 1992) and a low-storage fourth-order Runge-Kutta scheme (Carpenter and Kennedy, 1994).The time step is set by a Courant condition.The main simulation described in previous section runs for 18 900 time steps.Additional simulations described in Appendix A2 run for ∼ 3500 time steps.The pressure Poisson equation is solved using a Fourier decomposition along the periodic horizontal planes and a factorization of the resulting set of equations along the vertical coordinate (Mellado and Ansorge, 2012).
All the simulations were done in a cubic grid with 1024 3 points.The resolution parameter x/η sim is of the order of 2.0 or less, where x is the grid spacing.Using grid convergence studies (not shown), such a resolution has been proved to be enough for accuracies of the order of 2 % or better in the statistics discussed in this paper, using the numerical algorithm described above.Further details can be found in Mellado (2010).
Eulerian fields are extrapolated into the Lagrangian droplets' positions by using a trilinear interpolation from the eight closest grid points.The Lagrangian equations are integrated in time with the same low-storage fourth-order Runge-Kutta scheme (Carpenter and Kennedy, 1994) that we use in the Eulerian calculations.The Lagrangian calculations are parallelized.When using 10 9 Lagrangian droplets in 10 9 grid points, the calculation time in the Eulerian and Lagrangian part of the code is almost equal.The source code can be found in de Lozar (2016).

Trajectories
We first look at the trajectories of 20 Lagrangian cloud droplets during 1.2 t * , which were selected for visualization purposes.The tracking procedure is started after running the simulation for 5 t * , thus well after the initial transient.We classify the trajectories in three different categories (red, green, and orange) which are described in the following paragraphs.
Red trajectories correspond to cloud droplets that circulate quickly through the cloud, and spend a short time at cloud top.Figure 2 suggests that these droplets follow large convective eddies, which motivates us to call them convective droplets.In our simulations "convective droplets" recirculate in the cloud due to the stratification at cloud base, but we can expect that these droplets will exit the cloud and evaporate in well-mixed stratocumulus boundary layers.
Green trajectories correspond to cloud droplets that escape the downwards movements of the convective eddies, and stay at the cloud-top region for the whole tracing period.We use the name "long-resident droplets" for droplets that follow these trajectories.Long-resident droplets grow on average more than convective droplets, although they are less in number.Figure 2b shows that long-resident droplets do not stick to the inversion but feature a turbulent movement characterized by multiple length scales.This picture suggests that small-scale turbulence is important to escape the convective eddies.
Orange trajectories correspond to droplets that are driven out of cloud top by intermediate eddies, so that they keep circulating through the cloud bulk.Their near future is still undecided.They can either go back to the cloud top, and eventually become a long-resident droplet, or be driven downwards by a convective eddy.We use the name "erratic droplets" for those droplets that are neither convective nor long-resident.This is indeed poor naming because erratic droplets can have quite different dynamics, but it just shows how difficult it is to characterize the movement of all droplets in a turbulent flow.

Residence time
In order to quantify the different droplet dynamics we assign to each droplet a cloud-top residence time t res that is defined by the following rules: -The residence times are initialized 6.6t * after the start of the simulations, thus well after the initial transient.
-At each time step, the residence time advances only for droplets that are at a distance of less than 2λ from the maximum buoyancy gradient at the initial condition.The region where the residence time advances defines the cloud-top region.
-The residence time of cloud droplets that are distanced 10λ from the maximum buoyancy gradient is set to zero, i.e., we delete their memory.Since the cloudy domain comprises 20λ, this means that half of all cloud droplets lose their memory at each time step.
-The residence time of the cloud droplets whose normalized mass is below 10 −5 is set to 0. This condition eliminates the completely evaporated droplets from the residence time statistics.
The objective of this definition is to discriminate longresident droplets from convective ones.Convective droplets circulate through the mixed layer in a time comparable to the convective-eddy time, and their maximum residence time is a fraction of t * (Stevens and Feingold, 1996).On the other hand, long-resident droplets are able to escape the convective eddies and can reach considerably longer resident times.Erratic droplets are associated with very diverse dynamics and could have effectively almost any residence time, which complicates the discussion.
Figure 3a shows the evolution of the residence time histogram p(t res ).The histogram collapses for the short residence times, indicating that this statistic quickly reaches a quasi-stationary state.The sharp increase at the right end of the histograms shows that a considerable amount of droplets have the longest residence time of the distribution.This is an artifact of the finite time sampling because the longest residence time is always equal to the time elapsed between the measurement time and the residence time's initialization.As time advances, the longest residence time increases and the right part of the histogram develops.The right part of the histograms also collapses for different times, thus showing a quasi-stationary shape.
Figure 3a shows a plateau for t res < 0.6t * , which indicates that most droplets spend around this time at cloud top.This is consistent with the convective movements that drive most droplets to spend some time, comparable but smaller than the convective-eddy time, at cloud top.We identify the droplets with t res < 0.6t * with the previously defined convective droplets.According to this definition, 70 % of all  4).The colors stand for the time elapsed after initiating the residence times: t * , 2t * ,3t * ,4t * ,4.5t * .
droplets that have reached the cloud top behave as convective droplets.After t res ∼ 0.6t * , the residence time quickly decays.In order to investigate the quick decay of the residence times, we define a decay time, τ res , by: which, in case of an exponential decay, is constant and equal to the mean lifetime.The decay time in Fig. 3b is calculated by taking the derivative of the histogram presented in Fig. 3a, after using a Savitzky-Golay low-pass filter.We distinguish different regimes of τ res in Fig. 3b, which are discussed in the next paragraphs.
In the interval 0.6t * < t res < 0.9t * , the histogram and decay times rapidly decrease (1/τ res increases).Droplets with these residence times can be convective droplets that were driven by slow convective eddies, or erratic droplets that recirculate through the cloud.
In the interval 0.9t * < t res < 2.5t * , the histogram decays quasi-exponentially with a decay time that slowly increases from τ res t * /3.2 to τ res t * /2.In other processes like radiative decay, an exponential decay of the lifetime histogram means a random decay, and we use the same interpretation even when the definitions of lifetime and residence time are not identical.We explain the quasi-exponential decay of the histogram due to several droplets escaping the convective eddies that initially brought them into the cloud top.Once detached from the convective eddies, more unorganized cloudtop turbulence drives the droplets' movement, and their residence time becomes close to random.However, the small variation of the decay time in this part of the histogram suggests that these droplets can still feel the initial eddies, even when they are not directly driven by them.We identify these droplets with the long-resident droplets introduced in the previous section.Using this definition, we observe that 15 % of the droplets are long resident.
From t res > 2.5t * the number of droplets decays exponentially with τ res t * /2.These droplets become memoryless, in the sense that τ res does not vary for the later times investigated in our simulations.The memoryless condition implies that these droplets have completely lost their connection to the convective eddies that initially brought them to the cloud top.In fact, flow visualizations suggest that most convective eddies vanish after a time comparable to 2.5t * .We can think of droplets that wander through the cloud top region until being dragged downwards by some new eddy.This behavior can be identified as eddy hopping.Using this definition, we observe that 0.2 % of the droplets are memoryless longresident droplets.
Next, we look at where long-resident droplets can be found.In order to visualize our results, we first extrapolate the resident times to the Eulerian grid, for which we use the trilinear interpolation.Figure 4 shows the extrapolated residence times larger than 0.9t * .This value signals the presence of long-resident droplets, as explained in previous paragraphs.The plot is very patchy, as a result of the continuous mixing between long-resident and new convective droplets.Long-resident droplets can be found everywhere but they clearly prefer the downdraft regions of the flow.Updrafts are mainly composed by droplets that originate from the cloud bulk, which in general have shorter residence times.Within the downdraft regions in Fig. 4 there are also preferred locations with high densities of long-resident droplets.Those regions are probably characterized by broader DSDs, with higher probabilities for the collision-coalescence process.We can thus speculate that the generation of rain droplets is more likely in these regions than anywhere else in the cloud.

Impact of condensation due to radiative cooling for the cloud droplet size distribution
The different droplet dynamics discussed in previous sections should be reflected on the DSD.In particular we expect that long-resident droplets grow larger than average due to the combined action of radiative cooling and collisions, thus broadening the DSD.In this section, we investigate how the DSD broadens due to condensation induced by radiative cooling alone.Droplet evaporation is also included in the calculations, but we expect it to have little effect on the DSD be-  cause our Lagrangian algorithm considerably underestimates the evaporation of droplets (see Sect. 2.3).
Radiative and evaporative cooling in our simulation overcome the influx of heat at cloud top.As a result the cloud continuously cools down and the total liquid water linearly increases with time.This introduces a continuous shift in the DSD equal to Gt, where G is the mean growth of a droplet in the cloud.For a better comparison at different times, we subtract this mean growth from the DSDs presented in this section.The DSDs are normalized by the initial mean mass of the droplets M r .
Figure 5 shows the evolution of the DSD at the top 2λ in our simulations.We observe a significant broadening of the DSD, which does not seem to reach a steady state.At the last time step, 6.5 % of the droplets have doubled their initial mass.These are mostly long-resident droplets, as droplets need to grow for 1.25t * at maximum radiative cooling to double their mass in our setup.However, further growth becomes difficult, which can be explained by the exponential decay of the resident times.Only 5 in 1 million droplets triple their size, although this number seems to increase as time advances.
We compare the DSDs in our simulations to the in situ measurements from the VERDI campaign described in Sect.3.1.As with the simulations, we use only measurements close to the cloud top.The measured DSD has been rescaled in the vertical coordinate to properly compare the different sample sizes in observations and simulations.Figure 5 shows a good similitude between observations and simulations for the small sizes of the DSD ((M − Gt) 2.5M r ).We find this similitude rather surprising giving the large number of simplifications in our model, and that condensation due to radiative cooling is the only mechanism for the DSD broadening.Since our model is more appropriate for detached stratocu- muli, one possible reason for the good agreement is that the cloud in the observations was detached at cloud base, as suggested by the measured profile of liquid potential temperature.Unfortunately, we do not have observations below cloud base to verify this hypothesis.The agreement between observations and simulations is even more puzzling for the very small sizes of the DSD ((M − Gt) 0.5M r ), given the uncertainty in our evaporation calculations.For larger droplets ((M − Gt) > 2.5M r ) observations and simulations clearly differ, indicating that other processes, like collisions, are important for the further broadening the DSD.Even when considering that there is some degree of chance in the agreement between observations and simulations (such as in the final time of the simulations), the comparison strongly suggests that condensation due to longwave radiation was an important factor for broadening the DSD in the stratocumulus observed in the VERDI campaign.Since the radiative properties of this cloud are quite common, we conclude that longwave radiation is probably also relevant for the DSD evolution in many other radiatively driven stratocumuli.

Discussion and conclusions
Our simulations show cloud droplets that reside at the stratocumulus top for times considerably longer than the convective-eddy time, supporting some previous LES studies (Stevens and Feingold, 1996;Kogan, 2006;Magaritz et al., 2009).Here, we confirm that long-resident droplets are a result of the complex convective-turbulent dynamics, and are not an artifact of the turbulence model.This assertion is based on the resolution used in our simulations ( x = 73 cm), which is high enough to explicitly resolve a fully turbulent cloud top.The setup of our simulations mimics observations of Arctic stratocumuli during the VERDI campaign (Klingebiel et al., 2015).Turbulence is driven by longwave radiative cooling and to a lesser extent by evaporative cooling.Since the properties that determine the radiative forcing in the simulations are quite common for stratocumulus clouds, we expect that long-resident droplets can also be found in many other radiatively driven stratocumulus-topped boundary layers.
We observe that 15 % of the droplets partially escape the stratocumulus large-scale convective motions and reside at the cloud-top region for periods longer than the convectiveeddy time.This percentage is in rough agreement with cloudtop resident times in a drizzling stratocumulus LES (Kogan, 2006), and it is approximately twice as large as in nondrizzling stratocumulus LES (Stevens and Feingold, 1996).The comparison is not perfect because each study uses a different definition for the cloud-top residence time, but we find the approximate agreement of our simulations with LES with much coarser resolutions ( z = 25 m, x ≥ 50 m) rather surprising.
The large number of droplets in our simulations (10 9 ) allows for accurate quantification of extreme cloud-top residence times.We find that 0.2 % of the droplets spend more than 2.5t * at cloud top and completely decouple from the large-scale convective movements.The process of leaving cloud top becomes memoryless with a decay time τ res ∼ t * /2.Assuming that the separation of scales between the radiation length scale λ and the boundary layer depth z * is large enough, our results in the mixed-layer model can be roughly extrapolated to deeper well-mixed boundary layers: for a ∼ 1 km deep well-mixed radiatively driven stratocumulustopped boundary layer with t * = 20 min, we expect that 15 % of the droplets remain at cloud top for more than 20 min, and that 0.2 % of the droplets remain for more than 50 min and become memoryless with a decay time of 10 min.These residence times are probably long enough to allow for some long-resident droplets having multiple collisions at cloud top, as we discuss in the next paragraph.
Let us provide a simple argumentation to show how longresident droplets can experience multiple collisions and contribute to the formation of rain droplets at cloud top.The gravitational-settling collision rate for droplets of radius R = 10 µm with very small tracer droplets is p = 0.003 min −1 (Hall, 1980).For a typical convective-eddy time t * = 20 min, convective droplets that spend 10 min at cloud top have pt res = 0.03 collisions in mean value.Given that pt res 1, the probability of a single droplet colliding at cloud top is also 0.03.For memoryless long-resident droplets that spend 50 min at cloud top, the probability of colliding with a single droplet is 0.15, as it scales with time.This factor of 5 increase seems insignificant because it does not compensate for the fact that there are 500 convective droplets for each one that spends 2.5t * or longer at cloud top.The potential of long-resident droplets can be understood by looking at subsequent collisions.For simplicity let us assume that the probability of each subsequent collision is roughly equal to probability of the first collision.Although in reality, subse-quence collisions are more likely, the qualitative argumentation that follows is not affected by this simplification.Under this assumption, collisions become a Poisson process, in which the probability of a droplet having N collisions scales as ∼ (pt res ) N when pt res 1. Long-resident droplets are then 25 times more likely to have 2 collisions and this factor scales quickly with N. For N = 4, the enhancement factor is 625, which means that there are more long-resident droplets that had 4 collisions than convective ones.Cloud droplets require more than 4 collisions to become rain droplets, which are defined by a limiting radius R > 25-45 µm depending on the accretion model (Khairoutdinov and Kogan, 2000;Seifert and Beheng, 2001).We conclude that long-resident droplets can be relevant for rain formation.An important consequence is that rain models for cumulus clouds (where most droplets are convective), do not need to necessarily be valid for stratocumuli, in which rain could predominately originate from long-resident droplets.
We observe that long-resident droplets prefer to concentrate in the thin downdraft regions of the flow.This behavior agrees with the pattern of in-cloud residence times (different from the cloud-top residence times in this paper) found by Kogan (2006) inside drizzling stratocumuli.The flow in the LES of Kogan (2006) is characterized by thin updrafts and broad downdrafts, which is the opposite pattern of our simulations.The agreement suggests that the tendency of long-resident cloud droplets to be preferably placed in downdrafts is generic for droplets with relatively small sedimentation velocities (in the other extreme, long-lifetime ice crystals with much higher sedimentation velocities cluster in the updrafts regions Yang et al., 2015).As a consequence, longresident droplets are placed close to the small droplets generated in the entrainment process, which are also advected to the same convergence regions (as observed, for example, in Yamaguchi and Randall, 2012).Downdraft regions are thus characterized by a more polydisperse distribution than the updrafts, in which new droplets are brought to cloud top.As a consequence, the probability of collisions and the rate of formation of rain droplets are higher in the downdrafts than in the updrafts regions.We speculate that the resulting depletion of liquid in the downdrafts could help to explain the closedcell pattern which is commonly observed in stratocumulus clouds (Wood, 2012).This mechanism resembles the one introduced by Ovchinnikov et al. (2013), who consider that the closed-cell pattern is a consequence of convective droplets producing rain in a time comparable to the convective-eddy time.The main difference with Ovchinnikov et al. (2013) is that we hypothesize that long-resident droplets, and not convective droplets, are relevant for the generation of rain.Our hypothesis leads to a more frequent formation of the closedcell pattern.
Long-resident droplets also grow on average considerably larger than convective droplets due to condensation induced by longwave radiative cooling, which can speed up rain formation (Roach, 1976;Hartman and Harrington, 2005).In this paper, we isolate the effect of radiation, and investigate how condensation due to radiative cooling alone broadens the cloud-top DSD.We observe a significant broadening of the DSD: at the last time step 6.5 % of the droplets have doubled their initial mass.This broadening is probably wider for deeper stratocumulus-topped boundary layers with larger convective eddies.We compare the DSD in our simulations to the observations from the flight RF 11 from the VERDI campaign on which we based the simulations' setup.The measured DSD matches almost perfectly our measurements for droplets around the mean size at the last time step of the simulations.The perfect agreement is to some degree by chance, but it nevertheless strongly suggests that radiative cooling is an important factor for the broadening of the DSD for the smaller sizes.The radiative DSD broadening is similar to the broadening due to a wide aerosol size distribution, fluctuations at the activation/condensation process, or inhomogeneous mixing with the free atmosphere (Brenguier and Chaumat, 2000;Beals et al., 2015).At this point we cannot determine which effect is more relevant, and the answer might be different for diverse stratocumulus meteorological conditions.
It has been often noted that a better knowledge of turbulence is necessary to solve the size-gap problem in shallowclouds rain formation (Shaw, 2003;Devenish et al., 2012).In cumulus clouds small-scale turbulence (∼ 10 −3 -10 0 m) can help to bridge the size gap, as it roughly doubles the collision rates when compared to the laminar case (Grabowski and Wang, 2013;Onishi et al., 2015).However, the same smallscale turbulence amplification of collision rates is much weaker (∼ 1 %) for the droplet radius (∼ 10 µm) and turbulence dissipation rates (∼ 10 −3 m 2 s −3 ) that are typical for stratocumuli.This paper suggests that middle-scale turbulence (∼ 10 0 -10 2 m) can help to bridge the size gap in stratocumulus clouds by creating long-resident droplets at the cloud top.was varied by a factor of 4, Re 0 = 200, 400, and 800, which correspond to Kolmogorov length scales η sim = 38, 23, and 14 cm, respectively.We observe that changing the viscosity does not appreciably alter the cloud convective movements.The mean turbulence dissipation rate varies by ∼ 3 % between the different simulations, which is comparable to the variations introduced by the lack of statistical convergence.On the other hand, the integrated evaporative cooling decreases by ∼ 10 % when changing the viscosity from Re 0 = 200 to Re 0 = 800, which suggest that Re 0 = 200 is probably too low to study entrainment.The reason why changes in evaporative cooling do not significantly alter the convective movements is that the flow is mostly radiatively driven (the integrated evaporative-cooling buoyancy source is ∼ 70 % weaker than the integrated radiative-cooling buoyancy source).Since in this study we are mainly interested in the convective dynamics, we decided to use Re 0 = 200 for our reference case.This choice allows us to attain larger z * and longer simulation times at a reduced computational cost.
In Fig. B1 we compare the DSD of the Lagrangian droplets in a region 2λ below cloud top, 540 s after the start of the simulations.We find a strong similarity of all DSDs, indicating a small role of viscous effects.Small differences can be seen at the right side of the distribution, where the number of larger droplets increases with decreasing viscosity.This observation confirms that the cloud-top flow is turbulent.If the flow were laminar, a non-negligible fraction of droplets would stick to the cloud top and continuously grow due to radiative cooling.The number of sticky droplets would increase with viscosity, which is the opposite tendency of what was observed in our simulations.The absence of cloud-top laminar flow is expected, given that the Reynolds number based on the Taylor microscale is high enough at cloud top (Re λ ∼ 150).
If significant, the small increase of the large droplets' size with reducing viscosity can be attributed to different causes.One possible cause is that the radiative-forcing maximum decreases by 5 % when reducing the viscosity from Re 0 = 800 to Re 0 = 200, even when the integrated value of the radiative forcing does not change.As a consequence, the maximum possible growth is slightly higher in the lower viscosity simulations, which could lead to larger droplets.The reduction of the maximum is a consequence of a smoother Eulerian liquid water for the higher-viscosity simulations.Another possible cause is that small-scale turbulent fluctuations, which are stronger for the lower viscosities, broaden the DSD, thus leading to a larger number of large droplets.Even when we cannot quantify these different causes, we conclude that a too-high viscosity probably reduces the right tail of the DSD when the cloud top is fully turbulent, and therefore our results provide a lower limit for the DSD broadening.

Figure 1 .
Figure1.The left side represents a cloudy moist parcel and an unsaturated parcel before mixing.The right side represents the mixing processes in the Eulerian and Lagrangian formulations.In the Eulerian formulation (top) there is diffusion of liquid from the saturated into the unsaturated parcel.In the Lagrangian formulation (bottom) there is no diffusion, and the initially saturated parcel remains free of liquid.As a consequence, liquid water evaporates in both parcels in the Eulerian formulation, but only in the initially saturated parcel in the Lagrangian formulation.Consequently, the evaporation of in the Eulerian formulation is consistently higher.
Figure 2 shows these trajectories from two different viewing angles.The upper blue box in both plots represents the cloudtop region where radiative cooling is active, with a vertical extent of 2λ.The lower blue plane represents the cloud base.

Figure 2 .
Figure 2. Trajectories of Lagrangian cloud droplets from two different viewing points.The droplets have been tracked during 1.2 t * .The upper blue box in each figure represents the cloud-top region where the radiative cooling is strongest, and has a vertical extent of 2λ.The lower blue plane represent cloud base, which is placed 20λ = 300 m below cloud top.A strong stratification limits the flow at cloud base.The trajectories are sorted into three categories (green, red, and orange) that are described in the text.

Figure 4 .
Figure 4. Visualization of the stratocumulus cloud top at t = 11.1 t * .The white contours represent the liquid field from the Eulerian calculations.Red and green colors represent average residence times when extrapolating to the Eulerian grid.Only values larger than 0.9 t * convective-eddy times are shown, thus signaling longresident droplets that have escaped the convective movements.The color scale code goes from red to green, representing mean average residence times from 0.9 t * to 2 t * .

Figure 5 .
Figure 5. Evolution of the cloud-top DSD.The droplets in the simulations grow only by condensation due to radiative cooling.The points correspond to measurements from the VERDI campaign.

Figure B1 .
FigureB1.DSD dependence on viscosity (quantified by the reference Reynolds number Re 0 ) in a cloud-top mixing layer.The DSDs were measured 540 s after the start of the simulations.The overall convergence is good, pointing to a small role of the viscous effects and of the unresolved small scales.