MATRIX (Multiconfiguration Aerosol TRacker of mIXing state): an aerosol microphysical module for global atmospheric models

A new aerosol microphysical module MATRIX, the Multiconfiguration Aerosol TRacker of mIXing state, and its application in the Goddard Institute for Space Studies (GISS) climate model (ModelE) are described. This module, which is based on the quadrature method of moments (QMOM), represents nucleation, condensation, coagulation, internal and external mixing, and cloud-drop activation and provides aerosol particle mass and number concentration and particle size information for up to 16 mixed-mode aerosol populations. Internal and external mixing among aerosol components sulfate, nitrate, ammonium, carbonaceous aerosols, dust and sea-salt particles are represented. The solubility of each aerosol population, which is explicitly calculated based on its soluble and insoluble components, enables calculation of the dependence of cloud drop activation on the microphysical characterization of multiple soluble aerosol populations. A detailed model description and results of box-model simulations of various aerosol population configurations are presented. The box model experiments demonstrate the dependence of cloud activating aerosol number concentration on the aerosol population configuration; comparisons to sectional models are quite favorable. MATRIX is incorporated into the GISS climate model and simulations are carried out primarily to assess its performance/efficiency for global-scale atmospheric model application. Simulation results were compared with aircraft and station measurements of aerosol mass and number concentration and particle size to assess the ability of the new method to yield data suitable for such comparison. The model accurately captures the observed size distributions in the Aitken and accumulation modes up to particle diameter 1 μm, in which sulfate, nitrate, black and organic carbon are predominantly located; however the model underestimates coarse-mode number concentration and size, especially in the marine environment. This is more likely due to oversimplifications of the representation of sea salt emissions – sea salt emissions are only calculated for two size classes – than to inherent limitations of MATRIX.


Introduction
The impact of natural and anthropogenic aerosols on the climate system is the subject of numerous laboratory, experimental and theoretical studies (Ghan and Schwartz, 2007). These studies span spatial and temporal scales of several orders of magni-tool to estimate the integrative impacts of aerosols on climate and their role in climate change. The treatment of aerosol properties and processes in regional and global models is becoming increasingly complex (Ackermann et al., 1998;Adams et al., 2001;Binkowski and Shankar, 1995;Binkowski and Roselle, 2003;Easter et al., 2004;Gong et al., 2003;Herzog et al., 2004;Jacobson, 2001;Lauer et al., 2005;Riemer et al., 5 2003; Stier et al., 2005;Wilson et al., 2001;Wright et al., 2000Wright et al., , 2001Yu et al., 2003), taking into account more detailed microphysical and chemical interactions that underlie the formation of aerosol particles which then impact the climate system. More detailed representations of key aerosol properties such as size distributions, chemical composition (e.g., internal vs. external mixtures), hygroscopicity, optical properties, and cloud 10 activation are also being incorporated into global climate models. With this increased level of detail, the question arises as to which method is most appropriate.
The most common methods of including aerosol size distributions and aerosol mixing information into models are the sectional, modal, and moment based schemes. Sectional models (Adams et al., 2001;Gong et al., 2003;Jacobson, 2001), which divide 15 the size domain into intervals, or bins, calculate the evolution of the number concentrations in each size bin. This method has the disadvantage of being computationally expensive, because of the large number of bins required to accurately represent size and chemical mixing information. The modal method (Wilson et al., 2001;Iversen, 2002;Stier et al., 2005) approximates the size distribution by an assumed statistical 20 function, most commonly multiple lognormal functions. The quadrature method of moments (QMOM) (McGraw, 1997) provides a computationally efficient statistically-based alternative to modal and sectional methods for aerosol simulation, that does not make a priori assumptions about the shape of the size distribution. Numerous modules based on the flexible and efficient moments method tracking 2-6 moments have been imple-Introduction Interactive Discussion mass) that enter the covariance matrix of a principal component analysis are tracked instead of the distribution itself. The approach is flexible and highly efficient, yet can provide comprehensive representation of natural and anthropogenic aerosols and their mixing states. Until recently aerosols in the GISS climate model  were treated 5 in a mass-based scheme, in which the life cycles of sulfate, nitrate, black and organic carbon, sea salt, and dust were relatively independent from each other, except for surface reactions on mineral dust (Bauer and Koch, 2005;Bauer et al., 2007;Koch et al., 2006;Miller et al., 2006). Size information was included only for dust and sea salt (using a bin scheme). For other species, size distributions were specified, number 10 concentrations were not calculated, and water uptake effects were parameterized to depend upon relative humidity in the radiation and gravitational settling schemes. The implementation of MATRIX into the climate model will provide detailed aerosol characterization for formation, removal, and climate interactions to be calculated. In the current module each of several aerosol modes is treated by a two-moment method, 15 in which a set of mass and number concentrations is evolved. Efficient solvers based on analytic solutions to the moment equations are available (Binkowski and Roselle, 2003). This two-moment method is free of the concerns with moment consistency associated with larger moment sets, as number and mass concentrations are mathematically independent (unlike sets of three or more moments per mode), and solver 20 errors cannot generate inconsistent moment sets (Wright, 2007). The goal of this paper is to give a complete model description of the microphysical model MATRIX (Sect. 2), present applications of MATRIX as a box model (Sect. 3), and as part of the GISS global climate model (Sect. 4). The MATRIX model provides a wide set of possible model configurations of mode formulations and microphysical parame-but less than that of a sectional representation of multiple aerosol populations.
One goal of the model is to provide the number concentration and particle size information characterizing soluble particles for the treatment of cloud drop activation. This requires a representation of internal and external aerosol mixing processes. For example, insoluble particles can become soluble by acquiring soluble inorganic coatings 10 as a result of either condensation of soluble species or coagulation with populations containing soluble materials.
The total aerosol is represented by a user-selected set of distinct populations (or modes), along with a specification of the coagulation interaction between the modes, that we refer to as an "aerosol mechanism". Each mode has a distinct composition or 15 set of chemical components, independent of size. Modes may be primary (those that receive particle emissions), secondary (those formed by coagulation among primary modes or condensation of gaseous components onto primary modes); or the mixed mode including all aerosol constituents. When particles from two modes coagulate such that the resulting particle cannot be accommodated in one of the defined modes, Introduction Interactive Discussion dust coarse mode) represent soluble particles produced by condensation of inorganics onto DD1 (insoluble dust accumulation mode) and DD2 (insoluble dust coarse mode), respectively. Modes DD1 and DD2 can also acquire inorganics through coagulation with mode AKK or ACC. Particles are transferred from DD1 to DS1 or DD2 to DS2 when the volume fraction of inorganics in DD1 or DD2 exceeds a threshold (set to 5%).

5
The secondary mode DBC (dust-BC) receives particles solely through coagulation of any of the modes DD1, DD2, DS1, DS2 with any of the modes BC1, BC2, BC3, BCS. The mixed mode MXX receives the more complex aerosol mixtures and grows as the total aerosol approaches an internally mixed state. The transported species for the sea salt modes are the mass concentrations of dry sea salt (treated as NaCl) and the 10 total non-sea salt sulfate summed over both modes. The non-sea salt sulfate is apportioned between the SSA and SSC modes in proportion to the sea salt (NaCl) in each mode. Number concentrations are derived from mass concentrations using assumed lognormal distributions (D g,N , σ g,N ) given in Table 1a.
Mechanism 2, which also consists of 16 modes and 51 transported species, differs 15 from mechanism 1 only in omission of mode BC3 and inclusion of mode OCS (OCsulfate) resulting from coagulation of mode OCC with mode AKK or ACC, analogous to mode BCS. Mechanism 3 consists of 13 modes and requires 41 transported species. It differs from mechanism 1 in that (1) modes BC3 and BCS are omitted and black carbon aerosols are represented only by modes BC1 (insoluble) and BC2 (soluble), 20 and (2) mode DBC (mineral dust-BC) is omitted and these particles are now placed in the mixed mode MXX. Mechanism 4 consists of 10 modes and requires 34 transported species, described in Table 2. It differs from mechanism 1 in the omissions of modes AKK, BC3, BCS, DBC, and BOC, and the use of a single mode for sea salt (SSS). Mechanisms 5-8 25 differ from mechanisms 1-4 only in the representation of mineral dust by one insoluble and one soluble mode only. The numbers of modes are 14, 14, 11, and 8, and the number of transported species are 45, 45, 35, and 28, respectively. 9938 2.2 Equations for evolution of number and mass concentrations Equations for number and mass concentrations for each mode must be solved. These equations are similar to those solved by the Community Multiscale Air Quality (CMAQ) chemical transport model as described in Binkowski and Roselle (2003). Secondary particle formation from gaseous precursors, condensational growth, all coagulation 5 processes, and addition of sulfate formed by in-cloud oxidation are treated as a group without operator splitting. The partitioning of semi-volatiles between the gas-and particle-phases is done just prior to, and again subsequent to, the calculation of the primary dynamical processes just cited. Any reclassification of particles from one mode to another is done at the end of the time step. For the number concentration N i of where P i and L i are the production and loss rates, J i the secondary particle formation 15 rate (to be distinguished from the nucleation rate below), E num i the number concentration emission rate, R i the rate at which particles enter mode i due to coagulation among all other modes. K (0) kl is the mode-average coagulation coefficient (defined in Sect. 2.6 and Appendix B) for modes i and j , there are n modes, d i kl is unity if coagulation of mode k with mode l produces particles in mode i and zero otherwise; d i j is 20 unity if coagulation of mode j with mode i results in the removal of particles from mode i and zero otherwise. The coefficients a i , b i , and c i are defined for convenience and 9939 Introduction Analytic solutions are obtained by holding the coefficients constant over a time step ∆t=t−t 0 and are expressed using the quantities For c =0 the approximated solution of Eq. (4) is For c=0 and b =0 the solution is 10 and for c=0 and b=0 the solution is where sign errors in Eqs. (7) and (8) have been corrected from Binkowski and Roselle (2003). For the mass concentration Q i ,q of species q in mode i is defined through the second equality. Analytic solutions are again obtained by holding the coefficient constant over the time step. For f (3) i =0 the solution is and for f (3) i =0 the solution by an Euler forward step is The mass concentration production and loss terms treated in the aerosol module are where terms are included for particle emissions, secondary particle formation, coag-10 ulation, growth due to condensation and other gas-particle mass transfer, and incorporation of sulfate produced by in-cloud oxidation. Each of these terms and those for number concentrations in Eqs. (2) and (3) is discussed in a following section.

Secondary particle formation
The secondary particle formation production term is where m npf i ,q is the mass of species q in a newly formed particle added to mode i , and J p,i is the new particle formation (NPF) rate, determined from and less than the nucleation rate J. is directly acquired by mode AKK; water and ammonium uptake are determined by subsequent equilibrium calculations. The nucleation rate is the rate at which criticalsized clusters of ∼1 nm diameter are formed. The NPF rate J p,i is usually defined as the rate at which particles at the minimum detectable size are formed and is thus a function of measurement techniques, but in the present context J p,i is the rate at which 5 particles of a user-selected small ambient diameter D npf (3-20 nm) enter the Aitken mode. Explicit representation of particle dynamics at sizes less than D npf is avoided through parameterizations that convert the nucleation rate J to the NPF rate at size D npf , or directly give a NPF rate derived from field observations without reference to a nucleation rate. This section describes the MATRIX NPF model, which contains 10 several nucleation parameterizations. The following submodules of the NPF model are described in Appendix A: (A1) parameterizations to convert J to J p,i by implicitly treating the particle dynamics at sizes less than D npf , (A2) calculation of the mass of sulfate in a newly formed particle, and (A3) calculation of the steady-state concentration [H 2 SO 4 ] SS used in the nucleation and NPF parameterizations. 15 Five nucleation or direct NPF parameterizations are available, two for binary homogeneous nucleation of H 2 SO 4 -H 2 O (Jaecker-Voirol and Mirabel, 1989;Vehkamaki et al., 2002), one for ternary homogeneous nucleation of H 2 SO 4 -NH 3 -H 2 O (Napari et al., 2002), one for nucleation by ion-ion recombination (Turco et al., 1998), and one for formation of particles of 3-nm diameter derived from field observations (Eisele and 20 McMurry, 1997). The ion-ion recombination mechanism is taken as independent of the homogeneous mechanisms, and the ion-ion recombination contribution to the total nucleation rate is optionally included or set equal to zero.
-The Jaecker- Voirol and Mirabel (1989) parameterization is based on classical binary nucleation theory with hydrates, presented graphically, and designed for sim- 25 ple interpolation in T and RH. Each plot of log J vs. log [H 2 SO 4 ] for T= [223,248,273,298,323] K and RH=[20,40,60,80, 100]% were scanned, digitized, and fit to low-order polynomials for multi-linear interpolation and extrapolation to 0% RH. The range of rates covered by these fits is 10 −3 -10 5 cm −3 s −1 .

9942
- The Vehkamaki et al. (2002) parameterization is also based on classical binary nucleation theory with hydrates. Valid for T=230. RH=0. The conversion of the nucleation rate into a NPF and the calculation of the sulfate mass in a newly formed particle is presented in Appendix A.

Condensational growth and gas-particle mass transfer
The total production rate (averaged over the time step) of species q in mode i due to condensation of nonvolatile species and subsequent gas-particle mass transfer due to 5 equilibration of other species is given by where P growth i ,q is the total production rate of species q in mode i averaged over the time step, P kinetic q represents species treated kinetically (H 2 SO 4 ), P equil i ,q represents species for which equilibrium is assumed (H 2 O, NH 3 , HNO 3 ), h i is the fraction of the total 10 H 2 SO 4 condensation rate due to particles in mode i , and P kinetic i ,q is the total condensation rate of H 2 SO 4 summed over all modes, calculated as where k c is the total condensation sink. The condensation sink is the first-order rate constant for loss of [H 2 SO 4 ] onto particle surfaces (Kulmala et al., 2001) and is obtained 15 as k c = n i =1 k c,i with k c,i the contribution of mode i . For each mode the calculation of k c,i is based upon a single particle size, the diameter of average massD p,i , and an adjustment factor θ i = exp[−(l nσ g,i ) 2 ] to prevent excessive condensation due to the use of a monodisperse distribution of diameterD p,i (Okuyama et al., 1988). For mode i where D diff is the binary diffusion coefficient of H 2 SO 4 in air, D p is the particle diameter, K n=2λ/D p is the Knudsen number with λ the vapor (H 2 SO 4 ) mean free path in air, 9944 α the mass accommodation coefficient of H 2 SO 4 , n i the number size distribution of mode i , N i the particle number concentration in mode i , and β the transition regime correction to the condensational mass flux (Seinfeld and Pandis, 1998), given by A particle density is needed for each mode to obtain the diameter of average mass and 5 is calculated from ρ=m/V , where the volume concentration V is calculated as a sum of contributions for each chemical component based on the current mode composition. The condensation sink is implemented through lookup tables for k c,i /N i as a function of particle diameter calculated at representative heights for each host model vertical level and using the standard atmosphere to an appropriate characteristic temperature 10 and pressure to calculate D diff and λ. The condensation sink was calculated with the mass accommodation coefficient α for H 2 SO 4 set to 0.86 (Hanson, 2005), α=1 and λ=λ ai r .
2.4.1 Partitioning (P equil i ,q ) of semi-volatile species Two thermodynamic modules are included in MATRIX. EQSAM (Metzger et al., 2002b(Metzger et al., ,a, 2006 which is built upon a simplified non-iterative expression for the activity coefficients, is very efficient. ISORROPIA (Nenes et al., 1998) with ξ(RH)=r/r dry . This holds for the fractional relative humidity RHF above the crystallization value 0.45. The expression for ξ is ξ (RHF)=c SS [b SS +1/(1 − RHF)] 1/3 with ρ NaCl the density of NaCl (2.165 g cm −3 ) and with c SS =1.08 and b SS =1.2.

20
Hysteresis effects are represented by the treatment of Ghan et al. (2001). For RH between the crystallization and deliquescence RH, the total metastable aerosol water content (step 1+step 2) is calculated. If the actual current aerosol water content is greater than half this value, the aerosol is taken to be wet and the equilibrium water 9946 content is distributed over all modes in proportion to their sulfate and sea salt concentrations. Otherwise, the aerosol is taken to be dry. The RH-deliquescent point is set to a constant value of 80% (appropriate for ammonium sulfate, following Ghan et al., 2001). The critical RH is set to a constant value of 35% (appropriate for ammonium sulfate, following Ghan et al., 2001) for the hysteresis calculation. As the water associated specifically with sea salt is not a tracked variable, this treatment is for the aerosol phase water as a whole.

Incorporation of sulfate produced by in-cloud aqueous-phase oxidation
The sulfate produced by in-cloud oxidation of SO 2 is distributed over a selected set of modes such that each activating particle receives an equal portion of that sulfate: where ζ i is the fraction of particles in mode i that activate to form cloud droplets, N i the particle number concentration in mode i , and ∆[H 2 SO 4 ] cond is the concentration of sulfate produced in-cloud during the time step. MATRIX provides the host model cloud modules with either the number concentra-15 tion of soluble particles for each mode, an estimate of the number of activating particles based on constant activating fractions for each mode, or an estimate of the number of activating particles based on the aerosol activation parameterizations of Abdul-Razzak and Ghan (1998Ghan ( , 2000. For the number concentration of soluble particles, the modes treated soluble are selected by setting the soluble (or activating) fraction of particles κ i 20 for each mode. In this study the fraction was set to unity for modes ACC, DS1, DS2, BC2, BC3, BCS, OCS, SSA, SSC, SSS, and MXX, 0.7 for mode OCC, and zero for all other modes. The droplet activation parameterizations of Abdul-Razzak and Ghan (1998Ghan ( , 2000 treat multimodal and multicomponent aerosols. This parameterization provides the Interactive Discussion mode composition and the updraft velocity. The mode-average hygroscopicity parameters are key in this treatment. The composition of the soluble components of mineral dust is treated as in Ghan et al. (2001). Modes that are thought of as insoluble cores with soluble shells, such as the BC and dust modes, may have a small soluble mass fraction that is effective at enhancing activation, and bulk treatment of these parameter-5 izations may underestimate this enhancement. All modes are included in the activation calculation and each mode may have a nonzero activated fraction. The number concentration of activating particles for each mode is saved for use in the cloud modules, where aqueous chemistry, scavenging of interstitial particles (due to Brownian diffusion of particles to droplets), collision-coalescence and evaporation, aerosol resuspension, 10 and wet removal is calculated.

Coagulation
A derivation of the intermodal coagulation production terms P where g i kl ,q is unity if coagulation between modes k and l adds species q to mode i and is zero otherwise, K kl is a mode-average coagulation coefficient (defined in the Appendix B), δ kl is the Kronecker delta, and m j,q is the mean mass of species q per particle in mode j . Modes k and l are distinct modes, and mode i may be a third 20 distinct mode or the same mode as k or l . The loss terms are 9948 The coagulation coefficient K i j is the sum of five contributions: Brownian coagulation, the convective Brownian diffusion enhancement, gravitational collection, turbulent inertial motion and turbulent shear. The kernels are described in Jacobson (2005), and the values shown in Fig. 15.7 ofJacobson (2005) were reproduced by the MA-TRIX model. The collision efficiency for the gravitational collection was taken from 5 Eq. (12-78) of Pruppacher and Klett (1980). The dependence of K i j on temperature and pressure was examined by calculating K i j for all pairs of particle sizes for a set of diameters ranging from 0.003-30 µm at several values T and p. For temperatures of 288, 200, and 325 K at a pressure of 101325 Pa, the ratio K i j (200 K) / K i j (288 K) ranged from 0.60 to 1.4, and the ratio K i j (325 K)/K i j (288 K) ranged from 0.90 to 1.2, 10 showing a rather weak temperature dependence for K i j . The Brownian K i j is proportional to the sum of particle diffusion coefficients and thus depends on the Cunningham slip-flow correction to those coefficients, the particle Knudsen numbers, the mean free path, and thus the ambient pressure. For pressures of 101325, 10132.5, and 1013.25 Pa at a temperature of 288 K, the ratio K i j (10132.5 Pa) / K i j (101325 Pa) ranged from 15 0.41 to 8.0, and the ratio K i j (1013.25 Pa) / K i j (101325 Pa) ranged from 0.29 to 46, showing a significant pressure dependence for K i j . The T and p dependence on K i j is included through lookup tables for the mode-average coagulation coefficientsK i j with table dimensions for T , p, and each of the lognormal parameters D g,i , D g,j , σ g,i , and σ g,j for modes i and j . Table temperatures were 200, 260, 325 K and pressures 20 were 101325, 10132.5, and 1013.25 Pa. The D g,i and D g,j each spanned the interval [0.003, 30.0] µm with 81 evenly-spaced (log scale) values, and for σ g,i and σ g,j 1.6, 1.8, and 2.0. For each pair of coagulating modes, σ g,i and σ g,j are fixed, and T , p, D g,i and D g,j are each linearly interpolated from the tabulated values. Precise mass conservation is accomplished by rescaling all mass variables for a specie at the end 25 of each time step as needed to remove minor inaccuracies in the approximate analytic solutions when large time steps are used. These inaccuracies disappear as the time step or coagulation rates are reduced. As newly-formed particles grow by condensation and self-coagulation, the mean diameter of the Aitken (AKK) mode may approach that of the accumulation (ACC) mode.
To prevent the AKK mode from becoming too broad and containing excessively large particles, particles are transferred to the ACC mode as the AKK mode diameter ap-5 proaches that of the ACC mode. This can be done in MATRIX in one of three ways. The first method is to transfer the fraction F rac from AKK to ACC during the time step, where D p,AKK is the diameter of mean mass of mode AKK, D p,AKK,min is the minimum value of D p,AKK and is the pre-selected size at 10 which newly-formed particles appear in the model (10-20 nm is a good choice for our model), and D p,ACC is the diameter of mean mass of mode ACC. The exponent p is an adjustable parameter in the range 2≤p≤6; the larger its value, the more intermodal transfer is delayed until the AKK mode closely approaches the ACC mode in mean size. As long as mode AKK is near the size at which new particles are formed, there is 15 essentially no transfer; if mode AKK has grown to the size of mode ACC, the two modes are fully merged and the AK mode is emptied and ready to receive freshly-formed particles. When there are modest differences in mean particle size between the two modes, the transfer is spread over several time steps. This approach meets the criteria given in Whitby et al. (2002) (3) mode transfer is treated by the fractional transfer method (Eq. 23) with p=4.

Box model experiments: coagulation
Several test cases are examined focusing on coagulation of donor modes A and B to produce a distinct receptor mode C. This process is symbolized by the key interaction A+B → C, although there are six distinct interactions: A+C → C, B+C → C, A+B → C, and self-coagulation within all three modes. In all cases a 24-h period was simulated with a 0.5 h time step, with the discrete model subdividing the time step as needed to 15 avoid the attempt to remove more particles from grid point i than exist at point i . Initial lognormal distributions are indicated as [N(cm −3 ),D g (µm), σ g (1)].  fraction of accumulation mode dust was smaller than that of mode ACC, but increased over the 6-h period as dust particles acquired sulfate. Particles remaining in mode BC1 acquired enough sulfate to yield a small activating fraction (0.0006 with mechanism 1) after 6 hours. When mode BCS was present it showed an activated fraction 30%; when absent there was greater population of mode BC2. The total particle number 20 increased slightly due to a small amount of new particle formation. Mechanisms 5-8 do not represent distinct modes for accumulation and coarse mode dust, yet the total numbers of activating particles for these mechanisms was very similar to those for mechanisms 1-4; however, combining modes AKK and ACC to form a single sulfate mode (mechanisms 4 and 8) led to a significantly higher number of activated particles.

The global climate model
The Goddard Institute for Space Studies (GISS) General Circulation Model (GCM) climate modelE Hansen et al., 2005) , 2007) to the tropospheric gas-phase chemistry scheme (Shindell et al., 2003). The sulfate chemistry scheme provides the aqueous sulfate production rate for liquid clouds and the H 2 SO 4 concentration for the MATRIX scheme.

20
The term P emis i ,q (see Sect. 2.2) is the mass concentration emission rate of species q into mode i . Number concentration emission rates E num i are derived from mass concentration emission rates assuming lognormal size distributions and particle densities. Lognormal parameters D g,E and σ g,E (Tables 1a) are derived from Table 2  for the GISS interactive emission models. The emissions of the natural aerosols sea salt , dimethylsulfide (DMS) (Koch et al., 2006) and mineral dust (Miller et al., 2006) are calculated interactively in the model, depending on surface wind speed and other surface conditions. Sea salt emissions are calculated for two size classes, appropriate for the SSC and SSA mode. Mineral 5 dust emissions are calculated for 4 size bins, as described in (Miller et al., 2006), and are distributed into the DD1 and DD2 mode.
Most sulfur emissions are into gaseous phase sulfur dioxide (SO 2 ); we assume that 2.5% of the sulfur is emitted as sulfate particles. 99% of these direct sulfate emissions are assigned to the ACC mode, and 1% into the AKK mode. Industrial carbonaceous emissions are assigned to the BC1 and OCC mode, but carbonaceous emissions from biomass burning sources are assigned to the BOC mode.
Since here we wish to simulate present day climate, the trace gas SO 2 and ammonia (NH 3 ) (Bouwman et al., 1997) emissions are based on the anthropogenic emissions for 1995 from the Emissions Database for Global Atmospheric Research (EDGAR3.2) 15 (Olivier and Berdowski, 2001). Black (BC) and organic carbon (OC) emissions are from Bond et al. (2004). Biomass burning emissions for BC and OC are from the Global Fire Emission Database (GFED) model by van der Werf et al. (2003). Natural OM emissions are assumed to be derived from terpene emissions (Guenther et al., 1995), with a 10% production rate.

Transport and removal
Tracers, heat and humidity are advected using the highly nondiffusive Quadratic Upstream Scheme (QUS) of Prather (1986), which uses nine subgrid-scale spatial moments as well as the mean within each grid box. This increases the effective resolution of the tracer field and allows the GISS model to produce reasonable climate fields with 25 relatively coarse resolution.
Turbulent dry deposition is based on the resistence-in-series scheme described in Koch et al. (1999) and Chin et al. (1996) layer scheme of the GCM and depends on the aerosol mode mean diameter. Gravitational settling depends on the aerosol mode mean diameter and density and accounts for the effects of RH on density and size. The wet deposition schemes of the GISS modelE are described in Koch et al. (1999Koch et al. ( , 2006. The model has two types of clouds, convective and stratiform clouds. Tracer treatment in clouds follows the cloud pro-5 cesses, so that tracers are transported, dissolved, evaporated, and scavenged (with cloud water autoconversion and by raindrop impaction beneath clouds). This parameterizatons requires information about aerosol size and solubility, which are calculated for each aerosol mode by MATRIX. The solubility per mode is calculated by using a volume weighted approach, depending on the chemical composition of the aerosol 10 particles.

GCM simulations
The performance of MATRIX on the global scale is examined by comparing a simulation of the full MATRIX scheme (Mechanism 1) with a simulation that excludes the microphysical processes of nucleation, coagulation and condensation. This leads basically 15 to a simulation where all aerosols stay externally mixed as emitted, except particles that are emitted as internal mixtures (i.e. BC and OC from biomass burning). We compare these extreme cases, detailed mixed modes and externally mixed modes, to observations. As discussed in the previous chapters, MATRIX provides a wide choice of mode 20 configurations and parameterizations. In Sect. 3 all eight mode configurations were tested with the box model. On the global scale now we examine Mechanism 1, the most complex mechanism, using the following parameterizations: (1)  is simulated for 6 years under current climate conditions, and the mean over the last 5 years of the simulations will be analyzed in the following sections.

Mass simulation
Two experiments are carried out to examine the impact of micro-physical processes on the overall aerosol simulation. The experiment including the most complex mecha-5 nism 1 (BASE) will be compared with an experiment using exactly the same model setup, but excluding the micro-physical processes of nucleation, condensation and coagulation (NO-MIC). The NO-MIC experiment basically represents a mass only scheme, where only those aerosol modes are populated that experience chemical production or emission. For the current mechanism that will lead to population of the 10 following modes: ACC (will comprise all sulfate), BC1, BOC, OCC, DD1, DD2, SSA and SSC. Fixed size and solubility is assumed for each mode. Note that for NO-MIC there is no means to render the insoluble modes (BC1, DD1, DD2) soluble. Tables 3 and 4 present the budget per mass tracer and process for the BASE simulation. Figure 4 shows the annual mean column load concentrations per specie averaged 15 over the last 5 years of the simulations.
Including micro-physical processes leads to a global increase of sulfate mass by 22% in the BASE simulation compared to the NO-MIC case. NO-MIC is lacking sulfate nucleation, condensation and coagulation processes. When those processes are neglected less sulfate mass is calculated, which would mainly come from nucleation 20 and condensation (see Table 3). Nucleation adds small Aitken-mode particles into the atmosphere, leading to an overall longer lifetime of sulfate mass. This explains why the changes mainly occur in remote regions. Freshly nucleated particles tend to be engaged in the coagulation and condensation processes, which leads to a wide spread of sulfate material over aerosol mode and size ranges. 25 The changes seen in the nitrate concentrations simply complement the sulfate changes. More sulfate in the BASE case leads to less nitrate, as sulfate and nitrate compete for available ammonia in the system. Nitrate is reduced by 30% in the BASE 9956 important for reducing concentrations in remote regions. The dust load is increased by 5% in the BASE simulation. Including micro-physical processes, leads to decreases dust loads in remote areas and increases dust load close to the sources. The removal in remote regions is explained by the fact that aged dust has a higher solubility, and therefore has higher wet removal rates which 10 shortens lifetime and travel distance. The increase in the source region occurs because dust in the MXX mode has a smaller mean particle size than the pure DD1 or DD2 modes, which leads to less dry removal, especially gravitational settling of the very large particles in the source regions.
The changes in sea salt mass concentrations are rather small. Sea salt mass in- 15 creases by 6% without microphysics, again caused by the shift in size distribution. As number and mass are independent, the mean particle size is not constrained to appropriate or even physically reasonable values. Hence, these distributions of number mean diameter reveal model characteristics not to be taken for granted. For the coarse modes, MATRIX tends to transfer mass out of the DS2 and SSC modes too quickly, 20 but not number, with a tendency to predict mean particle sizes that are too small. As the SSC (and SSA) mean particle dry masses are fixed (so that number can be prescribed for these modes), the mean diameter distributions are sharply peaked for these modes. However, for mode DS2, there is a large spread in mean diameter extending to rather small sizes. This is not seen in mode DD2 as the mean diameter of this mode 25 is strongly influenced by emissions at a fixed mean particle size.
The evaluation of the mass concentrations for the BASE simulation is presented in Fig. 5 in Costa Rica (Schwarz et al., 2008) and two flights over Houston Texas in November 2004 (Schwarz and et al., 2006). The comparisons between model and observations of black carbon mass vertical profiles are presented in Fig 6. The observed concentrations show a sharp gradient in concentrations between boundary layer (PBL) and free tropospheric air over Costa Rica. The model predicts good boundary layer heights, but 25 underestimates BC in the PBL, especially in summer, and overpredicts BC in the free troposphere. Black carbon in Costa Rica mostly comes from biomass burning sources, which has large temporal variability. The model monthly average cannot capture the precise variability in the aircraft observations. The overestimated BC concentrations 9958 in the free troposphere over Costa Rica and Texas may be cause by inefficient cloud scavenging. The predominant aerosol mode in both regions is BOC. In both cases OC dominates the mode fraction, making those particles hydroscopic. However this still does not lead to sufficient wet removal. An interesting aspect of the BC aircraft measurements is that black carbon coatings 5 were also detected. Schwarz and et al. (2006) present (see Schwarz and et al. (2006), Fig. 6) the number fractions of internally mixed BC particles, ranging form 0.2 to 0.8. As shown in Fig. 6 Fig. 7. Total BC mass concentrations remain very similar to the original experiment, but the predominant aerosol mode is now BCS. This demonstrates the importance of the emitted aerosol size distribution.

20
The main mixing process among aerosols is coagulation, which strongly depends on particle size. Shifting OC particles into a larger size class leads to less mixing of BC and OC. Figure 8 shows the column load and zonal mean distribution of BC, as sum over all BC modes, for the BASE experiment, and the differences between BASE and Exp. 2. 25 On a global averaged, BC mass is decreased by 10% in Exp. 2. A decrease of BC mass is simulated all over the Northern Hemisphere and is most pronounced over Eastern Asia. In the BASE case, most BC is mixed with OC, whereas in Exp. 2 most BC is located in the BCS mode, therefore predominately mixed with sulfate.

Number and size simulation
Aircraft observations from various campaigns (flight tracks are displayed in Fig. 9) and station data from aerosol supersites of the Global Atmosphere Watch program (http://wdca.jrc.it), are used for the evaluation of modeled aerosol number and size distributions.

10
The PEM Tropics-B campaign Raper et al. (2001) took place over the South Pacific ocean in spring 1999. Aerosol size distributions were measured in 30 size bins by two aerosol probes on board the DC-8 plane. Aerosol number concentrations spanning from 0.1 to 3 µm particle diameters as measured by the PCASP, and for particles of diameters between 0.34 and 20 µm by the FSSP aerosol probe. In order to com- 15 pare those observations with the model, model data are interpolated to the flight tracks and monthly mean vertical profiles are calculated. Figure 10 shows the modeled and observed data for fine and coarse mode aerosols and indicates the simulated contributions from the individual modes to the total aerosol number concentration.
The differences between the two aerosol probes are quite large. The model simula-20 tion falls between the two measurements for particles sizes smaller than 1 µm. Here most particles come from the AKK, followed by the ACC, OCC, BOC and BCS modes. Close to the surface, the coated black carbon modes BC1, BC2 and BC3 show some small contribution to the total number concentration. Sulfate, nitrate and ammonium are the dominant species in this distribution, as they populate the AKK and ACC mode. 25 The two aerosol probes show very different number concentrations for particle larger than 1 µm. However the model underestimates the number concentrations of those larger particles, which mainly come from the MXX and SSA modes. 9960 The TRACE-P Jacob et al. (2003) campaign took place from February to April in 2001 over the North Pacific ocean. Vertical profiles of mass and number concentrations are shown in Fig. 11 and 12 as measured on the DC-8. In February the air craft flew between 150E and the US West coast sampling much cleaner air than during March and April, when Asian pollution dominated the shown profiles. Figure 11 shows an excellent comparison between modelled and observed SO 2 and a good comparison to sulfate mass concentrations. Nitrate and ammonia tend to be underestimated by the model. The measurement data show the same order of magnitude of sulfate and nitrate concentrations in the polluted profiles during March. Number size distributions during TRACE-P were measured in six size bins, spanning from 0.1 to 1550 µm. 10 Number concentration of particles lower than 1 µm are simulated reasonably well by the model (see Fig. 12) and the aerosol number concentrations in the coarse mode are underestimated. The precedence of modes is the same as over the South Pacific as seen from the comparisons to the PEM-Tropics-B campaign.
The campaigns discussed so far mainly sampled marine air masses. In contrast, 15 the INTEX-A (Singh et al., 2006) campaign took place over the North American continent and to a smaller extent over the North Atlantic and Pacific Oceans. The mass observations (Fig. 13) show a good simulation of nitrate close to the surface, but lack of modelled nitrate at higher altitudes. Model sulfate aerosol mass seems reasonable in the lower atmosphere, but overestimated in higher levels. The comparison to num-20 ber concentrations shows good agreement for particles lower than 2 µm and a lack of particles in the coarse mode, especially in the free troposphere. Observations are given in six size bins spanning from 0.3 to 2 µm. For these conditions, different modes are important. The BOC mode clearly dominates, followed by OCC and ACC in the boundary layer and the MXX mode in the free troposphere. The coarse mode shows 25 again mostly contributions from SSC and MXX.
Station measurements provide the unique possibility to compare model data to longterm observations. Unfortunately only few data sets providing detailed particle size resolution are available. Here we show comparisons to three European stations, pro-9961 Introduction vided by the Global Atmosphere Watch program (http://wdca.jrc.it/). Hohenpeißenberg is a mountain station at 988 m altitude located in southern Germany. Aerosol size distributions were measured from 0.1 to 6 µm. The year to year variability is very low, therefore we only compare data to the year 2002, as this was the most complete data set available. Figure 15 shows the annual mean size distribution 5 as observed at Hohenpeißenberg. The model data are interpolated into the observed size bins. The smallest size bin starts at 0.1 µm. That is the bin in which the highest number concentrations are observed, however even higher concentration might occur at smaller particle sizes. The model captures the size distribution fairly well, but overpredicts the number concentration of particles greater than 1µm. The two graphs in 10 the top panel of Fig. 15 show us the mode and species contributions to this distribution. Particles smaller than 1µm come from the AKK, and the coated black carbon modes (BC1,BC2,BC3). Dominating chemical species are sulfate and nitrate for the aitken mode particles and black carbon for particles around 0.05 µm size. The part of the size distribution that is observed and agrees well with the model, particle sizes 0.1 to 1 µm, have contributions from OCC, BOC, BCS and ACC, and to a limited extent contributions from the dust modes DS1 and DBC. The dominating chemical species in this part of the size distribution are sulfate, nitrate, organic and black carbon. Particles larger than 1µm belong to the mixed mode MXX, and comprise dust, sea salt mixed with sulfate and nitrate.

20
Pallas is a station located in northern Finland, at 560 m altitude. Size distribution measurements are available from 0.007 to 0.4 µm (see Fig. 15). The model captures the number concentrations at the lower end of the size spectrum. Those particles get formed through nucleation events, belong to the AKK mode and are comprised mostly of sulfate and nitrate. Particle number concentrations of around 0.05 µm size 25 are under-predicted by the model, which is primarily black carbon, and model and observations again agree pretty well for particles above 0.07 µm. Furthermore, the decreasing slope of the size distribution is captured by the model.
Alesund is located in Spitzbergen, even further north of the Finish station Pallas.

9962
The observations show somewhat smaller number concentrations in the nucleation mode, but a maximum in the distribution around 0.2 µm. The model underestimates the number concentrations at this remote station, where black carbon, sulfate and nitrate dominate the size distribution below 0.1 µm. Above that size organic carbon is present, and the agreement to observations is much improved.
5 Figure 16 presents the comparison of the three stations with model results from Exp. 2, where OC and BC emission sizes were doubled. It is very interesting to see how the prior visible gap in Fig. 15 around 1 µm particle sizes is now closed, at the expense of underestimating particles around 0.1 µm.

10
MATRIX has been successfully implemented into the GISS climate model and all 8 MATRIX mechanisms (see Sect. 2.) are functional in the global model. Although not shown here, the performances of mechanisms 1-4 show comparable results as the box model study. Especially the lag of the aitken mode AKK in mechanism 4 leads to significant changes in the overall model performance. Compressing the dust modes as 15 done in mechanism 5-8 further leads to a decrease in accuracy in the coarse mode model performance. Therefore from the global point of view mechanisms 1-3 are useful options for global applications, although model modes need further optimization.
The novelty of MATRIX is adding aerosol number, size and mixing to the mass concentrations in the GISS climate model. The strengths and weaknesses of the MATRIX 20 application on the global scale can be summarized as follows: -Aerosol size distributions are well represented by the model for Aitken and accumulation mode particle sizes. This was demonstrated by the comparisons to aircraft as well as the station observations. The evaluation of aerosol size and number concentrations demonstrated a very successful simulation regarding number 25 concentrations, below 1 µm. are important for aerosol mixing processes, as the coarse particles provide large surface areas for chemical and physical interaction. On the other hand the lifetime of coarse particles is very short, therefore their climate impact questionable.
-The most efficient process that leads to an internally mixed aerosol population is coagulation. Coagulation depends strongly on the size of the involved particles. 10 However, there is limited information in emission inventories about the size of the emitted particles. This information is extremely important for the results of microphysical models.
-Testing the design of the MATRIX mode concept on the global scale showed that the definitions of the modes can still be optimized. Some modes are sparsely 15 populated, for example the coated dust and black carbon modes, whereas fine aerosols tend to accumulate in the the BCS, the black organic carbon-sulfate mode and coarse mode aerosols accumulate in MXX, the mixed mode.
-Nitrate aerosols are calculated only as a mass based scheme, and the distribution over the modes is scaled relative to the sulfate distribution, as discussed in 20 section 2. (The same concept applies for ammonium and aerosol water). Therefore the condensation of nitric acid on particles surfaces is not yet explicitly taken into account in the model formulation. This can lead to a lack of nitrate on coarse aerosol particles and explains the poor nitrate correlations in the remote atmosphere and higher altitudes, as demonstrated in Bauer et al. (2007); Bauer and cesses. Furthermore sulfate material now spans the total size range of atmospheric particles. This will have an important impact on the calculation of the aerosol direct and indirect effects.

Conclusions
This manuscript gives a detailed description of the microphysical model MATRIX and 5 also describes its performance as box and global models, as part of the GISS climate model. The climate model including MATRIX takes 6h computing time to simulate 1 model year, using 15 processors on a Linux Networx cluster. Therefore this model system is very suitable to perform many sensitivity studies as well as perform transient climate simulations. MATRIX provides a wide set of possible mode configurations 10 and microphysical parameterizations. It is beyond the scope of this first publication to present, analyze and evaluate all options of the MATRIX model. However, this paper gives the complete model description, which will serve as reference for all future studies. The MATRIX code is written as a box model and should be readily adaptable within other regional or global model frameworks. MATRIX is based on the quadrature 15 method of moments (QMOM) scheme and the version presented in this paper includes two moments, number and mass, and one quadrature point. The two moment scheme was chosen as a starting point to develop the microphysical package, however our intention is to use higher aerosol moments in the future.
When designing a new model, decisions about what parameterizations to include or 20 how to represent the aerosols will have large impacts on the end product, the impact of aerosols on the earth system. MATRIX tries to be as flexible as possible at this point, by including various sub-models, and allowing a wide variety of different mode configurations. The presented box model results demonstrated the large impact mode configuration can have on cloud activation. However the simulations on the global scale showed that some modes are sparsely populated, due to the tendency for aerosols to mix and to include more than two chemical components. For the Aitken and accumula- tion mode particles, this leads to a dominance of the black, organic carbon and sulfate mixture, and in our formulation is allowed to include nitrate as well. In the coarse mode this favors the mixed mode, which allows all chemical species to be included in the aerosol. Our future work will focus on improving the simulation of aerosol mixtures, which will be an ongoing and difficult task. Detailed aerosol measurements are neces-5 sary to achieve this goal. The coupling of the MATRIX aerosol schemes to radiation and cloud schemes is work under progress. (2) A more physically-based expression taking into account the small-particle dynamics (described) is the analytic formula for F is derived 20 in Kerminen and Kulmala (2002) and Kerminen et al. (2004). It is given by where the time t' at which the new particles appear and is later than the time t at which nucleation occurred on account of the time required for growth, and the parameter η 9966 (nm) is calculated as γk c,α=1 /4πg r D H 2 SO 4 with γ (nm 2 m 2 h −1 ) a proportionality factor, k c,α=1 (s −1 ) the condensation sink (see Eq. 18) calculated with the mass accommodation coefficient α set to unity and the mean free path of air (λ air ), g d (nm h −1 ) the diameter growth rate d D p /d t, and D H 2 SO 4 (m 2 s −1 ) the diffusivity of H 2 SO 4 in air. The formula is applicable when the pre-existing aerosol and condensable vapor concen- Thus γ is nearly constant over the time step when these conditions are satisfied. If conditions are approximately steady-state, the difference in times t' and t is immaterial. Steady-state assumptions are discussed below. If η exceeds 56, F is set to 10 −17 .
The diameter growth rate g d of freshly nucleated particles is calculated as that due 15 to condensation of inorganic vapors (H 2 SO 4 , NH 3 , H 2 O) only (contributions from lowvolatility organic vapors in not yet represented) and is that of free molecular growth given by  15 and A NH 4 HSO 4 is the same as A (NH 4 ) 2 SO 4 on account of lack of data for ammonium bisulfate solutions. The first factor is due to the T dependence of surface tension of the solution (taken as that for water), the second factor is due to the dependence of the partial molar volume of water on h, and the third term is due to the dependence of the surface tension of the solution on h. Once χ S is obtained, then n H 2 O =(1−χ S )/χ S . A2 Calculation of the sulfate mass in a newly formed particle The sulfate mass in a newly formed particle of ambient diameter D npf as a function of RH and neutralization regime is now described. For each regime, and for D npf selected as either 3, 10, or 20 nm, the diameter ratio D npf /D npf , dry is calculated from polynomial fits and used to convert the ambient particle volume (πD 3 npf /6) to dry particle volume, 5 from which a sulfate (SO −2 4 ) mass can be calculated (using appropriate densities and molecular weights for each regime). The Kelvin effect is taken into account. The calculated sulfate mass at a given RH is stored in a lookup  (2002). This is a significant amount of aerosol, and hence smaller aerosol loadings and values of k c leading to τ>∆t will likely occur at a significant frequency in atmospheric models. When ∆t<n τ τ the steady-state assumption is not invoked and the equation and is used to calculate the nucleation rate and g d . This expression neglects the consumption of H 2 SO 4 through the formation of new particles and gives an upper limit 10 for the actual [H 2 SO 4 ] ∆t/2 . Eq. A7 shows that for t=∆ t=n τ τ, [H 2 SO 4 ](t) is within [100/(n τ +1)]% of its steady-state value when only production and condensation are considered, and since inclusion of NPF reduces τ and makes the approach to steadystate more rapid, the choice n τ =2 leads to application of the steady-state approximation under conditions for which the sulfuric acid concentration will rise to within 33%

Derivation of production and loss terms for coagulation
For the coagulation of mode k with mode l to produce particle number or mass (or neither) in mode i , there are three cases. The first case occurs when k=l =i , which 5 is self-coagulation within a single mode. The second case occurs when k = l = i , in which the distinct modes k and l coagulate to form a third distinct mode i . The third case occurs when either k=i or l =i (but not both) in which mass and number are lost from mode k or l and mass but not number is gained in mode i . Before discussing these cases, we define the mode-average coagulation coefficient or order n as for number size distributions n k and n l for modes k and l , respectively, and particle diameters D 1 and D 2 . Note that for n =0,K kl is not symmetric in k and l , and for n =0 the first index labels the donor mode (k or l ) providing concentration of diameter moment n to the receiving mode (i ). 15 The present model tracks number and mass concentrations only. The rate of change of number concentration N i in mode i is expressed as production and loss terms as and the rate of change of mass concentration Q i ,q of species q in mode i is similarly expressed as Case 2: For the second case, mode k coagulates with mode l to add new particles and mass of species q to a mode i that is neither mode k nor mode l . Each coagulation 5 event removes a particle from each of modes k and l and adds a particle to mode i . For the number concentration production rate in mode i The normalized third moment (m 3 aerosol per particle) for mode k can be expressed as 15 and combining Eqs. (A6) and (A7) The mass concentration Q k,q of species q (µg q m −3 air) in mode k is related to the third moment as with V k the total aerosol volume concentration (m 3 aerosol m −3 air) of mode k, ρ k,q the mass of q per volume of particle (µg q m −3 aerosol), m k,q the average mass of q per 5 particle (µg q particle −1 ) in mode k. In the parentheses following the third equality we have the ratio of units (µg q particle −1 ) / (m 3 aerosol particle −1 ), which is (µg q m −3 aerosol) as required. Combining Eqs. (B8) and (B9), for the mass concentration loss of q in mode k kl m k,q (B10) 10 and for likewise for mode l The production term for mode i then balances these losses l k m l ,q ) (B12) Case 3: For the third case, mode k (or l ) coagulates with mode i , reducing number 15 and mass concentrations in mode k and increasing the mass concentration of q (but not number concentration) in mode i . Each coagulation event removes a particle from mode k, and using previous results we have for the loss of number concentration in mode k and for the loss of mass concentration of q in mode k.
For the mass production of q in mode i The expressions for production and loss rates obtained in these three cases are now 5 combined and summed over modes as appropriate. For the production terms for number concentration in mode i , only the second case is nonzero, and summing over all pairs of modes k and l where d i kl is unity if coagulation of modes k and l is defined to produce particles in 10 mode i and zero otherwise, and there are n modes. Recall that in the second case mode i must be distinct from both k and l . The summand is symmetric in k and l , and the summation should include either k > l or l > k, but not both. For the loss terms for number concentration, there is the self-coagulation term and all pairs of modes whose intermodal coagulation results in the loss of particles from mode i where d i j (not symmetric in i , j ) is unity if coagulation of mode j with mode i results in the removal of particles from mode i and zero otherwise. For the production terms for mass concentration of q in mode i , only in the second and third cases are nonzero and are combined as where the Kronecker delta δ ki ( δ ki =1 for k=i , δ ki =0 for k = i ) omits the k term since if k=i , mode k does not add mass to mode i , and similarly for δ l i . For the first case above, δ ki = δ l i =0 ; for the second case, δ ki =0 and δ l i =1. Summing over all pairs of modes k and l , the production term for mode i is where g i kl ,q is unity if coagulation of modes k and l is defined as producing particles in mode i and either mode k or l is defined to contain species q; it is zero if either of these conditions is not met. For the loss terms for mass concentration of q in mode i where d i j (not symmetric in i , j ) is unity if coagulation of mode j with mode i results in 10 the removal of particles from mode i and zero otherwise (as defined above).