Vincenzo Bianco

Oronzio Manca

Sergio Nardini
email: sergio.nardini@unina2.it

Mario Roma

NUMERICAL INVESTIGATION OF TRANSIENT THERMAL AND FLUIDYNAMIC FIELDS IN AN EXECUTIVE AIRCRAFT CABIN

The objective of the present work is the numerical simulation of the thermal and fluid dynamic field in the cabin of an executive aircraft. A transient analysis is conducted on three and two-dimensional models of an executive aircraft cabin and a comparison of the two approaches is provided. In the model, a global thermal conductance was considered to take into account both the external environment and the fuselage material. The results, provided in terms of graphs and contours, refer to the three and two-dimensional models with inclined top inlets with respect to the vertical direction of an angle of 45°. A good

qualitative agreement between the results of the two and three dimensional models is detected.

Keywords: Airflow, aircraft air conditioning, CFD simulation, aircraft cabin, executive aircraft.

Nomenclature

Introduction

The reaction airplanes must fly to high altitude to furnish, in terms of fuel economy, the best performances. High altitudes are convenient in relationship to the fuel consumptions, but the same cannot be said for the external atmospheric conditions, that necessitate an isolation of the aircraft occupants from the external environment, characterized by low values of pressure and temperature. At an altitude of 11.000 m the air temperature is about -55 °C; the pressure is about one fifth of that with respect to the sea level and the damp is around zero. A system of environmental control is used for protection of passengers and crew in an aircraft from such extreme environmental conditions. The air distribution system is a fundamental component of the environmental control system of an aircraft because it should correctly distribute the conditioned air in the cabin, in order to realize a comfortable environment.

Higher occupant density, complex geometry and a lower supply of air rate per person when compared to buildings, make the design of a comfortable and healthy cabin environment for airplanes, especially for the executive class where top performances are required, very challenging.

In the last fifteen years, many numerical and experimental studies have been carried out, because of the importance of correct air distribution in the cabin. Mo et al. [START_REF] Mo | Application of particle image velocimetry for the measurement of the airflow characteristics in an aircraft cabin[END_REF] used the Particle Image Velocimetry (PIV) to measure the distribution of air in the cabin of an aircraft in the absence of seats. Dechow [START_REF] Dechow | Concentrations of selected contaminants in cabin air of airbus aircrafts[END_REF] and Waters et al. [START_REF] Waters | Measurements of indoor air quality on commercial transport aircraft[END_REF] studied the air quality in the cabin measuring the quantity of contaminants, particulate and the concentration of ozone without studying the distribution of the air in detail. Garner et al. [START_REF] Garner | CFD validation for contaminant transport in aircraft cabin ventilation flow fields[END_REF] conducted a campaign of measure in a Boeing 747 through two ultrasounds anemometers in absence of thermal loads and under conditions of static regime. Lately, Sun et al. [START_REF] Sun | Experimental characterization of airflows in aircraft cabins, Part 1, Experimental system and measurement procedure[END_REF] and Zhang et al.

[6] performed experimental measures through velocimetry with a tide of particles in a

Boeing 767-300 in presence of manikins. Recently, Kuhn et al. [START_REF] Kuhn | Experimental parametric study of forced and mixed convection in a passenger aircraft cabin mock-up[END_REF] experimentally investigated forced and mixed convection in a full scale passenger aircraft cabin mock-up.

Large scale particle image velocimetry (PIV) and temperature field measurements were conducted in a cross sectional plane of the mock-up. Their study demonstrates that the flow field in aircraft cabins is affected by various fluid mechanical phenomena, such as interaction between the supplied air jets, buoyancy forces and thermal plumes interactions.

All these parameters influence the flow field inside the cabin.

In comparison to the experimental investigations, the numerical ones are less expensive and it is possible to simulate situations similar to the reality, particularly for the transient state.

This is especially true for medium size firms which produce executive aircrafts, because it would not be economically viable for them to build a mock-up in an environmental chamber and engage an expert experimental team to perform wide experimental analysis on the environmental control system. Otherwise, for this class of aircraft, the best performances are required, thus the CFD represents a very powerful tool to improve design and performances. For these reasons, in the last ten years, many numerical studies have been published. Aboosaidi et al. [START_REF] Aboosaidi | Computational fluid dynamics applications in airplane cabin ventilation system design[END_REF] were among the first to use the CFD to study the air distribution in the cabin of a commercial aircraft without passengers and thermal loads. A similar study was conducted by Mizuno and Warfield [START_REF] Mizuno | Development of three dimensional thermal, airflow analysis computer program and verification test[END_REF] who also appraised the concentration of carbon dioxide. To take into account the presence of the passengers, Singh et al. [START_REF] Singh | Numerical simulation of airflow in an aircraft cabin section[END_REF] used some heated cylinders placed on the seats.

A section of the Boeing 767-300 was recently investigated by Lin et al. [START_REF] Lin | Numerical simulation of airflow and airborne pathogen transport in aircraft cabins-Part 1, Numerical simulation of the flow field[END_REF][START_REF] Lin | Numerical simulation of airflow and airborne pathogen transport in aircraft cabins-Part 2: Numerical simulation airborne pathogen transport[END_REF] using both the LES (Large Eddy Simulation) and RANS (Reynolds-average Navier-Stokes) models.

Lately Zhang and Chen [START_REF] Zhang | Novel air distribution systems for commercial aircraft cabins[END_REF] used FLUENT with a RNG (Renormalization Group) k -ε model to study the distribution of the air flow, temperature and the concentration of the CO 2 in a section of a commercial aircraft with different kinds of distributions. Bianco et al.

[14] investigated a two dimensional transient model of an aircraft section either with top inlets placed on the sides at an inclination of 45° with respect to the vertical or with top inlets placed in the middle section of the fuselage. The results show two vortex cells inside the cabin in both cases, that are unstable in the first case whereas they are stable in the second one.

Yan et al [START_REF] Yan | Experimental and CFD study of unsteady airborne pollutant transport within an aircraft cabin mock-up[END_REF] presented a study on the airborne pollutant transport within an aircraft cabin mock-up with both experimental and CFD approach. A steady airflow field was simulated first and then it was compared with the experimental data. The study demonstrated how CFD simulations can help to understand the pollutant transport mechanism. They provide more detailed information, which cannot be exclusively obtained by experiments.

Recently Zhang et al. [START_REF] Zhang | Experimental and numerical investigation of airflow and contaminant transport in an airliner cabin mockup[END_REF] investigated airflow and contaminant transport in an occupied cabin by both experimental measurements and numerical simulation with CFD. The numerical simulation used the Reynolds averaged Navier-Stokes equations based on the RNG k-ε model to calculate the air velocity, air temperature and concentration of gas contaminants. Their numerical results quantitatively agreed with the experimental data whereas some remarkable differences exist in airflow distribution.

Most papers present in the literature focus on the CFD airflow simulation and experimental investigation in large passenger aircraft cabins. The aim of the authors is to cover the lack of information regarding the executive class.

The purpose of the present work is the numerical visualization of the thermal and fluid dynamic field in the cabin of an executive aircraft, in order to collect information to improve the air conditioning design. A transient analysis is conducted on two and three dimensional models of the aircraft cabin and a comparison of the two approaches is performed, in order to estimate the trade off between accuracy and velocity of the solution.

In the model, a global thermal conductance is estimated and assigned to take into account both the external environment and the fuselage material. The results refer to the two and three dimensional models with the top inlets inclined with respect to the vertical axis of an angle of 45°. To the authors best knowledge, this represents the first paper which focuses specifically on this special class of aircraft.

It is important to remark that the present analysis is not completely able to capture the effects due to passengers and seats as well as the effect of a different overheads inclination.

Numerical model

Governing equations

The problem under investigation is a two and three dimensional simulation of an executive The flow is assumed incompressible and, depending on the model considered, two or three dimensional, a turbulence model is, however, used to take into account the viscous effects.

In the following equations, "D/Dt" is the material derivative and the symbol " " refers to time averaged quantities.

0    V (1)   T g T p Dt V D             ) ( (2)
  T q q Dt T D Cp      (3)
The governing equations of the fluid flow (1-3) constitute a highly non-linear and coupled partial differential equations system. In fact the density, ρ, is linked to energy equation being a temperature function. The well known Boussinesq approximation is considered to take into account temperature dependence of density.

The resulting equations system should be solved taking into account appropriate boundary conditions; and since the problem is turbulent, closure equations are required to arrive at the solution. Empirical data or approximate models are required to express the turbulent stresses and heat flux quantities of the related physical phenomenon.

In the present analysis, k -ε turbulent model in the RNG form is adopted, because it was broadly validated by other authors in different enclosed environments [START_REF] Zhang | Novel air distribution systems for commercial aircraft cabins[END_REF][START_REF] Zhang | Experimental and numerical investigation of airflow and contaminant transport in an airliner cabin mockup[END_REF][START_REF] Chen | Computational fluid dynamics for indoor environment modeling: past, present and future[END_REF][START_REF] Chen | Comparison of different kε models for indoor air flow computations[END_REF][START_REF] Zhang | Evaluation of various turbulence models in predicting airflow and turbulence in enclosed environments by CFD. Part-2: Comparison with experimental data from literature[END_REF] and, moreover, for the case analyzed in the present paper, to the best of our knowledge, experimental data, to validate the model, do not exist.

Boundary and initial conditions

For the considered aircraft section, inlet air velocity and temperature were fixed; the position of outlet sections was assigned and wall conductance and external temperature were assigned on the walls.

From the right and left overheads (Fig. 1), placed in the top part of the cabin, air is introduced at a velocity of 6.0 m/s and at a temperature of 318 K, while from the right and left side diffusers, positioned in the lower part of cabin, air is introduced at a speed equal to 0.4 m/s and at a temperature of 318 K. The air is expelled to the outside through the outlets corresponding to the corridor.

The other surfaces are walls subjected to convection at a temperature of 245 K with an imposed average conductance equal to 3 W/(m 2 K). The two closure faces are considered symmetrical surfaces, hence the derivatives of all the quantities are equal to zero.

In the case of the two-dimensional approximation, in order to balance the non perfect longitudinal symmetry, a lower inlet air temperature, equal to300 K, was considered. The two overheads in the 2D model were considered like a "strip" injection over all the longitudinal length, while in reality they are two punctual air injections, such that the quantity of hot air supplied to the cabin is higher. To avoid an over heating effect, the temperature is lowered.

The criteria used is to work with constant enthalpy, so that if the air flow rate is higher (strip injections), to balance the enthalpy, the air temperature is lowered. In this way, it is possible to obtain a 2D model which gives a good qualitative approximation of the 3D model, as will be shown in the results.

A time step of 0.1 s is adopted for the transient calculation and it was chosen in order to guarantee the solution convergence.

Initially the air inside the cabin is calm at a temperature equal to 295 K. At t=0 s, the air enters into the cabin through two overheads positioned in correspondence of passengers head, inclined of 45° in comparison to the vertical and from two side diffusers placed in proximity of the passengers feet.

The fuselage is composed of various materials assembled in parallel/series. Particularly, there are: "Microlite 0.4" and "Microlite 0.6" placed in series and both are parallel, with an aluminium frame, which has a structural function, as shown in Fig. 1c.

A further complication to the aforementioned problem has to do with the fact that aluminium and microlite conductivities differ by four magnitude orders, hence a simple series-parallel model is not applicable since the temperature field is not one-dimensional, but two-dimensional [START_REF] Incropera | Fundamentals of Heat Transfer[END_REF].

In order to calculate an average conductivity for aluminium, "Microlite 0.4" and "Microlite 0.6" (Fig. 1c), a finite volume numerical model was implemented. The model consists of a rectangular domain, composed by the aforementioned materials (as reported in Fig. 1c), and in two assigned temperatures imposed on the vertical sides, one cold and the other hot, with adiabatic horizontal sides. The average heat flux was determined by means of the numerical model and, consequently, the average thermal conductivity was calculated. The estimated value for the thermal conductivity (k avg ) is 0.16 W/(mK).

Once the average conductivity is estimated, it is possible to calculate the average conductance, which represents the parameter to impose in the aircraft cabin model.

In order to do this, the fuselage is considered as a hollow cylinder with length equal to the length of the elementary cell, Fig. 1c, composed of a part in aluminium and another in microlite, repeating itself in the space. The following equation is used to calculate the thermal resistance:

avg ext k L r r R      2 ln int (4)
where r ext and r int are equal to 0.96 m and 0.91 m respectively (Fig. 1c and2a), whereas L is equal to 0.185 m (Fig. 1c). Once the thermal resistance is obtained, the conductance can be calculated by making the reciprocal of R and dividing this value by the lateral area of the elementary cell:

L r A ext    2 (5)
This estimated conductance value, U, is equal to 3.0 W/(m 2 K).

Geometrical modelling and grid generation

The transient numerical investigation was carried out on two and three dimensional model of a section of an executive aircraft cabin. In this first study, the cabin was considered void, without taking seats and passengers into account.

The transversal section of the aircraft, shown in Fig. 1a and 2a, was modelled using Gambit

[17] and to obtain the three-dimensional model, it was extruded by 1237 mm along the longitudinal axis, Fig. 1b and2b.

Treated air enters the aircraft cabin from 4 inlets: two inlets on the top in proximity of the passengers head, inclined at an angle of 45° with respect to the vertical direction and two side diffusers in proximity of the passengers feet.

The computational grid was built using tetra elements, refined in proximity of the air inlets and outlets sections. The skew value is considered in order to measure the mesh quality; it varies between 0 and 1 and it should be the lowest possible in order to build a model that guarantees the convergence.

A grid independence analysis has been performed and, for this reason, different grids are tested. The mesh that introduces acceptable errors, with respect to the most refined grid tested, was selected in order to arrive at a correct compromise between time of calculation and accuracy of the results. The accuracy control parameter of the results was the transient average temperature on the cabin wall "Wall 1", Fig. 1b and Fig. 3.

The selected grid is the medium one, which is composed of 476.073 cells and introduces a maximum value of skewness equal to 0.75, while the difference in terms of medium temperature, expressed in degrees Celsius, with respect to the fine grid is about 5%, as reported in Fig. 3.

Numerical method

The system of governing equations, subject to their appropriate boundary conditions, was successfully solved by using the finite volume method. This method is based on the spatial integration of the conservation equations over finite control volumes for each time step, converting the governing equations to a set of algebraic equations.

The algebraic "discretized equations" resulting from this spatial integration process were sequentially solved throughout the physical and temporal domain considered. The computational analysis was performed by employing the commercial software FLUENT.

Conservation equations residuals are considered as convergence indicator and they were set in the following way: 10 -3 for mass conservation, 10 -4 for momentum equation and 10 -6 for energy conservation. During the iterative process the residuals were carefully monitored and the convergence was checked for each time step.

The solution was obtained on a XEON with two 32 bit processor of 2.40 GHz each and 2GB of RAM. It took about five days to simulate 800s.

Results

Results are presented in terms of velocity and temperature for two and three dimensional transient simulations. The quasi steady state in the case of three dimensional model is achieved after about 500 s, as shown in Fig. 3, while for the two dimensional model after about 30 s, as reported in Fig. 4b. The oscillating behaviour shown in Fig. 3 is due to the presence of vortex inside the cabin, as subsequently explained. The oscillations remain and the motion is therefore considered as quasi steady state.

Two-dimensional model results

For the two dimensional model, Fig. 4a shows the vertical component of velocity, v y , as a function of the horizontal coordinate for different times. The profiles are characterized by oscillation during the time. In fact, the absolute minimum continuously changes its position along the profile. This is due to the reciprocal impact of the two warm air jets coming from the overheads. One of them alternatively predominates on the other, such that one is recalled on the upper part of the cabin provoking a recirculation zone, while the other is able to penetrate deeper in the cabin, causing the main vortex. This oscillation causes the instability of the two main vortexes inside the cabin, as shown in Fig. 5. This behaviour is more pronounced in the two dimensional model, because the longitudinal direction is neglected. In the subsequent three dimensional model, it is possible to notice that the longitudinal dimension plays a part in limiting this phenomenon as the jets dissipate their energy longitudinally.

Fig. 4b shows that a quasi steady state, for the two dimensional case, is achieved after about 30 s. In fact, the temperature profiles at 60 s, 90 s and 120 s are practically overlapped.

Moreover, a slight increase of temperature along the profiles is detected in correspondence to maximum velocity magnitude. This happens in one fourth of the cabin width in central part of the section, as highlighted by the shaded areas of Fig. 4a and4b. Temperature is quite uniform inside the cabin, but there is a strong drop near the wall and probably one or more insulating substrates are required.

Three-dimensional model results

In Fig. 6a, temperature values for three dimensional model in overhead section along line 1, as a function of the time, are reported. Temperature decreases with time until the steady state around 500 s, as shown also in Fig. 3. The peak value, about 297 K close to the middle of the section, is due to the two jets which merge in this zone and then start to oscillate. In fact, the maximum value continuously changes its position. Moreover, considering the profiles at 600 s and 800 s, it is detected that, on the left of the maximum, for 600 s the temperature value is higher than the one for 800s, while on the right part, the opposite is observed. This is due to the left jet, at this time, which pushes the right one towards the cabin wall. This behaviour changes with time and consequently, temperature field is never symmetric, even if the geometry and boundary conditions are symmetric.

In Fig. 6b, temperature along line 1, for different sections at 800 s is given. In section 3, which is 0.089 m far from the overhead, the peak due to the jets is visible. The maximum value along the profile is attained at different points for different sections, as highlighted by the shaded area of Fig. 6b. It is concluded that there is jets oscillation also in the longitudinal direction, as clarified from the pathlines, which are omitted for the sake of brevity . However, the zone along the width is between about 0.7 m and 1.0 m.

Temperature is almost uniform inside all the cabin, reaching the value of about 294 K, and this is very important in order to guarantee the passengers thermal comfort. A strong temperature drop is detected just on the cabin wall.

Temperature fields, omitted for sake of brevity, confirm that the temperature is rather uniform in all the section, just a warmer vortex is in the central upper part of the section, which moves alternatively due to the two warm jets coming from the overheads.

The horizontal component of velocity, for different time steps, is presented in Fig. 7a. It is possible to notice some peak values at y=0.4 m with a value of about 0.4 m/s, which represents the jet coming from the side diffuser. For t=600 s, this maximum velocity value is not present, probably due to the vortex on the left side of the cabin which recalls the jet coming from the left side diffuser. Fig. 7b shows the horizontal velocity component at t=800 s along different sections. A peak in section 2, in the lower part of the cabin, is observed. This is due to the warm air jet coming from the left side diffuser, while in section 3, there is a peak value in the upper part of the cabin, which is probably due to overhead jets oscillation along the longitudinal plane, in agreement with Fig. 6b, where, in the same section, the maximum is detected. As for the other sections, it is possible to notice velocity fluctuations, probably due to the flow field instability, which is caused by the vortex movements. This guarantees a uniform temperature in order to ensure passenger comfort.

On the contrary, fluctuating velocity may determine discomfort conditions in the cabin.

In Fig. 7c When the jets meet, they also expand along the longitudinal direction, giving room to a motion in this direction, as shown by means of the path lines coming from the overheads at different times, which are not reported for sake of brevity. These path lines show that the motion field is quite disordered and characterized by many fluctuations, which affect the comfort of the cabin occupants.

However, the presence of these various motion fields tend to generate a more uniform temperature field, because of the high mixing effect.

Conclusions

In the present paper, thermal and velocity fields, in a section of an executive aircraft, were numerically investigated. The transient behaviour of the cabin is analyzed on two and threedimensional models and a global thermal conductance is considered to account for both the external environment and the fuselage materials.

Results of two-dimensional and three-dimensional models show a good qualitative agreement. In fact, all the main phenomena are predicted by the two dimensional model, irrespective of differences from the quantitative point of view. This is true, especially for the velocity field, that in the two dimensional model, does not consider the longitudinal direction and the kinetic energy is dissipated in a smaller domain, leading to higher values of the velocity. In terms of temperature the values are comparable because an enthalpy correction was adopted. It is important to remark that the indications given from the two dimensional model are very effective from a qualitative point of view, especially because they are immediately available. In fact, the numerical solution of a two dimensional model takes some minutes, while the three dimensional one takes some days. The aim is to perform a fast pre-optimization analysis on the two-dimensional model and then upgrade it to a three dimensional one.

Three dimensional model results outline that the thermal field seems adequate to guarantee the passengers comfort. In fact, there is an average temperature of about 294 K and just on the adjacent wall, there is a strong drop. Therefore, the wall insulation should be increased.

On the other hand, the velocity values have a correct magnitude to provide comfort conditions. They are in the order of ±0.1 m/s, except in proximity of the air inlets where they increase. However, they are characterized by fluctuations, which may cause passengers discomfort.

It is important to remark that, this is a first attempt to model thermal and velocity fields inside an executive aircraft cabin, therefore the presented results have limited practical application. Future improvement may be represented by the inclusion of seats and passengers, whose presence undoubtedly affect the flow configuration and by performing an optimization of the overheads inclination, in order to improve the significance of results.

A 2 k

 2 Lateral area of the elementary cell, m 2 Cp Specific heat of the fluid, J/(kg K) g Acceleration due to the gravity, m/s Thermal

  aircraft cabin in transient mixed convection. The transient problem was solved by the commercial code FLUENT 6.2 [17].

  , the vertical velocity component, for different times, is reported. The figure shows a minimum, which corresponds to the jets impact zone inside the cabin and fluctuations around x=0.8 m are observed at different times, as evidenced by the shaded area of Fig.7c. On the left and on the right of the minimum point, the velocity alternatively increases and decreases. Particularly, when the minimum is placed on the left part of the section, the average velocity value on the right part is higher and vice versa. This oscillating effect is due to the impact of the two jets: one jet tends to push the other. When the minimum velocity value moves on the left, it means that the right jet is prevailing on the left one and the opposite is the case when the minimum velocity value moves on the right. The vertical velocity component at t=800 s is reported in Fig.7dat different sections, where a minimum point at around x=0.8 m is observed, which corresponds to the jet penetration in the cabin. Moreover, this minimum point is present for section 2 and 3. This is due to the fact that the jets oscillate also in the longitudinal plane. In the other sections, the vertical velocity component results quite oscillating between -0.2 m/s and 0.1 m/s. The flow oscillations inside the cabin are highlighted by the velocity vectors fields at different times in Fig. 8. At t=200s the two jets merge into a quasi symmetrical flow in the middle section and four macro vortex cells are detected, two in each side. The merged flow moves toward the right side at t=400s, Fig. 8a, and the lower vortex cell in the left side expands. This effect should determine a depression in the left zone and the right jet push the left jet toward the left side. The merged flow moves toward the left part of the section, as shown in Fig. 8b. At t=600s and 800s, the right side has a lower pressure and flow moves toward the right side. Therefore, a sort of flip-flop flow is detected inside the cabin.

Figure 1 :

 1 Figure 1: (a) Two dimensional domain, (b) Three dimensional domain and (c) Schematic

Figure 2 :

 2 Figure 2: (a) References lines "Position 1" (P1) and "Line 1" (L1) for velocity and

Figure 3 :

 3 Figure 3: Average transient temperature on "Wall 1" used to assess grid independence.

Figure 4 :

 4 Figure 4: Velocity and temperature profiles for 2D model along line 1: (a) transient velocity

Figure 5 :

 5 Figure 5: Stream function fields in 2D model after: (a) 30s and (b) 60s.

Figure 6 :

 6 Figure 6: Temperature profiles in 3D model: (a) transient temperature behaviour in the

Figure 7 :

 7 Figure 7: Horizontal velocity component along position 1 (a) in section 2 at different time

Figure 8 :

 8 Figure 8: Velocity field in section 2 at different times (a) t=400 s and (b) t=800 s.


G Soldati

L Boschi

F Deschamps

D Giardini

Inferring radial models of mantle viscosity from gravity (GRACE) data and an evolutionary algorithm

Keywords: mantle rheology, inverse theory, viscous flow, gravity, tomography

Convective flow in the mantle can be thought of (and modeled) as exclusively driven by density heterogeneities in the mantle itself, and the resulting lateral variations in the Earth's gravity field. With this assumption, and a model of mantle rheology, a theoretical relationship can be found between 3-D mantle structure and flow-related quantities that can be measured on the Earth's surface, like free-air gravity anomalies. This relationship can be used to set up an inverse problem, with 1-D mantle viscosity as a solution. In the assumption that seismic velocity anomalies be of purely thermal origin, and related to density anomalies by a simple scaling factor, we invert the large-scalelength component of the above-mentioned measurements jointly with seismic observations (waveforms and/or travel times) to derive an accurate 5-layer spherically symmetric model of upper-and lower-mantle viscosity. We attempt to account for non-uniqueness in the inverse problem by exploring the solution space, formed of all possible radial profiles of Earth viscosity, by means of a non-deterministic global optimization method: the evolutionary algorithm (EA). For each sampled point of the solution space, a forward calculation is conducted to determine a map of gravity anomalies, whose similarity to GRACE is then measured; the procedure is iterated to convergence, according to EA criteria.

The robustness of the inversion is tested by means of synthetic tests, indicating that our gravity data set is able to constrain less than 6 radial layers, each with uniform viscosity. Independently of the tomographic model or the scaling factor adopted to convert seismic velocity into density structure, the EA optimization method finds viscosity profiles characterized by low-viscosity in a depth range corresponding to the transition zone, and relatively uniform elsewhere.

Viscosity is then typically estimated after solving coupled flow and gravitational potential equations for instantaneous deformation (flow, surface deformation, geoid) or time-dependent deformation (relative sea-level, plate motions); the former approach is sensitive to relative viscosity variations, while the latter allows also an estimate of absolute viscosity. In both cases, mantle viscosity is inferred by fitting modelled signals to various types of observations: relative sea level and variations in the Earth's rotational parameters for post-glacial rebound studies, dynamic topography, geoid and plate velocities for mantle convection analysis.

The surface observables of post-glacial rebound, geoid and dynamic topography have provided only first-order constraints on the radial viscosity structure of the mantle: while geoid/dynamic topography studies suggest that mantle viscosity increases by a factor of 30 or more from the basis of the lithosphere to the core-mantle boundary (e.g. Hager & Richards 1989), most post-glacial rebound studies [START_REF] Haskell | The motion of a fluid under a surface load 1[END_REF][START_REF] Peltier | Glacial isostatic adjustment II. The inverse problem[END_REF][START_REF] Peltier | Postglacial variations in the level of the sea: implications for climate dynamics and solid-earth geophysics[END_REF][START_REF] Mitrovica | Constraints on Mantle Viscosity Based Upon Post-Glacial Uplift Data From the Hudson Bay Region[END_REF][START_REF] Kaufmann | Glacial isostatic adjustment and the radial viscosity profile from inverse modeling[END_REF]) favour a moderate increase in viscosity at the upper-to-lower mantle discontinuity. These inferences are still subject of a contentious debate, and to reconcile convection-based and post-glacial reboundbased estimates, joint inversions of these two kinds of data have been performed, obtaining profiles with an overall increase in viscosity towards the lower mantle (e.g. [START_REF] Forte | New inferences of mantle viscosity from joint inversion of long-wavelength mantle convection and post-glacial rebound data[END_REF][START_REF] Mitrovica | A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic adjustment data[END_REF]).

Recent progress, including the introduction of compressibility [START_REF] Corrieu | Mantle dynamics and geoid Green functions[END_REF][START_REF] Panasyuk | Understanding the effects of mantle compressibility on geoid kernels[END_REF][START_REF] Forte | New inferences of mantle viscosity from joint inversion of long-wavelength mantle convection and post-glacial rebound data[END_REF], the evaluation of the performance of non-Newtonian rheology (Wu, 1992 Mantle circulation models that simultaneously predict seismic (P-and S-wave velocities) and geodynamical data (free-air gravity anomaly) have been shown 2 A c c e p t e d M a n u s c r i p t to be particularly good at fitting the latter (Forte et al., 1994). Here we model viscous flow in the Earth on the basis of a wide range of possible viscosity profiles, and attempt to identify the profiles for which the modeled viscous flow does the best job of predicting observed free-air gravity anomalies from GRACE [START_REF] Tapley | GGM02 -An improved Earth gravity field model from GRACE[END_REF]. We define a-priori density (ρ) models needed in mantle-flow calculations on the basis of seismic tomographic ones, scaled by a factor depending only on depth. This requires the assumption that the relative thermal and compositional contributions to seismic anomalies are the same everywhere at any given depth.

The interpretation of long-wavelength geoid/gravity anomalies in terms of mantle convection has a long history, starting from the pioneering works of [START_REF] Hager | A simple global model of plate dynamics and mantle convection[END_REF], [START_REF] Ricard | Geoid heights and lithospheric stresses for a dynamic earth[END_REF], [START_REF] Richards | Geoid anomalies in a dynamic Earth[END_REF], to the recent ones by [START_REF] Panasyuk | Models of Isostatic and Dynamic Topography, Geoid Anomalies, and Their Uncertainties[END_REF], [START_REF] Forte | High Viscosity Deep Mantle Flow and Thermochemical Structure Inferred From Seismic and Geodynamic Data[END_REF], [START_REF] Kaban | Importance of lateral viscosity variations in the whole mantle for modelling of the dynamic geoid and surface velocities[END_REF]. With few exceptions, authors have adopted a viscousflow theory which assumes mantle rheology to be represented in terms of an effective viscosity varying only with depth. Although mantle viscosity is likely to have lateral variations, 1-D viscosity profiles have been shown, to be an adequate representation of the horizontally-averaged mantle viscosity structure [START_REF] Moucha | Lateral variations in mantle rheology: implications for convection related surface observables and inferred viscosity models[END_REF]; conversely, the effect of lateral variations in viscosity is thought to be reflected almost exclusively in the small scale (high harmonic degree) component of the gravity field. The same is true of upper-mantle viscosity, while low-degree coefficients of the gravity field are more sensitive to the lower mantle (Richards & Hager, 1989;[START_REF] Forte | The kinematics and dynamics of poloidaltoroidal coupling in mantle flow: the importance of surface plates and lateral viscosity variations[END_REF]. Our goal is to identify a 1-D, whole-mantle viscosity model, and for this reason we neglect (except for a test in section 3.3) the high-degree component of gravity data and, consequently, of seismic models.

The goal of our contribution is not only to determine the viscosity profile of the mantle, but also to estimate the ability of gravity anomalies to resolve the radial distribution of viscosity as inferred from tomography and flow models.

We tested and used an evolutionary algorithm to invert various tomographic models for a radial profile of the mantle viscosity, and found that all the most likely viscosity profiles predict transition-zone (410-660 km depth) viscosity to be lower than in the uppermost and lower mantle.

Theory

Viscosity, gravity, and mantle flow

The relative radial variations in mantle viscosity can be determined from gravity measurements. An analytical theory of mantle flow [START_REF] Ricard | Geoid heights and lithospheric stresses for a dynamic earth[END_REF] A c c e p t e d M a n u s c r i p t [START_REF] Forte | Plate tectonics and aspherical Earth structure: The importance of poloidal-toroidal coupling[END_REF][START_REF] Forte | Viscous flow models of global geophysical observables: 1. Forward problems[END_REF] provides geoid kernels (G l , with l denoting the harmonic degree) given an average density profile (here, PREM), and a prescribed viscosity profile. The surface gravity anomalies (δg m l , with m harmonic order) are modeled by radial integration (from top of the outer core to the surface) of the ρ anomalies (δρ m l ), modulated by the geoid kernels.

For each harmonic, gravity anomalies are thus given by (e.g. [START_REF] Forte | Plate tectonics and aspherical Earth structure: The importance of poloidal-toroidal coupling[END_REF])

δg m l (θ, φ) = k l -1 2l + 1 r G l (ν/ν 0 , r)δρ m l (r, θ, φ) dr, (1)
where k is a constant that depends on the Earth's radius, surface gravity acceleration, and average mantle density, ν 0 is a reference value for viscosity [START_REF] Forte | Viscous flow models of global geophysical observables: 1. Forward problems[END_REF]) and the integration is carried out over the entire depth of the mantle. The multilayer approximation is employed, in which the viscosity is assumed to be constant within each layer and discontinuous at the layer boundaries; it should be noted that only the depth variation of relative viscosity ν/ν 0 is needed to compute the geoid kernels.

We scale ρ anomalies from shear-velocity (v S ) ones,

δρ m l (r, θ, φ) = ζ(r)δv S m l (r, θ, φ), (2)
where the scaling factor ζ is defined by

ζ(r) = δ ln ρ(r, θ, φ) δ ln v S (r, θ, φ) . (3)
The kernels G l are calculated following the approach of [START_REF] Forte | Viscous flow models of global geophysical observables: 1. Forward problems[END_REF], who expanded in terms of generalized spherical harmonics the constitutive equation, the conservation of mass and momentum, and solved for the poloidal flow using the method of propagator equations. The constraints arising from the observed geometry of rigid surface plates are included in a dynamically consistent manner by means of the buoyancy projection method [START_REF] Forte | Viscous flow models of global geophysical observables: 1. Forward problems[END_REF], in which the motions of surface plates are predicted (being coupled to the underlying mantle flow) rather than imposed. The plate geometries and corresponding projection operators are represented in terms of spherical harmonic basis functions up to degree l = 8, in order to reduce the effect of uncertainties in tomography (larger for higher spherical harmonic degrees), and because higher-degree geoid kernels are sensitive to heterogeneous structure in the upper mantle only, while we want to integrate ρ anomalies over the whole mantle.

Evolutionary algorithms

In this study we conduct a number of inversions based on different tomographic models and different values of the ratio between density and seismic velocity, and carry out a comparative evaluation of the resulting viscosity profiles. Owing to the above factors, and to nonlinearity (geoid kernels depend on viscosity itself), the problem of finding viscosity profiles from geophysical observables does not have a unique solution. We attempt to account for non-uniqueness in the inverse problem by exploring the solution space, formed of all possible radial profiles of Earth viscosity, by means of a stochastic optimization method:

the evolutionary algorithm (EA).

The first EAs, or optimization algorithms based on ideas from evolutionary theory, were conceived and implemented by [START_REF] Reichenberg | Evolutionsstrategie: Optimierung technisches Systeme nach Prinzipien der biologischen Evolution[END_REF]. Authors in earth sciences have already used EAs to find viscosity profiles from gravity data [START_REF] King | Radial models of mantle viscosity: results from a generic algorithm[END_REF][START_REF] Kido | Mantle viscosity derived by genetic algorithm using oceanic geoid and seismic tomography for wholemantle versus blocked-flow situations[END_REF], but our analysis represents an improvement in that it takes advantage of the increased power of modern calculators to explore in more detail the solution space. In particular, we have been able

to quantify the resolution limit of gravity data, finding that no more than 5 independent parameters (uniform layers) describing viscosity can be reliably constrained. In addition, a systematic test of the setting parameters of our EA A c c e p t e d M a n u s c r i p t lead to the conclusion that while the choice of a population of 100 individuals is appropriate for this problem, the number of generations over which King (1995) based his analysis is not sufficient to grant the stability of the solution;

Figure 1 shows that after 100 generations the fit of the solution model to the data may still be improving: we thus increased the number of generations from 100 to 500.

EAs use the idea of "survival of the fittest", to perform an iterative, multidimensional search for an optimal value of a given cost function. The mutation rate is a key parameter: if it is too low, the algorithm may converge prematurely to a local optimum, the EA failing to explore uniformly the space of parameters. In contrast, a high mutation rate may lead to slow or no convergence (an EA with high mutation is practically equivalent to a Monte Carlo algorithm). [START_REF] Charbonneau | A User's guide to PIKAIA 1.0[END_REF] suggest that a good compromise between allowing for new solutions and losing track of already identified ones is achieved by starting the EA run with a low mutation rate, and then allowing the mutation rate to grow as convergence is approached.

Although, ideally, the solution found by the EA should be independent of it, the choice of a specific fitness function might also play an important role in the speed and efficiency of the algorithm. We experimented with different cost functions (variance reduction, correlation), finding indeed rather similar solution models. To obtain maximum variability in best-fitting models, i.e., to best differentiate solution models with relatively similar fit as we progressively

-i=1 (δg i mod -δg i obs ) 2 i=1 (δg i obs ) 2 , ( 4
)
where δg i mod and δg i obs are the modeled and observed gravity anomalies, respectively, at the point i of a grid covering the Earth's surface. In a set of preliminary tests, we have verified that the cost-function as defined by eq. ( 4) results in the most effective convergence. Other cost functions that we have experimented with, including correlation and variance reduction (without exponential) did not allow to discriminate between close minima, beyond a certain refinement level.

We define an initial population, consisting of 100 randomly generated viscosity profiles. In most of our runs of the EA, convergence is achieved after roughly 100-300 generations (Figure 1a, red dots). The choice of the mutation rate adjustment (differential fitness) is reflected by the trend of blue dots in Figure 1a, representing the average fitness of the population for each generation:

average fitness achieves a maximum at 10-20 generations, then decreases with increasing fitness of the best solution.

We set the total number of generations to 500, corresponding to 50, 000 forward computations total; completing this task takes about 48 hours on a dual 2.7 Ghz PowerPC. Performance depends on the number of free parameters in the inverse problem (i.e. number of uniform viscosity layers), on the maximum considered harmonic degree, and on the precision chosen for the variables.

Since gravity data have little sensitivity to changes in viscosity larger than three orders of magnitude [START_REF] King | Radial models of mantle viscosity: results from a generic algorithm[END_REF], we used single precision (4 bytes) floating points variables, allowing approximately 7 digits of accuracy. We randomly generate an 'input' profile of mantle viscosity, and use our mantle flow model to predict the corresponding gravity anomaly map. We then use the resulting, 'synthetic' gravity anomaly map as the database to be inverted via the EA. The correlation between output and input model, shown in Figure 3 for 60 independent synthetic tests, is a measure of the accuracy and resolution of our method. We conducted 10 synthetic tests with 2-layer viscosity models, 10 with 4-layer models, and so on with 6-, 8-, 10-and 12-layer models. Figure 3 shows that the non-uniqueness of the problem grows quickly with the number of inversion parameters. If the unknown viscosity profile is parameterized in terms of more than 5 uniform layers, the chance of converging to a wrong solution is high. We use the EA, as described above, to identify the best-fitting 5-layer model of relative changes in mantle viscosity, and, again, scale tomography model S20RTS to define an a-priori density map.

We first explore the influence of the population size, the number of generations and the seed used to initialize the EA on the inversion results. Figure 4a shows the best-fitting viscosity profiles derived from runs of the EA with population sizes of 10, 50 and 100 individuals, evolved for 500 generations. Due to its stochastic nature, different runs of the EA inversion yield slightly varying results, but the important features (2-orders of magnitude viscosity jump at 410 km; smaller but significant jump at 1200 km) remain stable. Running the EA with 100 individuals for 100, 500 and 1000 generations, we obtain almost identical viscosity profiles (Figure 4b). Inverting the same data with same population size and number of generations, but different seed, we find approximately the same radial viscosity profile (Figure 4c). In all these cases, variance reduction (Figure 4d, 4e, 4f) converges to approximately the same maximum.

We invert, again, gravity anomalies from GRACE starting from v S model S20RTS and assuming a density-to-velocity scaling as in Figure 7b (red line).

We repeat the experiment varying the number of constant-viscosity layers from 2 to 12. The resulting viscosity profiles, shown in Figure 5, closely resemble the ones found in the other inversions of this Section, characterized by relatively low viscosity at depths corresponding to the mantle transition zone.

Concerning the fitness to the data, the gravity anomalies computed in five out of the six cases reduce the variance of about 45 to 50%. Conversely, we found no 2-layer model that reduces the variance at all. We infer that at least two viscosity discontinuities in the mantle are required to explain the gravity data in consideration, and 2-layer models can be rejected a-priori.

We run the EA with several different parameterizations, characterized by the same total number (5), but different depth ranges, of uniform viscosity layers.

We show the results in Figure 6. Independently of parameterization, solutions for the profiles in Figure 6a and 6b, suggesting that a radial viscosity structure with finer parameterization in the upper mantle is more consistent with the observed gravity field.

Effect of the short-wavelength component of tomography

The neglect of harmonic degrees >8 is justified by the goal of identifying a whole-mantle viscosity profile, while relatively high harmonic degrees are mostly, if not only, sensitive to the upper mantle. We test, however, the possibility that the whole-mantle viscosity profile we find be perturbed by accounting for degrees >8. We repeat the exercise described in section 3. 1997). Fewer authors found evidence that viscosity in the same region might be anomalously high [START_REF] Ricard | Mantle heterogeneities, geoid, and plate motion: a Monte Carlo inversion[END_REF][START_REF] Spada | Viscoelastic response of a hard transition zone: effects on post-glacial uplifts and rotational signatures[END_REF][START_REF] King | Radial models of mantle viscosity: results from a generic algorithm[END_REF].

Assumptions on the scaling factor

Since the choice of the velocity-to-density scaling factor may impart a bias to our results, and since the ones used here do not incorporate any mineral physics 5%. This result confirms our earlier findings (Figure 9b ).

Viscosity from inversion of gravity data and density models

Albeit commonly used [START_REF] Forte | Geodynamic Evidence for a Chemically Depleted Continental Tectosphere[END_REF][START_REF] Deschamps | Anomalies of temperature and iron in the uppermost mantle inferred from gravity data and tomographic models[END_REF], the procedure of estimating Earth's density via a depth-dependent scaling factor applied to seismic velocity models is, at least to some extent, inaccurate: lateral ρ anomalies directly observed from, e.g., normal-mode data are both uncorrelated with [START_REF] Resovsky | Using probabilistic seismic tomography to test mantle velocity-density relationships[END_REF], and too large with respect to [START_REF] Ishii | Normal-Mode and Free-Air Gravity Constraints on Lateral Variations in Velocity and Density of Earth's Mantle[END_REF][START_REF] Trampert | Probabilistic Tomography Maps Chemical Heterogeneities Throughout the Lower Mantle[END_REF]) seismic anomalies, for the scaling-factor approach to be valid. We replace the v S velocity models used so far with the ρ models provided by 2002) objected that density cannot yet be constrained in this way, because the sensitivity kernels for density are much smaller than those for velocities, and because the least-squares inversions conducted in this kind of studies require the use of a starting model, the choice of which is critical for the reliability of the results.

In Figure 12 we compare the best viscosity profiles found via EA on the basis of the ρ models, with those resulting from the corresponding v S models TRP246 and SPRD6. The difference with the profiles derived from velocity models (plotted for comparison in Figure 12a), is striking: viscosity increases almost monotonically with depth for the profile corresponding to ρ model TRP246, and the low-viscosity transition zone that characterizes all our solution viscosity models becomes much less pronounced in the profile derived from ρ model SPRD6. It is remarkable that variance reduction achieved by density with respect to velocity models drops from 57% to 10% for SPRD6 and becomes negative for TRP246, even though low even degrees of gravity data are appropriately fit by [START_REF] Ishii | Normal-Mode and Free-Air Gravity Constraints on Lateral Variations in Velocity and Density of Earth's Mantle[END_REF]. We explain this discrepancy in terms of the different approach used here with respect to the studies of Trampert et al. ( 2004) and [START_REF] Ishii | Normal-Mode and Free-Air Gravity Constraints on Lateral Variations in Velocity and Density of Earth's Mantle[END_REF], to establish a relationship between mantle flow and observations of gravity. We account for mantle flow explicitly (e.g., [START_REF] Richards | Geoid anomalies in a dynamic Earth[END_REF], while those authors do it by allowing for deflections of the internal boundaries. In the past, it has been assumed that the two approaches are equivalent, but we believe that this assumption must be reevaluated. We show in Figure 13 The corresponding variance reduction of the inverted gravity data amounts to 87%, much higher than achieved in earlier inversions.

CMB topography

Undulations of the CMB are generally believed to be the result of radial stresses generated by convective mantle flow induced, in turn, by lateral variations in density throughout the mantle. We compute here the topography of the CMB from the same v S and ρ models as in section 4.

The spherical harmonic coefficients δb m l of flow-induced CMB topography are related to density perturbations δρ m l by topography kernels B l via an equation similar to [START_REF] Anderson | Transition region of the Earth's upper mantle[END_REF],

δb m l (θ, φ) = 1 ∆ρ cm B l (ν/ν 0 , r)δρ m l (r, θ, φ) dr (5)
( [START_REF] Forte | Seismic-geodynamic determination of the origin of excess ellipticity of the core-mantle boundary[END_REF], where ∆ρ cm = -4.43 Mgm -3 is the density jump across the CMB according to PREM [START_REF] Dziewonski | Preliminary Reference Earth Model[END_REF], and the integration is done from the radius of the CMB to the Earth's surface. The topography kernels B l are calculated, as for the gravity anomalies, in the degree range l = 2 through 8 and, like geoid kernels, they implicitly depend on the (relative) viscosity profile of the mantle (ν/ν 0 ). As before, the harmonic coefficients δρ m l are found from a scaled velocity model.

We show in Figure 15 CMB topography computed directly from ρ models display a more complex pattern, and have amplitude three times bigger, though very similar to each other both in shape and in amplitude. Again, differences between the topography predicted by v S and ρ result from the fact that the two distribution are not correlated (e.g., [START_REF] Trampert | Probabilistic Tomography Maps Chemical Heterogeneities Throughout the Lower Mantle[END_REF].

The dynamic topography at the CMB is not directly observable from surface data, with the exception of the component δb 0 2 of the CMB topography, called excess or dynamic ellipticity, which can be inferred via VLBI measurements of the period of the Earth's free-core nutation. The most recent inferences [START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonridgid Earth and insights into the Earths interior[END_REF][START_REF] Mathews | What do nutations tell us about the Earth's interior?[END_REF] suggest a value closer to 0.4 km, rather than 0.5 km as determined in the earlier study by [START_REF] Gwinn | Geodesy by radio interferometry: studies of the forced nutations of the Earth. 2. Interpretation[END_REF].

The values of δb 0 2 we obtain (Table 1) on the basis of the different v S and ρ models are about three times larger, having absolute value bigger than 1.5

km. The poor fit to the observations of CMB ellipticity may be justified by the fact that our mantle flow models are constrained to only fit the free-air gravity data; other potential reasons are a poor velocity-to-density scaling in the lower mantle and a poorly resolved viscosity at depth. (i) Several high-pressure mineral-physics studies [START_REF] Smyth | Beta-Mg2SiO4: A potential host for water in the mantle[END_REF][START_REF] Kawamoto | Experimental evidence for a hydrous transition zone in the early Earth's mantle[END_REF][START_REF] Kohlstedt | The solubility of water in alpha, beta and gamma phases of (Mg,Fe)2SiO4[END_REF] Despite the many evidences of significant amounts of water in the transition zone, the actual content of water is still poorly constrained, being estimated to range between 0.1 wt% and 3 wt% [START_REF] Bercovici | Whole-mantle convection and the transitionzone water filter[END_REF]. Also, the effect of water on material properties is not clear, even if it is known that it controls the strength and deformation mechanism of minerals [START_REF] Kavner | Elasticity and strength of hydrous ringwoodite at high pressure[END_REF] and thus the rheology of rocks [START_REF] Karato | Plastic deformation of silicate spinel under transition-zone conditions of the Earths mantle[END_REF]. Since viscous deformation is a macroscopic form of creep depending on the presence of defects in the lattice structure, and since water increases the number of defects within a crystal, it enhances diffusion rates and this should decrease viscosity. While anomalously low values of viscosity in the transition zone are a robust result, some of our findings cast doubts on certain aspects of the approach we followed. In particular, we have illustrated in section 4.3 (Figures 12 and13)

the disagreement between density-based and velocity-based modeling results.

We have explained it as the consequence of a discrepancy between the approach followed here, where mantle flow is modeled explicitly (e.g., [START_REF] Richards | Geoid anomalies in a dynamic Earth[END_REF], and that of, e.g., [START_REF] Ishii | Normal-Mode and Free-Air Gravity Constraints on Lateral Variations in Velocity and Density of Earth's Mantle[END_REF] [97] Zhang, S., Christensen, U.R, 1993. 1 Predicted excess CMB topography (km) obtained from some v S and ρ models 25 M a n u s c r i p t

A c c e p t e d M a n u s c r i p t 1 Introduction

 1 The rheology of the Earth is of central importance for understanding both the Earth's transient deformation and long-term mantle dynamics. Presentday estimates of mantle viscosity are based on experimental studies of creep mechanisms in mantle minerals and on the analysis of geophysical observations of the Earth's response to surface and internal loading: mantle convection observables (timescale ∼ 10 6 yr), post-glacial rebound data (∼ 10 3 yr), and post-seismic relaxation (1-100 yr) following major earthquakes.

  ; Dal Forno et al., 2005), and the effects of laterally varying (3-D) viscosity structure (Richards & Hager, 1989; Zhang & Christensen, 1993; Kaufmann & Wu, 2002; Paulson et al., 2005, Moucha et al., 2007) have not clarified the question, and the only point of general agreement is that the lower mantle is more viscous than the upper mantle.

  A typical EA requires a genetic representation of the solutions (in general, as arrays of bits), which play the role of individuals in a population. The algorithm starts from a random population whose individuals are selected according to their fitness, and the best are used to form a new population, likely to be "more fit". Couples of parent chromosomes generate offspring by means of crossover and mutation. This procedure is repeated until a given maximum number of generations is reached, or convergence achieved. EAs are helpful because they can rapidly locate good approximate solutions to all types of problems, requiring no smoothness assumptions on the fitness function or its domain, and because of their robustness in finding global maxima in the presence of many local maxima. Furthermore, EAs are naturally parallel, thus allowing an easy optimization of machine resources.We use[START_REF] Charbonneau | A User's guide to PIKAIA 1.0[END_REF], freely available PIKAIA implementation of the EA. PIKAIA incorporates two basic evolutionary operators: uniform one-point crossover, and uniform one-point mutation. The mutation rate (i.e. the chance that a random variation in an individual's traits occurrs, independent of those of the parents) can be dynamically adjusted during the evolution, using either the linear distance in parameter-space or the difference in fit between the best and median solutions in the population.

6 A

 6 c c e p t e d M a n u s c r i p t refine our search, we chose to use as cost function the exponential of the variance reduction, or exp 1

Figure 2

 2 Figure 2 illustrates the sequence of steps in our algorithm: at each generation, tomographic anomalies are translated into ρ anomalies, which are then used to compute the gravity field at the Earth's surface associated with each viscosity profile in the population. The best-fitting viscosity profiles are then combined by the EA to identify a new, more fit population (a new generation), and the whole procedure is iterated.

A c c e p t e d M a n u s c r i p t 3

 3 Figure 7b (red curve). The scaling factor was obtained by inverting GRACE gravity data and S20RTS tomographic model with an input viscosity profile taken from Mitrovica & Forte (1997). The scaling factor ζ(r) is positive throughout the mantle, except for transition zones, where velocity anomalies are mainly compositional in origin. A thermal origin of anomalies, in fact, requires ζ to be positive, as, for any fixed composition, a perturbation in temperature causes perturbations of equal sign in density and v S .

8

 8

A c c e p t e d M a n u s c r i p t 3 . 2

 32 Testing the effects of different parameterization strategies and evolutionary regimesWe next replace synthetic data with true, free-air gravity anomalies from global Earth gravity model GGM02[START_REF] Tapley | GGM02 -An improved Earth gravity field model from GRACE[END_REF], based on the analysis of 363 days of GRACE (Gravity Recovery and Climate Experiment) in-flight data. Harmonic coefficients up to degree 160 are available, but we only consider degrees ≤ 8, consistently with our decomposition of the tomographic models.

  tend to be characterized by low viscosity in the second shallowest layer, and/or transition zone. The jump in viscosity found between 410 and 660 km depth ranges between 1 and 3 orders of magnitude. Variance reduction is maximum A c c e p t e d M a n u s c r i p t

4

  Figure7b), calculated from various tomographic models and an input viscosity profile selected fromMitrovica & Forte (1997). The three mantle viscosity profiles resulting, after running the EA, from the different tomographic models and scaling factors are shown in Figure7a. All profiles have approximately the same depth dependence, with important viscosity jumps at 410 and 660 km depth. For each run of the EA, corresponding to a certain tomography/density model, we also visualize in Figure7c-e the spread of the population, computing the mean and standard deviation of all viscosity profiles with fit above 53% (S20RTS); 40% (TRP246); 54% (SPRD6) (standard deviations are represented by gray intervals). The three thresholds have been chosen to always correspond to ∼ 10, 000 solution profiles. In all three cases, the best-fitting profiles are all very similar to each other, and different families of viscosity profiles fit the data equally well.

Figure 10 shows

 10 Figure10shows viscosity profiles averaged over models that fit the data best than a prescribed value (30, 40 and 50% in plots a, b, and c, respectively). All the profiles refer to the inversions of gravity data with the (scaled) tomographic model S20RTS. Even models with relatively low fit include a low-viscosity

  constraints, we try to use alternative scalings such as the ones based on laboratory experiments[START_REF] Karato | Origin of lateral heterogeneity of seismic wave velocities and density in the deep mantle[END_REF][START_REF] Cammarano | Inferring upper mantle temperatures from seismic velocities[END_REF]. Given that mineralogy-derived scalings between velocity and density are still subjected to a lot of uncertainties, we invert the gravity data (degrees 1 through 8) with v S velocity model SMEAN[START_REF] Becker | A comparison of tomographic and geodynamic mantle models[END_REF] and a scaling taken from Simmons et al. (2007) (see Figure 9a, blue line), selected among the ones proposed by Karato & Karki (2001) on the basis of the fit to a set of combined seismic and convection-related observables. The resulting viscosity profile is displayed in Figure 9b and does not differ significantly from the ones obtained with classical scalings, confirming the robustness of our results. The gravity anomalies computed with this viscosity achieve a variance reduction of 47.7%.We then attempt to account for the difference between sub-continental and sub-oceanic mantle, revealed by seismic tomography some 40 years ago (Jordan, 1975). The high-velocity roots below continents, absent below oceans (see[START_REF] Romanowicz | Global mantle tomography: Progress Status in the Past 10 Years[END_REF], for a review), are balanced by differences in the respective chemical composition. Here, we have computed radial models of ζ for oceans and continents separately. To define oceanic and continental areas, we have constructed a continent-ocean function derived from the 3SMAC tectonic regionalization (Nataf andRicard, 1996). The sub-continental and sub-oceanic scaling factors (red and green lines in Figure9a) are significantly different at A c c e p t e d M a n u s c r i p t depths shallower than 260 km, with the oceanic one negative at depths up to 80 km. Again we invert gravity data up to degree 8 based on v S velocity model, to find the viscosity profile of Figure 9c, with variance reduction of 46.

  how gravity anomalies computed on the basis of viscosity profiles from Figure 12 compare to GRACE data. While the A c c e p t e d M a n u s c r i p t v S -based viscosity profile of Figure 13a reproduces the data relatively well, the ρ-based results are in fact completely off. An alternative density model has been derived by Simmons et al. (2007), based upon seismic travel-time data, and geodynamic observations including dynamic topography, gravity, plate motions and CMB ellipticity. Simmons et al.'s (2007) approach also implicitly accounts for both thermal and compositional buoyancy effects on mantle flow. We have repeated our inversion experiment assuming density structure as mapped by Simmons et al. (2007). The resulting viscosity profile, shown in Figure 14, confirms the presence of a low-viscosity zone in the upper-to-lower mantle transition zone. It differs from those of the previous sections in the lower mantle, where it is characterized by lower values of relative viscosity, and by a pronounced viscosity jump at 2000 km depth.

  the total CMB topography obtained from v S models S20RTS, TRP246, SPRD6 (Figure15a, 15b, 15c), and that obtained from ρ models TRP246 and SPRD6 (Figure15d, 15e). The viscosity profiles implemented correspond, for each v S or ρ model, to the best-fitting profiles found in Sections 3 and 4 inverting that same model.

A c c e p t e d M a n u s c r i p t

  The v S -derived topographies are in close agreement with most published results (Morelli & Dziewonski, 1986; Forte et al, 1995; Obayashi & Fukao, 1997), characterized by a ring of depressions over the Americas, Eastern Asia, and Australia. CMB deflections based on v S model SPRD6 compares well with those found by Forte et al. (1995), but differ slightly in amplitude. Maps of

6

  Discussion and conclusionsWe have applied the evolutionary algorithm technique to identify a spherically symmetric model of viscosity in the Earth's mantle from global observations of free-air gravity anomalies in the degree range l = 2 through 8. We modeled perturbations in the Earth's gravity field induced by density heterogeneities via a viscous flow model, with no a priori barrier for the vertical flux at the 660 km discontinuity. This approach allows to derive the depth-dependence of relative viscosity, constraining its value uniquely in up to five uniform layers. The solutions we obtain on the basis of v S models S20RTS or SMEAN are consistent with classical estimates of the upper-to-lower-mantle viscosity jump. Additionally, they are characterized by a transition zone less viscous than the uppermost mantle by 2 or 3 orders of magnitude. This feature is parameterization-independent, and is shared by the viscosity profile we find based on the density model of Simmons et al. (2007). The found softening of transition zone minerals could be related to various processes: (i) high content of water (van der Meijde et al., 2003; Huang et al., 2005; Bolfan-Casanova, 2005); (ii) phase changes that occur at these depths, A c c e p t e d M a n u s c r i p t like transformation of pyroxenes into garnet, or olivine successively into wadsleyte and into ringwoodite; (iii) the extreme softening of a material as it undergoes a phase transition, known as transformational superplasticity (Sammis & Dein, 1974).

  have shown that transition zone minerals at average mantle temperatures have anomalously high water solubility compared to upper and lower mantle minerals, suggesting that the transition zone might act as a water reservoir. The potential presence of water in the transition zone, also revealed by the seismological analysis of Van der Meijde et al. (2003), could also explain the discrepancy between the velocity jump at 410 km observed seismically and the one expected for an olivine-rich (pyrolite) mantle (Duffy & Anderson, 1989; Anderson, 1989; Anderson & Bass, 1986). Smyth & Jacobsen (2006) proposed that lateral velocity variations in the transition zone may reflect variations in hydration rather that variations in temperature.

( 2 Page 17 of 38 A

 238 ii) The possible role of the dilution of pyroxenes into garnets (the major phase change at transition zone depths) could be enlightened by the knowledge of the creep laws for these two minerals. To date, available experimental data are sparse, but uniaxial compression and hot hardness tests[START_REF] Karato | Plastic deformation of garnets: systematics and implications for the rheology of the mantle transition zone[END_REF] demonstrated that the resistance to plastic deformation in garnets is significantly higher than most of the other minerals in the Earth's mantle; the pyroxene-garnet phase change, then, cannot explain the low viscosity in the transition zone. An alternative explanation could reside in the transformation of olivine into wadsleyte and then ringwoodite (Artem Oganov, personal communication, 2007).(iii) The third possible explanation for the soft transition zone might be attributed to a phenomenon known as transformational superplasticity, first pointed out by[START_REF] Sammis | On the possibility of transformational superplasticity in the earth's. mantle[END_REF], that consists in a dramatic reduction in effective viscosity observed during a phase transition in materials like metals and ceramics[START_REF] Poirier | Transformation plasticity, in Creep of Crystals[END_REF][START_REF] Maehara | Superplasticity in ceramics[END_REF][START_REF] Meike | A critical review of investigations into transformation plasticity[END_REF]. Panasyuk & Hager (1998) tested a model of transformational superplasticity for the upper mantle and estimated the degree of softening for mantle material at the phase change at 400 km depth: the viscosity decrease would be of 1-c c e p t e d M a n u s c r i p t orders of magnitude, consistent with what we found.

  s c r i p t[96] Wu, P., 1992. Deformation of an incompressible viscoelastic flat earth with power-law creep: a finite element approach, Geophys. J. Int., 108, 136-142.

Fig. 1 .

 1 Fig. 1. Example of the performance of the EA evolving from generation 1 to 500 (a). Blue and red dots represent the variance reduction (%) of the best-fitting model and the average fitness of the population, for each generation. The decrease in mean variance reduction at 10-20 generations corresponds to an adjustment of the mutation-rate parameter defined in Section 2.2. Fit of best model as a function of generation number (b) after 100 generations.

  A c c e p t e d M a n u s c r i p t

Fig. 2 .

 2 Fig. 2. After choosing velocity-to-density scaling relationship and tomographic model, we seek the mantle viscosity profile corresponding to the best fit of GRACE gravity data. Density anomalies and gravity data are related through sensitivity kernels, whose form in turn depends on the viscosity profile.

Fig. 5 .

 5 Fig. 5. Best-fitting mantle viscosity models resulting from inversions of GRACE free-air gravity anomalies and v S models S20RTS with different number of layers n.Only relative variations can be inferred from these models.

Fig. 6 .

 6 Fig.6. Best-fitting mantle viscosity models resulting from inversions of GRACE free-air gravity anomalies and v S model S20RTS, with same number of layers and different depths of viscosity discontinuities. Only relative variations can be inferred from these models. The fitness of gravity anomalies computed from these viscosity profiles to the data is displayed above each panel.

29

 29

  A c c e p t e d M a n u s c r i p t

Fig. 7 .

 7 Fig. 7. Best-fitting viscosity profile (a) obtained on the basis of some v S models, assuming n = 5 and converting v S anomalies to ρ anomalies via the scaling factors in (b). Frames(c-d-e) show average viscosity (red/green/blue lines) and standard deviation (gray intervals) of models with fit better than a given threshold.

Fig. 8 .

 8 Fig. 8. Mantle viscosity profile obtained on the basis of SMEAN v S model, assuming n = 5 and converting v S anomalies to ρ anomalies via the scaling factors in Figure (7b), red line. The harmonics are summed up to degree l = 16. The figure shows average viscosity (in red) and standard deviation (gray intervals) of models with fit better than a given threshold.

  A c c e p t e d M a n u s c r i p t

Fig. 9 .

 9 Fig. 9. Viscosity profile (b) obtained on the basis of SMEAN v S model, assuming n = 5 and converting v S anomalies to ρ anomalies via the scaling factor in (a), blue line (Karato and Karki, 2001). Frame (c) shows average viscosity derived using different velocity-to-sensity scaling factors for suboceanic mantle (red line) and subcontinental mantle (green line). Gray intervals represent the standard deviation of models with fit better than a given threshold.

  A c c e p t e d M a n u s c r i p t

Fig. 10 .

 10 Fig.10. Range of variability of the viscosity of the 5 layers (normalized to the value of viscosity of the shallowest layer) obtained inverting gravity data with density structure from v S model S20RTS (red curve Figure7).

Fig. 12 .

 12 Fig. 12. Best-fitting viscosity profiles based on v S models (panel a) TRP246 (solid line) and SPRD6 (dashed line) and ρ models (panel b) TRP246 (solid line) and SPRD6 (dashed line).

Fig. 14 .

 14 Fig. 14. Viscosity profile (b) obtained on the basis of the ρ model by Simmons et al. (2007). The red line represents the average viscosity, gray intervals correspond the standard deviation of models with fit better than a given threshold.

Fig. 15 .

 15 Fig. 15. CMB dynamic topography based on v S models S20RTS, TRP246, SPRD6 (frames a, b, c, respectively); and on ρ models TRP246 and SPRD6 (frames d, e). The scale for each map is ±5 km (a, b), ±8 km (c), ±15 km (d), ±12 km (e). Blue colors indicate elevation, red colors indicate depression.

  The mantle flow theory we employ takes into account many of the complexities of the real Earth, like sphericity, compressibility and self-

	gravitation, and additionally allows to compute surface dynamic topography,
	CMB deflections and plate motions via the surface divergence.
	This approach to modeling the surface gravity field has, however, several lim-
	itations: first, as seen before, only the ratios between viscosity values at dif-
	ferent depths can be constrained, rather than the absolute values of viscosity,
	thus neglecting the effects of toroidal flow and associated lateral viscosity vari-
	ations. The problem is also complicated by the existence of many, competing
	tomographic models of seismic velocity. The amplitude and pattern of seismic
	M velocities in the mantle are known only approximately, and despite the agree-ment at long wavelengths (Becker & Boschi, 2002), various tomographic im-ages differ in shape, depth extent, and amplitude of fine features (e.g., Becker & Boschi, 2002; Romanowicz, 2003; Boschi et al., 2007). Last, establishing an appropriate velocity-to-density scaling for the mantle is not straightforward. Growing evidence suggests that seismic velocity anomalies reflect both ther-mal and compositional heterogeneities (van der Hilst & Karason, 1999; Karato & Karki, 2001; Deschamps et al., 2001). Unlike that of temperature variations, the effect of compositional variations on seismic velocities and density is not yet well understood. A solution to this problem is to use a density model derived from a seismic-geodynamic inversion which implicitly includes both a n u s c r i p t thermal and compositional effects on buoyancy (Simmons et al., 2007).
	A c c e p t e d

  Some effects of lateral viscosity variations on geoid and surface velocities induced by density anomalies in the mantle, Geophys. J. Int., 114, 531-547.

	Model	Type of anomaly Dynamic ellipticity (km)
	S0RTS	δv S	1.5
	TRP246	δv S	1.8
	SPRD6	δv S	1.8
	TRP246	δρ	4.8
	SPRD6	δρ	1.5
	Table		M a n u s c r i p t
	A c c e p t e d

Acknowledgments. This research was partially supported by the NERIES program. We are grateful to Adam Dziewonski and Giorgio Spada for fruitful discussions, to Dave Yuen for continuing encouragement, and to Miaki Ishii, Jeroen Ritsema and Jeannot Trampert for providing their tomographic models. We thank Rob Moucha and Mark Jelinek for useful comments that improved a first version of this paper. The implementation of the EA technique used in this analysis is made via the PIKAIA Fortran routine, developed by Charbonneau and Knapp (1995) at the High Altitude Observatory in Boulder, Colorado, and available electronically on the Observatory ftp archive. Figures were prepared with GMT by Wessel and Smith (1991).