28 Aug-1st Sept., 2006.

de Viron, O., I. Panet, and M. Diament (2006).
Extracting low frequency climate signal
from GRACE data, *e-Earth*, **1**, 9-14.

Duquenne H. (2007a). Modelling the vertical gravity gradient for gravity measurements reduction. Proc. of the 1st International Symposium of the IGFS. Forthcoming.

Duquenne H. (2007b). A data set to test geoid computation methods. Proc. of the 1st International Symposium of the IGFS. Forthcoming.

Francis O., van Dam T., Amalvict M., Andrade da Sousa M., Bilker M., Billson R., d'Agostino G., Desogus, Falk R., Gremak A., Gitlein O., Jonhson D., Klopping F., Kostelecky J., Luck B., Mäkinen J., McLaughlin D., Nunez E., Origlia C., Palinkas V., Richard P., Rodriguez E., Ruess D., Schmerge D., Thies S., Timmen L., van Camp M., van Westrum D., Wilmes H. (2005). Results of the International Comparison of Absolute Gravimeters in Walferdange (Luxembourg) of November 2003, IAG Symposia, vol. 129, Springer.

Froger, J-L., Remy, D., Bonvalot, S., Legrand, D. (2006). Dynamic of long term multi-scale inflations at Lastarria-Cordon del Azufre volcanic complex, central Andes, revealed from ASAR-ENVISAT interferometric data. *Earth and Planetary Science Letters* (doi: 10.1016/j.epsl. 2006.12.012).

Gabalda, G., Bonvalot, S., Hipkin, R. (2003).
CG3TOOL : An interactive computer

program for Scintrex CG-3M gravity
data processing. *Computers &
Geosciences*, 29, 2, pp. 155-171.

Gattacceca T. (2005). Campages de mesures gravimétriques. Internal report IGN/SGN CR/G 193, Institut Géographique National.

Gébelin A., Martelet G., Brunel M., Faure M. and Rossi Ph. (2004). Late Hercynian leucogranites modelling as deduced from new gravity data: the example of the Millevaches massif (Massif central, France), *Bull. Soc. Géol. Fr.*, 175, 239-248.

Gébelin A., Martelet G., Chen Y., Brunel M., and Faure M. (2006). Structure of late Variscan Millevaches leucogranite massif in the French Massif Central: AMS and gravity modelling results, *Journal of Structural Geology*, 28, 148-169.

Greff-Lefftz, M., L. Metivier, and H. Legros (2005). Analytical solutions of love numbers for a hydrostatic ellipsoidal incompressible homogeneous Earth, Celestial Mechanics & Dynamical Astronomy, 93 (1-4), 113-146.

Guillen A., Joly.A., Calcagno.P., Courrioux.G., Fitzgerald.D., McInerney.P., Seikel.R. (2006). Coupling 3D modelling and forward - Inverse modelling of potential field data (gravity and magnetic data and tensor components). Oral presentation, IAMG'2006.

Gunawan, H. (2005). Etudes gravimétrique et microgravimétrique appliquées à la
volcanologie. Exemples de la Soufrière
de Guadeloupe (Antilles) et du Merapi
(Java). Ph. Thesis, IPGP.

Hetenyi, G., R. Cattin, J. Vergne, and J.L. Nabelek (2006). The effective elastic thickness of the India Plate from receiver function imaging, gravity anomalies and thermomechanical modelling, *Geophys. J. Int.*, **167** (3), 1106-1118.

Hinderer, J., Sedighi, M., Bayer R., Ghazavi, K., Amalvict, M., Luck, B., Nilforoushan, F. Masson F., Peyret, M., Djamour, Y. and Kouhzare, A. (2003). The absolute gravity network in Iran: an opportunity to analyse gravity changes caused by present-day tectonic deformation, in Proc. IMG-2002 (Instrumentation and Metrology in Gravimetry) *Cahiers du Centre Européen de Géodynamique et de Séismologie*, Luxembourg, vol. 22, pp. 137-141.

Ihde J., Baker T., Bruyninx C., Francis O., Amalvict M., Kenyeres A., Makinen J., Shipman S., Simek J. and Wilmes H. (2005). First Results of the Development of a European Combined Geodetic Network (ECGN), IAG Symposia, vol. 129, Springer.

Jacob T., Bayer R., Le Moigne N., Maia M., Moreau F., Bour O., Davy P., Hinderer J., Luck B., Lalancette M.F., Batany C. (2006). Compte-rendu des Mesures Gravimétriques à Ploemeur, Présentation orale, 3-15 Mai 2006,

Colloque du GDR-G2, La Rochelle, 22-24 novembre 2006

Jamet, O., and M. Diament, On the estimation of local spectral representations of gravity anomalies, Proc. Analytical Representation of Potential Field Anomalies for Europe (AROPA), *Cahiers du Centre Européen de Géodynamique et de Séismologie*, 20, 71-77, 2003.

Joly A., Chen Y., Faure M., Martelet G. and Truffert C. (2006). A multidisciplinary study on the Carboniferous (Namurian) Montmarault granitic pluton and its tectonic bearings on the activity of the Sillon Houiller Fault (Variscan Massif Central, France). Oral presentation, EGU General Assembly, Vienna, Austria.

Lequentrec-Lalancette M.F., Rouxel D. & Lenamouric C. (2006). Géoïde européen: contribution du SHOM, *Report n° 449 EPSHOM/CH/GG/NP*.

Llubes, M., Florsch N., N., Hinderer, J., Longuevergne, L., Amalvict, M. (2004). Local hydrology, Global Geodynamics Project and CHAMP / GRACE / GOCE perspective: Some case studies, *Journal of Geodynamics*, **38**, 355-374.

Llubes M. and the loading team (2006). A multitechnique monitoring of ocean loading in North of France » submitted to C.R.A.S.

Louis P. (ed.) (1996). Bilan et prospective en gravimétrie. CNFGG.

Mäkinen J., Amalvict M., Shibuya K., Fukuda Y. (2007). Absolute gravimetry in
Antarctica: Status and prospects,
*Journal of Geodynamics*, in press

Martelet G., Calcagno P., Gumiaux C., Truffert C., Bitri A., Gapais D. and Brun J.P. (2004). Integrated 3D geophysical modelling of the Hercynian Suture Zone in the Champtoceaux area (south Brittany, France), *Tectono-physics*, 382, 117-128.

Martelet G. and Debeglia N. (2006). Banque gravimétrique de la France : travaux réalisés en 2006, tâches à réaliser en 2007, BRGM/RP-55009-FR.

Martelet G., and Sarrailh M. (2006) Gravity data control for the determination of an accurate regional geoid, Ext. Abstr., in 1st symposium of IGFS, 28/08-1/09/2006, Istambul, Turkey.

Masson F. (2006) Active tectonics studied by absolute gravimetry: the AGTO project in Taiwan. Oral presentation, Taiwan, November 2006

Mazabraud Y., N. Béthoux and S. Deroussi (2005). “Characterisation of the seismological pattern in a slowly deforming intraplate region: central and western France.” *Tectonophysics* 409(1-4): 175-192.doi:10.1016/j.tecto.2005.08.021.

Métivier, L. (2004). Influence des variations latérales de densité et de paramètres rhéologiques sur la déformation de la Terre. Ph. Thesis, IPGP.

Métivier, L., M. Greff-Lefftz, and M. Diament (2005). A new approach to computing accurate gravity time variations for a realistic Earth model with lateral

heterogeneities, Geophys. J. Int. (doi: 10.1111/j.1365-246X.2005.02692.x).

Métivier, L., M. Greff-Lefftz, and M. Diament (2006). Mantle lateral variations and elasto-gravitational deformations: (1) numerical modeling., Geophys. J. Int., 167, 1060-1076.

Mikhailov, V., A. Galdeano, M. Diament, A. Gvishiani, S. Agayan, S. Bogoutdinov, E. Graeva, and P. Sailhac (2003). Application of artificial intelligence for Euler solutions clustering, *Geophysics*, 68 (1), 168-180.

Mikhailov, V., S. Tikhotsky, M. Diament, I. Panet, and V. Ballu (2004). Can tectonic processes be recovered from new gravity satellite data ?, *Earth Planet. Sci. Lett.*, 228, 281-297.

Mikhailov, V.O., and M. Diament (2006). Some aspects of interpretation of tensor gradiometry data, *Izvestiya, Physics of the Solid Earth*, 42 (12), 971-978.

Moreaux, G; M.-F. Lequentrec-Lalancette, D. Rouxel, F. Ponchaut (2004) On the Parameterization of Planar Models. Oral presentation, *Geoid Gravity and Space Missions Meeting, Sept 2004*, Porto.

Moreaux, G; Lequentrec-Lalancette MF, Rouxel D. (2005). Marine geoid from gravity and gradiometry. Oral presentation, *ECGS/GOCINA Workshop*, April 13-15, Luxembourg.

Nicolas J., Nocquet J.-M., Van Camp M., Boy J.P., Hinderer J., Gegout P., Calais E., and Amalvict M. (2006). Seasonal effect on Laser, GPS, and Absolute
Gravimetry vertical positioning at the
OCA geodetic station, Grasse, France,
*Geophys. J. Int.*, 10.1111/j.1365-246X.
2006.03205.x

Pajot, G., Mikhailov, V., LeQuentrec-Lalancette, M.F. & Diament, M., Gravity gradiometry: enhanced methods for enhanced data. Oral presentation, EGU 2006, Vienna, Austria, 2006.

Panet, I., O. Jamet, M. Diament, and A. Chambodut (2004). Modelling the Earth's gravity field using wavelet frame, in *IAG Symposia*, vol. 129, pp. 48-53, Springer.

Panet, I. (2005). Les ondelettes sphériques en gravimétrie spatiale. Applications en Polynésie Française et à l'étude du séisme de Sumatra. Ph. Thesis, IPGP, 2005.

Panet, I., A. Chambodut, M. Diament, M. Holschneider, and O. Jamet (2006). New insights on intraplate volcanism in French Polynesia from wavelet analysis of Grace, Champ and sea-surface data., *J. Geophys. Res.*, 111, doi:10.1029/2005JB004141.

Pavez, A. (2005). Structure et déformations du volcan LASCAR à partir d'observations par satellites et au sol : apport à la connaissance et la surveillance des volcans andésitiques. Ph. Thesis, IPGP.

Pavez, A., Remy, D., Bonvalot, S., Diament, M., Gabalda, G., Froger, J.L., Legrand, D., Julien, P., Moisset, D. (2006). Insight into ground deformations at Lascar Volcano (Chile) from InSAR, Photogrammetry and GPS data:

implications on volcano dynamics and
future space monitoring. *Remote
Sensing of Environment*, 100, 3, 307-
320.

Petit, C., Déverchère, J. (2006). Structure and evolution of the Baikal rift, a synthesis, G3, Vol 7, 11.

Peyrat, S., DeChabalier, J-B., Nercessian, A., Bouin, M-P., Campos, J., Legrand, D., Bonvalot, S., Vilotte, J-P. Source process of the June 13, 2005 Tarapaca earthquake (Chile) (2006a) *American Geophysical Union (AGU) Fall meeting, San Francisco*. Geophysical Research Abstracts, Vol. 8, 04951, Dec. 2006.

Peyrat, S., J. Campos, J.B. de Chabalier, A. Perez, S. Bonvalot, M.-P. Bouin, D. Legrand, A. Nercessian, O. Charade, G. Patau, E. Clevede, E. Kausel, P. Bernard, J.-P. Vilotte (2006b). The Tarapaca intermediate-depth earthquake (Mw 7.7, 2005, Northern Chile): A slab-pull event with horizontal fault plane constrained from seismologic and geodetic observations. *Geophysical Research Letters*, V33, LL22308, doi:10.1029/2006GL027710.

Rémy, D. (2005). Analyse et inversion de séries interférométriques et micro-gravimétriques temporelles sur les volcans actifs : apport à la quantification d'effets de sites et à la compréhension de la dynamique volcanique. Ph. Thesis, IPGP.

Sato T., Okuno J.,. Hinderer J., MacMillan DS, Plag HP, Francis O., Falk R. and Fukuda Y. (2006). A geophysical
interpretation of the secular
displacement and gravity rates observed
at Ny-Ålesund, Svalbard in the Arctic—
effects of post-glacial rebound and
present-day ice melting, *Geophys. J. Int.*
165, 729–743.

Savanier, D., R. C. Maury, G. Guille, C. Legendre, P. Rossi, H. Guillou, S. Blais, S. Deroussi (2005). Carte géol. France (1/50 000), feuille de Nuku Hiva, Polynésie française, notice explicative. Orléans: BRGM.

Sepúlveda S., Lahsen A., Bonvalot, S., Cembrano J., Alvarado A., Letelier P. (2005). Morpho-structural evolution of the Cordón Caulle geothermal region, Southern Volcanic Zone, Chile: insights from gravity and 40Ar/39Ar dating. *Journal of Volcanology and Geothermal Research*, Volume 148, Issues 1-2, 15 October 2005, pp 165-189

Talbot, J.Y. (2003). Apport des études ASM et gravimétriques des plutons cévenols à la caractérisation structurale de l’évolution tardi-hercynienne du Massif Central. Ph. Thesis, Université d’Orléans.

Talbot J.Y., Martelet G., Courrioux G., Chen Y. et Faure M. (2004). Emplacement in an extensional setting of the Mont Lozère-Borne granitic complex (SE France) inferred from comprehensive AMS, structural and gravity studies, *Journal of Structural Geology*, 26, p. 11-28.

Talbot J.Y., Faure M., Chen Y. and Martelet G. (2005). Pull-apart emplacement of the Margeride granitic complex (French

Massif Central). Implications for the
late evolution of the Variscan orogen,
*Journal of Structural Geology*, 27,
1610-1629.

Tiberi, C., M. Diament, J. Déverchère, C. Petit-Mariani, V. Mikhailov, S. Tikhotsky, and U. Achauer (2003). Deep structure of the Baikal rift zone revealed by joint inversion of gravity and seismology, *J. Geophys. Res.*, 108 (B2), 2109, doi:10.1029/2002JB001880.

Tiberi, C., Ebinger, C., Ballu, V., Stuart, G., Oluma, B. (2005). Inverse model of gravity data from the Red Sea-Aden-East African rifts triple junction zone, *Geophysical Journal International*, 163, doi: 10.1111/j.1365-246X.2005.02736.x

Tirel, C., Gueydan, F., Tiberi, C., Brun, J.-P. (2004). Aegean crustal thickness inferred from gravity inversion. Geodynamical implications, *Earth and Planetary Science Letters*, 228, 267–280.

Tiwari, V.M., M. Diament, and S.C. Singh (2003). Analysis of satellite gravity and bathymetry data over Ninety-East Ridge: variation in the compensation mechanism and implication for emplacement process., *J. Geophys. Res.*, 108 (B2), 2109, doi: 10.1029/2000JB000047.

Verdun, J., E.E. Klingele, R. Bayer, M. Cocard, A. Geiger, and H.G. Kahle (2003). The alpine Swiss-French airborne gravity survey, *Geophysical Journal International*, 152 (1), 8-19.

Verdun, J., and E.E. Klingele (2005). Airborne
gravimetry using a strapped-down
LaCoste and Romberg air/sea gravity
meter system: a feasibility study,
*Geophysical Prospecting*, **53** (1), 91-
101.

Vitushkin, L., M. Becker, Z. Jiang, O. Francis, T. M. van Dam, J. Faller, J.-M. Chartier, M. Amalvict, S. Bonvalot, N. Debeglia, S. Desogus, M. Diament, F. Dupont, R. Falk, G. Gabalda, C. G. L. Gagnon, T. Gattacceca, A. Germak, J. Hinderer, O. Jamet, G. Jeffries, R. Käker, A. Kopaev, J. Liard, A. Lindau, L. Longuevergne, B. Luck, E. N. Maderal, J. Mäkinen, B. Meurers, S. Mizushima, J. Mrlina, D. Newell, C. Origlia, E. R. Pujol, A. Reinhold, Ph. Richard, I. A. Robinson, D. Ruess, S. Thies, M. Van Camp, M. Van Ruymbeke, M. F. de Villalta

Compagni, S. Williams (2003). Results of the Sixth International Comparison of Absolute Gravimeters ICAG-2001. *Metrologia*, 39, 5, pp. 407-424.

Widiwijayanti, C., V. Mikhailov, M. Diament, C. Deplus, R. Louat, S. Tikhotsky, and A. Gvishiani (2003). Structure and evolution of the Molucca Sea area: constraints based on interpretation of a combined sea-surface and satellite gravity dataset, *Earth Planet. Sci. Lett.*, 215, 135-150.

Widiwijayanti, C., C. Tiberi, C. Deplus, M. Diament, V. Mikhailov, and R. Louat (2004). Geodynamic evolution of the northern Molucca Sea area (Eastern Indonesia) constrained by 3-D gravity field inversion, *Tectonophysics*, 386, 203-222.
# Time-Varying Gravity with Superconducting Gravimeters

J. HINDERER, J.-P. BOY, S. ROSAT, Y. ROGISTER,

M. AMALVICT, P. GEGOUT & C. DE LINAGE

IPGS (UMR 7516 CNRS/ULP), 5, rue Descartes 67084 Strasbourg Cedex, France
Jacques.Hinderer@eost.u-strasbg.fr

**Abstract** This paper is devoted to the superconducting gravimeter (SG) and its broad spectral range of applications in Earth Sciences. We present the recent scientific achievements using SGs from the Global Geodynamics Project (GGP). We show how the extremely low noise in the long-period seismic band helps to retrieve weak seismic normal modes. Special attention is paid to the search for the Slichter translational mode of the solid inner core. Atmospheric loading in gravity is investigated, as well as solid Earth tidal effects including the Free Core Nutation resonance in the diurnal band. Ocean loading can also be investigated through the detection of non-linear ocean tides and the observation of long-period tides. The signature in SG data of non-tidal ocean circulation, continental hydrology and Earth rotation effects all demonstrate the ability of these instruments to study long period geodynamical processes. Some specific applications related to tectonics (post-glacial rebound, coseismic effects) are shown. We end on the possibility of using surface gravity measurements with a regional network of SGs (like in Europe) to validate satellite-derived gravity missions such as GRACE.

**Key words** superconducting gravimeter, GGP, time-varying gravity, normal modes, Slichter triplet, Earth tides, ocean tides, atmospheric pressure, hydrology, polar motion, tectonics, GRACE, satellite gravimetry

## Variations Temporelles du Champ de Pesanteur à l'Aide de Gravimètres Supraconducteurs

**Résumé** Cet article est consacré aux gravimètres à supraconductivité et leur large utilisation en Sciences de la Terre. Nous présentons les résultats scientifiques récents
les plus marquants obtenus à partir des gravimètres supraconducteurs du réseau international GGP (Global Geodynamics Project). Le niveau de bruit très faible dans la bande sismique longue période permet de restituer des modes propres d'amplitude faible. Nous nous intéressons aussi à la recherche du mode de Slichter de translation de la graine solide. Les effets de surcharge atmosphérique et de marée solide incluant la résonance du noyau liquide dans la bande diurne sont décrits. La surcharge océanique peut être étudiée par l'intermédiaire des ondes de marée non-linéaires et des ondes à longue période. La signature dans les données des gravimètres supraconducteurs des effets de circulation océanique, d'hydrologie continentale et de rotation de la Terre démontre la capacité de ces instruments à restituer les phénomènes géodynamiques de longue période. Quelques applications tectoniques (rebond post-glaciaire, effets co-sismiques) sont considérées. Nous terminons sur la possibilité d'utiliser les données de variations de pesanteur au sol d'un réseau régional de gravimètres supraconducteurs (exemple de l'Europe) afin de valider les missions de gravimétie spatiale de type GRACE.

***Mots clefs*** gravimètre supraconducteur, GGP, variation temporelle de gravité, modes propres, triplet de Slichter, marées terrestres et océaniques, pression atmosphérique, hydrologie, rotation, tectonique, GRACE, gravimétrie spatiale

# I. INTRODUCTION

The Global Geodynamics Project (GGP) (http://www.eas.slu.edu/GGP/ggphome.htm l) is an international research effort that was launched in 1991. In 2003 it became an Inter-Commission Project of IAG (International Association of Geodesy) and reports to Commission 2 (The Gravity Field) and Commission 3 (Earth Rotation and Geodynamics). The GGP consists of a worldwide network of superconducting gravimeters (SGs), currently about 25 instruments. Vertical time-varying gravity acceleration data in a standard form, basically untouched raw data decimated to 1

minute samples, are sent at the end of every month to a database (http://ggp.gfz-potsdam.de). Phase I of the project was the period (1997-2003) and we are currently in phase II (2003-2007). The geographical distribution of the GGP stations is shown in Fig. 1. Current operational stations are in yellow, temporarily stopped stations are in red and the planned stations (most of them to be installed in 2007) are in green.

The scientific objectives of the GGP cover geophysical phenomena throughout the wide period range of applications of the instruments (from a few seconds to several
years), covering topics such as normal modes, mantle rheology, tides, solid Earth-oceans-atmosphere interactions, hydrology, and Earth rotation. There are recent review papers on the scientific achievements using

SGs (Hinderer and Crossley 2004; Meurers 2004; Crossley & Hinderer 2005). In the following, we present the recent results achieved in the various fields of SG applications.

<figure><img src="image_65.png" /><figcaption>**Figure 1:** Geographical locations of the current and planned GGP stations. In yellow: current operational stations; in red: temporarily stopped sites; in green: planned stations (most of them to be installed in 2007).</figcaption></figure>

## II. NORMAL MODES

Investigation with SGs of the seismic normal modes excited by large earthquakes has led to some new and impressive results because of the very low noise of these instruments (Widmer-Schniedrig 2003) below 1 mHz compared to the long-period seismometers; see also Van Camp et al. (2005) for a noise comparison between AG

and SG. Of recent large earthquakes, the 2001 Peru event of magnitude $M_w = 8.4$ and the $M_w = 9.3$ Sumatra huge event on December 26, 2004 strongly excited the long period seismic modes, and significant observations were made using SG data of the GGP network. In particular, the fully split $^0S_2$ multiplet (with 5 individual
singlets) that has been rarely visible on a
single instrument, could be fully analyzed
from the Strasbourg C026 instrument (Rosat
et al. 2003a) (see Fig. 2) but was also
present at other SG stations. The most
important new result was the detection of
the overtone $_2S_1$ elastic mode never
observed before. Another example of the
strong signal to noise ratio in the modes
generated by the Sumatra quake can be
found in (Rosat et al. 2005) and in (Roult et
al. 2006). The frequency and damping
analysis of these low-frequency seismic
modes put some new constraints on the
density structure inside the Earth.

A recent investigation on SG data
has also shown the fundamental interest of
using the amplitudes and phases of low
frequency vertical components of normal

modes to retrieve the source parameters of
large earthquakes such as rupture length and
duration (Lambotte et al. 2006).

One major goal of GGP was to
detect the translational motion of the solid
inner core (the Slichter triplet $_1S_1$) the
period of which is determined primarily by
a gravitational restoring force; therefore its
period(s) would bring a new constraint on
the density contrast at the inner core
boundary (ICB) (Rogister 2003; Rosat et al.
2006) and possibly also on viscosity just
above the ICB. There is obviously no clear
observational evidence in the stacking of
SG data of any of the theoretically predicted
triplets (Rosat et al. 2003b and 2004).
Precise atmospheric and hydrologic
corrections are required in the search for
this translational motion.

<figure><img src="image_66.png" /><figcaption>**Figure 2:** Splitting of the fundamental spheroidal mode $_0S_2$ after the June 2001 Peru earthquake of magnitude M=8.4.</figcaption></figure>
<figure><img src="image_67.png" /><figcaption>**Figure 3:** Bayesian inversion in superconducting gravimetry (up) and in VLBI vertical displacement (down) (Florsch et al., 2006).</figcaption></figure>

### III. ATMOSPHERIC AND HYDROLOGIC EFFECTS

In order to have the most accurate atmospheric corrections for the new gravity satellite missions, CHAMP and GRACE, 2D pressure loading computations have been extended to 3D modeling where pressure, temperature and humidity data varying with height are also taken into account in addition to the surface data (see

Boy and Chao 2005; Neumeyer et al. 2004a). It was recently shown for the GGP station Medicina (Italy) with the help of balloon radio sounding that there is a seasonal vertical air mass change in the atmosphere without ground pressure changes which is caused by warming and cooling that leads by attraction to a non-
negligible gravity effect of the order of one
µGal.

Early studies in gravity focused mainly on tides, oceanic and atmospheric effects but, in the last decade, more and more attention has been paid to the hydrology signature in gravity, especially in SG data where the drift is not a large issue. Most of the studies are restricted to modeling local effects by trying to find correlations between gravity residuals and a relevant hydrological parameter such as water table level, rainfall or soil moisture (see Harnisch & Harnisch 2006a). Active experiments, where known amounts of water were added to specific areas in the gravimeter vicinity, were conducted at Moxa Observatory and helped significantly in the validation and parameterization of hydrology models (Kroner and Jahr 2006). More recently a study was devoted to the seasonal changes in SGs (Boy and Hinderer 2006). These gravity changes can be linked to global hydrology models. In fact, for

more than half of the 20 analysed SG, there is a good correlation. For the other stations, the discrepancies may be associated with local hydrology effects, especially when a station is partly underground like Moxa or Strasbourg. Very recently a paper by Van Camp et al. (2006) investigated in details the hydrogeological contributions to gravity for the station Membach in Belgium. Finally several papers (Romagnoli et al. 2003; Richter et al. 2004) carefully examined both height and gravity changes from continuous GPS and SG (C023) observations at the station Medicina (Italy). These studies provide a convincing interpretation of the seasonal signal from the combined loading contributions of air pressure, ocean circulation and hydrology (surficial water table) illustrating once more the need to have simultaneous height and gravity measurements in addition to the monitoring of various environmental parameters in the close vicinity of the station.

## IV. EARTH AND OCEAN TIDES

Earth and ocean tides are by far the largest components of surface gravity changes. The SGs have brought two areas of improvement to tidal studies. First, the high sensitivity of these meters which enables them to detect small amplitude tidal signals previously hidden in noise (e.g. non-linear

ocean tides) and to retrieve with better precision larger tidal signals. With this high precision, Xu et al. (2004a) revisited the question of the possible latitude dependence of tidal gravimetric factors. Using 19 GGP stations they found that the discrepancy of the four principal waves ($O_1$ and $K_1$ in the
diurnal band, $M_2$ and $S_2$ in the semi-diurnal band) between observations and theoretical models is less than 0.2 %. This means there is no significant latitude dependence. Second, the much lower instrumental drift of SGs versus mechanical spring meters permits more precise studies of long period tides ($M_f, M_m, S_{Sa}, S_a$) (Ducarme et al. 2004; Boy et al. 2006).

The largest tidal components occur in both semi-diurnal and diurnal frequency bands. Within the diurnal band, some waves are affected by a resonance that occurs due to the Nearly Diurnal Free Wobble (NDFW), also called Free Core Nutation (FCN), which is a differential rotation of the fluid outer core with respect to the mantle. We do not discuss here the numerous papers using GGP data to derive the FCN parameters and refer the reader to previous reviews on this subject (Hinderer and Crossley 2004). A more global analysis by Xu et al. (2004) who used simultaneously tidal gravity observations at 19 GGP stations leads to a FCN period of 429.9 (427.2 - 432.7) sidereal days and a $Q$ value about 20 000 (12 000-72 000). The FCN period is in good agreement with space geodetic studies, but the gravity-derived $Q$ factor is often much smaller than that obtained using VLBI. An artifact of the analysis is sometimes even a negative $Q$ (Sun et al. 2004a) due to an inappropriate application of the least squares method

usually used for the determination of the resonance parameters (obviously a linearized form applied to a non-linear problem). The classical least squares method is inadequate because it implies Gaussian statistics that are not correct for this problem.

An appropriate Bayesian method (Florsch and Hinderer, 2000) is required to solve for the FCN parameters due to the non-standard form of the probability distribution for $Q$. Fig. 3 shows the Bayesian approach applied to the Strasbourg SG data by plotting the joint probability distributions for the eigenperiod $T$, quality factor $Q$ and resonance strength $a$ (Florsch et al. 2004, 2006). The mismodelling of the ocean tides that are corrected in SG data before determining the resonance parameters can lead to underestimated $Q$ values (Florsch and Hinderer, 2000).

The impact of the ocean and atmospheric loading contributions in gravity is now better understood (Boy et al. 2003, 2004).

Observations from the GGP network can be compared to computations from different global ocean tide models, some of them being partly hydrodynamic, or assimilating data derived directly from satellite altimetry (Baker and Bos, 2003, Boy et al., 2003 and Bos and Baker, 2005). The present situation, from the
comprehensive treatment by Bos and Baker (2005) is that the class of recent ocean tide models, specifically GOT00 (Ray 1999), NAO99b (Matsumoto et al. 2000) and FES99 (Lefèvre et al. 2002), is in better agreement with SG observations than earlier models such as the classical Schwiderski's model proposed in 1980. When differences are carefully examined most models are in agreement and there is no 'best' model valid in all areas of the world.

A corollary to this investigation is the limits we can place on Earth's inelasticity from the gravimetric amplitude factors and phase delays after correction for ocean tidal loading. This was pointed out by Baker and Bos (2003) where they compared for some waves the gravimetric amplitude factor and phase delay to some reference elastic or slightly inelastic values.

We have pointed out that the SGs have a very small instrumental drift compared to classical spring meters and this is why the investigation of long period tides is particularly suited to SG data. In a recent paper, Boy et al. (2006) analyzed long series from 18 GGP stations to estimate the ocean loading for the monthly ($M_m$) and fortnightly ($M_f$) tides. The available models were an equilibrium tide and three recent hydrodynamical models with satellite altimeter data assimilation NAO99b (Matsumoto et al. 2000), FES99 (Lefèvre et al. 2002) and TPXO.6 (Egbert and

Erofeeva, 2002). They concluded that the uncertainty for $M_m$ is still too large to discriminate between the newest models. On the other hand the hydrodynamical models for $M_f$ are clearly closer to the SG observations than the equilibrium tidal model or the older model proposed by Schwiderski (1980).

Non-linear tides can also be seen in SG records. At CA (Cantley) in Canada it is possible to see clearly these small abnormal waves originating from non-linearities in the ocean tidal response in the Bay of Fundy (Merriam, 1995). The sensitivity inferred from the SG measurements is equivalent to 1 mm of open ocean tide, which means that these observations are a unique tool for validating these small ocean tidal waves that exist only near the coastlines.

Some years ago, small signals of a few nGal were identified in the SG residuals with periods around 6 hours from Strasbourg and Cantley stations. A more systematic study by Boy et al. (2004) on various GGP stations demonstrated that these signals are indeed due to non-linear tides in the quadiurnal frequency band. Due to the integrative properties of the gravity loading Green's functions, inland SG observations act as a complementary large scale validation tool to point-like tide gauge observations of non-linear ocean tides. The latter are generated only in shallow water
near the coasts and hence mostly escape
detection by satellite altimeters like

TOPEX/POSEIDON or JASON.

<figure><img src="image_68.png" /><figcaption>**Figure 4:** GGP/GRACE gravity comparison in Europe.</figcaption></figure>

## V. NON-TIDAL OCEAN EFFECTS

In addition to tidal oceanic contributions, non-tidal effects related to the general oceanic circulation are also detectable in SG measurements. Virtanen and Makinen (2003) investigated the loading effect of the Baltic Sea on the Metsähovi (Finland) instrument (T020) which is located only 15

km away from the open sea. They found a good correlation between SG residuals and sea level changes from the nearby Helsinki tide gauge. Zerbini et al. (2004) investigated also non-tidal ocean effects both in gravity and height observations at the Medicina station in Italy.

## VI. EARTH ROTATION

At the long period end of the observable spectrum of gravity, there are two isolated

theoretical periods of the Earth's normal mode spectrum. One is the inner core
wobble (ICW), whose period is of the order
of years and which is undetected in gravity.
The other mode is much more accessible,
and is of course the Chandler Wobble, with
a period of about 435 days in the mantle
reference frame. It is important to use a long
data set when trying to infer the transfer
function of the polar motion; clearly more
than 6.5 years of continuous data are
required to separate the annual component
(whatever its origin) from the 435 day term.
Harnisch and Harnisch (2006b) used data
from 12 GGP stations with lengths varying
between 4 and 18 years to investigate the
polar motion contribution to gravity. They
found that in general the gravimetric
amplitude factors for the Chandler wobble

are close to the nominal value of 1.16 and
phase lags of a few degrees. Xu et al.
(2004b) also made a similar study on 5 SGs
from the GGP network. Combining the
results leads to a δ factor of 1.16 ± 0.07 and
to a phase delay of -1.30° ± 1.33°. Finally a
data set of 9 SGs was analyzed for the polar
motion response by Ducarme et al. (2006)
who found a mean δ factor of 1.179 ±0.004.

A value slightly larger than the
nominal value of 1.16 is to be expected
when including mantle inelasticity and/or
ocean pole tide contribution. The variability
on the phase delays between different
stations is however still large and
unexplainable in terms of Earth's
rheological properties.

VII. TECTONICS

There is another class of long-term gravity
changes which is not periodic, but secular
and primarily due to active tectonics or
post-glacial rebound (PGR). Until recently,
most studies have been done using AGs
rather than SGs. One reason is that the SGs
have a small instrument drift that must be
accounted for as part of any estimate of the
secular gravity trend. The potential of long
term gravity measurements to solve tectonic
problems is however significant, as clearly
demonstrated by a recent study (Lambert et
al. 2006) devoted to PGR in North America.
In particular, they could show that the

admittance of gravity to height changes (in
µGal/mm) deduced from combined GPS
and AG measurements is indeed close to the
theoretically predicted values. In fact,
combining gravity and height changes
measurements is very efficient in
discriminating PGR effects from present-
day ice melting contributions as shown by a
recent investigation of the long term gravity
changes in Svalbard (Sato et al. 2006).
Another study based on 7 years of
collocated gravity measurements of SG and
AG at Membach (Belgium) by Francis et al.
(2004) indicates that there is a small
decrease in gravity connected to uplift seen
by GPS at the same location.

Such comparisons of AG and GPS
measurements are much easier to do when
there are continuous SG measurements at
the same location. In fact the continuity of
the SG record is often important to check
the integrity of the AG measurements and to
model effects such as hydrology that
increase the AG scatter. Amalvict et al.
(2004, 2006) analyzed the long term gravity
changes of the Strasbourg SG with regular
AG measurements and collocated GPS and
interpreted the results in terms of hydrology
and tectonics of the Rhine graben.

Several recent SG results point the
way to future studies involving gravity. One
of these is the study by Imanishi et al.
(2004) of coseismic displacements that
accompanied a large Mw = 8.0 earthquake

in Japan. Using 3 SGs in a linear array they
showed that gravity offsets occurred that
matched the theoretical prediction of static
displacements that accompany a fault
rupture process. The significance of this
result can be transferred to another context,
that of subduction induced silent
earthquakes that have been identified off the
coast of Japan and also with the Cascadia
subduction zone. Episodic tremor that
accompanies the silent slip events can last
for several days, but there is no identifiable
earthquake. Such events are identified
primarily as horizontal displacements on
GPS arrays, but Lambert et al. (2006) have
shown that AG measurements also reveal
coincident offsets in gravity of several
µGal, in between which there is a secular
drift (in this case a negative trend).

VIII. GROUND/ORBIT VALIDATION

There is a major international effort in the
present decade to measure variations in the
Earth's global gravity field using low orbit
satellites. The first satellite CHAMP
(Challenging Mini-satellite Payload) was
launched in 2000 and was followed two
years later by GRACE (Gravity Recovery
and Climate Experiment). In the near future,
there will be a third mission called GOCE
(Gravity field and steady-state Ocean
Circulation Explorer) that will orbit even

closer to the ground and hence be even
more sensitive to smaller scale gravity
changes. The primary goal of these missions
is to use the temporal changes of the Earth's
gravity field to infer changes in regional and
continental water storage, and ocean
circulation.

A major concern with satellite
measurements of time varying gravity is
how to calibrate and validate such
observations. In addition to comparisons
with models (the primary technique used to date), there are several possibilities using actual measurements at the ground (GPS, gravity) or at the ocean bottom (water pressure). This problem is important because if successful, it would produce an independent method of validation that does not rely solely on modeling. We will show hereafter that surface measurements from the GGP network provide a useful additional constraint on space gravity data. The validation signal is related to seasonal gravity changes that are coherent on length scales appropriate to satellite altitudes (typically a few hundred of km).

A first study directly comparing CHAMP data to 6 SG ground observations was done by Neumeyer et al. (2004b) and has led to satisfactory results for all the stations in the one year analysis period. Neumeyer et al. (2006) recently extended this study to GRACE data pointing out again the partial agreement between surface and satellite-derived gravity at specific locations. Before a detailed comparison can be made, however, one has to remember that ground gravity measurements include necessarily a contribution from the vertical motion of the instrument through the ambient gravity field. This signal does not affect the orbiting satellite and hence there is a difference in the gravity changes as seen at (moving) ground level and by the satellite (Hinderer et al. 2006a).

We note in these studies that the comparison of single station results with the large-scale satellite solutions is problematic due to the completely different error budgets involved. GRACE data for example are good to 1 µGal only over length scales longer than 500-1000 km, whereas SGs are good to the same accuracy (or better) at a single point. In order to average SG measurements and reduce local effects, there have been attempts to assemble a network solution from nearby SG stations rather than doing the above single station comparison. Within the existing rather sparse GGP network, Europe is obviously the best place to try such an approach.

The approach was first initiated using one year of SG data and extended to longer data sets by Crossley et al. (2003; 2004, 2005). This approach was further extended to a 21 month time interval to inter-compare surface data (GGP European sub-network), satellite data from GRACE, and theoretical predictions for two hydrology models LaD (Milly and Shmakin 2002) and GLDAS (Rodell et al. 2004) (Andersen et al. 2005a; Hinderer et al. 2006a).

The results show the existence of an annual signal that is coherent over Europe with a few µGal amplitude mostly due to the seasonal loading from continental hydrology (soil moisture + snow) according to recent global models. There is even a
possibility to detect in GRACE data
interannual signals (Andersen and Hinderer
2005) and, in particular, there is a clear
evidence that GRACE has been affected by
the heat wave that occurred in summer 2003
in Europe (Andersen et al. 2005b). As
indicated in Fig. 4, there is a decrease in the
soil water content of about 6 cm in summer
2003 with respect to summer 2002. The
Wettzell (Germany) and the Medicina

(Italy) SG ground data seem to confirm this
point.

A proposal is underway to combine
an SG with AG and portable Scintrex-type
gravimeters to calibrate/validate GRACE
satellite measurements in West Africa
(Hinderer et al. 2005; 2006b). One profile
covers the extremely dry regions in Central
Sahara desert down to the coast of Benin in
the monsoon region of very high rainfall.

IX. CONCLUSION

The necessity to precisely correct time-
varying gravity data from the atmospheric
and hydrologic effects has been highlighted,
particularly for the comparison between
surface and satellite-derived gravity
variations, and also for the search of weak
signals associated with the internal
dynamics of the Earth (seismic normal
modes and Slichter mode). The ocean
loading corrections still need to be
improved in the diurnal tidal band if we
want to improve the retrieval of the free
core nutation parameters.

We have also emphasized the use of
precise tidal determinations with SGs for
validating the newest ocean tide models
(linear and non-linear tides). In the long
period part which can be investigated
thanks to the very low instrumental drift of
SGs, a number of interesting applications in
hydrology, Earth rotation or tectonics were
shown. Finally the most recent progress
achieved in the ground validation of the
GRACE mission with a regional network of
SGs has been pointed out.

X. REFERENCES

Amalvict, M., Hinderer, J., Makinen, J., Rosat,
S., and Rogister, Y., 2004. Long-term
and seasonal gravity changes and their
relation to crustal deformation and
hydrology, J. Geodynamics, 38, 343-353.

Amalvict, M., Hinderer, J., and Rozsa, S., 2006.

Crustal vertical motion along a profile
crossing the Rhine graben from the
Vosges to the Black Forest Mountains:
Results from absolute gravity, GPS and
levelling observations, J. Geodynamics,
41, 358-368
Andersen, O., and Hinderer, J., 2005. Global inter-annual gravity changes from GRACE: early results, Geophys. Res. Lett., 32, L01402, doi: 10.1029/2004GL020948.

Andersen, O., Hinderer, J., and Lemoine, F. G., 2005a. Seasonal Gravity Field Variations from GRACE and Hydrological Models, in (eds) Jekeli, Bastos and Fernandes, Gravity, Geoid and Space Missions. IAG symposia vol. 129, 316-321, Springer Verlag.

Andersen O., Seneviratne, S., Hinderer, J., and Viterbo, P., 2005b. GRACE-derived terrestrial water storage depletion associated with the 2003 European heat wave, Geophys. Res. Lett., 32, L18405, doi: 10.1029/2005GL023574.

Baker, T., and Bos, M., 2003. Validating Earth and ocean tide models using tidal gravity measurements, Geophys. J. Int., 152, 468-485.

Bos, M. S., and Baker, T. F., 2005. An estimate of the errors in gravity ocean tide loading computations, J. Geod. 79, 50-63.

Boy, J.-P., and Chao, B. F., 2005. Precise evaluation of atmospheric loading effects on Earth's time-variable gravity field, Journal of Geophysical Research, 110, B08412, doi: 10.1029/2002JB002333.

Boy, J.-P., and Hinderer, J., 2006. Study of the seasonal gravity signal in superconducting gravimeter data, Journal of Geodynamics, 41, 227-233.

Boy, J.-P., Llubes, M., Hinderer, J., and Florsch, N., 2003. A comparison of tidal ocean loading models using superconducting gravimeter data, J. Geophys. Res., 108, B4, 2193, doi: 10.1029/2002JB002050.

Boy, J.-P., Llubes, M., Ray, R., Hinderer, J., Florsch, N., Rosat, S., Lyard, F., and Letellier, T., 2004. Non-linear oceanic tides observed by superconducting gravimeters in Europe, J. Geodynamics, 38, 391-405.

Boy, J.-P., Llubes, M., Ray, R., Hinderer, J., and Florsch, N., 2006. Validation of long-period oceanic tidal models with superconducting gravimeters, Journal of Geodynamics, 41, 112-118.

Crossley, D. and Hinderer, J., 2005. Using SG arrays for hydrology in comparison with GRACE satellite data, with extension to seismic and volcanic hazards. Korean J. Remote Sensing, 21 (1) 31-49.

Crossley, D., Hinderer, J., Llubes, M., and Florsch, N., 2003. The potential of ground gravity measurements to validate GRACE data, Advances in Geosciences, 1, 1-7.

Crossley, D., Hinderer, J., & Boy, J.-P., 2004. Regional gravity variations in Europe from superconducting gravimeters, J. Geodynamics, 38, 325-342.

Crossley, D., Hinderer, J., and Boy, J.-P., 2005. Time variation of the European gravity field from superconducting gravimeters, Geophys. J. Int., 161, 257-264.

Ducarme, B., Venedikov, A., Arnoso, J., and Vieira, R., 2004. Determination of the long period tidal waves in the GGP superconducting gravity data, Journal of Geodynamics, 38, 307-324.
Ducarme, B., Venedikov, A., Arnoso, J., Chen, X.-D., Sun, H.-P., and Vieira, R., 2006. Global analysis of the GGP superconducting gravimeters network for the estimation of the pole tide gravimetric amplitude factor, Journal of Geodynamics, 41, 334-344.

Florsch, N. and Hinderer, J., 2000. Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data, Phys. Earth Planet Int., 117, 21-35.

Egbert, G. D., and S. Y., Erofeeva, 2002. Efficient inverse modeling of barotropic ocean tides, Journal of Atmospheric and Oceanic Technology, 19 (2), 183-204.

Florsch, N., Hinderer, J., Petrov, L., Boy, J.-P., & Llubes, M., 2004. *Combinaison bayésienne de données de gravimétrie et de déformations VLBI pour la détermination des paramètres de la résonance du noyau*, Colloque GDR G2, Le Mans, novembre 2004.

Florsch, N., Llubes, M., Petrov, L., Boy, J.-P., & Hinderer, J., 2006. Combined gravimetric and VLBI free core nutation parameter Bayesian inversion, Combined Workshop of the special study group on 'Analysis of Environmental Data for the Interpretation of Gravity Measurements', the Global Geodynamics Project (GGP), and of the new working group on 'High Precision of Tidal Predictions', Jena, Germany, March 2006.

Francis, O., Van Camp, M., van Dam, T., Warnant, R., and Hendrickx, M., 2004. Indication of the uplift of the Ardenne in

long-term gravity variations in Membach (Belgium), Geophys. J. Int., 158, 346-352.

Harnisch, G., & Harnisch, M., 2006a. Hydrological influences in long gravimetric data series, Journal of Geodynamics 41, 1-3, 276-287.

Harnisch, M., & Harnisch, G., 2006b. Study of long-term gravity variations, based on data of the GGP co-operation, Journal of Geodynamics, 41, 1-3, 318-325.

Hinderer, J., and Crossley, D., 2004. Scientific achievements from the first phase (1997-2003) of the Global Geodynamics Project using a worldwide network of superconducting gravimeters, J. Geodynamics, 38, 237-262.

Hinderer, J., de Linage, C., and Boy, J.-P., 2005. Issues in the ground validation of satellite-derived gravity measurements: a proposal to validate GRACE in Africa from the Sahara to the equatorial monsoon, EOS Trans. AGU, 85 (52), Fall Meet. Suppl., Abstract G33A-0028.

Hinderer, J., Andersen, O., Lemoine, F., Crossley, D., and Boy, J.-P., 2006a. Seasonal changes in the European gravity field from GRACE: A comparison with superconducting gravimeters and hydrology model predictions, Journal of Geodynamics, 41, 59-68.

Hinderer, J., de Linage, C., & Boy, J.-P., 2006b. How to validate satellite-derived gravity observations with gravimeters at the ground?, Bull. Inf. Marées Terr., 142, 11433-11441.

Imanishi, Y., Sato, T., Higashi, T., Sun, W., and
Okubo, S., 2004. A network of superconducting gravimeters detects submicrogal coseismic gravity changes, Science, 306, 476-478.

Kroner, C., and Jahr, T., 2006. Hydrological experiments around the superconducting gravimeter at Moxa Observatory, Journal of Geodynamics, 41, 268-275.

Lambert, A., Courtier, N., and James, T.S., 2006. Long-term monitoring by absolute gravimetry: Tides to postglacial rebound, Journal of Geodynamics, 41, 307-317.

Lambotte, S., L. Rivera, and J. Hinderer, 2006. Rupture length and duration of the 2004 Aceh-Sumatra earthquake from the phases of the Earth's gravest free oscillations, Geophys. Res. Lett., 33, L03307, doi: 10.1029/2005GL024090.

Lefèvre, F., F. H. Lyard, C. Le Provost and E. J. O. Schrama, 2002. FES99: A global tide finite element solution assimilating tide gauge and altimetric information, Journal of Atmospheric and Oceanic Technology, 19, 1345-1356.

Matsumoto, K., T. Takanezawa and M. Ooe, 2000. Ocean tide models developed by assimilating TOPEX/Poseidon altimeter data into hydrodynamical model: A global model and a regional model around Japan, Journal of Oceanography, 56, 567-581.

Merriam, J., 1995. Non-linear tides observed with the superconducting gravimeter at Cantley, Quebec, Geophys. J. Int., 123 (2), 529-540.

Meurers, B., 2004. Investigation of temporal gravity variations in SG-records, Journal of Geodynamics 38, 423-435.

Milly, P. C. D., and A. B. Shmakin, 2002. Global modeling of land water and energy balances. Part I: The land dynamics (LaD) model, Journal of Hydrometeorology, 3 (3), 283-299.

Neumeyer, J., Hagedoorn, J., Leitloff, and Schmidt, R., 2004a. Gravity reduction with three-dimensional atmospheric pressure data for precise ground gravity measurements, J. Geodyn., 38, 437-450.

Neumeyer, J., Schwintzer, P., Barthelmes, F., Dierks, O., Imanishi, Y., Kroner, C., Meurers, B., Sun, H.-P., and Virtanen, H., 2004b. Comparison of superconducting gravimeter and CHAMP satellite derived temporal gravity variations. In: Reigber, Ch., Lühr, H., Schwintzer, P., Wickert, J. (Eds.), Earth Observations with CHAMP Results from Three Years in Orbit, pp. 31-36.

Neumeyer, J., Barthelmes, F., Dierks, O., Flechtner, F., Harnisch, M., Harnisch, G., Hinderer, J., Imanishi, Y., Kroner, C., Meurers, B., Petrovic, S., Reigber, C., Schmidt, R., Schwintzer, P., Sun, H.-P., & Virtanen, H., 2006. Combination of temporal gravity variations resulting from Superconducting Gravimeter (SG) recordings, GRACE satellite observations and global hydrology models, J. of Geodesy, doi: 10.1007/s00190-005-0014-8.

Ray, R., 1999. A global ocean tide model from TOPEX/Poseidon altimetry: GOT99.2, NASA Technical Memorandum, 209478, 58 pp.
Richter, B., Zerbini, S., Matonti, F., and Simon, D., 2004. Long-term crustal deformation monitored by gravity and space techniques at Medicina, Italy and Wettzell, Germany, Journal of Geodynamics, 38, 281–292.

Rodell, M., P. R. Houser, U. Jambor, J. Gottschalck, K. Mitchell, C.-J. Meng, K. Arsenault, B. Cosgrove, J. Radakovich, M. Bosilovich, J. K. Entin, J. P. Walker, D. Lohmann and D. Toll, 2004. The Global Land Data Assimilation System, Bulletin of the American Meteorological Society, 85 (3), 381-394.

Rogister Y. (2003): Splitting of seismic free oscillations and of the Slichter triplet using the normal mode theory of a rotating, ellipsoidal Earth, Phys. Earth planet. Int., 140, 169-182.

Romagnoli, C., Zerbini, S., Lago, L., Richter, B., Simon, D., Domenichini, F., Elmi, C., and Ghirotti, M., 2003. Influence of soil consolidation and thermal expansion effects on height and gravity variations, J. of Geodynamics, 521-539.

Rosat, S., Hinderer, J. & Rivera, L., 2003a. First observation of $_{2}$S$_{1}$ and study of the splitting of the football mode $_{0}$S$_{2}$ after the June 2001 Peru earthquake of magnitude 8.4, Geophys. Res. Lett., 30, 21, 2111.

Rosat, S., Hinderer, J., Crossley, D. & Rivera, L., 2003b. The search for the Slichter Mode: Comparison of Noise Levels of Superconducting Gravimeters and Investigation of a Stacking Method, Phys. Earth planet. Int., 140, 183-202.

Rosat, S., Hinderer, J., Crossley, D., and Boy, J.-P., 2004. Performance of superconducting gravimeters from long-period seismology to tides, J. Geo-dynamics, 38, 461-476.

Rosat, S., Sato, T., Imanishi, Y., Hinderer, J., Tamura, Y., McQueen, H., & Ohashi, M., 2005. High-resolution analysis of the gravest seismic normal modes after the 2004 M$_{w}$ = 9 Sumatra earthquake using superconducting gravimeter data, Geophys. Res. Lett., 32, L13304, doi: 10.1029/2005GL023128.

Rosat S., Rogister Y., Crossley D. & Hinderer J., 2006. A search for the Slichter triplet with superconducting gravimeters: Impact of the density jump at the inner core boundary, Journal of Geodynamics, 41, 296-306.

Roult, G., Rosat, S., Millot-Langet, R., Clevede, E., & Hinderer, J., 2006. First determination of *Q* quality factors and new measurements of frequencies for the whole set of singlets of the Earth's normal modes $_{0}$S$_{0}$ , $_{0}$S$_{2}$ , $_{0}$S$_{3}$ and $_{2}$S$_{1}$ using Superconducting Gravimeter data from the GGP network, J of Geodynamics, 41, 345-357.

Sato, T., Okuno, J., Hinderer, J., MacMillan, D. S., Plag, H.-P., Francis, O., Falk, R. and Fukuda, Y., 2006. A geophysical interpretation of the secular displacement and gravity rates observed at Ny-Ålesund, Svalbard in the Arctic- Effects of the post-glacial rebound and present-day ice melting, Geophys. J. Int., 165, 729-743.

Schwiderski, E. W., 1980. On charting global
ocean tides, Reviews of Geophysics and
Space Physics, 18 (1), 243-268.

Sun, H.-P., Jentzsch, G., Xu, J.-Q., Hsu, H.-Z.,
Chen, X.-D., and Zhou, J.-C., 2004.
Earth's free core nutation determined
using C032 superconducting gravimeter
at station Wuhan/China, Journal of
Geodynamics 38, 451-460.

Van Camp, M., Williams, D. P., Francis, O., 2005. Uncertainty of absolute gravity measurements, Journal of Geophysical Research, Vol. 110, B05406, doi: 10.1029/2004JB003497.

Van Camp, M., M. Vanblooster, O. Crommen, T. Petermans, K. Verbeeck, B. Meurers, T. van Dam and A. Dassargues, 2006. Hydrogeological investigations at the Membach station, Belgium, and application to correct long periodic gravity variations, J. Geophys. Res., 111, B10403, doi: 10.1029/2006JB004405.

Virtanen, H., and Mäkinen, J., 2003. The
effect of the Baltic Sea level on gravity at

the Metsähovi station. Journal of
Geodynamics, 35, 4-5, 553-565.

Widmer-Schnidrig, R., 2003. What can superconducting gravimeters contribute to normal mode seismology, Bull. Seis. Soc. Am., 93 (3), 1370-1380.

Xu, J.-Q., Sun, H.-P., and Ducarme, B., 2004a. A global experimental model for gravity tides of the Earth, Journal of Geodynamics, 38, 293-306.

Xu, J.-Q., Sun, H.-P., and Yang, X.-F., 2004b. A study of gravity variations caused by polar motion using superconducting gravimeter data from the GGP network, Journal of Geodesy, 78: 201-209, doi: 10.1007/s00190-004-0386-1.

Zerbini, S., F. Matonti, F. Raicich, B. Richter and T. van Dam, 2004. Observing and assessing nontidal ocean loading using ocean, continuous GPS and gravity data in the Adriatic area, Geophys. Res. Lett., 31, L23609, doi: 10.1029/2004GL021185.
SISMOLOGIE ET PHYSIQUE
DE L'INTÉRIEUR DE LA TERRE
# Ionospheric Seismology

P. LOGNONNÉ¹, G. OCCHIPINTI¹, A. KHERANI¹, P.E. GODET¹,
J. ARTRU², R. GARCIA³, T. FARGES⁴, H. HÉBERT⁴, V. BAZIN⁵,
F. CRESPON⁶, G. MOREAUX⁶, J. HELBERT⁶ & E. JEANSOU⁶.

¹Institut de Physique du Globe de Paris, 4 Avenue de Neptune, 94100 Saint Maur des Fossés
lognonne@ipgp.jussieu.fr

²Centre national d'Etudes Spatiales, Toulouse, France

³Dynamique Terrestre et Planétaire, Observatoire Midi-Pyrénées, Toulouse, France

⁴Commissariat à l'Energie Atomique/DASE, Bruyères le Chatel, France

⁵Office National d'Etudes et de Recherches Aérospatiales/DEMR, Palaiseau, France

⁶NOVELTIS, Toulouse, France

**Abstract** The sounding of the ionosphere with GPS, Doppler sounder or Radar allows the detection of acoustic waves or gravity waves generated by quakes or tsunami. These waves are induced by the vertical displacement of the Earth surface. They propagate almost vertical, with an amplification with altitude associated to the exponential decay of the atmospheric density. For ground displacement of a few mm, amplitudes of a few tens to hundred meters are achieved at 250-300 km of altitude. The first observations of these waves in the ionosphere were done in the 60th. The technological development in the last 10 years has lowered the detection threshold and movies can now be done of the signals, allowing seeing the horizontal propagation of the perturbations.

Even if this new way in seismology is making the first steps, several perspectives seem to be possible: the first one is related to the tsunami detection. The tsunami signal appears about 15 minutes after the uplift of the see-level and is completely developed about 30 minutes after. It is therefore possible, with optimized radars, like over-the-horizon radar, to perform a regional monitoring of the ionosphere. A second perspective is related to the precise imaging of the wave front of the signals in the ionosphere. This might detect not only the surface waves, but also the body waves generated by large quakes. This might be an interesting alternative to the deployment
of dense oceanic bottom seismic network, even if this will detect only the largest quakes.

**Key Words** Seismology, ionosphere, remote-sensing, tsunami

## Sismologie de l'ionosphère

**Résumé** Le sondage de l'ionosphère, par l'intermédiaire de réseaux GPS, de sondeurs Doppler ou de Radar, permet de détecter les ondes acoustiques ou les ondes de gravités générées par les séismes ou tsunami. Typiquement, ces ondes sont générées par les déplacements verticaux de la surface terrestre. Elles se propagent alors quasi-verticalement, en étant amplifiées avec l'altitude, en raison de la décroissance de la densité de l'atmosphère. Pour des déplacements du sol de l'ordre du mm, des amplitude de quelques dizaines à centaines de mètre sont donc atteinte vers 250-300km d'altitude. D'une manière générale, les premières observations de ces ondes dans l'ionosphère remonte aux années 1960-1970. Les développements technologiques des dix dernières années permettent cependant d'abaisser les seuils de détection et de réaliser des véritables films qui permettent de suivre la propagation horizontale de ces ondes dans l'ionosphère.

Même si cette nouvelle approche en sismologie est encore à ses premiers pas, plusieurs perspectives semblent déjà se dégager: la première est associé à la détection et au suivi des tsunami. Le signal ionosphérique apparait en effet 15 minutes après l'initiation du tsunami et est complètement développé 30 minutes environ. Il est donc possible, avec des radars adaptés, comme les radar trans-horizon, de réaliser une surveillance régionale des zones à risques. Une seconde perspective a trait à l'imagerie fine des fronts d'ondes des ondes de surface, voir des ondes de volume. Ceci permettrait ainsi d'avoir une alternative au déploiement de réseaux sismiques denses, particulièrement difficile pour ce qui est des réseaux fond de mer.

**Mots clefs** Seismologie, ionosphère, télédétection, tsunami

# I. INTRODUCTION

Recent developments in ionosphere remote sensing, in particular techniques using the Global Positioning System (GPS), Doppler HF sounder or even over-the-horizon radar

provide unprecedented capabilities for monitoring the reaction of the ionosphere to seismic waves. Three types of signals can be addressed by a routine survey and
monitoring of the ionosphere. The first
ones are produced by an acoustic coupling
between the solid Earth and the ionosphere
near the source and lead to post-seismic
acoustic signals. The second are related to
the generation of acoustic waves
associated to surface waves and even body
waves, and can be detected remotely far
from the source. For quakes with moment

magnitude larger than about 7, the
detection is worldwide. Finally, the third
are related to coupling of oceanic gravity
waves (i.e. tsunami) with atmospheric and
ionospheric gravity waves.

We present briefly in this paper the
state of the art in the detection of these
signals by several groups in France and
address the perspectives.

<figure><img src="image_69.png" /></figure>

**Figure 1:** Seismic surface waves after the Mw = 7.6 Chi-Chi earthquake (Taiwan, September 20, 1999) as measured on a ground seismometer (bottom panel) at the Geoscope station SSB (Saint-Sauveur, France) and on the CEA ionospheric Doppler sounding network (Francourville, France), corresponding to the vertical motion of ionospheric layers at altitudes 168 and 186 km. All traces show the vertical velocity perturbation in the 1-50 mHz frequency band. An amplification of 4.10⁴ is observed between the ground and the ionosphere. The ~8 minutes delay between the ground and the ionosphere corresponds to the propagation time of the acoustic wave. Due to this delay, body waves are expected to arrive in the ionosphere at 18h17, 18h23 for PP and S waves respectively, while surface waves arrive at about 18h40. SV waves, due to SV-P conversion, are detected.
## II. THEORETICAL AND HISTORICAL BACKGROUND

After an earthquake, seismic waves or tsunami generate surface motions at the Earth's surface. By continuity of the vertical displacement, the atmosphere is forced to move with a vertical velocity equal to the surface vertical velocity, and the induced perturbation propagates upward as an atmospheric wave. The propagation is done primarily by acoustic processes, for signals with frequencies higher than the Brünt-Vaïssalla frequency, and by gravity processes, for frequencies above.

These atmospheric vibrations produce adiabatic pressure and temperature variations in the neutral atmosphere. At higher altitude, the velocity of the neutral atmosphere can be transferred to ions, through collision processes, and forces the later to oscillate. Due to the neutrality of the ionosphere, these ions perturbations generate electron density perturbations. Both the electrons and ions in the processes are sensitive to the Lorentz forces related to the magnetic and electric fields.

In the 1960s, the first published observations were performed with ionospheric sounders for $M_s=8$ quakes [1-4]. The development of Doppler sounders in the last thirteen years has lead to a reduction of the detection threshold, as

we will see later.

The Global Positioning System (GPS), with the development of the associated ground networks, has opened new possibilities in the detection and imaging of these signals. As a by-product of geodetic measurements, phase and code measurements from GPS stations can indeed be used to calculate the electron density of the ionosphere between the receiver and the satellite. The primary data provide the integrated value of the electron density, or Total Electronic Content (TEC) [5]. After the first detection following the 1994 Northridge earthquake [6], seismic signals have been detected for many other quakes and have been interpreted as the result of a shock wave generated by the supersonic Rayleigh wave [7].

A detailed theory has been developed [8] to take into account the coupling between the solid Earth, the ocean and the atmosphere. In the later, the boundary conditions of the elasto-dynamic operator at the solid Earth - atmosphere interface is integrated in the normal modes theory. A radiative boundary condition simulates the escape of acoustic and gravity atmospheric waves in the upper ionosphere, where no refraction of waves is observed. This theory allows the computation of normal modes, and those can be used to compute
not only seismograms for spherical
models, but also barograms and neutral
velocities. The dissipation related to
viscosity was later incorporated [9]. This
technique provides the neutral density

perturbation and velocities of the upper
atmosphere. The later acts as neutral wind
and generates ionospheric perturbations
[10].

III. TOWARD REMOTE SENSING SEISMOLOGY?

During more than 4 decades, the detected
signals described above were however
more or less considered as some “funny”
or “exotic” observation in seismology,
unable to provide new valuable
informations, either on the source or on the
internal structure of the Earth. However,
we are now facing, with the development
of new technologies in ionospheric
sounding, or with the dense GPS networks
progressive changes, which put a new light
on these researches and start to point out
possible seismological interests and
applications.

Following the pioneering works
done with analog Doppler sounder, new
generations of sounders have been
developed such as the Doppler sounder
operated by CEA/DASE in France and
working continuously since August
1999[11]. The French instrument is
detecting most of the earthquakes with Ms
greater than 6.5 [12]. Doppler data are
very similar to seismograms in the sense
that they measure directly the vertical
motion of a ionospheric layer.

The performances of Doppler
sounder are moreover such that both the
surface waves and body waves are
detected in the ionosphere, including SV
waves (see Fig 1). However they remain
limited to a small number of point
measurements and cannot resolve the 3D
structure of the perturbation. New studies
are therefore done in order to use Over-
The-Horizon radars, which might provide
maps of the ionospheric vertical
displacements [13]. Preliminary results
show that the signal to noise ratio of these
instruments is probably comparable to
those obtained by Doppler sounders and
that these instruments could therefore be a
way to provide dense measurements of the
seismic wavefront, with sampling as low
as 1 measurement per 25 km2 over surface
of several 106 km2.

A next step was performed
following the 2002 Denali earthquake
[14,15]. The dense California GPS
networks was used to map the ionospheric
perturbations and to compute the group
velocity with a high spatial resolution
above the Northwestern US coast, with a
precision compatible with the
identification of the signature associated to
lateral variation in the Earth's lithosphere.
The 3D structure of the seismic
ionospheric signal was then characterized
[16] and with such approach, the
comparison of signals from identical

altitude can be performed. Figure 2a shows
a vertical cut of the signal detected above
Hokkaido, after the Hokkaido Tokacho-
Oki earthquake of September 25, 2003
while Figure 2b shows signals with
increasing epicentral distances, all
corresponding to a 300 km altitude.

<figure><img src="image_70.png" /></figure>

**Figure 2a:** Vertical cut of the 3D Rayleigh waves impact in the ionosphere for the Tokacho-Oki event. The Total Electronic Content amplitudes observed are typically 0.1 TECU peak-to-peak but 3D local variations reach a few $10^9$ e/m$^3$. No wavefront is observed with a north or northwest propagation direction, due to a poor coverage of the GPS satellite in these directions.

These results were obtained thank
to a collaboration with Japan GSI. In the
near field, these signals are associated to
the near-field seismic waves and especially
to infrasonic waves generated near the
source and propagating in the atmospheric

wave-guide. It has been shown that the
ionospheric total electronic content
perturbations are sensitive to the focal
mechanism and can therefore be used for
inverting the source dynamic and
geometry [17]. Note that 2D ionospheric
maps over Europe are now available for either seismic or other applications, in the framework of an academic-SME

collaboration with the Toulouse based SME Noveltis [18]. See also [19] for a review on GPS seismic observations.

<figure><img src="image_71.png" /></figure>

**Figure 2b:** TEC ionospheric signals at about 300 km of altitude after 3D tomography reconstruction. Note in the far field the signal associated to the Rayleigh waves (dashed line) and in the near field those associated to the acoustic pulse.

## IV. IONOSPHERIC TSUNAMIS WARNING

As for surface waves, early theoretical works in the 1970s predicted that atmospheric gravity waves are generated in the wake of a tsunami [20]. About 30 minutes are needed for the gravity wave to develop its first maximum perturbation in the ionosphere (versus ~10 minutes for seismic-acoustic waves). But after this delay the ionospheric perturbation follows the tsunami front and, as for the seismic waves, the atmospheric oscillations are amplified with altitude. It should be noted moreover that, due to their much shorter

wavelength and period, the surface noise of ocean swell does not produce significant upward propagating waves in the atmosphere: the atmosphere acts as a filter, enhancing the long wavelength tsunami perturbation over other sources. Figure 3 shows the result of simulation, where the tsunami first generates an atmospheric gravity wave which is then generating, through collisions between neutral atmosphere and ions, perturbations in the electronic density.

The first observation had however
to wait almost 30 years. It was performed
after the Peru, June, 2001 tsunami [21].
The tsunami arrival was observed on
Japanese tide gauges between 20 and 22
hours after the earthquake, with wave
amplitudes between 10 and 40 cm (open
ocean amplitude were estimated to be of 1-
2 cm) and dominant periods of 20 to 30
minutes. Shortly after, a large ionospheric

perturbation was detected through a
specific processing of data from the
continuous GPS network in Japan
(GEONET). The arrival time, orientation,
wavelength, velocity of the wave packet
observed are consistent with what is
expected for a tsunami-induced
perturbation.

<figure><img src="image_72.png" /></figure>

**Figure 3:** Coupling between the neutral atmosphere gravity wave induced by a tsunami and the ionosphere. The tsunami amplitude has a 0.5 meter amplitude and about 13min period, corresponding to the amplitude of the Sumatra, 2004 tsunami. From top to below are the normalised neutral wind, and the absolute and relative electron density. This shows that large perturbations, reaching 10%, are generated by such tsunamis.

The gigantic and dramatic Sumatra
tsunami of December, 2004 confirmed the
possibilities of observing tsunami
ionospheric signals, and signals were
detected on the Total Electronic Content
(TEC) measurement on-board the

TOPEX/Poseidon and JASON satellites,
and on the GPS stations in Indonesia and
in the India Ocean [19]. A few GPS
stations were quite close from the
epicenter and have shown that the
ionospheric perturbation starts about 15
minutes after the event and reach a first maximum less than 30 minutes after the signal was performed, and both the waveform and the amplitude observed by Jason and Topex has been reproduced [22]. These results confirm the interest of a

real-time monitoring of the ionosphere, which could be carried out either with active microwave radar or by optical systems detecting the airglow associated with the ion recombination in the ionosphere.

<figure><img src="image_73.png" /><figcaption>**Figure 4:** Long period vertical atmospheric oscillations, for a $10^{18}$ Nm quake ($M_w=5.9$) and for period larger than 100 sec on Venus. Due to the difference in the acoustic coupling at the ground, ionospheric signals at 150 km of altitude are about 100 stronger on Venus for the same magnitude and altitude than on Earth.</figcaption></figure>

# V. EXPORTING REMOTE SENSING SEISMOLOGY ON VENUS?

Although on Earth this technique would never provide the same quality of seismic data as a seismic network, it can be the unique way to obtain seismic data on planets too hostile for the deployment of long lived seismic stations. Venus is the

best example [23,24]. In addition, the coupling strength is proportional to the acoustic impedance of the atmosphere, equal to $\rho_c$ where $\rho$ is the density and $c$ the acoustic speed. As the atmospheric density at the surface of Venus is about 60 kg/m$^3$
and the acoustic velocity is slightly higher (410 m/s) than on Earth, this leads to an acoustic impedance about 60 times greater than on Earth, where the atmospheric density is 1.2 kg/m³.

Moreover, at 50 km of altitude, where the Venus pressure is comparable to Earth ground pressure, the decrease by almost 2 order of magnitude of the density

leads already to an amplification of 10 of the waves. Consequently, Venus quakes will generate atmospheric infrasonic waves with amplitudes much larger than on the Earth surface [Figure 4]. This profitable effect give an unique opportunity for a future Venus quakes detection by a satellite sounding the Venus ionosphere.

## VI. CONCLUSION

Advance in the monitoring of small-scale perturbations of the ionosphere have allowed the detection of atmospheric Rayleigh waves as well as tsunami-induced gravity waves with both ground systems based on GPS, and ionospheric sounding performed by TOPEX and JASON and Doppler sounders. Doppler sounders are also sensitive to body waves. These new data open exciting prospects in seismology such as the remote sensing of the Rayleigh seismic wave fronts, especially over the ocean, where the deployment of dense seismic networks is the most challenging. These techniques might also provide in a future a high resolution picture of the wave front of body waves. These prospects are also very

exciting for tsunamis because they are extremely difficult to observe in the open ocean, but their associated gravity waves have a clear impact on the ionosphere and can be detected by remote sensing systems. The monitoring of the ionosphere by joint ground/space techniques, such as continuous GPS networks, over-the-horizon radar or even by a future dedicated space system, might improve our understanding of tsunami propagation in the open ocean and possibly the efficiency of the future tsunami warning systems. Other applications of this technique are also found in planetology, especially with interesting prospects in the remote sensing of quakes on Venus.

## VII. REFERENCES

[1] Bolt, B.A., *Nature*, **202**, 1095-1096, 1964

[2] Davies, K., and D.M. Baker, *Journal of*

*Geophysical Research*, **70**, 1251-1253, 1965.
[3] Yuen, P. C., Weaver, P. F., Suzuki, R. K., and Furumoto, A. S., *J. Geophys. Res.*, **74**, 2256-2264, 1969.

[4] Weaver, P. F., Yuen, P. C., Prollss, G. W., and Furumoto, A. S., *Nature*, **226**, 1239-1241, 1970.

[5] Mannucci, A. J., *Radio Science*, **33**, 565-582, 1998.

[6] Calais, E., & Minster, J. B., *Geophys. Res. Lett.*, **22**, 1045-1048, 1995

[7] Afraimovich, E.L., Perevalova, N.P., Plotnikov, A.V., Uralov, A.M., *Annales Geophysicae*, **19**, 395-409, 2001

[8] Lognonné, P., C. Clévédé and H. Kanamori, *Geophys. J. Int.*, **135**, 388-406, 1998.

[9] Artru J., P. Lognonné et E. Blanc, *Geophysical Res. Lett.*, **28**, 697-700, 2001.

[10] Kherani A., N. Kamath, F. Crespon, R. Garcia and P. Lognonne, submitted to *Geophys. J. Int.*

[11] T. Farges, J. Artru, P. Lognonné, A. Le Pichon, *Choc*, **26**, CEA pub., 6-18, 2002.

[12] Artru, J. T. Farges, P. Lognonné, *Geophys. J. Int.*, **158**, 1067-1077, 2004.

[13] G. Occhipinti, Thèse de l'Institut de Physique du Globe de Paris, 2006.

[14] Ducic, V. J. Artru and P. Lognonné, *Geophys. Res. Lett.*, **30**, 1951, 2003.

[15] Ducic, Thèse de l'Université de Paris VII, 2003.

[16] Garcia, R., F. Crespon, V. Ducic, P. Lognonné, *Geophys. J. Int* **163**, 1049-1064, 2005.

[17] Crespon, thèse de l'Institut de Physique du Globe de Paris, 2007.

[18] More information on the Spectre project can be found at http://www.noveltis.fr/spectre

[19] Lognonné, P., Juliette Artru, Raphael Garcia, François Crespon, Vesna Ducic, Eric Jeansou, Giovani Occhipinti, Jérôme Helbert, Guilhelm Moreaux, Pierre-Emmanuel Godet, *Planet. Space. Science*, **54**, 528-540, 2006

[20] Peltier, W. R and Hines, C.O., *J. Geophys. Res.*, **81**, 1995-2000, 1976.

[21] J. Artru, V. Ducic, H. Kanamori, P. Lognonné and M. Murakami, *Geophy. J. Int.*, **160**, 840-848, 2005.

[22] G. Occhipinti, P. Lognonné, E. Kherani, H. Hébert, *Geophys. Res. Lett*, **33**, L20104, 2006.

[23] Garcia, R., P. Lognonné and X. Bonnin, *Geophys. Res. Lett.*, **32**, L16205, 2005

[24] Lognonné P. and C. Johnson, *Treatise in Geophysics*, 2007.
# Plaques Plongeantes et Panaches Dans le Manteau Inférieur

E. STUTZMANN¹, J. BESSE¹, A. DAVAILLE¹ & Y. REN¹

¹ Institut de Physique du Globe de Paris, 4 place Jussieu, 75252 Paris Cedex 5,
France
stutz@igp.jussieu.fr

**Résumé** La tomographique sismique fournit des images des structures du manteau profond. Les anomalies de vitesses rapides sont associées aux plaques froides qui plongent tandis que les anomalies de vitesses lentes correspondent aux panaches chauds qui remontent. Mais la tomographie donne accès uniquement à la position actuelle des grandes structures du manteau et pour comprendre la localisation et la morphologie des plaques plongeantes et des panaches il faut également considérer leur histoire temporelle. Celle-ci est accessible en considérant les observations géologiques de surface, la cinématique des plaques et les expériences de mécanique des fluides. En combinant ces observations, il est alors possible d'interpréter les images tomographiques en terme de convection thermochimique du manteau. Nous présentons ici l'interprétation des images tomographiques du manteau dans deux régions, sous les Amériques et dans la région indo atlantique.

Sous les Amériques, nous avons mesuré 16.000 temps d'arrivée différentes par cross corrélation en sélectionnant des séismes et des stations le long d'un corridor reliant les Aléoutiennes à l'Amérique du Sud. Un modèle de vitesse des ondes S a ensuite été obtenu par inversion de ces temps d'arrivée. La morphologie des plaques plongeantes a été déterminée en comparant les images tomographiques et la reconstruction de la position des frontières de plaques au cours des derniers 100 Ma. Vers 800 km de profondeur, notre modèle tomographique montre quatre anomalies de vitesses rapides que l'on peut associer aux subductions des plaques Nazca et Cocos sous l'Amérique du Sud et Centrale, à la plaque Farallon sous l'Amérique du Nord et à la plaque Pacifique sous les Aléoutiennes. A la base du manteau, la grande anomalie rapide, peut être associée avec la subduction, au Crétacé, de la plaque Farallon sous les Aléoutiennes et sous l'Amérique du Nord et Centrale et peut être avec la plaque
Phoenix sous l'Amérique du Sud. En remontant dans le manteau on peut suivre la fragmentation progressive de la plaque Farallon d'abord en deux puis en trois plaques. Cette interprétation est validée par le bon accord entre le volume de plaque estimé à partir des images tomographiques et à partir de la cinématique des plaques.

Dans la région indo atlantique, nous interprétons conjointement les images tomographiques, les observations géochimiques, géologiques et paléomagnétiques et les expériences de convection mantelliques pour retrouver l'histoire de la convection au cours des derniers 260Ma. Le « boite » indo atlantique contient a) une anomalie de vitesse lente à la base du manteau qui a la même forme que la Pangée il y a 250 Ma et qui se divise en plusieurs branches en remontant dans le manteau inférieur, b) un « superswell » sous l'Afrique du Sud, c) au moins 6 « points chauds primaires » reliés par une longue trace à des traps, d) de nombreux petits points chauds. Au cours des derniers 250 Ma la boite indo atlantique a été le siège de 10 épisodes de traps dont 7 sont relié à la cassure d'un continent. Les expériences de laboratoires montrent que les « superswells », les points chauds avec trace et les traps sont les différentes stades de la déstabilisation périodique d'une couche limite thermique chaude et chimiquement hétérogène à la base du manteau. On observe dans le manteau terrestre 7 instabilités et une récurrence d'environ 150 Ma ce qui est en accord avec les prédictions de la mécanique des fluide. Ainsi les remontées mantelliques sont actives et correspondent à des instabilités thermochimiques provenant de la base du manteau.

***Mots clefs*** Tomographie sismique, plaques plongeantes, panaches, manteau terrestre, convection

## Sinking Slabs and Plumes in The Lower Mantle

**Abstract** Tomographic models provide images of the deep mantle structures. Fast seismic anomalies can be associated with sinking slabs of subducted lithosphere whereas slow anomalies correspond to hot rising plumes. Tomographic images only give access to the present day location of slabs and plumes and if we want to understand their location and morphology, it is necessary to consider their evolution with time. To get access to the time variations, we have considered plate motion history, geological observations and fluid mechanic experiments. Combining these observations, it has been possible to reinterpret the tomographic images in terms of thermo-chemical mantle convection. In this paper we review results obtained in two
areas: beneath the Americas where slabs have been subducting for at least the last 120 Myr and in the indo-atlantic "box" which does not host any large scale present day subduction.

Beneath the Americas, we constructed 3D models of S wave velocity anomalies through travel time inversion. Earthquakes and receivers located along a wide corridor from Alaska to South America were selected. 16.000 differential travel times were measured by broadband waveform cross-correlation and then inverted to obtain a high resolution S wave tomographic model. We interpret the slab structures and unravel the subduction history by comparing our Vs tomographic images with the reconstructed plate motion from present-day up to 100 Ma. In the upper part of the lower mantle we recognize the imprint of the most recent phase of the plate history. Around 800 km depth, four distinct fast anomalies can be associated with subduction of the Nazca and Cocos plates beneath South and Central America respectively, the Juan de Fuca plate beneath North America, and the Pacific plate beneath the Aleutian islands. In the lowermost mantle, the large fast anomalies can be associated with Late Cretaceous subduction of the Farallon plate beneath the Aleutian islands and most of the Americas, and perhaps with the Phoenix plate beneath the southern part of South America. Going upwards, at 2000 km depth, the images record the fragmentation of the proto Farallon plate into the Kula plate in the north and into the Farallon plate in the north-east since 80 Ma. Around 1000 km depth, we observe separated fast anomalies interpreted as the Kula-Pacific, the Juan de Fuca, and the Farallon slabs. This interpretation is validated by the good agreement between the volume of slabs estimated from the tomographic images and from plate history reconstruction.

In the Indo-Atlantic box, we interpret the tomographic images of the mantle together with geochemical, geological and paleomagnetic observations and with fluid mechanic experiments. We attempt to unravel the pattern of convection in the Indo-Atlantic "box" as well as its temporal evolution over the last 260 Myr. The « box » presently contains a) a broad slow seismic anomaly at the CMB which has a similar shape as Pangea 250 Myr ago, and which divides into several branches higher in the lower mantle, b) a "superswell" centered on the western edge of South Africa, c) at least 6 "primary hotspots" with long tracks related to traps, and d) numerous smaller hotspots. In the last 260 Myr, this mantle box has undergone 10 traps events, 7 of them related to continental break up. Several of these past events are spatially correlated with present-day seismic anomalies and/or upwellings. Laboratory
experiments show that superswells, long-lived hotspot tracks and traps may represent three evolutionary stages of the same phenomenon, i.e. the episodic destabilization of a hot, chemically heterogeneous thermal boundary layer, close to the bottom of the mantle. When scaled to the Earth's mantle, the recurrence time of such phenomenon should be on the order of 100-200 Myr. Also, at any given time, the Indo-Atlantic box should contain 3 to 9 of these instabilities at different stages of their development, in agreement with observations. The return flow of the downwelling slabs, although confined to two main « boxes » (Indo-Atlantic and Pacific) by subduction zone geometry, may therefore not be passive, but rather take the form of active thermochemical instabilities.

***Key words*** Seismic tomography, sinking slabs, plumes, mantle, convection

# I. INTRODUCTION

La convection dans le manteau est caractérisée par des plaques froides qui plongent depuis la surface de la terre et des panaches chauds qui montent. Pour comprendre ce que deviennent les plaques qui plongent dans le manteau ou pour localiser les panaches montants, la tomographie sismique est un outil indispensable car elle apporte des images de l'intérieur de la terre. Si on suppose que les anomalies de vitesse sont dues principalement à des variations de température, on peut interpréter les images tomographiques dans le manteau profond en associant les anomalies de vitesses rapides à des structures froides qui plongent et les anomalies de vitesses lentes à des structures chaudes qui montent. On obtient alors une image instantanée des grandes structures du manteau. Ensuite,

pour relier ces structures à la convection mantellique, il faut considérer l'évolution temporelle du manteau. On obtiendra cette information à partir des observations paléomagnétiques et géologiques. On peut alors confronter les images tomographiques et les observations de surface avec les modélisations de mécanique des fluides pour comprendre le mode de convection dans le manteau et ainsi déterminer le rôle respectif des plaques plongeantes et des panaches mantelliques.

Les zones de subduction ont bougé au cours des 400 derniers millions d'années mais elles sont restées situées, pour la plupart, autour du Pacifique (Storey, 1995). Dans le manteau inférieur, on observe des anomalies positives de vitesse autour du Pacifique qui sont
interprétées comme des fragments de plaques océaniques qui auraient subducté dans le passé. Cet anneau de vitesses rapides a été observé par Masters et al. (1982) dans l'un des premiers modèles tomographiques. Les anomalies du géoïde présentent également un fort degré 2 qui s'explique bien par la localisation des plaques dans le manteau terrestre essentiellement autour du Pacifique (Ricard et al., 1993).

Les modèles tomographiques plus récents montrent que certaines plaques plongent dans le manteau inférieur tandis que d'autres semblent stagner dans la zone de transition entre 400 et 670 km de profondeur (e.g. Grand et al., 1997 ; van der Hilst et al., 1997 ; Fukao et al., 2001 ; Ren et al., 2007). Ainsi, les plaques ne plongent pas partout de la même manière et il reste à comprendre pourquoi. La cinématique des plaques, la migration des fosses de subduction, la température et la composition des plaques, ou encore le saut de viscosité à 670 km de profondeur sont autant de paramètres dont il faut tenir compte pour comprendre comment et jusqu'à quelle profondeur les plaques plongent dans le manteau ?

L'anneau d'anomalies positives de vitesse autour du Pacifique semble séparer le manteau en deux entités, la boîte Pacifique et la boîte Indo-Atlantique dans lesquelles on observe des anomalies de

vitesse négatives (e.g. Grand et al., 1997 ; Mégnin and Romanowicz, 2000 ; Ritsema and van Heijst, 2000). Les modèles tomographiques actuels montrent que ces deux régions lentes ne sont pas homogènes mais présentent des structures complexes qui sont encore mal connues. Ces anomalies de vitesse négatives sont interprétées comme des panaches chauds qui montent mais dont la morphologie des courants montants dans le manteau reste à définir. On s'interroge notamment sur l'origine, soit purement thermique, soit thermochimique, des courants montants. Pour répondre à cette question, il faut être capable de déterminer quelle est la part des hétérogénéités de température et de composition dans les images tomographiques du manteau (e.g. Saltzer et al., 2004). Une autre question ouverte concerne le lien entre les anomalies de vitesse négatives dans le manteau profond et les points chauds en surface. Certains auteurs pensent que tous les points chauds ont une origine asthénosphérique (Anderson, 2000) tandis que d'autres observent des anomalies lentes de vitesse sous la plupart des points chauds au moins jusqu'à la zone de transition et dans de nombreux cas jusqu'à la limite noyau-manteau (Montelli et al., 2004). Les deux cas co-existent sans doute dans la terre mais cela reste à démontrer.

C'est dans ce contexte que nous
avons étudié les plaques et les panaches
dans le manteau profond en nous focalisant
sur deux régions: sous les Amériques et
dans la région indo-atlantique. Sous les
Amériques, nous avons recherché dans le
manteau profond les anomalies de vitesse
rapide qui pourraient correspondre à des
morceaux de plaques issus de la
fragmentation de la grande plaque Farallon
qui a subducté dans cette région au cours
des 120 derniers Ma. Dans la région Indo-

Atlantique, nous avons étudié les
structures chaudes qui remontent depuis la
base du manteau en interprétant
conjointement les images tomographiques,
les observations de surface ainsi que les
modélisations de convection. Avant de
présenter ces deux études, nous faisons un
bref résumé des principales
caractéristiques de la convection
mantellique.

<figure><img src="image_74.png" /></figure>

**Figure 1:** Carte des séismes (étoiles) et des stations (triangles) utilisés pour obtenir le modèle tomographique. Les points chauds sont représentés par des points blancs. Les deux figures dans le coin gauche en bas représentent les positions des plaques durant les périodes 10-25 Ma et 84-94 Ma (d'après Lightow-Bertelloni and Richard, 1995).
## II. LA CONVECTION MANTELLIQUE

Le manteau terrestre est un système thermique complexe. La surface de la terre est à la température moyenne de 2° C et refroidit le manteau par le haut. L'existence d'un champ magnétique issu du noyau terrestre implique que celui-ci convecte et perd de la chaleur. Il réchauffe le manteau par le bas et la température à la limite noyau-manteau est estimée environ 3800° K. Si on considère un manteau terrestre de composition homogène et isovisqueux de la surface à la limite noyau-manteau, le nombre de Rayleigh associé au manteau est compris entre $10^7$ et $10^9$. Cela correspond à une convection thermique vigoureuse du manteau avec des plaques froides qui plongent depuis la surface et des panaches chauds qui montent depuis la base du manteau (Bercovici et al. 2000).

Entre les deux couches limites thermiques à surface et à la base du manteau, ce dernier est chauffé de l'intérieur par les éléments radiogéniques. La proportion de chauffage interne par rapport au chauffage par le bas est encore mal connue et on estime que 80 à 90% de la chaleur du manteau provient du chauffage interne. A la surface de la terre, la couche limite thermique froide doit évacuer la chaleur provenant du chauffage

par le bas et du chauffage interne. A la base du manteau, la couche limite thermique chaude transporte uniquement la chaleur provenant du noyau. Ainsi, le gradient de température dans la couche limite thermique supérieure est bien plus important que celui de la couche limite inférieure et c'est pourquoi les plaques plongeantes sont le moteur de la convection mantellique (Davies and Richards, 1992). En plus de l'incertitude sur la proportion entre le chauffage interne et le chauffage par le bas, d'autres paramètres ont également une influence sur la convection mantellique. Un paramètre important est la dépendance en température de la viscosité du manteau et notamment le saut de viscosité de plusieurs ordres de grandeur à la surface de la terre et d'un facteur 30 à 100 entre le manteau supérieur et inférieur. Il faut également tenir compte de la présence de changements de phase au niveau de la zone de transition ou encore de la sphéricité de la terre (voir Bercovici et al., 2000 pour une review). Enfin la présence possible d'hétérogénéités de composition dans le manteau implique que la convection mantellique est thermochimique et non uniquement thermique (e.g. Davaille, 1999).
<figure><img src="image_75.png" /></figure>

**Figure 2 :** *Gauche:* Cartes des perturbations de vitesse d'onde S dans les couches 750-850 km (a), 1000-1150 km (b), 1300-1450 km (c), 1750-1900 km (d), 2650-2890 km (e). *Droite:* cartes des frontières de plaque aux ages 10-25 Ma (f), 48-56 Ma (g), 56-64 Ma (h), 74-84 Ma (i), 84-94 Ma (j) (figure d'après Ren et al., 2007).
<figure><img src="image_76.png" /></figure>

**Figure 3:** Coupes à travers le modèle de perturbation de vitesse d'onde S le long des trajets noirs. Les lignes rouges et vertes représentent les positions passées des frontières entre les plaques de l'Océan Pacifique et les continents américains au cours des derniers 120 Ma.

De nombreux modèles du manteau sont proposés mais aucun n'est pour l'instant entièrement satisfaisant. Il existe des modèles de convection à une couche qui englobe l'ensemble du manteau, à deux couches avec un manteau supérieur séparé du manteau inférieur, à une couche et demi avec des passages intermittents de matière entre le manteau supérieur et inférieur (e.g. Machetel et Weber, 1990). Ces dernières années, plusieurs auteurs ont proposé des modèles dans lesquels la limite entre les deux couches se trouverait non pas dans la zone de transition mais

plutôt dans le manteau inférieur, voire à la base du manteau inférieur et les couches correspondraient à des régions de composition (de quoi ?) distinctes (Davaille, 1999, Kellog et al., 1999). Une synthèse des arguments en faveur de l'un ou l'autre des modèles pourra être trouvée dans Albarède et van der Hilst (1999).

Pour étudier la convection mantellique, nous allons dans la suite de cette article considérer d'abord les plaques plongeantes sous les Amériques puis les panaches dans la boîte indo-atlantique.
### III. LES PLAQUES PLONGEANTES SOUS LES AMÉRIQUES

Les zones de subduction du pourtour du Pacifique existent depuis plusieurs centaines de millions d'années. Elles se déplacent lentement au gré de l'ouverture des Océans Indien et Atlantique, en réduisant progressivement la superficie du Pacifique.

Nous avons choisi d'étudier les plaques plongeantes sous les Amériques car cette région présente une forte

concentration en séismes et stations, ce qui permet d'obtenir des modèles tomographiques à haute résolution. Le modèle présenté ici sera utilisé pour déterminer s'il est possible de corréler les anomalies de vitesse rapides dans le manteau terrestre avec l'histoire du mouvement des plaques à la surface de la terre.

#### III.1 Les données

Nous avons sélectionné les couples de séismes et stations le long d'un large corridor reliant les Aléoutiennes au Sud de l'Amérique du Sud en passant par l'Amérique du Nord et l'Amérique centrale (Figure 1). Nous avons utilisé des sismogrammes large-bande de 616 séismes enregistrés entre 1990 et 2002 par un total de 841 stations provenant des réseaux permanents GEOSCOPE, IRIS ainsi que d'expériences temporaires. Nous avons mesuré des temps d'arrivée différentiels

par cross corrélation. Pour chaque séisme, nous avons mesuré les temps d'arrivée différentiels des phases directes P et S entre les différentes paires de stations qui avaient enregistré le séisme. Nous avons également mesuré les temps d'arrivée différentiels PcP-P et ScS-S respectivement sur les composantes transverses des sismogrammes. Nous avons ainsi obtenu environ 18.000 mesures P1-P2 et S1-S2 et 1000 mesures PcP-P et ScS-S.

#### III.2 La Tomographie

Le manteau a été paramétré en cellule de 2° x 2° x 200 km. Les résidus de temps d'arrivée différentiels ont été inversés simultanément en utilisant l'algorithme de gradient conjugué itératif LSQR (Paige et

Saunder, 1982) pour déterminer la perturbation de lenteur de chaque cellule par rapport au modèle de référence 1D ak135 (Kennett et al. 1995). Deux paramètres de régularisation ont été
introduits dans l'inversion pour contrôler
l'amplitude des perturbations de lenteur
par rapport au modèle de référence ainsi
que le lissage du modèle. Ces paramètres
de régularisation ont été déterminés en
considérant la réduction de variance et le

paramètre χ2 ainsi qu’un grand nombre de
tests synthétiques. Une description
complète des données, de la méthode
d’inversion, des tests synthétiques et du
modèle peut être trouvée dans Ren et al.
(2007).

<figure><img src="image_77.png" /></figure>

**Figure 4.** Carte bathymétrique de la région Indo-Atlantic présentant les points chauds actifs (hexagones noirs et blancs) et les traps (gris). Les points chauds primaires (hexagones noirs) sont reliés à leur trap par une ligne continue. Le « superswell » africain est entouré d'une ligne blanche. JM : Jan Mayen ; I : Iceland ; Ei : Eifel ; Az : Azores ; Ber : Bermuda ; GM : Great Meteor ; Can : Canary ; Ho : Hoggar ; Tib : Tibesti ; CV : Cape Verde ; Da : Darfur ; Af : Afar ; Cam : Cameroon ; Fer : Fernando ; As : Ascension ; Com : Comores ; Td : Trinidad ; SH : St Helena ; Re : Reunion ; Tris : Tristan ; Ve : Vema ; Dis : Discovery ; Crz : Crozet ; Met : Meteor ; Bou : Bouvet ; Ma : Marion ; Ker : Kerguelen. (d'après Davaille et al., 2005)

### III.3 Identification des plaques dans le modèle tomographique.

Différentes plaques océaniques ont
subducté sous les Amériques durant les
derniers 150 Ma (Engebretson, 1985 ;
Lightow-Bertelloni and Richard, 1995).
Durant l'ère Mésozoïque, l'ancien océan
Pacifique était constitué des plaques
Izanagi, Farallon, Phoenix et Pacifique

(Figure 1). Les plaques Pacifique et
Farallon existent encore aujourd'hui tandis
que les autres ont entièrement disparu dans
le manteau. Au cours de sa subduction
sous les Amériques, la grande plaque
Farallon s'est progressivement cassée pour
former plusieurs plaques plus petites. Vers
80 Ma, une ride a séparé la plaque Farallon en deux, et la partie nord est devenue la plaque Kula qui a entièrement disparu de la surface de la terre il y a environ 48 Ma. Puis vers 30 Ma, la partie restante de la plaque Farallon s'est à nouveau fragmentée pour créer les plaques Nazca, Coco et Juan de Fuca dont nous continuons à observer aujourd'hui la subduction sous les Amériques.

Au cours des derniers 120 Ma, les plaques de l'Océan Pacifique sont arrivées à la fosse de subduction de l'ouest des Amériques avec une vitesse moyenne à la surface d'environ 8cm/an. Supposons que les plaques aient plongé dans le manteau supérieur à la vitesse de leur convergence à la surface puis qu'elles aient été ralenties d'un facteur 4, lorsqu'elles ont pénétré dans le manteau inférieur plus visqueux. Dans ce cas simplifié, une plaque qui plongerait verticalement atteindra la base du manteau en 120 Ma. Les plaques plus anciennes se trouveraient alors à la base du manteau et les anomalies de vitesse positives plus superficielles que l'on observe dans les modèles tomographiques pourraient être corrélées avec des épisodes de subduction plus récents.

Nous avons comparé les anomalies de vitesse positives de notre modèle tomographique avec les positions des frontières de plaques au cours des derniers 120 Ma en utilisant la reconstruction de

l'histoire des plaques de Lithgow-Bertelloni et Richard (1995) (Figure 2). Dans la partie supérieure du manteau inférieur, nous distinguons plusieurs anomalies de vitesse rapides que nous pouvons associer avec les subductions récentes des plaques Nazca et Coco sous les Amériques du Sud et Centrale et de Juan de Fuca et Pacifique sous l'Amérique du Nord et sous les Aléoutiennes (Figure 2a et 2f). L'anomalie de vitesse positive au niveau des Caraïbes correspond à la subduction de la plaque Atlantique sous l'arc des Antilles.

Pour retrouver l'histoire des plaques plus anciennes, il est plus simple de partir de la base du manteau et de remonter vers la surface pour tenter de suivre la fragmentation de la plaque Farallon. A la base du manteau, les anomalies de vitesse d'onde S sont essentiellement positives avec deux exceptions entre 0° et 15°S et entre 30° et 45°N (Figure 2e). Ces structures rapides pourraient correspondre aux plaques Phoenix au Sud et Farallon au Nord ou à toute autre plaque plus ancienne qui aurait subducté sous les Amériques. Entre 120 et 80 Ma, la plaque Farallon se déplaçait vers le Nord Nord-Est et subductait sous l'Amérique du nord, l'Amérique centrale et sous les Aléoutiennes (Figure 2j).

Vers 80 Ma, la plaque Farallon s'est brisée en deux plaques (Figure 2i). La
partie nord est devenue la plaque Kula qui a continué sa progression vers le nord tandis que le mouvement de la plaque Farallon s'est progressivement réorienté vers le Nord Est autour de 65 Ma. La localisation de l'ancienne ride entre les plaques de Kula et Farallon est mal connue et Engebretson et al. (1985) l'ont placée dans la région hachurée de la figure 2i entre le Mexique et le centre de l'Amérique du Nord. La carte d'anomalie de vitesse d'onde S dans la couche 1750-1900 km de profondeur (figure2d) montre un coin sous l'Amérique du Nord entre deux anomalies de vitesse positives alignées suivant des directions différentes. Si ce coin marque la frontière entre les deux plaques, nous pouvons associer l'anomalie rapide au Nord à la plaque Kula et l'anomalie au Sud à la plaque Farallon. Dans ce cas, nous avons retrouvé la signature de l'ancienne ride: elle se trouverait à l'extrémité nord de la zone hachurée.

La frontière entre les plaques Kula et Farallon s'est ensuite déplacée vers le Nord-Ouest durant le période 65-56 Ma et nous pouvons la suivre dans notre modèle tomographique en remontant dans le manteau à la profondeur 1300-1450 km (figure 2c et 2h). Nous retrouvons également le coin entre les deux anomalies positives mais il est situé plus au Nord

Ouest qu'à la profondeur précédente, ce qui confirme notre interprétation.

Enfin vers 48 Ma, la dorsale entre la plaque Kula et Pacifique a cessé d'être active et la plaque Kula a été intégrée à la plaque Pacifique (Engebretson et al., 1985) qui se déplaçait vers le nord. A la même période, la plus grande partie de la plaque Farallon avait déjà subducté sous les Amériques et la plaque Pacifique est entrée en contact avec le continent Américain au niveau de la côte ouest du Canada et du Golfe de Californie (Figure 2g). La partie restante de la plaque Farallon entre ces deux régions s'appelle maintenant la plaque Juan de Fuca et continue à subducter sous l'Amérique du Nord dans la zone de subduction des Cascades. Nous observons sur notre modèle tomographique à la profondeur 900-1200 km de profondeur plusieurs anomalies de vitesse positives séparées par des anomalies de vitesse négatives que nous interprétons comme les anciennes plaques Pacifique/Kula au Nord, Farallon/Juan de Fuca sous le Canada, et Farallon sous la partie est de l'Amérique du Nord et sous l'Amérique centrale (Figure 2b).

Enfin, à la profondeur 750-850 km, nous retrouvons les anomalies rapides associées aux plaques Juan de Fuca, Cocos et Nazca qui sont les derniers fragments de la grande plaque Farallon (Figure 2a et 2f)
dont nous avons parlé au début de l'interprétation de notre modèle.

Nous observons une très bonne corrélation entre le mouvement des plaques au cours du temps et notre modèle tomographique à différentes profondeurs. Pour valider quantitativement cette corrélation, nous avons comparé le volume de plaque observé dans nos images tomographiques avec celui estimé à partir du mouvement des plaques. Le long de différents profils, nous avons déterminé le mouvement des plaques de l'actuel océan Pacifique par rapport aux Amériques en utilisant, d'une part le référentiel des points chauds et d'autre part le circuit de plaques passant par l'Antarctique. Le détail des reconstructions se trouve dans l'article de Ren et. al (2007). Nous avons montré que les deux méthodes de reconstruction donnent des résultats similaires pour les 80 derniers millions d'années. Pour les périodes plus anciennes, les vitesses de convergence calculées par la méthode de circuit de plaques sont systématiquement supérieures à celles obtenues par la méthode utilisant les référentiels de points chauds fixes. La

différence de vitesse entre les deux méthodes peut atteindre 8 cm/an. Ces reconstructions nous ont permis de faire des coupes dans notre modèle tomographique qui sont perpendiculaires à la position présente et passée de la fosse de subduction en tenant compte notamment du brusque changement de direction de convergence vers 80 Ma (Figure 3). En supposant que les plaques ont une épaisseur d'environ 100 km, nous avons estimé la longueur totale de plaques dans le manteau à partir des anomalies de vitesse rapides de nos images tomographiques le long des 3 coupes de la figure 3. Ces longueurs sont en bon accord avec les longueurs prédites à partir des vitesses de convergence des plaques, soit 16.000 km sous l'Amérique du Nord, 12.000 km sous l'Amérique Centrale et 5000 km sous l'Amé- rique du Sud (coupe A-A', B-B' et C-C' de la figure 3). Ce bon accord confirme que l'on a effectivement retrouvé dans notre modèle tomographique les différents morceaux de plaque qui résultent de la fragmentation de la plaque Farallon au cours des derniers 90-120 Ma.

## IV. LES PANACHES MANTELLIQUES

Si les plaques qui subductent étaient le seul moteur de la convection mantellique, on devrait observer dans le manteau des

vitesses rapides associées aux plaques froides à proximité des zones de subduction et le reste du manteau devrait
être uniformément plus lent car il correspondrait à un courant de retour passif. Les modèles tomographiques récents (e.g. Grand et al., 1997; Mégnin and Romanowicz, 2000; Ritsema and van Heijst, 2000) montrent que dans la région Indo-Atlantique et sous l'océan Pacifique le manteau n'est pas uniformément plus lent et on observe différents extrema dans

les régions d'anomalies de vitesse négatives qui pourraient correspondre à des panaches montant depuis le noyau. Par ailleurs, les anomalies de bathymétrie et de topographie associées aux « superswells » sud-africain et du Pacifique (Mc Nutt, 1998) ne peuvent pas non plus être expliquées par une convection uniquement contrôlée par la subduction.

<figure><img src="image_78.png" /></figure>

**Figure 5:** Reconstruction paléomagnétique à différentes périodes de la position des continents par rapport à l'axe de rotation de la terre et position des traps (d'après Besse et Courtillot, 2002).
<figure><img src="image_79.png" /></figure>

**Figure 6** : Modèle tomographique TXBW de Grand (2000). Sur la carte à la base du manteau, les 4 extrema de vitesse sont en vert. L'anomalie lente en forme de banane est en blanc. Les coupes 2-5 sont le long des trajets en couleurs. 2) coupe le long de l'Atlantique, 3) coupe du Cap Vert à l'Afar, 4) coupe de l'Afrique du Sud à l'Afar, 5) coupe de l'Afar à Kerguelen.

Nous avons étudié le manteau dans
la région Indo-Atlantique en comparant la
localisation des anomalies de vitesse
négatives des modèles tomographiques

globaux du manteau avec les observations
géologiques, l'histoire du mouvement des
plaques à la surface et les modèles de
convection, dans le but de comprendre les
paramètres qui contrôlent la morphologie
et la localisation des anomalies de vitesse
négatives dans le manteau de cette région.
Nous présentons ici les principaux

résultats obtenus. L'étude complète est
publiée dans l'article de Davaille et al.
(2005).

IV.1. les observations de surface

La boîte Indo-Atlantique est définie
comme la région qui est séparée de la boîte
Pacifique par les subductions du pourtour
du Pacifique. Cette région est couverte par
50% de continents et 50% d'océans. Elle
se raccourcit le long des zones de

convergence des Alpes et de l'Himalaya et
s'agrandit le long des rides des océans
Indien et Atlantique. La variation de
surface de cette région n'a été que de 0.6
km²/an au cours des derniers 200 Ma.

<figure><img src="image_80.png" /></figure>

**Figure 7** : A) Position de la Pangée il y a 250 Ma superposée à la carte tomographique à la base du manteau (modèles TXBW). B) Localisation des traps au moment de leur éruption projetée sur la carte tomographique du modèle TXBW à la base du manteau. Le cercle en tirets verts représente la localisation d'un événement qui s'est produit il y a 145 Ma, qui a créé le rifting de l'Ouest et du Centre de l'Afrique et qui est peut être associé à la trace du point chaud de St Helen (e.g. Wilson 1992).
Nous trouvons dans la boîte Indo-
Atlantique 25 points chauds
essentiellement localisés au Nord Ouest du
continent africain et près des rides de
l'Atlantique et de l'ouest de l'océan Indien
(Figure 4). Ils se trouvent tous au dessus
ou à la périphérie d'une anomalie de
vitesse lente à la base du manteau. Dans
cette région, Courtillot et al. (2003)
définissent six points chauds « primaires »
qui proviennent de la couche D'' à la
limite manteau-noyau et sont associés à un
volcanisme basaltique de type traps. Ces
points chauds primaires sont l'Islande,

Tristan, Marion, Kerguelen, la Réunion et
l'Afar. Les traps de CAMP, de Sibérie et
d'Emeishan, par contre, ne sont pas
associés à des points chauds connus. Une
autre particularité de la région Indo-
Atlantique est l'existence d'un bombement
topographique au sud de l'Afrique du Sud
qui s'étend dans l'océan Atlantique où la
bathymétrie est anormalement élevée
(Nyblade and Robinson, 1994). Ce
« superswell » est interprété comment un
soulèvement actif provenant de structures
qui remonterait du manteau inférieur
(Lightow-Bertelloni and Silver, 1998).

IV.2 Cinématique des plaques et des continents

La reconstruction de la cinématique des
plaques et des continents est détaillée dans
l'article de Besse et Courtillot (2003). Le
super continent de la Pangée est resté
assemblé entre 320 Ma et 200 Ma (Figure
5a). Vers 180 Ma, l'ouverture de
l'Atlantique central a divisé la Pangée en
deux continents, la Laurasia au nord et le
Gondwana au Sud (Figure 5b, c). Vers 160
Ma, le Gondwana s'est lui aussi divisé en
deux pour créer un bloc Afrique-
Amérique-du-Sud et un bloc Madagascar-
Australie-Inde-Antarctique (Figure 5c),
suivi à partir de 130 Ma de l'ouverture de
l'Atlantique Sud séparant l'Amérique du
Sud de l'Afrique. Vers 110 Ma, l'Inde puis
progressivement l'Australie et la Nouvelle
Zélande se sont détachés de l'Antarctique.

L'Inde a ensuite rapidement migré vers le
Nord depuis 84 Ma en laissant sur place
Madagascar et les Seychelles (Figure
5c,e). L'Atlantique Nord s'ouvre en
plusieurs phases à partir de 62 Ma
séparant d'abord le Groenland de
l'Amérique du Nord, puis de l'Europe
(Figure 5e). L'ouverture des océans
Atlantique et Indien se poursuit et
l'épisode le plus récent est l'ouverture de
la mer rouge et du Golfe d'Aden depuis 30
Ma.

La séparation d'un super continent
en entités plus petites est généralement
expliquée par l'action de courant de
convection se développant sous une
lithosphère continentale isolante (Guillou
et Jaupar, 1995). Les grandes cassures
continentales sont souvent associées à du volcanisme de type trap provenant du manteau profond. La région Indo Atlantique a été le siège de 10 traps au cours des derniers 260 Ma, dont 7 sont reliés à la cassure d'un continent (Figure 5). CAMP est associé à l'ouverture de l'Atlantique central vers 200 Ma. Le volcanisme associé aux traps de Karoo, Ferrar et Tasman a précédé la cassure du Gondwana Est. Vers 133 Ma, les traps de Parana et Etandeka provenant du volcan Tristan précédent l'ouverture de l'Atlantique Sud. Les traps du point chaud de Kerguelen, datés à environ 118 Ma, sont situés dans la région de Rajmahal en

Inde et de Bundury en Australie; elles ont permis la séparation de l'Inde, de l'Australie et de l'Antarctique. Dans cette même région, le point chaud de Marion est entré en éruption vers 88 Ma en créant les traps de Madagascar et le volcan de la Réunion a produit les traps du Deccan vers 65 Ma. Durant cette période, l'Inde s'est séparée des Seychelles, au cours de sa migration vers le Nord. Les traps de BTIP et du Groenland associés à l'Islande se sont mis en place vers 62 Ma, juste avant l'ouverture de l'Atlantique Nord. Enfin le point chaud de l'Afar est le plus récent, les traps d'Ethiopie et du Yemen étant datés à 30 Ma.

### IV.3 Les modèles tomographiques

Nous avons considéré trois modèles tomographiques globaux récents TXBW (Grand et al., 1997), SAW24B16 (Mégnin and Romanowicz, 2000) et S20RTS (Ritsema and van Heijst, 2000) qui ont été obtenus à partir de jeux de données différents (temps d'arrivée ou forme d'onde). Les modèles ont été paramétrés différemment (splines ou blocs) et les théories reliant les données aux modèles sont également différentes (théorie des rais ou noyaux de sensibilité 2D). Nous avons considéré uniquement les structures de grande échelle présentes dans les trois modèles.

Les modèles tomographiques globaux montrent sous la région Indo-Atlantique une large anomalie de vitesse lente juste au-dessus de la couche D'' de forme allongée entre le nord de l'Islande et le sud-ouest de l'océan Indien, en passant sous l'Atlantique Est et une grande partie de l'Afrique (Figure 6). Dans cette anomalie lente, on distingue quatre extrema de vitesse lente (A, B, C, D sur la figure 6) et en remontant vers la surface, la large anomalie se divise en plusieurs branches que l'on peut suivre jusqu'à des profondeurs différentes suivant les régions.
<figure><img src="image_81.png" /></figure>

**Figure 8:** Interprétation en terme de mécanique des fluides du modèle de tomographique TXBW. 1) Les zones d'influences des instabilités thermo-chimiques, sont représentées par des disques de diamètre 2000 et 3500 km (centrés sur A, B, C, D, les Afars, la Réunion, l'Islande), superposé au modèle TXBW à 2700 km de profondeur. 2) coupe à travers l'Atlantique ; 3) coupe du Cap Vert à l'Afar ; 4) coupe de l'Afrique du Sud à l'Afar ; 5) coupe de l'Afar à Kerguelen en passant par la Réunion. Les instabilités contenant un fluide légèrement hétérogène en composition sont en rouge et celles fortement hétérogènes sont en vert.
## IV.4 Modélisation de la convection

Les expériences de convection analogique d'Anne Davaille lui ont permis de montrer que différentes formes de remontée mantellique peuvent coexister dans le manteau terrestre, depuis les petits panaches d'origine purement thermique, jusqu'aux gros dômes d'origine thermochimique qui pourraient venir d'une couche de composition différente à la base du manteau et correspondrait au volcanisme de type trap (Davaille, 1999). Suivant le contraste de densité de cette couche à la base du manteau, il pourrait

également y avoir des structures chaudes dans le manteau profond trop dense pour atteindre le manteau supérieur mais qui pourraient créer de la topographie dynamique à la surface. Les remontées mantelliques sont des phénomènes épisodiques et Lebars et Davaille (2004) ont établi des lois d'échelle pour le manteau terrestre. Celles-ci prédisent pour le manteau terrestre de la région Indo-Atlantique la présence de 3 à 9 instabilités qui monteraient dans le manteau avec une périodicité de 100 à 200 Ma.

## IV.5 Convection du manteau dans la boîte-Indo-Atlantique

Nous pouvons maintenant comparer les observations de surface, l'histoire des plaques et des continents de la région Indo-Atlantique et les images tomographiques aux prédictions de la mécanique des fluides. Nous remarquons tout d'abord que la forme et la localisation de l'anomalie lente à la base du manteau ressemblent de façon remarquable à celles de la Pangée il y a 250 millions d'années (Figure 7a) comme si le continent avait servi de couvercle en dessous duquel le manteau serait resté plus chaud. Nous observons également que la plupart des épisodes de traps qui ont eu lieu au cours des 250 derniers millions d'années se situent à l'aplomb des bords de cette large

anomalie lente (Figure 7b).

Les modèles tomographiques montrent quatre instabilités chaudes sous la forme d'anomalies de vitesse négatives visibles dans l'ensemble du manteau ou au moins le manteau inférieur sous les points chauds primaires d'Islande, d'Afar, de Kerguelen et de la Réunion (Figure 8). Nous observons également deux anomalies lentes de forme cylindrique qui remontent de la base du manteau jusqu'à la zone de transition sous la région du Cap-Vert-Canaries-Açores et sous les points chauds de Marion et Crozet. Enfin, nous observons dans le manteau profond une anomalie cylindrique sous l'Afrique équatoriale à l'aplomb du point C (figures
8.1 et 8.4) qui ne peut être reliée à aucun
point chaud connu. Au total, nous
dénombrons 7 instabilités (centres des
cercles bleus de la figure 8.1).

Ce nombre est en accord avec les
prédictions de la mécanique des fluides qui
prévoit entre 3 et 9 instabilités. Ainsi les
anomalies de vitesse lente observées dans
le manteau pourraient correspondre à de
gros dômes ancrés dans une couche fine,
plus visqueuse et légèrement plus dense
(0,1-0,6%) située à la base du manteau. La
mécanique des fluides prédit également un
diamètre entre 2000 et 3500 km pour ces
dômes (disques bleus sur la figure 8.1). Ce
diamètre qui est en accord avec celui des
dômes que nous voyons dans le manteau à
partir des modèles tomographiques
(figures 8.2 à 8.4) ainsi qu'avec la distance
entre les différents dômes (figure 8.1).

Nous observons également à la
base du manteau sous l'Afrique du Sud
une anomalie de vitesse lente allongée en
forme de banane qui pourrait ressembler à
une ride (en vert sur la figure 8). Cette
structure pourrait correspondre à la
remontée de matériau chaud mais dont le
contraste de densité d'origine chimique

serait plus élevé (δp~0,6-1,4%) qui
stagnerait dans le manteau inférieur mais
pourrait néanmoins être à l'origine de la
topographie dynamique nécessaire pour
expliquer le « superswell » africain.

En comparant la localisation des
traps au moment ou ils sont arrivés à la
surface, nous remarquons que CAMP âgé
de 200 Ma est arrivé à la surface de la terre
là où se trouvent aujourd'hui les points
chauds d'Atlantique central âgés de 40
Ma. Plus au Sud, les traps de Karoo
coïncidaient parfaitement au moment de
leur éruption, il y a 183 Ma, avec le
« superswell » Africain de 30 Ma. Ainsi on
semble observer une récurrence de 150 Ma
entre deux remontées successives depuis la
base du manteau. Cela est en accord avec
la mécanique des fluides qui prédit des
temps de récurrence entre 100 et 200 Ma.

Le bon accord entre les
observations et les prédictions de la
mécanique des fluides montrent que les
anomalies lentes observées dans le
manteau sont la signature d'une
déstabilisation périodique d'une couche
limite thermique chaude et chimiquement
hétérogène à la base du manteau.

V. CONCLUSIONS

L'amélioration de la résolution des
modèles tomographiques du manteau et
leur comparaison avec les observations en

surface, la cinématique des plaques et la
mécanique des fluides permettent de
comprendre de mieux en mieux la
convection mantellique. Nous avons ainsi
pu observer dans le manteau les différents
morceaux de plaques qui ont résultats de la
fragmentation de la plaque Farallon au
cours des derniers 90-120 Ma. Ce type
d'étude va permettre par exemple
d'apporter des contraintes sur la
localisation de frontières de plaque qui ont
disparu de la surface de la terre.

L'amélioration de la résolution des
modèles va également permettre de
déterminer la nature des anomalies de
vitesse rapide associées aux slabs. Les
plaques océaniques sont froides à la
surface de la terre et elles se réchauffent
progressivement quand elles plongent dans
le manteau. Si les anomalies de vitesse
positives autour du Pacifique observées à
la base du manteau correspondent à des
plaques qui ont subducté il y a 100 voire
200 millions d'années, on peut supposer
que ces anomalies de vitesse positives ont
une origine chimique et non plus
uniquement thermique. Cela reste à
démontrer et des modèles tomographiques
de la vitesse des ondes P et S avec une
haute résolution sont indispensables pour
cela. Ces modèles permettront également
des études comparatives des différentes
plaques en subduction dans le but de
comprendre le rôle des différents
paramètres qui contrôlent leur descente
dans le manteau.

Par rapport à l'étude des plaques

plongeantes, l'identification des courants
chauds ascendants dans le manteau
profond à partir des anomalies de vitesse
négatives est moins avancée et plus
polémique. Les points chauds, les traps, les
« superswells » sont la signature en surface
de courants chauds ascendants mais
l'origine en profondeur, asthénosphérique
ou profonde, des courants chauds est un
sujet très débattu (e.g. Anderson, 2000,
Montelli 2005). Par exemple, malgré des
expériences sismologiques dédiées à
l'étude du panache sous Hawaii, une
anomalie de vitesse lente dans l'ensemble
du manteau n'a pas été observée. Mais
cette observation est instantanée tandis que
les panaches mantelligues sont des
phénomènes transitoires qui ne sont pas
tous au même stade de leur évolution
temporelle. S'ils ne sont plus alimentés par
la couche limite thermique dont ils sont
issus, le conduit reliant le point chaud et
cette couche peut avoir disparu
partiellement ou complètement (e.g.
Silveira et al., 2006).

Dans la région Indo-Atlantique,
nous avons considéré à la fois les images
actuelles du manteau et le volcanisme en
surface qui nous a donné accès aux
variations temporelles. Nous avons montré
que les anomalies de vitesses négatives
observées dans cette région peuvent être
interprétées comme la signature de la
déstabilisation périodique d'une couche
limite thermique chaude et chimiquement
hétérogène à la base du manteau.

Nous avons pour l'instant identifié
des gros dômes du manteau profond. Mais
suivant la nature thermique ou thermo-
chimique des hétérogénéités du manteau, il

peut exister toute une zoologie de courants
chauds ascendants dans le manteau
(Davaille, 2003) de formes très diverses.
Ces courants chauds ascendants sont à
identifier -le cas échéant- dans les modèles
tomographiques.

VI. BIBLIOGRAPHIE

Albarède, F. and R. van der Hilst, 1999.
On the nature and scale of mantle
convection. EOS Trans. Amer.
Geophys. Union 80, 535-539.

Anderson D. L., 2000. Thermal State of
the Upper Mantle; No Role for
Mantle Plumes, Geophys. Res.
Lett., 27, 22, 3623-3626

Bercovici, D., Y. Ricard, and M.
Richards, 2000. *The relation
between mantle dynamics and plate
tectonics: a somewhat biased
primer*, AGU monography, 121,
Richards, Gordon and van der
Hilst, eds, 5-46.

Besse J. and V. Courtillot, 2002. Apparent
and true polar wander and the
geometry of the geomagnetic field
over the last 200 Myr J. Geophys
Res, 107, B11 2300, doi:10.1029
/2000JB000050.

Courtillot V., A. Davaille, J. Besse and J.
Stock, 2003. Three distinct types of
hotspots in the Earth's mantle,
Earth Planet. Sci. Lett., 205 295-

308.

Davaille A., 1999. Simultaneous
generation of hotspots and
superswells by convection in a
heterogeneous planetary mantle,
Nature, 402, 756-760.

Davaille A., E. Stutzmann, G. Silveira, J.
Besse and V. Courtillot, 2005.
Convective Pattern under the Indo-
Atlantic box, Earth Planet. Inter.
Let., 239, 233-252.

Davies G. and M. Richards, 1992. Mantle
convection, J. Geology, 100, 151-
206.

Dziewonski and Anderson, 1981.
Preliminary reference Earth model,
Phys. Earth Planet. Inter., 25, 297-
356.

Engebretson D., A. Cox and R. Gordon,
1985. Relative motions between
oceanic and continental plates in
the pacific basin., Spec. Pap. Geol.
Soc. Am., 2006, 1-59.

Fukao Y. S. Widiyantoro and M.
Obayashi, 2001. Stagnant slab in
the mantle transition zone, Review of Geophysics, 39,3, 291-323.

Guillou L. and C. Jaupart, 1995. On the effect of continents on mantle convection," J. geophys. res., 100, B12, 217-238.

Grand S., R. van der Hilst, and S. Widiyantoro, 1997. Global seismic tomography: a snapshot of convection in the earth. GSA Today 4, 1-7.

Kellogg L., B. Hager and R. van der Hilst, 1999. Compositional stratification in the deep mantle, Science, 283,1881-1884.

Kennett B.L.N. and Engdahl E.R., 1991. Travel times for global earthquake location and phase identification, Geophys J Int., 105, 429-465.

Kennett B.L.N., Engdahl E.R. and Buland R., 1995. Constraints on seismic velocities in the Earth from travel times, Geophys J Int, 122, 108-124.

LeBars M. and A. Davaille, 2004. Large interface deformation in two-layer thermal convection of miscible viscous fluids, J. Fluid Mechanics 499, 75-110.

Lithgow-Bertelloni C. and M. Richards, 1995. Cenozoic plate driving forces, Geophys. Res. Lett. 22 1317-1320.

Lithgow-Bertelloni C. and P. Silver, 1998. Dynamic topography, plate

driving forces and the African superswell, Nature, 395, 269-272.

McNutt M., 1998. Superswells, Rev. Geophys., 36, 211-244.

Machetel P. and P. Weber, 1991. Intermittent layered convection in a model mantle with an endothermic phase change at 670 km, Nature, 350, 55-57.

Masters G., T. Jordan, P. Silver and F. Gilbert, 1982. Aspherical Earth structure from fundamental spheroidal-mode data, Nature 298, 609-613.

Mégnin C. and B. Romanowicz, 2000. A model of shear velocity in the mantle from the inversion of waveforms of body, surface and higher mode waveforms, Geophys. J. Int. 143 709.

Montelli R., G. Nolet, F. Dahlen, G. Masters and E. Engdahl, S.-H. Hung, 2004. Finite-frequency tomography reveals a variety of plumes in the mantle, Science 303, 338-343.

Morgan W., 1971. Plate motions and deep mantle convection, Nature, 230, 42-43.

Nyblade A. and S. Robinson, 1994. The African Superswell, Geophys. Res. Lett., 21,765-768.

Paige C. and M. Saunders, 1982a. LSQR: An Algorithme for Sparse Linear
Equations and Sparse Least Squares, ACM Transactions on Mathematical Software, 8, 43-71

Paige C. and M. Saunders, 1982b. ALGORITHM 583 LSQR: Sparse Linear Equations and Sparse Least Squares Problems ACM Transactions on Mathematical Software, 8, 195-209

Ren Y. E. Stutzmann, R. van der Hilst and J. Besse, 2007. Understanding seismic heterogeneities in the lower mantle beneath the Americas from seismic tomography and plate tectonic history, J. Geophys. Res., in press.

Ricard Y., M.A. Richards, C. Lithgow – Bertelloni & Y. Le Stunff, 1993. A geodynamical model of mantle density heterogeneity, J. Geophys. Res., 98, 21895-21909.

Ringwood A., 1962. A model for the upper mantle J. Geophys. Res. 67, 857-867

Ritsema J. and H. van Heijst, 2000. Seismic imaging of structural heterogeneity in Earth's mantle:

Evidence for large-scale mantle flow, Science Progress, 83, 243-259.

Romanowicz B., 2003. Global mantle tomography: progress status in the last 10 years, Annu. Rev. Geoph. Space Phys, 31, 1, 303

Saltzer R., E. Stutzmann, and R. van der Hilst, 2004. Poisson's ratio beneath the Alaska from the surface to the core mantle boundary, J. Geophys. Res. (B6): B0630110.1029/2003JB0027 12.

Silveira G., E. Stutzmann, A. Davaille, J.-P. Montagner, and L. Mendes-Victor, 2006. Azores hotspot deep signature in tomographic models., J. Volcano. Geotherm. Res., 156, 23-34.

Storey B., 1995. The role of mantle plumes in the continental break up: case hisotries from Gondwanaland, 377, 301-308.

Van der Hilst R., S. Widyantoro and E. Engdahl, 1997. Evidence for deep mantle circulation from global tomography, Nature, 386, 578-584.
# Seismic Imaging Using Ambient Noise

L. STEHLY<fnref n="1" />,<fnref n="2" />, F. BRENGUIER<fnref n="1" />, P. ROUX<fnref n="1" />,

M. CAMPILLO<fnref n="1" /> & Y. COLIN DE VERDIERE<fnref n="3" />

**Abstract** In the present paper we show why and how cross-correlation of noise records can be used for seismic imaging. We first define theoretically the conditions required to observe the emergence of the Green functions between two stations A and B from the cross-correlation of noise records only. In a second part we study the origin of the Earth seismic noise in the period band 5-20 s. Our result indicates that the Earth seismic noise is spread over large surface and can actually be used for seismic imaging. We test the accuracy of the travel time measurements of surface wave performed on the correlation of seismic noise by analyzing 11 years of continuous noise records at three Californian stations. We show that the Green functions emerging from the correlation are robust and can even be used to routinely detect instrumentals clock errors at the stations. We then present two applications of noise correlation to seismic imaging using surface waves at the regional scale in Western Europe as well as at the local scale in the Piton de la Fournaise volcano. Finally, we conclude on the evidence of P-waves extracted from noise correlation at small scale in the Parkfield seismic network.

**Key words** Seismology, surface waves, tomography, seismic noise, data control quality, correlation

## Tomographie à Partir de Bruit de Fond Sismique

**Résumé** Les méthodes de tomographie sismique habituellement utilisées se basent sur des enregistrements de séismes. De là découlent leurs principales faiblesses : la répartition inhomogène des tremblements de terre limite la résolution qu'il est

<footnote n="1">¹ Université Joseph Fourier, CNRS, LGIT BP 53, 38041 Grenoble, France
[lstehly@obs.ujf-grenoble.fr](mailto:lstehly@obs.ujf-grenoble.fr)</footnote>

<footnote n="2">² CEA/DASE, BP 12 91680 Bruyères-le-Châtel, France</footnote>

<footnote n="3">³ Institut Fourier, CNRS-UJF BP 74, 38041 Grenoble, France</footnote>
possible d'atteindre, en particulier dans les régions asismiques. La faible occurrence des séismes interdit en outre d'étudier des changements intervenant dans le milieu, par exemple dans des volcans ou à proximité des failles.

Nous présentons ici, une approche différente de la tomographie utilisant uniquement des enregistrements de bruit de fond. L'idée principale est que, sur des temps longs, le bruit de fond peut-être considéré comme un champ d'onde aléatoire. Dans ce cas, la corrélation du bruit enregistré en deux points A et B est la fonction de Green du milieu entre ces deux points. La fonction de Green est l'enregistrement qu'on aurait obtenu en A si on avait placé une source en B ou vice-versa. Cette idée a été utilisée avec succès dans différents domaines, en océanographie, en acoustique, ou encore en héliosismologie.

La transposition à la sismologie n'est pas immédiate. En effet cette méthode repose sur une hypothèse forte : celle que le bruit de fond puisse être considéré comme un champ d'onde aléatoire. C'est le cas par exemple si on a des sources de bruit blanc en chaque point du milieu. Dans le cas de la terre, le bruit sismique en dessous de 1 Hz, est générée essentiellement par l'activité océanique, par un mécanisme encore mal connu.

Dans le présent article, nous présentons tout d'abord une approche théorique générale, montrant que pour un milieu quelconque la corrélation du bruit enregistrée est la fonction de Green du milieu, en supposant que les sources de bruit sont réparties de façon homogène en chaque point du milieu.

Dans la seconde partie, nous étudions les sources de bruit sismique de périodes entre 5 et 20s. En utilisant des réseaux situés en Californie, sur la côte Est des Etats-Unis, en Europe et en Tanzanie, nous montrons que le bruit est généré dans les océans. Entre 10 et 20s, il provient essentiellement de l'hémisphère nord en hiver et de l'hémisphère sud en été. Ce résultat est important car il signifie que moyenné sur un an, le bruit sismique est quasi homogène.

Dans la troisième partie, nous évaluons quelle est la précision des mesures de temps d'arrivée d'ondes de surface effectuées sur les corrélations, en utilisant 13 ans de données continues enregistrées par 3 stations californiennes. Nous montrons que les mesures sont suffisamment précises pour détecter et corriger des erreurs de temps présentes dans les données de l'ordre du dixième de seconde.

Nous présentons ensuite deux applications à la tomographie en ondes de surfaces, tout d'abord à l'échelle régionale en Europe de l'Ouest, puis à plus petite échelle sur le
volcan du Piton de la Fournaise. Enfin dans la dernière partie, nous montrons que des
ondes P peuvent également être visibles sur les corrélations.

***Mots clefs*** Sismologie, onde de surface, tomographie, bruit sismique, corrélation, control de la qualité des données

# I. INTRODUCTION

Traditional observational methods in
seismology are based on earthquake
records which result in two main
shortcomings:

1) Most techniques are based on waves emitted by earthquakes that occurred only in geologically active area, mainly plate boundaries. This results in a limited resolution in all other area where earthquakes are not present. In particular at stations far away from the source region all the high frequency information is lost due to the attenuation of the medium.

2) The occurrence of earthquakes is too low, preventing from studying “real time” change of active structures such as volcanoes or faults.

Here we explore an alternative way
of probing the earth interior using noise
records only. The main idea is to consider
seismic noise as a random source field
when averaged over long time series. In
this particular case, cross-correlation
between two stations yields the Green

function between these two points. As the
noise is mainly generated at the subsurface
by oceanic and atmospheric activity, the
surface wave part of the Green function is
mostly extracted from the cross-correlation
process.

In this paper, we first show theoretically the relationship between noise cross-correlation and the Green function [1,2]. Then we study the origin of the earth noise in the 5-20 s period band [3]. Our results indicate that the source area of seismic noise is spread over large surfaces, and can actually be used for seismic imaging. In section 4, we evaluate the accuracy of the surface wave velocity measurement performed from noise correlation [4]. Section 5 is devoted to seismic imaging. We present two examples of surface wave tomography at the regional scale in western Europe, as well as a more local scale at the “Piton de la Fournaise” volcano [5]. Finally, section 6 shows the evidence of P-waves locally extracted from noise correlation in the Parkfield network at the San Andreas Fault [6].
## II. BACKGROUND AND THEORETICAL APPROACH

The Green function of a medium between two points A and B represents the record we would get at A if we put a source at B. The use of random noise to reconstruct the Green function has already been applied successfully in various fields of physics such as helioseismology [7,8], acoustic [9] or oceanography [12].

In the case of a completely random wavefield, the cross-correlation of the wavefield recorded between two points converges to the complete Green function of the medium, a property that is reminiscent of the fluctuation-dissipation theorem [10]. By complete Green function we mean it includes all reflection, scattering and propagation modes. To demonstrate this result and to define more precisely under which assumption it is valid, various experimental, numerical and theoretical approaches have been developed. Historically speaking, helioseismology was the first field where, ambient noise cross-correlation performed from recordings of the sun surface random motion was used to produce coherent images of the sun subsurface [7,8]. More recently, a seminal paper was published in 2001 by Weaver and Lobkis [9] that showed how diffuse thermal noise recorded and cross-correlated at two

transducers fasten to one face of a duralumin sample provided the complete Green's function between these two points. Weaver theoretically interpreted this result using the equipartition of the modes excited in the duralumin sample. Derode et al [11] proposed to interpret the Green's function reconstruction in terms of a Time-Reversal analogy and showed that correlation of multiply scattered waves could be used for passive imaging in acoustic [9,14,15]. Similar results were then obtained at a much larger scale in shallow underwater acoustics where both direct and reflected wavefronts were retrieved from ambient noise cross correlation [12]. At the same time, the cross correlation process was investigated from the perspective of multiple uncorrelated sources that were considered as noise sources. By summing the contribution of all sources to the correlation, it has been shown numerically that the correlation contains the causal and acausal Green function of the medium [13]. Following this approach, the precedent result was generalized to the case where randomization is not produced by the distribution of the sources, but is provided by multiple scattering that takes place in heterogeneous medium [15].
Experimental evidences demonstrated the feasibility of passive imaging in acoustics [9,15,16] and oceanography [17]. In seismology Campillo and Paul [18] used this idea to retrieve the Green's function between two seismic stations from a collection of earthquake.

Multiple theoretical approaches quickly followed the first successful experimental demonstrations. The convergence of the noise correlation function toward the Green function in an unbounded medium is now classically interpreted through the stationary phase theorem [19,20]. The use of Green's theorem [15] or the fluctuation-dissipation approach [21] yields very rigorous theoretical approaches to interpret experimental results. Even the case of non-reciprocal (presence of a flow, for example) or inelastic media have now been theoretically investigated [13,21].

In seismology, Aki [22] proposed a long time ago to use seismic noise to retrieve the dispersion properties of the subsoil. For the problem of elastic waves, it has been theoretically showed that the convergence of noise correlation to the Green function was bonded by the equipartition condition of the different components of the elastic field [23,24]. In other words, the emergence of the Green function is effective after a sufficient self-averaging process that is provided by

random spatial distribution of the noise sources when considering long time series as well as the scattering [26,27]. Shapiro and Campillo [25] reconstructed the surface wave part of the Green function by correlating seismic noise at stations separated by distance of hundreds to thousands of kilometers, and measured their dispersion curves at periods ranging from 5 to about 150 seconds. This method led to the first application of passive seismic imaging in California [28,29] with a much greater spatial accuracy than classical active techniques. Note also that Larose et al [30] used noise cross-correlation at small distances on the moon.

In the following, we summarize the theoretical approaches using mathematical tools that allow a global view of the correlation problem in any damped propagation medium. Cross-correlation of noise recorded at two distant station A and B yields the Green function, assuming only that the noise is a white noise distributed everywhere in the medium.

We consider any medium X, homogeneous or not, where the wave propagation equation is controlled by a damped equation that can be written as:

$$ \frac{\partial^2 u}{\partial t^2} + 2a \frac{\partial u}{\partial t} - Lu = f \quad (1) $$

Here $a > 0$ is a constant that corresponds to the attenuation of the medium, $f(\vec{r}, t)$ is the source field (i.e. the noise field) and
$u(\vec{r}, t)$ denotes the displacement field. $L$ can be the Laplace operator or any other negative self-adjoint elliptic differential operators. If we choose $L = c^2(\vec{r})\Delta$, we recognize the usual wave equation.

The Green function $G(t, \vec{r}_1, \vec{r}_2)$ is the solution of the wave equation (1) if the source function $f$ is a Dirac. It represents the displacement that would be recorded at $\vec{r}_1$ if a source $\delta(\vec{r}, t)$ acts at $\vec{r}_2$. In the case of a noise field acting in medium $X$ between time $0$ and $\infty$, the resulting displacement at $\vec{r}_1$ is:

$$u(\vec{r}_1, t) = \int_0^\infty ds \int_X G(s, \vec{r}_1, \vec{r}) f(t-s, \vec{r}) d\vec{r} \quad (2)$$

We define $P = \sqrt{-L-a^2}$, and we note by $\|X\|$, the integral kernel of the operator $X$. By solving the differential equation (1), we obtain:

$$G(t, \vec{r}_1, \vec{r}_2) = Y(t)e^{-at} \left\| \frac{\sin tP}{P} \right\| (\vec{r}_1, \vec{r}_2) \quad (3)$$

with $Y$ the Heaviside function.

We now assume that $f$ is a white noise distributed everywhere in the medium $X$, acting at anytime $t$. In the frequency domain, a white noise contains all the frequency with a random phase. In the time domain, a white noise is a random signal such that the displacement field at the same point $\vec{r}$ at two different times are completely uncorrelated as well as the wavefield measured at the same time $t$ at two different points $\vec{r}_1$ and $\vec{r}_2$. We define

the time-domain normalized cross correlation process as:

$$C_{1,2}(\tau) = \int_0^T u(\vec{r}_1, \tau) u(\vec{r}_2, t+\tau) dt / \sqrt{\int_0^T u^2(\vec{r}_1, t) dt \int_0^T u^2(\vec{r}_2, t) dt} \quad (4)$$

The denominator in (4) is a normalization factor that practically mitigates local station effects and helps in determining the actual coherence of the noise field between the stations. Indeed, the auto-correlation $C_{1,2}(\tau)$ has a maximum of 1 whereas noise signals strongly decorrelated (i.e. when stations are too far away or when the incoming noise is incoherent) do not exhibit any peak in the normalized noise correlation function but a random temporal signal with a variance of $1/(2T\Delta\omega) \ll 1$, when dealing with frequency-limited noise sources of bandwidth $\Delta\omega$ [31, 32].

In case of a white noise, the law of large number allows us to replace the large $T$ limit by an ensemble average. We then get the following explicit expression for the correlation (without any normalization) between the wavefields recorded at $\vec{r}_1$ and $\vec{r}_2$:

$$C_{1,2}(\tau) = \frac{\sigma^2 e^{-a|\tau|}}{4a} \left[ (P^2 + a^2)^{-1} \left( \cos \tau P + a \frac{\sin |\tau| P}{P} \right) \right] (\vec{r}_1, \vec{r}_2) \quad (5)$$

It is interesting to consider the small $a$ behavior of the Eq(5):

$$C_{1,2}(\tau) \approx \frac{\sigma^2 \cos \tau \sqrt{-L}}{4a(L - L)} \quad (6)$$

The time derivative of Eq.(6) is expressed in terms of the Green function, giving the
more familiar expression:

$$ \frac{d}{d\tau} C_{1,2}(\tau) = \frac{\sigma^2 e^{-d|\tau|}}{4a} [G(\tau, \vec{r}_1, \vec{r}_2) - G(-\tau, \vec{r}_1, \vec{r}_2)] \quad (7) $$

This means that for any medium, the time-derivative of the cross-correlation computed between the wavefields

recorded at two stations $\vec{r}_1$ and $\vec{r}_2$ is the Green function of the medium, provided that the damping coefficient is small enough and that noise sources behave as white noise acting everywhere in the medium.

### III. SEISMIC NOISE ORIGIN IN THE 5-20 SECONDS PERIOD BAND

We have seen in the previous section that the time cross-correlation function of random wavefields computed between a pair of distant stations contains, at least partially, the actual Green function between the two stations. This provides us with a possibility to retrieve the propagation properties of deterministic seismic waves along long paths by analyzing seismic ambient noise only (in the following, we use the term noise for the microseism which actually have no relation with instrumental noise). Understanding to what extents seismic noise can be used for seismic imaging requires to know the spatial and temporal distribution as well as the nature of its source.

Ambient seismic noise is mostly made of surface waves [33, 34]. Noise sources are likely close to the Earth's surface and are believed to be loads caused by pressure perturbations in the atmosphere and the ocean. At relatively

short periods (< 20 s) the noise spectrum is dominated by two peaks around 7 s and 14 s, i.e. the primary and the secondary microseismic peaks. They are believed to be related to the interaction of the sea waves with the coast [35]. The primary microseism periods are similar to the main swell (10-20 s), while the strongest secondary microseism peak in the noise spectrum originates from the non-linear interaction between direct and reflected swell waves that results in half period (5-10 s) pressure variations [36]. This interaction leads to pressure variations at the sea bottom that do not exhibit the rapid exponential decay with depth expected for primary gravity waves.

Seismic ambient noise records are dominated by a few energetic sources [33] and cannot be considered as a random field. However the contribution of the most energetic noise sources can be reduced by disregarding completely the amplitude considering 1-bit signal only
[18,25,37]. Alternative whitening
techniques were used by Sabra et al [29].

This “normalized background
seismic noise” is assumed to be a random
field when considering long time series

[25,32]. Another process contributing to the noise randomization is the scattering of seismic waves on heterogeneities within the Earth that is significantly strong at periods less than 40 s.

<figure><img src="image_82.png" /></figure>

<figure><img src="image_83.png" /></figure>
<figure><img src="image_84.png" /><figcaption>**Figure 1:** Average normalized amplitude of cross correlations of noise versus azimuth, for various network a) During winter between 5-10 s, b) Summer 5-10s, c) Winter 10-20s, d) Summer 10-20s.</figcaption></figure>

Optimization of seismic imaging based on noise correlation requires to understand the origin of the seismic noise and the spatial and temporal distribution of

its sources [38,39]. Here we study the origin of the background seismic noise averaged over long time to determine the seasonally-averaged azimuthal distribution
of the Normalized Background Energy Flow (NBEF) through several networks located in California, on the East coast of the USA, in Western Europe and in Tanzania. The term NBEF is used to remind that we don't study the energy flow

of the “real” seismic noise, but of the noise whose amplitude has been normalized by one-bit coding which is the relevant analysis in the context of the Green function reconstruction by noise cross-correlation.

### III.1 Using asymmetry of cross-correlation to study the origin of the seismic noise

In the case of a spatially homogeneous distribution of noise source, the noise cross-correlation function (NCF) between

two stations 1 and 2 is expected to be symmetric in time and amplitude (see Eq. (7)).

<figure><img src="image_85.png" /></figure>
<figure><img src="image_86.png" /><figcaption>**Figure 2:** Comparison between seasonal variations of the location of seismic noise sources and significant wave height. (a) and (c): geographical distribution of the mean back-projected normalized amplitude of the 10-20s noise cross-correlations (see text for detailed explanation of the back-projection) during the winter and the summer. (b) and (d): Global distribution of the square of wave height measured by TOPEX/Poseidon during the winter and the summer.</figcaption></figure>

The positive time of the NCF corresponds to the causal Green function of the medium (ie wave going from 1 to 2), and the negative time to its anti-causal counterpart (i.e. waves going from 2 to 1). The amplitude of the Green function reconstructed by the NCF directly depends on the energy flux of the waves traveling

from one station to the other [40,41]. This means that in the case of an isotropic distribution of noise sources, the energy flux between two stations is the same in both directions and the resulting NCF between these stations is symmetric. On the other hand, if the spatial source density is asymmetric, the amounts of energy
propagating in both directions are
different.

The asymmetry of the NCF
computed between several pairs of stations
of a network can be used to measure the
main direction of the energy flux across an

array of seismic stations. Performing such
measurements at different arrays allows us
to determine the location of the coherent
seismic noise sources in the bandwidth of
interest.

### III.2 Data processing and results

We analyze the noise origin at several
networks located in the East coast of USA,
in Western Europe, in California and
Tanzania. Each of these networks contains
about 20 broadband seismic recorders
separated by distances ranging from one
hundreds to a few hundreds of kilometers.

At each of these networks, one year
(2003) of continuous vertical records has
been analyzed. The data are processed in
the following way:

* All noise records are corrected
  from the instrumental response.

* The records are then filtered in the
two period bands we are studying: 5-10 s
and 10-20 s.

* To reduce the contribution of the
most energetic source of noise, we
completely disregard the amplitude of the
records by considering 1-bit signal only.

* The NCF between each pair of
station within each network are computed.

* We estimate the normalized
amplitude of the NCF at arrival times
corresponding to Rayleigh wave
propagation time with respect to the

azimuth between each station pair.

Figure 1 shows the average
normalized amplitude of the NCF versus
azimuth for various networks during
winter and summer and for the two period
bands 5-10 s and 10-20 s. The NBEF of
the 5-10 s noise remains stable over the
whole year at all networks. This
observation is compatible with the idea
that the generation of the secondary
microseism is mostly controlled by the
bathymetry and the coastlines geometry.
The noise in the 5-10 s period is
dominated by the activity in a region close
to the station since absorption limits the
contribution of distant source zones. It
explains the stability of the distribution
during the year.

The behavior of the 10-20 s noise is
very different (Figure 1). Indeed the effect
of attenuation is weaker for longer periods
and sources at the global earth scale are
contributing. In this period band, the
regions where the noise is originating are
different during summer and winter. A
first order observation from Figure 1 is
that the major part of the noise seems to
come from the North during winter and

from the South during the summer.

### III.3 Comparison of the origin of the 10-20 s noise map with ocean wave height map

In spite of their limited number and
coverage, the considered networks
illuminate a large part of the Earth surface.
In the following, we use them for a rough
identification of the source regions that
generate the 10-20 s noise consistently
recorded by the four networks.

We use a simple back projection
procedure. For every considered station
pair of a given network, we trace the great
circle connecting the stations and divide it
in two equal parts corresponding to two
waves propagating in opposite directions.
Amplitudes measured from positive and
negative parts of the NCF computed
between these stations are attributed to the
corresponding half great circle.

We then define a 2° x 2° grid on the
Earth surface and take the mean value of
the normalized correlation amplitude
associated with all great circles crossing
every cell. The results for the four
networks in the 10-20 s band are presented
in Figure 2. During the northern
hemisphere winter, the observed average
noise is apparently dominated by a source
zone in the northern Atlantic ocean. Two
weaker sources are visible in the northern
Pacific, one close to Alaska and the other

close to Japan. During the northern
hemisphere summer, most of the noise
originates from the southern Indian ocean
and the southern Pacific. It has to be noted
that the resolution of these patterns is
limited by the small number of networks.
An example is the artifact in the Pacific
Ocean which is due to the same great
circle shared by the California, eastern US
and Western Europe networks.

These results clearly show that,
contrary to the secondary microseism, the
10-20 s noise is not generated locally but
is excited by sources acting at the global
scale and having a clear seasonal
variability. An attempt to locate these
sources suggests they have to be searched
in the central parts of the oceans, although
we cannot exclude that part is generated on
the coast or in shallow water areas with
our limited resolution. Specifically, with
the form of normalization and averaging
used here, we do not contradict the
observation that the strong storm arrivals
along the coasts are associated with high
amplitude microseisms.

The hypothesis that the average
background noise energy originates in the
zones of storm activity is further
confirmed by strong similarities between
the maps of the apparent sources of the 10-
20 s seismic noise and the global wave
height maps obtained from the radar
altimeter data collected by the satellite
TOPEX-POSEIDON (Figure 2). Both
maps show clear seasonal variations with

maxima located in deep oceans in the
southern and the northern hemispheres
during the summer and the winter,
respectively. A similar seasonal variability
was recently reported for the very long
period (> 150 s) Earth hum [42].

IV. TRAVEL TIME MEASUREMENTS FROM NOISE
CORRELATION:

STABILITY AND DETECTION OF INSTRUMENTAL ERRORS
USING THE PRINCIPLE OF TIME SYMMETRY

We have seen in the previous section that
although the “normalized background
seismic noise” is not fully isotropic, it
seems to be enough randomized to enable
us to retrieve Green function by cross-
correlation. Our goal here is to assess the
accuracy of the Green function
reconstructed by cross-correlation of
ambient seismic noise and to determine to
what extents the noise-based travel time
measurements can be used to detect and to
quantify station instrumental errors and to
monitor changes in the physical properties
of the medium.

To that goal, we analyze the
temporal fluctuations of surface waves
travel times measured on the cross-
correlation of seismic noise between three
Californian stations GSC, PAS and PFO.

As previously mentioned, the surface wave
arrival time should be identical on the
positive and negative part of the NCF.
However, three main factors may lead to
fluctuations of the NCF between two
stations:

1) A physical change in the medium would result in either a faster or slower travel time measured in both positive and negative cross-correlation time.

2) A clock error in one of the two stations would produce a time shift of the whole cross-correlation resulting in a larger travel time in the positive time and a smaller apparent travel time in the negative time or vice-versa. A change in the phase response of one of the
sensors would have the same
effect.

3) A change in the spatial distribution of the source of the noise should affect the positive and negative correlation time independently, since the positive and negative time are sensitive to noise sources

located in different regions.

Studying simultaneously how the
apparent travel time of the reconstructed
surface waves change with time in the
positive and negative part of the cross-
correlation allows us to evaluate these
three main causes of fluctuations.

<figure><img src="image_87.png" /><figcaption>**Figure 3:** (a) Variation of the surface waves travel time measured with one month stack for GSC-PAS in the 5-10 s period band. (b) Measurement of the variation of the surface waves travel time in the 5-10 period band expressed in second, after correction of the instrumental errors.</figcaption></figure>

IV.1 Data and measurements

We analyze the fluctuations of the
apparent travel time of Rayleigh and
Loves waves reconstructed by cross
correlating seismic ambient noise, by
using 11 years (1991-1996 and 1999-2005)
of continuous records on the three
components at three Californian stations:
GSC, PFO and PAS. We first define a

"Reference Green Function" (RGF) for
every path in two period bands 5-10 s and
10-20 s, and for all components (ZZ, TT,
RR, ZR and RZ) by computing the one-bit
cross-correlation functions of 6 years
(1999-2005) of continuous seismic noise
records.

Horizontal components records
were rotated to radial and transverse
directions assuming the propagation along
the inter-station great circle. Then we
computed cross-correlations for each
thirty-day-period during 11 years (1991-
1996 and 1999-2005). For every time
window, we compare the arrival time of

the surface wave of the current month with
the RGF for both the causal and anti-
causal NCF. The travel time variations are
measured from the phase of the cross-
spectrum computed between the RGF and
the NCF.

<figure><img src="image_88.png" /><figcaption>**Figure 4:** Paths where 16 s Rayleigh wave group velocity measurements were obtained from cross correlations of ambient seismic noise. White triangle shows stations used for this study.</figcaption></figure>

IV.2 Stability of the travel time measurements between 1999 and 2005

We first consider period 1999-2005, for
which we did not detect any obvious
instrumental errors. Figure 3a shows the
apparent travel time variations of surface
waves for paths GSC-PAS in the 5 to 10 s

period band, averaged over all components
of the correlation tensor, that is including
both Love and Rayleigh waves. Results
from the positive time of the correlation
shown with a blue line correspond to
waves going from GSC to PAS. The red line shows results from the negative part of the NCF corresponding to waves traveling from PAS to GSC.

The path GSC-PAS is nearly perpendicular to the Californian coastline. The negative time of the GSC-PAS NCF is sensitive to sources located in the Pacific Ocean [3] and exhibits smaller fluctuations

than the positive time which is sensitive to sources in the Atlantic Ocean. In the 5-10 s period band the fluctuations on the negative time of GSC-PAS never exceed 0.05 s and has a variance of only 0.0007 s which represents respectively 0.08% and 0.0012% of the Rayleigh waves travel time. This shows that the reconstruction of the Green function is extremely robust.

<figure><img src="image_89.png" /><figcaption>**Figure 5:** Rayleigh wave group velocity map at 16 s period constructed by cross-correlating one year of seismic noise records between each pair of station.</figcaption></figure>
<figure><img src="image_90.png" /></figure>

**Figure 6:** (a) Map of the Piton de la Fournaise volcano. Seismic stations are represented as inverted triangles. The gray zone indicates the limits of the rift zone. The thin dashed rectangle corresponds to the limits of the presented tomographic images. Geographic coordinates are Gauss-Laborde kilometric coordinates (Transverse Mercator). Contour lines of elevation are spaced every 100m. (b) Two hours of ambient seismic noise (ANR). (c) Causal and acausal reconstructed Rayleigh waves (positive and negative times, dominant period=4 s) between station RMR (not shown on the map) and the rest of the network. The trace envelopes are represented as thin gray curves.
IV.3 Evidence of instrumental errors between 1992 and 1996

We now consider the period 1992-1996, for which instrumental errors were detected. Variations of the arrival time for the path GSC-PAS for the 5-10 s period bands during years 1991-1996 are displayed in Figure 3a. Averaged time delays measured from the positive and negative times of the NCF are clearly correlated:

* In 1992, the time delay exhibits a peak reaching 2 s, that can be observed both for positive and negative correlation time with the same polarity.

* A positive time shift is observed during all the considered period on the positive and negative time.

* There is a small positive trend between 1992-1996, which means that the apparent arrival time becomes larger on the positive time and smaller in the negative time.

These observations cannot be related to changes of the physical properties of the crust that would result in fluctuations of opposite polarity for the positive and negative correlation time (i.e. shorter or larger travel time for both positive and negative correlation time). It is also unlikely that such almost perfectly even variations are due to the variations in the distribution of the microseism sources because the positive and negative NCF are not sensitive to sources located in the same region. Moreover, all these features are observed on the two period bands 5-10 s and 10-20 s while noise sources are not expected to be the same (see section 2 and [3]). In conclusion, these even and period-independent time fluctuations strongly suggest a clock drift or instrumental errors at one of the two stations.

IV.4 Correction of the instrumental errors

Apparent travel time fluctuations due to instrumental errors are isolated by averaging the measurement performed in the positive and negative correlation time in the two period bands 5-10s and 10-20s. We distinguish two kind of robust features:

* Rapid fluctuations with a high magnitude: The main peak of 1992, a minor peak in 1994, the offset of the curve

between 1991-1996 and 1999-2005.

* Slowly evolving errors: the average values of the measured instrumental error change slightly with time.

By averaging the rapid and slowly evolving instrumental errors measured on the two period bands 5-10 s and 10-20 s we obtain an evaluation of the instrumental errors that can be used to correct the data.
<figure><img src="image_91.png" /></figure>

**Figure 7:** 3D S-wave velocity model. We show 6 horizontal slices extracted from the 3D model at different depths. Average S-wave velocity is shown in white boxes on the bottom of corresponding slices. Black dashed line at depth -0.5 km shows the limits of the rift zone. We also plot a 3D view of the model. The 3D blue patch delimits the iso-velocity perturbation surface corresponding to 2.5% of velocity perturbation. The encircled numbers show the horizontal positions of the vertical velocity profiles shown in Figure 6c.
Figure 3b presents the travel time variations re-measured for the GSC-PAS path between 5 and 10 s after removing the estimated instrumental errors. As expected, the peaks of 1992 and 1994 have completely disappeared, as well as the time shift between the 1991-1996 and 1999-2005 period and the small trend between 1991 and 1996 (Figure 3b).

Travel time fluctuations after the correction of instrumental errors have a variance of 0.012 s and 0.0022s on the positive and negative time, respectively. This represents 0.02% and 0.003% of the travel time of Rayleigh waves. This shows that the Green Functions reconstructed by

correlating ambient seismic noise are accurate enough to detect and to correct small instrumental errors and to measure inter-station travel times with a precision of a few tenths of second for stations 200 km apart.

This measure of instrumental error does not require any assumption on the wave velocity and relies only on the fundamental properties of time symmetry. It could be used routinely to detect instrumentals problems. Finally note that similar synchronization corrections have been successfully performed between underwater acoustic arrays from ocean ambient noise [43].

<figure><img src="image_92.png" /><figcaption>**Figure 8:** Topographic map of the Parkfield area (an 11-km large square) showing stations (triangles) and SAF (blue).</figcaption></figure>

# V. APPLICATION TO IMAGING

We have seen in the section 3 that the spatial distribution of the background

seismic noise in the 5-20 s period band covers a large surface when integrated
over long time. This allows us to retrieve
the Green function between two stations
by correlating background seismic noise
records. The emerging signal of the NCF
is dominated by surface waves, since the
background seismic noise is mainly made
of surface waves. The reconstructed Green
function are stable over time (see section
4) and robust enough to measure surface
wave propagation time with a precision of
a few tenth of a second, independently of
the azimuth of the considered station pair

path.

Passive imaging from seismic noise
and Rayleigh wave group velocities has
been first used by Shapiro et al and Sabra
et al [28, 29] who provided images of the
Californian crust. Here, we present two
examples of seismic noise processing to
produce Rayleigh waves group velocity
map at the regional scale in Western
Europe and to compute 3D S-velocity
model at the Piton de la Fournaise volcano
[5].

<figure><img src="image_93.png" /></figure>

**Figure 9:** Angular-speed distribution of pre-processed incoming noise on the Parkfield network averaged over one month. Plane wave beamforming is summed incoherently over 100 frequencies from (a) 0.2 to 0.5 Hz and (b) 0.6 to 1.3 Hz. The x-axis corresponds to noise directivity θ (North is 0°), the y-axis to speed c. Units are arbitrary.
<figure><img src="image_94.png" /><figcaption>**Figure 10:** Range-time representation of the Z-Z component of the noise correlation tensor averaged over one month in two frequency bands (a) [0.1-1.3 Hz], and (b) [0.6-1.3 Hz]. Each plot has been normalized by its own maximum.</figcaption></figure>

## V.1 Surface wave tomography of Europe

One year of continuous records from October 2004 to October 2005 at 190 3-component broadband stations spreads over Western Europe has been used with a particular emphasis in the Alps (Figure 4). For each station pairs, we computed the one-bit cross-correlation components that include the Rayleigh wave part of the Green function (ZZ, RR,ZR and RZ). The

spectrum of the noise records were whitened in the 5-150 s period band prior to the correlations.

Rayleigh waves dispersion curves are evaluated from the emerging Green function using frequency-time analysis [44,45] for the 17,000 inter-station paths. For each path we get eight evaluations of the dispersion curves since we consider
four components of the correlation tensor (ZZ, RR, RZ and ZR) and both the positive and the negative part of the NCF.

We reject all waveforms (1) with signal to noise ratio lower than seven, (2) with group velocities measured on the positive and negative correlation time differing by more than 5 percents, and (3) with path shorter than two wavelengths.

This result in about 3000 paths over the initial 17000 paths inter-stations

paths at 16 s (Figure 4). We then apply a tomographic inversion following [46] to this data set to obtain a group velocity map on $100*100 = 10,000$ 25x25km cells across Europe (Figure 5). Several geological features can be seen on the resulting group velocity map, in particular sedimentary basin corresponding to low velocity anomalies such as the Po basin (North Italy), the Rhone Basin (South-East France) and the Pannonia basin.

## V.2 3-D S-wave tomography of the Piton de la Fournaise volcano

Continuous seismic noise, recorded during 18 months (Jul. 1999 to Dec 2000) from 21 short period stations were collected by the Observatoire Volcanologique du Piton de la Fournaise (Figure 6a).

An example of a noise record at one of the stations (ANR) is shown in Figure 6b. All noise records are first bandpass filtered between 1 to 5 s and their spectral amplitude are whitened in order to avoid strong dominant spectral peaks in the background noise.

For each available station pair, the one-bit NCF is computed day per day. Some of the NCFs are rejected upon a signal to noise ratio criterion. The remaining traces are stacked over 18 months. For each path, group velocity dispersion curve are estimated using an automatic Frequency-Time Analysis (FTAN, [44,45]). We manually select

dispersion curves according to group velocity limits and for the station-to-station distance shorter than one wavelength. We finally obtain 75 reliable dispersion curves from which group velocities are extracted for periods equal to 2, 2.5, 3, 3.5, 4 and 4.5 s.

2-D Rayleigh wave group velocity map are obtained from tomographic inversion of the arrival-time measurements at each period using the algorithm described by Barmin et al [46]. Our 2-D model involves $22 \times 28 = 616$ 1 x 1 km cells. Because of the sparse ray coverage and the low resolution of the data set, we chose to apply a strong smoothing to the tomographic inversion which results in a ~4-km spatial resolution. The inversion results are thus robust and show a moderate variance reduction varying from 38 to 18 % with increasing periods (from 2
to 4.5 s).

Dispersion curves are then extracted from the Rayleigh-wave group-velocity maps for every model cell. We fit these curves by polynomial functions in a least-squares sense and invert them using a Monte-Carlo algorithm [47], the synthetic dispersion curve being calculated using a method by Herrmann [48]. We thus obtain a Vs-versus-depth velocity profile for each cell. We present six horizontal slices as well as a 3-D view of the 3-D smoothed model on Figure 7. The results clearly show the presence of a high velocity anomaly which moves westward with depth (+1.3 to -1.1 km above sea level).

This structure is surrounded by a low-velocity ring interpreted as effusive products associated to the construction of the Piton de la Fournaise volcano on the flank of the older Piton des Neiges volcano. This high velocity anomaly has also been detected by a previous earthquake and active P-wave tomography on the Piton de la Fournaise volcano [49]. Recent works also imaged the presence of a high velocity chimney on different volcanoes [50-54]. We interpret this anomaly as a solidified intrusive magma body. The high velocity anomaly is also well correlated with the rift zone at sea level (z=-0.5 km).

# VI. P-WAVES EXTRACTION FROM SEISMIC NOISE CROSS-CORRELATION

A main issue in the convergence of the correlation process to the transfer function is the influence of variations in the temporal and spatial distribution of the noise sources. From the temporal point of view, the noise spectrum defines the frequency bandwidth over which the impulse response can be retrieved. When receivers are widely separated, the coherent propagating noise must have sufficient amplitude to be recorded on both receivers despite geometrical spreading and attenuation. This explains why the slowly-attenuated Rayleigh waves have

dominated the impulse response obtained so far from correlation of seismic noise.

In this section, we present results from the correlation function of seismic noise recordings among pairs of stations in the dense Parkfield network, California. When performed on many station pairs at short ranges, the noise correlation function (NCF) demonstrates the presence of both a P-wave and a Rayleigh wave in the NCF.

We processed data recorded on the dense temporary seismic network installed in the Parkfield area between July 2001 and October 2002. One month of seismic
noise recordings were cross-correlated between each pair of 30 broadband 3-component seismic stations located in an 11-km square (Fig 8). This network has been extensively used to monitor and image the San Andreas Fault (SAF) using both man-made explosions and earthquakes as part of the San Andreas Fault Observatory at depth (SAFOD) project. Inversion results have confirmed the spatial heterogeneity of P-wave velocity across the Fault up to 6 km in depth [55-58].

At first, frequency-incoherent beamforming is performed using the $N=30$ stations of the network to determine the average velocity $c$ and direction $\theta_0$ of the seismic noise (Fig. 9). Beamforming is performed in two bandwidths of interest [0.2-0.5 Hz and 0.6-1.3 Hz] on data segments of one day of seismic noise, recorded in February 2002, as:

$$B(\theta, c) = \frac{1}{\Delta\omega} \int_{\omega_c - \Delta\omega/2}^{\omega_c + \Delta\omega/2} \sum_{i=1}^{N} \tilde{S}_i(\omega) \exp \left[ i \frac{\omega}{c} (x_i \sin\theta + y_i \cos\theta) \right] d\omega \quad (8)$$

where $\omega_c$ is the central noise frequency and $\Delta\omega$ the frequency bandwidth, $\tilde{S}_i(\omega)$ is the complex Fourier component at frequency $\omega$ of the noise signal $S_i(t)$ on the $i^{th}$ seismic station ($i \in [1, N]$), and $(x_i, y_i)$ are the longitude/latitude coordinates of station $i$. Working with a dense seismic network having a small coverage area

allows to determine an average apparent incoming velocity for this zone using plane wave beamforming. In the two frequency bandwidth, the noise field clearly originates from the Pacific Ocean with an incident direction $\theta_0=55^\circ$ on the Parkfield network. On the other hand, the low frequency beamformer exhibits an apparent velocity of 2.8 km/s compatible with a Rayleigh wave (Fig. 9a) while the upper frequency beamformer shows an apparent velocity of 5 km/s (Fig. 9b).

Since Fig. 9 reveals a strong directivity in the seismic noise, only the station pairs aligned with the noise main direction $\theta_0$ are chosen to compute the point-to-point NCF. Restricting the cross-correlations to those pairs ensures that 1) the travel time of the main peak of the NCF is not biased and corresponds to the actual travel time of the Green's function between the stations [19,20] and 2) the signal-to-noise ratio (SNR) of the NCF is maximized. Practically speaking, 145 station pairs are selected in the Parkfield network whose angles $\theta_{ij}$ are included in a directivity angle $\Delta\theta = |\theta_{ij} - \theta_0| \le \sqrt{c/R_{ij}\omega_c}$ dependent on the distance $R_{ij}$ between stations and the frequency characteristics of the seismic noise field [12].

The NCF is computed for each selected station pairs as in Eq.(4) and coherently averaged over 30 days to
further increase the SNR. Figure 10a is a
display of the NCFs obtained for the 145
selected pairs sorted by ascending range
$R_{ij}$. A propagating wavefront clearly
appears at high frequency (Fig. 10b) which
corresponds to a ~5 km/s velocity wave. A
polarization study between the Z-R and Z-
Z components of the correlation tensor
confirmed the P-wave nature of this
wavefront [6].

Going back to Fig. 9b, we note that
the apparent velocity of the P-wave
corresponds to the P-velocity at the turning
point. Recent inversions of the P-wave

velocity profile reveal a strong velocity
gradient at the surface, the 5km/s speed
being reached at no more than 1.5 km in
depth on the west side of the SAF. This
confirms that the noise sources that excite
P-waves are local and can't be confused
with deep incident waves that would hit
the seismic array with a much higher
apparent velocity. One realistic hypothesis
is that P-waves are locally generated by
conversion of incident Rayleigh waves
coming from the Pacific by local
heterogeneities at the Earth's surface or
subsurface [59].

VII. REFERENCES

[1] Y. Colin de Verdiere (2006a), http://fr.arxiv.org/abs/math-ph/0610043/

[2] Y. Colin de Verdiere (2006b), http://fr.arxiv.org/abs/math-ph/0610044/

[3] Stehly et al (2006), JGR Vol 111.

[4] Stehly et al (2007), submitted.

[5] Brenguier et al (2007), GRL 34.

[6] Roux et al (2005), GRL 32, L19303.

[7] Duvall et al (1993), *Nature*, 362, 430-432

[8] Gilles, et al (1997), *Nature*, 390, 63-64.

[9] Weaver, R. L., and O. I. Lobkis (2001),
*Physical Rev. Lett.*, 87(13).

[10] Kubo, R. (1966), *Rep. Prog. Phys.* 29,
255-284.

[11] Derode et al (2003), *J. Acoust. Soc. Am.*,
113, 2973-2976.

[12] Roux et al (2004), *J. Acoust. Soc. Am.*,
116, 1995-2003.

[13] Wapenaar, K. (2004), *Phys. Rev. Lett.* 93,

254301.

[14] Derode et al (2003), *Appl. Phys. Lett.*, 83,
3054-3056.

[15] Lobkis, O. I and R. L Weaver (2001), *J. Acoust. Soc. Am.*, 2001a.

[16] Larose et al (2004), *Journal of Applied Physics*, 95, 8393-8399.

[17] Roux, P., and W. A. Kuperman (2003), *J. Acoust. Soc. Am.*, 116, 1995-2003.

[18] Campillo, M., and A. Paul (2003),
*Science*, 299, 547-549.

[19] Snieder, R. (2004), *Physical Review E*,
69.

[20] Roux et al (2005), *J. Acoust. Soc. Am.*,
117, 79-84.

[21] Godin, O. A (2007), submitted.

[22] Aki, K., 1957. *Bull. Earthq. Res. Inst.*, 35, 415-456

[23] Sanchez-Sesma et al (2006), GRL, 33.
[24] Sanchez-Sesma, F., and M. Campillo
(2006), *BSSA in press*.

[25] Shapiro, N., M. Campillo (2004), *GRL*, 31, doi: 10.1029/2004GL019491.

[26] Campillo, M. (2006), *Pure App. Geo-phys.*, 163, 475-502.

[27] Larose et al (2006), Geophysics, Vol 71, issue 4.

[28] Shapiro et al (2005), *Science*, 307, 1615-1618.

[29] Sabra et al (2005), *GRL*, 32, doi: 10.1029/2005GL023155.

[30] Larose et al (2005), *GRL*, 32 (16), L16201, doi: 10.1029/2005GL023518.

[31] Weaver, R.L., and O.I. Lobkis (2005), *J. Acoust. Soc. Am.*, 117, 3432-3439.

[32] Sabra et al (2005), *GRL*, 32.

[33] Friedrich et al (1998), *J. of Seismology*, 2, 47-48

[34] Ekström, G. (2001), *JGR*, 106, 26,483-26,494, doi: 10.1029/2000JB000086.

[35] Gutenberg, B. (1951), *American Meteorological Society*, pp. 1303-1311.

[36] Longuet-Higgins, M. (1950), *Philos. Trans.R.Soc.Lond.*, 243, 137-171.

[37] Derode et al (1999), *J. of Applied Physics*, 89(9), 6343-6352.

[38] Pedersen et al (2006), submitted.

[39] Schulte-Pelkum et al G3, 5, 3004, doi: 10.1029/2003GC000520.

[40] Van Tiggelen, B. (2003), *PRL*, 91(24).

[41] Paul et al (2005), *JGR*, 110, 8302, doi: 10.1029/2004JB003521.

[42] Rhie, J., and B. Romanowicz (2004), *Nature*, 431, 552-554.

[43] Sabra et al (2005), *IEEE J. Of Oceanic Engineering*, 30(2), pp 338-347.

[44] Levshin et al (1989), Seismic Surface Waves in a Laterally Inhomogeneous Earth.

[45] Ritzwoller, M., and A. L. Levshin (1998), *JGR*, 103(4839-4878).

[46] Barmin et al (2001), *Pure and appl. Geophys.*, 158, 1351-137

[47] Shapiro et al (1997), *GJI*, 128, 151-166.

[48] Herrmann R.B. (1987), St Louis University, St Louis, MO.

[49] Lankar V. (1997), PhD, *Institut de Physique du Globe de Paris*, France.

[50] Laigle et al (2000), *JGR*, 105(B9), 21.633-21.646.

[51] Tanaka et al (2002), *GRL*, 29, 59-62.

[52] Zollo et al (2002), *GJI*, 151, 566-582.

[53] Sherburn et al (2006), *GJI*, 166, 957-969.

[54] Patan et al (2003), *Science*, 313, 821-823.

[55] Ben-Zion, Y and P. Malin (1991), *Science*, 251, 1592-1594.

[56] Catchings et al, (2002), *Bull. Seismol. Soc. Am.*, 92, 2493-2503.

[57] Thurber et al (2004), *GRL*, 31, L12S02.

[58] Hole et al (2006), *GRL*, 33, L07312.

[59] Aki, K., and P. G. Richards (1980), W. H. Freeman & Co, 487-498.
VOLCANOLOGIE ET CHIMIE
DE L'INTÉRIEUR DE LA TERRE
GÉOMAGNÉTISME ET AÉRONOMIE
# New Perspectives on Mars' Crustal Magnetic Field

BENOIT LANGLAIS¹ & YOANN QUESNEL¹, ²

1 CNRS UMR 6112, Nantes Atlantique Université, Laboratoire de Planétologie et Géodynamique, Faculté des Sciences et Techniques, 2 rue de la Houssinière, 44300 Nantes cedex, France

2 Now at GFZ Potsdam, Section 2.3 Earth's Magnetic Field, Telegrafenberg, 14477 Potsdam, Germany

Benoit.Langlais@univ-nantes.fr

**Abstract** The planet Mars lacks today a planetary, dynamic magnetic field. Instead, the Mars Global Surveyor NASA spacecraft detected strong, intense, localized magnetic fields of lithospheric origin. On the average, these magnetic fields are one to two orders of magnitude larger than the terrestrial lithospheric fields at comparable altitudes. The current distribution of magnetic fields results from past magnetization and demagnetization processes. The magnetic minerals of the crust acquired a remanent magnetization while cooling under their Curie temperature in the presence of a magnetic dynamo. This dynamo subsequently shut down, and the remanent magnetization was altered by surface or deep events, such as impact craters or volcanic activity. This crude scenario can be more accurately specified, by interpreting global and local models of the current magnetic field of Mars in terms of magnetization intensities and directions. Some specific areas are studied, including the intensely magnetized Terra Sirenum, as well as the magnetic anomaly associated with Apollinaris Patera. Magnetic minerals could be of primary and/or secondary origin; this latter would imply an early hydration of a basaltic crust. A scenario, in which Mars experienced a major polar wander due to the Tharsis bulge, prior to the cessation of its dynamo, is proposed and discussed.

**Key words** Mars, Lithospheric Magnetic Field, Paleodynamo, Remanent Magnetization
# Nouvelles Perspectives sur le Champ Magnétique Crustal de Mars

**Résumé** La planète Mars n'a actuellement pas de champ magnétique dynamique global. En fait, la mission de la NASA Mars Global Surveyor a mis en évidence un champ magnétique lithosphérique très intense et très localisé. En moyenne, ces champs magnétiques sont un à deux ordres de grandeur plus forts que ce que l'on observe sur la Terre. La répartition de ce champ magnétique lithosphérique est le résultat de processus anciens d'aimantation et de désaimantation. Les minéraux magnétiques de la lithosphère ont acquis leur aimantation en présence d'une dynamo martienne, lorsque ils ont refroidi sous leur température de Curie. Cette dynamo martienne s'est ensuite arrêtée, et l'aimantation lithosphérique rémanente a été modifiée par différents phénomènes en surface ou en profondeur, comme des impacts ou des volcans. Ce scénario d'évolution très simpliste peut être affiné, grâce à l'interprétation de modèles globaux et locaux du champ magnétique de Mars, en termes de direction et d'intensité de l'aimantation. Des zones particulières sont étudiées, dont la région très aimantée de Terra Sirenum, et l'anomalie magnétique associée au volcan Apollinaris Patera. Les minéraux magnétiques peuvent avoir une origine primaire ou secondaire, ce qui implique alors une hydratation précoce de la croûte martienne. Un scénario, dans lequel l'axe de rotation de Mars s'est fortement déplacé à cause de la mise en place du complexe volcanique Tharsis, est présenté et discuté.

*Mots clefs* Mars; Champ Magnétique Lithosphérique, Paléodynamo, Aimantation Rémanente

## I. INTRODUCTION

Mars and the Earth differ in several aspects, including size, color, atmosphere, water, and topography. However, their internal structures are very alike, as it has been inferred from the mean density and the moment of inertia, complemented by meteorite analyses. The martian crust is on average 50-km thick, as it has been

inferred from geoid and topography analyses [Wieczorek and Zuber, 2004]. Beneath this crust, there is a silicate mantle, between 1300 and 1800 km thick. At the center, there is a metallic, iron-rich core. The exact state of the core remains unknown, even if recent observations concluded that it has to be at least partially
liquid [Yoder *et al.*, 2003].

The existence of an active magnetic dynamo on Mars was debated until 1996, when the NASA mission Mars Global Surveyor (MGS) was launched. This spacecraft was designed to study the past and present geological and geophysical evolution of Mars. Together with Mars Pathfinder, MGS was the first successful mission to Mars since the Viking 1 and 2 missions in 1975. MGS was also the first mission devoted to the study of the magnetic field of Mars, thanks to the MAGnetometer and Electron Reflectometer (MAG/ER) package. Prior to MGS, very few spacecrafts had carried magnetometers, and only Mariner 4 [Smith *et al.*, 1965], and Mars 2, 3 and 5 [Dolginov *et al.*, 1973] made some measurements at large distances of the planet. Conclusions were poor about the existence of a dynamic (i.e. reversing) or fossil magnetic field. Experiments on Viking were also inconclusive [Hargraves *et al.*, 1977], as well as the results from the Phobos 2 spacecraft [Riedler *et al.*, 1989]. The measurements made on board MGS were thus eagerly awaited.

MGS was inserted around Mars on Sept. 2nd 1997 on a very elliptical orbit. Its low periapsis made it possible to measure the magnetic field at very low altitude (as low as 90 km) for a short period of time.

Such measurements were performed as early as during the 4th orbit. These first measurements led to the conclusion that Mars lacks today a global, planet-scale, magnetic field [Acuña *et al.*, 1998].

The magnetic field of Mars is characterized by localized, sometimes very intense, magnetic field anomalies. The term ‘anomalies’ here refers to anomalous with respect to the global (zero) field of Mars, by analogy with the terrestrial lithospheric magnetic field anomalies. This came as a major surprise, as it was not envisaged that a remanent magnetic field could be that intense (MGS measured fields as large as 1500 nT at 90 km altitude, above the very cratered Terra Cimmeria and Terra Sirenum areas). Such remanent magnetic fields are the relics of an ancient, now extinct, dynamic magnetic field, modified by several alteration processes (volcanoes, impact craters, ...). Any interpretation consequently relies on a careful analysis of magnetic measurements, but also on the understanding of the different processes that affected the remanent magnetization during Mars’ history.

In the following, we discuss only the internal magnetic field of Mars. We first present the lithospheric magnetic field of Mars as it is modeled from MGS magnetic measurements. We then focus on
local anomalies, using different modelings
approaches. These local analyses allow
some interesting properties of the ancient

magnetic field of Mars to be investigated,
which bring new and important constraints
on Mars' formation and evolution.

<figure><img src="image_95.png" /></figure>

**Figure 1:** Magnetic field (radial component) predicted at 200-km altitude, from the Equivalent Source Model of [Langlais *et al.*, 2004], superposed onto a shaded relief of Mars. Abbreviations denote Hellas, Syrtis Major, Isidis, Elysium, Apollinaris Patera, Olympus Mons, Tharsis, and Argyre. The three areas correspond to the ones studied by Quesnel, 2006 and Quesnel *et al.*, 2007. See text for details.

## II. A GLOBAL MAP OF THE MARTIAN MAGNETIC FIELD
MGS was in orbit around Mars between Sept. 97 and Nov. 2006. Several orbital phases were originally planed, namely a short AeroBraking (AB) period, followed by a Mapping Orbit (MO) cycle period, MGS being then supposed to serve as a radio relay for other missions. But because a structural problem occurred very early in the mission, the orbit circularization was slowed down and the MO phase was delayed. Due to solar conjunction, a Science Phase Orbit (SPO) was intercalated between what was called AB-1 and AB-2. During this 6-month phase, MGS was on an elliptical orbit (175x17800 km) during which the MAG/ER experiment stayed on. MGS entered its MO cycles in March 1999, and flew during more than 34000 orbits.

MGS measured the three components of the magnetic field, $B_r$ (radial field, positive upward), $B_\theta$ (horizontal North-South, positive South), and $B_\phi$ (horizontal East-West, positive East). The accuracy is of the order of 1 nT for each component [Acuña, 2003], resulting of a combination of sensor uncertainties and spacecraft-generated magnetic fields. Mainly because the altitude is constantly changing along orbit, it is difficult to directly interpret magnetic field measurements in terms of correlations with geologic structures, or in

terms of magnetization contrasts of the lithosphere.

On the Earth, magnetic field measurements made at different altitudes are generally combined into spherical harmonics models [Cain et al., 1989]. However, this method is very sensitive to incomplete data coverage [Schmitz et al., 1989]. In contrast, equivalent source dipole (ESD) [Purucker et al., 1996] approach is less sensitive to geographical data distribution. In addition, it gives insights on what could be the magnetization contrasts in the source region.

Using this method, Langlais et al. [2004] produced a global map of the martian magnetic field at a constant altitude. The model is based on 4840 equivalent sources, located below Mars' surface. Sources are placed onto an equisurface grid, using an icosahedral network (spherical triangles). The mean resolution of the model is relied to the mean space interval between adjacent sources, namely 170 km, representing 2.92 degrees at the equator.

One of the major differences with terrestrial studies is that the magnetization is not induced, i.e. aligned onto the core, main magnetic field. As a consequence, there are no a priori constraints on the magnetization directions; all three components of the magnetization were
thus solved for. The magnetic field computed at a given location is the sum of all dipole contributions. In order to limit the size of the inverse problem, only dipoles within 1500 km of the observation are considered. Measurements of the three MGS phases were used. Local time selection was applied on SPO and MO datasets, in order to avoid the expected larger external magnetic field of the day time. Outliers were removed, by comparing measurements to an earlier model prediction [Purucker et al., 2000]. Measurements were subsequently averaged onto 1°x 1°x10 km, resulting in 75380 bins for the MO dataset and 103996 bins for the joint AB and SPO datasets. Whenever possible, associated variance (with respect to the mean value) was computed and used during the inversion. Several models were computed, using either one or the other data subset, and based on different ESD mean spacing. An iterative scheme was performed, and inversions were stopped when the root mean square (rms) differences between observations and model predictions did not evolve by more than 1% between successive iterations. Final model rms are of the order of 3 nT for each components for the MO dataset, and of the order of 20 nT for the joint AB/SPO dataset.

This discrete model can be used to compute the magnetic field of Mars at a

constant altitude, for altitudes larger than the model resolution. We show on Figure 1 the radial component of the magnetic field, computed at a constant altitude of 200 km. The magnetic field is relatively weak, ± 10 nT, everywhere North of the topographical dichotomy. Two features can nonetheless be identified, close to (185°E, 45°N) and (30°E, 65°N). These two magnetic anomalies reach amplitudes of 60 nT at spacecraft altitude. The case for the southern hemisphere is very different. There are very intense magnetic fields, up to ± 650 nT above Terra Cimmeria and Terra Sirenum. The largest volcanic provinces, Tharsis, Elysium, Olympus, are not associated with magnetic anomalies. The largest impact craters (Hellas, Argyre, Isidis) also lack magnetic signals.

These observations imply that these destructive events took place after the Martian magnetic dynamo had shut down. As a consequence, only some parts of the lithosphere kept its remanent magnetization, acquired while the dynamo was still active. The scenario in which the dynamo started after these events [Schubert et al., 2000] is unlikely, as magnetized terrains are older than the de- or non magnetized ones. The observed magnetic field in the early history of Mars could arise from two kinds of dynamo generation: thermal or compositional
(Stevenson, 2001). An early thermal dynamo, however, is favored by thermal evolution models [e.g. Hauck and Phillips, 2002; Breuer and Spohn, 2003, 2006; Williams and Nimmo, 2004]. Assuming an initially hot core, with a 14% sulphur composition [Breuer and Spohn, 2003, 2006], the growth of the inner core would not have started before 3 billion years (Gyr), which is not compatible with the magnetization of the ancient, > 4 Gyr old units of the southern hemisphere. Furthermore, the thickening of the lithosphere would progressively turn the convecting regime into a stagnant lid regime, thus facilitating the dynamo shutdown. [Breuer and Spohn, 2003]. Any thermal convection in the core does indeed rely on the ability of the mantle to evacuate the core heat to the surface of the planet. However, because of the limited lateral resolution of the MGS-based magnetic models, it is difficult to ascertain the timing of the dynamo and of its shutdown [Langlais, Thébault and

Quesnel, Large impact demagnetization on Mars, in prep.] (denoted as [Langlais et al., in prep.] in the following).

Lower altitude magnetic fields remain unknown. Downward continuations, based on spherical harmonics models [Cain et al., 2003] predict fields as large as 12000 nT. However, this value has to be seen as a minimum one, since models do not take into account the short, very energetic wavelengths of the magnetic field. It is possible that surface magnetic fields of lithospheric origin are of similar amplitude than the Earth' main magnetic field. Such large intensities imply that the magnetic mineral content is higher on Mars, and/or that the ancient, dynamic magnetic field, was more intense. Although non-unique, the magnetization model of [Langlais et al., 2004] predicts magnetization as large as 12 A/m for a 40-km thick lithosphere. In order to improve the resolution of this Martian magnetization global model, more regional studies are needed.

### III. LOCAL STUDIES

Remanent magnetic anomalies trace both the ancient magnetic field and the processes that locally modified the magnetization since the dynamo shut down. However, since it is very difficult to study Mars' surface as a whole, local

studies are easier to performed. This can be done using forward approaches, using simple geometry magnetized bodies (spheres, surface disks, vertical prisms, ...) [Blakely, 1995; Telford et al., 1976]. In this case, magnetization and location
parameters have to be systematically
explored, and the set of parameters
associated with the closest-from-the-
observations predictions are kept.
However, it is not possible to make an
exhaustive exploration of the different
parameters. In addition, some parameters
are too correlated with each other to allow
an unique solution to be found. This
approach is limited, as it implies simple
geometry bodies, and isolated magnetic
anomalies.

### III.1 Semi-inverse approach

The semi-inverse method [Quesnel, 2006]
is a two-step method. It can be applied on
any area. First, a raw, forward approach is
applied, based on homogeneously
magnetized spheres. This first step gives
some hints on the source characteristics
(number of sources, approximated location
and magnetization components). These
characteristics are used is the second step,
as a priori knowledge of the solution in
the inverse problem [Tarantola and
Valette, 1982]. This two-step method was
proven to be efficient on terrestrial studies,
by comparing its results to other
geophysical information on a metamorphic
complex South of the Armorican Massif
[Quesnel, 2006]. On synthetic datasets, the
efficiency of the method was also shown
to be improved when the inversion is
based on multi-altitude datasets, such as

It is also possible to apply an
inverse scheme. Two approaches have
been developed. A local inversion, using
equivalent sources and a priori
information, and a local version of the
equivalent source dipole approach. The
main difference between these methods is
that the second one uses sources that are
located onto an equi-surface array, while
in the first one, both magnetization
parameters and source locations are
investigated.

the one MGS provided [Quesnel *et al.*,
2006].

This method was applied on three
distinct Martian areas [Quesnel, 2006],
representative of the martian magnetic
anomalies. A first area is located within
the heavily crater Terra Sirenum. Magnetic
field is as intense as 1300 nT at 100 km
altitude. While the anomalies appear to be
scattered at low altitude, they coalesce at
high altitude. The semi-inverse method is
particularly adapted in this case. The
second area is located between Hellas (to
the South) and Syrtis Major (to the North).
This area is also cratered, but the magnetic
field is weaker, close to 400 nT at 100-
km altitude. The last area is chosen in the
northern hemisphere, within Vastitas
Boralis, where an anomaly was clearly
identified [Hood and Zacharian, 2001].
This anomaly is only 125 nT large at low altitude. It becomes very weak at high altitude, close to 10 nT.

Many tests were performed, in order to assess the accuracy of the method and its sensitivity with respect to the input parameters. The relationship between magnetic field anomalies and magnetized bodies is non unique. This is the reason why a priori information is required in this approach. A priori data errors, as well as a priori parameter errors, are also used. Large parameter errors will allow a faster convergence through a solution, while a small one will prevent the code to converge to a solution. Resulting models are a posteriori tested versus their associated χ²:

$$ \chi^2 = \sum_{o}^{N_{obs}} \left( \frac{d_m - d_o}{\sqrt{C_{d_o d_o}}} \right)^2 \frac{1}{N_{obs}} $$

where $d_o$ and $d_m$ are the observed and modeled field, $C_{d_o d_o}$ is the covariance associated with the measurements, and $N_{obs}$ is the number of observations, respectively. This value has to be as low as possible. Another criterion that is used is the distribution of the differences between model predictions and observations: it has to be as close to a Gaussian distribution as possible. If necessary, possible outliers are removed, and a priori parameter errors are modified in order to optimize these two quality criteria.

We detail the results associated with the first area [Quesnel et al., 2007]. Data were first selected with respect to their altitude, which leads to two distinct datasets. Source parameters were then investigated using either one or the other or both datasets. Three magnetized bodies were chosen for this area. Models based on a single low- or high-altitude dataset failed to correctly predict high- or low-altitude measurements. The use of both datasets strengthened the results.

Despite the relative proximity of the three magnetic sources, they are associated with different magnetization parameters. They however have a similar and intene magnetic moment, of the order of 2.10¹⁶ A.m². Their depths range between 30 and 60 km. Assuming the magnetic sources to tangent-to-the-surface spheres, their magnetization would be 54 and 32 A/m for the deepest ones, and up to 143 A/m for the shallowest one. Such values are not met on the Earth on such large scales (> 10⁵ km³). The magnetization directions also vary, from a steep (-50°) inclination, to a sub horizontal (~20°) inclination for the last one. The 3-body model predicts more than 57% of the raw signal. Correlation coefficients exceed 0.95 for Br and B₀ (B₀ is more polluted by external fields, and as a consequence, more difficult to model with lithospheric sources). Rms differences are of the order