G Vermot Des Roches

E Balmes
email: balmes@sdtools.com

R Lemaire
email: remi.lemaire@fr.bosch.com

T Pasquet

Design oriented simulation of contact-friction instabilities in application to realistic brake assemblies

Keywords: brake squeal, reduction, component mode, time simulation, instability, contact-friction

This paper presents advances in non-linear simulations for systems with contact-friction, with an application to brake squeal. A method is proposed to orient component structural modifications from brake assembly simulations in the frequency and time domains. A reduction method implementing explicitly component-wise degrees of freedom at the system level allows quick parametric analyses giving modification clues. The effect of the modification is then validated in the time domain where nonlinearities can be fully considered. A reduction method adapted for models showing local non-linearities is purposely presented along with an optimization of a modified non linear Newmark scheme to make such computation possible for industrial models. The paper then illustrates the importance of structural effects in brake squeal, and suggests solutions.

INTRODUCTION

Automotive brake design is nowadays oriented towards an optimized weight/performance ratio which tends to generate noisy systems. High friction coupling happening at the pad/disc interface is responsible for self-sustained instabilities in the auditive frequency range. The resulting noise can attain 120dB in the brake vicinity and is known as squeal between 1 and 16kHz or moan under 1kHz.

Unlike the case of low frequency vibrations, the brake performance is not altered by squeal, which happens mostly at low pressure, low speed conditions. The perceived quality is however altered, as the driver's feeling is disturbed and the environmental nuisance is not welcome. Methods to design silent brakes are mainly empirical and difficult to control, due to modelling issues (e.g. contact complexity) or implementation difficulties.

Classical design methods for brake vibration are set in the frequency domain, following experimental results on unstable mode lock-in patterns [START_REF] Massi | Dynamic and tribological analysis of brake squeal[END_REF] as function of global parameters such as the friction coefficient or the braking pressure. The system is linearized around a working point to apply Lyapunov theorem in which the system stability is related to damping of complex modes. This method shows great limitations and is not sufficient as only potentially unstable modes computed for a fixed contact state are found. In particular, recent work [START_REF] Lorang | Instabilité vibratoire des structures en contact frottant: application au crissement des freins de TGV[END_REF] showed that more unstable complex modes can be found than there actually exist.

In the time domain the system with its full non linearities can be simulated, which gives a clear view of brake stability. Such implementation raises many issues as a direct simulation on a full industrial model would require prohibitive computational costs. Contact handling requires relatively small time steps around 10 -6 s which makes long (100ms) simulations difficult to handle. Using a non linear implicit Newmark scheme on a 600,000 DOF system would actually generate over 1TB of data in over 700 hours.

These issues are dealt with model reduction techniques detailed in section 1.3 and an adaptation of the Newmark scheme to non linear penalized contact vibrations in section 2. Simulation costs thus become affordable yielding from 500MB to 5GB of data in 12 hours. The computational cost of time simulation is still too important to allow quick parametric studies, besides; from the design point of view, time simulation results must be exploited in a way allowing to close the loop with classical design parameters. It is then necessary to consider components and system characteristics, mostly known in the frequency domain. This paper thus aims at implementing in section 1.2 the Component Mode Tuning (CMT) method presented in [START_REF] Vermot Des Roches | Using component modes in a system design process[END_REF] to obtain quick design parametric studies from which design directions can be deduced.

A specific treatment of time simulations, based on shape identification presented in section 2.3 allows making the parallel between frequency and time methods. A design method is eventually suggested and illustrated for an industrial brake model provided by Bosch in section 3.

MODEL REDUCTION STRATEGIES

General concept

Simulation of large industrial models in the design and validation processes requires specific needs which do not match actual computational performance. Indeed, design studies are comfortable with very quick reanalysis simulations from a nominal model, in the hour scale. Validation simulations can take longer, with classical limit of acceptability around 10h (a night).

Although such systems as the industrial brake can be nowadays computed on workstations, and new real eigenvalue solvers like Automated Multi Level Solvers (AMLS) [START_REF] Kaplan | Implementation of automated multilevel substructuring for frequency response analysis[END_REF] allow to solve systems over a few million DOF, the time required scales in hours, not in minutes. The full real modes at the nominal state are however accessible, which opens the way to new reduction methods more adapted than the traditional Component Mode Synthesis (CMS) method [START_REF] Craig | A review of time domain and frequency domain component synthesis method[END_REF].

CMS was indeed based on the assumption of component independence, static solution capability and explicit boundary coordinates. Interaction information is thus a priori ignored, such that the full finite element basis of their interface is kept. The target application was however different from the aim of this study, which is to reproduce dynamic vibrations of an automotive brake working at a steady state.

Two strategies are presented in the following, which can be combined in the future. Section 1.2 features a reduction method introducing component-wise DOF at the assembly level so that the design loop between structural modification and reanalysis is optimized. Section 1.3 presents a reduction method adapted to models with local non-linearities, which is the case for brake squeal assemblies under the modelling considered for time simulations.

Explicit use of component level DOF

The approach considered here assumes that the components are physically disjoint and that coupling occurs through interfaces that have a physical extent. Since the components are disjoint, the associated DOF q ci are always distinct. The setup of such model is presented in figure 1.

[ ] [ ] [ ]       +       = + =
I I I I el el I el K K K K K K K K K
Figure 1 shows the matrix connectivities for a sample 3 component/2 interface model. The block structure of the component and interface models is indicated through colours. In the application of section 3, the interface matrices correspond to tangent contact/friction coupling models.

Since DOF associated with components are distinct, model reduction is simply performed by generating an assumed basis (classical Rayleigh Ritz reduction) T ci for each component response, thus yielding a reduction basis T R on the whole model [START_REF] Lorang | Instabilité vibratoire des structures en contact frottant: application au crissement des freins de TGV[END_REF] To obtain accurate predictions, one now seeks to retain exact modes for both the system and components. For each component i, one thus uses a reduction basis that combines component modes in free/free conditions φ ci and the trace on the component of exact system modes Φ |ci [START_REF] Vermot Des Roches | Using component modes in a system design process[END_REF] This method, fully described in [START_REF] Vermot Des Roches | Using component modes in a system design process[END_REF] allows the generation of very compact systems, as the size of each component block in the reduced matrices, is the size of the number of modes kept in T ci , that is to say a few hundred. The final brake assembly is thus reduced to around 1,300 DOF.

The reduction basis employed yields a reduced system showing identical dynamics, characterized by the real modes.

The matrix decomposition detailed in equation ( 1) once reduced on the basis described by equations ( 2) and (3) shows a very interesting topology, as the elastic matrix is now diagonal, each term corresponding to a component mode frequency. These frequencies are thus explicit in the reduced system and become a design parameter. Since the model is very compact, reanalyses from the reduced model as function of component modes is very quick (less than a minute for real mode computations).

Reduction of models showing local non linearities

Brake squeal dynamics rely on the quality of the pad/disc interface, as it is the location of the main instability. The idea for time simulations is then in a first approach to consider the remaining of the system as linear, based on the pseudo-periodic initial state.

All DOF in the vicinity of this contact area are thus kept explicitly, noted as q c . The remaining of the system DOF are noted q i .

To achieve the accuracy objective of exact dynamic behaviour, the same method as applied in section 1.2 is employed for the part to reduce, using the trace of the exact real modes of the assembly as assumed shapes. As only the q i DOF will be reduced, the trace (or restriction) of the Rayleigh-Ritz basis on this part is only considered. The reduction basis, illustrated in figure 2, is then expressed as

(4)
The pad/disc section is then kept unreduced, while all remaining parts are reduced in a superelement.

In the process, the interface DOF are implicitly reduced on the system real modes.

Figure 2 Industrial brake and model reduction strategy for non linear time computations

To cinematically generate a coherent coupling between the finite element part and the superelement, the reduction basis must in addition contain the pseudo-periodic stationary condition q 0i . The final time model features then 30,700 DOF, with a reasonable sparsity. This reduction method was detailed in [START_REF] Vermot Des Roches | Time simulation of squeal phenomena in realistic brake models[END_REF][START_REF] Vermot Des Roches | Design oriented time/frequency analysis of contact friction instabilities in application to automotive brake squeal[END_REF].

TIME SIMULATIONS FOR LARGE INDUSTRIAL SYSTEMS WITH CONTACT-FRICTION

Contact -friction formulations

Contact/friction modelling is commonly split into two formulation strategies, depending on the resolution strategy of the Signorini/Coulomb conditions, giving contact-frictions forces f n , and f t depending on the friction coefficient µ.

Contact can either be formulated as a structural constraint using a Lagrange formulation, or regularized using a penalization formulation. Both formulation strategies are based on the gap g between a defined set of contact points selected at

      = 2 1 0 0 c c R T T T [ ] . Orth ci ci ci T Φ = φ [ ]             Φ =       c R c i i c i q q I
q q q 0 0 0 6 th European Conference on Braking JEF 2010, Lille, France the interface, through the notion of gap between two surfaces, which is defined as the relative displacement along the contact normal N with a possible offset g 0 .

The penalization method chosen here aims at relaxing the contact constraint of Signorini by authorizing a controlled level of interpenetration. Practically, a relationship between the gap and contact forces is established to account for an approached contact constraint. where p is the contact pressure, g is the gap, p 0 and λ are parameters to define depending on the interface properties.

Friction implementation follows the definition of the Coulomb law, which for two solids relates the sliding velocity and the friction forces. A basic regularization, plotted in figure 3 and considered in the study, penalizes low sliding velocities (noted w S ) through the introduction of a parameter k t , such that if

n t S f k w ⋅ ≤ ⋅ µ (6)
else The pairing strategy and numerical implementation details can be found in [2,6,7,8],

Choice of a time integration scheme

The most common methods come from the Newmark family which implementation is well documented [START_REF] Géradin | Mechanical vibrations, Theory and application to structural dynamics[END_REF]. Simulation of instabilities requires using a well controlled damping strategy [START_REF] Hughes | Stability, convergence and growth and decay of energy of the average acceleration method in nonlinear structural dynamics[END_REF]. Modal damping is used in the frequency range of interest, which is detailed in section 2.3. Difficult to assess numerical damping implementations are avoided through the choice of an implicit non linear Newmark scheme [START_REF] Géradin | Mechanical vibrations, Theory and application to structural dynamics[END_REF] with scheme coefficients set to β=.25 and γ=0.5.

A non linear Newmark integration scheme is presented in figure 4. The classical method would recompute and factorize the Jacobian [J] at each time step for each correction increment. The most efficient solvers still need around 10s to factorize a matrix of a size of 30,700, which would take over 600 hours to perform, assuming one correction increment per step. The ability to pass such computation in a reasonable amount of CPU time lies in the use of a fixed Jacobian, which is more detailed in [START_REF] Vermot Des Roches | Time simulation of squeal phenomena in realistic brake models[END_REF]. To take contact forces into account, the Jacobian is set up with a mockup linear bilateral contact law which stiffness k c is parameterized, so that the Jacobian is written [START_REF] Vermot Des Roches | Design oriented time/frequency analysis of contact friction instabilities in application to automotive brake squeal[END_REF] where M, C, K are respectively the mass, damping and stiffness matrices, K C a linear coupling matrix of the normal DOF at the contact using a uniform contact law of unitary stiffness.

Design oriented identification of time response

Most common system characterizations are known in the frequency domain, which makes critical the ability to link non linear time simulations of a system to its modal properties -this is especially true for design applications. This section aims at creating this link using the system real modes, which can be used as a parameter through a time domain modal damping strategy or by identification in the response.

The most common and easy to use damping strategy in the frequency domain is modal damping. The concept comes from the assumption [START_REF] Géradin | Mechanical vibrations, Theory and application to structural dynamics[END_REF] that the structure modes are separated enough

g e p g p λ - = 0 ) (    ⋅ = = n t S t t f f w k f µ C c K k K C h M h J + + + = β γ β 2 1
Design oriented frequency/time simulation of industrial brake squeal not to interact with each other, ensuring a relatively good independence in terms of damping. [START_REF] Vermot Des Roches | Design oriented time/frequency analysis of contact friction instabilities in application to automotive brake squeal[END_REF][START_REF] Bianchi | Using modal damping for full model transient analysis[END_REF] have shown that the modal damping C M can be expressed by its application to a velocity vector v as [START_REF] Vermot Des Roches | Time frequency analysis of contact-friction instabilities[END_REF] which relates first the projection of the velocity on the real mode k by φ κ T Mv, the application of a damping ratio 2ζ k ω k then the restitution on the finite element basis by the projection Mφ k .

This concept can be extended to any shape ψ defined on the model DOF, provided it is mass normalized: ψ T Mψ = 1, which is always possible to obtain.

The modal damping strategy presented is based on the notion of projection of space fields on the system real modes, here the velocity. The same is possible for the displacement so that the projection of the time response on each real mode is possible and the response can be expressed as the sum of each real mode response and an orthogonal part [START_REF] Vermot Des Roches | Design oriented time/frequency analysis of contact friction instabilities in application to automotive brake squeal[END_REF]. Modal participations are computed by projecting the response q(t) as [START_REF] Géradin | Mechanical vibrations, Theory and application to structural dynamics[END_REF] It is then possible to compute the mechanical energy associated to each real mode E mk in the time response by computing [START_REF] Hughes | Stability, convergence and growth and decay of energy of the average acceleration method in nonlinear structural dynamics[END_REF] This response allows to evaluate and rank each mode response in the time response thus identifying the main ones. The computation of modal mechanical energies through equation ( 10) yields the notion of modal sensors.

INDUSTRIAL APPLICATION

The concepts presented in sections 1 and 2 are here applied to an industrial brake assembly provided by Bosch, presented in figure 2.

The application of the CMT allows to detect a priori unstable modes of the system along with their stability sensitivity as function of the component modes. The time simulation allows to detect modes which are the most unstable in realistic condition (i.e. regarding the full non linearity). This makes both methods complementary.

A cross application to evaluate the effect of a design change on a component is then suggested. In the following full system modes are designed with the prefix R for real modes and C for complex modes followed by the mode number, no prefix is used for component free/free real modes.

Stability predictions in the frequency domain

At 12 Bar, the brake system provided features several unstable modes, following the stability diagram of figure 5.

Figure5. Brake stability diagram at 12 Bar

Some unstable modes can be highlighted, in particular complex modes C44 and C51, which shapes are plotted in figure 6.

Figure 6. Shapes of unstable complex modes C44 (left) and C51 (right). Colors from blue to red ranking in ascending elastic strain energy per component

From the CMT, the strain energies by components can be extracted, so that the participations can be coarsely observed. Unstable mode C44 shows in particular a disc/outer-pad/caliper interaction.

More detailed information can be obtained by analyzing the sensitivity of the real part of each complex mode to component mode frequencies [START_REF] Vermot Des Roches | Using component modes in a system design process[END_REF]. Figure 7 plots the sensitivity results for mode C44 shown in figure 6.

∑ = = N k T k k k M Mv M v C 1 )) ( 2 ( φ ω ζ φ ) ( ) ( ) ( 2 2 2 2 t t t E k k k mk α ω α + =    = = ) ( ) ( ) ( ) ( t q M t t Mq t T k k T k k φ α φ α
The sensitivity amplitude shown in figure 7 allows to sort component modes thus obtaining the most sensitive ones. Reanalysis studies are preferred to perturbation studies as they are more precise, but require more computational time. Sensitivity computations then give an a priori guess to target specific modes in the analyses performed below.

Focusing on mode C44, it can be seen that the system is mostly sensitive to the caliper. This is a counter-intuitive observation as the friction induced instabilities happen at the pad/disc interface, thus highlighting structural effects.

More detailed reanalyses can then be performed by focusing on the caliper, especially on modes 12 to 14 that seem more relevant.

The shapes of caliper modes 12 to 14 are close, which will make difficult a physical single modification of mode 13. Its effect can however be studied numerically.

Figure 8 shows the evolution of mode C44 as function of caliper modes frequencies (13 or 12 to 14), from -50 to +50% It can be observed that mode C44 can be stabilized, but the instability is in fact transferred to mode C42, so that modifications of caliper mode 13 alone do not yield a better overall behaviour. Altering caliper modes 12 to 14 does not fully stabilize mode C44 but avoids the instability transfer to mode C42, which makes it a more reasonable solution.

Similar results can be obtained by predicting the effect of a variation of caliper mode damping in the assembly. Figure 9 and 10 show the stability diagrams as function of caliper mode 13 damping from 0 to 50% and of caliper modes 12 to 14 from 0 to 50%. Altering caliper mode 13 alone destabilizes mode C42 while full stability can be obtained for a given level of caliper mode damping when working on modes 12 to 14, as shown in figure 10.

Care must be taken when trying to stabilize a given mode as destabilization can occur at other frequencies. Stabilizing a specific mode can let another unstable one become predominant in the response.

Figure 10. Stability diagram for a variation of caliper modes 12 to 14, damping from 0 to 50%

Squeal simulation

The nominal squeal simulation is performed without damping in the frequency band of interest for 100ms and a time step of 10 -6 s. Figure 11 shows the simulation result by plotting the braking torque (which is a synthesis of the friction forces). Such simulation is taken as a limit cycle, although no attempt of mathematical characterization is performed. The fact that an exponential divergence occurs after which a saturation pattern appears in which the signal remains bounded is the criterion retained here.

Initialization strategies can be discussed to obtain such cycle in an optimal computation time [START_REF] Vermot Des Roches | Time frequency analysis of contact-friction instabilities[END_REF].

The limit cycle obtained is however robust to the initialization strategy, provided an initial perturbation is set.

The use of modal sensors allows to observe the contribution of each real mode to the response in terms of mechanical energy. The study presented in [START_REF] Vermot Des Roches | Design oriented time/frequency analysis of contact friction instabilities in application to automotive brake squeal[END_REF] showed that the limit cycles obtained are non linear, so that classical characterizations are difficult. It is in particular interesting to consider new space decompositions, as presented in [START_REF] Vermot Des Roches | Design oriented time/frequency analysis of contact friction instabilities in application to automotive brake squeal[END_REF].

Design oriented frequency/time simulation of industrial brake squeal

Cross study providing design information

Section 3.1 provided design information regarding the brake stability predictions, while section 3.2 briefly illustrated the information that can be output from an actual time simulation. A natural design loop thus appears where rapid predictions can be obtained in the frequency domain with precise component modification indications, which can be validated in the time domain.

The effect of altering caliper modes on the brake assembly has been predicted as efficient to avoid initial instabilities. The use of shape damping in the frequency domain is then used to validate the design direction. Figure 13 indeed shows the residual displacements on the outer pad.

Figure 13. Residual displacement at the caliper/outerpad interface (right). MAC correlation between the original caliper mode and its restitution from the projection on the time reduction basis (left)

The projection remains however relevant as the caliper mode projected on the time reduction basis and restituted on the full model has perfect MAC correlation with the initial caliper mode shape.

The validation is performed by applying shape damping on the caliper modes considered, starting at the squeal simulation end of figure 11. The results are presented in figure 14 for respectively damping caliper mode 13 and modes 12 to 14 from 2 to 10%.

Damping caliper mode 13 only generates a transient evolution from the limit cycle to another one. Increasing the damping ratio yields another divergence due to an instability transfer to mode C42. This was predicted in section 3.1 and is verified in figure 15, which plots the mechanical energy contributions modes R42 and R45. The vibration amplitude is decreased for low damping values, with no stabilization process. For larger damping ratios, the instability is transferred. Adding damping to mechanical assemblies does indeed not necessarily stabilize a system. Such observations were also related in [START_REF] Fritz | Investigation of the relationship between damping and mode coupling patterns in case of brake squeal[END_REF]. Figure 16 however shows that damping caliper modes 12 to 14 has a stabilizing effect.

Low values are not sufficient to obtain an actual stability, however, a threshold exists for which the energy associated to unstable mode C44 is dissipated. Another unstable mode then becomes preponderant, as figure 16 shows that mode C51 starts to express itself in the response. Further computations are needed to evaluate actual vibration levels obtained by mode C51 instability regime. A method such as suggested in [START_REF] Vermot Des Roches | Time frequency analysis of contact-friction instabilities[END_REF] would also allow some predictions.

CONCLUSION

The present paper provided a method to perform brake design regarding vibration levels. In the case of brake squeal, instabilities due to the friction coupling at the pad/disc interface generate complex movements at the level of the brake assembly.

The method presented is based on both frequency and time simulations, through different model reduction strategies. The use of the CMT (Component Mode Tuning) method allows generation of very compact tangent models of the order of a thousand DOF allowing very quick parametric studies, in which 20 design points can be computed in half an hour. In addition, the reduction implements component parameters explicitly at the assembly level, so that local modification can be quickly evaluated.

Time simulation of brake squeal are based on a hybrid model, in which most of the brake is reduced on the trace of its real modes, while the pad/disc interface, critical to the instability, is kept unreduced, thus allowing the full expression of the contact-friction non linearity. The optimization performed on the time integration scheme allows to pass a 100ms simulation in 12h for a system of approximately 30,000 DOF. The computation is then performed for validation purposes once design directions have been obtained in the frequency domain by the CMT. Actual vibration level can indeed be obtained in the time domain.

The application considered highlighted the criticity of structural effects in a brake assembly as the caliper was identified as the most sensitive parameter to a pad/disc generated instability.

Further work to the development presented will focus on evaluating design directions more extensively and suggesting a realistic structural modification to stabilize the brake system in practice.

The importance of structural effects is of enough importance that vibrations in the whole assembly are believed to reach levels at which the linearization of other interfaces is arguable. The study of non linear effects at other component interfaces is thus a clear perspective to the work presented.

Figure 1 .

 1 Figure 1. A general mechanical assembly of component with disjoint interfaces and the stiffness matrix, split into a component part and a coupling part Matrices associated with full or reduced components are non zero for a single component, while interfaces can have non-zero values coupling multiple components. Their diagonal blocks are non null on interface component DOF. In the case of two components one thus has

Figure 3 .

 3 Figure 3. Sample contact (left) and friction (right) laws.

Figure 3

 3 Figure3plots several adapted formulations of penalized contact laws. The contact law retained by Bosch is of the exponential type defined at each contact point by[START_REF] Craig | A review of time domain and frequency domain component synthesis method[END_REF]

Figure 4 .

 4 Figure 4. A modified non linear Newmark scheme

Figure 7 .

 7 Figure 7. Real part sensitivity ranking for unstable modes C44 as function of the component modes (sensitivity computed for a 5kHz frequency variation)

Figure 8 .

 8 Figure 8. Stability evolution of modes C42 and C44 for caliper modes frequency variatons from -50 to +50%.Left. Caliper mode 13, Right modes 12 to 14.

Figure 9 .

 9 Figure 9. Stability diagram for a variation of caliper mode 13 damping from 0 to 50%.

Figure 11 .

 11 Figure 11. Nominal squeal simulation. Braking torque.

Figure 12 .

 12 Figure 12. A free/free caliper mode and its projection on the time model. Color code:Left, strain energy. Right, displacement.

Figure 12

 12 Figure 12 presents a caliper mode and its projection on the time model. Since the time model reduction basis does not exactly feature caliper modes, residual displacements at the caliper interface with other components exists. Figure 13 indeed shows the residual displacements on the outer pad.

Figure 14 .

 14 Figure 14. Braking torque results from the initial squeal simulation continued with caliper modes damping 13 (Top) or 12-14 (Bottom)

Figure 16 .

 16 Figure 16. Mechanical energy contribution of modes R45 and R51 when damping caliper mode shapes 12 to 14 from 2 to 10%

th European Conference on Braking JEF 2010, Lille, France


C Reimann

B Flem

A Arnoldussen

P Englmaier

T E Finne

F Koller

Ø Nordgulen

Gaël Le Roux
email: gael.le.roux.pro@gmail.com

Jeroen E Sonke

Christophe Cloquet

Dominique Aubert

François De Vleeschouwer

Comment on ''The biosphere: A homogeniser of Pb-isotope signals" by

 [START_REF] Reimann | The biosphere: a homegeniser of Pb-isotope signals[END_REF]recently published a study on Pbisotope signature along a 120 km long transect cutting the city of Oslo. Based on concentration but also isotope data, they misinterpret Pb concentration of the biosphere in rural places and explain these large enrichments of Pb as being due to natural processes. The study ignores numerous previous studies either on local, regional or global scales (see reviews by Shotyk and Le Roux, 2005;[START_REF] Callender | Heavy metals in the environment-historical trends[END_REF], and references therein), which clearly demonstrate that anthropogenic Pb emitted in the atmosphere from different sources (leaded gasoline, coal burning, metallurgy, etc.) was and is dispersed worldwide. The study also ignores work on Norway by the Steinnes and colleagues group [START_REF] Harmens | Temporal trends (1990-2000) in the concentration of cadmium, lead and mercury in mosses across Europe[END_REF]Steinnes et al., 2005a, b;Åberg et al., 2004), and measurements and modelling by the EMEP network (www.emep.int/, EMEP, 2005). The study also neglects numerous works on preanthropogenic Pb deposition rate and isotopic signature using continental archives of atmospheric deposition like peat bogs [START_REF] Shotyk | History of atmospheric lead deposition since 12,370 14 C yr BP from a peat bog, Jura mountains, witzerland[END_REF][START_REF] Klaminder | Isotopic trends and background fluxes of atmospheric lead deposition in N Europe: analyses of three ombrotrophic bogs from south Sweden[END_REF][START_REF] Kylander | Refining the pre-industrial atmospheric Pb isotope evolution curve in Europe using an 8000-year old peat core from NW Spain[END_REF][START_REF] Le Roux | Recent atmospheric Pb deposition at a rural site in southern Germany assessed using a peat profile and snowpack, and comparison with other archives[END_REF]. These studies have shown that preanthropogenic Pb atmospheric deposition rate and its Pb isotopic signature is regionally defined, but also that those signals are negligible compared to past 2 ka and recent Pb atmospheric fluxes (Table 1). [START_REF] Reimann | The biosphere: a homegeniser of Pb-isotope signals[END_REF] argue that the Pb isotopic signatures of their biogenic samples ( 206 Pb/ 207 Pb of 1.14-1.17) are possibly explained by mass-dependent and mass-independent stable isotope fractionation processes during biogeochemical cycling, rather than by source mixing of 20th century anthropogenic Pb and local geogenic background Pb. In doing so, they ignore basic knowledge on Pb-isotope geochemistry. In contrast to Zn, Fe or Hg isotopes, 3 Pb-isotopes ( 206 Pb, 207 Pb and 208 Pb) are radiogenic and are endmembers of U-Th decay chains. Therefore the variability of Pb isotopic compositions due to combined U-Th radioactivity and U-Th/Pb elemental fractionation is up to 200 times larger than the variability due to mass-dependent equilibrium and kinetic isotope fractionation [START_REF] Faure | Principles of Isotope Geology[END_REF]. In addition, Reimann et al.'s argument that the recently proposed mass-independent nuclear field shift mechanism [START_REF] Schauble | Role of nuclear volume in driving equilibrium stable isotope fractionation of mercury, thallium, and other very heavy elements[END_REF] could induce large natural variations in the very heavy stable isotopes, such as those of Hg and Pb, is similarly flawed. The nuclear field shift effect (also known as nuclear volume effect) is a theoretical concept invoked to explain experimental U isotope fractionation [START_REF] Biegeleisen | Nuclear size and shape effects in chemical reactions. Isotope chemistry of the heavy elements[END_REF]. [START_REF] Schauble | Role of nuclear volume in driving equilibrium stable isotope fractionation of mercury, thallium, and other very heavy elements[END_REF] theoretical predictions of nuclear field shift fractionation factors for Hg and Tl are <1‰ per atomic mass unit. Although nuclear field shift fractionation has potentially been observed at a very limited level (<1‰) in experimental liquid Hg evaporation [START_REF] Estrade | Mass independent fractionation of Hg isotopes during evaporation and condensation processes[END_REF], it remains to be detected in natural samples. Rather, recent work on Hg isotopes notes the notorious absence of the nuclear field shift effect in mass-independent Hg isotope fractionation in biological samples [START_REF] Bergquist | Mass-dependent and -independent fractionation of Hg isotopes by photo-reduction in aquatic systems[END_REF][START_REF] Laffont | Mercury stable isotope variations in a Bolivian watershed: evidence for mass-independent fractionation[END_REF][START_REF] Epov | Measurement of Hg isotopic ratios on molecular level using GC-MC-ICP-MS hyphenation[END_REF]. Finally, the predicted mass dependent and mass-independent isotope fractionation of Pb and neighbouring heavy elements of <1‰ per amu, is smaller than the long-term analytical reproducibility of Reimann et al.'s sector field ICP-MS Pb-isotope analyses (2‰, 2SD), and therefore undetectable. Moreover, [START_REF] Reimann | The biosphere: a homegeniser of Pb-isotope signals[END_REF] point out that Fe and Zn isotopes are fractionated doi:10.1016/j.apgeochem.2008.05.017 in biological systems but the variation is again within a few ‰ per amu (e.g: [START_REF] Cloquet | Variation in the isotopic composition of zinc in the natural environment and the use of zinc isotopes in bigeosciences: a review[END_REF][START_REF] Guelke | Fractionation of stable iron isotopes in higher plants[END_REF][START_REF] Weiss | Isotopic discrimination of zinc in higher plants[END_REF][START_REF] Viers | Evidence of Zn isotopic fractionation in a soil-plant system of a pristine tropical watershed (Nsimi, Cameroon)[END_REF]. A lower or in the best case a similar variation would be expected for Pb which will change the 206 Pb/ 207 Pb ratio in the third decimal. These observations lead to the same conclusion that such variations could not be detected and are in any case a possible explanation for Pb ''homogenisation".

In summary, stable isotope fractionation cannot, by any physico-chemical means, be responsible for the terrestrial variation in stable Pb-isotope ratios, nor in the variation in Reimann et al.'s samples. Reimann et al.'s conclusion that biological activity fractionates and homogenizes stable Pb-isotope signatures, and that therefore the ''20th century airborne pollution signal as defined by Bindler et al. (1999) ( 206 Pb/ 207 Pb ratio of c.a 1.15) ... is not related to anthropogenic sources." is therefore based on a misconception.

Another hypothesis by Reimann et al. is a preferential leaching of some minerals, which will have another isotope composition than the bulk soil. Lead will be taken up by plants from this soil solution and therefore have a different isotope composition than bulk soil. That is right and this was verified for other isotopic systems (i.e. 86 Sr/ 87 Sr), however this is negligible in comparison with atmospheric deposition of anthropogenic Pb [START_REF] Klaminder | Isotopic trends and background fluxes of atmospheric lead deposition in N Europe: analyses of three ombrotrophic bogs from south Sweden[END_REF].

This was also clearly shown by [START_REF] Bacon | Lead grass in the Scottish uplands: deposition or uptake?[END_REF], when they investigated the isotopic composition of grass in a field enriched in 207 Pb. They showed that in rural places in Scotland, >80% of Pb in grass was derived from atmospheric deposition. In soils, Pb, even in rural places, mainly comes from atmospheric deposition of anthropogenic Pb apart from specific geological settings. Lead soil inventories, which are closely related to atmospheric precipitation, are even used in forested and mountainous areas as indicator of total atmospheric deposition [START_REF] Weathers | The effect of landscape features on deposition to Hunter Mountain, Catskill Mountains[END_REF][START_REF] Weathers | Empirical modeling of atmospheric deposition in mountainous landscapes[END_REF].

Once it is concluded that plants do not ''homogenise" the Pb-isotope signal, the data from [START_REF] Reimann | The biosphere: a homegeniser of Pb-isotope signals[END_REF] can be looked at differently and indeed using the same results, inverse conclusions supported by numerous recent studies (i.e. Lahd [START_REF] Lahd Geagea | Baseline determination of the atmospheric Pb, Sr and Nd isotopic compositions in the Rhine valley, Vosges mountains (France) and the Central Swiss Alps[END_REF][START_REF] Cloquet | Atmospheric pollutant dispersal around an urban area using trace metal concentrations and Pb isotopic compositions in epiphytic lichens[END_REF]) can be adopted. In fact, [START_REF] Reimann | The biosphere: a homegeniser of Pb-isotope signals[END_REF] ignore their own data. The discussion about Fig. 5 in [START_REF] Reimann | The biosphere: a homegeniser of Pb-isotope signals[END_REF] is especially short and only points out that ''completely different processes determine Pb concentrations and isotope ratios in the geosphere and in the biosphere": for sure, because the biosphere and the surface soil layers are influenced by atmospheric deposition of anthropogenic Pb. If only the O-Horizon samples are considered, it appears clearly that there is a relationship between Pb concentration and Pb-isotope composition (Fig. 1). This allows definition of the the isotope composi- tion of an end-member enriched in Pb between 1.14 and 1.16 for 206 Pb/ 207 Pb. The raw results by [START_REF] Reimann | The biosphere: a homegeniser of Pb-isotope signals[END_REF] are also interesting in order to validate or develop bioindicators of atmospheric deposition. Mosses are often used as bioindicators of atmospheric deposition, however in some areas, they can be rarely found. A comparison between the concentrations in mosses and birch leaves (Fig. 2) from the data of [START_REF] Reimann | The biosphere: a homegeniser of Pb-isotope signals[END_REF] show that birch leaves would possibly be good alternatives, even if they represent a shorter time of aerosol interception. The isotopic shift (Fig. 3) of leaves leads to a more radiogenic component compared to mosses. This can be perhaps explained by an additional uptake of Pb from the soil via the tree roots to the leaves (Fig. 4) if it is considered that mosses represent the atmospheric Pb isotopic signature. Further investigations with more precise Pb-isotope measurements and adapted protocol are needed to explain this shift. Other bioindicators such as wood or bark (Fig. 5) may be interesting, especially in order to understand the Pb biogeochemical cycle on a longer time-scale, but then it can be expected that the Pb contamination in the wood is dependent on the history of the tree and its environment. Recently analyses of tree bark using a dedicated sampler and protocol have shown, however, that they can be used to monitor atmospheric pollution and corresponds to 2-8 a of accumulation [START_REF] Lahd Geagea | Baseline determination of the atmospheric Pb, Sr and Nd isotopic compositions in the Rhine valley, Vosges mountains (France) and the Central Swiss Alps[END_REF]. Here (Fig. 5), the data of [START_REF] Reimann | The biosphere: a homegeniser of Pb-isotope signals[END_REF] again show a relationship between concentration and isotopic composition.

To conclude, the set of data from Reimann et al. ( 2008) is inappropriately interpreted and should be revised in order to integrate results and conclusions from previous studies on environmental archives. In rural and remote places from western Europe far from any point source, the isotopic chronology of atmospheric Pb deposition and consequently ecosystem contamination displays a similar history mainly dominated by the dispersion of submicron aerosols carrying Pb from leaded gasoline before its recent banishment.

Consequently to that and to other environmental improvements, Pb emissions in the atmosphere have decreased over 10 a in western Europe along with Pb deposition in rural places. This is causing other previous ''minor" sources of Pb to influence the isotopic composition of Pb in the environment. Will it spatially differentiate this isotopic composition of Pb? This should be verified in the near future using a similar approach to [START_REF] Reimann | The biosphere: a homegeniser of Pb-isotope signals[END_REF], combining isotopic measurement and GIS but without ignoring previous studies specifically on environmental archives.

Åberg, G., Abrahamsen, G., Steinnes, E., Hjelmseth, H., 2004

Fig. 1 .

 1 Fig. 1. Pb isotope composition vs. Pb concentration À1 in O-horizons.

Fig. 5 .

 5 Fig. 5. Pb isotope composition vs. Pb concentration in birch bark.

Table 1

 1 Pb atmospheric accumulation rate (AR) estimated using peat bogs in Europe for ''preanthropogenic" periods (adapted from[START_REF] Kylander | Refining the pre-industrial atmospheric Pb isotope evolution curve in Europe using an 8000-year old peat core from NW Spain[END_REF][START_REF] Le Roux | Recent atmospheric Pb deposition at a rural site in southern Germany assessed using a peat profile and snowpack, and comparison with other archives[END_REF] and recent measured wet deposition

	Study	Country	Pre-anthropogenic AR lg m À2 a À1	206 Pb/ 207 Pb signature	Period
	Le Roux et al. (2005)	Southern Germany	30 ± 14 (1-90, n = 51)	1.18-1.19 (n = 10)	5800-1300 cal. BC
	Shotyk et al. (1998)	Northern Switzerland	10 ± 2 (5-14, n = 18)	1.199-1.21	7000-4100 cal. BC
	Klaminder et al.(2003)	Southern Sweden	1-10	1197-1212	4000-1500 cal. BC
	Kylander et al. (2005)	Galicia, Spain	32 ± 30 (n = 11)	1.25470 ± 0.02575	6100-3200 cal. BC
	Shotyk et al. (2005) and unpublished data	Faroe Islands	0.2-10 (n = 12)		8000-3500 cal. BC
	Study	Country	AR lg m À2 a À1	206 Pb/ 207 Pb signature	Period
	Le Roux et al. (2005)	Southern Germany	2500 ± 1500	1.14-1.16	2003
	EMEP (2005)	Southern Germany	1500		2003
	EMEP (2005)	Southern Norway	800-2300		2003
	EMEP (2005)	Southern Sweden	700-1100		2003

in moss mg kg -1 Pb in Birch leaves mg kg -1

  . Utilization of bark pockets as time capsules of atmospheric-lead pollution in Norway. Atmos. Environ. 36-38, 6231-6237.Fig. 2. Pb concentrations in birch leaves vs. Pb concentrations in mosses.Fig. 3. 206 Pb/ 207 Pb in birch leaves vs. 206 Pb/ 207 Pb in mosses.Fig. 4. ''Isotopic shift ''of birch leaves compared to the atmospheric signal (mosses) vs. 206 Pb/ 207 Pb in birch wood.

		0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1	0	2	4 Pb (8.81 ; 7.2) 8 6	10	206 Pb/ 207 PbBirchLeaves -206 Pb/ 207 PbMosses	-0.005 -0.001 0.003 0.007 0.011 0.015 1.12	1.14	1.16	1.18	1.2	(1.193; 0.033) 1.22 1.24
													206 Pb/ 207 Pb Birch Wood
											1.17	Bark
										Pb	1.16
										206 Pb/ 207	1.14 1.15
											1.13
											1.12
											0		1	2	3	4	5
													1/Pb
		1.175
	Birch leaves	1.165 1.17
	206 Pb/ 207 Pb	1.16
		1.155
		1.15
			1.15	1.155	1.16	1.165	1.17	1.175
					206 Pb/	207 Pbmoss


Maud Leriche
email: maud.leriche@aero.obs-mip.fr

R Lyana Curier

Laurent Deguillaume

Dimitri Caro

Karine Sellegri

Nadine Chaumerliac

Maud Leriche $,£

Numerical quantification of sources and phase partitioning of chemical species in cloud: Application to wintertime anthropogenic air masses at the Puy de Dôme station

Keywords: aerosol particles, cloud chemistry, droplets nucleation, Puy de Dôme station

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

For several years, significant efforts were made by the scientific community to understand tropospheric chemistry. However, these investigations have been mainly focused on homogeneous gas phase chemistry, which is currently well documented [START_REF] Atkinson | Evaluated kinetic and photochemical data for atmospheric chemistry: Volume I -gas phase reactions of O x , HO x , NO x and SO x species[END_REF]2006). At a global scale, cloud chemistry has a potential impact on the atmospheric oxidizing capacity and on climate change via the indirect effect of aerosols on the radiative budget of the earth's system, which is still poorly understood and is a major source of uncertainties in climate forecasting [START_REF] Ramaswamy | Radiative forcing of climate change[END_REF]. Some modeling studies have attempted to assess the overall impact of cloud chemistry on global atmospheric composition and climate with a particular focus on the ozone chemistry [START_REF] Lelieveld | The role of clouds in tropospheric photochemistry[END_REF][START_REF] Liang | Effect of aqueous phase cloud chemistry on tropospheric ozone[END_REF][START_REF] Barth | Effect of marine boundary layer clouds on tropospheric chemistry as analyzed in a regional chemistry transport model[END_REF] and on the sulfur cycle [START_REF] Dentener | Aqueous phase reaction of HNO 4 : The impact on tropospheric chemistry[END_REF][START_REF] Kreidenweis | Modification of aerosol mass and size distribution due to aqueous-phase SO 2 oxidation in clouds: Comparisons of several models[END_REF]. However, the details of in-cloud processes are still poorly understood and some new aspects of cloud chemistry have recently been put forward concerning for instance its potential contribution to the formation of secondary organic aerosol particles [START_REF] Blando | Secondary organic aerosol formation in cloud and fog droplets: a literature evaluation of plausibility[END_REF][START_REF] Gelencsér | Evaluation of the atmospheric significance of multiphase reactions in atmospheric secondary organic aerosol formation[END_REF][START_REF] Carlton | Link between isoprene and secondary organic aerosol (SOA): Pyruvic acid oxidation yields low volatility organic acids in clouds[END_REF]. More efforts are needed to improve our understanding of the cloud impact at local scale critical for radiative and health effects, in particular via integrated modeling and in-situ studies.

A cloud is a multiphase system comprising a gaseous phase, a particulate phase, and an aqueous phase. To assess the role of cloud chemistry at a local scale, interactions between gases, cloud droplets and aerosol particles have to be taken into account in numerical cloud chemistry models. The interactions between aerosol particles and cloud droplets are usually not considered in cloud chemistry models [START_REF] Ervens | CAPRAM2.4 (MODAC mechanism): An extended and condensed tropospheric aqueous phase mechanism and its application[END_REF][START_REF] Williams | The influence of cloud chemistry on HOx and Nox in the moderately polluted marine boundary layer: a 1-D modelling study[END_REF]. Aerosol particles contribute to the chemical composition of cloud droplets by nucleation and impaction scavenging followed by dissolution of soluble compounds. During the cloud lifetime, the cloud chemistry leads to the formation of new chemical species with relatively low volatility such as inorganic and organic acids, which can modify the physico-chemical properties of aerosol particles after the cloud dissipates [START_REF] Feingold | Cloud processing of aerosol as modeled by a large eddy simulation with coupled microphysics and aqueous chemistry[END_REF][START_REF] Yin | Vertical transport of aerosols in a mixed-phase convective cloud and the feedback on cloud development[END_REF] and lead to secondary organic aerosols formation [START_REF] Gelencsér | Evaluation of the atmospheric significance of multiphase reactions in atmospheric secondary organic aerosol formation[END_REF]. For some chemical species, aerosol particles dissolution is the only source in cloud droplets; an example is transition metal ions, and in particular iron, which is well known to play a major role in the oxidizing capacity of clouds [START_REF] Deguillaume | The role of transition metal ions on HO x radicals in clouds: a numerical evaluation of its impact on multiphase chemistry[END_REF]2005). Study of such complex interactions needs laboratory experiments for kinetics and thermodynamics database development and modelling efforts integrating in-situ measurements.

In this framework, the Model of Multiphase Cloud Chemistry model (M2C2) [START_REF] Leriche | Modeling study of strong acids formation and partitioning in a polluted cloud during wintertime[END_REF][START_REF] Deguillaume | The role of transition metal ions on HO x radicals in clouds: a numerical evaluation of its impact on multiphase chemistry[END_REF] has recently been developed to simulate aerosol particles temporal evolution and their activation in cloud droplets. This cloud chemistry model includes aerosol particles, cloud microphysics and chemistry in both aqueous and gas phases. While measurements do not allow the determination of the origin of the different sources of chemical species in the cloud water and to identify these sources, the model provides this information as well as the main chemical production pathways and the impact of cloud microphysical processes on chemistry. In this study, the M2C2 model is used to simulate typical wintertime anthropogenic air masses at the Puy de Dôme station. At this station, continuous measurements (http://wwwobs.univbpclermont.fr/observ/en/index.htm) as well as intensive field campaigns are performed leading to an interesting and explicit multiphase dataset including physico-chemical characterisation of air mass, aerosol particles and clouds encountered at the high altitude site. The present study aims at evaluating the role of the chemistry versus both the gas scavenging and the aerosol particles nucleation scavenging in the cloud water chemical composition. Then, the role of the chemistry in the phase partitioning can be determined and the various sources of chemical species inside the cloud, either from gas phase or from particulate phase are quantified.

Firstly, the model, its new developments and the initial conditions of the simulation are presented. Secondly, numerical concentrations in cloud water are compared to experimental values to assess the quality of the simulation. Then, the analysis of the different sources of chemical species in cloud water is performed. Finally, the partitioning of chemical species obtained from measurements and from simulation are compared and discussed.

Model description

2-1 Model formulation

The Model of Multiphase Cloud Chemistry (M2C2) results from the coupling [START_REF] Leriche | Coupling quasi-spectral microphysics with multiphase chemistry: A case study of a polluted air mass at the top of the Puy de Dôme mountain (France)[END_REF] between a multiphase chemistry model described in [START_REF] Leriche | A model for tropospheric multiphase chemistry: Application to one cloudy event during the CIME experiment[END_REF] and a two-moment warm microphysical scheme module predicting the number concentration of cloud droplets and raindrops and the mixing ratio of cloud water and rainwater categories using log-normal distributions [START_REF] Caro | A two-moment parameterization of aerosol nucleation and impaction scavenging for a warm cloud microphysics: description and results from a two-dimensional simulation[END_REF]. The dynamical framework of the model is an air parcel [START_REF] Gérémy | A modelling study on the activation of small Aitken-mode aerosol particles during CIME 97[END_REF]. The chemistry included in the chemical module is explicit. The exchange of chemical species between the gas phase and the aqueous phase is parameterized following the mass transfer kinetic formulation developed by Schwartz (1986). The aqueous phase chemistry includes the detailed chemistry of HO x , chlorine, carbonates, NO y , sulphur, the oxidation of organic volatile compounds (VOCs) with one carbon atom [START_REF] Leriche | Modeling study of strong acids formation and partitioning in a polluted cloud during wintertime[END_REF], the chemistry of transition metal ions for iron, manganese and copper [START_REF] Deguillaume | The role of transition metal ions on HO x radicals in clouds: a numerical evaluation of its impact on multiphase chemistry[END_REF] and the oxalate oxidation and interaction with iron chemistry (Table 1). Oxalate represents a strong organic iron complexing agent and the highest fraction of dicarboxylic acids in Table 1 aerosol particles [START_REF] Legrand | Atmospheric year-round records of dicarboxylic acids and sulfate at three French sites located between 630 and 4360 m elevation[END_REF]. The value of the accommodation coefficient for oxalic acid is estimated to be 0.05 [START_REF] Leriche | A model for tropospheric multiphase chemistry: Application to one cloudy event during the CIME experiment[END_REF] because no measurement is available. The model considers calculations of photolysis frequencies both in the gas phase and in cloud droplets. The pH is calculated at each time step by solving the electroneutrality equation.

The activation of aerosol particles into cloud droplets is parameterized according to [START_REF] Abdul-Razzak | A parameterization of aerosol activation, 1, Single aerosol type[END_REF] and [START_REF] Abdul-Razzak | A parameterization of aerosol activation, 2, Multiple aerosol types[END_REF]. Aerosol particles are represented by the sum of four log-normal distributions as follows: Where N ap is the total number concentrations of aerosol particles, N api is the number concentrations of aerosol particles for each mode; a mi is the median radius of the lognormal distribution for each mode; i is the geometric standard deviation of the lognormal distribution of each mode and a ap is the radius of aerosol particles.

The parameterization used to allow for nucleation of cloud droplets comes down to find an expression for the maximum of the supersaturation on each model time step [START_REF] Abdul-Razzak | A parameterization of aerosol activation, 1, Single aerosol type[END_REF][START_REF] Abdul-Razzak | A parameterization of aerosol activation, 2, Multiple aerosol types[END_REF] considering the physico-chemical properties of aerosol particles and based upon the Köhler theory [START_REF] Köhler | The nucleus in the growth of hygroscopic droplets[END_REF][START_REF] Pruppacher | Microphysics of Clouds and Precipitation (2 nd Edition)[END_REF]. For each mode, the activated fraction is calculated and leads to the number of droplets newly formed by nucleation and to the associated cloud water mixing ratio.

The introduction of the cloud droplets nucleation process in the model allows for the monitoring of the evolution of aerosol particle size distribution with time. The microphysical scheme for aerosol particles is a two-moment scheme in which the prognostic variables are the number concentration of aerosol particles and the corresponding volumetric concentration [START_REF] Caro | A two-moment parameterization of aerosol nucleation and impaction scavenging for a warm cloud microphysics: description and results from a two-dimensional simulation[END_REF]. At each time step, the evolution of aerosol particles size distribution is computed from the number of nucleated droplets and the associated activated aerosol mass. Until now, the impaction scavenging is neglected in the model, which is not very restrictive as nucleation scavenging in non precipitating clouds represents more than 90% of the total aerosol particles scavenged by the cloud droplets [START_REF] Flossmann | The scavenging of two different types of marine aerosol particles using a two-dimensional detailed cloud model[END_REF].

The impact of aerosol particles on cloud chemistry is considered in the model via the activated mass of aerosol particles. The fraction of soluble chemical species in aerosol particles is used to compute the chemical concentrations in cloud droplets of species coming from the activated particles. The mass fraction of soluble chemical species in each aerosol mode is prescribed from measurements.

Initialization of the model: Anthropogenic air mass at the Puy de Dôme station

The Puy de Dôme station is located in the centre of France at the summit of the Puy de Dôme Mountain (1465 m above sea level) and is well recognized in the scientific community [START_REF] Van Dingenen | A European aerosol phenomenology -1: physical characteristics of particulate matter at kerbside, urban, rural and background sites in Europe[END_REF][START_REF] Putaud | A European aerosol phenomenology -2: chemical characteristics of particulate matter at kerbside, urban, rural and background sites in Europe[END_REF] In order to perform a simulation of a typical wintertime anthropogenic air mass at the Puy de Dôme station with the M2C2 model, a scenario based upon the air mass classification of Sellegri et al. (2003a), which defines a typical anthropogenic influenced free tropospheric air mass encountered at the station, was elaborated using the available measurements.

The air trajectory used to initialize the parcel comes from a three dimensional simulation of a typical dynamical situation corresponding to wintertime polluted air masses at the Puy de Dôme station coming from the North-North-East. This three dimensional simulation was performed with the meso-scale non-hydrostatic Clark model [START_REF] Clark | Multi-domain simulations of the timed dependent Navier-Stokes equations: benchmark error analysis of some nesting procedure[END_REF] using three nested grids. The trajectory used in M2C2 corresponds to the back-trajectory computed by the three dimensional model arriving at noon local time at the summit of the Puy de Dôme. The trajectory is shown in Figure 1 together with the corresponding vertical wind speed. This trajectory at noon corresponds to the maximum of the photolysis leading to the highest intensive photochemistry. The beginning of the lifting of the air parcel is also indicated in Figure 1. At this moment, the M2C2 simulation run was started for a duration of twelve minutes.

Figure 1 The gas phase chemical species are initialized using measurements from [START_REF] Voisin | Scavenging of acidic gases (HCOOH, CH 3 COOH, HNO 3 , HCl, and SO 2 ) and ammonia in mixed liquid-solid water clouds at the Puy de Dôme mountain (France)[END_REF], performed during the CIME experiment in winter 1997-1998. During this campaign, a large amount of chemical species was measured kin addition to continuous measurements of, e.g., ammonia and nitric acid. An event corresponding to an anthropogenic air mass arriving at the Puy de Dôme station as defined in Sellegri et al.

(2003a) was selected in this database. Initial concentrations taken from this event for gas phase chemical species are shown in Table 2.

The initialization of the physico-chemical properties of aerosol particles are derived from Sellegri et al. (2003a) and [START_REF] Sellegri | Etude du processus d'activation des gouttelettes de nuage: Implications en chimie multiphases[END_REF] for anthropogenic influenced free Figure 2 Figure 3 particles in each mode was derived from the aerosol particle measurements database following [START_REF] Sellegri | Etude du processus d'activation des gouttelettes de nuage: Implications en chimie multiphases[END_REF]. The mass fraction of iron is estimated from Particle Induced X-Ray Emissions (PIXE) data (Sellegri et al., 2003a) and assuming an iron solubility of 15% [START_REF] Deguillaume | Transition metals in atmospheric liquid phases: sources, reactivity and sensitive parameters[END_REF]. The oxalic acid is assumed to represent 60% of total measured organic acids and formic acid, 40%. This partitioning followed the ratio observed between dicarboxylic and mono-carboxylic acids (Sellegri et al., 2003a) if we assume that mono-carboxylic acids are mainly represented by formic acid and dicarboxylic acids by oxalic acid.

Results and discussion

Comparisons between simulated and observed cloud concentrations for main chemical species

The simulated concentrations and pH in cloud water are compared with observed values as shown in Table 4. Numerical results are averaged values on the whole simulation.

Observed values are taken from previous studies using the database of chemical measurements available in cloud water at the Puy de Dôme station coming from the ORE-BEAM project mentioned above [START_REF] Marinoni | Cloud chemistry at the Puy de Dôme: variablility and relationships with environmental factors[END_REF][START_REF] Parazols | Speciation and role of iron in cloud droplets at the puy de Dôme station[END_REF], which used air mass classification similar to Sellegri et al. (2003a). The maximum value of the simulated cloud liquid water content is 0.4 g m -3 at the top of the mountain and the mean value is 0.16 g m -3 . For anthropogenic influenced free tropospheric air masses, the measured cloud liquid water content is around 0.3 g m -3 and hence there is a reasonable agreement between model and measurements. Simulation and in-situ measurement indicate both acidic pH values for the wintertime anthropogenic air masses at the Puy de Dôme station. In general, results from the model compare favorably with the measured concentrations. However, formate and oxalate are overestimated by the model and iron is underestimated.

Table 4 For both carboxylic acids (formate and oxalate), the overestimation is about a factor 2.

However, the amplitude of this factor is quite uncertain because of the low accuracy for carboxylic acids concentrations determined by ion chromatography due to their low concentrations close to the detection limit. Overestimation of formic acid has also been obtained by the SPACCIM model applied on a continental event of the FEBUKO field campaign, which took place at the Schmücke Mountain (937m above sea level) located in central Germany in October 2001 [START_REF] Wolke | SPACCIM: A parcel model with detailed microphysics and complex multiphase chemistry[END_REF]. However, the SPACCIM chemical mechanism takes into account organic acids with up to four carbon atoms and the modelled acetic acid concentrations compared well with those obtained from the FEBUKO measurements. This suggests that the observed overestimation of formic acid is not solely due to weighing formic acid as the source of mono-carboxylic acid.

Measurements show that oxalic acid is the main dicarboxylic acid in aerosol particles.

Moreover, recent studies have shown that, in addition to aerosol particles, chemical reactivity is an important source of oxalic acid in cloud droplets, especially in rural areas influenced by biogenic emissions [START_REF] Ervens | A modeling study of aqueous production of dicarboxylic acids: 1. Chemical pathways and speciated organic mass production[END_REF][START_REF] Sorooshian | Oxalic acid in clear and cloudy atmospheres: Analysis of data from International Consortium for Atmospheric Research on Transport and Transformation 2004[END_REF].

Thus, no satisfactory reason is available to explain the simulated overestimation of oxalic acid by the model.

The only source of iron in cloud droplets is from aerosol particles, hence, its underestimation is explained by: the uncertainty on measured iron concentration in aerosol particles which are often close to the detection limit; the fact that measurements of elemental species in aerosol particles were not available for all anthropogenic events used to estimate physico-chemical properties of initial aerosol particles; and the two assumptions made to estimate the iron soluble part in the particulate phase (Table 3).

Measurements of iron in particles have been made using PIXE analysis which gives the elemental species concentration, the two assumptions made to estimate soluble iron mass fraction are: all iron in the particulate phase is iron oxide (FeO) and the solubility of iron is 15%. The first assumption could lead to an underestimation of iron mass fraction in the particulate phase depending on the real chemical species of iron which occur in various forms: iron oxides, iron hydroxides, adsorbed iron, etc… [START_REF] Deguillaume | Transition metals in atmospheric liquid phases: sources, reactivity and sensitive parameters[END_REF]. For the second assumption, available values of iron solubility are in the range of 0.01 to 35% depending on the occurrence of iron in a crystalline or in a carboneous matrix [START_REF] Desboeufs | Dissolution and solubility of trace metals from natural and anthropogenic aerosol particulate matter[END_REF] and the dissolution is instantaneous in the model although it is well known that the dissolution is a kinetic process [START_REF] Deguillaume | Transition metals in atmospheric liquid phases: sources, reactivity and sensitive parameters[END_REF]. Thus, the underestimation of iron by the model is still an open question.

In the model, the cloud water chemical composition is initialized by chemical species coming from the activation of aerosol particles. In this study, the whole coarse mode and accumulation mode 2 are activated, as well as a large fraction of the accumulation mode 1, and only a small part of the Aitken mode. Cloud water concentrations do not only reflect the chemical composition of activated aerosol particles which would lead to similar concentrations for nitrate and sulphate as shown by the chemical composition of aerosol particles (Figure 2). Indeed, the nitrate concentration in cloud droplets is higher than the concentration of sulphate due to the impact of chemical reactivity on cloud water concentrations as will be demonstrated later. To quantify the contribution of reactivity, the sources of chemical species in cloud water have to be determined.

Quantifications of chemical species sources in cloud water

While measurements do not allow determining the origin of the different chemical species in cloud water without artefacts, the model can provide this information.

Analysis of the sources of chemical species (particulate, gaseous, or chemical reactivity in cloud water) is shown in Table 5 as relative contributions to their final concentration in the liquid phase. These results show that, for inorganic compounds, the scavenging of Table 5 aerosol particles predominates due to the complete activation of the accumulation mode 2 and the coarse mode with the largest contribution for sulphate and ammonium.

In the cloud droplets, 72.5% of the ammonium concentration comes from the dissolution of activated aerosol particles. The 27.5% remaining comes from the gas phase, as ammonia is a very soluble gas, and it is nearly completely scavenged by the cloud water. Similar ratios are observed for sulphate, with 73% of its cloud droplets concentration coming from the dissolution of activated aerosol particles whereas the remaining 27% are chemically produced in the aqueous phase. This aqueous phase production is due to the scavenging of gaseous sulphur dioxide in cloud droplets where it is oxidized into sulphate by pernitric acid and hydrogen peroxide. Despite the high iron concentrations in cloud droplets, as [START_REF] Deguillaume | The role of transition metal ions on HO x radicals in clouds: a numerical evaluation of its impact on multiphase chemistry[END_REF] showed, the catalysis of sulphur autoxidation by iron is negligible in the formation of sulphate.

For nitrate, the situation is quite complex: 56% of the nitrate-containing cloud droplets come from activated aerosol particles, 15.5% come from the gas phase initial concentration, 23% come from a gas phase in cloud production and the remaining 5.5%

come from an aqueous phase in cloud production. As nitric acid is a very soluble gas, its initial gas phase concentration and its gas phase production are completely and rapidly scavenged in cloud droplets in a similar way to ammonia. The simulated air mass has anthropogenic characteristics, hence the nitrogen oxides concentration levels are high (Table 2) and lead to a significant production of nitric acid in the gas phase and of nitrate in the aqueous phase. The main chemical pathway producing nitrate in the aqueous phase is the oxidation of sulphur dioxide by pernitric acid, which also contributes to the sulphate production [START_REF] Leriche | Modeling study of strong acids formation and partitioning in a polluted cloud during wintertime[END_REF]. This analysis of nitrate sources in cloud droplets explains the higher concentration of nitrate than sulphate in cloud droplets mentioned above (Table 4). Indeed, as the nitrate mass fraction in activated particles is as high as those of sulphate (Table 3) and its main aqueous phase chemical production pathway is the same as for sulphate, its higher concentration in cloud droplets may be explained by its efficient gas phase chemical production followed by its scavenging. High NO x levels lead to this efficient nitrate production, which confirms the hypothesis of Sellegri et al. (2003b) of a favoured nitrate production in cloud over sulphate under the atmospheric conditions typical of the Puy de Dôme site.

In cloud droplets, iron, which is a major OH precursor [START_REF] Deguillaume | Transition metals in atmospheric liquid phases: sources, reactivity and sensitive parameters[END_REF], come only from the dissolution of activated particles as seen before, as well as the oxalic acid which can inhibit the role of iron in the OH production due to its complexing ability.

Introducing aerosol particles in the M2C2 model allows for correctly representing the in-cloud source of these two chemical species, and thus, the in-cloud OH production.

Analysis of the relative contribution of the main OH production pathways in cloud droplets shows that iron reactivity contributes 30% to the OH production in droplets via the photolysis of iron(III) and the Fenton reaction (Table 6) while its scavenging is the main pathway contributing 57%. In order to study the impact of oxalate on OH production in cloud droplets, an additional simulation was performed neglecting oxalate in aerosol particles. Results show that, when oxalate chemistry is not considered, the contribution of iron in OH aqueous phase production is higher, 57%, and the scavenging is less important (Table 6). Oxalate chemistry leads to a decrease in the contribution of iron in OH production as expected by its inhibiting role due to its high complexing power on iron. However, these results show that, even if oxalate chemistry is considered and despite a relatively low photochemical activity due to the hibernal situation, iron contributes significantly to the OH radical production inside the cloud.

All these results show the interest of such an approach, which allows for quantifying the relative importance of gas and aqueous phase chemistry against aerosol particles as Table 6 source of soluble chemical species in cloud droplets by considering explicitly the chemical source coming from activated particles.

Partitioning of chemical species among cloud phases: droplets, interstitial aerosol particles and gas

To get a global view of chemical species in cloud, the phase partitioning between the gas phase, the particulate phase and the liquid phases is calculated for nitric acid, sulphuric acid, ammonium and formic acid. Results of this analysis are shown in Figure 4, which can be directly compared to the estimated phase partitioning from measurements shown on Figure 5 (adapted from Sellegri et al., 2003b) for anthropogenic influenced free tropospheric air masses sampled during the 2001 Puy de Dôme campaign in winter 2001. The experimental procedure is explained in detail in Sellegri et al. (2003b), but we will briefly summarize its principle. During a cloud event, interstitial and residual phases are separately sampled with specific inlets: Round Jet Impactor (RJI) and Counterflow Virtual Impactor (CVI, [START_REF] Noone | Design and calibration of a counterflow virtual impactor for sampling of atmospheric fog and cloud droplets[END_REF]. The gas phase is sampled using mist chambers, while the particulate phase is sampled using cascade impactors, both in the residual and interstitial phases. Moreover, the bulk liquid concentrations were measured with a CDI device (Cloud Droplet Impactor, [START_REF] Kruisz | A cloud water samples for high wind speeds[END_REF]. While the contributions of gaseous, particulate and liquid concentrations to the total atmospheric concentration of a given species can be experimentally evaluated within the experimental uncertainties, it is difficult to evaluate the particulate and gaseous origin of the liquid concentration by an experimental procedure. Sellegri et al. (2003b) attempted to do this by neglecting the in-cloud chemical liquid production. The experimental gas and particulate originating liquid concentrations were deduced by combining the results from the CVI device and the results from the CDI device. The strong hypothesis in this calculation was that the fraction of gaseous species which were Figure 4 Figure 5 evaporating from cloud droplets was equivalent to the one which had absorbed from the original gas phase and dissolved in the liquid phase during the formation of the cloud.

In the present study, it is possible to evaluate this hypothesis by modelling the fractions of gas-originating and particle-originating compounds present in the liquid phase as well as the contribution of chemical reactivity in liquid phase concentrations.

First, the model can be compared to measurements regarding the three phase apportionment (interstitial gas, interstitial particles and liquid phase). The main difference between measurements and simulation lies in the fraction of the interstitial particulate phase. This fraction is close to zero for numerical results (in black in Figure 4) whatever the chemical species considered whereas it ranges from 1% to 12.4% of the total concentration for in-situ data. Measurements showed bulk mass scavenging efficiencies for anthropogenic air mass of 0.91 for sulphate, 0.96 for nitrate and 0.88 for ammonium, no data were available for formic acid (Sellegri et al., 2003c) whereas numerical results show bulk mass scavenging efficiencies close to 1 for all species as almost all aerosol mass is activated. However, the in-situ number scavenging efficiency is between 0.4 and 0.48 which is similar to the simulated value (0.41) highlighting the ability of the model to simulate the cloud droplets nucleation process. The discrepancy in the mass scavenging efficiencies indicates an overestimation of the activated mass by the model. This drawback could be improved in the future by the implementation of the effect of organic compounds on the surface tension and on the bulk hygroscopicity of the particles following Abdul-Razzak and Ghan (2004).

With respect to measurements, in addition to gas and aerosol particles scavenging, the simulation allows to identify the chemical production sources in the liquid phase as detailed above. For an improved comparison with the partitioning obtained from measurements, the contribution of the chemical gas phase production in the liquid phase concentrations has been included in gas-originating liquid. For nitric acid, the fraction present in the gas phase agrees between simulation and measurements. For the liquid phase, the simulation shows a small but appreciable amount of nitrate originating from chemical production. The amount of gas and particulate originating nitrate in the liquid phase as determined from the model simulation and the measurements are within experimental uncertainties. For sulphuric acid, a direct comparison between measured and simulated results is not possible because the phase partitioning obtained from measurements is for the system SO 2 /SO 4 2-

. The experimental findings that a large amount of sulphate in cloud droplets originates from particles are confirmed by the simulation. For ammonium, measured and simulated partitioning show completely different patterns: while the measured interstitial gas fraction is about 56 %, the simulated fraction is 1% of the total atmospheric ammonium concentration. Moreover, for measurements, the amount of gas-originating ammonium in the liquid phase is close to the particles-originating concentration, while it is three times more for the simulation.

These two discrepancies could be related to high ammonium emission around the Puy de Dôme Mountain not considered in the simulation and leading to observed subsaturation of the liquid phase as underlined by Sellegri et al. (2003b). For formic acid, the simulated and measured partitioning compare well with the major contribution coming from the gas phase. Moreover, the simulation shows a negligible production in the liquid phase, which is consistent with the observations.

Finally, such a modelling approach leads to a global view of the partitioning of chemical species in clouds including the identification of their various sources. This approach is globally in quite good agreement with the measured partitioning, which makes it a valuable tool for evaluating liquid concentrations in clouds, without developing a heavy instrumental set-up.

Conclusion

Partitioning of chemical species between the liquid, gaseous and particulate phases was simulated for typical wintertime anthropogenic air masses at the Puy de Dôme research station using the M2C2 model. The M2C2 model has been extended to include the aerosol particles evolution due to nucleation scavenging which is parameterized using the Abdul-Razzak and Ghan ( 2000) formulation. The model simulates the chemical and the microphysical evolution of clouds including the interactions between aerosol particles and chemical reactivity in both liquid and gaseous phases. Modelled chemical concentrations could be compared to measurements performed at the station during extensive field campaigns.

The comparison of the simulated and measured concentrations in cloud water shows a good agreement for major species found in droplets, i.e. for sulphate, nitrate and ammonium ions. The model and measurements agree to a smaller extent for organic acids which are overestimated and iron which is underestimated by the model. The reasons for the overestimation of organic acids are not clear and the discrepancy could be partly due to larger uncertainties in the experimental determination of organic concentrations in cloud water, compared to major ions. The underestimation of iron can be explained by the uncertainties on its measured concentration in the particulate phase and by the assumptions made for estimating its soluble fraction in aerosol particles depending on its chemical form into the particle matrix and on its solubility.

The simulated concentrations in cloud water correctly reflect the chemical composition of the activated aerosol particles. Indeed, the analysis of the sources of the chemical species in cloud water shows an important contribution from nucleation scavenging of particles which prevails for nitrate, sulphate and ammonium. However, the gas scavenging is the main source for formic acid in cloud water. The chemical reactivity is a significant source for nitrate and for sulphate. In addition, the calculation of OH production in cloud droplets shows the significant contribution from iron even if the oxalate inhibiting role is considered.

Finally, the model could further estimate which fraction of the liquid droplets was due to particulate scavenging, gas dissolution, or chemical reaction in the liquid phase.

Again, results could be compared to experimental measurements performed by Sellegri et al. (2003b). Numerical results show an underestimation of the interstitial particulate phase which could be due to an overestimation of the activated mass by the model.

However, the numerical number scavenging efficiency of particles is in agreement with measured value. Another important discrepancy is found for ammonium partitioning:

while measurements show a predominant fraction of the species present in the interstitial gas, the simulation prevails almost all of ammonium in the liquid phase. An explanation for this discrepancy may be the large emissions of ammonium around the Puy de Dôme as suggested by measurements. In general, other results show good agreement with the partitioning obtained from measurements. Moreover, the simulation provides the quantification of chemical in-cloud production which can not be obtained directly from measurements.

In conclusion, this study shows the ability of the M2C2 model to study the complex interactions between cloud microphysics, aerosol particles and photo-chemical reactivity. In future work, such an approach will be extended to the study of organic compounds chemistry in cloud including photochemistry and biologic pathways using

  tropospheric aerosols without marine input sampled at the Puy de Dôme station. The chemical composition for these aerosol particles is shown in Figure2for each four identified modes: Aitken, accumulation 1, accumulation 2 and coarse mode. Usually, at the Puy de Dôme station, the accumulation particles are spread in two distinct modes with different chemical signature as shown in Figure2. As expected, the dust contribution is significant for the coarse mode whereas the organic contribution decreases with the particle size. The coarse mode appears to be acidic with a low contribution of ammonium. The initial number distribution of aerosol particles corresponding to this chemical composition are represented in Figure3together with the four log-normal distributions used to fit it. Starting from the left-hand side, the first fitted distribution represents the Aitken mode, the second and the third ones the two accumulation modes and the fourth one the coarse mode. Table3presents the physicochemical properties of aerosol particles and the fitting parameters of the number distribution used in the simulation. The density, the number of dissociated ions, the mass fraction of soluble materials, the mass fraction related to soluble mass of nitrate, sulphate, ammonium, formate, oxalate and soluble iron, and the molar mass of aerosol Table2

  Figure captions

Figure 2 .

 2 Figure 2. Chemical composition of the four modes of aerosol particles for anthropogenic air mass sampled at the top of Puy de Dôme mountain, adapted from Sellegri et al. (2003a).

Figure 3 .

 3 Figure 3. Fit of the observed aerosol mass distributions with four lognormal distributions for anthropogenic air mass sampled at the top of Puy de Dôme mountain, adapted fromSellegri et al. (2003a). The solid grey line is the distribution obtained from measurements (obtained by weighting), the dotted grey line is the fitted distribution, the solid black line is the first lognormal distribution for the Aitken mode, the dot-dashed black line is the second lognormal distribution for the accumulation mode 1, the dotted black line is the third lognormal distribution for the accumulation mode 2, and the dashed black line is the fourth lognormal distribution for the coarse mode.

Figure 4 .

 4 Figure 4. Simulated phase partitioning between the gas, particulate and liquid phases in cloud. The contributions of activated particles, gas scavenging and chemical production are reported for the liquid phase.

Figure 5 .

 5 Figure 5. Estimated phase partitioning from measurements for anthropogenic influence free tropospheric air masses between the gas, particulate and liquid phases in cloud

  Intensive field campaigns are regularly organized at the station like the European CIME experiment (Wobrock et al., 2001) as well as various national campaigns as the 2001 Puy de Dôme campaign during winter 2001.

	Atmosphérique et Modifications anthropiques) project database (http://wwwobs.univ-
	bpclermont.fr/atmos/orebeam/orebeampage1.htm) including measurements of iron (total
	content and oxidation state), ionic chromatography, dissolved organic carbon;
	conductivity, redox potential, pH and some biological parameters (see Amato et al.,
	2005).
	. The station is equipped
	with continuous measurements since June 1995 for meteorological variables and ozone;
	since February 1998 for radiation fluxes (UV, visible, diffusive fluxes and J(NO 2 ));
	since May 2000 for number concentration of aerosol particles; since September 2000 for
	black carbon content; since November 2000 for CO 2 ; since April 2002 for CO; since
	January 2003 for SO 2 and NO x ; since April 2003 for NO y ; and since June 2003 for drop
	surface area, liquid water content and effective radius when a cloud is present on the site
	(free access to data from http://wwwobs.univ-bpclermont.fr/observ/en/index.htm). In
	addition, detailed measurements of the cloud liquid phase are made since February 2001
	for selected cloud events mainly during winter and spring in the framework of the ORE-
	BEAM (Observatoire de Recherche en Environnement Biophysicochimie de l'Eau

Table 1 .

 1 Henry's law constant and equilibrium constants for oxalic acid, and oxalate oxidation and iron-oxalate chemistry in M2C2 model. No temperature dependence (i.e. activation energy) is available for these equilibriums and reactions.

			Henry's law constant	H (M atm -1 )	References
			COOHCOOH(g)	5.0 10 8	Saxena and Hilderman, 1996
			Equilibrium constant	K (M)	Reference
			H 2 C 2 O 4	H + + HC 2 O 4	-	5.6 10 -2	Smith and Martell, 1976
			HC 2 O 4	-	H + + C 2 O 4 2-	5.42 10-5	Smith and Martell, 1976
					Reactions	k 298 (M -n+1 s -1 )	References
	HC 2 O 4	-+ OH C 2 O 4	-+ H 2 O	(A1)	4.7 10 7	Getoff et al., 1971
	C 2 O 4	2-+ OH C 2 O 4	-+ OH -	(A2)	7.7 10 6	Getoff et al., 1971

Table 2 .

 2 Initial gas phase concentrations from[START_REF] Voisin | Scavenging of acidic gases (HCOOH, CH 3 COOH, HNO 3 , HCl, and SO 2 ) and ammonia in mixed liquid-solid water clouds at the Puy de Dôme mountain (France)[END_REF]

	Species	Gaseous concentration (ppbv)
	N 2	7,8.10 8
	O 2	2,1.10 8
	H 2 O	5,7.10 7
	O 3	37
	NO	0.9
	NO 2	6.2
	CH 4	1600
	CO	100
	HNO 3	0.15
	HNO 2	0.054
	H 2 O 2	0.07
	HCHO	6
	CH 3 OH	2.25
	SO 2	0.45
	HCOOH	0.28
	NH 3	0.13
	HCl	0.07

Table 3 .

 3 Initial physico-chemical properties of aerosol particles; the first three parameters are defined in equation (1), api is the density of aerosols for mode i, i is the number of dissociated ion for mode i, i is the mass fraction of soluble materials in mode i, M api is the molar mass of mode i. The mass fraction of chemical species is related to soluble mass.

	Mode	N api (# cm -3 )	a mi (µm)	log i	api (g cm -3 )	i	i (%)	M api (g mole -1 )
	1	111.9	0.038	0.255	1.75	2	31	175
	2	4.2	0.205	0.278	1.74	2	40	274
	3	1.5	0.330	0.041	1.74	2	61.5	248
	4	0.026	1.300	0.301	1.75	2	65	272
	Mode	Nitrate (%) Ammonium (%) Sulphate (%) Formate (%) Oxalate (%)	Iron (%)
	1	30	30.8	26.4	2.2		3.3	0
	2	29.6	18.2	35	1.2		1.8	0.1
	3	30	14.8	39.5	0.8		1.2	0.2
	4	43	2.7	19	2.5		3.7	0.8

Table 4 .

 4 Comparison of observed and simulated averaged pH value and concentrations in cloud water in µg L -1 , except for iron in µM. Observed values are adapted from[START_REF] Marinoni | Cloud chemistry at the Puy de Dôme: variablility and relationships with environmental factors[END_REF] from Parazols et al. (2006) for anthropogenic air mass sampled at the Puy de Dôme station.

		pH	Nitrate	Ammonium	Sulphate	Formate	Oxalate	Iron (µM)
	Observed 3.7 ± 0.4 10702 ± 6631 3074 ± 934 5460 ± 1584 412 ± 304 155 ± 114 5.3 ± 3.2
	Simulated	3.6	14880	2520	6940	920	360	1.9

Table 5 .

 5 Contribution to the concentrations of nitrate, sulphate, ammonium and formic acid in the liquid cloud droplet at the end of the simulation from activation of aerosol particles, from scavenging of gases and from chemical reactivity .

		Aerosols scavenging	Gas scavenging	Chemistry Aqueous phase production Gas phase production
	Nitrate	56	15.5	23	5.5
	Ammonium	72.5	27.5	0	0
	Sulphate	73	0	0	27
	Formate	20	75	0	5

Table 6 .

 6 Contribution of gas scavenging, iron reactivity and other pathways in OH cloud droplets production for the reference simulation including iron-oxalate chemistry and a simulation neglecting oxalate chemistry.

		OH scavenging	Iron reactivity Iron(III) + h Fenton	Other pathways
	Iron and oxalate	57%	26%	4%	13%
	Iron	36%	55%	2%	7%

line figure Click here to download line figure: Figure1.ppt

				Aitken mode		accumulation mode 1
	NH 4	other ions 5%	dust 0%	H 2 O 3%	BC 21%	7% other ions 2% NH 4	dust 1%	H 2 O 3%	BC 11%
	9%
	SO 4 8%					SO 4 14%
	NO 3
	9%					NO 3
	OA					12%				OM-OA
	2%					OM-OA 43%	OA 1%			49%
			accumulation mode 2				coarse mode
	dust other ions	H 2 O	BC				H 2 O	BC
	NH 4	2%	2%	3%	10%		dust		3%	2%	OM-OA
	9%						14%			16%
										OA
						OM-OA				4%
						29%
						other ions
	SO 4					17%
	25%
						NH 4
						OA	2%			NO 3
				NO 3	1%	SO 4			29%
				19%		13%

line figure Click here to download line figure: Figure5.ppt

				HCOOH/HCOO -		HCl/Cl -			HNO 3 /NO 3	-
	88.6%			76.8%
										39.5%
		11.3%	11.0%	22.2%	21.3%	6.9%	90.4%	82.1%	3.4% 1.0%	95.6%
	0.1%					2.6%			56.1%
			0.4%	1.0%		0.9%				8.3%
	SO 2 /SO 4	2-		H 2 C 2 O 4 /C 2 O 4	2-		NH 3 /NH 4	+
										Interstitial gas
	3.2%		30.4%	23.3%		33.7%	55.6%			15.9%	Interstitial particles
	7.5%	89.3%	59.0%	64.3%	30.6%		40.6%	24.6%	Liquid phase
				12.4%						gas-originating Liq.
										part-originating Liq.

Acknowledgments

The authors are very grateful to Wolfram Wobrock for providing the three dimensional simulation with the Clark model. This work was supported by the "Programme National de Chimie Atmosphérique" (PNCA) of the INSU (Institut des Sciences de l'Univers) and by the ORE-BEAM ("Observatoire de Recherche en Environnement Biophysicochimie de l'Eau Atmosphérique et Modifications anthropiques") project of the French Ministry for Research and Education. Computer resources were provided by I.D.R.I.S. (Institut du Développement et des Ressources en Informatique Scientifique, project n° 000187).

focus on dicarboxylic acids chemistry representing a crucial scientific issue in organic aerosol particles chemistry [START_REF] Fuzzi | Critical assessment of the current state of scientific knowledge, terminology, and research needs concerning the role of organic aerosols in the atmosphere, climate, and global change[END_REF].