Methods and Published by Copernicus Publications Data Systems on behalf of the European Geosciences Data Systems Union.

Abstract. A cloud-resolving model (CRM) coupled to a new intermediate-complexity bulk aerosol scheme is used to study aerosol–boundary-layer–cloud–precipitation interactions and the development of pockets of open cells (POCs) in subtropical stratocumulus cloud layers. The aerosol scheme prognoses mass and number concentration of a single lognormal accumulation mode with surface and entrainment sources, evolving subject to processing of activated aerosol and scavenging of dry aerosol by clouds and rain. The CRM with the aerosol scheme is applied to a range of steadily forced cases idealized from a well-observed POC. The long-term system evolution is explored with extended two-dimensional (2-D) simulations of up to 20 days, mostly with diurnally averaged insolation and 24 km wide domains, and one 10 day three-dimensional (3-D) simulation. Both 2-D and 3-D simulations support the Baker–Charlson hypothesis of two distinct aerosol–cloud "regimes" (deep/high-aerosol/non-drizzling and shallow/low-aerosol/drizzling) that persist for days; transitions between these regimes, driven by either precipitation scavenging or aerosol entrainment from the free-troposphere (FT), occur on a timescale of ten hours. The system is analyzed using a two-dimensional phase plane with inversion height and boundary layer average aerosol concentrations as state variables; depending on the specified subsidence rate and availability of FT aerosol, these regimes are either stable equilibria or distinct legs of a slow limit cycle. The same steadily forced modeling framework is applied to the coupled development and evolution of a POC and the surrounding overcast boundary layer in a larger 192 km wide domain. An initial 50% aerosol reduction is applied to half of the model domain. This has little effect until the stratocumulus thickens enough to drizzle, at which time the low-aerosol portion transitions into open-cell convection, forming a POC. Reduced entrainment in the POC induces a negative feedback between the areal fraction covered by the POC and boundary layer depth changes. This stabilizes the system by controlling liquid water path and precipitation sinks of aerosol number in the overcast region, while also preventing boundary layer collapse within the POC, allowing the POC and overcast to coexist indefinitely in a quasi-steady equilibrium.

accumulation mode with surface and entrainment sources, evolving subject to processing of activated aerosol and scavenging of dry aerosol by cloud and rain.
The LES with the aerosol scheme is applied to a range of steadily-forced simulations idealized from a well-observed POC case. The long-term system evolution is explored with extended two-dimensional simulations of up to 20 days, mostly with 10 diurnally-averaged insolation. One three-dimensional two-day simulation confirms the initial development of the corresponding two-dimensional case. With weak mean subsidence, an initially aerosol-rich mixed layer deepens, the capping stratocumulus cloud slowly thickens and increasingly depletes aerosol via precipitation accretion, then the boundary layer transitions within a few hours into an open-cell regime with scattered 15 precipitating cumuli, in which entrainment is much weaker. The inversion slowly collapses for several days until the cumulus clouds are too shallow to efficiently precipitate. Inversion cloud then reforms and radiatively drives renewed entrainment, allowing the boundary layer to deepen and become more aerosol-rich, until the stratocumulus layer thickens enough to undergo another cycle of open-cell formation. If mean sub- 20 sidence is stronger, the stratocumulus never thickens enough to initiate drizzle and settles into a steady state. With lower initial aerosol concentrations, this system quickly transitions into open cells, collapses, and redevelops into a different steady state with a shallow, optically thin cloud layer. In these steady states, interstitial scavenging by cloud droplets is the main sink of aerosol number. The system is described in a re- 25 duced two-dimensional phase plane with inversion height and boundary-layer average aerosol concentrations as the state variables. Simulations with a full diurnal cycle show similar evolutions, except that open-cell formation is phase-locked into the early morning hours.
The same steadily-forced modeling framework is applied to the development and evolution of a POC and the surrounding overcast boundary layer. An initial aerosol perturbation applied to a portion of the model domain leads that portion to transition 5 into open-cell convection, forming a POC. Reduced entrainment in the POC induces a negative feedback between areal fraction covered by the POC and boundary layer depth changes. This stabilizes the system by controlling liquid water path and precipitation sinks of aerosol number in the overcast region, while also preventing boundarylayer collapse within the POC, allowing the POC and overcast to coexist indefinitely in 10 a quasi-steady equilibrium.

Introduction
Marine stratocumulus clouds cover broad swaths of the world ocean, exerting a strong net radiative cooling effect on climate due to their high albedo (Hartmann et al., 1992). The turbulent circulations maintaining marine stratocumulus are on the order of the 15 boundary layer thickness (1 km), far below the spatial scale resolved by global climate models (GCMs). Marine stratocumulus clouds are thin (100-500 m deep), and are usually capped by a sharp, strong temperature inversion. These sharp gradients are not resolved by GCMs and a challenge to parameterize. They are maintained by strong feedbacks between turbulence, cloudiness, radiation, aerosols and precipitation; 20 in a GCM these are parameterized processes involving substantial subgrid variability. These factors combined make marine stratocumulus a key challenge for simulating climate, including two central uncertainties in simulating anthropogenic climate change -boundary cloud feedbacks (Bony and Dufresne, 2005) and cloud-aerosol interaction (Quaas et al., 2009). 25 Cloud microphysical properties are modulated by natural and anthropogenic sources of cloud condensation nuclei (CCN). Their influence on net cloud radiative forcing is 18145 usually discussed in terms of the first and second aerosol indirect effects. The first indirect effect (Twomey, 1977) refers to the increase of cloud albedo for constant liquid water path (LWP) when cloud droplet number concentration (N d ) is increased. The second indirect effect (Albrecht, 1989) describes a broader class of feedbacks relating to the impact of changes in N a on N d and albedo via cloud macrophysical properties 5 such as fractional cloud cover, LWP, and lifetime by means of interactions between microphysical perturbations and cloud dynamics. Cloud microphysics and dynamics also affect the aerosol. For instance, the droplet collision-coalescence that forms precipitation also reduces aerosol number concentration N a , increasing cloud droplet size and further enhancing precipitation in a positive 10 feedback loop. Baker and Charlson (1990) suggested that this feedback process could lead a cloud layer to evolve toward one of two very different aerosol concentrations ("multiple equilibria") for a single set of large-scale forcings. Their simple steady state model was based on a mixed layer capped by a stratocumulus cloud of fixed liquid water content (LWC). They predicted marine boundary layer (MBL) aerosol concentra- 15 tion given a specified surface aerosol source, treating aerosol sinks due to coagulation, collision-coalescence, and drizzle fall-out. Over a commonly-realized range of aerosol source strengths, their model predicted that an MBL with initial N a below a threshold value would evolve over a few days toward a drizzling state with aerosol concentration ∼ 10 cm −3 and small cloud albedo, while an MBL with initial N a above the threshold 20 value would evolve toward a non-precipitating state with very high aerosol concentration ∼ 1000 cm −3 and much larger cloud albedo. Ackerman et al. (1994) approached the same problem using a column model including a cloud layer which responded to the aerosol, parameterized turbulent vertical transports, sophisticated microphysics and radiation. Unlike Baker and Charlson (1990), they found no evidence for bistability, 25 obtaining instead a smooth increase in equilibrium aerosol concentration with source strength, because the total aerosol sink term in their model was an increasing function of N a . The relevance of aerosol bistability to real stratocumulus cloud regimes remains an open question, to be further discussed in this paper.
An observable manifestation of a positive aerosol-cloud feedback and perhaps bistability is the formation of pockets of open cells, or POCs (Stevens et al., 2005), lowalbedo regions of cumuliform, open-cellular convection embedded within a sheet of high-albedo closed-cell stratocumulus convection. Early aircraft observations of POCs during the EPIC (Bretherton et al., 2004) and DYCOMS-II (Stevens et al., 2003) field 5 campaigns (Comstock et al., 2005(Comstock et al., , 2007van Zanten and Stevens, 2005) documented the contrasts between POCs and surrounding stratocumulus regions. Subsequent observations indicated that the aerosol concentration within POCs is also highly depleted, especially in an "ultraclean layer" just below the inversion (Petters et al., 2006;Wood et al., 2011a). Satellite observations showed POCs form preferentially in the early 10 morning, when stratocumulus cloud is typically thickest and drizzliest, and persist once formed (Wood et al., 2008). The VOCALS Regional Experiment (REx; Wood et al., 2011b) in September-October 2008, targeting the massive and persistent stratocumulus deck off the west coast of Chile, was designed in part to systematically document POC structure. VOCALS Research Flight 06 comprehensively sampled a mature POC 15 one day after its formation, providing an invaluable dataset for modeling and analysis (Wood et al., 2011a). Four other POCs were also sampled by the NCAR/NSF C-130 during VOCALS-REx (Wood et al., 2011b). However, observations alone cannot definitively explain how POCs are formed and maintained -this requires a model, and the challenges of resolution and parameterization uncertainties make GCMs nonideal for 20 this task.
Our study focuses on large-eddy simulation (LES) of the MBL as a means to understand the complexities of stratocumulus-aerosol-precipitation interaction and their role in POC formation and maintenance. LES of MBL turbulence and clouds explicitly simulate the larger, most energetic eddies that dominate turbulent fluxes. Like GCMs, they 25 require parameterizations of moist thermodynamics, cloud microphysics, radiation, and subgrid turbulence, but they do not require the elaborate formulations of subgrid variability that make these parameterizations complex and uncertain in GCMs. Early LES studies of aerosol impacts on cloud dynamics included Ackerman et al. (2003) and 18147 Ackerman et al. (2004), which examined the response of a small-domain simulation of marine stratocumulus to specified changes in cloud droplet concentration. Their simulations found a low-aerosol, heavy drizzle regime in which the cloud cover increased and thickened with more aerosol and a higher-aerosol, nearly nondrizzling regime in which the cloud thinned with more aerosol. 5 Recent LES modeling studies in domains with sizes on the order of tens of km have examined the influence of aerosol changes on the ubiquitous mesoscale cellular variability in stratocumulus cloud regimes (Wood and Hartmann, 2006). Some studies focused on cloud response to specified N d or aerosol. Savic-Jovcic and  and Xue et al. (2008) found increased open-cellular organization with decreas-10 ing cloud droplet or aerosol concentration, though the simulations didn't exhibit as marked a decrease in cloud cover as found in POCs. Berner et al. (2011) simulated the VOCALS-REx RF06 POC case, imposing observationally based N d differences between the POC and the surrounding overcast region, and reasonably reproduced the observed contrasts in cloud, precipitation and boundary layer structure across the 15 POC boundary. Their simulations also revealed that the overcast region entrained much more strongly than inside the POC, yet the mean inversion height across the domain remained essentially level, with mesoscale circulations compensating the reduced entrainment in the POC.
A full understanding of the role of aerosol-cloud interactions in the climate sys-20 tem requires simulation of the feedback of cloud processes on aerosols. For this purpose, several models have added a simple CCN budget including a number sink due to precipitation-related processes. Mechem and Kogan (2003) used this approach in a mesoscale model (not an LES) with horizontal resolution of 2 km. They simulated transitions from aerosol-rich stratocumulus layers to aerosol-poor, precipitating layers with partial cloud cover. Mechem et al. (2006) also included a prescribed surface CCN source and varying free troposphere (FT) CCN concentration. They identified entrainment of FT CCN as an important buffering mechanism for MBL CCN, because the surface CCN source is often too weak to balance the collision-coalescence sink.
Subsequent studies have used LES with simple CCN budgets in mesoscale-size domains. Wang and Feingold (2009a) used 300 m (barely eddy-resolving) resolution to look at the evolution of a stratocumulus layer with three initial CCN concentrations, finding sensitivities similar to Mechem and Kogan (2003). Wang and Feingold (2009b) used a 60 km × 180 km domain with an initial linear gradient in CCN, showing that the 5 development of open-cell organization smoothly increased as initial CCN decreased. Wang et al. (2010) showed that POCs could be rapidly triggered by reducing the initial aerosol concentration in a mesoscale region within a solid stratocumulus layer, and that the resulting aerosol perturbations and POC structure persist for the 8 h length of their simulation. They also showed that depending on the initial aerosol concentra-10 tion, a stratocumulus-capped boundary layer with the same initial cloud characteristics could either quickly transition into open cells that further reduced aerosol concentration, or persist in a high-aerosol, nearly nonprecipitating state for the length of a 36 h simulation. Several LES studies have also included more complete models of aerosol processes 15 in cloud-topped boundary layers. Feingold et al. (1996) coupled bin aerosol microphysics to a 2-D LES and analyzed the role of aqueous chemistry. Feingold and Kreidenweis (2002) explored the effects of different initial aerosol distributions and aqueous chemistry on cloud dynamics over periods up to eight hours, and coined the term "run-away precipitation sink" for the increasingly efficient removal of aerosol via pre-20 cipitation processes at lower values of N d . Ivanova and Leighton (2008) implemented a three mode, two moment bulk aerosol scheme within a mesoscale model based on the approximation that aerosol mass within rain or cloud water droplets is proportional to their water mass, but coagulates into one aerosol particle per droplet, as suggested by Flossmann et al. (1985). This approach allows for the inclusion of relatively complete 25 aerosol microphysics with a minimum of additional advected scalars. Kazil et al. (2011) coupled detailed aerosol and gas chemistry into an LES to simulate open-cell convection within the VOCALS RF06 POC, obtaining a realistic simulation of the vertical distribution of aerosol and of the ultraclean layer, and simulating an episode spontaneous 18149 nucleation of new aerosols within the ultraclean layer. Arguably, this is the most realistic LES depiction of coupled aerosol-stratocumulus cloud-precipitation interactions to date.
In this paper, we build on prior aerosol-cloud-precipitation studies to analyze the multiday evolution and equilibrium states of the coupled stratocumulus cloud-aerosol 5 system subject to a range of constant forcings and initial conditions. For this purpose, we couple a new intermediate-complexity aerosol model to our LES. We recognize that in reality, MBL air is always advecting over a changing SST and subject to changing synoptic forcings and free-tropospheric conditions. However, it is still reasonable to ask whether one can define preferred "regimes" through which the aerosol-cloud system tends to evolve, and if so, whether the system can be expected to evolve smoothly between regimes or whether there are circumstances favoring rapid transitions, or in which the system evolution is sensitive to small changes in the external forcings or initial conditions, and lastly, how such results bear on the aerosol bistability controversy.
Our aerosol model is slightly more complete than schemes which only predict CCN 15 number. It prognoses mass as well as number for a single log-normal accumulation mode, a simplified form of an approach used by Ivanova and Leighton (2008) that is easily coupled into the warm rain portion of the Morrison et al. (2005) double-moment microphysics package used in our LES. Feingold (2003) showed that aerosol geometric mean radius affects the activated fraction of aerosol and thereby aerosol indirect ef- 20 fects. In addition to the processes usually included in simple CCN-predicting schemes (parameterized surface source, entrainment from the FT, collision-coalescence, and precipitation fall out), we include interstitial scavenging by cloud and rain, which we find to be an important number sink in less heavily precipitating stratocumulus. Unlike a full-complexity aerosol scheme like that of Kazil et al. (2011), our scheme is still com- 25 putationally efficient enough for the extended runs used for the present study, at least when run in two spatial dimensions (2-D).
Our study is organized as follows: we develop a single-mode, double-moment bulk aerosol scheme inspired by Ivanova and Leighton (2008). Section 2 describes the LES and the aerosol scheme. Section 3 describes how our mostly 2-D simulations are initialized and forced, and Sect. 4 gives an overview of the different simulations. In Sect. 5, we discuss a control simulation which undergoes a transition from closed cell stratocumulus to more cumuliform drizzle cells that is qualitatively consistent with observations from the VOCALS field campaign, as well as prior simulations of Mechem and Kogan 5 (2003) and Wang et al. (2010). We also show that a identically forced two-day 3-D run has a qualitatively similar aerosol-cloud evolution as the 2-D control run. Section 6 demonstrates that a small increase in subsidence forcing leads to a stable, overcast equilibrium state instead of the control run evolution, and sensitivity runs test the role of the availability of FT aerosol and the diurnal cycle of insolation. In Sect. 7, we examine 10 a set of identically-forced simulations initialized with different MBL aerosol concentrations with a phase plane analysis similar to Bretherton et al. (2010). We use this phase plane to discuss the concept of regimes and their transitions, and we show that there can be multiple equilibria as a function of the initial aerosol state. Finally, in Sect. 8, we use runs initialized with an initial gradient of aerosol concentration to examine the 15 formation of a POC and the subsequent evolution and long-term stability of the coupled system comprised of two cloud regimes -a POC and the surrounding overcast stratocumulus. We find that this system can achieve long-term stability due to entrainment feedbacks. 20 The simulations in this paper are performed using version 6.9 of the System for Atmospheric Modeling (SAM; Khairoutdinov and Randall, 2003). SAM uses an anelastic dynamical core. An f -plane approximation with Coriolis force appropriate for the specified latitude is used. In our simulations, the clouds are liquid phase only, and SAM partitions water mass into vapor mixing ratio q v , cloud water mixing ratio q c (drops 25 smaller than 20 micron radius), and rain water mixing ratio q r (drops larger than 20 micron radius). These mixing ratios are separately advected in our version, along with 18151 liquid-ice static energy s li = c p T + gz − Lq l (neglecting ice). Here c p is the isobaric heat capacity of air, g is gravity, z is height, L is the latent heat of vaporization, and the liquid water mixing ratio q l is the sum of q c and q r . Microphysical tendencies are computed using two-moment Morrison microphysics (Morrison and Grabowski, 2008;Morrison et al., 2005), which requires that we also prognose a cloud water number concentra-5 tion N d and a rain number concentration N r ; saturation adjustment is used to repartition between q v and q c each time step. We have modified SAM to advect scalars using the selective piecewise-parabolic method of Blossey and Durran (2008), which is less numerically diffusive than SAM's default advection scheme. The 1.5 order TKE scheme of Deardorff (1980) is used as the sub-grid turbulence closure (SGS). The SGS length scale is chosen as the vertical grid spacing, as this inhibits unrealistically large mixing on the highly anisotropic grids (large horizontal relative to vertical spacing) needed to efficiently resolve both the inversion and mesoscale structure in stratocumulus LESs. Surface fluxes are computed in each column from Monin-Obukhov theory. Radiation is calculated every 15 s using the Rapid Radiative Transfer Model (RRTM; Mlawer et al., 15 1997); the solar zenith angle is set for the VOCALS RF06 POC location at 17.5 • S 79.5 • W. Drizzle has been included in the radiation calculation by specifying a combined cloud-drizzle water effective radius within each grid cell, computed to give the sum of the optical depths due separately to cloud and drizzle drops in that grid layer. 20 We have implemented a new, computationally efficient, single mode, double moment aerosol scheme for warm clouds which tightly couples with Morrison microphysics, including surface fluxes, entrainment, collision-coalescence, evaporation, and scavenging of interstitial aerosol. The aerosol is described by a single log-normal distribution with prognostic mass and number concentration. All the aerosol is assumed to be equally hygroscopic, with the properties of ammonium sulfate. Where condensate (either cloud or rain water) is present, we assume that the "wet" fraction of the distribution that resides in the condensate particles will correspond to the largest (and hence most easily activated) aerosol particles. The prognosed aerosol should loosely be thought of as corresponding to the accumulation mode. Since we do not represent the Aitken mode, some processes such as spontaneous nucleation of new aerosol particles, or the growth of a population of small 5 aerosol particles via deposition of gasses such as SO 2 , are not represented. However, the single-mode approximation does capture a variety of important aerosol-cloud interaction processes, as we will see, and its simplicity and efficiency make it especially attractive for more idealized simulations.

Single-mode aerosol scheme
A schematic representation of the various source, sink, and transfer terms for mass 10 and number, considered in the MBL mean, is shown in Fig. 1. We follow the approximation of Ivanova and Leighton (2008) that activated aerosol mass q aw is affected by moist processes in proportion to the cloud water mass q c , and similarly for rain aerosol mass q ar vs. rain water mass q r : where | mp denotes tendencies due to microphysics. Including the effects of aerosol in this way requires prognosing and advecting only four additional scalars, namely dry aerosol number concentration (N ad ) and mass mixing ratio (q ad ) , activated aerosol 20 mass mixing ratio in cloud (q aw ), and aerosol mass mixing ratio in rain (q ar ). We assume that the aerosol is fully soluble within condensate and produces a single particle upon evaporation of a cloud or rain drop. At any time, the total aerosol concentration and mass that define the log-normal distribution are the sum of components from the dry 18153 (unactivated) aerosol, cloud droplets and rain drops: q a = q ad + q aw + q ar .
Ivanova and Leighton included three modes in their scheme: unprocessed aerosol, 5 cloud-processed aerosol, and a coarse mode resulting from evaporated precipitation. By simplifying their approach to a single accumulation mode, we keep the number of auxiliary scalar fields to a minimum. We carry the aerosol mass mixing ratio in addition to number concentration for compatibility with the droplet activation parameterization of Abdul-Razzak and Ghan (2000) used by the Morrison microphysics scheme (which requires these aerosol parameters), and compatibility with a newly developed scavenging parameterization for interstitial aerosol. In the interstitial scheme, described further in the Appendix, scavenging tendencies for q ad and N ad are computed using the cloud and raindrop size spectra from the Morrison scheme, together with approximate collection kernels for convective Brownian diffusion, thermophoresis, diffusiophoresis, turbulent coagulation, interception, and impaction. While interstitial scavenging has been ignored in several other studies, cloud droplets have a reasonably high collection efficiency for unactivated accumulation-mode aerosol (Zhang et al., 2004). Surface fluxes are computed using a modified version of the wind-speed dependent sea salt parameterization of Clarke et al. (2006), where we have refit the size-resolved 20 fluxes with a single, log-normal accumulation mode. As we are concerned mainly with particles at sizes where they will be viable CCN, we choose to center the source distribution about the geometric radius of 0.13 µm. To include the number and mass contributions from the smaller and most numerous portion of the coarse mode, as well as the smallest end of the accumulation mode that may be active CCN at higher super- 25 saturations, we choose a width parameter σ g = 2. We then choose an aerosol number source that gives 50 % of the total integrated number flux in the distribution given by Clarke et al. (2006); the remaining particles being assessed to be too small to act as CCN.
The mass flux is only 0.5 % of that given by Clarke et al. (2006), for which the main 5 mass flux is in large coarse-mode aerosols that again lie outside the range of sizes across which our fit is optimized. Figure 2 plots the size-resolved mass and number fluxes at a windspeed of 9 m s −1 for the Clarke et al. (2006) parameterization against our unimodal approximation. Aerosol in the free troposphere can be brought into the cloud layer through model-10 simulated entrainment. Observations suggest that over remote parts of the oceans, the FT generally has a substantial number concentration of aerosol particles that can act as CCN (e. g. Clarke, 1993;Clarke et al., 1996;Allen et al., 2011), but these particles usually have diameters significantly smaller than 0.1 microns, and must grow via coagulation and gas-phase condensation in the boundary layer before they activate into 15 cloud droplets. Because we can not accurately represent this process, we short-circuit it by specifying a free-tropospheric aerosol number concentration comparable to measured CCN concentrations, but by choosing these particles to already have a mean size of 0.1 micron and distribution width σ g = 2, similar to what we assume for the surface aerosol source. Functionally, this is like assuming entrained aerosols instantaneously 20 grow to this mean size by coagulation and gas-phase condensation when they enter the boundary layer, a process which in reality may take hours to days (Clarke, 1993). This assumption also applies to the contribution of number from the surface source at smaller sizes, which we consider as viable CCN when refitting the Clarke et al. (2006) parameterization to our single accumulation mode. Inclusion of nucleation from the gas 25 phase and other chemistry is outside the scope of our idealized formulation.

18155
The full system of equations governing the bulk aerosol moments in each grid cell are: Here subscript Srf denotes a surface flux, Act is CCN activation, ScvCld is interstitial scavenging by cloud droplets, ScvRn is interstitial scavenging by rain droplets, Evap is evaporation, Auto is autoconversion, Accr is accretion, SlfC is self-collection of rain, and NMT are non-microphysical terms (advection, large scale subsidence, and sub-15 grid turbulent mixing). The equations for N d and N r are identical to the standard implementations within version 3 of the Morrison microphysics implemented in SAM v6.9, though the aerosol mass and number used within the Abdul-Razzak and Ghan activation parameterization are now the local values of q a = q ad + q aw and N ad + N d , respectively. The equation for 20 N ad also includes cloud and rain evaporation terms from the Morrison schemes, as well as contributions from the surface flux, interstitial scavenging, and advection/turbulence. The microphysical mass sources and sinks are derived from the corresponding Morrison cloud and rain mass sources and sinks using the Ivanova-Leighton approximations (1) and (2). One modification was made to the Morrison scheme. The default scheme assumes no loss of cloud droplet number due to evaporation until all water is removed 5 from a grid box (homogeneous mixing). This assumption has been replaced with the heterogeneous mixing assumption that cloud number is evaporated proportionately to cloud water mass, which theory and field measurements suggest may be more appropriate for stratocumulus clouds (Baker et al., 1984;Burnet and Brenguier, 2007).
The discretized system preserves aerosol mass and number budgets within the do-10 main. The domain-integrated mass source and sink terms are due to surface flux of q ad , fallout of q ar , and mean vertical advection. Aerosol number has a more complex budget. One unforeseen complication was the limiter tendencies within the Morrison microphysics code necessary to prevent unphysical rain and cloud droplet size distributions. To conserve mass, the Morrison scheme adds or subtracts rain and cloud 15 number to keep the droplet size distributions within observationally derived bounds. We maintain a closed aerosol budget under these conditions by shifting aerosol number between N r or N d and N ad . A spurious source of total number (which we keep track of) can still result from the rare case when the required droplet number source exceeds the available dry aerosol number.  The thermodynamic sounding and winds used in this study are loosely based on the VOCALS-RF06 derived profiles of Berner et al. (2011). Changes include a reduction of the inversion height to 1300 m from 1400 m and initial boundary layer q t reduced to 10 7.0 g kg −1 from 7.5 g kg −1 . This results in a thinner and less drizzly initial cloud layer that does not deplete a large fraction of the initial aerosol during the spin-up of the simulations. The wind is forced using a vertically uniform geostrophic pressure gradient and the initial sounding is tuned to minimize inertial oscillations. The free-tropospheric moisture and temperature are nudged to their initial profiles (or downward linear extrap- 15 olations thereof, if the inversion shallows more than 150 m from its initial specification) on a one hour timescale in a layer beginning 150 m above the diagnosed inversion height. The sounding and winds averaged over hour three of the control case are depicted in Fig. 3.

20
Most simulations use diurnally averaged insolation and an insolation-weighted solar zenith angle appropriate for the latitude and time of year of the VOCALS-RF06 POC observation. Section 6.2 presents a sensitivity study using a diurnal cycle of insolation.

Subsidence
Mean subsidence is assumed to increase linearly from the surface to a height of 3000 m and to be constant above that. Throughout this paper, the subsidence profile is indicated using its value at a height of 1500 m, which is between 4.75-6.5 mm s −1 in the simulations presented. These values are somewhat larger than our best guess of 5 the actual mean subsidence at 1500 m during VOCALS-RF06 (2 mm s −1 ; Wood et al., 2011b). This subsidence range was chosen because it exhibits an interesting range of cloud-aerosol-precipitation interaction and long-lived regimes under steady forcing. A simulation with the observed subsidence produces an initial evolution qualitatively similar to the first case discussed below, but with more rapid initial cloud deepening 10 and onset of drizzle. In the VOCALS study region, there is a significant diurnal cycle of subsidence, but we have not included it here for simplicity, since in this paper our focus is on cloud-aerosol regimes which evolve mainly in response to the average forcing over longer periods of time.

15
The aerosol number concentration in the FT is set to 100 mg −1 , typical of FT CCN concentrations over the remote ocean observed in VOCALS (Allen et al., 2011). A sensitivity test with no FT aerosol is discussed in Sect. 6.1. Within the MBL, cases with vertically-uniform initial aerosol concentrations ranging from 100 mg −1 to as low as 10 mg −1 are considered, inspired by VOCALS-observed ranges within the MBL over 20 the remote ocean (Allen et al., 2011). Table 1 lists the simulations discussed in this paper and summarizes our naming convention. Slight random noise in the initial temperature field is used to initiate eddy motions. Before discussing the results in detail, we provide a brief synopsis. In the base-25 18159 line W5/NA100 run, the cloud layer deepens, thickens, and starts to drizzle, transitions to open cellular convection via strong precipitation scavenging, and then collapses; this evolution is also found in a three-dimensional version of this run, as well as in a larger-domain 2-D simulation. If W is raised to 6.5 mm s −1 , the cloud layer deepening is suppressed and a nearly nonprecipitating steady state stratocumulus layer develops.

5
The early development of these runs is reminiscent of prior simulations of Mechem and Kogan (2003) and Wang et al. (2010). A similar bifurcation between cloud thickening and collapse vs. development of a steady state is seen when the diurnal cycle is included. For the W6.5 case, we also consider identically forced runs with different initial aerosol concentrations (NA50, NA30, NA10). The NA10 run evolves through an 10 collapsing open-cell regime into a different thin-cloud equilibrium than the other runs. Finally, we consider "POC" simulations in a larger 2-D domain where the initial MBL aerosol includes a horizontal variation in initial MBL aerosol concentration. These runs produce a low-aerosol open-cell state that develops from the initially lower N a region, but in the 100 : 50 run, the surrounding region remains a well-mixed, nearly nonprecip- 15 itating stratocumulus topped boundary layer with higher aerosol concentrations; this POC-like combined state persists indefinitely given steady forcing.
The model behavior agrees qualitatively with observations in a number of important respects that increase our confidence in its applicability. In the VOCALS campaign, POCs were observed to form rapidly, almost always in the early morning hours when 20 LWP reaches its maximum (Wood et al., 2011b(Wood et al., , 2008). In the model, stable stratocumulus decks can persist with diurnally averaged LWPs up to around 100 g m −2 , but larger LWPs result in precipitation sinks of aerosols that cannot be balanced by reasonable source strengths. This agrees well with the LWP climatology of Wood and Hartmann (2006) for subtropical stratocumulus. The modeled transition from stratocumulus to 25 open cells occurs via strong precipitation feedbacks, again in agreement with observations (van Zanten and Stevens, 2005;Comstock et al., 2005;Wood et al., 2011b), and in an observationally reasonable period on the order of ten hours. The model also produces a post-transition vertical structure for aerosol in good agreement with obser-vations, with a surface layer N ad on the order of 20-30 mg −1 and a decoupled upper layer with much lower values, including an "ultra-clean layer" as observed in VOCALS RF06 (Wood et al., 2011b).

Evolution through multiple cloud-aerosol regimes
We begin our analysis with run W5/NA100. Figure 4 depicts time-height plots of horizontally-averaged total (dry plus wet) aerosol number concentration and liquid water content, as well as time series for important domain-averaged MBL meteorological variables and for individual tendency terms from the MBL-averaged aerosol number budget. To compute the budget, terms are calculated between the surface and the time-varying horizontal-mean inversion height z i , determined as the height at 10 which the domain mean relative humidity goes below 50 %. The entrainment source is w e (N aFT − N aMBL )/z i . Here w e is the entrainment rate, diagnosed from the difference between the z i tendency and the mean vertical motion at the inversion height. This approach is approximate, but the residual that it induces in the aerosol number budget is only a few percent of the dominant terms. The scavenging of interstitial aerosol num-15 ber by rain is not shown on the aerosol number budget plot, because it never exceeds 1 mg −1 day −1 , which is much smaller than the terms shown.
The plots reveal three distinct regimes: an initial period with a well-mixed stratocumulus-topped boundary layer, a transitional period with decoupled vertical structure, reduced cloud, and heavy precipitation, and a collapsed boundary layer state 20 with continual weak precipitation and sharply reduced aerosol concentrations.
Over the first two days, the inversion slowly deepens, the stratocumulus layer thickens and its LWP increases. There is a 25 % decrease in N a despite relatively negligible surface precipitation. Figure 4d shows that over the first day, the dominant aerosol loss term is interstitial scavenging by cloud, with smaller and roughly comparable losses from autoconversion and accretion, and a very small term due to aerosol number nonconservation from limiters within the microphysics scheme acting to keep the cloud size 18161 distribution within defined bounds. These losses are balanced primarily by the surface source, with a slowly increasing contribution from entrainment. As noted by Mechem et al. (2006), the FT may act as a reservoir for MBL CCN when the MBL-average CCN number concentration falls below that of the FT. During the second day, accretion losses begin to rise more sharply, while interstitial scavenging by cloud diminishes. This 5 occurs due to the improved collection efficiency for drizzle drops for constant q r and diminishing N r , while interstitial scavenging becomes less efficient as cloud droplets grow larger. During day three, accretion losses grow to dominate the number budget, while Fig. 4c shows entrainment sharply declining as the domain-averaged surface precipitation 10 rate climbs above 0.5 mm day −1 . This is an example of what Feingold and Kreidenweis (2002) called a "runaway precipitation process". The decrease in entrainment is due to a rapid decrease in turbulence near the top of the boundary layer, because of drizzle-induced stabilization of the boundary layer (Stevens et al., 1998) and reduced destabilization by boundary-layer radiative cooling as cloud cover reduces. The result 15 is a crash of MBL N a and a drastic reduction in LWP. Figure 4d shows the sink terms diminish rapidly at the end of day three, primarily because the vast majority of aerosol has been removed from the cloud layer. The transition from well-mixed stratocumulus to showery, cumuliform dynamics is abrupt. Figure 5 shows x-z snapshots of q l and N a at three times spanning a six-20 teen hour period (day 2.625 to day 3.295). At the first time, cloud cover is essentially 100 %, with q l maxima near 1 g kg −1 and a slight amount of cloud base drizzle under the thickest clouds. Total aerosol concentration N a has only a slight vertical gradient and is about 40 mg −1 in the cloud layer. At the second time, eight hours later, the cloud is thinning and breaking towards the center of the domain, and in a smaller section 25 of the domain there is surface precipitation beneath a strong drizzle cell with cloudtop q l in excess of 1.5 g kg −1 . The decoupled structure of the boundary layer inhibits aerosol transport from the surface layer into the remaining cloud. Hence, aerosol concentration has become quite vertically stratified, with cloud layer values depleted near or below 10 mg −1 , except for higher values around the strong updraft of the drizzle cell. At the final time, cloud cover has fallen below 50 %; one weak drizzle cell remains, but the cloud layer q l is nearly totally depleted. A narrow band of highly depleted aerosol concentrations sits 100-300 m below the inversion, with concentrations falling below 5 mg −1 (an "ultra-clean layer"), while the surface layer concentrations remain between 5 20-30 mg −1 . Referring back to Fig. 4d, the net MBL aerosol number sink rate during this 16 h period exceeds 60 mg −1 day −1 , allowing for near complete aerosol depletion in less than 24 h. The showery, cumuliform conditions following the sharp transition from stratocumulus persist for approximately two days. During this period, Fig. 4c shows that strong drizzle events occur periodically with spikes in surface precipitation up to 4 mm day −1 , cloud cover oscillates around 60 % (much of which is optically thin), and domain-mean LWP oscillates between 20-40 g m −2 . Entrainment is negligible, so the boundary layer continually shallows and the cumulus cloud layer thins. During the sixth day, surface precipitation becomes lighter and more continuous at around 1 mm day −1 , while liquid 15 water path falls to between 10-20 g m −2 , as the cloud layer becomes too thin to support episodic cumulus showers. Figure 4a also suggests there is less vertical stratification of the boundary layer aerosol. The boundary layer maintains this state while shallowing due to subsidence from a depth of 700 m down to 300 m. This slow boundary layer collapse is reminiscent of results from a one-dimensional turbulence closure modeling 20 study of Ackerman et al. (1993). After six days of slow collapse, the MBL is sufficiently shallow and the inversion is weak enough that thin cloud forms near the top of the boundary layer and reinvigorates entrainment. The sudden influx of entrained aerosol into the cloud decreases cloud droplet sizes, reducing precipitation efficiency and cloud processing sinks within 25 the boundary layer. This results in a positive entrainment-cloud-aerosol feedback that rapidly replenishes boundary layer aerosol, allowing the layer to sustain radiatively active cloud and begin deepening again via entrainment. While cloud type is initially thin Sc, as the layer deepens the thin cloud is unable to develop sufficient turbulence to 18163 mix the full layer depth, and a decoupled structure takes over with around day 10.5 with lower cloud fraction, slightly reduced entrainment, and a larger vertical aerosol gradient.
The boundary layer continues to deepen and moisten over the next few days with a similar boundary layer structure. On day 15, entrainment begins to strengthen again 5 due to larger LWP and radiatively driven turbulence; full cloud cover is achieved during day 16. While Fig. 4 only extends 20 days, similar simulations we have run past 22 days exhibit limit cycle behavior, collapsing again due to the runaway precipitation sink of aerosol. In a sensitivity test using a 96 km wide domain (run W5/NA100/LD), mesoscale variability allows a region of stratocumulus to form, begin replenishing MBL aerosol, 10 and restart inversion deepening earlier in the collapse process, when the inversion is still at a height of 700 m. This indicates the collapsed regime is fairly delicate in this model.

Comparison with 3-D results
To check the robustness of our 2-D results, we performed a two-day, 3-D simulation 15 W5/NA100/3-D identical in forcing and initialization to the W5/NA100 case, except using a domain of 24 km × 24 km in horizontal extent. The evolution of this run compared with the 2-D control run is depicted in Fig. 6. A few differences are apparent. Entrainment in the 3-D simulation is less efficient than in the 2-D run; the boundary layer first shallows during spin-up, then slowly deepens to a maximum of 1270 m after one day 20 of evolution, while the 2-D run is nearly 50 m deeper at this point. This results in less entrainment drying of the boundary layer, and a peak LWP of 175 g m −2 is reached in the 3-D case as compared to a maximum of 155 g m −2 in the 2-D case. The transition to a collapsing state begins after only 24 h in the 3-D run, while taking 54 h to reach in the 2-D case. This is due partly to the larger LWP accelerating precipitation losses, as well as the greater strength of cloud interstitial scavenging sink, which is nearly twice as large in the 3-D case. The combination of these effects leads to accelerated scaveng-ing of aerosol within the boundary layer and faster onset of the runaway precipitation sink.
While the strengths of the aerosol number source and sink terms differ quantitatively between the 2-D and 3-D configurations, their relative roles are similar. Interstitial scavenging is initially the largest sink of N a and the surface flux the largest source.
As multiple processes act to reduce N a , the entrainment source strengthens, but is ultimately unable to compete with the accretion losses as drizzle becomes significant. Cloud fraction decreases in both runs following the development of significant surface precipitation and the rainout of q l , demonstrating qualitatively similar dynamics. As the 3-D configuration is nearly 200 times as computationally expensive to run, we use the 10 2-D framework for the remainder of our simulations to enable the long run times necessary to examine aerosol-cloud regimes and equilibrium states.

Stable equilibrium and sensitivity to forcing
In the W5/NA100 run, cloud breakup occurs when LWP increases enough to induce a runaway drizzle-aerosol loss feedback. This suggests that this transition might be 15 suppressed if LWP remains sufficiently low. With this in mind, Case W6.5/NA100 increases 1500 m subsidence to 6.5 mm s −1 , limiting LWP by inhibiting the boundary layer from deepening. This case, shown in Fig. 7, reaches a steady, well-mixed, stratocumulus-topped equilibrium state with an inversion height of 1100 m. LWP settles to a mean of approximately 65 g m −2 after eight days, with oscillations of ±5 g m −2 about 20 this value thereafter. The boundary layer N a settles around 88 mg −1 , with 100 % cloud cover, less than 0.1 mm day −1 of cloud base drizzle and negligible precipitation reaching the surface. The steady state aerosol budget is dominated by an MBL-averaged surface N a source of 20 day −1 and an interstitial scavenging N a sink of 12 day −1 , with smaller contributions from other processes. 25 18165

Sensitivity to FT aerosol
Case W6.5/NA100-0FT is identical to Case W6.5/NA100, except that the FT aerosol concentration is set to zero. In this run, entrainment dilution is always a sink term for MBL-mean N a . After a gradual depletion of aerosol, the run-away precipitation sink transitions the system into a collapsed MBL. The simulation was run out to 20 days and 5 remained in a collapsed state, with the inversion continually sinking slowly throughout (no plots shown).

Sensitivity to diurnal cycle of insolation
Case W5/NA100/DC (Fig. 8) is configured similarly to the control run W5/NA100, except with a diurnal cycle of insolation. Surprisingly, the inversion does not initially deepen as fast as in the control run, because the daily-mean entrainment rate is lower. This occurs because the cloud dissipates during the daytime, leading to a diurnallyaveraged reduction in longwave cooling and turbulence compared to the control case. As a result, the cloud layer evolves towards a steady diurnal oscillation rather than thickening until it undergoes the runaway precipitation feedback. The cloud LWP exceeds 15 200 g m −2 for a few hours at night, but this is not long enough for accretion losses to build up before the cloud thins during the daytime and aerosol regenerates. A curious feature of the case is a nearly 2 m s −1 oscillation in simulated surface wind speed, driving a large simulated diurnal cycle in surface aerosol number flux(which is proportional to the cube of the wind speed). This effect is attributed to reduced downward turbu-20 lent mixing of momentum to the surface layer during the daytime. Observational data show a much weaker diurnal cycle in windspeed of a few tenths of a meter per second (Dai and Deser, 1999); we speculate that our 2-D approach is artificially amplifying the simulated wind and surface aerosol flux oscillation, although we don't think this has a major impact on our results. The case was rerun with the 1500 m subsidence rate re-25 duced 5 % to 4.75 mm s −1 , and behavior similar to the control case reappears (Fig. 9). The rapid reduction in cloud LWP and cloud fraction occurs in the early morning when precipitation is most active, as in the formation of observed POCs (Wood et al., 2008). The sensitivity of the long-term behavior to this slight change in an external parameter is an indicator of a positive-feedback system. The diurnal modulation of the terms in the aerosol number budget does not seem to alter their relative importance in each cloud regime compared to the control case. 5 7 A reduced-order phase-plane description of the aerosol-cloud system Schubert et al. (1979) discussed characteristic timescales on which a stratocumuluscapped mixed layer adjusts to a sudden change in boundary conditions and forcings. They pointed out that there is a quick (few hours) thermodynamic adjustment of the MBL, followed by a slower adjustment timescale (several days) for the MBL depth to adjust into balance with the mean subsidence. Bretherton et al. (2010) elaborated these ideas using both mixed-layer modeling and LES of stratocumulus-capped boundary layers with fixed cloud droplet concentrations. They showed that with fixed boundary conditions, for any initial condition, the MBL evolution converged after thermodynamic adjustment onto a "slow manifold" along which the entire boundary-layer thermody-15 namic and cloud structure was slaved to a single slowly-evolving variable, the inversion height (z i ). With a cloud droplet concentration of 100 cm −3 they found there were two possible slow manifolds, a "decoupled" manifold evolving toward a steady state with small cloud fraction and a shallow inversion, and a "well-mixed" manifold evolving toward an solid stratocumulus layer with a deep inversion. In both equilibria, precipita-20 tion was negligible. Simulations initialized with well-mixed boundary layers capped with a cloud layer that was optically thick but non-drizzling converged onto the well-mixed manifold; simulations in which the initial cloud layer was either optically thin or so thick as to heavily drizzle converged onto the decoupled manifold.
In this section, we explore the use of similar concepts for our cloud-aerosol sys-25 tem. Given the MBL-average aerosol N a , the clouds and turbulence will pattern the aerosol sources and sinks to set up the vertical structure of N a within the MBL within 18167 its turbulent overturning time of a few minutes (for a coupled boundary layer) to a few hours (for a decoupled boundary layer). Hence, the combination of N a and z i provides a reduced-order two-dimensional phase-space to describe the "slow manifold" evolution of the cloud-aerosol system on timescales of a day or longer. We start by viewing the control W5/NA100 case in this way. Figure 10 shows the 5 position in N a -z i phase space averaged over sequential 12 h periods, colored by the LWP during that period. The first 12 h adjustment period of the vertical structure of the boundary layer and the aerosol is shown in light grey, and the second in dark grey. The spacing between successive points indicates how fast a run is evolving. Different "regimes" through which the system evolves are labeled. Here, a cloud-aerosol regime is defined as a part of phase space with qualitatively similar cloud and aerosol characteristics, and hence different balances of terms in the aerosol budget. For instance, the open cell regime is characterized by low aerosol, low LWP, low entrainment and efficient precipitation (accretion) removal of aerosol, while the thick cloud regime is characterized by high aerosol, high LWP, little precipitation, high entrainment, and 15 a balance between cloud scavenging and surface/entrainment aerosol sources. The regimes grade into each other as N a and z i change; they need not have sharp boundaries in phase space. The shading indicates two regions of phase-space in which the aerosol concentration evolves comparatively rapidly, either due to runaway precipitation feedback, or due to rapid "runaway" entrainment of aerosol when the shallow boundary 20 layer with open-cell convection redevelops thin inversion-base stratocumulus clouds. Overall, the phase space trajectory has converged onto a limit cycle in which it will indefinitely move between the thick-cloud, open-cell, and thin-cloud regimes. It is even more interesting to examine the W6.5 case in this framework. The inverted triangles in Fig. 11 show the W6.5/NA100 evolution in N a -z i phase space. This run 25 shows a slow decrease in z i , with an initial decrease of N a , later changing to a slight increase during the final approach to a "thick-cloud" steady state at z i = 1100 m and N a = 88 mg −1 .
We performed a sequence of runs identical to W6.5/NA100, but with the initial MBL aerosol concentration varied to 10, 30, and 50 mg −1 ; results are also plotted on Fig. 11.
The simulations with initial N a values of 30 and 50 mg −1 (five and six pointed stars, respectively) evolve to the same thick-cloud equilibrium as the NA100 case. The initial z i drop is faster for lower N a , because increased drizzle inhibits entrainment. After 5 a few days, both simulations have settled to nearly their equilibrium z i but still slowly drift toward higher N a as they approach the steady state. The run with initial N a of 10 mg −1 (diamonds) immediately enters a runaway precipitation feedback as it spins up, transitioning after the first 12 h to a low-LWP collapsing boundary layer. The "elbow" in the lower left corner of the phase space indicates a time 10 at which the MBL eventually recovers a thin layer of stratocumulus and starts to rapidly re-entrain aerosol, as in the W5/NA100 case. In contrast to that case, the stratocumulus layer never thickens enough to support vigorous turbulence, and the entrainment is only able to slightly deepen the boundary layer, settling into a second stable "thincloud" equilibrium that is different than for the other initial conditions, with high N a and 15 no drizzle.
The labelled regimes occupy the same parts of phase space as for W5, because the aerosol loss rate is not directly affected by the instantaneous mean subsidence rate. However, increased subsidence causes individual simulations to evolve through phase space differently than for W5. On Fig. 11, we have sketched in some additional 20 features whose exact structure we can only guess at. The grey and red dashed curves indicates the upper and lower boundaries of attractor basin for the thick-cloud steady state, i. e. all points in phase space from which the slowly-evolving system converges to that steady state. Although the NA50 and NA30 runs appear to start outside this attractor basin, note that it is only after the initial thermodynamic adjustment, which 25 takes about a day, that the system behavior locks into this slowly-evolving phase-plane description. By that time these runs are well within the notional boundary of the attractor basin. The sketched arrow indicate the rough direction in which the system evolves through different parts of phase space. Between the two steady states, the red dashed 18169 curve separates trajectories for which the cloud thickens and the inversion deepens toward the thick cloud equilibrium, and ones for which the cloud thins and the inversion shallows toward the thin cloud equilibrium. This is an example of sensitive dependence of the long-term system evolution to its initial state. From the simulations we have performed, the open-cell regime is quite robust, so we anticipate that if N a is slightly 5 larger than on this trajectory, it will adjust downward onto a similar trajectory as the NA10 simulation. For large z i > 1300 m, the W5/NA100 simulation suggest that LWP is very large, increasing accretion loss of N a . Thus we hypothesis that both the attractor boundary and the "runaway precipitation" area of very rapid aerosol depletion move rightward as z i increases. With an expanded set of runs that varies the initial z i in addition to N a , the phase plane evolution could be filled out more confidently and precisely.
The two equilibria of the W6.5 phase plane are strikingly similar to those found by Bretherton et al. (2010). This suggests that there may be parts of phase space with two possible slow manifold behaviors, depending on the initial state of the system. 15 However, unlike in Bretherton et al. (2010), initial states with very thick cloud will quickly remove their aerosol through precipitation accretion, forcing them into a low-aerosol open-cell state before returning back to the thin cloud equilibrium via runaway aerosol entrainment.
Unlike Baker and Charlson (1990), we do not find a high-aerosol and a low-aerosol 20 steady state for W6.5. However, the open-cell regime is quite long-lived, and without free-tropospheric aerosol, the runaway-entrainment transition to the thin cloud, highaerosol equilibrium would not occur. Thus, our simulations are consistent in spirit with their hypothesis. In our work, however, interstitial aerosol removal is much more efficient, and entrainment dilution is active, regulating the high-aerosol state to number 25 concentrations of around 100 mg −1 , about ten times smaller than in their simple model.

POC simulations
So far, we have considered the evolution of an aerosol-cloud system in a small domain through multiple regimes over time, with the possibility of multiple steady states for a given forcing. We now examine the possibility of simultaneously supporting adjacent regions in different aerosol-cloud regimes in an idealized representation of a POC. We extend our domain to 192 km and specify an initial gradient in aerosol concentration, with a 42 km region initialized to N a + , followed by a 12 km half-sine wave transition to 84 km initialized at N a − , then another half-sine wave transition back to 42 km initialized at N a + . POC runs are initialized with identical thermodynamic profiles and forcings to the W5/NA100 case. The idea is that the lower N a cloud within the center of the domain 10 (the incipient POC) will start drizzling and transition via runaway precipitation feedback to open cell structure. This will reduce entrainment in those regions. Because the inversion height in the open cell region must stay similar to that in the overcast region, reduced entrainment in part of the domain will slow or reverse inversion deepening over the entire domain, as discussed in Berner et al. (2011). This can prevent the further 15 growth of LWP and hence drizzle in the overcast region, and thereby suppress the transition there, allowing a "coupled slow manifold" behavior hypothesized by Bretherton et al. (2010) of this section documents this remarkable behavior. Figures 12 and 13 show x-z cross-sections of N a and q l at days 0.475 and 1.205. In Fig. 12, the boundary layer is vertically well-mixed in N a and topped by an over-18171 cast stratocumulus layer. No significant drizzle is present, and no strong differences other than the N a gradient are visible across the domain. Figure 13, 17.5 h later, shows a picture strongly reminiscent of observations from a POC, with several strong drizzle cells present in the central region and broken cloud. N a is highly depleted in places, falling below 15 mg −1 in the upper layer of the POC, and below 2 mg −1 in the bright ma-5 genta regions, nicely capturing the observed ultra-clean layer (Wood et al., 2011b). The surrounding overcast region maintains full cloud cover without substantial surface precipitation, and N a remains above 75 mg −1 in the majority of the outer region. Another model feature resembling observations of the VOCALS RF06 POC is the location of the horizontal N a gradient, with the bulk of the aerosol decrease located within the overcast 10 region (Wood et al., 2011b). The qualitative similarity to key features in the observations suggests the simple model captures at least some of the important mechanisms in the real system. After transition, the POC and overcast regions evolve jointly over the next several days, linked together by circulations acting to maintain an essentially level domain- 15 wide inversion against large differences in entrainment and thermodynamic profiles between the overcast and POC regions. Figure 14 depicts Hovmöller plots of N a and LWP. The spatial evolution of the aerosol and LWP fields in Fig. 14a and b shows how the POC initially widens rapidly following transition. The different character of the regions is apparent in Fig. 14b, in that the domain maximum and minimum LWP values 20 are contained within the POC, while the overcast region is much more homogeneous. Figure 14a shows the degree to which the aerosol gradient is sharpened and maintained by cloud processing. Interestingly, after its initial expansion, the POC reaches a maximum of 70 % of the domain on day 2.4, then shrinks back into a quasi-steady equilibrium with the overcast by the time the run ends on day 9, occupying 25-30 % of 25 the total area.
In order to understand the apparent equilibration, we plot time series of several variables from averages over 24 km wide stripes at the center of each region in Fig. 15. In Fig. 15a, the transition to open cellular convection within the POC at the end of the first day is evident, with a sharp decrease of N a there and corresponding loss of cloud fraction (Fig. 15f); within the overcast region, N a diminishes less rapidly. Inversion heights fall in both regions after the transition, as depicted in Fig. 15b. This is a response to the reduced domain-mean entrainment resulting from the growth of the POC, shown in Fig. 15c. The falling inversion thins the stratocumulus clouds, reducing 5 LWP (Fig. 15d) and cloud processing of aerosol in the overcast region, in turn allowing N a to increase. This further limits precipitation sinks of both aerosol and cloud water, and overcast LWP begins to recover soon after surface precipitation ceases on day three. The inversion height and LWP then increase in unison, and as the cloud thickens, stronger cloud processing and buffering by entrainment dilution (Mechem et al., 2006) arrests the growth of overcast N a . This suggests a negative feedback loop, where the coupling between POC area, domain mean entrainment, and LWP as moderated by z i stabilizes the system. Figure 16 shows the nearly equilibrated behavior on day 9.58. Strong latent heating in the upper region of the POC contrasts with continuous radiative cooling of the well 15 mixed overcast. We hypothesize the 50 m difference in inversion height between the POC and surrounding overcast represents a balanced response to this gradient in vertically integrated buoyancy. While the POC has diminished in area, drizzle processes are still very active. The cloudbase smoothly drops from the edges of the overcast, reflecting the cooler, decoupled POC surface layer, and only the most active drizzle cells 20 have cloud top heights near the inversion. Without the surrounding overcast region, the inversion in the POC region would have collapsed several hundred meters due to the lack of entrainment, while that in the overcast region would have deepened due to the large entrainment there and undergone runaway precipitation feedback, as in Fig. 4. That is, each region holds the other in balance. 25 Interpreting this system in a reduced-complexity slow-manifold framework, the quasiequilibrium evolution of each region is controlled its mean N a and by the coupled evolution of inversion height, i. e. domain-mean entrainment, which is in turn set, to first order, by the areal fraction of the domain occupied by the overcast region vs. 18173 the POC. It is unclear exactly what controls the areal fraction of POC to overcast, but reduced POC entrainment couples to domain-wide boundary layer depth, LWP, and aerosol tendencies in a manner that creates a negative feedback on the horizontal expansion of the POC. 5 We have implemented a new single log-normal mode, double-moment bulk aerosol scheme coupled to the Morrison microphysics parameterization for the SAM LES. The scheme follows the general approach of Ivanova and Leighton (2008), though we have used only a single dry mode for simplicity and computational efficiency. We use this framework to examine the evolution of the cloud-aerosol system through different 10 regimes, sensitivity to forcing and initial conditions, and the formation and equilibrium of POCs.

Conclusions
As in some past studies, transition from closed-cell stratocumulus to more cumuliform, open-cellular convection occurs via the sharp enhancement of the accretion sink as LWP increases and N d decreases. If LWP does not get sufficiently high, this tran-15 sition never takes place. With diurnally-varying insolation, the transition occurs in the early-morning hours, as with observed POC formation.
A series of identically forced runs varying the initial boundary layer aerosol concentration demonstrated multiple equilibria, with higher N a runs settling into a deeper, well-mixed state and a run initialized with N a of 10 mg −1 transitioning to a collapsed 20 boundary layer before recovering to a shallower layer with very thin cloud and little liquid water content and high N d . A 2-D phase plane representation of the evolution of the cloud-aerosol system in terms of boundary-layer mean aerosol concentration and inversion height provides a convenient way to identify cloud-aerosol regimes and transitions between them, and basins of attraction to the respective equilibria. In spirit, our 25 results support the bistability hypothesis of Baker and Charlson (1990), though modi-fied in important ways by entrainment feedback and interstitial aerosol scavenging by cloud droplets. Another series of runs used a larger domain with a horizontal aerosol gradient to explore the formation of POCs. When initialized with an outer domain N a of 100 mg −1 and an inner N a of 50 mg −1 , the inner region transitioned to open-cellular cumuliform 5 convection while the outer region remained well-mixed, closed-cell stratocumulus. The transition of the inner region reduced domain averaged entrainment, reducing the overall inversion height and reducing cloud thickness and LWP in the outer region, resulting in increased N a there and providing a negative feedback on growth of the open-cellular region. The resulting POC-overcast system appears to be a stable steady state under 10 constant forcing. This is consistent with the observed robustness of POCs, which are rarely observed to close up once formed (Wood et al., 2008). Future work will involve continued exploration of the phase space in an effort to more fully map the behavior of the system, including longer three dimensional runs. Three-dimensional POC simulations will allow for the development of more complete 15 regional aerosol budgets, particularly characterizing the role of horizontal transport of aerosol from the overcast into the POC. Another area of potential research involves the sensitivity of nascent and fully developed POCs to FT aerosol perturbations, as the question of whether it is possible to close a POC remains open.

18175
Appendix A

Interstitial Scavenging
A log-normal form f (A) is assumed for the aerosol size spectrum. The number and mass scavenging tendencies are: where A is the aerosol diameter, ρ a is aerosol density, and γ(A) is the scavenging coefficient (e.g. Seinfeld and Pandis, 1998) (Khvorostyanov and Curry, 2002). The collector drop size distribution is assumed to follow a generalized gamma distribution F (D) = N 0 D ν exp(−λD), where N 0 is the intercept parameter, ν the shape parameter, and λ is the slope parameter (e.g. Ulbrich, 1983).
With the above simplifications, the scavenging coefficient becomes: In Morrison et al. (2005) and other double moment schemes, N 0 and λ may be calculated from the total number (N d for cloud, N r for rain) and water mass (q c for cloud, q r for rain) where ρ w is the density of water. Collision efficiency is given as the effective collision kernel over the geometric colli- , where K (D, A) is the sum of the individual process kernels. For cloud droplets, the included processes are Brownian diffusion, thermophoresis, diffusiophoresis, and turbulent coagulation (using Pruppacher and Klett, 1997, Eqs. (11-59), (17)(18)(19)(20)(21)(22)(23)(24)(25), (17-33), 15 and (11-77), respectively). The turbulent kernel currently assumes a fixed value for dissipation in convective clouds of = 4.22 × 10 −5 m 2 s −3 ; this will be updated to use the local model values of epsilon in future versions. For rain, K (D, A) is the semi-empirical formulation of Slinn (1983). Appropriate values of sulfate aerosol thermal conductivity for these kernels may be found in Seinfeld and Pandis (1998, p. 481). Figure 17 plots 20 18177 logE (D, A) for aerosol diameter vs. hydrometeor diameter at a pressure of 900 mb, temperature of 282 K, and relative humidity (RH) of 100.5 %. There is some sensitivity to log E (D, A) around RH values near 100 % for aerosol larger than 0.1 µm and cloud droplets between 5-20 µm; as the version of SAM used here does not include prognostic supser-saturation, we choose an in-cloud value of 100.5 %. The discontinuity at D = 80 µm shows the transition between cloud and rain kernels. Since collision efficiencies for larger hydrometeors and accumulation mode (sub-micron diameter) aerosol are small (mostly below 1 %), this discontinuity is of little concern. The number and mass scavenging rates (A1) and (A2) are integrated numerically in a method generalizing the approach of Berthet et al. (2010). 10 Look-up tables of γ(A) for cloud and rain as a function of temperature, pressure, RH, and cloud/rain precipitation fluxes are generated on initialization of the scheme. This is done by integrating A3 across a default gamma distribution for cloud and the Marshall-Palmer distribution for rain using Gauss-Laguerre quadrature. Gauss-Hermite quadrature is then used to integrate over the aerosol size distribution using the current values 15 of N a , q a , and look-up values of γ(A) to calculate the scavenging tendencies of N a and q a . For efficiency, the tendencies are updated once every three timesteps and are set to zero after all liquid or aerosol has been removed from a gridbox between updates. Future versions of the code will use the prognostic Morrison cloud and rain distributions when numerically integrating (A3) and include an option to calculate tendencies using 20 an approximate analytical approach.  Clarke et al. (2006) size-resolved sea-salt aerosol number flux and unimodal CCN approximation used in this study, plotted for u 10 of 9 m s −1 . Blue vertical line indicates cut-off for viable CCN at 0.07 µm. All aerosol particles in this study are assumed to have the properties of ammonium sulfate.     Fig. 10. Phase-plane plot for slow evolution of control run W5/NA100. Symbols represent 12 h box averages and are colored by LWP. First and second 12 h averages is shown in light and dark grey, respectively, because the vertical profile of cloud and LWP is still undergoing an initial fast adjustment.    x-z sections from W5/POC-100-50 on day 9.58. Panels as in Fig. 12.  Fig. 17. Plot of collision efficiency as a function of aerosol and hydrometeor diameters at a pressure of 900 mb, temperature of 282 K, and relative humidity of 100.5 %. The rain collision efficiency is plotted for hydrometeor diameters above 80 µm; the cloud collision efficiency is plotted below this threshold.