Development of Two-Moment Cloud Microphysics for Liquid and Ice within the NASA Goddard Earth Observing System Model ( GEOS-5 )

This work presents the development of a two-moment cloud microphysics scheme within the version 5 of the NASA Goddard Earth Observing System (GEOS-5). The scheme includes the implementation of a comprehensive stratiform microphysics module, a new cloud coverage scheme that allows ice supersaturation and a new microphysics module embedded within the moist convection parameterization of GEOS-5. Comprehensive physically-based descriptions of ice nucleation, 5 including homogeneous and heterogeneous freezing, and liquid droplet activation are implemented to describe the formation of cloud particles in stratiform clouds and convective cumulus. The effect of preexisting ice crystals on the formation of cirrus clouds is also accounted for. A new parameterization of the subgrid scale vertical velocity distribution accounting for turbulence and gravity wave motion is developed. The implementation of the new microphysics significantly improves the 10 representation of liquid water and ice in GEOS-5. Evaluation of the model shows agreement of the simulated droplet and ice crystal effective and volumetric radius with satellite retrievals and in situ observations. The simulated global distribution of supersaturation is also in agreement with observations. It was found that when using the new microphysics the fraction of condensate that remains as liquid follows a sigmoidal increase with temperature which differs from the linear increase assumed 15 in most models and is in better agreement with available observations. The performance of the new microphysics in reproducing the observed total cloud fraction, longwave and shortwave cloud forcing, and total precipitation is similar to the operational version of GEOS-5 and in agreement with satellite retrievals. However the new microphysics tends to underestimate the coverage of persistent low level stratocumulus. Sensitivity studies showed that the simulated cloud properties are robust to 20 1 https://ntrs.nasa.gov/search.jsp?R=20140013016 2018-09-04T17:47:57+00:00Z


Introduction
Cloud microphysical schemes in global circulation models (GCMs) have evolved from directly prescribing cloud properties (i.e., particle size and number, cloud amount and Published by Copernicus Publications on behalf of the European Geosciences Union.D. Barahona et al.: Cloud microphysics in GEOS-5 concentration of condensate) (e.g., Manabe et al., 1965) to explicit representation of the microphysics involving the formation, evolution, and removal of cloud droplets and ice crystals (e.g., Gettelman et al., 2010;Lohmann, 2008;Sud et al., 2013;Quaas et al., 2009).The development of sophisticated cloud schemes allows a more realistic description of the variability and interdependence of cloud properties, and will likely improve model predictions of climate (Lohmann and Feichter, 2005).However, their increased complexity has also brought about new challenges in the description of small-scale dynamics, cloud particle nucleation and growth, and the generation of precipitation.Most models rely on simplified representations of such processes.
Current GCMs typically use either single-(e.g., Del Genio et al., 1996;Bacmeister et al., 1999) or two-moment cloud microphysics schemes (e.g., Gettelman et al., 2010;Sud et al., 2013;Lohmann et al., 2008).More detailed schemes have also been developed; however, their computational costs make them unsuitable for climate studies (Khain et al., 2000).The advantage of two-and higher-moment schemes is that the characteristics of the cloud particle size distribution are explicitly calculated and allowed to interact with radiation and influence the evolution of cloud properties.Some schemes also allow for supersaturation with respect to the ice phase, required to model ice nucleation explicitly (e.g., Gettelman et al., 2010;Wang and Penner, 2010).When coupled to an appropriate aerosol activation parameterization, multi-moment microphysics schemes are capable of modeling the modification of cloud properties by aerosol emissions, a key component of anthropogenic climate change (IPCC, 2007;Lohmann and Feichter, 2005; Kaufman and Koren, 2006).
Mounting evidence suggests that aerosols, both natural and anthropogenic, play a key role in many atmospheric processes.For example, the presence of ice in clouds at temperatures above 235 K depends on the presence of waterinsoluble ice nuclei (IN) (Pruppacher and Klett, 1997).IN in turn act as precipitation-forming agents in convective systems and mixed-phase clouds (Ramanathan et al., 2001;Rosenfeld and Woodley, 2000).Although IN originate mostly from natural sources (i.e., dust and biogenic material), anthropogenic emissions can modify the natural IN concentration.The effect of aerosols on clouds has also been associated with planetary radiative perturbations from the modification of clouds by anthropogenic aerosol emissions (Twomey, 1977(Twomey, , 1991;;Lohmann and Feichter, 2005).Emissions of cloud condensation nuclei (CCN) may also lead to the modification of the precipitation onset in convective cumulus, decreasing the average size of cloud droplets (Rosenfeld et al., 2008).Recent studies suggest that the interplay between CNN and IN plays a significant role in the maintenance of Arctic clouds (Morrison et al., 2012;Lance et al., 2011).Accurate representation of these effects in atmospheric models is critical for reliable climate prediction, yet difficult due to their complexity and gaps in the understanding of CCN and IN activation and the interactions between clouds and radiation (Stevens and Feingold, 2009).
A recent simulation of the non-hydrostatic implementation of the NASA Goddard Earth Observing System at 14 km spatial resolution demonstrated that as the spatial resolution increases, the parameterized convective transport of moisture plays a weaker role in the generation of cloud condensate.At high resolution, the simulated cloud properties are controlled by the cloud microphysics (Putman and Suarez, 2011).For typical GCM resolutions (∼ 2 • ), the parameterization of convective precipitation strongly impacts the simulation of the hydrological cycle and the distribution of cloud tracers in the atmosphere (e.g., Arakawa, 2004).Most GCMs use singlemoment schemes to describe the microphysics of convective systems; two-moment microphysical schemes have also been developed for convective clouds, mostly based on ideas originally developed for stratiform clouds (e.g., Lohmann, 2008;Song and Zhang, 2011;Sud et al., 2013).
The NASA Goddard Earth Observing System, version 5 (GEOS-5) is a system of models integrated using the Earth System Modeling Framework (ESMF) (Rienecker et al., 2008).The operational version of GEOS-5 is regularly used for decadal predictions of climate, field campaign support, satellite data assimilation, weather forecasts and basic research (Rienecker et al., 2008(Rienecker et al., , 2011;;Molod, 2012).GEOS-5 uses a single-moment cloud microphysics scheme to parameterize condensation, sublimation, evaporation, autoconversion and sedimentation of liquid and ice (Bacmeister et al., 2006).This single-moment approach captures the main climatic features related to the formation of stratocumulus decks and tropical storms (Reale et al., 2009;Putman and Suarez, 2011), but prevents the explicit linkage of aerosol emissions to cloud properties and omits subgrid variability in cloud properties.In this work, we develop a new microphysical package for GEOS-5 that addresses these issues.The new microphysics scheme presented here explicitly predicts the mass and number of cloud ice and liquid, and links the number concentration of ice crystals and cloud droplets to aerosol through the processes of cloud droplet activation and ice crystal nucleation.

Operational GEOS-5
The cloud scheme in the operational version of GEOS-5 considers a single phase of condensate; however, the removal and evaporation of cloud water from detrained convection and in situ condensation are treated separately.The fraction of condensate existing as ice is assumed to increase linearly between 273 and 235 K. Processes of autoconversion, evaporation/sublimation and accretion of cloud water and ice are treated explicitly (Bacmeister et al., 2006).Moist convection is parameterized using the relaxed Arakawa-Schubert (RAS) scheme (Moorthi and Suarez, 1992).Generation and evaporation of convective, anvil and stratiform precipitation are parameterized according to Bacmeister et al. (2006).Longwave radiative interactions with cloud water, water vapor, carbon dioxide, ozone, N 2 O and methane are treated following Chou and Suarez (1994).The Chou et al. (1992) scheme is used to describe shortwave absorption by water vapor, ozone, carbon dioxide, oxygen, cloud water, and aerosols and scattering by cloud particles and aerosols.Cloud particle effective size is prescribed and tuned to adjust the radiative balance at the top of the atmosphere.The current version of GEOS-5 also accounts for the radiative effect of precipitating rain and snow according to Molod et al. (2012).Aerosol transport is calculated interactively using the Goddard Chemistry, Aerosol, Radiation, and Transport model, GOCART (Colarco et al., 2010), a global aerosol transport model that considers dust, sea salt, black and organic carbon, and sulfate aerosols.Scavenging of aerosol mass is based on a convective mass flux approach; however, it is not explicitly linked to droplet and ice crystal nucleation (Colarco et al., 2010).
The calculation of large-scale condensation and cloud coverage in GEOS-5 follows a total-water probability distribution function (PDF) approach (Smith, 1990;Rienecker et al., 2008;Molod, 2012).The total-water PDF is assumed to follow a uniform distribution characterized by the critical relative humidity, based on the formulation of Slingo (1987).Anvil cloud fraction is parameterized following Tiedtke (1993).

New cloud variables
The cloud microphysical scheme in GEOS-5 was augmented to consider the evolution of the mass and number of ice crystals and cloud droplets.Four new prognostic variables were added to GEOS-5: q l , q i , n l and n i , representing the gridaverage mass and number mixing ratio of liquid and ice, respectively.The evolution of a given tracer, η, is described by where the terms on the right-hand side of Eq. ( 1) represent the tendency in η due to advective and turbulent transport and large-scale and convective cloud processes, respectively.Advective transport and turbulent transport in GEOS-5 are described in Rienecker et al. (2008).∂η ∂t ls refers to the change in η from non-convective cloud processes (i.e., anvil and stratus clouds), whereas ∂η ∂t cv describes the change in η from processes occurring within convective cumulus.

Microphysics of stratiform and anvil clouds
In GEOS-5, clouds are classified as stratiform (cirrus, anvils and stratocumulus) and convective.Stratiform clouds are formed by in situ condensation and anvil detrainment.The stratiform scheme of Morrison and Gettelman (2008, hereafter MG08) was incorporated into GEOS-5 as part of the new cloud scheme.Since MG08 allows for ice supersaturation, and accounts for activation of aerosols based on a subgrid vertical velocity, other aspects of the cloud scheme were updated.The calculation of cloud fraction and large-scale condensation was modified to account for supersaturation with respect to ice and microphysical processing (Sect.2.3.1).The new scheme uses the CCN activation and ice nucleation parameterizations of Fountoukis and Nenes (2005) and Barahona and Nenes (2009b), respectively (Sects.2.3.2 and 2.3.3).A parameterization of subgrid vertical velocity, w sub , was also developed (Sect.2.3.4), and MG08 was modified to account for the effect of preexisting ice crystals on cirrus formation (Sect.2.3.5).A new microphysical scheme for convective clouds that explicitly considers CCN and IN activation was implemented (Sect.2.4).These modifications represent a complete overhaul of the cloud microphysics of GEOS-5.
The MG08 scheme includes prognostic equations for the mass and number mixing ratio of cloud ice and liquid, and diagnostically predicts the vertical profiles of rain and snow.The detailed mass and number balances leading to ∂n l ∂t ls , ∂q l ∂t ls , ∂q i ∂t ls and ∂n i ∂t ls are presented in Morrison and Gettelman (2008).The MG08 scheme is used to describe the microphysics of convective detrainment and stratiform condensate.
The size distribution of cloud droplets, rain, ice and snow is assumed to follow a gamma distribution; i.e., where the subscript y is used to represent a hydrometeor species, and N 0,y and λ 0,y are the intercept and slope parameters of n y (D), calculated as in Morrison and Gettelman (2008) (cf. Eq. 3).D y and µ y are the sphere-equivalent diameter and the size dispersion of the y species, respectively.A Marshall-Palmer distribution (Marshall and Palmer, 1948) is assumed for rain and snow; i.e., µ y = 0.The version of MG08 implemented in GEOS-5 follows closely the description of Gettelman et al. (2010), with some modifications as follows.MG08 also uses an exponential approximation of the size distribution of ice crystals; i.e., µ i = 0. Theoretical considerations however suggest that n i (D i ) in recently formed clouds is better represented by log-normal and gamma functions in which the concentration of ice crystals decreases steeply for very small sizes (Barahona and Nenes, 2008).Since this behavior cannot be reproduced using an exponential distribution, setting µ i = 0 may lead to underestimation of λ 0,i and overestimation of crystal size.This assumption is relaxed in GEOS-5, and µ i is calculated as a function of T following the correlation of Heymsfield et al. (2002), obtained from extensive measurements in cirrus clouds.Following Heymsfield et al. (2002), it is assumed 1736 D. Barahona et al.: Cloud microphysics in GEOS-5 that 0.5 < µ i < 2.5.The critical size for ice autoconversion was set to D cs = 400 µm.The sensitivity of cloud ice water to µ i and D cs is analyzed in Sect. 4.
The droplet autoconversion parameterization in MG08 (Khairoutdinov and Kogan, 2000) was replaced by the formulation of Liu et al. (2006).The latter was preferred because of its greater flexibility in representing the effect of cloud droplet dispersion on the autoconversion rate.The liquid water content exponent in Liu's parameterization was set to 2.0 (Liu et al., 2006).Following Liu et al. (2008), the cloud droplet size dispersion, µ l , was parameterized in terms of the grid-scale mean droplet mass.
Other modifications to MG08 include the calculation of the nucleated droplet number and ice crystal concentration and the parameterization of the subgrid-scale vertical velocity (Sects.2.3.2 to 2.3.4).Partitioning of total condensate accounts for the Bergeron-Findeisen process following Morrison and Gettelman (2008) and Gettelman et al. (2010).Ice and liquid cloud fraction are however not discriminated, and the total cloud fraction is calculated using the probability distribution function (PDF) of total water (Sect.2.3.1).

Stratiform condensation and cloud fraction
Cloud fraction, f c , plays a crucial role in cloud and radiative processes, and is intimately tied to the in-cloud number and mass mixing ratios.In GEOS-5 it is calculated using a totalwater PDF scheme; i.e., where P q (q) is the normalized total-water PDF in the nonconvective part of the grid cell, and f cn is the convective detrainment mass fraction.P q (q) is assumed uniform with width equal to q (Appendix A); q * is the weighted saturation mixing ratio between liquid and ice, given by where f ice is the mass fraction of ice in the total condensate, and q * l and q * i are the saturation specific humidities for liquid and ice, respectively.The total condensate is given by where q c,det is the contribution of convective detrainment to the total condensate.The term S crit in Eqs.
(3) and ( 5) is termed the critical saturation ratio.S crit controls the minimum level of supersaturation required for cloud formation within a model grid cell.Equation (3) implies that regions within the grid cell for which q t > q * S crit are covered with cloud (Appendix A).The total water in the non-convective part of the grid, q t , is calculated assuming water saturation for the convective detrainment.Note that Eqs.
(3) to ( 5) are coupled through the energy balance (not shown), and must be solved simultaneously.
In the operational version of GEOS-5, it is assumed that S crit = 1 for all conditions.In this work, the same assumption is used for mixed-phase and liquid clouds.However, for ice clouds, linking S crit to ice nucleation processes increases the minimum relative humidity required for cloud formation, allowing for supersaturation with respect to ice.S crit is thus controlled by the subgrid-scale dynamics and the aerosol properties.In cirrus clouds, S crit is calculated by the ice nucleation parameterization (Sect. 2.3.3,Eq. 13).
To make an initial estimate of f c , the width of P q (q) is prescribed and parameterized in terms of a critical relative humidity (Molod et al., 2012).This is fully diagnostic, since the width does not depend on state variables.However, the convective contribution to f c is fully prognostic and depends on the detrained mass flux parameterized using a Tiedke-style approach (Tiedke, 1993).Using this assumption, an initial estimate of f c is calculated in the form (Eq. A2) where q mx = q t + 0.5 q and q are the upper limit and the width of P q (q), respectively.Similarly, for total condensate (Eq.A3), where accounts for changes in q * due to latent heating during condensation.Equation ( 6) may lead to a reduction in f c if q t < S crit q * , even if q t > q * (i.e., the grid cell is on average supersaturated), which may lead to inconsistency between ice crystal growth and total condensate.This is resolved by assuming a proportional increase in f c with water vapor deposition onto preexisting ice crystals.Cirrus clouds thus persist in supersaturated grid cells (however, is not created) even if q t < S crit q * .Evaporation, water vapor deposition and condensation, and sedimentation processes can modify f c .Microphysical processes can also alter q t and P q (q) via the formation of precipitation.Fully prognostic schemes parameterize these effects by assuming some proportionality between changes in f c and microphysical rates (e.g., Del Genio et al., 1996;Sud and Walker, 1999;Tompkins, 2002;Kärcher and Burkhardt, 2008).Here, an alternative approach, maintaining the form of the PDF, is proposed as follows.Assuming that the totalwater PDF (i.e., anvil and stratiform) after microphysical processing follows a uniform distribution, an equation similar to Eq. ( 7), but without an explicit contribution from detrainment, can be written for the total condensate (Eq.A5).Since total water, q t , and total condensate, q c , are known after the microphysics, then a new width, q , consistent with the new state of the system, can be calculated, as detailed in Appendix A. Using q in Eq. ( 3), a new cloud fraction corrected for microphysical processing can be written in the form (Eq. A8) In practice, an initial estimate of f c (Eq. 6) is used to calculate q c and q t at the beginning of the time step.Then assuming that microphysical processes proceed at a constant cloud fraction, q c and q t are calculated and introduced into Eq.( 8) to calculate f c .This procedure has the limitation that microphysical processes are calculated using an initial estimate of f c instead of its final value; however ensures consistency between f c and q c at the end of the time step.

Cloud droplet activation
CCN activation into cloud droplets is parameterized following the approach of Fountoukis and Nenes (2005) (hereafter FN05).FN05 give an analytical solution of the equations of an ascending cloudy parcel using the method of population splitting (Nenes and Seinfeld, 2003).Sulfates, hydrophilic organics and sea salt are considered CCN active species.Aerosol number concentrations were derived from the predicted mass mixing ratio for each species using size distributions obtained from the literature (Table 1).Sulfate and organics are considered internally mixed, and five separate bins are used to describe dust.Aerosol composition is parameterized in terms of the hygroscopicity parameter (Petters and Kreidenweis, 2007): κ was set to 0.65, 0.2 and 1.28 for sulfate, hydrophilic organics, and sea salt, respectively.The water uptake coefficient was set to 1.0 (Raatikainen et al., 2013).In this work, the adiabatic version of the FN05 parameterization is employed.However, FN05 can readily be extended to include dust activation (Kumar et al., 2009b), entrainment (Barahona and Nenes, 2007), and giant CCN (Barahona et al., 2010b).The contribution of CCN activation in stratiform clouds to the droplet number concentration is given by where N l and N l,act are the in-cloud preexisting and activated droplet number concentrations, respectively.N l,act is calculated at w+ sub = w + 0.8σ w (Peng et al., 2005;Fountoukis and Nenes, 2005), w and σ w being the mean and standard deviation of the subgrid distribution of vertical velocity (Sect.2.3.4), and w+ sub the vertical velocity averaged over the positive side of the distribution.This approximation is valid for w σ w , and may introduce up to 20 % non-systematic discrepancy in N l,act when compared to the direct solution of the integral in Eq. ( 14) (Morales and Nenes, 2010); however, it is justified for its computational efficiency.

Ice nucleation
The ice nucleation parameterization implemented in GEOS-5 was developed by Barahona andNenes (2008, 2009a, b) (hereafter BN09), and is summarized in Barahona et al. (2010a).The parameterization of BN09 is derived from the analytical solution of the governing equations of an ascending cloud parcel, and considers the dependency of the ice crystal concentration on cloud formation conditions, subgridscale dynamics, and aerosol properties.At cirrus levels (T < 235 K), both homogeneous and heterogeneous ice nucleation and their competition are considered; i.e., where N s i,nuc is the ice crystal concentration nucleated in a single parcel ascent, N hom and N het the ice crystal concentrations produced by homogeneous and heterogeneous ice nucleation, respectively, and S i,max the maximum saturation ratio reached within the cloudy parcel.In BN09, S i,max is explicitly calculated, accounting for the competition between water vapor deposition onto ice crystals and supersaturation generation by expansion cooling.S i,max (hence N s i,nuc ) thus depends on dynamics, temperature and the concentration of ice nuclei; i.e., S i,max = S i,max (w sub , T , N het ) (Barahona and Nenes, 2009b).Since homogeneous freezing quickly depletes supersaturation, S i,max is limited, so that S i,max ≤ S hom , S hom being the saturation threshold for homogeneous freezing (Ren and Mackenzie, 2005;Koop et al., 2000).For T > 235 K and S i,max < S hom , only heterogeneous ice nucleation takes place.
N hom is determined by the homogeneous ice nucleation rate of sulfate solution droplets, parameterized in terms of the water activity following Koop et al. (2000).Heterogeneous ice nucleation is described through a generalized ice nucleation spectrum, N het = N het (S i , T , m 1...n ), so that N het = N het (S i,max ), with S i being the saturation ratio with respect to ice, and m 1...n the moments of the aerosol number distribution.N het depends on the aerosol composition, and in principle can have any functional form (Barahona and Nenes, 2009b).The usage of N het (S i , T , m 1...n ) also obviates the need for prescribing fixed nucleation thresholds, which may carry uncertainty (Barahona, 2012).In this work, N het is described using the formulation of Phillips et al. (2013) (hereafter Ph13), considering immersion and deposition ice nucleation on dust, black carbon, and soluble organics.In simplified form, the Ph13 spectrum can be written as where N x , D g,x , σ g,x , and sp,x are the total number concentration, the geometric mean diameter, the geometric size dispersion, and the mean particle surface area of the x aerosol  (Lance et al., 2004).D g (µm) and σ g are the geometric mean diameter and dispersion, respectively.i is the particle number fraction in mode i.The "polluted" size distribution parameters for sulfate and organics are used when the total aerosol mass exceeds 5.0 µg m −3 .
Aerosol species, respectively, and ϕ x (S i , T , sp,x ) is the number of active nucleation sites per particle (Phillips et al., 2013(Phillips et al., , 2008)).
The summation in Eq. ( 11) is carried out over five lognormal modes for dust, and single log-normal modes for black carbon and organics (Table 1).Primary biological particles are not predicted by GEOS-5 and are not considered in this work.However, on a global scale, their effect on ice cloud formation may be small (Hoose et al., 2010).Since dust and soot aerosol are typically irregular aggregates rather than spherical particles, sp,x was obtained from the mean sphere-equivalent particle volume, assuming a bulk surface area density of 10 m 2 g −1 for dust (Murray et al., 2011) and 50 m 2 g −1 for soot (Popovitcheva et al., 2008).The BN09 parameterization also allows the calculation of S crit for cirrus (Eq.6, Sect.2.3.1).According to BN09, the ice saturation ratio at which most IN freeze in a polydisperse aerosol population, S het , is given by (Barahona and Nenes, 2009b) If no IN are present, then S het approaches the saturation threshold for homogeneous freezing, S hom (Barahona and Nenes, 2009b).S het and S hom represent the minimum saturation ratio required for cloud formation by heterogeneous and homogeneous freezing, respectively.They thus have the same meaning as the critical saturation ratio of Eq. ( 6). S crit is then calculated as where f het = N het /(N hom +N het ) is the fraction of ice crystals produced by heterogeneous ice nucleation.The grid cell averaged nucleated ice crystal concentration, N i,nuc , is calculated by weighting N s i,nuc over the distribution of updrafts within each grid cell (Sect.2.3.4): where φ( w, σ 2 w ) is the normal distribution and w max = w + 4σ w (Sect.2.3.4).The latter is used as an upper limit to the integral to avoid numerical instability.Note that, for ice nucleation, using the approximation N i,nuc ≈ N i,nuc w+ sub may introduce a much larger bias than for cloud droplet activation (Sect.2.3.2).This is because the competition between homogeneous and heterogeneous nucleation introduces strong nonlinearity in N i,nuc (Barahona and Nenes, 2009a).The characteristic value of w sub for N i,nuc therefore generally differs from the average vertical velocity.PDF averaging is also applied for S crit and S i,max .
The contribution of ice nucleation in cirrus to the ice crystal number concentration is given by The factor P q (q t > S crit q * i ) accounts for the probability of finding an air mass leading to cloud formation within the grid cell.This term was proposed by Barahona and Nenes (2011) to account for the effect of prior nucleation events on current cloud formation.
For the mixed-phase regime (T > 235 K), Eq. ( 11) is applied directly by assuming saturation with respect to liquid water, to find the contribution of deposition and condensation heterogeneous nucleation to N i .In this regime, cloud droplet freezing by immersion and contact ice nucleation contribute to the ice crystal population.These are explicitly treated as follows.The tendency in N i from immersion freezing of cloud droplets is parameterized in the form where is the average cooling rate (Sect.2.3.4).N x and n s,x are the number concentration and the active site surface density for species x, respectively.The latter is calculated according to Niemand et al. (2012) for dust and Murray et al. (2012) for black carbon.
Contact ice nucleation is parameterized as the product of the collection flux of aerosol particles by the cloud droplets and the ice nucleation efficiency in contact mode.Young (1974) suggested that phoretic effects and Brownian motion are responsible for collection scavenging of ice nuclei.Baker (1991) however showed that Brownian motion is the dominant factor, although phoretic effects may be significant in deep convective clouds (Phillips et al., 2007).Considering only Brownian collection, the contribution of contact ice nucleation to the ice crystal formation tendency is written in the form where dN x dt Brw is the Brownian collection flux of the x aerosol species (Young, 1974).Consistent with laboratory studies (e.g., Fornea et al., 2009;Ladino et al., 2011), the active site density in the contact mode is assumed to be the same as for immersion freezing shifted towards a higher temperature; i.e., T cont ≈ T − 3 K.
The in-cloud contribution of ice nucleation in mixed-phase clouds to the ice crystal number concentration tendency is given by where the subscripts cont, imm, and dep refer to contact, immersion, and deposition/condensation ice nucleation, respectively.The term N d t is used to limit the nucleated ice crystal concentration to the existing concentration of cloud droplets.

Subgrid-scale dynamics
Besides information on the aerosol composition and size, parameterization of cloud droplet and ice crystal formation requires knowledge of the vertical velocity, w sub , on the spatial scale of individual parcels (typically under 100 m), which is not resolved by GEOS-5.w sub depends on radiative cooling (Morrison et al., 2005), turbulence (Golaz et al., 2010), gravity wave dynamics (e.g., Barahona and Nenes, 2011;Kärcher and Ström, 2003;Jensen et al., 2010;Joos et al., 2008) and local convection.To account for these dependencies, we employ a semi-empirical formulation as follows.
In situ measurements (e.g., Peng et al., 2005;Bacmeister et al., 1999;Conant et al., 2004) suggest that w sub is approximately normally distributed.The mean vertical velocity of the distribution is written as (Morrison et al., 2005) where w ls is the grid-scale vertical velocity, c p is the heat capacity of air, g is the acceleration of gravity, and ∂T ∂t rad is the diabatic heating due to radiative transfer.Variance in w sub for stratiform clouds results from subgrid-scale eddy motion, σ 2 w,turb , and gravity wave dynamics, σ 2 w,gw ; i.e., A first-order closure is used to diagnose σ 2 w,turb (Morrison and Gettelman, 2008): where K T is the mixing coefficient for heat (Louis et al., 1983) and l m is the mixing length.MG08 prescribed a fixed l m = 300 m.To account for the spatial variation of l m , the formulation of Blackadar (1962) is used instead; i.e., where k is the von Kármán constant, z is the altitude and λ m is the value of l m in the free troposphere (Blackadar, 1962).The latter is estimated as 10 % of the boundary layer height from the previous time step (Molod, 2012).This approach also takes into account the vertical variation of l m within the planetary boundary layer (PBL).The minimum value of σ 2 w,turb is set to 0.01 m 2 s −2 within the PBL.Small-scale gravity waves strongly affect the formation of cirrus and mixed-phase clouds (e.g., Haag and Kärcher, 2004;Jensen et al., 2010;Joos et al., 2008;Barahona and Nenes, 2011;Dean et al., 2007).In situ measurements suggest that the dynamics of the upper troposphere are characterized by the random superposition of gravity waves from different sources (e.g., Jensen and Pfister, 2004;Bacmeister et al., 1999;Sato, 1990;Herzog and Vial, 2001).Random wave superposition results in a Gaussian distribution of vertical velocities (e.g., Bacmeister et al., 1999;Barahona and Nenes, 2011).Using this, a semi-empirical parameterization for σ 2 w,gw is derived in the form (Eq. B5) www.geosci-model-dev.net/7/1733/2014/Geosci.Model Dev., 7, 1733-1766, 2014 where τ 0 is the surface stress (Lindzen, 1981), U the horizontal wind, ρ a the air density, N the Brunt-Väisälä frequency, and L c the wavelength of the highest-frequency waves in the spectrum, also referred to as the characteristic cirrus scale (here assumed to be 100 m).Equation ( 23) is obtained by relating |τ 0 | to the equivalent perturbation height at the surface.This is scaled to obtain the maximum wave amplitude at each vertical level (Joos et al., 2008;McFarlane, 1987) and then used to compute σ 2 w,gw (Barahona and Nenes, 2011).This approach parameterizes orographically generated gravity waves.In practice, both the background and the orographic surface stress are used in Eq. ( 23) to avoid underestimation of σ 2 w,gw in marine regions.The second term in brackets on the right-hand side of Eq. ( 23) limits σ w,gw to account for wave saturation and breaking (Eq.B3).The derivation of Eq. ( 23) is detailed in Appendix B. Only activation processes are modified by subgrid vertical velocity variability; i.e., φ( w, σ 2 w ) is assumed to be uncorrelated with the subgrid distribution of condensate.

Preexisting ice crystals
Ice nucleation can be inhibited by water vapor deposition onto preexisting ice crystals (i.e., ice crystals present in the grid cell from previous nucleation events).Their impact on cirrus properties may be significant at low temperatures, where ice crystals are small and have low sedimentation rates (Barahona and Nenes, 2011).The effect of preexisting crystals on ice nucleation can be parameterized by reducing the vertical velocity for ice nucleation in cirrus by a factor dependent on the preexisting ice crystal concentration and size (Eq.C5); i.e., w sub,pre = (24) where N i,pre is the preexisting ice crystal concentration, λ 0,i,pre is the slope of the size distribution of preexisting ice crystals, c is a shape factor (here assumed to be equal to 1), ρ i is the bulk density of ice, and A i , α and β are temperaturedependent parameters (Appendix D). w sub,pre represents a corrected vertical velocity accounting for the effect of preexisting ice crystals limiting expansion cooling.Water vapor deposition onto preexisting crystals acts against the increase in ice supersaturation from expansion cooling.Thus, by enhancing water vapor deposition within cloudy parcels, preexisting ice crystals tend to decrease S i,max (Eq.10), leading to a reduction in N i,nuc .To account for this, S i,max is calculated using w sub,pre instead of w sub ; i.e., S i,max = S i,max (w sub,pre , T , N het ).A similar approach was proposed by Kärcher et al. (2006), who used a numerical integration technique instead of the analytical approach presented here.The derivation of Eq. ( 24) is detailed in Appendix C. The effect of preexisting ice crystals on cirrus properties is analyzed in Sect. 4.

Microphysics of convective cumulus
While all the main features of RAS are preserved in the new scheme, the removal of condensate is reformulated to account for the effect of IN and CCN emissions on the generation of convective precipitation.RAS calculates the convective cloud condensate and mass flux at each model level by averaging over an ensemble of ascending parcels, each one lifted from the top of the PBL (Molod et al., 2012;Rienecker et al., 2008).Each ascending parcel is characterized by its detrainment level and entrainment rate (Moorthi and Suarez, 1992), and saturation adjustment is used to find the amount of condensate present in each parcel.In the current RAS implementation in GEOS-5, a single parcel detrains at each model level, so that the tendency of the tracer η due to cloud convective processes is given by where D is the detrainment rate and W the convective mass flux.In the operational GEOS-5, a prescribed fraction of condensate is assumed to precipitate from each parcel before reaching the cloud top.The remaining condensate is then linearly partitioned between ice and liquid as a function of T and detrained at the neutral buoyancy level.Each convective parcel is assumed to develop independently, and the detrained condensate from different parcels is weighted by the convective mass flux.The subscript "cp" in the following equations refers to processes occurring within each parcel.A detailed description of the GEOS-5 convective scheme is presented elsewhere (Moorthi and Suarez, 1992;Rienecker et al., 2008).The balance of a tracer, η, within a convective parcel is written as where dη dt cp is the rate of change in η from microphysical processes occurring within convective parcels, w cp is the parcel vertical velocity, λ is the per-length entrainment rate, and η is the value of η in the cloud-free environment.Detrainment of condensate is assumed to occur only at the cloud top.
The rate of change in η from microphysical processes occurring within a convective cloud parcel is given by where the subscript "source" refers to nucleation, condensation and deposition processes, "precip" to precipitation, and "freezing" to phase transformation.Equation ( 26) is integrated for each parcel from cloud base to cloud top, at which all remaining condensate detrains into the anvil; i.e., 1 = Dη.The initial condition for Eq. ( 26) depends on the tracer.At cloud base, the concentration of ice crystals and the ice mass mixing ratio are assumed to be zero, whereas the activation of cloud droplets at cloud base is explicitly considered (Sect.2.4.2).
Solution of Eq. ( 26) requires knowledge of the vertical velocity within each parcel, w cp , which is also necessary for driving the droplet activation and ice nucleation parameterizations.This is calculated by solving (Frank and Cohen, 1987) 1 2 where γ = 0.5 (Sud and Walker, 1999), T v and T v are the virtual temperature of the cloud and the environment, respectively, and q cn is the mixing ratio of total condensate in the convective parcel.The first term on the right-hand side of Eq. ( 28) represents the increase in the convective parcel's kinetic energy by buoyancy, whereas the second and third terms account for the entrainment of stagnant air into the parcel and the drag from the weight of the condensate, respectively.Equation ( 28) is forward-integrated from the level below cloud base to cloud top using w cp,in = 0.8 m s −1 as an initial condition (Guo et al., 2008;Gregory, 2001); the vertical profile w cp is not very sensitive to this assumption (Sud and Walker, 1999).Note that w cp,in differs from the vertical velocity used for cloud droplet activation.The latter depends on the local buoyancy; i.e., w cp,cloudbase = w cp,in + dw cp dz z base , where z base is the model layer thickness at the cloud base.

Partitioning of convective condensate
Total condensate is partitioned between liquid and ice as follows.Nucleated ice crystals are assumed to grow by accretion of water vapor in an environment saturated with respect to liquid water.That is, the coexistence of liquid water favors a high concentration of water vapor available for deposition onto the ice crystals, and the ice and liquid phases remain in quasi-equilibrium within the convective parcel.Hydrometeor species are assumed to follow a gamma distribution (Eq.2).The growth rate of ice crystals within convective cumulus is given by (Pruppacher and Klett, 1997;Korolev and Mazin, 2003) where dq cn dt is the rate of generation of total condensate calculated by the convective parameterization, c is a shape factor (assumed equal to 1), ρ i the bulk density of ice, A i is a temperature-dependent growth factor (Appendix D), n i,cp is the ice crystal concentration within the convective parcel, λ 0,i is the slope parameter of the ice size distribution within the convective parcel, and S i,wsat is the value of S i at saturation with respect to liquid water.Using Eq. ( 29), and since q cn = q l + q i , the source term for liquid water within convective cumulus is given by where dq l dt cond is the rate of generation of liquid water within convective parcels.

Droplet activation and ice crystal nucleation in convective cumulus
Explicit activation of CCN into cloud droplets is only considered at cloud base and used as an initial condition to Eq. ( 26) (Sect.2.4).Entrained aerosols (sulfate, sea salt, and organics) above cloud base are assumed to activate instantaneously as they enter the cloud parcel.Dust and soot IN lead to the heterogeneous freezing of cloud droplets in the immersion and contact modes, described using Eqs.( 16) and ( 17).Since soot and dust particles would likely adsorb water within convective parcels (Wiacek et al., 2010;Kumar et al., 2009a), ice nucleation in the deposition mode within convective cumulus is not considered.Cloud droplets freeze homogeneously at 235 K. Frozen droplets rapidly quench supersaturation within convective cumulus.The homogeneous nucleation of deliquesced sulfate, which requires high supersaturation (S i ∼ 145-170 %, Koop et al., 2000), is thus not likely to occur within convective parcels.Homogeneous freezing of interstitial aerosol is therefore not considered within convective cumulus.

Generation of convective precipitation
Precipitation is generated within each convective parcel and assumed to reach the surface during each time step.The remaining condensate is then detrained into anvil clouds following Eq.( 25).Ice water in convective cumulus is likely to exist as graupel, snow and ice crystals with different size distributions and falling velocities, affecting the formation of precipitation.Following Del Genio et al. ( 2005), a simplified treatment is proposed, where total ice is partitioned between ice and snow (assumed as a single species) and graupel.The two species are differentiated by their terminal velocity.This partitioning is prescribed as a function of temperature and used to calculate the formation of ice precipitation within convective clouds.For ice crystal growth and detrainment a single ice species is assumed.Droplet-to-rain autoconversion is calculated according to Liu et al. (2006), and all autoconverted water is assumed to be lost as surface precipitation within one time step.The size dispersion of the droplet population, µ l , follows the formulation of Liu et al. (2008).Evaporation of convective precipitation is parameterized according to Bacmeister et al. (1999).Total ice water within convective parcels is assumed to partition between ice/snow (taken as a a single species) and graupel, and differentiated by their terminal velocity (Table 2).The fraction of total ice existing as graupel is approximated by (Del Genio et al., 2005) The particle sizes of ice/snow and graupel are assumed to follow an exponential distribution (µ g = µ i/s = 0.0) (McFarquhar and Heymsfield, 1997).The number precipitation rate of ice/snow within convective parcels is given by the number flux across the critical size for the ice/snow species, D c,i/s (Seinfeld and Pandis, 1998), where n i/s = (1 − f gr )n i,cp is the number concentration of ice/snow particles, and λ 0,i/s is the slope parameter of the ice/snow size distribution.The mass precipitation rate of ice/snow is calculated as where q i/s = (1 − f gr )q i,cp is the mixing ratio of ice/snow within the convective parcel, and + 6λ 0,i/s D c,i/s + 6] is the ratio of the volume to number fraction above D c,i/s in the size distribution of ice/snow.The term ξ i/s is introduced to account for the preferential precipitation of the largest particles of the population, which tends to enhance the mass over the number precipitation rate.The critical size for precipitation, D c,i/s , is calculated by equating the hydrometeor terminal velocity, w term , to w cp (Table 2).
Equations ( 32) and (33) assume that ice and snow grow mainly by diffusion within the convective parcel.The same assumption cannot be applied to graupel since it also grows by collection of cloud droplets.The precipitation rate of graupel is therefore calculated by removing the fraction of the size distribution above the graupel critical size, D c,g , at each model level (Ferrier, 1994) dn gr dt precip,cp = where n gr = f gr n i,cp is the graupel number mixing ratio, λ 0,g the slope parameter of the graupel size distribution and t L = z w−1 cv is the time spent by the parcel in a given model layer.Similarly for q gr , dq gr dt precip,cp = (35) where q gr = f gr q i,cp is the graupel mass mixing ratio.

Model evaluation
Model evaluation is carried out by comparing cloud properties against satellite retrievals and in situ observations.Satellite data sets included level 3 products from the NASA MODIS (http://modis.gsfc.nasa.gov/)combined TERRA and AQUA data product (Platnick et al., 2003), and the International Satellite Cloud Climatology Project (ISCCP) (Rossow and Schiffer, 1999) and CloudSat (Li et al., 2012(Li et al., , 2014) ) projects.Level 3 MODIS monthly output for the years 2003-2009 was used.CloudSat data spanned over the years 2007-2008 and a climatology for the years 1983-2008 was used for ISCCP data (Rossow and Schiffer, 1999).When possible, the Cloud Feedback Model Intercomparison Project Observation Simulator Package, COSP (Bodas-Salcedo et al., 2011), was used to compare model output against satellite retrievals.COSP uses the model output to simulate the retrieval of satellite platforms, minimizing in this way errors from the sampling of the model output when comparing against satellite observations.Global cloud radiative properties were obtained from the Clouds and Earth's Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) level 4 data product (Loeb et al., 2009) and the NASA Earth Radiation Budget Experiment (ERBE; Barkstrom, 1984).Total precipitation was obtained from the Global Precipitation Climatology Project (GPCP) data set (Huffman et al., 1997) and the CPC merged analysis of precipitation (CMAP) (Xie and Arkin, 1997).CERES, GPCP and CMAP data were available over the entire period of simulation.A climatology spanning the years 1985-2003 was used for the ERBE data.Runs were performed for a period of 10 years starting on 1 January 2001, with a spin-up time of one year using a c48 cubed sphere grid (about ∼ 2 • spatial resolution) and 72 vertical levels.Sensitivity studies (Sect.4) were performed by running the model for two years at the same resolution.Test runs showed that 19.3 0.37 Locatelli and Hobbs (1974) two years were enough to elucidate the first-order effect of variation in microphysical parameters on cloud properties.All simulations were forced with observed sea surface temperatures (Reynolds et al., 2002).Initial conditions were obtained from the Modern-Era Retrospective Analysis for Research and Applications (MERRA; Rienecker et al., 2011).
The aerosol concentration was calculated interactively using the GOCART model (Colarco et al., 2010), with emissions as described in Diehl et al. (2012).Results obtained with the operational version of GEOS-5 and using the new microphysics are referred to as the CTL and NEW runs, respectively.

Cloud fraction
The parameterization of cloud fraction in GEOS-5 was modified to account for the effect of microphysical processing on P q (q) (Sect.2.3.1) and allow supersaturation with respect to the ice phase.Figure 1 shows the effect of these modifications on the low (CLDLO), middle (CLDMD), and high (CLDHI) cloud fraction in GEOS-5.In general the CTL and NEW simulations present similar distributions of cloud fraction.However, in NEW, f c tends to be higher and in better agreement with ISCCP retrievals.The new cloud fraction scheme resulted in higher CLDLO in the remote Atlantic and Pacific oceans and reduced the cloud bias over South America and Asia.CLDLO associated with the low-level stratocumulus decks on the western coasts of North and South America and South Africa is still underpredicted in the NEW simulation.This feature is common in climate models (Kay et al., 2012); in GEOS-5 it is likely caused by the absence of an explicit shallow cumulus parameterization.The overprediction of CLDLO at the high latitudes of the NH in CTL is also significantly reduced in the NEW simulation.Overall, the global mean bias in CLDLO is significantly lower in NEW (−3 %) than in CTL (−5 %).
The global mean bias in CLDMD is also lower in NEW (−9 %) than in CTL (−15 %).The overestimation of CLDMD at the low and middle latitudes of the Southern Hemisphere (SH) and the Northern Hemisphere (NH) in CTL is largely removed in NEW, which results from a more realistic distribution of ice water content in NEW than in CTL (Sect.3.6).The underestimation in CLDMD at the middle latitudes (∼ 30 • ) of NH is also smaller in NEW than in CTL, particularly over land.However, NEW tends to increase the overestimation in CLDMD at the high latitudes of the SH.Similarly, although the CTL and the NEW simulations present similar distributions of high-level clouds (CLDHI), in general CLDHI tends to be overestimated at the marine high latitudes and underestimated over the continents.The NEW simulation also tends to underpredict CLDHI over the Tropical Warm Pool.The global mean bias in CLDHI is about 1 and 4 % in the CTL and NEW run respectively.Biases in CLDHI and CLDMD at the high latitudes (above 60 • ) of the SH and the NH tend to be more pronounced in NEW than in CTL.Although the source of these biases is not clear, they may be related to a low value of q * (Eq.4) in mixed phase clouds.Note that ISCCP retrievals tend to be uncertain in those regions as well (Rossow and Schiffer, 1999).

Supersaturation over ice
Restricting cloud formation to S i > S crit implies that supersaturation must be built before new ice clouds can form.The term P q (q t > S crit q * i ) in Eq. ( 15) also restricts ice nucleation to supersaturated regions and reduces the nucleated ice crystal concentration and the water vapor relaxation time scale.Furthermore, MG08 allows for supersaturation within cirrus since it does not apply saturation adjustment for ice clouds.These factors lead to sustained supersaturation at cirrus levels (T < 235 K).
Cloud formation and ice crystal nucleation are controlled in part by S crit , which provides an internal link between ice nucleation, f c and q i .S crit depends on T and on the local vertical velocity at the scale of individual cloudy parcels (∼ 100 m to 1 km).S crit is also determined by the availability of IN: in general, high IN concentration leads to low S crit (Barahona and Nenes, 2009b).The global distribution of S crit for T < 235 K (Fig. 2, right panel) presents two characteristic modes, showing predominance of heterogeneous (S crit ∼ 120 %) and homogeneous (S crit ∼ 150 %) ice nucleation.The peak at 150 % and the highest S crit values correspond to low T regions with high vertical velocities and low aerosol concentration, common around the tropopause (Fig. 2, left panel).Values of S crit as low as 105 % are also not uncommon, and are associated with high concentrations of active IN (e.g., dust).These are often located around T ∼ 230-240 K, where deposition/condensation IN are active and abundant enough to impact supersaturation (Sect.3.5).For lower T , the concentration of active IN is too low to decrease supersaturation substantially, and S crit increases towards S hom (Fig. 2, left panel).
The global mean value of S crit (∼ 144 %) is close to S hom , which would in principle indicate a strong predominance of homogeneous nucleation (Fig. 2, left panel).This however depends on whether a cloud is actually formed under those conditions.Although high values of S crit are very frequent for p < 50 hPa, most cirrus clouds form between 100 and 300 hPa (Sect.3.6), where S crit ∼ 110-130 %.At these vertical levels, S crit is relatively high (∼ 130 %) in the Southern Hemisphere, but lower in the Northern Hemisphere.Homogeneous freezing would thus tend to be more predominant in the Southern Hemisphere.This behavior is further analyzed in Sect.3.5.
The distribution of clear sky saturation ratio, S i,c = (q v − f c q * )/(1.0−fc ), is shown in Fig. 3. In-cloud S i is assumed to be 100 %.In reality, supersaturation relaxation may be slow in cirrus clouds, particularly at low T (Krämer et al., 2009;Barahona and Nenes, 2011).However, it is expected that for p > 200 hPa most supersaturation is relaxed inside clouds over the time step of the simulation (∼ 1800 s) (Barahona and Nenes, 2008).Figure 3 also shows data from the AIRS (Atmospheric Infrared Sounder) (Gettelman et al., 2006) and MOZAIC (Measurement of ozone and water vapor by Airbus in-service aircraft) (Gierens et al., 1999) projects.The uncertainty in the retrieval increases with S i,c .However, both MOZAIC and AIRS data show an exponential decrease in Figure 3. Global frequency distribution of clear sky saturation ratio with respect to ice obtained from 6 h instantaneous GEOS-5 output over a 3-year subset (2002)(2003)(2004) of the NEW run (left panel, black dots).Filled areas correspond to the frequency distributions from AIRS (solid area) satellite retrievals (Gettelman et al., 2006) and the MOZAIC (hatched area) data set (Gierens et al., 1999), respectively, for the years 2002-2004.Uncertainty in the observations was calculated as one standard deviation around the mean value within a 2 • × 2 • grid cell and introducing a 10 % perturbation in S i along the x axis.The center and right panels show the zonal mean frequency (%) of clear sky supersaturation from GEOS-5 and AIRS, respectively.

P (S i,c
) with increasing S i,c (Fig. 3, left panel).GEOS-5 also shows this exponential decrease and is in agreement with AIRS and MOZAIC data.The peak P (S i,c ) in the model is shifted towards S i,c ∼ 100 % since retrievals tend to avoid zones with S i,c ∼ 100 % near the cloud edges (Gettelman and Kinnison, 2007).The frequency of S i,c > 101 % in GEOS-5 distributes almost symmetrically around the tropics (Fig. 3, middle panel), with a slightly higher probability of supersaturation in SH than in NH.This is in part due to lower IN concentrations in SH (Fig. 7), although differences in the dynamics of SH and NH also play a significant role.In agreement with AIRS data (Fig. 3, right panel), GEOS-5 predicts about 10 % supersaturation frequency in the upper tropical levels.GEOS-5 seems to slightly overpredict P (S i,c > 100 %) above 300 hpa at the high latitudes of the NH and SH and near the TTL; however, the uncertainty in the retrieval in these regions is also high (Gettelman and Kinnison, 2007).

Subgrid-scale vertical velocity
The nucleation of ice crystals and cloud droplets is strongly influenced by the subgrid-scale vertical velocity, w sub .φ( w, σ 2 w ) in stratocumulus and anvils is mainly determined by σ w , whereas w is typically small (∼ 10 −2 m s −1 ).For convective clouds w cp is explicitly calculated by solving Eq. ( 28).In general the eddy contribution to σ 2 w is significant near the surface and negligible above 500 hPa.At 900 hPa, where mostly liquid clouds are formed, σ w ranges between 0.1 and 0.7 m s −1 and is typically lower over the ocean than over land (Fig. 4).High σ w is however found in the storm track regions of the Southern and Northern hemispheres.At this vertical level σ w is the lowest in the Arctic region (∼ 0.1 m s −1 ).The range of σ w shown in Fig. 4 is in good agreement with in situ measurements of vertical velocity at   38).Data for latitudes higher than 60 • have been excluded from the analysis.cloud base in marine stratocumulus (Peng et al., 2005;Guo et al., 2008), and continental regions (Fountoukis et al., 2007;Tonttila et al., 2011), showing σ w mostly between 0.2 and 1 m s −1 .However, global measurements of σ w have not been reported.Compared to similar schemes (e.g., Golaz et al., 2010) Eq. ( 21) results in higher velocities within the PBL since the characteristic length decreases near the surface, consistent with the vertical momentum balance within the PBL (Blackadar, 1962).σ 2 w thus rarely hits the prescribed minimum (∼ 0.01 m s −1 ) within the PBL.
Gravity wave motion dominates the global distribution of σ w at 500 and 150 hPa, being typically larger over land than over the ocean (Fig. 4).Air flowing over orographic features produces high-frequency waves that propagate to the free troposphere (Bacmeister et al., 1999;Herzog and Vial, 2001).σ w is thus highest over the mountain ranges of Asia, South America, and the Antarctic.At 500 hPa, σ w is about 0.1 m s −1 over land, and may reach up to 0.5 m s −1 over mountain ranges.These values are in good agreement with in situ measurements (Gayet et al., 2004).A similar distribution of σ w is found at 150 hPa, with values over land slightly higher than at 500 hPa.Over the ocean, σ w is typically larger at 150 hPa than at 500 hPa, particularly over the tropics, since gravity waves in these regions can reach larger amplitudes before breaking.Figure 4 shows that σ w in the upper troposphere varies by up to three orders of magnitude around the globe.Such variability has important implications for the effects of IN emissions on cloud formation (Sect.3.5).

Cloud droplet number concentration
Comparison of cloud droplet number concentration against satellite retrievals is typically challenging.Retrieval algorithms generally introduce assumptions on the droplet size distribution that may bias the cloud droplet number concentration.To compare satellite retrievals and model data over the same basis, we take advantage of the output generated by the COSP MODIS simulator to obtain a "model retrieved" column-integrated droplet concentration, N l,cum , in the form (Han et al., 1998) where τ is the liquid cloud optical depth, b = 0.193 (Han et al., 1998), and R eff,liq is the effective radius of cloud droplets.To apply Eq. ( 38), R eff,liq and τ are obtained either from the GEOS-5 COSP output or from the MODIS retrieval.This procedure does not aim to produce an accurate retrieval of N l,cum , but rather to compare GEOS-5 and MODIS data equally.Equation ( 38) is applied between 60 • S and 60 • N, where the MODIS retrieval is more reliable (Platnick et al., 2003).
Figure 5 shows the global distribution of N l,cum from GEOS-5 (NEW run, FN05) and MODIS.GEOS-5 is able to capture the high N l,cum found in regions of high sulfate emissions i.e., Europe, Central and Southeast Asia and the eastern coast of North America.There is also agreement between MODIS and GEOS-5 in regions with high biomass burning emissions like Subsaharan Africa and South America.However, the model tends to slightly underpredict N l,cum in the remote Atlantic and Pacific Oceans.There is also underprediction of N l,cum off the western coasts of North and South America and Africa.This is due to underprediction of shallow stratocumulus in GEOS-5 (Fig. 1) and because w sub tends to be small in these regions (Fig. 4).The global mean N l,cum in GEOS-5 (1.68 cm −2 ) is slightly lower than with MODIS results (1.96 cm −2 ).
Droplet concentration is influenced by the CCN activation parameterization and the aerosol size distribution.The GO-CART model uses a single moment aerosol microphysics, and some uncertainty may result from assuming a fixed size distribution to obtain the aerosol number concentration.The impact of this assumption is discussed in Sect. 5.The sensitivity of N l,cum to the CCN activation parameterization was studied by implementing the Abdul-Razzak and Ghan (2000) activation parameterization (Fig. 5, middle plot) and is analyzed in Sect. 4.

Ice crystal number concentration
At any given T , N i varies by up to four orders of magnitude, although mostly within a factor of 10 (Fig. 6a).The mean N i peaks around 200 L −1 at 225 K, decreasing to ∼ 20 L −1 at 190 K, and below ∼ 1 L −1 at 180 K.For T > 245 K N i remains mostly below ∼ 10 L −1 .Global mean N i is around 66 L −1 for all clouds and around 166 L −1 for cirrus (T < 235 K). Figure 6 shows agreement of GEOS-5 values with in situ measurements of N i over the whole T interval (Krämer et al., 2009;Gultepe and Isaac, 1996).There is good agreement of GEOS-5 with field campaign data at T < 200 K, where most models show a large positive bias (e.g., Barahona et al., 2010a;Salzmann et al., 2010;Gettelman et al., 2012).This results from the proper consideration of the effect of prior nucleation events on ice crystal nucleation (Section 2.3.3).N i is also influenced by the presence of preexisting ice crystals; their effect is analyzed in Sect. 4.
The relative contribution of different mechanisms to the source of N i is shown in Fig. 6.To facilitate comparison against in situ measurements, integrated variables, instead of number tendencies, are used.Thus, the ice crystal concentration from ice nucleation in the deposition and condensation modes, N dep , is calculated using Eq. ( 11) and the BN09 parameterization.N i from immersion freezing, N imm , is calculated by integration of Eq. ( 16) over the model time step.The concentration of detrained ice crystals, N cnv , is given by the ice crystal concentration at the cloud top calculated by Eq. ( 26).
N dep varies mostly within the range from 0.1 to 50 L −1 , and is largest around 240 K, where the aerosol concentration is large enough to result in significant IN concentration (Fig. 6b).There is however large variability in N dep around the globe.Most deposition IN come from dust, although the concentration of black carbon IN may be significant, reaching 2 L −1 at T ∼ 230 K (not shown).A few deposition IN (∼ 1 L −1 ) are found at T as high as 260 K, mostly in regions of large dust concentration.
N imm reaches up to 40 L −1 around 240 K, but decreases rapidly for lower T , where it is prevented by the homogeneous freezing of cloud droplets (Fig. 6c).In agreement with in situ observations of mixed-phase clouds (e.g., DeMott N cnv remains below 50 L −1 for T > 240 K, characteristic of heterogeneous ice nucleation.For T > 250 K, N cnv reaches up to 10 L −1 mostly from immersion and contact freezing of supercooled droplets within the convective cumulus (Fig. 6d).Homogeneous freezing of cloud droplets is evident in the strong increase in N cnv around T ∼ 240 K, which in some instances may reach up to 10 cm −3 .Such very high N ncv is responsible for the highest values of N i in Fig. 6.Along with immersion freezing, detrainment from convective cumulus determines N i for T > 240 K.
Figure 7 (left panel) shows the spatial distribution of ice crystal concentration nucleated in cirrus (T < 235 K) and weighted by cloud fraction.The spatial distribution (also  weighted by cloud fraction) and zonal mean of the contribution of heterogeneous ice nucleation to N i,nuc are shown in the middle and right panels of Fig. 7, respectively.Globally, about 70 % of the production of ice crystals in cirrus proceeds by homogeneous freezing, with a clear contrast between the Northern Hemisphere (NH) and the Southern Hemisphere (SH).Homogeneous freezing is most prevalent in the SH, and only on the western coasts of South America and Africa is the contribution of heterogeneous freezing significant (∼ 30 %; Fig. 7, middle panel).By contrast, most of the NH is influenced by IN emissions, which in some cases dominate the ice crystal production.
Part of the higher predominance of heterogeneous ice nucleation in NH than SH is explained by the greater abundance of dust in NH.However, comparison of Figs. 4 (right planel) and 7 (left panel) also reveals a marked effect of σ w on N i .Low σ w tends to enhance the effect of IN on N i because of the greater residence time of the heterogeneously frozen ice crystals in the parcel before the onset of homogeneous freezing, and the lower rate of increase of supersaturation Figure 9. Zonal mean non-convective liquid water mass mixing ratio (mg kg −1 ) (upper panels) and total liquid condensate (water and rain, bottom panels) for non-convective clouds from the CTL and NEW runs and the CloudSat retrieval (Li et al., 2014).Model results span over 10 years of simulation, whereas CloudSat retrievals are plotted for the period 2007 to 2008.(Barahona and Nenes, 2009a).Heterogeneous freezing thus tends to dominate ice crystal production in regions of low σ w and low N i,nuc like Sub-Saharan Africa, the Arctic, and the western coast of North America, even though these regions are not characterized by high emission rates of IN (Fig. 7,left panel).This result is also consistent with the study of Cziczo et al. (2013) who found predominance of heterogeneous ice nucleation in these regions.However, Fig. 7 (middle and right panels) shows that in most other regions, and globally, homogeneous ice nucleation tends to dominate the global production of ice crystals.This suggests that variability in σ w plays a significant role in defining the effect of IN emissions on cirrus formation.

Cloud liquid and ice water
The implementation of the new microphysics resulted in significant improvement of the representation of ice and liquid water content in GEOS-5.Figure 8 shows the zonal mean ice mass mixing ratio, q i , from the NEW and CTL simulation compared to the CloudSat retrieval for non-convective, non-precipitating ice (Li et al., 2012).The global distribution of q i in the NEW simulation is in better agreement with the satellite retrieval than that obtained in CTL.The excessive freezing around T = 240 K, characterized by the bulls-eye pattern around 600 hPa in the CTL run, is not present in the NEW simulation.In absolute terms, q i in the NEW and CTL runs is generally lower than CloudSat data, although mostly within the intrinsic error of the retrieval, about a factor of 2 (Li et al., 2012;Eliasson et al., 2011).Including snow in the comparison (Fig. 8, bottom panels) still results in lower ice and snow concentration than in CloudSat, although within the error of the retrieval.
Figure 9 shows the zonal mean liquid mass mixing ratio, q l , from GEOS-5 for the CTL and NEW runs compared against the CloudSat retrieval for non-convective, nonprecipitating liquid water (Li et al., 2014).There is far lower q l in the NEW than in the CTL run, particularly over the tropics and the subtropics of the NH.Above 900 hPa, the spatial distribution of q l in the NEW run is in better agreement than CTL.In absolute terms q l in NEW is closer to Cloud-Sat than in CTL.However, this must be taken with caution as CloudSat may not retrieve liquid water close to the ground (Devasthale and Thomas, 2012).The NEW and CTL simulations however show that most liquid water is held below the 850 hPa level in GEOS-5.The bottom panels of Fig. 9 also suggest that the rain mass mixing ratio is lower in NEW than in the CTL simulation and CloudSat.Still, the spatial distribution of the concentration of liquid and rain from NEW and from the CloudSat retrieval shows similar characteristics.
The spatial distribution of the liquid water path (LWP) (Fig. 10) in the NEW simulation is similar to that observed by CloudSat, although in general LWP is larger in the NEW simulation that in CloudSat, particularly over marine regions.Comparison against other retrievals reveals uncertainty in experimental observations of LWP.Annual average LWP from MODIS is 144 g m −2 , about twice as much as the GEOS-5 output when using COSP to simulate the MODIS retrieval (60 g m −2 ).MODIS however tends to predict higher LWP in polar regions than in the tropics, pointing to an artifact of the retrieval (Platnick et al., 2003).SSMI data (Spencer et al., 1989) is also typically used for model evaluation, although it is restricted to oceanic regions.Annual mean LWP from SSMI is about 84 g m −2 , which is higher than predicted by GEOS-5 over the ocean (∼ 48 g m −2 , not shown).
Figure 10 shows the annual mean IWP (non-precipitating, non-convective) from GEOS-5 and CloudSat (Li et al., 2012).In general there is reasonable agreement in IWP between CloudSat and GEOS-5, with and slightly higher ) for non-convective, non-precipitating clouds from GEOS-5 output using the new microphysics and from the CloudSat retrieval (Li et al., 2012(Li et al., , 2014)).
IWP in GEOS-5 (27.1 g m −2 , NEW run) than in Cloud-Sat (25.8 g m −2 ).There is also uncertainty in IWP obtained by different retrievals; however, a recent intercomparison showed agreement between the ISCCP and CloudSat retrieved IWP (Eliasson et al., 2011).GEOS-5 is able to capture the high IWP observed in the Tropical Warm Pool, Central Asia, and over the mountain ranges of Africa, and North and South America.The high IWP of the latter regions results in part from strong ice crystal production over mountain ranges (Sect.3.5).GEOS-5 however underestimates IWP in the tropical western Pacific Ocean.The spatial distribution of the total-water path (liquid and ice) is similar to that obtained with CloudSat, although the global mean TWP is higher in GEOS-5 (∼ 64 g m −2 ) than in the retrieval (∼ 49 g m −2 ) due to the larger LWP in GEOS-5.

Supercooled cloud fraction
Figure 11 shows the supercooled cloud fraction (e.g., the fraction of cloud condensate present as liquid, SCF = 1−f ice , in mixed-phase clouds for the CTL and NEW simulations.In the CTL simulation the total condensate is linearly partitioned into liquid and ice between 235 and 270 K (Bacmeister et al., 2006).In the NEW simulation partitioning of the condensate is carried out taking into account the activity and concentration of IN and the Bergeron-Findeisen process.In CTL most values of SCF below 260 K follow the prescribed linear tendency.Variability in SCF increases strongly above 260 K due to the freezing of condensate at 273 K and iceenhanced precipitation (Fig. 11).The tendency of SCF with T in NEW shows different features than in CTL following a sigmoidal instead of a linear tendency.This behavior has been observed in satellite retrievals and field campaigns (Choi et al., 2010;Hu et al., 2010) and is characteristic of immersion freezing mediated mainly by dust (e.g., Murray et al., 2011;Marcolli et al., 2007).The region of maximum SCF frequency in Fig. 11 however expands about 10 K, which results from variation in particle size and concentration, the presence of black carbon IN, enhanced precipitation in mixed-phase clouds, and variation in σ w .There is also a higher frequency of SCF > 0.4 for T < 255 K in the NEW than in the CTL simulation, which results from a higher fraction of supercooled liquid in the convective detrainment in NEW than in CTL.
Compared with CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) retrievals (Fig. 11, solid lines) (Choi et al., 2010), SCF in NEW is shifted by about 6 K towards higher T , which implies that clouds tend to glaciate at higher T in the model than observed by the satellite.This would indicate higher IN activity (i.e., a higher dust concentration or more active dust) in GEOS-5 than implied by the CALIOP data.However, CALIOP is sensitive mostly to cloud-top properties, and SCF may be biased low in deep convective clouds, where most of the supercooled liquid is found below the cloud top (Hu et al., 2010).The influence of these factors on SCF requires more investigation, and will be undertaken in a future study.The sigmoidal increase of SCF with T in both GEOS-5 and the satellite retrieval still indicates that SCF is significantly influenced by the presence of IN.

Cloud droplet and ice crystal effective radii
The annual mean droplet effective radius R eff,liq from the NEW simulation (14.3 µm) is in agreement with MODIS retrievals (14.8 µm) (Fig. 12).This is higher than the prescribed mean for the CTL run and simulated by other models also using the MG08 stratiform microphysics (∼ 9-11 µm) (Gettelman et al., 2008;Salzmann et al., 2010) but similar to the one obtained in Sud et al. (2013) in GEOS-5.The results presented in Fig. 12 benefit from using the COSP package that accounts for the preferential cloud-top sampling of MODIS (Bodas-Salcedo et al., 2011).Other studies (Gettelman et al., 2008;Salzmann et al., 2010) however did not use COSP for comparison.In agreement with the MODIS retrieval the spatial distribution of R eff,liq in the NEW run shows a clear ocean-land contrast (Fig. 12).R eff,liq is overestimated in the western coasts of South America, Africa, and to a lesser extent, North America, due to low N l over these regions.Over the land R eff,liq is underestimated in southern and central Asia, Europe and the western coast of North America, likely due to the high concentration of cloud droplets predicted by GEOS-5 in these regions (Sect.3.4).
The global distribution of ice effective radius, R eff,ice , for the NEW run is presented in Fig. 13 along with MODIS retrievals.The global mean value of R eff,ice in the NEW simulation (26.2 µm) is in good agreement with the satellite (24.2 µm).GEOS-5 is able to reproduce the low R eff,ice seen by MODIS over most of the large mountain ranges, e.g., over the Andean and Himalayan regions, although it tends to underestimate R eff,ice over northeastern Asia.Low R eff,ice is caused by strong homogeneous freezing events with N i > 1 cm −3 in high orographic uplift (Fig. 4), although local convection may also have an effect on R eff,ice as detrainment from deep convection tends to increase N i (Sect.3.5).There is some contrast in R eff,ice between land and ocean in the MODIS retrievals, which is captured by GEOS-5.However, the model tends to overestimate R eff,ice in the subtropical regions of NH and SH, which may be caused by low σ w leading to low N i .
There may be some uncertainty in the retrieval of R eff,ice , particularly for optically thick clouds (Chiriaco et al., 2007).To corroborate the GEOS-5 results further, in situ observations of the volumetric ice crystal radius, R vol,ice = , are used.Figure 14 shows R vol,ice as a function of T along with a composite of in situ data from several field campaigns (Krämer et al., 2009;McFarquhar and Heymsfield, 1997).There is reasonable agreement between the field data and the model, particularly for T < 230 K, where both show a decrease in R vol,ice with decreasing T .Around T ∼ 230 K the model tends to predict slightly higher R vol,ice than the observations, although mostly within the spread of the data.The discrepancy may also be a result of crystal shattering in ice crystal probes, which tends to increase measured N i decreasing R vol,ice (Krämer et al., 2009).The smooth   transition in R vol,ice at 235 K indicates that both homogeneous and heterogeneous ice nucleation significantly contribute to ice crystal formation at this temperature (Sect.3.5).In agreement with observations (McFarquhar and Heymsfield, 1997) R vol,ice increases steadily for T > 235 K, which results from increasing vapor deposition rates and decreasing N i as T increases (Sect.3.5).

Annual mean diagnostics
Table 4 and Fig. 15 show the summary of the annual mean cloud properties analyzed in this work.Annual mean LWP for the NEW run is 37.1, and 60 g m −2 if the MODIS COSP simulator is used.LWP in NEW is higher than the CloudSat retrieval (23.0 g m −2 ) (Li et al., 2014) mostly from higher LWP at the middle latitudes of the SH, and lower than MODIS retrievals (∼ 100 g m −2 ).Ocean-only LWP is also lower than SSMI output by about a factor of two (not shown).LWP in GEOS-5 refers only to non-convective (anvil and stratiform) clouds and it is likely that the discrepancy with SSMI and MODIS originates from the consideration of convective clouds in the retrievals.IWP in NEW (27.2 g m −2 ) is in better agreement with CloudSat (25.8 g m −2 ) (Li et al., 2012), although GEOS-5 tends to overestimate IWP at the middle latitudes of the SH and the NH.Including snow in the comparison does not affect IWP in the tropics, but results in larger subtropical IWP in NEW than in CloudSat.Global mean LWP in CTL is higher (60.0 g m −2 ) and IWP slightly lower (19.0 g m −2 ) than in NEW.
The prescribed R eff,liq and R eff,ice in CTL are generally smaller than those retrieved by MODIS with a global mean bias of about −5 and −4 µm for R eff,liq and R eff,ice , respectively.R eff,liq and R eff,ice in NEW are closer to MODIS with a global bias of about −0.5 and 2 µm, respectively (Table 4), well within the intrinsic error of the retrieval (King et al., 2003).Zonal mean R eff,liq is however overestimated in the Northern Hemisphere from underestimation of N l in oceanic regions (Sect.3.4).
Global mean cloud fraction in the NEW simulation is higher than in CTL but still lower than ISSCP retrievals (Rossow and Schiffer, 1999).The higher f c in NEW results from higher cloud coverage over continental regions (Sect.3.1).There is good agreement between NEW and IS-CCP cloud fraction at the continental middle latitudes and most of the underestimation in NEW originates in marine regions.However, in these regions both the NEW and CTL simulations show agreement with the MODIS retrieval.The reason for the better agreement of GEOS-5 with MODIS than with ISCCP in these regions is however not clear but may be related to differences in the cloud masks of ISCCP and MODIS (Pincus et al., 2012).Global annual mean precipitation, P tot , is lower in the NEW (2.72 mm d −1 ) than in the CTL (2.85 mm d −1 ) simulation and in better agreement with GPCP (Huffman et al., 1997) and CMAP (Xie and Arkin, 1997) observations (∼ 2.6 mm d −1 ), although both simulations tend to overestimate P tot in the tropics.In SH the NEW simulation tends to predict P tot higher than CMAP and lower than GPCP, whereas CTL is in better agreement with GPCP data.In NH, P tot in the NEW and CTL simulations is closer to GPCP than to CMAP data, although in NEW it tends to be lower than the GPCP observations.
The global top of the atmosphere (TOA) net radiative balance is about +0.95 W m −2 in the NEW simulation.The slight radiative imbalance in NEW results in part from the negative bias in stratocumulus cloud coverage in the NEW simulation (Sect.3.1).The liquid cloud optical depth in NEW however agrees with MODIS data (Fig. 15), particularly over the tropics.In CTL liquid clouds tend to be optically much thicker than MODIS observations (Fig. 15), which results from larger LWP and smaller R eff,liq than the observations (Sects.3.6 and 3.8).The higher optical depth in CTL leads to a more negative SWCF (−52.1 W m −2 ) than in CERES and to a higher net radiative imbalance W m −2 .Longwave cloud effect (LWCF) is similar in the CTL and NEW runs (∼ 25.0 W m −2 ) and in agreement with CERES data (26.2W m −2 ).Compared to MODIS ice cloud optical depth is however overestimated in CTL and underestimated in NEW.In NEW the low bias in ice optical depth is compensated by a positive bias in the high-level cloud fraction (Sect.3.1).

Sensitivity studies
Tables 3 and 4 present a summary of the sensitivity of GEOS-5 to different microphysical parameters.To study the sensitivity of cloud properties to the description of CCN activation, the parameterization of Abdul-Razzak and Ghan (2000) (hereafter, ARG) was implemented (run ARGACT, Table 3 and Fig. 5, middle panel).ARG is based on a fit to the numerical solution of the equations of an ascending parcel written in terms of dimensionless parameters.Compared to the NEW run, the usage of ARG resulted in slightly higher N l than with the FN05 formulation, particularly over marine regions (Fig. 5, middle panel).The ARG parameterization also predicts higher droplet concentration in regions of high aerosol emissions like Southeast Asia and southern Africa.Global mean R eff,liq was lower for ARG than for FN05 by about 0.7 µm leading to about 2 W m −2 more negative SWCF (Table 4).LWP and cloud fraction remained almost the same as in NEW suggesting that the change in SWCF was driven by modification of cloud albedo.
The sensitivity of cloud properties to the characteristic cirrus scale, L c , was also investigated.L c is associated with the wavelength of the highest-frequency waves leading to cirrus formation (Eq. 23), and impacts the subgrid vertical velocity variability in the upper troposphere.Increasing L c from 100 to 400 m reduced global N c by about a factor of two (run LC400) due to a reduction in σ w and a decrease in the rate of ice crystal formation.The global mean R eff,ice increased by about 3 µm and LWCF decreased by 2 W m −2 .Global mean σ w for L c = 400 m is about 0.07 m s −1 and 0.11 m s −1 at 500 and 150 hPa, respectively, about half the value obtained in the NEW simulation (Fig. 4).These values are still within the observed values in field campaigns (e.g., Gayet et al., 2004), and further observations are needed to better constraint L c .Table 4 however shows that GEOS-5 results are robust to moderate changes in σ w .
The effect of the dispersion in the ice crystal size distribution, µ i , on ice cloud properties (Table 4) was analyzed by setting µ i = 0.0 (run MUIZERO) instead of using a temperature dependent parameterization for µ i (Sect.2.3).This led to about a factor of two lower IWP and R eff,ice than in NEW, which resulted from an increase in autoconversion and accretion of ice by snow at low T (not shown).Despite the lower IWP, the lower ice crystal size increased the ice cloud optical depth and resulted in slightly higher LWCF and SWCF than in the NEW simulation.Because of this compensating effect the radiative properties of ice clouds are robust to moderate changes in the ice crystal size distribution.Decreasing the critical size for ice autoconversion from 400 to 200 µm (run DCS200) also increased ice autoconversion leading to lower IWP than in NEW.R eff,ice was also reduced, although to a lower extent than in MUIZERO.The net radiative effect of reducing D cs was thus a decrease of about ∼ 6 W m −2 in LWCF.
Several studies were performed to investigate the sensitivity of GEOS-5 to the description of heterogeneous ice nucleation.In NOBC and NOGLASS the effect of black carbon and glassy IN, respectively, was switched off.These runs suggested that black carbon and glassy IN only have a subtle effect on global climate (Table 4), although their local effects may be significant.In particular black carbon IN tend to increase LWCF in regions of high aerosol emissions like East Asia and the eastern coast of North America.In the same regions glassy IN tend to reduce N i at low T (Fig. 16).The global TOA radiative imbalance due to black carbon and glassy IN amounts to −0.05 and −0.18 W m −2 , respectively.Although these values are comparable to other published studies (Gettelman et al., 2012), they must be taken with caution, since they are based on limited results.A comprehensive description of the aerosol indirect effect in GEOS-5 will be addressed in future studies.
In the PDA08 run the Phillips et al. (2008) (hereafter Ph08) ice nucleation spectrum was used.Ph08 was employed in previous studies to study the effect to the ice nucleation spectrum on N i (Barahona et al., 2010a;Morales Betancourt et al., 2012;Liu et al., 2012).Ph08 accounts for the effect of both soluble and insoluble organic material acting as IN, whereas in Ph13, only soluble organics are considered to be IN.Using the Ph08 parameterization reduced N i increasing R eff,ice by about 1 µm, slightly decreasing LWCF.This resulted in part from the effect of organic IN inhibiting homogeneous freezing in cirrus clouds.Other cloud properties remained similar as in NEW.
The effect of preexisting ice crystals on ice crystal formation was analyzed in NOPREEX, where it was assumed that N i,pre = 0.For this run, the global mean N i was 359 L −1 , about twice that in NEW, with the greater increase occurring between 200 and 240 K (Fig. 16), and mostly in the tropics (not shown), indicating that the presence of ice crystals from convective detrainment tends to inhibit new ice nucleation events.Mean R eff,ice was reduced by about 6 µm, increasing LWCF by 5 W m −2 .
In NOCNV the generation of precipitation in cumulus convection was described by a single-moment approach (Bacmeister et al., 2006).Some studies (e.g., Gettelman et al., 2008;Salzmann et al., 2010) did not consider explicitly the freezing and activation of aerosol particles in convective cumulus.It is thus important to study how this assumption would affect GEOS-5 results.In NOCNV the contribution of convective detrainment to ice crystal and droplet number concentration was approximated by assuming a fixed droplet size of 10 µm for droplets and using the correlation of McFarquhar and Heymsfield (1997) to obtain the ice crystal size as a function of T .Compared to NEW, the single-moment approach resulted in enhanced precipitation rates, particularly over the Tropical Warm Pool.SWCF and LWCF were lower than in NEW by about 3 W m −2 , which was in part the result of a lower detrainment flux of condensate in the tropical upper troposphere.R eff,liq decreased by about of 1 µm due  (c,  d) Shortwave (SWCF) and longwave (LWCF) cloud forcing against CERES-EBAF retrievals (Loeb et al., 2009).(e) Liquid water path against CloudSat (black lines) and MODIS (black circles) retrievals.(f) Non-convective, non-precipitable ice water path against Cloud-Sat retrievals (Li et al., 2012(Li et al., , 2014)).Also shown is the total (ice and snow) non-convective ice water path (red circles) from GEOS-5 using the new microphysics.(g) Total cloud fraction from COSP output against MODIS (black lines) and ISCCP (black circles).(h) Total precipitation against GPCP data (Huffman et al., 1997).Also shown are data from the CMAP data set (Xie and Arkin, 1997) (black circles).(i, j) Liquid and ice optical depth (COD) from COSP output against MODIS retrievals.to an increase in droplet number concentration.Mean R eff,ice only changed by about 0.5 µm; however, N i was slightly increased, particularly at low T (Fig. 16).
Finally it is important to analyze the effect of microphysical parameters on N i at low T .Figure 16 shows the temperature dependency of N i for the runs of Table 4.All curves of Fig. 16 show the same characteristics, increasing N i with decreasing T to a maximum around 210 K and then decreasing to values typically below 10 L −1 at 185 K.The only exception to the latter is the NOCNV run in which mean N i is about 140 L −1 at 185 K, resulting from the lower detrained  4.
N i acting as preexisting ice crystals at low T .The maximum N i is around 300 L −1 for most runs, and only for the NO-PREEX run does it increase up to 800 L −1 .The fact that in all runs N i decreases for T below 200 K indicates that as the T decreases, N i becomes more dependent on S crit (Sect.3.2).This indicates that parcel history plays a primary role in determining N i at low T , whereas preexisting ice crystals and IN only play a secondary role.

Summary and conclusions
A new cloud microphysics scheme was developed for the NASA GEOS-5 global atmospheric model.The main features of the new microphysics are -A comprehensive two-moment microphysics description for stratiform clouds (Morrison and Gettelman, 2008).
-Consistent coupling of the cloud fraction and stratiform condensation with the microphysics.The stratiform condensation scheme was modified to allow supersaturation in ice clouds.
-A two-moment microphysics scheme embedded within the RAS convective parameterization.The new scheme explicitly treats the formation of droplets and ice crystals, the partitioning of condensate between ice and liquid, and the generation of precipitation within convective cumulus.
-A comprehensive description of cloud droplet activation and ice nucleation in stratiform and convective clouds, linked to the aerosol physicochemical properties.The description of ice formation considers homogeneous freezing of cloud droplets and interstitial aerosol as well as heterogeneous ice nucleation on ice nuclei.Competition between homogeneous and heterogeneous ice nucleation, and between different ice nuclei is explicitly treated.Immersion, contact, condensation and deposition ice nucleation modes are considered.
-Explicit calculation of the critical saturation ratio for ice formation considering aerosol properties, temperature and subgrid-scale dynamics.
-Explicit parameterization of the effect of preexisting ice crystals on ice nucleation.
-Explicit parameterization of the distribution of subgridscale vertical velocity in stratiform clouds, accounting for the effect turbulence and gravity wave motion on the vertical velocity variance.A new parameterization in terms of large-scale variables was developed for the latter.
The new microphysics was evaluated against satellite retrievals and field campaign data.Usage of the COSP satellite simulator greatly facilitated the comparison with satellite observations, reducing the uncertainty in the sampling of the model results.In general, cloud microphysical fields like ice water, liquid water content and droplet and ice crystal size were in much better agreement with observations than when obtained with the operational version of GEOS-5.The model performance in reproducing the observed total cloud fraction and longwave and shortwave cloud forcings is also improved, and is in reasonable agreement with satellite observations.
In the new microphysics ice and cloud droplet nucleation are tightly linked to the evolution of the cloud properties.Cloud droplet number impacts the formation of precipitation.Precipitation decreases total water, which in turn feeds back into the cloud fraction through modification of P q (q) (Sect.2.3.1).The link between N i , f c , and q i is stronger, since the production of condensate is controlled in part by S crit , which depends on the presence of IN (Eq. 13).The linkage between cloud micro-and macro-physical variables in the model emphasizes the internal consistency of the new cloud scheme.
A new cloud coverage scheme was developed to allow supersaturation with respect to the ice phase.The frequency and spatial distribution of supersaturation simulated by the model was in good agreement with satellite and in situ observations.It was shown that supersaturation is controlled in part by ice crystal nucleation and the value of S crit .The latter dictates the minimum water vapor threshold required for cloud formation.S crit is highly variable over the globe, and dependent on aerosol concentration and temperature.Models that assume a single threshold for ice cloud formation are thus inherently biased.
The variation of supercooled cloud fraction with temperature in the new microphysics followed a sigmoidal tendency.This is in agreement with CALIOP data (Choi et al., 2010) and differs from the typical linear increase of SCF with T assumed in most GCMs.There are no temperature-based constraints to the occurrence of the Bergeron-Findeisen process nor to the partition of total condensate between ice and liquid in the new microphysics.The sigmoidal tendency in SCF resulted from explicit consideration of homogeneous, immersion and contact freezing in the model.This suggests that rather than temperature alone, the presence of IN greatly influences the frequency of supercooled liquid in mixed-phase clouds.
A new approach was proposed to parameterize the distribution of subgrid-scale vertical velocity in cirrus and stratocumulus, which takes into account turbulence and gravity wave motion.Although no observational studies have been reported on the global distribution of σ w , the parameterization results were within reported values in field campaigns.Since the parameterization proposed here focuses on surface and orographic stresses, which are higher over land, σ w may be underestimated in the upper troposphere in oceanic regions.The ability to predict σ w as a function of large-scale variables still points in the right direction to reduce one of the main sources of uncertainty in the modeling of the effect of aerosol emissions on climate.It was also shown that the variability in σ w is a determining factor defining the effect of IN emissions on cirrus formation.
The simulated ice crystal concentration was in agreement with field campaign data, even at very low T , where most models tend to overestimate N i (e.g., Barahona et al., 2010a;Salzmann et al., 2010;Hendricks et al., 2011).In GEOS-5, the decrease in N i with decreasing T results from an increase in S crit (Fig. 2), which limits P q (q t > S crit q * i ) at low T , decreasing the probability of homogeneous freezing events.The term P q (q t > S crit q * i ) in Eq. ( 15) provides a link between current cloud formation and prior ice nucleation events (Barahona and Nenes, 2011).This suggests that a statistical rather than a single-parcel approach (e.g., Jensen et al., 2012;Spichtinger and Cziczo, 2010) is required for the correct modeling of low-temperature cirrus.
A new parameterization of the effect of preexisting ice crystals on ice cloud formation was developed.It was shown that their effect is more pronounced for T around 200 K, typically reducing N i .However, preexisting ice crystals alone can not explain the low ice crystal concentration at low T .The effect of organic glassy IN on cloud formation was also analyzed and it was found that it tends to reduce N i at low temperatures.Although these factors alone cannot explain the tendency of N i at T < 190 K, they are still necessary for reproducing the observed N i in the upper troposphere.In fact, it was found that the observed values of ice crystal concentration in the upper troposphere result from a combination of several factors: parcel history, IN concentration, convective detrainment and subgrid dynamics.
Effective cloud droplet size simulated with GEOS-5 was in agreement with the MODIS retrieval.There was however a slight underestimation in R eff,liq over the land, and overestimation over the tropical marine regions.This points to the need for a more sophisticated description of aerosol microphysics in GEOS-5.Sensible assumptions were made regarding the aerosol size distribution; however, there is high variability in the aerosol properties around the globe, which may affect CCN activation.The inclusion of a more comprehensive aerosol microphysics in GEOS-5 will be addressed in a future study.The simulated cloud droplet number concentration also showed some sensitivity to the parameterization of CCN activation, which in turn influences the cloud albedo.
There was good agreement in the global mean ice effective radius between GEOS-5 and the MODIS retrieval.The decrease in R vol,ice as T decreases, a common feature of in situ observations (Krämer et al., 2009) was also captured by GEOS-5.The model was able to capture key features of the spatial distribution of R eff,ice , as for example the predominance of low R eff,ice near mountain ranges.This was a result of the explicit consideration of ice nucleation and of the spatial variation of σ w,gw .R eff,ice was however overestimated in marine regions, particularly in the Southern Hemisphere.The parameterization of σ w,gw developed in this work may underestimate σ w over the ocean.Other IN sources like biological particles (Burrows et al., 2013) and sea salt (Wise et al., 2012) were not considered in this study but may enhance ice nucleation in marine environments.Some uncertainty may be introduced by the single-moment approach used for the aerosol microphysics in GEOS-5 ice nucleation, although ice nucleation is less dependent on aerosol size than CCN activation.There is also uncertainty in the formulation of the heterogeneous ice nucleation spectrum, since factors like mixing of dust/soot with sulfate, which may lead to IN deactivation/activation, are not taken into account in Ph13.The role of the uncertainty in the satellite retrieval must also be accounted for when comparing R eff,ice against model results.All of these factors require further investigation.Nevertheless, the approach proposed here results in a realistic and reasonable spatial distribution of R eff,ice .
It was shown that the cloud radiative fields modeled in GEOS-5 with new microphysics are in good agreement with observations, although local biases may be significant.GEOS-5 tends to underestimate the optical depth of persistent stratocumulus decks, which leads to a negative radiative bias in the western Pacific.Reducing such bias requires an explicit representation of shallow cumulus condensation in GEOS-5.The long-term and large-scale climatic response of GEOS-5 with the new microphysics will be analyzed in a future study.
A simple approach was assumed to describe the cloud microphysics in convective clouds.The description of precipitation within convective cores is highly complex due to the interplay of several cloud species (e.g., graupel, hail, rain, ice and snow).Some authors have developed more comprehensive microphysical packages for convective clouds, including processes of autoconversion, aggregation, collection and accretion (e.g., Song and Zhang, 2011;Sud and Walker, 1999;Lohmann, 2008).To be effective, a detailed description of microphysics in convective clouds requires prognostic prediction of the vertical profiles of rain and snow, which is not implemented in most GCMs.Also, collection and aggregation rates depend on the vertical profiles of rain and snow, which are not known in advance.The advantages of a complex representation of the microphysics of convective cores must thus be weighted against the uncertainty introduced in accommodating such descriptions within the diagnostic integration schemes of the convective parameterizations in GCMs.
The model results were quite robust to variation in microphysical parameters.The largest differences from the base configuration were found for a decrease in the size dispersion parameter of the ice crystal size distribution and in the critical size for ice autoconversion.Both changes lead to a reduction in R eff,ice and IWP and modified the long wave cloud forcing.The high sensitivity of R eff,ice and IWP to the value of µ i suggests that more attention must be put on its correct parameterization in GCMs.
The implementation of the comprehensive microphysics developed in this work resulted in a more realistic simulation of cloud properties in GEOS-5.Mounting evidence suggests that the explicit description of processes of droplet and ice crystal nucleation and precipitation is necessary for the correct representation of clouds in Earth system models.The new microphysics would likely result in improved and more realistic climate simulations in GEOS-5.The new parameterizations developed here may also help to improve our understanding of the role of microphysics and aerosol emissions on the evolution of clouds.Within the larger picture, the further development of the microphysics GEOS-5 will help to understand the role of clouds on climate and eventually reduce the uncertainty in their prediction.scale of gravity wave motion leading to the formation of clouds is typically smaller than the scale of the GCM grid cell.A spectrum of vertical velocities rather than a single wave may thus be a more realistic representation of the subgrid dynamics in the upper troposphere.Still, surface perturbations are likely to determine the maximum w sub in the spectrum (Joos et al., 2010;Barahona and Nenes, 2011).Using this concept, a semi-empirical parameterization for σ w,gw can be developed as follows.
The mean vertical momentum flux at the surface (McFarlane, 1987) is given by where δh s is the vertical displacement at the surface caused by the orographic perturbation, N s the Brunt-Väisälä frequency at the surface and U s the surface wind (taken as the geometrical mean between the meridional and zonal components), and k is the horizontal wave number.Equating τ with the mean surface stress, τ 0 , and scaling δh according to McFarlane (1987) i.e, δh = δh s ρ a U s N s /ρ a U N 1/2 , the mean wave amplitude, δh, at any height can be written as where U N is the saturation wave amplitude (Dean et al., 2007).The maximum vertical velocity in the gravity wave spectrum is related to δh by (Joos et al., 2008) w max = kU δh. (B3) In a spectrum of randomly superimposed gravity waves, w max can be empirically related to σ w,gw by (Barahona and Nenes, 2011) σ w,gw ≈ 0.133w max (B4) making k = 2π L c and combining Eqs.(B2) to (B4), we obtain σ 2 w,gw = 0.0169 min where L c is a characteristic scale associated with the wavelength of the highest-frequency waves in the spectrum, typically between 50 and 500 m (Bacmeister et al., 1999), although considered a free parameter.

Appendix C: Parameterization of the effect of preexisting crystals on ice nucleation
Water vapor deposition onto ice crystals from previous nucleation events decreases supersaturation and may reduce N i , particularly at low temperatures (Barahona and Nenes, 2011).To account for this effect, the local rate of change of S i in a cloudy parcel with preexisting crystals is written in the form (Barahona and Nenes, 2011) dS Ice crystal nucleation in cirrus occurs over small S i intervals (Barahona and Nenes, 2008;Kärcher and Lohmann, 2002).Therefore, to a good approximation, the size of preexisting ice crystals can be considered constant during ice nucleation.With this assumption, Eq. ( C3) can be reorganized as dS i dt = αw sub S i 1 − N i,pre πβcρ i A i (S hom − 1) 2λ 0,i,pre αw sub S hom (C4) where it was assumed that S i −1 S i ≈ S hom −1 S hom .If N i,pre = 0 then Eq. (C4) reduces to the saturation balance of a parcel with no preexisting crystals present (Barahona and Nenes, 2008).The effect of preexisting crystals on ice nucleation can thus be accounted for by redefining the cloud scale vertical velocity in the form w sub,pre = (C5) w sub max 1 − N i,pre πβcρ i A i (S hom − 1) 2λ 0,i,pre αw sub S hom , 0 .
Equation (C5) shows that the effect of water vapor deposition onto preexisting crystals can be understood as a reduction in the rate of increase of supersaturation by expansion cooling.Since w sub is typically an input to ice cloud formation parameterizations, Eq. (C5) also provides a simple way of accounting for the effect of preexisting ice crystals on ice cloud formation, independently of the ice nucleation parameterization employed.Graupel number concentration and mixing ratio, respectively n i/s , q i/s Ice plus snow number concentration and mixing ratio, respectively N het Ice nucleation spectrum N imm Ice crystal concentration produced by immersion freezing N 0,y Intercept parameter of n y (D) n s,x Immersion active site surface density for the x species N x Aerosol number concentration of the x species n y (D) Size distribution of the y species p Pressure P q (q) Probability distribution of total cloud condensate p s,w , p s,i Liquid water and ice saturation vapor pressure, respectively P tot Total precipitation q * Weighted saturation mixing ratio between liquid and ice q c Total condensate mixing ratio q c,det Detrained condensate mixing ratio www.geosci-model-dev.net/7/1733/2014/D. Barahona et al.: Cloud microphysics in GEOS-5

Figure 2 .
Figure 2. Annual zonal mean (left panel) and global frequency distribution (right panel) of the critical saturation ratio, S crit (%), for the cirrus regime (T < 235 K), obtained from 6 h instantaneous GEOS-5 output over a 3-year subset (2002-2004) of the NEW run.Solid bold lines (left panel) represent the annual mean tropopause pressure.

Figure 4 .
Figure 4. Annual mean subgrid vertical velocity standard deviation, σ w , for the NEW run.

Figure 6 .
Figure 6.Global frequency of in-cloud ice crystal number concentration as a function of temperature from 6 h instantaneous GEOS-5 output over a 3-year subset (2002-2004) of the NEW run.(a) Ice crystal concentration, N i .Solid lines represent the 25 and 75 % quantiles from the field campaign data analysis of Krämer et al. (2009).Solid-dotted lines represent the typical range of mean N i found in mixed-phase clouds (Gultepe and Isaac, 1996).(b) Ice crystal concentration from deposition/condensation ice nucleation, N dep .(c) Ice crystal concentration from immersion ice nucleation, N imm .(d) Ice crystal concentration from convective cumulus detrainment, N cnv .

Figure 7 .
Figure 7. Annual mean ice crystal concentration nucleated in cirrus (T < 235 K) weighted by cloud fraction for the NEW run (left panel).Also shown are the weighted average (center panel) and zonal mean (right panel) fractions of ice crystal production by homogeneous freezing in cirrus.

Figure 8 .
Figure8.Zonal mean non-convective ice water mass mixing ratio (mg kg −1 ) (upper panels) and total ice condensate (ice and snow, bottom panels) for non-convective clouds from the CTL and NEW runs and the CloudSat retrieval(Li et al., 2012).Model results span over 10 years of simulation, whereas CloudSat retrievals are plotted for the period 2007 to 2008.

Figure 11 .
Figure 11.Global frequency of supercooled cloud fraction (SCF) from GEOS-5 for the CTL and NEW runs.The most frequent SCF value for each temperature is marked * .The solid lines represent the range of SCF (mean plus and less one standard deviation) derived from the CALIOP satellite retrieval for the years 2006-2007(Choi et al., 2010).

Figure 14 .
Figure 14.Global frequency of ice volumetric radius as a function of temperature from GEOS-5, NEW run.Solid lines represent the 25 and 75 % quantiles from the field campaign analysis of Krämer et al. (2009).Filled circles were calculated using the correlation obtained by McFarquhar and Heymsfield (1997) from field measurements in mixed-phase and cirrus clouds.

Figure 15 .
Figure 15.Annual zonal means from the GEOS-5 model for the CTL (blue lines) and NEW (red lines) runs compared against different observations (black lines).(a, b) Liquid (R eff,liq ) and ice (R eff,ice ) effective radius from COSP output against MODIS.(c,d) Shortwave (SWCF) and longwave (LWCF) cloud forcing against CERES-EBAF retrievals(Loeb et al., 2009).(e) Liquid water path against CloudSat (black lines) and MODIS (black circles) retrievals.(f) Non-convective, non-precipitable ice water path against Cloud-Sat retrievals(Li et al., 2012(Li et al., , 2014)).Also shown is the total (ice and snow) non-convective ice water path (red circles) from GEOS-5 using the new microphysics.(g) Total cloud fraction from COSP output against MODIS (black lines) and ISCCP (black circles).(h) Total precipitation against GPCP data(Huffman et al., 1997).Also shown are data from the CMAP data set (Xie and Arkin, 1997) (black circles).(i, j) Liquid and ice optical depth (COD) from COSP output against MODIS retrievals.

Figure 16 .
Figure 16.Annual mean ice crystal concentration as a function of temperature for the different runs of Table4.

Table 1 .
Log-normal size distribution parameters used in this study

Table 3 .
Description of sensitivity runs performed with GEOS-5 using the new microphysics.

Table 4 .
Annual mean model results and observations.The experimental data sets are described in Sect.3.CTL and NEW refer to runs with the operational version of GEOS-5 and with the implementation of the new microphysics, respectively.Sensitivity studies are described in Table3and Sect. 4.
i dt = αw sub S i − β

Table D1 .
Continued.Ice crystal concentration nucleated in a single parcel ascent N i,nuc Ice crystal concentration nucleated in cirrus N l,act Activated cloud droplet number concentration N l,cum Column-integrated droplet number concentration n l , N l Grid mean and in-cloud droplet number concentration, respectively n i , N i Grid mean and in-cloud ice crystal number concentration, respectively n i,cp Ice crystal number concentration within convective parcels N cnv Detrained ice crystal concentration N dep Ice crystal concentration produced by deposition and condensation nucleation n gr , q gr