Exploring Impacts of Size‐Dependent Evaporation and Entrainment in a Global Model

While most observations indicate well‐buffered clouds to aerosol perturbations, global models do not. Among the suggested mechanisms for this discrepancy is the models' lack of connections between cloud droplet size and two processes that can contribute to reduced cloudiness when droplets become more numerous and smaller: evaporation and entrainment. In this study, we explore different implementations of size‐dependent evaporation and entrainment in the global atmospheric model CAM5.3‐Oslo. We study their impact on the preindustrial‐to‐present day change in liquid water path ( LWPPD‐PI ) and the corresponding aerosol indirect effect ( AIEPD‐PI ). Impacts of the 2014–2015 fissure eruption in Holuhraun, Iceland, are also presented. Our entrainment modifications only have a moderate effect on AIEPD‐PI (changes from − 1.07 W m −2 to − 0.98 W m −2 ), and a small impact on the signal from the Holuhraun eruption compared to other suggested compensating mechanisms. Simulations with added size‐dependent evaporation in the top of the stratiform clouds also show small evaporation differences between PI and PD. Moderate changes in AIEPD‐PI were achieved when also including an entrainment feedback to the evaporation changes, mixing air between the cloudtop layer and the layer above. These changes were not associated with the size dependency, but changes in the cloud susceptibility to aerosols in both PI and PD when adding evaporation. We find that increased evaporation of smaller droplets at stratiform cloud tops can reduce LWP , but can increase LWP in some areas due to enhanced shallow convection caused by destabilization.


Introduction
Aerosol-cloud interactions are singled out by the Intergovernmental Panel on Climate Change (IPCC) as the largest contributor to the uncertainty in present-day global radiative forcing Myhre et al., 2013). An increase in the atmospheric aerosol concentration usually favors more numerous, but smaller cloud droplets that scatter more of the incoming solar radiation back to space (Twomey, 1977). The uncertainty in this so-called "Twomey effect," "albedo effect," or "first aerosol indirect effect," has been reduced during the last decades, and there is strong confidence that it is cooling the climate on a global scale (Bellouin et al., 2013;Bréon et al., 2002;Ghan et al., 2016;Feingold et al., 2003;Lebsock et al., 2008;McCoy et al., 2017;Quaas et al., 2008). However, estimates of the total magnitude of radiative effects by aerosol-cloud interactions remain highly uncertain, mainly due to the modulation of clouds through fast cloud adjustments that follow a reduction in cloud droplet size. This effect is often referred to as the "second aerosol indirect effect." Since smaller cloud droplets decrease the efficiency of precipitation formation by suppressing the collision-coalescence process, aerosols may increase the water content and the lifetime of clouds (Albrecht, 1989;Pincus & Baker, 1994). This is illustrated in the left branch of Figure 1, which shows how different cloud processes can affect LWP when the cloud droplet size is reduced. This so-called "Albrecht effect" or "lifetime effect" has, for a long time, been the only cloud process included in global models that can affect how aerosol perturbations impact cloud water content or extent through cloud physics (Kirkevåg et al., , 2018Neale et al., 2012;Walters et al., 2014;Zhang et al., 2012). As a consequence, estimates by these models of the total effective radiative forcing through aerosol-cloud interactions (ERF aci , including both the first and the second aerosol indirect effect) are not complete . to aerosol perturbations are questionable due to the difficulty of measuring aerosols and clouds simultaneously. Toll et al., 2017Toll et al., (2017Toll et al., , 2019 compared satellite measurements of clouds downwind of known aerosol sources, like cities, volcanoes, shiptracks, and fires, to unpolluted clouds nearby, showing a weak average decrease in LWP in the polluted clouds. Gryspeerdt et al. (2019) found both positive and negative relations between N c and LWP, with a domination of reduced LWP in clouds with high N c . McCoy et al. (2018) used a combination of satellite measurements and aqua planet simulations to show that high values of LWP in midlatitude cyclone clouds are associated with high values of N c . Malavelle et al. (2017) studied satellite retrievals of LWP during the 2014-2015 fissure eruption in Holuhraun, finding no significant increase in LWP compared to previous years. Thus, the sign of the second aerosol indirect effect is under debate, which makes it worth investigating cloud processes that can contribute to a reduction in cloudiness associated with increased aerosol loadings.
Sensitivity tests turning off the dependency of precipitation on droplet number in global models have shown better agreement with observations (Gettelman, 2015;Gettelman et al., 2013;Malavelle et al., 2017;Quaas et al., 2009;Wang et al., 2012), but cannot explain situations where LWP decreases with increased droplet concentrations. Another suggested reason for observed constant or decreasing LWP with increasing aerosol loading is the lack of representation of processes in the models that buffer cloud responses to aerosol perturbations (Stevens & Feingold, 2009). Some of these processes are seen in Figure 1. The middle branch shows that more numerous, but smaller droplets result in more evaporation. This size-dependent evaporation mechanism is not represented in most global models (Ackerman et al., 2004;Malavelle et al., 2017;Zhou & Penner, 2017). As for the growth rate of cloud droplets by condensation, the rate of evaporation also varies with the cloud droplet size. The characteristic phase relaxation time, ∝ (Nr) −1 (Squires, 1952a), where N is the cloud droplet number concentration and r is the mean radius of the cloud droplets, is used as a measure of the conversion rate between the liquid and the vapor phase. If a clean cloud (N 1 , r 1 ) is perturbed to a polluted cloud (N 2 > N 1 , r 2 < r 1 ) and the total liquid water content (M c ) is conserved, the phase relaxation time is reduced (see supporting information, Text S1), so that This change in is due to the larger total surface area compared to the total volume of the more numerous but smaller cloud droplets in the polluted cloud.
As seen in Figure 1, enhanced evaporation of cloud water will not only directly result in decreased LWP, but also trigger other fast feedback processes. When water at the cloudtop evaporates, the temperature drops, which promotes sinking of air masses, enhancing the turbulent mixing and entrainment of air from above, which again can result in more evaporation. Results from fine-scale modeling have shown examples of decreased cloud fraction and LWP caused by increased aerosol loading, probably due to enhanced evaporation and evaporation-entrainment feedbacks when the droplets become smaller (Ackerman et al., 2004;Altaratz et al., 2008;Feingold et al., 2006;Hill et al., 2009;Xue et al., 2008;Zhou & Penner, 2017). Small et al. (2009) presents observations from a flight campaign of nonrecipitating clouds that support these results. Bretherton et al. (2007) also points out that sedimentation of droplets away from the entrainment zone is suppressed when droplets are smaller, which can result in more liquid water available for evaporation in the entrainment zone, more evaporation and thus enhanced entrainment efficiency and reduced LWP. This is seen in the right branch of Figure 1.
The aim of this study is to implement size-dependent evaporation and size-dependent entrainment into a global model and see how that affects the second aerosol indirect effect. The hypothesis is that enhanced evaporation and entrainment in a polluted cloud will counteract the increase in water content due to suppressed autoconversion. This will make a perturbed cloud less reflective compared to a similar cloud without size dependency. The orange boxes indicate where size dependency is included in the modified simulations of this study, which will be described in detail in section 3.
Size-dependent evaporation and entrainment are implemented separately into the model, and both setups are described in section 3.  . Recent volcanic eruptions are, on the other hand, better documented and well observed by both satellites and ground-based instruments, and can serve as testbeds when modeling aerosol-cloud interactions. In Section 4.3, we investigate how the inclusion of size dependency affects the modeled changes in cloud properties caused by a volcanic eruption. We also compare the impact of implementing size dependency to the impact of corrections of other processes that have been suggested to affect modeled aerosol indirect effects.

General Description
This study uses the global model CAM5.3-Oslo, which is the atmospheric component of the Norwegian Earth System Model (NorESM). While CAM4-Oslo was used in the first version of NorESM Iversen et al., 2013;Kirkevåg et al., 2013) and CAM6-Oslo (under development) will be the atmospheric component of the second version of NorESM, CAM5.3-Oslo is an intermediate version of the atmospheric component, documented by Kirkevåg et al. (2018). It is based on the Community Atmospheric Model version 5.3 (Liu et al., 2016;Neale et al., 2012). The main difference between CAM5.3 and CAM5.3-Oslo is that the latter has its own aerosol module, OsloAero. This module is neither modal or sectional, but so-called "production tagged". In OsloAero, background tracers form lognormal distrubutions, but process-tracers resulting from condensation, coagulation and cloud processes change the shape of the distribution. Lookup tables produced by a sectional model are applied to obtain optical properties. More detailed description of the module is found in Kirkevåg et al. (2018). The stratiform clouds are treated by the double moment bulk microphysics scheme MG1.5, which is almost the same as MG1 , but with cloud droplet activation moved before the rest of the microphysical process rate calculations (Gettelman, 2015). The scheme applies autoconversion parameterizations based on Khairoutdinov and Kogan (2000). The treatment of shallow convective clouds and moist turbulence are based on work done at the University of Washington , while the parameterizations of deep convective clouds are based on Zhang and McFarlane (1995). Aerosols are activated following Abdul-Razzak and Ghan (2000). The macrophysics is described in Park et al. (2014) and includes saturation adjustment.
The blue boxes and the arrows in between in Figure 2 show the order of some of the most important processes in CAM5.3-Oslo that affect the cloud water mass, q c , and droplet number, N c . The aerosols are only impacting the stratiform clouds, which means that aerosol-cloud interactions lack from the convective clouds. Detrained convective condensate is added to the stratiform clouds, with a number concentration estimated by assuming a constant mean volume cloud droplet radius of 10 and 8 μm for shallow and deep convective clouds, respectively. The implementation of size dependency, indicated by the orange boxes in Figure 2, are described later.
Emissions of aerosols and aerosol precursor gases, and the concentrations of greenhouse gases in the atmosphere follow the standard of the Coupled Model Intercomparison Project number 5 (CMIP5) (Lamarque et al., 2010). The atmospheric component is coupled to a satellite phenology version of the Community Land Model version CLM4.5 (Oleson et al., 2013) and an ocean with prescribed climatological sea surface temperatures and sea ice extent. Methods by Ghan (2013) are used for calculating the effective radiative forcing of aerosols. For more detailed descriptions of the model, see Kirkevåg et al. (2018) and Karset et al. (2018).

Configurations
The model was configured with 30 levels in the vertical and a horizontal resolution of 0.9 • (latitude) by 1.25 • (longitude). Nudged meteorology was applied to constrain the natural variability (Kooperman et al., 2012). This meteorology was generated by the model itself in an earlier simulation, applying all present-day conditions (emissions, oxidants, land use, sea surface temperatures, etc.). Equation (2) describes how the horizontal wind components and the surface pressure were nudged from the old modeled value X m,old to the new value X m,new by nudging to the prescribed value X p from the meteorology previously produced by the model.
A relaxation time scale, , of 6 hr was used, corresponding to a nudging intensity of 8.3% when the time step of the model, Δt, is 30 min. The air temperature was freely evolving (Zhang et al., 2014), as were all other variables.

Experimental Setup
To find out how the clouds in the model respond to aerosol perturbations, two simulations with the default model setup, N, were carried out. The simulations applied emissions of aerosol precursor gases, aerosols, and prescribed oxidants from present day (PD, year 2000) and preindustrial times (PI, year 1850), respectively. The simulation length is 6 years, where the last four are analyzed. Simulations from a previous study, applying a similar setup, showed a standard error of only 0.01 W m −2 for the total aerosol indirect effect, and the same result, with no drift in the signal, when extending the simulation length to 11 years and analyzing the last 10 years (Karset et al., 2018). The results from the nudging simulations of that study were also in the uncertainty range of the simulated total aerosol indirect effect with free meteorology, running for 50 years, where the last 30 were analyzed. The impact of including size dependency on the evaporation and the entrainment processes are investigated separately in the following sections.

Size-Dependent Evaporation
An additional evaporation is implemented in the top layer of clouds without ice. The argument of adding extra evaporation, on top of the one already calculated by the model itself, is that we assume the model is underestimating the amount of evaporation and entrainment, thus the reduction in cloud water, due to the lack of an evaporation-entrainment feedback parameterization. Size-dependent evaporation occurs in both warm and mixed-phased clouds, but we choose to only focus on the warm clouds to better understand the physical mechanisms behind the results, without possible impacts by complex processes including ice. In a later study, mixed-phase clouds should also be included. Figure 3a shows the effective radius of the cloud droplets in the top layer of stratiform clouds without ice in the default PD-simulation, while Figure 3b shows the difference in the effective radius between PI and PD. The effective radius of cloud droplets in CAM5.3-Oslo is slightly lower than that of MODIS (Moderate Resolution Imaging Spectroradiometer) (Rausch et al., 2017), which we find acceptable due to the positive biases associated with the r e -products coming from MODIS (Liang et al., 2015;Painemal & Zuidema, 2011). Changes in r e between PI and PD are mostly negative, owing to the Twomey effect when aerosol concentrations increase. Positive areas over land are associated with reduced emissions from forest and grass fires since PI, while positive areas over remote oceans are due to oxidant differences between PI and PD, resulting in a shift in the location of secondary aerosol formation (Karset et al., 2018).
As a first approach to the question of how size-dependent evaporation will affect cloud responses to aerosol perturbations, 5% of the liquid in the cloudtop layer is evaporated if the cloud droplets are smaller than 4 μm. The choice of a lower limit of 4 μm is not observational based, and smaller values have been tested without changing our conclusions. Since the evaporation changes can occur every time step of 0.5 hr, it results in an evaporation rate in the cloud top layer of up to 10% hr −1 . As we do not have good global constraints by observations of this quantity at this point, a higher evaporation factor of up to 50% is also tested. This corresponds to an evaporation rate in the cloud top layer of up to 100% hr −1 .
The rate of evaporation depends on the surface area. If the total volume of water in a cloud, V, is constant, the total surface area, A, of all the cloud droplets is given as Since A is proportional to r −1 , we let the additional evaporation also vary with r −1 (Rogers & Yau, 1989) for droplets larger than 4 μm. This results in an evaporation rate seen in Figure 4a. Note that r in equation (3) refers to the volume mean droplet radius, r v , while we in Figure 4a and in the modified model setups apply the effective radius, r e , which is an area weighted droplet radius. Since r e and r v usually are linearly proportional (Freud & Rosenfeld, 2012;Martin et al., 1994;Wood, 2000), and we are not interested in using the absolute value, but scale the processes by it, our conclusions are not changing when switching between the two. We simply applied r e since it is more available in all parts of the modified model code. The horizontal distribution of the evaporation rate from the PD-simulation is seen in Figure 4b with its global mean value in the heading.
The standard evaporation calculations in the model is a part of the macrophysics (see Figure 2). The evaporation changes in this study are applied to the model code right after the microphysics (see Figure 2). One reason why we do not add it right after the standard evaporation in the macrophysics is the way activation is calculated. The activation process makes assumptions based on saturation-adjustment, which is not the case if we add extra evaporation. If trying to evaporate cloud droplets between the macrophysics and activation, the scheme will activate back the number, but not the mass, resulting in an unphysical reduction in cloud droplet size. Another reason why the code changes are applied where they are is because a the residual  condensation term in the microphysics scheme that condenses all excess vapor before entering radiation. With our implementation, excess vapor will condense in the next timestep, but not before the condensation calculations in the macrophysics scheme. This allows our changes to affect both radiation, turbulence, entrainment, convection and detrainment (see Figure 2).
When water evaporates, vapor is released and the temperature drops. The change in vapor can be implemented in the model in different ways, and several options were tested in this study by different sets of PDand PI-simulations. The main part of this study will focus on the three setups listed in Table 1, while others will be studied and presented in section 4.1.4. In EVAP_NOMIX and EVAP_NOMIX_HIGH, the released vapor is kept inside the same gridbox as the evaporation of cloud water occurred. These setups indicate no vertical mixing between the cloudtop layer (CT) and the layer above (CT+1). In EVAP_MIX, the evaporation is assumed to either stem from mixing of air from above the cloudtop or that it initiates downdrafts at the cloudtop due to evaporative cooling, which also contributes to enhanced turbulence, entrainment, and mixing of air between the cloudtop layer and the layer above (see Figure 1). The released vapor is therefore moved one layer up. The released vapor is given as q c · evap rate , where q c is the liquid water content in the cloudtop layer. When this mass moves, the air around also moves. Due to conservation of mass, the same amount of air moves down from the layer above the cloudtop and mixes into the cloudtop layer. This is taken into account in EVAP_MIX. The mixed variables are mixing ratio of water vapor (q v , [kg kg −1 ]) and potential temperature ( , [K]). The amounts are calculated as follows in equations (4)-7.
CF is the cloud fraction in the layer the evaporation occurred, m v is the mass of water vapor (kg m −3 ), a is the air density (kg m −3 ) and is a factor converting from absolute temperature, T, to potential temperature, . = The mixing of cloudy air with drier air from outside the edges of a cloud can take place in different ways, resulting in different impact on the cloud droplet number concentration after the following evaporation (Baker et al., 1980;Latham & Reed, 1977). The two extremes are commonly referred to as homogeneous and extreme inhomogeneous mixing Lehmann et al., 2009;Pinsky et al., 2015). Homogeneous mixing occurs on a much shorter time scale than evaporation. The drier air has time to be mixed with all the cloudy air, which results in some evaporation of all droplets. In this scenario, none of the droplets evaporate completely, keeping the cloud droplet number concentration constant. This is illustrated in the For the homogeneous mixing case (left column), all droplets evaporate to some extent, but the number stays the same. For the extreme inhomogeneous mixing case (right column), some droplets evaporate completely, while the others remain unchanged. Exteme inhomogeneous mixing are also referred to as heterogeneous mixing in this paper. left column of Figure 5. On the other hand, extreme inhomogeneous mixing will evaporate some droplets completely, keeping the sizes of the remaining droplets the same due to the much slower mixing than evaporation. This is illustrated in the right column of Figure 5. A number of observational based studies have documented that both types of mixing, and also a mix between the two, can occur in liquid clouds (Freud et al., 2011;Korolev et al., 2015;Painemal & Zuidema, 2011). It is not well understood when and under which conditions one or the other will take place, although there are some studies linking the mixing type to the relative humidity above the cloud top, finding more homogeneous mixing when the air is humid (Lu et al., 2011;Yeom et al., 2017). We use a 50% mix between the two, but results from simulations with purely homogeneous and purely extreme inhomogeneous mixing (from now referred to as heterogeneous mixing) are also shown in section 4.1.4.

Size-Dependent Entrainment
The entrainment rate at the top of stratiform clouds in the model is calculated within the moist turbulence scheme, based on Bretherton and Park (2009), which again is based on the entrainment closure of Nicholls and Turton (1986). Its strength is, among others, given by an entrainment efficiency factor, A, which is affected by the evaporative cooling at the cloudtop. A is given by equation (8), where L v is the latent heat of vaporization, q l is the liquid water content at the cloudtop, and Δs is the jump of mean liquid static energy across the entrainment interface. a 2 is a tuning parameter that is commonly used to tune global models with a CAM-based atmospheric component to retrieve radiative balance and a realistic cloud deck. Bretherton and Park (2009) uses a 2 = 15, based on observations by Stevens and Feingold (2009) and Caldwell et al. (2005), while Nicholls and Turton (1986) bases their choice of a 2 = 60 on aircraft measurements. Bretherton and Park (2009) suggest tuning of a 2 within the range of 10-100. The default value of a 2 in CAM5.3-Oslo is 30. In this study, size-dependent entrainment is implemented in the model by letting a 2 vary with the effective radius of the cloud droplets in the cloudtop layer. Since the model is lacking size-dependent evaporation-entrainment feedback, scaling a 2 with the cloud droplet size is an easy, first attempt to implement this effect. As a first approach, the values in equation (9) are applied. a 2 = 100, r e < 4 μm This is also seen in Figure 6a. The r −1 relation is applied for the same reason as for the evaporation case. The horizontal distribution of a 2 and A from the PD-simulation is seen in Figures 6b and 6c.

Size-Dependent Evaporation 4.1.1. PD-Changes in LWP
Before looking at changes in LWP between PI and PD caused by differences in the additional size-dependent evaporation, we first examine how LWP in the PD-simulations of EVAP_NOMIX, EVAP_NOMIX_HIGH, and EVAP_MIX (see Table 1) differ from that of the original default model setup, N (Figure 7a). Results from EVAP_NOMIX in Figure 7b shows that evaporating up to 10% hr −1 of the liquid in the cloudtop layer has a negligible effect on LWP when keeping all changes inside the gridbox of evaporation. In EVAP_NOMIX_HIGH, the maximum evaporation rate is tuned up. Figure 7c shows that this results in a change in global mean LWP of −3.7 g m −2 (−6.9%). Since the simulations are carried out with the use of nudging, the signal vary little from year to year, with a standard error of 0.1 g m −2 for LWP. When allowing the evaporated liquid and the surrounding air masses to be mixed with the layer above, LWP changes as Description of the setups are found in Table 1. shown in Figure 7d. The global mean value in Figure 7d is negative, but there are areas with large positive values, especially outside the west coast of Africa.

How Can LW P PD Increase With Increased Evaporation?
When adding evaporation, reduced LWP is expected. The result from EVAP_MIX in Figure 7d shows that the opposite can occur. A closer look at the mechanisms behind this signal indicates that this is linked to changes in shallow convection. Figure 8 shows that the frequency of occurence of shallow convection is enhanced in many of the same areas as LW P PD increase in EVAP_MIX. Enhanced shallow convection can be caused by reduced CIN (convective inhibition), which can occur when the gradient in virtual potential temperature is affected when the evaporation changes dry the cloudtop layer and moisten the layer above . Figure S1 in the supplementary shows that this drying and moistening occur in the simulations of EVAP_MIX, especially in the regions where the vapor differences between the cloudtop layer and the layer above is large. Figure S2a in the supporting information shows reduced CIN in the same areas. When shallow convection is enhanced, more liquid cloud water is detrained to the stratiform clouds ( Figure S3).
While cloud top evaporation seems to enhance shallow convection through its impact on the stability, the Pearson's correlation coefficient in Figure 8b is only weakly positive, indicating other mechanisms involved. Figure S2b shows that the estimated inversion strength (EIS), which is a good indicator for the stratiform cloud cover (Wood & Bretherton, 2006), decreases in the same areas. While unstable regimes favor convective clouds, stable regimes favor stratiform cloud formation, trapping the moisture under the inversion. When the inversion strength is reduced, the cloud can break up, and LWP is reduced. The combination of reduced LWP due to evaporation, enhanced LWP due to shallow convection, and reduced LWP due to reduced EIS, give us a complex picture of the LWP response to enhanced evaporation of small droplets in EVAP_MIX seen in Figure 7d. Figure 9 shows how the total aerosol indirect effect (AIE PD−PI ) and the LWP change between PI and PD are affected when size-dependent evaporation is added. The left column shows results from the default model setup. The middle column shows results from the simulation EVAP_NOMIX_HIGH (see Table 1). Figure 9b shows a negligible impact on AIE PD-PI , even though we look at the case with high additional evaporation rate. Since the implementation of this effect resulted in a reduction in LWP of 6.9% (see Figure 7c) in the PD-simulation, and a corresponding change in radiative balance on the top of the atmosphere of 1.09 W m −2 (not shown here), the lack of result on AIE PD-PI is not due to weak evaporation. It is also not caused by PD-PI differences in compensating processes due to for example residual condensation or feedbacks on other processes, such as convection and mixing, since the evaporation modifications are implemented right before the radiation calculations in the code. The lack of a result stems from the r −1 -relation giving a small difference in the amount of evaporated liquid between PI and PD, resulting in the small change in ΔLW P PD-PI seen in Figure 9e. Even though the global mean evaporation rate for the additional evaporation in the cloudtop layer of the PD-simulation of EVAP_NOMIX_HIGH is as high as 51% hr −1 , it is only 1.4% hr −1 lower in the corresponding PI-simulation. As described in section 3.1, the r −1 -relation stems from surface area to volume relation of a collection of cloud droplets, and the assumption of the strength of the evaporation being proportional to the total cloud droplet surface area. Going away from the r −1 -relation means going away from this assumption.

Changes in Aerosol Indirect Effects due to Evaporation
When also including vertical mixing between the cloudtop layer and the layer above after the evaporation (case EVAP_MIX), AIE PD-PI gets less negative (Figure 9c) , which was our expected change when implementing size-dependent evaporation. This occurs even though we only allow for up to 10% hr −1 additional evaporation rate for the smaller droplets. But, when looking at Figure 9f, which shows how the change in LWP between PI and PD (ΔLW P PD-PI ) differs compared to the default model simulations, N, it is clear that the change in AIE PD-PI does not correspond everywhere with the change in ΔLW P PD-PI , at least not south of 30 • N. The hypothesis before this study was that the higher aerosol concentration in PD would evaporate more liquid in PD compared to that in PI due to the larger total cloud droplet surface area, dampening the positive signal in ΔLW P PD-PI in Figure 9d caused by suppressed autoconversion, giving negative patterns in Figure 9f. This could have been the reason for the signal north of 30 • N, where ΔLW P PD-PI is negative and AIE PD-PI is positive, but we question this link since EVAP_NOMIX showed almost equal dampening in the change in LW P PD-PI , without any impact on AIE PD-PI . Since analysis in section 4.1.2 showed that adding evaporation to the PD-simulation EVAP_MIX did not result in reduced LWP everywhere, but also increased LWP some places due to enhanced shallow convection, positive signals in Figure 9f are expected. Positive patterns in Figure 9f means that LWP changes more between PI and PD compared to that of the default setup. This should result in more negative AIE PD−PI in these regions than what seen in Figure 9 Figure 9c must stem from changes in the susceptibility of the clouds to aerosol perturbations when adding extra evaporation, rather than effects from the size dependency. Section 4.1.2 showed that adding more evaporation at the cloud tops facilitates shallow convective clouds rather than stratiform clouds in some areas. The former is associated with reduced horizontal extent, promoting a reduction in the total cloud fraction, giving a global mean absolute change in total cloud fraction of −1.8% from that of the default model setup in PD. Clouds with reduced horizontal extent no longer have the same potential to impact the radiative balance when imposing an aerosol perturbation through their reduced cloud weighted susceptibility (Alterskjaer et al., 2012;Karset et al., 2018). With the same imposed aerosol perturbation in the different setups, this change in cloud weighted susceptibility could explain the dampening of the total aerosol indirect effect seen in Figure 9c.

Sensitivity Tests: Evaporation
Five sets of simulations are carried out to study how sensitive our results are to the treatment of exchange of mass between the cloudtop layer and the layer above followed by the evaporation, and to the type of mixing. Table 2 gives an overview of the simulations. It also includes the three simulations discussed earlier in this paper for comparison. When it comes to the exchange of mass between the layers, the similar results between EVAP_MIX and EVAP_MIX3, where moist air is mixed between the layers, clearly shows that the previously discussed signals in ΔAIE PD-PI and ΔLW P PD come when not only mixing the released vapor and temperature (EVAP_MIX2) or the released vapor, temperature and dry air (EVAP_MIX4), but also the moist air (EVAP_MIX3). Assuming purely homogeneous mixing slightly reduces the impact on ΔLW P PD , while purely heterogeneous mixing has the opposite effect. This could be due to their different impact on the cloud droplet size, resulting in different autoconversion rates. Homogeneous mixing results in smaller droplets, which suppresses the autoconversion rate, contributing to a less negative signal in ΔLW P PD for EVAP_NOMIX_HIGH_HOMO compared to the other setups. In our study, the mixing assumption does not affect ΔAIE PD-PI , but this is probably due to the small difference in evaporation factor between PI and PD. With larger differences in the evaporation between PI and PD, using purely heterogeneous mixing would probably enhance the signal due to the same reasons as for ΔLW P PD . Figure 10 shows how the total aerosol indirect effect and the change in LWP between PI and PD are affected when adding size dependency on the entrainment rate at the top of the clouds. As seen in Figure 6a, the entrainment is more efficient for the smaller droplets in PD compared to those in PI (Figure 3b), resulting in a less negative indirect effect in the case named ENTRAIN compared to that of the default model setup (ΔAIE PD-PI = 0.09 W m −2 , 8.4%), as seen in Figure 10. Since the simulations are carried out with the use of nudging, the standard error is low (0.01 W m −2 ). Parts of this result are caused by a reduction in the increase in LWP between PI and PD (ΔLW P PD-PI = −0.19 g m −2 , 7.3%). While enhanced evaporation resulted in various complex responses in LWP, the signal due to enhanced entrainment is more as expected, with reduced LWP in all areas with enhanced entrainment (not shown here). As for the result in EVAP_MIX, dampened indirect effect could be due to changes in the susceptibility rather than due to the difference in cloud droplet size between PI and PD, but one of the sensitivity tests in the following section shows that this is not the case.

Sensitivity Tests: Entrainment
Four sets of simulations are carried out to study how sensitive our results are to the entrainment modifications. Table 3 gives an overview of the simulations. It also includes the simulations discussed in the previous section for comparison. As seen in Figure 6b, the global mean value of a 2 in ENTRAIN_PD is 44.8, which is higher than that of the default model setup (30.0). A higher value of a 2 means stronger entrainment efficiency, which eventually results in a reduction in the cloud water content. This can change how susceptible clouds are to aerosol perturbations, meaning that the result seen in Figure 10 could be a result of this rather than size dependency. To check if this is the case, the sensitivity test ENTRAIN_CONST is carried out applying the same constant high a 2 -value of 44.8 in both PI and PD. These simulations only gave a change in ΔAIE PD-PI of +0.02 W m −2 compared to that of the default model setup, indicating that most of the signal in Figure 10a stems from the size dependency rather than susceptibility changes.
Even though some observations suggest a 2 -values ranging from 10-100 (see section 3.2), large uncertainties are still associated with this parameter. In ENTRAIN_HIGHLIM, larger values of a 2 are allowed. The results in Table 3 show that this increases ΔAIE PD-PI slightly, but not much, probably due to its similar high impact in PI and PD. Going away from the r −1 -relation could create larger differences between the eras. This means allowing for evaporation or entrainment differences caused by other reasons than the surface area difference. One such reason could be that smaller droplets often are associated with more liquid water due to dampened autoconversion, thus even larger total surface area. In ENTRAIN_HIGHSENS, the sensitivity of cloud droplet size is increased, and the relation r −1.5 is applied. The results from ENTRAIN_HIGHSENS show almost the same ΔAIE PD-PI as ENTRAIN_PD, probably due to the reduction in global mean value of a 2 in both PI and PD. By combining higher sensitivity to cloud droplet size and higher upper limit if a 2 , ENTRAIN_HIGHSENSLIM gives the largest change in ΔAIE PD-PI among our setups, but even with the doubling of the upper limit of a 2 , and the much higher sensitivity of the entrainment rate to cloud droplet size, ENTRAIN_HIGHSENSLIM does not give more than a moderate impact on ΔAIE PD-PI of 0.13 W m −2 .

Size Dependency Applied to the Holuhraun Eruption
The 2014-2015 fissure eruption in Holuhraun, Iceland, ejected large amounts of the aerosol precursor gas sulfur dioxide (SO 2 ) into the troposphere . When SO 2 is oxidized to sulfate (SO 4 ), the SO 4 particles can act as cloud condensation nuclei, increasing the cloud droplet number concentration. Satellite retrievals show that this resulted in smaller cloud droplets during the eruption, thus a first aerosol indirect effect as expected (Malavelle et al., 2017;McCoy & Hartmann, 2015). Second aerosol indirect effects were not observed since cloud amounts and liquid water path did not change significantly (Malavelle et al., 2017). Simulations by four global models, including CAM5.3-Oslo, replicated the effect on the sizes of the cloud droplets when modeling the eruption in the study of Malavelle et al. (2017), but they failed in modeling the second indirect effect, reporting large increases in LWP that were not clearly observed by satellite. We want to check if adding size dependency to the model code could make the models agree better with observations. Since size differences had a small impact on the different amounts of evaporated liquid between small and large droplets, we focus on implementing our size-dependent entrainment.

Including Size-Dependent Entrainment
Two simulations were carried out with the default model setup, with and without the eruption. The simulations without the eruption (NOHOLU_N) starts on January 1st 2013, and last for two years. Concerning the nudging, the model-produced meteorology applied for the PD-PI aerosol indirect effect is replaced by meteorology from ERA-interim (Berrisford et al., 2009). The simulation including the eruption (HOLU_N) restarts from NOHOLU_N on 1 January 2014, but include the emissions from the eruption from 31 August 2014, as described in Malavelle et al. (2017). Similar simulations are carried out, including size-dependent entrainment, resulting in the two cases HOLU_ENTRAIN and NOHOLU_ENTRAIN.
The resulting impact on the October 2014 changes in LWP are found in Figure 11. While the change in LWP, averaged over the domain 44-80 • N, 60 • W to 30 • E, was 6.3 g m −2 (8.2%) in the simulations with the default model setup, including size-dependent entrainment dampens the signal only slightly to 5.8 g m −2 (7.7%). As for the PD-PI indirect effect, adding size-dependent entrainment described in section 3.2 is not very efficient in dampening the LWP-response to aerosol perturbations, even though we are using the outer most recommended values for the tuning parameter a 2 to enhance the entrainment efficiency for the smaller droplets. Our implementation of size-dependent entrainment cannot alone make the simulated change in LWP correspond with the observed lack of positive change (Malavelle et al., 2017). This result can indicate that the effect of size-dependent entrainment on the second aerosol indirect effect in this case study is small, but it can also just be a result of our model not being able to represent this process properly, for example due to too crude vertical resolution.

Including Other Processes
When it comes to simulating aerosol-cloud interactions, other weaknesses of global models, beside the lack of size dependency, are already known. Since the impact of adding size-dependent entrainment could differ in more updated model versions due to these weaknesses, additional sensitivity experiments, using the Holuhraun eruption as a testbed, are carried out. One of the known weaknesses is the lack of indirect effects in clouds other than the stratiform clouds (Kirkevåg et al., 2018). Another is the too strong sensitivity of the autoconversion rate to cloud droplet number concentration (Gettelman et al., 2013;Malavelle et al., 2017;Quaas et al., 2009;Wang et al., 2012). We carry out three sets of simulations, with and without the eruption, modeling the impact adjustments of these weaknesses can have on the aerosol indirect forcing alone (SCONVAIE and AUTOLOW), and together combined with the inclusion of size-dependent entrainment (ENTRAIN_SCONVAIE_AUTOLOW). In the setup SCONVAIE, detrained droplets from shallow convection are not all set to 10 μm, as in the default model setup, but they differ in size based on the background aerosol number concentration. If the column integrated aerosol number concentration, N a , is larger than 9 · 10 12 m −2 , the mean volume radius of the droplets, r vol is set to 8 μm. If N a is smaller than 3 · 10 12 m −2 , r vol is set to 12 μm, while a linear relationship is applied for concentrations in between. The N a -limits are chosen to optimize the effect of the size dependency, without the need of increasing or decreasing the r vol -limits to unrealistic values. When more clouds are susceptible to aerosol perturbations, it is likely that both the impact on the radiative forcing due to the eruption, and the effect of adding size-dependent entrainment is enhanced. In the setup AUTOLOW, we reduce the sensitivity of autoconversion rate to cloud droplet number concentration. Like many other models, CAM5.3-Oslo applies the autoconversion scheme of Khairoutdinov and Kogan (2000). This scheme is given by Equation (10), where q r and q c are the mixing ratios of rain and cloud liquid water, E is a constant representing subgrid variability, k, and are constants and N c is the cloud droplet number concentration.
In the default scheme, k is 1,350, is 2.47 and is 1.79. In the simulations with the setup AUTOLOW, and are changed to 1.67 and 0.67, respectively, based on a suggested relation between the two in the study of Feingold et al. (2013). k is tuned to 9.5 · 10 −2 so that the LWP in the simulation without the eruption is similar to the LWP in NOHOLU_N.
Allowing for indirect effects also for the detrained water from shallow convective clouds only increases the signal slightly from 6.3 g m −2 (8.2%) in the default model setup to 6.7 g m −2 (8.8%). Dampening the dependency on N d and q c on the autoconversion rate almost halves the signal to 3.6 g m −2 (4.6%). The combination of all three cases, SCONVAIE_AUTOLOW_ENTRAIN, results in a change in LWP caused by the Holuhraun eruption of 2.9 g m −2 (3.8%). Based on the individual cases, the large dampening from the default model setup is dominated by AUTOLOW, with a smaller contribution from size-dependent entrainment. Figures of the horizontal distribution of LWP changes caused by the eruption with the different setups are found in Figure S4 in the supplementary information.

Summary and Conclusions
We have explored different implementations of size-dependent evaporation and entrainment in the global atmospheric model CAM5.3-Oslo. The impacts on the total aerosol indirect effect and the second indirect effect, using the change in LWP caused by an aerosol perturbation as a measure, are studied. Our main findings are the following: We also show that the stability changes can have the opposite effect on LWP through altering EIS. (section 4.1.2).
While previous studies have highlighted the lack of buffering effects when modeling aerosol-cloud interactions with global models as a big reason why models overestimate the LWP-response to aerosol perturbations (Malavelle et al., 2017;Stevens & Feingold, 2009;Wang et al., 2003;Zhou & Penner, 2017), our simulations show a more complex picture of the cloud response when including these mechanisms. When previous studies have mentioned the lack of processes that buffer cloud responses to aerosol perturbations, enhanced evaporation of smaller cloud droplets are associated with a reduction in cloud water. Our study highlight that the response in LWP is more complex due to competing effects acting in opposite directions. Even though enhanced evaporation at the top of the stratiform clouds reduces their liquid content and extent, the following response by shallow convective clouds can act in the opposite way due to reduced stability (reduced CIN). On the other hand, reduced stability (reduced EIS) in regions dominated by stratiform clouds captured under inversions are associated with reduced cloudiness. This could mean that global mean values of AIE PD-PI and ΔLW P PD-PI may be results of compensating effects in different regions, especially with buffering processes as clouds are more or less shallow convective of sensitive to EIS. These mechanisms may not be captured by regional LES or cloud resolving models, which usually do not cover large areas with different cloud regimes and may not capture the interaction between them, since the modeled domain usually is nested one way with large scale dynamics and boundary conditions from reanalysis data or a coarser model. Our results indicate some buffering by entrainment effects, but not comparable to the magnitude of the potential of dampening the response in LWP to aerosol perturbations through altering the autoconversion process. With the large amounts of various sensitivity simulations giving small impact on aerosol indirect effects, our study indicates that discrepancies between global model estimates of LWP responses to aerosol perturbations are dominated by overestimations in the dependency of autoconversion to cloud droplet number rather than the lack of size dependency in the evaporation or the entrainment process.