Development of an online radiative module for the computation of aerosol optical properties in 3-D atmospheric models: validation during the EUCAARI campaign

Obtaining a good description of aerosol optical properties for a physically and chemi-cally complex evolving aerosol is computationally very expensive at present. The goal of this work is to propose a new numerical module computing the optical properties for complex aerosol particles at low numerical cost so that it can be implemented in 5 atmospheric models. This method aims to compute the optical properties online as a function of a given complex refractive index deduced from the aerosol chemical composition and the size parameters corresponding to the particles. The construction of look-up tables from the imaginary and the real part of the complex refractive index and size parameters will also be explained. This approach is vali- 10 dated for observations acquired during the EUCAARI campaign on the Cabauw tower during May 2008 and its computing cost is also estimated. These comparisons show that the module manages to reproduce the scattering and absorbing behaviour of the aerosol during most of the ﬁfteen-day period of observation with a very cheap computationally cost.


Introduction
While the greenhouse effect on global warming is quite well understood and leads to a quantification of global temperature increases, the effects of aerosol particles on the radiative budget of the atmosphere are still modelled only roughly (Forster et al., 2007). The direct, indirect and semi-direct effects of aerosols may have opposite impacts on rect forcing and feedback. The current methods used to quantify the optical properties are not accurate enough to take the high variability of aerosol composition and size distribution into account.
At those resolutions, it is imperative to consider the ageing of the aerosol and so, the evolution of the size distribution and the chemical composition of the aerosol particles. 15 For example, the absorptivity of primary absorbing carbonaceous particles is due to the coating of carbonaceous aerosols by secondary hydrophilic particles (Jacobson, 2000;Mikhailov et al., 2006). In parallel, the increase of the particle size also impacts the scattering and absorbing power of the aerosol (Seinfeld and Pandis, 1997). In that sense, it is necessary to have a good description of both chemical and physical 20 aspects. Modelling the chemical evolution of aerosols and computing the aerosol optical properties as prognostic variables is numerically very expensive. Because of this high cost, which results from an iterative computation of particle optical properties, the computation of aerosol optical properties as a function of its size distribution and chemical composition is never performed within meso-scale atmospheric models nor within 25 weather models. In several models, the impacts of aerosols on radiation is parameterized either by climatology tables or as a function of the main compounds concentration deduced from the literature (Kinne et al., 2006). In Grini et al. (2006), the effect of dust particles on radiation is studied with an elvolving aerosol size distributions, but the Introduction aerosol is supposed to be only one kind of chemical composition. We propose a more complex approach, with a number of aerosol compounds taken into account higher than in usual methods (Lesins et al., 2002). In this context, the goal of the present work is to set up a computationally cheap module depending on a complex physical and chemical aerosol description to compute 5 the optical properties at six short wavelengths.
The EUCAARI field experiment (Kulmala et al., 2009) which took place during May 2008 near Cabauw (Netherlands) provided us with a complete data set of aerosol measurements for the validation of the work. The Cabauw tower was fully equipped with a set of instruments measuring various aerosol properties. The location of the 10 Cabauw site between the North Sea and the industrialized area of Rotterdam allowed different aerosol types, from polluted to maritime air masses, to be observed.
The first section will describe the bases used to compute the aerosol optical properties, and highlight the sensitivity of those parameters as a function of both physical and chemical aerosol description. Then the building of look-up tables in order to save 15 time during the computation will be explained, and the resulting optical properties will be evaluated in regard of Mie computations. Finally, the comparisons between observations acquired during the EUCAARI campaign at the Cabauw tower and computations from the module will be presented.

Computation of optical parameters 20
In order to compute the aerosol optical properties, we applied the Mie theory (Mie, 1908) assuming that the aerosol was made up of spherical particles. Relative to other theories that exist to compute the optical properties of a particle, the Mie theory remains the method that is the most commonly used at present. For a given complex refractive index K =k r +i ·k i and a given effective radius r at a wavelength λ, the Mie 25 theory computes the extinction, absorption and scattering, efficiency of the particle which are noted Q ext,abs,sca (r,λ), respectively, and the Legendre polynomials of the phase function. Because of the high cost of the explicit fluxes computation, the radiative models commonly use the following three optical parameters: the Extinction coefficient or b, the Single Scattering Albedo (ratio of scattering to extinction) or SSA, and the asymmetry factor, g, giving the radiative direction information. Considering the aerosol size distribution n(r), the optical parameters can be computed for a given radius r as: One of the first steps to apply the Mie theory to an aerosol particle with a complex chemical composition, is to define a refractive index for the whole aerosol.

10
It is important to consider an evolving refractive index for the aerosol. Indeed, the main objective of this study is to perform an online computation of the aerosol optical properties considering an aerosol evolving through its size distribution and its chemical composition. The purpose of this study is to set up a module for a large set of aerosol compositions. As discussed in Chylek et al. (2000), the Maxwell-Garnett mixing rule 15 suits the best for situations with many insoluble particles suspended in solution. Yet, the small scale studies about aerosol particles and with a high spatial variability usually take place in urban areas, which precisely show those conditions.
In order to define a refractive index corresponding to an aerosol particle composed of different chemical components, the Maxwell-Garnett equation (Maxwell-Garnett, 1904) 20 as defined in Tombette et al. (2008) allows us to link the chemical composition of the aerosol to a refractive index and then to take the particle size distribution into account. Introduction

The Maxwell-Garnet equation
This approach considers the aerosol as being made up of an inclusion and an extrusion.
The inclusion is composed of the primary and solid parts of the aerosol, whereas the extrusion is composed of the secondary and liquid parts of the aerosol. Then the 5 effective refraction index is: where i are the complex effective dielectric (square root of the refractive index) constants in which subscripts 1 and 2 stand for the inclusion and the extrusion.
f is the volumic fraction of inclusion.
In the computation it is assumed that aerosol particles are only composed of primary Organic Carbon (OCp), Black Carbon (BC), Dust, nitrates (NO − 3 ), sulphates (SO 2− 4 ), ammonium (NH + 4 ), water (H 2 O) and Secondary Organic Aerosols (SOA). The inclusion or core of the aerosol is then composed of the OCp, the BC and the dust whereas the extrusion or shell is composed of NO − 3 , SO 2− 4 , NH + 4 , H 2 O and SOA. 15 The refractive index for each component considered is as defined by Krekov (1993) andTulet et al. (2008) and shown in Table 1.

Sensitivity of optical parameters
Here, we conducted various sensitivity tests based on variations of the aerosol size distribution (Sect. 2.3.1) and chemical composition (Sect. 2.3.2). Although these re-Introduction proposed parametrization. In order to show the importance of considering a consistent size distribution and chemical composition for the aerosol, Mie computations were performed with different refractive indexes and different size distributions.

Sensitivity to size distribution
The size distributions were represented by a lognormal function and described by a me-5 dian radius and a geometric standard deviation. The effective radius, representing the predominant radius with respect to radiation, is expressed as r eff =r 3 v /r 2 s with r v , r s standing respectively for the volume mean radius and surface mean radius depending on the geometric standard deviation. Figure 1 shows that for a given refractive index of 1.55−0.1i, the increase of the 10 median radius and the geometric standard deviation has a direct impact on the mass extinction efficiency distribution. For the three geometric standard deviations considered: 1.55, 1.75 and 1.95, the evolution of the mass extinction efficiency is non-linear with respect to both the median radius and the effective radius, and can show a difference of 40%. The maximum value for each geometric standard deviation considered 15 stands for a different median radius and can reach 3.6 m 2 g −1 for a geometric standard deviation of 1.55 whereas it reaches only 2.5 m 2 g −1 for a geometric standard deviation of 1.95. Because aerosol size distributions are usually expressed according to the median radius and also because the non-linearity of the parameters stand as a function of both 20 effective and median radius, the results will be shown as a function of the median radius in the rest of the study, for ease of understanding. Figure 2 shows the extinction efficiency Q ext , drawn in black and the three size distributions corresponding to the same median radius of 0.045 µm and the three geometric standard deviation are overlaid in red blue and green. The mass extinction efficiency results from the integration of the 25 extinction efficiency multiplied by the size distribution (and a factor of πr 2 ). Figure 2 shows that the non-linearity of the extinction efficiency versus the median or effective radius mainly comes from the high fluctuations of the extinction efficiency up to 2 µm.

Sensitivity to chemical composition
The influence of the chemical composition of the aerosol is shown in Fig. 3. The aerosol is considered to have a fixed median radius of 0.15 µm, a fixed geometric standard deviation of 1.6, and a fixed total mass concentration of 10 µg m −3 . We calculated the aerosol SSA for different refractive indexes by using the Maxwell-

5
Garnett equation (see Sect. 2.2), considering an aerosol composed of 20 to 80% of core (primary aerosol), and 80 to 20% of a shell of secondary species composed with the following mass ratio: 0.75 NH + 4 , 1 H 2 O, 1.35 SO 2− 4 , 0.9 NO − 3 and 1 SOA to ensure the electroneutrality of the secondary solution.
The resulting SSA is computed for 3 different core compositions depending of the

Methodology for the building of the look-up tables
An analytical solution was used, employing a look-up table of aerosol optical properties and a mathematical analytical function approximating the Mie computation.
In order to minimize the number of stored terms, the construction of the look-up tables was adjusted in two different ways described thereafter.

Description of the polynomial interpolation method
The first way to minimize the number of stored terms was by approximating the optical parameters with a double 5th degree polynomial interpolation according to the the median radius. Considering the median radius as a stored input parameter would be consistent for no more than a hundred cases as pointed out by Grini et al. (2006). To 10 avoid this constraint, the computation and storage of the fifth degree polynomial coefficients that best fitted the optical parameters evolving according to the median radius was taken. In this way, the input terms of the look-up tables are the complex refractive index and the geometric standard deviation. However, in several cases as shown in Fig. 4, the shape of the Mie resulting parameter evolution was poorly reproduced by the 15 best fitted mono polynomial interpolation. The mean relative errors were around 200% and showed a difference of 1.5 m 2 g −1 at the maximum of mass extinction efficiency.
For this reason the interpolation of the curve was divided into two different curves defined by the evolution from the first radius considered to radius corresponding to the maximum of the curve r cut , and from r cut to the last radius considered. Then as shown Introduction To summarize, the polynomial coefficients stored corresponding to the polynomial P (r) follow the equation: P (r) = a 10 + a 11 r + a 12 r 2 + a 13 r 3 + a 14 r 4 + a 15 r 5 if r < r cut a 20 + a 21 r + a 22 r 2 + a 23 r 3 + a 24 r 4 + a 25 r 5 if r > r cut where r cut is defined as the median radius corresponding to the maximum value of the optical parameter, and the 12 a i ,j , with i = 1,2, j = 0 to 5 are stored in the look-up table.

5
The coefficients a i ,j were computed by a least square approach for the fitting of the two parts of the curve.

Range of parameters chosen
The input terms of the look-up tables are then the imaginary and real part of the complex refractive index of the aerosol, and the geometric standard deviation of the size 10 distribution. The second way to minimize the stored terms was by choosing of input terms in the look-up tables so as to obtain pseudo linearity between adjacent pairs of resulting optical properties. Tests were performed to select a minimum amount of stored terms (not shown here), and six k i , eight k r and eight σ were chosen because of the pseudo linearity between 15 pairs of terms. This selection method allowed us to interpolate between the results of each optical property corresponding to the two stored values. The stored terms of the input parameters as well as the 6 wavelengths are reproduced in Table 2. As an example, if the value of the geometric standard deviation is 1.52, the corresponding aerosol optical properties are the weighted interpolations between the optical properties corre-20 sponding to the stored geometric standard deviations of 1.45 and 1.65.
As in most cases, if we position ourselves between two stored k i , k r and σ we compute the eight polynomials corresponding for the eight different combinations and we interpolate with the appropriate weight for obtaining the weighted optical properties. Although the computation is performed at 6 wavelengths which may be considered Introduction as input parameters, there is no possible interpolation between the wavelengths as performed for the other input parameters.

Evaluation of the method in regard of direct Mie computations
In order to evaluate the dual polynomial interpolation method, a comparison between the resulting and the direct Mie computed aerosol optical properties was performed.

5
This comparison was done at the 6 previously explicited λ, for 100 considered values of median radii, and for 20k r , 20k i and 18σ equally spaced in order not to be located at the stored values, leading to 4 320 000 comparison points for the 3 optical parameters. Because the computed parameters stands for values of k r , k i and σ which require an interpolation by the module (they do not match to the stored values), this comparison 10 allows to evaluate the module in regard of both polynomial interpolation method and choice of input stored terms. Figure 6 shows the result of the comparison for the three optical properties between the dual polynomial interpolation module and the direct Mie application. On the top (respectively middle and bottom) is represented the mass extinction efficiency, single 15 scattering albedo and asymmetry factor, computed by the dual polynomial interpolation as a function of the same parameter computed from the direct Mie application. Due to the very high number of scatter plots (4 320 000 for each optical parameter), the comparison is shown as a 2-D density function computed over 100 equally splitted bins. For each bin, the 2-D density function is computed as a function of the number of 20 sampling included in the bin. The first coloured contour (dark blue) stands for 99% of the optical parameters computed by the dual polynomial interpolation included within the area. The 1-D density function of the Mie computed optical parameters are represented by the continuous red lines and shows to which optical parameter values do the most cases stand for. The 1-D density function is computed as a function of the total 25 number of sampling.
First, Fig. 6 shows that for the three optical parameter computations, the dual polynomial interpolation method manages to reproduce the values computed by the direct 745 Introduction Mie application. Indeed, almost all of the parameters computed by the module follow the linear curve x=y represented by a continuous black line. The correlation coefficients value standing for the three optical parameters are respectively 0.9992, 0.9946 and 0.9994. The most large differences stand for extreme cases with very low median radius or very high geometric standard deviation. As an example, the highest values of mass extinction efficiencies which can reach more than 15 g m −2 stand for very unusual combinations of parameters such as a geometric standard deviation below 1.4, a real part of the refractive index of 1.80. These extreme cases all occur on the first wavelength. We can also notice that for the single scattering albedo, the highest dispersion stand for very low values, meaning that the scattering is mainly absorbing. Moreover, 10 the most computed cases, represented by the continuous red line, occur where the density function is the highest, and at this points there is a good linearity between the two methods.

Note to obtain the look-up tables
The very large number of parameters inside the tables does not allow to write them in 15 the present article. However, the tables are available on request to the corresponding authors (benjamin.aouizerats@cnrm.meteo.fr, pierre.tulet@meteo.fr).

Validation during the EUCAARI campaign
For the validation of the module, a full set of measured data for the inputs and of the outputs was required. Concerning inputs, a description of the aerosol size distribution 20 giving the median radius and the geometric standard deviation for each mode, associated with the aerosol chemistry giving the aerosol composition was needed. Concerning outputs, the measurements of the aerosol optical properties was also necessary. During the EUCAARI field experiment, a complete set of instruments were available at the Cabauw tower. Between

Period and instruments
The data were processed to give one average value over 30 min for each parameter during the 15 days.
-The aerosol size distribution was deduced after merging the SMPS (model TSI 5   3034) and APS (model TSI 3321) size spectral observations: the SMPS measures the aerosol size distribution between 10 and 500 nm, and the APS between 500 nm and 20 µm. The measurements were considered to have been made in dry conditions (RH<50%). Then, a lognormal fit that best approached the observed distribution was found.
At each time step, a least square approach gave the median radius, the geometric standard deviation and the mass concentration for each aerosol mode observed by both instruments. To link the number of aerosol directly measured to the mass of aerosol, the hypothesis of spherical particles with a density of 2.5 g cm −3 was done.

15
-The chemical composition of the aerosol was deduced from the observations of an aerosol mass spectrometer (AMS, cTOF type) and a multi angle absorption photometer (MAAP model 5012 550 nm by measuring the Angström coefficient with an aethalometer. These measurements were also made in dry conditions.
A summary of the instruments and the associated aerosol properties is presented in Table 3.

Methodology and hypotheses
5 Some assumptions were made from the observations. The choice was to split the POM between primary organic carbon and secondary organic carbon with a 60-40% distribution. A second hypothesis was to split the aerosol compounds equally between modes. We also multiplied the composition concentrations observed by the AMS by a factor deduced from the ratio of masses, itself deduced from the fit of the SMPS+APS 10 size spectrum and the total of masses observed by the AMS and the MAAP. It was noticed that the AMS had a cut-off diameter of 500 nm. Then, a checking was also processed in order to see whether the masses observed by the AMS+MAAP were the same as those observed by the SMPS+APS integrated from 0.01 to 0.5 µm. 15 The 15-day period of the study showed different aerosol concentration level due to different regimes. As presented in Fig. 7, the aerosol composition gives us several items of information: from 15 to 17 May, the aerosol shows a continental composition; from 17 to 21 May, the low level of aerosols is characteristic of a scavenging period with a clean air mass and fresh aerosols; from 21 to 29 May, the aerosol comes again from 20 a continental air flux and the secondary fraction of the aerosol, in particular the inorganic species, also grows in proportion. Finally, the 29 May air flux was different and dusts were present during this period. The total measured mass concentration is also represented in Fig. 7 as a black line. Also one can note that the ratio of primary aerosol to secondary aerosol remains quite stable during the fourteen days of measurements.  Figure 8 shows the evolution of the aerosol size distribution during the fourteen days of measurement. This figure is the result of the fit applied to the SMPS+APS data and used as inputs in the module.

GMDD
The Fig. 8 shows the evolution of all the size distribution parameters versus time. The median diameter evolves as the geometric standard deviation. The maximum of 5 mass also evolves consequently as already shown in Fig. 7. Finally, the combination of data measuring size description and chemical composition gives a good representation of the period. Figure 9 shows the single scattering albedo at 550 nm as calculated form the neph-10 elometer and the MAAP measurements (continuous red line) and computed by the module (dashed blue line). Figure 10 shows the mass extinction efficiency at 550 nm as calculated from observations (continuous red line) and computed by the module (dashed blue line). There were no instruments able to measure the asymmetry factor during the EUCAARI cam-15 paign. Therefore, there is no comparison possible with the computed values.

Comparison of computed optical properties and observations
First, the main trend shows that both mass extinction efficiency and single scattering albedo are correctly reproduced by the module during the 15 days. The correlation coefficient for the 15-day period and for the mass extinction efficiency is 0.94, and 0.89 for the single scattering albedo.

20
Concerning the single scattering albedo, the module mainly manages to reproduce the highest as the lowest values observed, even with high temporal variability. The main error is the overestimation of the single scattering albedo on 20 May when the observations reach 0.5 whereas the computed value is 0.6.
Concerning the mass extinction efficiency comparison, again there is a good agree-25 ment between the observations and the computations. However, it is noticeable that from 20 to 24 May, the computations overestimate the mass extinction efficiency by 10 to 80%. To understand why there is a less good correlation especially between 20 and 749 Introduction  Fig. 11 shows the evolution of the difference between the observed and computed mass extinction efficiency (in red) during the 15-day period, and the evolution of the relative error on the integrated mass up to 500 nm observed by the SMPS+APS and by the AMS+MAAP. The discrepancy between observed and computed mass extinction efficiency between 20 and 24 May, may be directly linked to a major difference of masses calculated from by the two different experimental systems, the SMPS+APS for the aerosol size distribution and AMS+MAAP for the aerosol chemical composition. Between 20 and 24 May, the relative error between both systems on the aerosol integrated mass up to 500 nm can reach more than 100% and the major differences between the observed 10 and computed mass extinction efficiency occur at these same maxima. Thus, the main differences observed may probably come from a discrepancy in the instruments in term of mass measurement, leading to an inconsistency between the inputs and outputs of the module.

15
In order to quantify the direct and semi-direct effect of an aerosol, it is necessary to compute the aerosol optical properties in atmospheric models depending on an evolving and complex aerosol particle. To study the aerosol particles interaction with radiation in highly temporal and spatial variable locations, such as urban areas, it is necessary to consider a large number of chemical species.

20
Nevertheless, the high computing cost of the classically applied Mie theory generally leads to a climatological parametrization of aerosol optical properties in regionalclimate models.
The numerical computation was performed for six wavelengths using a vectorial computer (NEC SX8 type), in order to compare the time cost of this module with the other Introduction The comparison was made between two standard simulations, with and without computation of the aerosol optical properties. The standard simulation required a cpu time per point per time-step of 98.0 µs versus the standard simulation with the aerosol optical properties computation has a cpu time per point per time-step of 99.7 µs. These results show that the numerical cost of the module was no more than 1.7% of the stan-5 dard simulation total cost. As an example, the chemistry solving cost of these same simulations is 35.6%.
We can then consider that the previously described module is numerically very cheap relative to other processes and is consequently affordable for most atmospheric modelling. 10 This work presents a new, computationally cheap module dedicated to the online computation of optical properties according to the particle chemical composition and size distribution.
To minimize the computing cost, the module is founded on look-up tables built by a dual fifth degree polynomial interpolation of the parameters depending on the me-15 dian radius of the aerosol size distribution. The parameters are the geometric standard deviation of the lognormal size distribution of the aerosol mode considered, and the imaginary and real part of the complex aerosol refractive index corresponding to a chemical composition deduced from the Maxwell-Garnett equation.
The module was then evaluated by using observations acquired during the EUCAARI 20 campaign (May 2008). The numerical cost was also tested. The comparisons between optical properties modelled by the module and those acquired at the Cabauw tower during fifteen days of measurements shows a good correlation. Furthermore, the numerical cost of the module is shown to be very low, allowing its implementation in atmospheric models treating aerosol size distribution and chem-25 ical composition. This optical module is already coupled with the ORILAM aerosol scheme (Tulet et al., 2005(Tulet et al., , 2006 implemented into the atmospheric research model Meso-NH (Lafore et al., 1998). This development will be used in future works to investigate feedbacks of polluted aerosols on radiation and urban climate (especially impacts GMDD GMDD Tombette, M., Chazette, P., Sportisse, B., and Roustan, Y.: Simulation of aerosol optical properties over Europe with a 3-D size-resolved aerosol model: comparisons with AERONET data, Atmos. Chem. Phys., 8, 7115-7132, doi:10.5194/acp-8-7115-2008Phys., 8, 7115-7132, doi:10.5194/acp-8-7115- , 2008. 739 Tulet, P., Crassier, V., Cousin, F., Shure, K., and Rosset, R.: ORILAM, a three moment lognormal aerosol scheme for mesoscale atmospheric model.    Fig. 11. Relative error between the mass measured by the AMS+MAAP and that deduced from the SMPS+APS up to 500 nm (blue line) and difference between mass extinction efficiency observed and computed (red line).