To model aerosol–cloud interactions in general circulation models (GCMs), a prognostic cloud scheme of cloud liquid water and amount is expanded to include droplet number concentration (Nd) in a way that allows them to be calculated using the same large-scale and convective updraft velocity field. In the scheme, the evolution of droplets fully interacts with the model meteorology. An explicit treatment of cloud condensation nuclei (CCN) activation enables the scheme to take into account the contributions to Nd of multiple aerosol species (i.e., sulfate, organic, and sea-salt aerosols) and to consider kinetic limitations of the activation process. An implementation of the prognostic scheme in the Geophysical Fluid Dynamics Laboratory (GFDL) AM2 GCM yields a vertical distribution of Nd with a characteristic maximum in the lower troposphere; this feature differs from the profile that would be obtained if Nd is diagnosed from the sulfate mass concentration based on an often-used empirical relationship. Prognosticated Nd exhibits large variations with respect to the sulfate mass concentration. The mean values are generally consistent with the empirical relationship over ocean, but show negative biases over the Northern Hemisphere midlatitude land, perhaps owing to the neglect of subgrid variations of large-scale ascents and inadequate convective sources. The prognostic scheme leads to a substantial improvement in the agreement of model-predicted present-day liquid water path (LWP) and cloud forcing with satellite measurements compared to using the empirical relationship.
The simulations with preindustrial and present-day aerosols show that the combined first and second indirect effects of anthropogenic sulfate and organic aerosols give rise to a steady-state global annual mean flux change of −1.8 W m−2, consisting of −2.0 W m−2 in shortwave and 0.2 W m−2 in longwave. The ratios of the flux changes in the Northern Hemisphere (NH) to that in Southern Hemisphere (SH) and of the flux changes over ocean to that over land are 2.9 and 0.73, respectively. These estimates are consistent with the averages of values from previous studies stated in a recent review. The model response to higher Nd alters the cloud field; LWP and total cloud amount increase by 19% and 0.6%, respectively. Largely owing to high sulfate concentrations from fossil fuel burning, the NH midlatitude land and oceans experience strong radiative cooling. So does the tropical land, which is dominated by biomass burning–derived organic aerosol. The computed annual, zonal-mean flux changes are determined to be statistically significant, exceeding the model’s natural variations in the NH low and midlatitudes and in the SH low latitudes. This study reaffirms the major role of sulfate in providing CCN for cloud formation.
Clouds, which cover a significant fraction of the earth’s surface, play a critical role in affecting the radiation balance by partly reflecting the incoming shortwave sunlight back to space and by absorbing infrared radiation emitted by the surface. As a highly dispersed system made up of tiny droplets, a water cloud can be characterized in terms of liquid water content (LWC) and droplet number concentration (Nd). The corresponding cloud optical depth (τ) is approximately proportional to LWC2/3N1/3d (Twomey 1977), indicating that the distribution of water mass among droplets has a direct bearing on the radiative properties of clouds. Such a mechanism is the origin of aerosol–cloud interactions and indirect effects of aerosols.
Soluble aerosol particles can be activated into droplets as cloud condensation nuclei (CCN) through dissolution, allowing water condensation to occur at a relatively low supersaturation. As a variety of natural and anthropogenic sources increase CCN concentrations, a higher Nd not only makes clouds more reflective by increasing cloud albedo (the first indirect effect), but also lengthens cloud lifetime by suppressing precipitation (the second indirect effect). Both effects may also lead to a weakening of the hydrological cycle (Ramanathan et al. 2001). The complexity and our relatively poor understanding of aerosol–cloud interactions contribute to the persisting uncertainties associated with the estimated radiative impacts of indirect effects (Ramaswamy et al. 2001). For example, Lohmann and Feichter (1997) reported that the first and second indirect effects together can cause a global annual mean top-of-the-atmosphere (TOA) flux change ranging from −1.4 to −4.8 W m−2, depending on the parameterizations of cloud cover and autoconversion.
Any effort to incorporate aerosol–cloud interactions into general circulation models (GCMs) involves the task of relating Nd to aerosol properties (e.g., mass concentration, size distribution, chemical composition, solubility, and surface activity) important for determining CCN efficiency. As one of the early attempts, Boucher and Lohmann (1995, hereafter BL95) proposed an empirical sulfate–Nd relationship by fitting measured/inferred Nd to sulfate mass concentrations. The relationship was adopted by a number of GCM studies of indirect effects (e.g., Lohmann and Feichter 1997; Rotstayn 1999; Kiehl et al. 2000; Ming et al. 2005), and led to remarkably similar results. Nonetheless, the large scatter in the measured Nd used in BL95 is an indication of controlling factors other than the sulfate mass concentration, implying that BL95 may oversimplify the aerosol–Nd linkage.
The more recent and physically based approaches used in Ghan et al. (1997) and Lohmann et al. (1999) explicitly accounted for the meteorological and microphysical processes involved in droplet evolution as source and sink terms in the prognostic scheme of Nd in GCMs. As the single most important source for Nd, the rates of CCN activation were calculated using the parameterizations of Abdul-Razzak and Ghan (2000) and Chuang and Penner (1995), respectively. These factored in the influence of aerosols and updraft velocity on the highly nonlinear activation process by fitting parcel model simulations. These studies assumed that most sink processes remove cloud liquid water and droplets proportionally. Hence the fractional removal rates for Nd were the same as those for liquid water, a common prognostic variable in GCMs.
In the atmosphere, cloud liquid water and droplets are produced and dissipated simultaneously, as they describe different aspects of the same process. No matter which meteorological process (e.g., large-scale ascent, radiative cooling, turbulent cooling, and convection) causes cloud formation at a particular time and location, the same driving force must be responsible for the production of droplets as well as cloud liquid water. This seemingly plain real-world observation poses a fundamental constraint of self-consistency on any prognostic scheme used for simultaneously tracking cloud liquid water and droplets (in some cases, cloud amount as well). For example, all the prognostic variables should share the same source and sink terms. If cloud water detrained from convective updrafts contributes to stratiform clouds, the droplets contained in detrainment ought to be accounted for accordingly. Also, all the variables should be driven by exactly the same processes. Thus, for example, it is inadvisable to use different updraft velocities for cloud liquid water and for droplets.
This paper proposes a self-consistent prognostic scheme for cloud liquid water, amount and droplets. The predicted geographical and vertical distributions of prognosticated Nd, obtained from an implementation of the scheme in the Geophysical Fluid Dynamics Laboratory (GFDL) AM2 GCM, are compared with those diagnosed in the same model using the BL95 empirical relationship (Ming et al. 2005). The predicted cloud parameters are compared with satellite measurements. The combined aerosol first and second indirect effects are quantified using the difference in TOA radiation flux between preindustrial (PI) and present-day (PD) conditions. Statistical significance of the difference is assessed by comparing with the model’s unforced variations.
2. Description of scheme
a. Prognostic equation for cloud droplet number concentration
An earlier version of the GFDL AM2 GCM (Ming et al. 2005) employed a variant of the cloud scheme of Tiedtke (1993) that was adapted for predicting cloud liquid water (l) and amount (a) in GCMs, while diagnosing Nd from the sulfate mass concentration. The new cloud scheme largely retains the treatment of l and a, and adds the airmass-based droplet number concentration (n) as a third prognostic variable (cloud ice is a separate prognostic tracer). Note that n is defined as the number of droplets per unit air mass (cloudy and clear together), as opposed to Nd that is defined on the basis of unit cloudy air volume (in-cloud), and is often used in reporting in situ measurements. For a grid box, n can be converted to Nd following
where ρa is air density.
By explicitly accounting for the processes that contribute to droplet formation and dissipation, we arrive at the governing equation for n:
(∂n/∂t)|ad and (∂n/∂t)|tu represent the advective and turbulent transport of droplets across the boundaries of model grid boxes. The convective and large-scale sources and sinks are denoted by (∂n/∂t)|cv and (∂n/∂t)|ls, respectively.
b. Convective source
The GFDL AM2 GCM uses the relaxed Arakawa–Schubert (RAS) scheme (Moorthi and Suarez 1992; GFDL Global Atmospheric Model Development Team 2004) to calculate the contribution of water detrainment from convective updrafts to cloud condensate. Detrainment from each updraft member j is assumed to occur only at cloud top. Consequently, one can write in analogy to the liquid source term:
where n*jtp is the droplet number concentration in the convective updraft j at cloud top if it happens to detrain at layer i (i.e., cloud top is at layer i). Here, ni is the droplet number concentration in the environment at layer i; MjI is the corresponding convective mass flux; and p and g are pressure and the gravitational constant, respectively. The mass detrainment rate, Dj, can be calculated as
where Mjtp and ptp are the mass flux and pressure at cloud top.
The assumption of no detrainment at layers beneath cloud top and conservation of droplet number leads to:
where ncv-ac,ba represents the number concentration of droplets activated at the base of a convective cloud, while njcv-ac,i is for those generated from activating CCN entrained at a particular layer i above cloud base. Also, rj represents the fraction of droplets that is removed through precipitation, and is assumed to be a function of detrainment pressure (see appendix).
An assumption of Eq. (5) is that both clear and cloudy air can be entrained with probability of (1 − a) and a, respectively, consistent with the assumptions applied to cloud liquid and amount. Equation (5) also implies that CCN that are not activated at the time of being entrained will remain unactivated all the way toward cloud top. This is a hypothesis that can be justified on the basis that experimental and modeling evidence suggests that typically the supersaturation within a cloud arrives at a maximum fairly close to the base and then decreases monotonically with height (see section 2d; Pruppacher and Klett 1997). So, immediately following the initial activation, the supersaturation decreases as an air parcel is transported vertically in a convective updraft, thus precluding further activation.
c. Stratiform source/sinks
Large-scale ascent, turbulent and radiative cooling may give rise to stratiform clouds by saturating air. Following the formulation of Tiedtke (1993), cloud droplet number shares the same large-scale sources and sinks as liquid water and amount. The term (∂n/∂t)|ls in Eq. (2) is written as
where the ∂n/∂t terms with the subscripts ls-ac, ev, er, au, ac, fr, ri, and bf refer to gain/loss in droplet number due to large-scale activation, evaporation, erosion, autoconversion, accretion, freezing, riming, and the Bergeron–Findeisen process, respectively. Note that the term evaporation specifically refers to evaporation of droplets due to large-scale descent and warming, as opposed to erosion that describes evaporation caused by horizontal mixing of saturated and unsaturated air.
Since CCN activation occurs only in the newly formed cloudy fraction of a grid box, the contribution of large-scale activation can be expressed as
where nls-ac is the droplet number concentration of new or incremental cloud amount. The time derivative of a[(∂a/∂t)|c] due to large-scale condensation is calculated following Tiedtke (1993).
This study treats the impact of large-scale descent and warming on droplets differently from that on cloud water. Parcel model simulations (not shown here) suggest that droplets dissipate not nearly as fast as cloud water. So, we assume that large-scale descent can partly evaporate cloud water, but cannot dissipate droplets unless a cloud completely evaporates as a result. We also assume that cloud erosion and all microphysical processes—namely autoconversion, accretion, freezing, riming, and the Bergeron–Findeisen process—remove cloud droplets and liquid water proportionally. Hence, an existing parameterization that is formulated for cloud water in terms of fractional removal can be applied to droplets without modification. Because the homogeneous freezing of cloud water is treated as occurring instantaneously at temperatures below −40°C, it is reasonable to assume proportional removal for the freezing process. Note that the Bergeron–Findeisen process involves evaporation of droplets. The treatment here is inconsistent with large-scale evaporation. However, sensitivity studies (not shown here) indicate that the assumption of proportional removal for the other processes including the Bergeron–Findeisen process and autoconversion causes negligible difference in the model results as compared to alternative assumptions.
d. Parameterization of activation and other microphysical processes
It is important to distinguish CCN activation in the absence of preexisting droplets from that in their presence in the case of convective clouds. Initial activation at cloud base belongs to the former case. Then, when the droplets produced at cloud base are carried upward in a convective system, they compete for water with entrained CCN in a way that suppresses supersaturation, and subsequently affects the concentration of newly formed droplets (the latter case).
Without preexisting droplets, the concentration of activated droplets is a function not only of aerosol size distribution and chemical composition, but also of updraft velocity, which varies from a few meters per second for convective clouds to as low as a fraction of a centimeter per second for stratiform clouds. The tremendous range of updraft velocity poses a serious challenge for any effort to parameterize the highly nonlinear activation process, especially at low updraft and/or high aerosol concentrations (i.e., conditions favorable for kinetic limitations). Given that, this study utilizes the parameterization proposed by Ming et al. (2006) because of its satisfactory performance over a wide range of updraft velocity to calculate ncv-ac,bc in Eq. (5) and nls-ac in Eq. (7). Because the RAS scheme does not resolve the vertical profile of updraft velocity, the average updraft velocity in an ensemble of convective clouds is approximated as 0.5CWFj, where cloud work function (CWFj) is the integral of buoyancy along the path of convection for the jth cloud type. By doing so, convective updraft velocity is linked to cloud type. Because the vertical velocity is used to calculate activation near cloud base, a factor of 0.5 is used rather than the maximum possible, 2, to approximately account for the fact that a convective updraft may accelerate along the path. The rates of turbulent and radiative coolings are converted to equivalent updraft velocities using the moist adiabatic lapse rate. These are added to the actual model-resolved grid-mean velocity of large-scale ascent outside of convective updrafts and downdrafts to form a composite updraft velocity that is used to calculate large-scale condensation as well as droplet activation.
Though most droplets in convective clouds are activated at cloud base, we still need an approach to quantifying activation above cloud base. For an entraining air parcel rising at a velocity of V, the time derivative of the supersaturation s obeys
where WL is the concentration of cloud liquid water [unit: kg(water) m−3(air)]. The coefficients α′ and γ depend on the concentration of water vapor and temperature of the cloud (WV and T, respectively) and of the environment (W ′V and T ′, respectively) as well as the entrainment rate (μ):
The physical parameters involved in Eqs. (9) and (10) are the universal gas constant R, the molecular weight of water Mw, the latent heat of water ΔHυ, the heat capacity of air Cpa, the molecular weight of air Ma, and the saturated vapor pressure of water ps.
It is a good assumption that the amount of water condensing onto droplets is significantly larger than that onto interstitial particles after the initial activation at cloud base occurs. Thus the condensation rate can be expressed as
where Dd and n are the average diameter and number concentration of droplets, respectively, seq is the equilibrium supersaturation of droplets. The growth coefficient G depends on droplet size, and ρw is water density.
A common characteristic of the profiles of supersaturation simulated with a parcel model (Ming and Russell 2004) and plotted in Fig. 1 is that, after reaching a maximum at cloud base (the bottom 10 m or so), it levels off rapidly, and keeps on decreasing at a pace as slow as 10−5 s−1, thus entering a pseudoequilibrium state and making it plausible to set ds/dt to zero numerically. This is based on an assumption of homogeneous mixing of entrained air. By inserting Eq. (11) into Eq. (8) and setting its left-hand side to zero, we can solve for the pseudoequilibrium s. It is further assumed that, among the newly entrained CCN, only those whose critical supersaturation is lower than the calculated s will be activated, so that njcv-ac,i in Eq. (5) can be calculated. Though this study suggests that the contribution to Nd of activation above cloud base, as modeled here, is negligible for shallow convection that gives rise to liquid water clouds, it may be significant for deep convection; this approach could be readily modified to account for its impact in a study of indirect effects involving ice clouds.
Cloud erosion follows the same expression as in Tiedtke (1993), and the rate constant is higher in convective and turbulent conditions (see appendix). The parameterization of the microphysical processes largely follows Rotstayn (1997) and Rotstayn et al. (2000). As the only exception, autoconversion is calculated according to the work of Khairoutdinov and Kogan (2000), which considers the continuous effect of Nd on autoconversion rate.
e. Implementation of scheme in GCM
The new cloud scheme is implemented in the GFDL AM2 GCM. A detailed overview of model structure and validation is provided by the GFDL Global Atmospheric Model Development Team (2004). The prescribed sea surface temperature (SST) simulations described in this paper are run at a horizontal resolution of N45 (2° latitude by 2.5° longitude) and with 24 vertical layers, which lie mostly in the troposphere. All the runs discussed in this study cover a duration of 6 yr with the first year as spinup. The aerosol climatology used in the GCM is prescribed by monthly mean 3-dimensional distributions of mass concentrations calculated offline with a chemistry transport model with prescribed global emissions of aerosol precursors. The aerosols thus do not interact with the model meteorology. Globally uniform assumptions of aerosol physical and chemical properties have to be made since not enough measurements are available to constrain them on the global scale. Aerosol species are assumed to be externally mixed based on field measurements (e.g., Svenningsson et al. 1992; Berg et al. 1998). All the species share the same size distribution, consisting of two lognormal modes (N1:N2 = 17:3, Dg,1 = 0.01 μm, σ1= 1.6, Dg,2 = 0.07μm, σ2 = 2.0; Whitby 1978), which is represented by 40 logarithmically spaced equal size bins in the GCM. This is a simplification as numerous chemical and physical processes may affect aerosol sizes. Sulfate and sea-salt aerosols are treated as (NH4)2SO4 and NaCl, respectively, while organic aerosol is modeled as a mixture of four soluble compounds including malic acid (48% by mass), citric acid (22%), glucose (4.8%), and fructose (4.7%), and a generic insoluble species (20.5%). This composition represents organic aerosols in both clean marine and polluted continental environment well in terms of hygroscopic growth and CCN efficiency (Ming and Russell 2001, 2004). Surface-active organic aerosol can lower droplet surface tension below that of pure water, and the relationship between surface tension and total concentration of organic species in droplets used in the parameterization of CCN activation is derived from the measurements of Facchini et al. (2000). Dust and black carbon are not assumed to be CCN here. The appendix lists the values of model parameters and describes the numerical method used in integrating the governing equation for the stratiform cloud source/sinks [Eq. (6)].
a. Column burdens
The monthly mean aerosol climatology of sulfate and organic carbon (OC) is generated using the Model for Ozone and Related Chemical Tracers (MOZART; Tie et al. 2001; Horowitz et al. 2003; Tie et al. 2005; Horowitz 2006). The total mass concentration of coarse and fine sea salt is related to Special Sensor Microwave Imager (SSM/I)-observed monthly mean surface wind speed (Haywood et al. 1999). Since fine particles account for most CCN in terms of number, this study assumes that 10% of the total concentration is from fine sea salt (O’Dowd et al. 1997). The transport of fine sea salt over land is not considered. The PI and PD geographical distributions of column burdens in January and July are plotted in Fig. 2.
Largely resulting from burning sulfur-containing fossil fuels, the highest PD burdens of sulfate aerosol center over the Northern Hemisphere (NH) midlatitude industrial regions, namely East Asia, Europe, and United States. Long-range transport of continental emissions considerably elevates the concentrations over the NH oceans, namely North Pacific and North Atlantic. The seasonal variations are strong with the burden in July significantly higher than in January. The biogenic nature of PI sulfate causes higher burdens in lower latitudes and during summertime. The PI and PD global annual mean sulfate burdens are 0.22 and 0.81 Tg, respectively. Both fossil fuel and biomass burning gives rise to PD OC. Though O’Dowd et al. (2004) showed evidences for oceanic sources of OC, they are not considered in the simulations in the absence of a global inventory. The three industrial source regions of sulfate have some of the highest OC burdens, the rest of which occur over the regions characteristic of heavy biomass burning, namely India, central Africa, and South America. The global annual mean OC burden is 0.17 Tg for PI and 1.36 Tg for PD. An organic matter (OM)/OC ratio of 1.67 is used to convert the burden of OC to that of OM. The PI and PD burdens of naturally occurring sea salt are assumed to be the same, and are highest over the oceanic regions with high wind speed. The global annual mean burden of fine sea salt is 0.02 Tg.
b. Distributions of Nd
Figure 3 depicts the geographical distribution of in-cloud Nd at 904 mb from a simulation with the prognostic cloud scheme driven by the PD aerosol climatology. Note that Nd values are linearly averaged in this study. The corresponding tendencies of convective and large-scale sources and sinks of n are plotted in Fig. 4. In January, relatively high Nd (over 100 cm−3) is seen over the tropical and subtropical oceans and land, especially the regions off the coast of East and South Asia, Africa, and South America, where droplets form mainly from convection (Fig. 4). In contrast, Nd is generally less than 30 cm−3, much lower than 500–100 cm−3 diagnosed from sulfate mass concentrations using BL95 (Ming et al. 2005). The large-scale processes over the NH midlatitude land are generally weaker, producing Nd less than 50 cm−3. These are smaller than 100–150 cm−3 diagnosed using BL95. Because the measurements presented in Leaitch et al. (1992) showed strong seasonal variations and were used as the only dataset in fitting the sulfate–Nd relationship over land by BL95, Nd in January diagnosed using BL95 may have positive biases. In July, the convectively active regions in the Tropics and subtropics have Nd frequently in excess of 100 cm−3 over ocean, which are often higher than those over land. This is contrary to what BL95 suggests and gives rise to a land–ocean contrast. A detailed discussion on the discrepancies is given in section 4b.
A closer look at Figs. 3 and 4 reveals that although strong tendency of sources usually leads to higher Nd, the opposite occurs occasionally. For example, despite the strong convection over the SH subtropical oceanic region off the west coast of Africa in July, the equilibrium Nd is less than 10 cm−3. The explanation lies in the second indirect effect. Because lower Nd has a tendency to accelerate cloud dissipation and shorten lifetime, droplets form more frequently than at higher Nd. So, even though Nd is low due to the lack of CCN, the enhanced frequency of cloud formation/dissipation gives rise to strengthening of source tendency.
The vertical distribution of zonal-mean Nd is shown in Fig. 5. As droplet evolution interacts with the model meteorology, Nd simulated with the prognostic scheme exhibits a vertical structure distinct from diagnosed Nd, for which the signature of vertical transport of ground-based emissions with maxima adjacent to the surface is clearly borne out. As an overall trend, prognosticated Nd increases with altitude starting from the surface, and then decrease after reaching maxima at around 800–900 mb. The vertical distributions of zonal-mean tendencies of sources and sinks are plotted in Fig. 6. Convective detrainments are much stronger than large-scale sources at the pressure levels where the maxima in Nd are located. As compared to large-scale updraft, convective updraft is much stronger, rendering Nd from convective detrainments higher. This suggests that the main feature of the vertical distributions of Nd (maxima in the lower troposphere) results from convective detrainments. Nd gradually decreases to less than 10 cm−3 for the water fraction of mixed-phase clouds above the boundary layer. In contrast, the diagnosed Nd is not constrained by the boundary layer, and could be as high as 100 cm−3 for mixed-phase clouds. Since the field measurements on which the empirical relationship of BL95 is based were taken for low-altitude warm clouds, its applicability to mixed-phase clouds is uncertain, especially given the differing meteorological conditions of the two cloud types.
c. Comparison with satellite measurements
We use the GCM with prognostic treatment of Nd to simulate the PD liquid water path (LWP) over ocean. The results in January are characteristic of two major peaks at 30°N and 45°S and two others close to the equator, a pattern similar to that in the Special Sensor Microwave Imager (SSM/I) measurements (Greenwald et al. 1995; Weng and Grody 1994; Fig. 7). The pattern of simulated LWP is generally in good agreement with the observations in both hemispheres, but underestimates in the high latitudes. In January, the model’s peaks in LWP are displaced relative to those observed. The model accurately catches the increase in LWP in the NH midlatitudes in July, while underestimating LWP in the Tropics and in the SH midlatitudes. In contrast, the use of diagnosed Nd yields much higher LWP compared to prognosticated Nd, and worsens the agreement with measurements. Most notably, the model with diagnosed Nd overestimates LWP almost by a factor of 3 in the NH midlatitudes in July. On the other hand, as discussed later in section 4b, the mean values of prognosticated Nd are within a factor 2 of diagnosed Nd. Autoconversion is particularly efficient in removing cloud water under conditions characteristic of low Nd as the influence of Nd on the autoconversion rate is highly nonlinear. Though the average values of Nd from the prognostic scheme are comparable to BL95, the former leads to more autoconversion than the latter because of the large variations of prognosticated Nd. This is the main reason why the prognostic scheme yields lower LWP than the diagnostic scheme.
The GCM-simulated PD cloud optical depth is compared with the International Satellite Cloud Climatology Project (ISCCP) measurements in Fig. 8. Note that logarithmic (as opposed to linear) averaging is applied in computing the monthly mean cloud optical depth since the change in irradiance measured by satellites and used in deriving cloud optical depth is less sensitive to optically thick clouds than thin clouds as a result of photon saturation effect. The measured cloud optical depth is relatively uniform in the low latitudes both in January and in July. The model predictions are typically twice as high as the measurements. The general trend of the ISCCP data in the midlatitudes is that cloud optical depth increases with latitude from the equator, opposite to what the model suggests except in the NH summertime. Overall the agreement is within a factor of 2. The overestimation caused by using diagnosed Nd in the NH midlatitude in July is even higher, by almost a factor of 4.
Except in the high latitudes, the GCM-predicted total cloud amount largely repeats the latitudinal pattern of the ISCCP measurements (Fig. 9). The model underestimates total cloud amount by less than 10% in the Southern Hemisphere (SH) midlatitudes in January, and in the NH midlatitudes and in the SH low latitudes in July. The different schemes of Nd cause little change in the simulated total cloud amount.
The PD shortwave (SW) and longwave (LW) cloud forcings diagnosed from the GCM are shown in Figs. 10 and 11, respectively. The prognostic cloud scheme achieves a satisfactory agreement between the model predictions and Earth Radiation Budget Experiment (ERBE) measurements in terms of shortwave cloud forcing except underestimations in the SH midlatitudes in January and in the NH midlatitudes in July. In contrast, an overestimation of 16 W m−2 in SW cloud forcing occurs at 60°N in July when diagnosed Nd is used. This is not nearly as significant as in the case of LWP in terms of percentage increase, largely due to the nonlinear relationship between radiation and LWP. The model’s underestimation of total cloud fraction also partly compensates for the overestimated LWP. A global annual mean net radiation of 1.4 W m−2 arises from an absorbed SW radiation of 236.7 W m−2 and an outgoing LW radiation of 235.3 W m−2 when diagnosed Nd is used. The global annual mean absorbed SW and outgoing LW radiations are 240.0 and 236.2 W m−2, respectively, resulting in a net radiation of 3.8 W m−2 when prognosticated Nd is used. It is noted that the GCM has not been retuned in any manner in the wake of the physics implementation here.
d. Flux changes
The difference in TOA radiation flux between two simulations with the prognostic cloud scheme driven separately by PI and PD aerosol climatology is employed to assess the sign and scale of indirect effects (designated as the BASE case). The aerosol climatology used for clear-sky calculation remains the same at the PD level in both runs, and thus neither direct nor semidirect effect is included. Note that both runs allow the model meteorology to evolve, so the calculated flux change is a measure of combined first and second effects, and includes the impact of model feedbacks (Ming et al. 2005). As seen in Fig. 12, owing to the response of the model to the change in Nd, regional radiative warming appears in addition to the radiative cooling; this is scattered across the globe, especially over the SH oceans and Australia mainly as a result of the model’s natural variability (radiative warming or cooling implies, respectively, deficit or surplus of column radiation when referring to flux change). Nonetheless, radiative cooling typically ranging from −3 to −10 W m−2 is found over the NH midlatitude industrial regions because of high burdens of sulfate and organic aerosols. Strong radiative cooling is widespread over the NH oceans as a result of long-range transport of continental aerosols. Mainly because of high burdens of biomass burning OC, the tropical parts of South America and Africa undergo radiative cooling ranging from −5 to −20 W m−2, which is even stronger than that over the sulfate-dominant industrial regions. The statistical significance of the geographical distribution of flux changes will be discussed later. The global annual mean flux change amounts to −1.8 W m−2, including an SW component of −2.0 W m−2 and an LW component of 0.2 W m−2 (Table 1). The ratios of the flux change in NH to that in SH (NH/SH) and of the flux change over ocean to that over land (ocean/land) are 2.9 and 0.73, respectively. A recent review by Lohmann and Feichter (2004) stated that the average value of the total indirect effect estimated in some latest studies is approximately −1.7 W m−2 with one standard deviation of 0.6 W m−2, and the average ocean/land ratio is 0.78 with one standard deviation of 0.53. The estimates here are very close to the average values.
The geographical distributions of the flux changes in SW and LW are plotted in Fig. 13. The SW radiative cooling is stronger and more widespread than when the net radiation is considered. The LW components of flux changes are generally smaller than SW, and are largely due to model noise.
Higher prognosticated Nd incurs a percentage increase in annual mean LWP of 19% (from 47 g m−2 for PI to 56 g m−2 for PD), while the increase in total cloud amount is 0.6% (from 64.6% for PI to 65.0% for PD). The different responses of cloud water (l) and amount (a) can be attributed to the ways in which they are modeled. The time derivative of a is proportional to (1 − a)2, which means that the increase in a is limited as a approaches 1. In contrast, the time derivative of l is proportional to a, which means that l can keep on increasing even if the grid box is completely cloudy. The cloud amount in the preindustrial case has already been frequently limited, and could not be significantly higher than in the present-day case. Another note is that though an increase of 0.6% in cloud fraction appears to be small, it is not so for radiation. A first-order calculation shows that if clouds cover 60% of the earth and have an average cloud forcing of 30 W m−2, the increased cloud fraction causes an increase of reflection of 0.3 W m−2, which is small compared to the total SW indirect effects, but by no means negligible.
As shown in Fig. 14, the largest increases in LWP occur over the NH midlatitude industrial regions and oceans, as well as over the tropical biomass burning regions as higher Nd tends to prolong cloud lifetime. The geographical pattern is very similar to that of flux changes (Fig. 12), indicating that enhanced LWP is the underlying mechanism for radiative cooling caused by aerosol indirect effects. When diagnosed Nd is used, the global annual mean LWP changes from 70 g m−2 for PI to 78 g m−2 for PD, representing an increase of 11%. In this case, the increased LWP centers predominantly in the sulfate-dominant NH midlatitudes since BL95 uses sulfate as a surrogate for all aerosol species.
e. Role of OC in aerosol–cloud interactions
Numerous experimental and modeling studies (e.g., Cruz and Pandis 1997; Corrigan and Novakov 1999; Russell et al. 2000) have indicated that soluble and/or surface-active OC can be efficient CCN. The prognostic cloud scheme employs a parameterization of CCN activation that has been shown to be able to satisfactorily represent OC (Ming et al. 2006). In a sensitivity case designed to evaluate the role of OC in aerosol–cloud interactions on the global scale (the no-OC case, or NOOC), the PI and PD prescribed SST simulations are rerun with only sulfate and sea-salt aerosols. The size distributions remain unchanged as OC is assumed to be externally mixed with other aerosols. The calculated flux changes are plotted in Fig. 15. The omission of OC considerably weakens the radiative cooling over the tropical parts of central Africa and South America, where high burdens of biomass burning OC are located (Fig. 2). In the NH midlatitudes, the radiative cooling is also reduced but to a lesser extent because indirect effects are largely due to sulfate aerosol. The resulting global annual mean flux change is −1.1 W m−2 as compared to −1.8 W m−2 in BASE and −2.3 W m−2 calculated from Nd diagnosed with BL95 (Ming et al. 2005; Table 1). As shown in Fig. 15, the geographical pattern of the flux changes calculated using diagnosed Nd is generally similar to the NOOC case. Both cases yield comparable radiative cooling over the NH midlatitude land, but the values over the oceans in NOOC are not nearly as strong as suggested by using diagnosed Nd. This is because the relative increases in Nd with respect to the sulfate mass concentration are comparable over the NH midlatitude land for both cases, while prognosticated Nd shows no pronounced dependence on sulfate over the oceans (see section 4b). Note that the conclusions drawn here depend on the assumed properties of organic aerosol. If OC is less soluble than assumed or internally mixed with other species, its impacts would not be as important as suggested by this sensitivity study.
f. Statistical significance
The model’s natural variability is responsible for inherent fluctuation of radiation flux, which may affect the statistical significance of the flux change resulting from perturbation in aerosol climatology. The model’s zonal natural variations are represented by the standard deviations estimated from a total of 10-yr simulation with aerosol climatology unchanged at PI level. For this particular version of the GFDL GCM, relatively large zonal-mean variations can be found near to the equator and in the high latitudes, while the variations are generally less than 1 W m−2 in the low and midlatitudes (Fig. 16). The calculated zonal-mean flux changes in BASE exceed the standard deviations to various extents in the NH low and midlatitudes and in the SH low latitudes, implying that the probability of these flux changes being forced perturbations as opposed to natural variations is at least greater than 84%. In the NH midlatitudes, where the strongest radiative cooling is located, the flux changes are approximately 3 times greater than the corresponding standard deviations, effectively ruling out the possibility that they are caused by the model’s natural variability. Similar to BASE, the NOOC case yields unambiguous flux changes in the NH midlatitudes, but considerably weakens flux changes in the SH low latitudes to levels comparable to or less than standard deviations. This confirms the significant role played by OC in aerosol–cloud interactions over the biomass burning regions. A similar statistical analysis is performed for the geographical distribution of flux changes (Fig. 12). The strong radiative cooling over the NH midlatitude industrial regions (particularly East Asia and Europe) and oceans and over the tropical parts of South America and Africa all exceeds at least one standard deviation, and is considered statistically significant. In contrast, the scattered radiative warming mostly in the SH is usually within one standard deviation, and largely results from the model’s natural fluctuation.
a. Subgrid variability of updraft velocity
Owing to convection and large-scale processes, the upward motion of moist air is the fundamental driving force behind cloud formation. Critical for prognostic determination of cloud water as well as droplet number in GCMs, updraft velocity (w) can be thought of as the summation of grid-mean (w) and subgrid variation (w′). In reality, it is w that drives cloud formation, but only grid-mean w can be resolved from GCMs. In the earliest study of prognostic scheme of Nd, Ghan et al. (1997) argued that since w′ is necessary for accounting for the impact of subgrid variability of w on CCN activation, its probability distribution can be diagnosed either from the turbulent kinetic energy (TKE) or from the vertical eddy diffusivity (K), and then relating the standard deviation of w (σw) to TKE. A simpler method used in Lohmann et al. (1999) assumed w = w + cTKE with c empirically set at 0.7. Ghan et al. (1997) did not distinguish convective and large-scale sources of cloud droplets, and Lohmann et al. (1999) effectively assumed negligible aerosol impact on convective source by considering a constant radius of 10 μm for these droplets.
If approaches that use grid-mean w to drive cloud formation are applied to the GFDL GCM, the inconsistency in updraft velocity poses a fundamental problem in that the same updraft field should drive both cloud droplet growth and liquid water increase simultaneously, whether subgrid variability is explicitly considered or not in a GCM. In other words, if subgrid components of w are taken into account in determining Nd, the impact on l ought to be fully accounted for. If one uses w for driving cloud liquid and w′ for droplets, the combination of a zero w and a nonzero w′ (as is the case for stratocumulus) would allow droplets to be formed in the absence of cloud liquid water, projecting an unphysical picture.
This study approaches the issue of subgrid variability by first drawing upon observations of the real atmosphere and recognizing that two major types of updraft (i.e., convective and large-scale) operate. Rather rapid ascents take place in convective systems, and w can reach a few meters per second, while relatively slow large-scale updrafts are on the order of centimeter per second. Thus convection could be an important contributor to subgrid variability if present. Approaches developed for modeling convection (e.g., the RAS convection scheme used here) treat ensembles of convective clouds and enable subgrid variability of w resulting from convection to be fully accounted for. In the case of large-scale stratiform clouds, the choice of grid-mean w for calculating CCN activation allows us to apply the same w to forming cloud water and droplets. The GCM simulations described here show that this treatment can produce enough droplets to sustain LWP at levels comparable to measurements. Our approach is able to account for the impact on Nd of subgrid variability of w to some extent without having to diagnosing subgrid variability empirically. Therefore, it represents a reasonable alternative to the previous formulations of Ghan et al. (1997) and Lohmann et al. (1999), who must impose an artificial lower bound on w; this is avoided in the present study. If a large-scale cloud scheme of cloud water and droplets formulated explicitly on subgrid-scale w can be devised, that would be the best tool for studying indirect effects.
Because the RAS scheme is unable to resolve the vertical structure of updraft velocity, this study has to estimate it as 0.5CWF for an ensemble of convective clouds. This may lead to an underestimation of updraft velocity and Nd. Nonetheless, such a method differs fundamentally from deriving subgrid variations from grid-mean TKE and applying them to driving large-scale stratiform clouds. The limitation posed by diagnosing convective updraft from CWF can be readily overcome since an explicit vertical distribution of updraft is well within the reach of an improved convection scheme (Donner 1993).
b. Comparison with other studies
Numerous studies, including the most recent one by Ming et al. (2005), used the empirical relationship of BL95 relating Nd to mass concentrations of sulfate (a surrogate for all aerosols) to simulate indirect effects. Ming et al. (2005) reported a global annual mean flux change of −2.3 W m−2, stronger than −1.1 W m−2 in the NOOC case of this study (we believe that NOOC is more appropriate for comparison than BASE since the BL95 relationship does not account for OC explicitly). In light of the sharp contrast between the simplistic nature of BL95 and the complexity surrounding the prognostic scheme, the large differences in flux change (1.2 W m−2 in global annual mean), despite similar geographical patterns (Fig. 15), are not unexpected, but still prompt a further investigation of the practicality of BL95, with the assistance of the prognostic cloud scheme used here.
Scatter plots of predicted Nd and corresponding sulfate concentrations at the grid points located from 20° to 60°N in BASE are given in Fig. 17. These latitudes cover all major industrial regions (i.e., East Asia, Europe, and the United States), and are characteristic of relatively high sulfate burdens (Fig. 2). Because BL95 was fitted from field data gathered in North America and North Atlantic, it presumably reflects sulfate–Nd correlation in these sulfate-dominant regions, if there is such a correlation. It can be seen from Fig. 17 that a large scatter is associated with prognosticated Nd for sulfate concentrations, spanning almost two orders of magnitude from 10 to 1000 cm−3. This indicates that factors besides sulfate concentrations (e.g., updraft velocity and chemical composition) contribute to the determination of Nd. Nonetheless, Nd over the NH midlatitude land generally increases with the sulfate mass concentration in the range of 0.8–4 μg m−3 that covers a large fraction of the data points, reaffirming the role of sulfate in providing CCN for cloud formation. On the other hand, the slope of mean Nd, which dictates the relative increases in Nd due to anthropogenic sulfate aerosol, and subsequently the flux changes, are rather similar to BL95 in that range. This may be partly responsible for the similar indirect effects over the NH midlatitude land calculated with prognosticated Nd as well as with diagnosed Nd. In contrast, mean Nd over ocean does not show strong dependence on the sulfate mass concentration except, implying that marine clouds may not be as susceptible to aerosols as suggested by BL95. This explains why the radiative cooling over the NH midlatitude oceans in NOOC is much weaker than when diagnosed Nd is used. Figure 18 is for OC-abundant 20°S to 20°N regions. Here, Nd is generally higher than in the midlatitudes as a result of stronger convective source (Fig. 4). Mean Nd even decreases with respect to the sulfate mass concentration when it is above 1 μg m−3, as the presence of large amount of OC that provides a significant fraction of CCN over this region renders the dependence of Nd on sulfate less significant. So, BL95 may not be as applicable to the tropical and subtropical regions as to the midlatitudes.
An observation of Fig. 17 is that the absolute values of prognosticated Nd are generally lower than BL95 over the NH midlatitude land. Though this may hint at potential negative biases of the scheme, it is worth noting that the measurements, from which the BL95 relationship over land was fitted, were mostly taken over central Ontario (Leaitch et al. 1992). It is for the most part unclear how well BL95 is representative of, and hence applicable to other continental regions especially given the large variations among such empirical relationships based on different field measurements (e.g., Twohy et al. 2005). This prompts us to make an effort to compare prognosticated Nd with BL95 over central Ontario at the pressure level of 844 mb, which is chosen in light that most samples in Leaitch et al. (1992) were collected at altitudes from 1 to 2 km. As shown in Fig. 19, the model predicts no liquid clouds over central Ontario in January, rendering it impossible to make a direct comparison with BL95. The prognosticated Nd in July are in the range of 150–200 cm−3. Though somewhat lower than 200–300 cm−3 diagnosed from BL95, the prognosticated Nd are largely within the natural scatters inherent in the Nd measured by Leaitch et al. (1992) ranging approximately from 100 to 500 cm−3, which partly result from the variance in updraft velocity (Penner et al. 2001). This comparison proves inconclusive for possible signs of underestimation.
A comprehensive database of measured Nd of stratiform clouds can be found in Miles et al. (2000). However, the limited number of samples does not allow us to compare them with the GCM predictions at a particular site in a statistically significant fashion in light of the large variations in measured Nd even for the same location and season. Nonetheless, prognosticated Nd does lie at the low end of measured Nd in a global mean sense. For continental clouds, Miles et al. (2000) reported a global mean Nd of 288 cm−3 with a standard deviation of 159 cm−3, while prognosticated Nd is frequently below 100 cm−3. The scheme uses grid-mean large-scale updraft velocities for driving droplet activation in stratiform clouds, and thus may be subject to negative biases by neglecting subgrid variations of large-scale ascents. This shortcoming particularly precludes adequately accounting for the presence of stratocumulus even though the convection scheme in the model partly compensates for the problem by generating shallow convective clouds over regions where stratocumulus could be important. In addition, the convection scheme appears not to generate enough convective clouds over land (Fig. 4). Both factors may contribute to the scheme’s possible negative biases as well as the strong land–ocean contrast in Nd (Fig. 3). Note that the assumed aerosol properties, namely mixing state, size distribution, and composition may also bias the model, but the signs and magnitudes of potential biases cannot be reliably quantified due to a lack of measurements on the global scale. A previous study by Ming et al. (2006) on CCN activation suggests that droplet numbers are probably more susceptible to updraft velocity than to reasonably assumed aerosol properties. Hence the assumed aerosol properties may not be as important sources of model biases as the factors contributing to underestimation of updraft velocity.
The agreement between prognosticated Nd and BL95 is much better over ocean than over land; mean Nd over ocean are mostly within a factor of 2 of BL95. They also generally fall into the range of measured Nd at several oceanic sites (Ramanathan et al. 2001) at low sulfate concentrations (<1 μg m−3), but are close to the upper bound from 20°S to 20°N when the concentration is below 0.3 μg m−3. As another contributing factor to the land–ocean contrast in Nd, the relatively higher Nd over the tropical and subtropical oceans result from shallow convection. So, any tendency of the RAS scheme to overestimate shallow convection will inevitably give rise to too high Nd. Research is currently under way to better understand the roles of different convection schemes on modeling aerosol–cloud interactions and indirect effects.
The prognostic schemes of Nd developed by Ghan et al. (1997) and Lohmann et al. (1999) have been utilized to estimate indirect effects of sulfate aerosol (Ghan et al. 2001; Lohmann et al. 2000). The differences in the treatment of subgrid variability, parameterizations of activation and autoconversion, and other aspects of the schemes make it impossible to make a clean comparison of this study with the two previous ones. Instead, we reiterate some findings of Ghan et al. (2001) that are useful for putting this study into perspective. First, the large difference in the calculated indirect effects due to externally mixed anthropogenic sulfate aerosol between Lohmann et al. (2000) (−0.1 W m−2) and Ghan et al. (2001) (−2.4 W m−2) is attributable to the arbitrary lower limits on aerosol and droplet numbers. Note that this study avoids using such limits and differs from Ghan et al. (1997) (a lower limit on σw of 0.1 m s−1) and Lohmann et al. (1999) (a background Nd value of 40 cm−3) in this regard. This adds to the usefulness and novelty of the prognostic scheme of this study, as it eliminates the need for such artifacts. Second, Ghan et al. (2001) examined the sensitivity of indirect effects to aerosol numbers diagnosed with a fixed normalized distribution, as is the case for this study, and found that the flux change increases from −2.4 W m−2 with predicted aerosol numbers to −3.1 W m−2 with diagnosed aerosol numbers. The NOOC case of this study obtains a global annual mean SW flux change of −1.5 W m−2, much weaker than what Ghan et al. (2001) found when aerosol numbers are diagnosed. The underlying reasons for the differences can be fully explained only by a carefully designed model intercomparison and is beyond the scope of this study.
This paper describes an effort to expand a prognostic cloud scheme of cloud liquid water and amount to include droplet number concentration. It does so in a way that allows them to be calculated using the same large-scale and convective updraft velocity field. The evolution of droplets fully interacts with the model meteorology. Explicitly treated as a major source of droplets, CCN activation is parameterized using the microphysically based approach proposed by Ming et al. (2006), which allows the scheme to take into account the impacts on Nd of multiple aerosol species (i.e., sulfate, OC, and sea salt), decreased droplet surface tension due to OC, and kinetic limitations. An implementation of the scheme in the GFDL GCM yields a vertical structure of Nd characteristic of maxima at around 800–900 mb, as opposed to being adjacent to the surface for Nd diagnosed from sulfate mass concentrations using BL95. This can be attributed to detrainment of convection at higher altitudes. The agreement between model diagnostics and satellite-measured LWP, cloud amount and forcing is improved compared to using diagnosed Nd.
The simulations with PI and PD aerosols show that the combined first and second indirect effects of anthropogenic sulfate and organic aerosols give rise to a global annual mean flux change of −1.8 W m−2, consisting of −2.0 W m−2 in SW and 0.2 W m−2 in LW, with the model response altering the cloud fields, and subsequently LW radiation. Specifically, LWP and total cloud amount increase by 19% and 0.6%, respectively, from PI to PD conditions. The radiative cooling is centered over the NH midlatitude industrial regions (i.e., East Asia, Europe, and United States) and oceans (i.e., North Pacific and North Atlantic) mainly due to sulfate aerosol, and over the tropical parts of South America and Central Africa mainly due to organic aerosol. The NH/SH and ocean/land ratios are 2.9 and 0.73, respectively. The calculated annual zonal-mean flux changes are determined to be statistically significant, exceeding the model’s natural variations in the NH low and midlatitudes and in the SH low latitudes. Sensitivity analyses show that anthropogenic sulfate alone gives rise to an annual mean flux change of −1.1 W m−2, and the absence of organic aerosol significantly weakens the radiative cooling over the tropical land.
The contribution of S. Klein was under the auspices of the U.S. Department of Energy at the University of California Lawrence Livermore National Laboratory under Contract W-7405-Eng-48.
Model Parameters and Numerical Integration
The threshold value of grid-mean relative humidity above which stratiform clouds start to form is set at 80%. This study assumes that convective precipitation removes water condensate and droplets proportionally. As a function of detrainment pressure, the removal fraction rj in Eq. (5) is 0.5 below 800 mb, 0.975 above 500 mb, and linearly interpolated at intermediate pressure levels. The cloud erosion rate constant is set at 4.7 × l0−6 s−1 in convective conditions, 5.0 × l0−5 s−1 in turbulent conditions, and 10 × 10−6 s−1 otherwise. The mass and heat accommodation coefficients used in the parameterization of CCN activation are both 1 (Laaksonen et al. 2004).
Because the time step of the GCM cloud scheme (30 min) is much longer than the characteristic time scale of Eq. (6), analytical integration is necessary. If large-scale activation is treated as constant and the other microphysical processes as linear to n, Eq. (6) turns into
Corresponding author address: Yi Ming, Geophysical Fluid Dynamics Laboratory, Princeton, NJ 08542. Email: Yi.Ming@noaa.gov