Anthropogenic Air Pollution Delays Marine Stratocumulus Breakup to Open Cells

Marine stratocumulus cloud (Sc) decks with high cloud fraction typically breakup when sufficient drizzle forms. Cloud breakup leads to a lower cloud radiative effect due to the lower cloud amount. Here we use realistic Lagrangian large eddy simulations along a 3‐day trajectory, evaluated with satellite observations, to investigate the timing of Sc breakup in response to aerosol conditions. We show that the timing of the breakup is strongly modulated by the diurnal cycle and large‐scale meteorology but varies systematically with the initial aerosol concentration: the more polluted the clouds, the later the breakup. This indicates that the cloud radiative effect via cloud cover adjustments is not saturated, in contrast to the effect of aerosol on cloud albedo at fixed cloudiness, which weakens with increasing aerosol levels. The results also show that the cloud radiative impact of anthropogenic aerosol is strongest far from its origin over land.


Introduction
The effect of atmospheric aerosol particles on cloud albedo and coverage is one of the largest sources of uncertainty in climate projections (Boucher et al., 2013). For hypothetical, fixed cloud amounts, the aerosol effect on cloud albedo becomes smaller when the reference background is more polluted (Twomey, 1977). It has therefore been proposed that adding aerosol to present day, already high aerosol levels would result in small changes to cloud albedo (Carslaw et al., 2013;Stevens, 2013). However, the response of changes in cloud amount-broadly speaking, cloud fraction (CF) and liquid water path (LWP; referred to as rapid adjustments; Boucher et al., 2013)-does not necessarily exhibit a saturated response to aerosol perturbations. Marine stratocumulus clouds (Sc) are highly sensitive to such adjustments as typical cloud cover varies between being either fully overcast (closed cells) or broken (open cells, with a typical CF≤65 %) (Atkinson & Zhang, 1996;Goren & Rosenfeld, 2014;Muhlbauer et al., 2014). The large change in CF upon transition between the cloud cover regimes exerts a large cloud radiative effect (CRE) (Goren & Rosenfeld, 2014), with a relatively small contribution from the cloud albedo changes at fixed cloud amount (Chen et al., 2014;Goren & Rosenfeld, 2014). (Atkinson & Zhang, 1996;Muhlbauer et al., 2014). Observational (Sharon et al., 2006;Stevens et al., 2005;Wood & Hartmann, 2006;Wood et al., 2011) and modeling (Berner et al., 2013;Kazil et al., 2014;Savic-Jovcic & Stevens, 2008;Xue et al., 2008;Wang & Feingold, 2009;Yamaguchi et al., 2017) studies have shown that the onset of precipitation is the primary factor in the chain process that leads to the breakup of the fully overcast closed cells: Once in the precipitating state, the Sc system transitions from a primarily top-down driven closed-cell circulation to a bottom-up driven cumulus cloud state (Berner et al., 2011;Feingold et al., 2010;Wang et al., 2010). As the precipitating cumulus cloud elements grow vertically and become more dominant, the overcast layer dissipates and is rapidly replaced by precipitating open cells. Thus, the cloud breakup is the last stage of the transition process. Aerosol, among other factors such as sea surface temperature, boundary layer deepening, entrainment, and free tropospheric (FT) humidity along the cloud trajectory, can promote or inhibit precipitation formation and therefore also the timing of the Sc breakup. Untangling the role of aerosol and large-scale conditions on cloud adjustments, as well as investigating their covariability, is the focus of recent studies (Gryspeerdt et al., 2016;Mülmenstädt & Feingold, 2018).
The prevalence of the Sc decks downwind of the major continents expose them to the influence of continental and anthropogenic aerosol (Adebiyi et al., 2015;Adebiyi & Zuidema, 2018;Allen et al., 2011;George et al., 2013;Goren & Rosenfeld, 2015;Schwartz et al., 2002;Saide et al., 2012;Yang et al., 2011). Using a chemistry and aerosol transport model, colocated with satellite observations, Goren & Rosenfeld (2015) (hereafter GR15) related a long-lived extensive overcast Sc deck to air pollution originating in Western Europe. They showed that the complete breakup of the overcast closed-cell cloud deck occurred 3 days after its formation and only after the cloud drop concentration (N d ) was reduced sufficiently to allow precipitation. It was hypothesized that the higher level of anthropogenic aerosol inhibited precipitation over these 3 days, thereby delaying the closed-cell breakup and inducing a large local radiative forcing, mainly due to the larger CF. Had the air mass not been influenced by anthropogenic air pollution, would the overcast closed cell have broken up sooner?
In the current work, we use realistic Lagrangian large eddy simulation (LES), evaluated by observations, to address the above question and to explore the underlying processes in closed-to open-cell transitions and their radiative impact.

Simulation
We developed a new methodology to simulate boundary layer clouds. The approach, realistic Lagrangian LESs, allows for enhanced realism of simulated cloud properties and cloud state evolution relative to observations over periods of days. We employ the System for Atmospheric Modeling (Khairoutdinov & Randall, 2003) as a LES model. We operate the model so that the simulation domain follows a given air mass in the boundary layer, along its Lagrangian trajectory. The model simulates the evolution of the cloud deck embedded in the airmass. Trajectories are determined with the HYSPLIT model (Hybrid Single Particle Lagrangian Integrated Trajectory Model Stein et al., 2015;Rolph et al., 2017) from the fifth generation of the European Centre for Medium-Range Weather Forecasts atmospheric reanalyses (ERA5) reanalysis data (Hersbach, 2018). The trajectories are located 500 m above sea level, at some distance from shear effects near the surface and the inversion. The ERA5 reanalysis provides the mean state of the atmosphere as well as sea surface temperature along the trajectory. We use Newtonian relaxation to nudge mean FT temperature and water vapor in the LESs towards the ERA5 reanalysis, at altitudes 300 m above the simulated inversion, with a nudging time constant of 1,800 s. The mean horizontal wind speed is nudged towards ERA5 at all levels.
The horizontal (vertical) domain size is 76.8 × 76.8 km (5,000 m). A sufficiently large horizontal domain size allows capturing mesoscale organization of the open-  and closed-cell SC cloud state . The horizontal grid spacing is 200 m. The vertical grid spacing for the simulation presented here is 20 m between 0 and 890 m, 10 m between 890 and 1,890 m, and coarsens from 10 m by 10%̇per level between 2,100 and 5,000 m. For additional description of the model setup we refer to section 1 in the supporting information (Iacono et al., 2008;Mlawer et al., 1997;Monin & Obukhov, 1954;Yamaguchi et al., 2011Yamaguchi et al., , 2017. We process the simulated thermodynamic and microphysical quantities for evaluation with satellite data. Cloud water path, rain water path (RWP), and LWP are calculated by vertically integrating cloud water, rain water, and liquid water, respectively, over locations with the associated optical depth ≥3. Cloud optical depth, rain optical depth, and liquid water optical depth are calculated by sampling locations with the respective optical depth ≥3. The satellite optical depth is sampled in a similar fashion. The hydrometeor effective radius is calculated as the ratio of the third to second moments of the cloud and rain drop size distribution. The third and second moments are integrated over locations with cloud and rain optical thickness ≥3, and over one optical depth below cloud top, before their ratio is taken. N d at cloud top is calculated as the average over locations with cloud and rain optical thickness ≥3 and over one optical depth below cloud top. The in-cloud N d is calculated as the average over locations with a cloud water content ≥0.01 g kg −1 .

Satellite Observations
Observed cloud properties were retrieved from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) instrument on board the geostationary Meteosat Second Generation. SEVIRI provides observations with high temporal resolution of 15 min with spatial resolution of 3.5×5.5 km 2 in the northeast Atlantic. We used SEVIRI-based cloud products (CF, cloud optical depth c , drop effective radius r e , and LWP) generated with the Meteosat Second Generation-Cloud Physical Properties algorithm developed at the Royal Netherlands Meteorological Institute available via http://data.knmi.nl (Roebeling et al., 2006). Moderate resolution imaging spectroradiometer (MODIS) true-color images were provided by the National Aeronautics and Space Administration Earth Observing System Data and Information System available at https:// worldview.earthdata.nasa.gov/. In order to reduce uncertainty in cloud property retrieval, pixels having c and r e smaller than 3 are excluded from the analysis as these retrievals are less reliable (Sourdeval et al., 2015). All cloud retrievals were obtained from single-layer liquid clouds. Cloud droplet number concentration is retrieved using the observed r e and c (Szczodrak et al., 2001) and by assuming an adiabatic fraction of 0.8 and cloud top temperature and pressure of 275 K and 950 hPa, respectively. Potential uncertainties in N d arise from the high sensitivity of N d retrieval to r e (Grosvenor et al., 2018;Szczodrak et al., 2001), especially when r e is relatively small. Observations of N d therefore should be considered rather qualitative. The cloud properties along the trajectory were retrieved from an area covering approximately 100×100 km 2 , to be comparable with model domain size.
The HYSPLIT back trajectories deviate from actual boundary layer air mass trajectories. Because the clouds of interest are located close to the interface between the clean and polluted air mass (see Figure 1 top panel), the deviation can locate the area used for averaging of the satellite data partly over the clean air mass. For this reason we constructed back trajectories visually from SEVIRI imagery by tracking cloud features by eye to include only clouds in the polluted, continental air mass. Figures 1a-1f show the observations for both ERA5/HYSPLIT and the visually-determined trajectory. Figures S4 and S11 compare the ERA5/HYSPLIT back trajectories with the visually determined trajectories. Sampling along the visually determined back trajectory captures the time evolution of the simulated cloud properties better, while the overall trend is similar for both sampling approaches.
The LWP c was calculated as follows: an adiabatic calculation with prescribed cloud base temperature and pressure (5 • C and 950 hPa, respectively; estimated from the simulated trajectory) was used to generate a liquid water content profile with vertical resolution of 10 m and assuming 80 % adiabaticity. Given the simulated N d at every time step, N d t , and for each simulation, a vertical profile of the volume effective radius, r v , was calculated as: , with subscripts i and t indicating the vertical level and LES simulation time step, respectively, and w being the water density. For converting r v to r e we assumed that r e = 1.08r v . Having derived the vertical profile of r e , we determined LWP c to be the LWP at the level where r e = 14 m.

Results
The simulations were run along a 3-day trajectory over the northeast Atlantic Ocean downwind of the coasts of Western Europe, on which closed cells transitioned to open cells (GR15). Figure 1, top panel, shows MODIS true-color images of the Sc for each day of the 3-day event, as well as the trajectory location. Figure 1, lower panel, shows the simulated cloud properties (CF, cloud optical thickness c , effective radius r e , N d , LWP, and surface precipitation) along the trajectory for varying initial aerosol concentration, N a i (NA125 to NA300 in intervals of 25 cm −3 ). Also shown are the satellite-observed cloud properties along the trajectory. The satellite retrievals of cloud properties are limited to daytime only due to the dependence on solar radiation. The time evolution of the cloud properties for the simulations with N a i of 200-300 cm −3 tracks the observations reasonably well (Figure 1a-1e), providing confidence that the underlying cloud dynamics and microphysics are captured by the model. In the following, we refer to the simulations with N a i < 200 cm −3 as representing Sc clouds with lower degrees of pollution, from which assessments of the Sc breakup times, as well as fundamental features of the system, are derived.

Cloud Breakup Time Dependence on Initial Aerosol Concentration
We diagnose the closed-to open-cell breakup time based on the CF. As can be seen in Figure 1a, the cloud breakup occurs later for higher N a i , implying that indeed, under cleaner conditions, the closed cells would breakup sooner. The trajectory follows closed cells that were observed to breakup on Day 3, while its surroundings are already broken (Figure 1; upper panel and GR15). On Day 1 the clouds in the location of the trajectory had elevated N d compared to the surrounding clouds (supporting information Figure S1). Therefore, simulations with lower N a i can represent trajectories of the nearby less polluted clouds, which breakup sooner.
The closed-cell breakup and emerging spatial organization of the open-cell state is shown in the animated sequence of cloud optical depth (Animation S1). The clearing shown in Animation S1 and Figure 1a proceeds rapidly once surface precipitation begins and accelerates for rates exceeding 0.1 mm day −1 (Figures 1f and 2). The formation of surface precipitation occurs simultaneously with r e , reaching about 12-14 m, both in the observations and in the simulations (Figure 1c). An r e of 12-14 m is a good indicator of significant drizzle (Gerber, 1996;Rosenfeld et al., 2012) and has been shown to be associated with transitions from closed to open cells (Goren & Rosenfeld, 2014;Rosenfeld et al., 2006).

Cloud Transition Time Scales
Prior to cloud breakup, N d exhibits two regimes of removal ( Figure 1d): (1) a slow decay of N d followed by (2) a rapid decrease that leads to precipitating open cells and very clean conditions. The existence of the two removal regimes can be explained by the nonlinearity of drizzle formation processes (Gunn & Phillips, 1957); as long as cloud droplets are small, collision-coalescence processes are inefficient and N d decreases slowly. As the droplets grow and form rain drops (accompanied by the simultaneous rapid increase in r e ; Figure 1c), collision-coalescence accelerates in a positive feedback that accelerates reduction in N d . This rainout-cleansing process leads eventually to a very clean state with N d on the order of 10 per cm −3 (Figure 1d). As the boundary layer deepens ( Figure S2), entrainment of cleaner, FT air dilutes the aerosol concentration, N a , and can also serve as a significant factor in lowering N d (Berner et al., 2013;Mechem et al., 2006). Entrainment rates and FT aerosol concentrations can therefore control to some extent the rate of the slow decay regime of N d , thereby affecting the cloud breakup time (Mechem et al., 2006;Wood et al., 2012;. The maximum scavenging rates, defined as d 2 N d ∕dt 2 =0, commensurate with a rapid buildup of rain water mass, are marked in Figure 1d for each of the simulations (black line) and replotted against N a i in Figure 2. Also shown in Figure 2 are the times of selected key physical processes leading to cloud breakup. These illustrate the chronological order in which the cloud transition occurs: (1) initial formation of rain drops; (2) rapid scavenging of N d ; (3) surface precipitation, which is associated with the beginning of the cloud breakup (Figure 1a and 1f); and (4) formation of open cells, defined here as CF≤65 % (a typical CF for open cells Goren & Rosenfeld, 2014;Muhlbauer et al., 2014). This emphasizes that the breakup becomes possible only when clouds start to precipitate heavily enough (surface precipitation rates of 0.1 mm day −1 ; Figure 2), and N a become sufficiently low, in agreement with previous studies (Berner et al., 2013;Stevens et al., 2005;Wang & Feingold, 2009). Figure 2 introduce the delay in transition onset and transition duration time scales, respectively, which together determine the cloud breakup time. The delay in transition onset is the additional time it takes to begin buildup of rain water, here defined as RWP c >1 g m −2 (blue symbols in Figure 2), from one simulation to its adjacent one with larger N a i . The transition duration is defined here as the time between RWP c and the formation of surface precipitation exceeding 0.1 mm day −1 (blue and orange symbols, respectively, in Figure 2). A systematic delay in transition onset, and the commensurate cloud break-up, can be seen with higher N a i . This can be explained based on the dependence of drizzle formation on both N d and LWP (Berner et al., 2013;Comstock et al., 2004;Mechem et al., 2012;Pawlowska, 2003;Terai et al., 2012;.

The green and black vertical lines in
with q l being the adiabatic cloud water content at cloud top. Using with LWP adi the adiabatic LWP (see supporting information) and combining (1) and (2) gives By assuming a critical r e of 12-14 m for initiation of drizzle (Gerber, 1996), r e c , one can assign an equivalent theoretical critical LWP for breakup, LWP c .
Equation (4) expresses a quadratic relationship between LWP c and N d . Thus for a given LWP, polluted clouds would require a higher LWP c and form precipitation later than cleaner clouds Konwar et al., 2012), in accordance with our simulations. Figure 3a illustrates the microphysical processes (manifested by changes in N d ) and large-scale meteorology (manifested by changes in LWP) that together determine the time when the actual LWP = LWP c (equation 4). Theoretical curves of LWP c were calculated based on the simulated N d at every time step (colored lines in 3a). They represent processes that are responsible for reducing N d (e.g., slow scavenging rate and dilution with FT air), which slowly lower LWP c with time. Also shown is a time series of the LWP for NA300 (blue line in 3a), which is the last to breakup and therefore provides an envelope LWP that is similar in all simulations up to the onset of the transition processes (see also Figure 1e). The LWP associated with NA300 thus represents large-scale meteorology and diurnal cycle. When the LWP c and the LWP associated with NA300 intercept, surface precipitation forms (roughly when the simulated r e is about 12 m) and the breakup takes place. Equation (4) thus successfully identifies the cloud breakup time based on the meteorological and aerosol state of the system and may serve as a diagnostic for Sc breakup to open cells in global climate models (GCMs).

Meteorology Versus Aerosols in Controlling Cloud Breakup Time
Evidence for the strong role of meteorology is the rapid increase in LWP on Day 3 (Figure 1e), which is associated with a FT moist layer making contact with the deepening boundary layer (Figures S3 and S2). This allows the LWP of NA250-NA300 to reach LWP c faster and to form precipitation that eventually breaks up the clouds. An additional set of simulations was performed along a second trajectory of transitioning closed cells originating nearby and having roughly similar N a on Day 1 but which takes a different path and breaks up a day sooner. A similar systematic response of breakup time versus N a i was found but with shorter delay times due to a more dramatic water vapor intrusion from the FT. Analysis of this case can be found in the supporting information (section S2 and Figures S7-S11). Figure 3b shows the NA300 LWP tendency, in which changes in LWP driven by the diurnal cycle are embedded within changes in LWP driven by the large scale meteorology. Whether the LWP increases or decreases, Figure 3. Illustration of the microphysical processes (manifested by changes in N d ) and large-scale meteorology and diurnal cycle (manifested by changes in liquid water path [LWP]) in codetermining the cloud state evolution. (a) Time series of LWP of the NA300 simulation and the theoretical critical LWP for each initial N a i simulation, LWP c (equation 4), for which r e = 14 m, when precipitation often exists and cloud breakup is anticipated. Symbols represent the local time when the simulated r e = 12 m is reached and precipitation forms (each symbol and color represents each case, e.g., red points represent NA300). An r e in the range of 12-14 m is an indicator of significant drizzle (Gerber, 1996;Rosenfeld et al., 2012). (b) The tendency of NA300 LWP (black line) with symbols representing rain water path (RWP) c , the transition onset time (similar convention as in panel a). The quadratic relationship between LWP c and N d (see text) explains the systematic delay in transition onset. Decreases in LWP (i.e., negative tendencies) can lengthen the delay in transition onset and transition duration (e.g., see illustration for NA150), while rapid increases (i.e., positive LWP tendencies) make them shorter (e.g., see illustration for NA275). Gray shading in a and b indicates nighttime. (c) Correlation between delay in transition onset and LWP tendency (R 2 =0.69). (d) Correlation between transition duration and LWP tendency for a period of 1 hr after RWP c is reached (R 2 = 0.76).
hastens or dampens precipitation formation processes, thereby affecting the transition onset time and duration and the resultant breakup time. The symbols in Figure 3b represent the transition onset (RWP = RWP c ) for each simulation. The delay in transition onset can be inferred from the times at which adjacent simulations (in terms of N a i ) first form an RWP c . The transition duration corresponds to the period between RWP c in Figure 3b to the time at which r e reaches 12 m in Figure 3a. The two time scales are illustrated in Figures 3a and 3b for NA150 and NA275, respectively. The delay in transition onset and the transition duration were found to be anticorrelated with the LWP tendency (Figures 3c and 3d, corresponding to R 2 of 0.69 and 0.76, respectively), which greatly modulates the breakup time due to N a i . Nevertheless, even within the variability of LWP, greater N a i systematically delays the cloud breakup (Figure 2, red lines).

Summary
In this study, we have identified causal relationships between initial aerosol concentration N a i , meteorology, and cloud breakup time in LES. We have shown that the delay of cloud breakup due to additional aerosol is modulated by changes in LWP driven by large-scale meteorology and the diurnal cycle. Together, they determine the onset of the transition process, as well as its duration. The breakup time varies systematically with N a i so that the greater N a i , the later the breakup occurs. This means, that for the event studied here, cloud cover adjustments to aerosol are not saturated, in contrast to the cloud albedo response to aerosol at fixed cloud amount, which weakens as the background aerosol levels are increased. This is shown in Figure 4, where the daily mean shortwave CRE in the LES is plotted as a function of N a i for the last 2 days of the simulations. The CRE is nearly saturated on Day 2 for N a i > 200 cm −3 , in line with the logarithmic nature of the cloud albedo response to aerosol (Twomey, 1977). On Day 3 on the other hand, CF responds strongly to additional N a i , driving a commensurate response in the CRE.
Our results imply that the cloud radiative forcing due to cloud cover adjustments to aerosol in polluted Sc decks may well be found over the westernmost parts of the overcast decks, unintuitively far from anthropogenic aerosol sources. Assessing the global significance requires quantifying the frequency of occurrence of interactions between anthropogenic aerosol and Sc susceptible to transition. In order to do so, developing observational methods and strategies is required. GCMs could also be used; however, shallow clouds and their adjustments to aerosol are poorly represented in current models (Ghan et al., 2016;Malavelle et al., 2017;Neubauer et al., 2014). Improving shallow cloud representation in GCMs, and in particular Sc cloud cover adjustments to aerosol, is therefore vital to reducing uncertainty in climate projections.