permet non seulement d'estimer de manière plus précise la géométrie des structures mantelliques, mais aussi les incertitudes sur cette estimation.

Résumé

Cette thèse vise à proposer de nouvelles méthodes permettant de reconstruire la circulation dans le manteau terrestre et l'évolution de la tectonique de surface pour les deux cents derniers millions d'années. Nous utilisons des modèles numériques de convection mantellique dans lesquels la dynamique de surface est comparable à la tectonique terrestre. En combinant ces modèles avec des reconstructions de la tectonique des plaques il est possible d'estimer la structure et l'évolution du champ de température dans le manteau. Jusqu'à présent, l'inclusion des reconstructions de la tectonique des plaques se faisait en imposant des conditions aux limites du modèle (équilibre des forces, vitesses imposées...). Ces techniques, bien que permettant de tester la validité de différents scénarios tectoniques alternatifs, n'autorisent pas de rétroaction dynamique de la convection mantellique sur la tectonique de surface.

Dans ce travail, nous avons développé des techniques d'assimilation de données permettant d'intégrer les reconstructions de la tectonique des plaques dans un modèle numérique tout en laissant se développer de manière auto-cohérente cette rétroaction.

Les techniques développées permettent également de prendre en compte les incertitudes associées aux reconstructions de la tectonique des plaques et de calculer les erreurs sur l'estimation finale de la circulation mantellique.

Dans un premier temps, nous avons développé un filtre de Kalman suboptimal qui permet d'estimer la structure et l'évolution de la circulation mantellique la plus probable à partir d'un modèle numérique de convection et d'une série temporelle d'observations de surface, ainsi que de leurs incertitudes respectives.

Ce filtre a été testé sur des expériences synthétiques. Celles-ci consistent à tenter de retrouver une évolution témoin à partir d'une série temporelle de données issues de cette évolution. Ces expériences ont montré qu'il était possible, en principe, de reconstruire la structure et l'évolution de l'ensemble du manteau à partir d'observations de vitesses et de flux de chaleur à la surface.

Dans un second temps, nous avons développé un filtre de Kalman d'ensemble. Ce filtre

Introduction

The Earth is an active planet, as evidenced by the regular occurrence of earthquakes, volcanic eruptions and the presence of large regions of highly deformed rocks. This activity is the surface expression of a dynamic system comprised of the two outermost layers of the Earth: the crust (of variable thickness between 0 up to 60 km) and the mantle (extending from the base of the crust to 2900 km in depth). To date, the Earth science community uses two different theories to describe the dynamics of the mantlecrust system at a global scale: plate tectonics theory describes the motion at the surface of the Earth while mantle convection theory describes its internal dynamics.

The theory of plate tectonics is a kinematic theory. Its original formulation [START_REF] Mckenzie | The north pacific: an example of tectonics on a sphere[END_REF][START_REF] Morgan | Rises, trenches, great faults, and crustal blocks[END_REF][START_REF] Le Pichon | Sea-floor spreading and continental drift[END_REF] states that the Earth's surface is made of a jigsaw of internally rigid blocks, the plates, each of them moving with respect to the others. The relative motion of these plates is accommodated by deformation along the lines separating plates, called plate boundaries. This theory is successful in predicting the pattern of surface velocities and explaining the repartition of mountain ranges, earthquakes and volcanism (Figure 1 shows for example the fit between seismicity and plate boundaries). It also gives a global framework to integrate various types of data (magnetic, bathymetric, geodetic, paleontologic, geologic...) and provide models of the motions at the surface of the Earth at present (Geologically current plate motion models, e.g. [START_REF] Demets | Geologically current plate motions[END_REF], and in the past (classical plate reconstruction models, e.g. [START_REF] Torsvik | Plate tectonics and net lithosphere rotation over the past 150my[END_REF][START_REF] Seton | Global continental and ocean basin reconstructions since 200ma[END_REF]Müller & Wessel, 2015). FIGURE 1: Synthetic world map taken from Müller & Wessel (2015). This map shows: in white, the plate boundaries from [START_REF] Coffin | Present-day plate boundary digital data compilation[END_REF], in red dots the shallow seismicity (0 -70 km depth) from [START_REF] Engdahl | Global teleseismic earthquake relocation with improved travel times and procedures for depth determination[END_REF], and with areas of lighter illumination the regions of diffuse plate boundaries based on Figure 1 of [START_REF] Gordon | Diffuse oceanic plate boundaries: strain rates, vertically averaged rheology, and comparisons with narrow plate boundaries and stable plate interiors, The history and dynamics of global plate motions[END_REF]. Plate abbreviations: B, Borneo; AN, Antarctica; AR, Arabia; AU, Australia; CA, Caribbean; CAP, Capricorn; CL, Caroline; CO, Cocos; EU, Eurasia; I, Indo-China; IN, India; JF, dynamics of the Earth. Originally, the possibility of mantle convection was studied [START_REF] Holmes | Radioactivity and earth movements[END_REF][START_REF] Pekeris | Thermal convection in the interior of the earth[END_REF][START_REF] Hales | Convection currents in the earth[END_REF]Griggs, 1939b;[START_REF] Runcorn | Convection currents in the earth's mantle[END_REF]) so as to provide a viable driving mechanism for the theory of continental drift of [START_REF] Wegener | The origin of oceans and continents, Trans. JGA Skerl. Originally published as Die Entstehung der Kontinente und Ozeane[END_REF].

The formulation of plate tectonics theory provided a more detailed description of the dynamics of the Earth's surface and it was soon thereafter recognized that the plates correspond to the thermal boundary layer in mantle convection theory [START_REF] Turcotte | Mantle convection and the new global tectonics[END_REF], or in short that plates are mantle convection. Since then, two of the central questions of mantle convection studies have been: how does plate tectonics arise from mantle convection, and what are the dynamical interactions between plates and the rest of the mantle?

The physical description of mantle convection is based on the derivation of partial differential equations governing the flow [START_REF] Turcotte | Mantle convection and the new global tectonics[END_REF]. These equations are highly nonlinear, so the computation of a solution generally requires the development of numerical methods [START_REF] Mckenzie | Convection in the earth's mantle: towards a numerical simulation[END_REF]. The development of models of mantle convection has followed two complementary lines of inquiry to investigate the relationship between surface tectonics and mantle dynamics (see [START_REF] Lowman | Mantle convection models featuring plate tectonic behavior: An overview of methods and progress[END_REF], for a review).

The first consist in using plate reconstruction models, either past or present, as boundary conditions of mantle convection models. In doing so, they provide estimations of the structure and evolution of the flow in the Earth's mantle. We will call these estimations mantle circulation reconstructions. [START_REF] Hager | Kinematic models of large-scale flow in the earth's mantle[END_REF] pioneered this type of studies and showed how the presence of plates at the surface affects the large scale flow of the mantle. Mantle circulation reconstructions have provided valuable insights into the dynamics and evolution of the deep earth mantle structures [START_REF] Bunge | Time scales and heterogeneous structure in geodynamic earth models[END_REF][START_REF] Mcnamara | Thermochemical structures beneath africa and the pacific ocean[END_REF][START_REF] Davies | Reconciling dynamic and seismic models of earth's lower mantle: The dominant role of thermal heterogeneity[END_REF][START_REF] Bower | Lower mantle structure from paleogeographically constrained dynamic earth models[END_REF], the evolution of mantle plumes and their relationship to hotspots [START_REF] Hassan | A rapid burst in hotspot motion through the interaction of tectonics and deep mantle flow[END_REF], changes in the Earth's rotation axis [START_REF] Steinberger | Changes of the earth's rotation axis owing to advection of mantle density heterogeneities[END_REF], sea-level [START_REF] Moucha | Dynamic topography and long-term sea-level variations: There is no such thing as a stable continental platform[END_REF] or dynamic topography [START_REF] Flament | A review of observations and models of dynamic topography[END_REF], for example.

The second line of research consist of implementing the physics necessary to generate plates in convection codes, so that plate tectonics arises naturally from mantle convection models [START_REF] Bercovici | The generation of plate tectonics from mantle dynamics[END_REF]. Tackley (2000a) described this aim as the "Holy Grail" of mantle convection modellers: it would describe physically the dynamics of both the mantle and the surface, and, in doing so, unify the theories of plate tectonics and mantle convection. 16 years later, this aim is still not reached, but the progresses in this direction have led to the emergence of numerical models of mantle convection which display a surface dynamics similar to that of the Earth, at first order [START_REF] Moresi | Mantle convection with a brittle lithosphere: thoughts on the global tectonic styles of the earth and venus[END_REF][START_REF] Moresi | Plate tectonics and convection in the earth's mantle: Toward a numerical simulation[END_REF][START_REF] Richards | Role of a lowviscosity zone in stabilizing plate tectonics: Implications for comparative terrestrial planetology[END_REF][START_REF] Stein | The effect of rheological parameters on plate behaviour in a self-consistent model of mantle convection[END_REF][START_REF] Bercovici | The generation of plate tectonics from mantle convection[END_REF]Tackley, 2000b;[START_REF] Van Heck | Planforms of self-consistently generated plates in 3d spherical geometry[END_REF]. Although the mechanisms of plate generation are still not clearly understood, these models are realistic enough to study the dynamics of seafloor spreading [START_REF] Coltice | Dynamic causes of the relation between area and age of the ocean floor[END_REF][START_REF] Coltice | Convergence of tectonic reconstructions and mantle convection models for significant fluctuations in seafloor spreading[END_REF], the long-term motion of continents [START_REF] Rolf | Statistical cyclicity of the supercontinent cycle[END_REF][START_REF] Yoshida | Preliminary three-dimensional model of mantle convection with deformable, mobile continental lithosphere[END_REF] or the repartition of plate sizes and the dynamics of specific tectonic features [START_REF] Mallard | Subduction controls the distribution and fragmentation of earth's tectonic plates[END_REF].

In this dissertation, we propose to take advantage of the predictive power of such models to develop a new generation of mantle circulation reconstructions, where observations of past and present plate motions are not imposed as boundary conditions, but assimilated in mantle convection models by an inverse method.

The problem of reconstructing mantle circulation is indeed an inverse problem: past and present motions at the surface of the Earth are the consequence of the structure and evolution of the flow in the mantle. Hence, we aim at determining the cause (the structure and evolution of the flow in the mantle) from its consequences (observations of surface motions). From the standpoint of inverse problem theory (see e.g. [START_REF] Wunsch | Discrete inverse and state estimation problems: with geophysical fluid applications[END_REF], for a general introduction to inverse methods), the problem posed by mantle circulation reconstruction corresponds to the general question: how to estimate the state evolution of a system from a time series of data and a dynamical model? Methods solving this type of problem are grouped together under the generic term data assimilation methods.

Data assimilation methods were originally developed in the context of numerical weather prediction [START_REF] Daley | Atmospheric data analysis[END_REF], and have since then found numerous applications, especially in studies of geophysical fluids (i.e. the atmosphere [START_REF] Bengtsson | Dynamic meteorology: data assimilation methods[END_REF][START_REF] Daley | Atmospheric data analysis[END_REF][START_REF] Kalnay | Atmospheric modeling, data assimilation, and predictability[END_REF], the ocean Chassignet & Verron, 2006 and the liquid outer core [START_REF] Fournier | An introduction to data assimilation and predictability in geomagnetism[END_REF]. These areas of research provide extensive resources on methods to solve data assimilation problems. We build on the experience and advances of data assimilation in those fields to design algorithms for the mantle circulation reconstruction problem.

This dissertation is organized as follows:

Chapter 1 presents the mantle circulation direct problem. We derive the partial differential equations governing mantle dynamics and the associated boundary and initial conditions. We also discuss constitutive laws chosen to describe mantle material, with a special focus on the rheology of the mantle, which plays a crucial role in the modelling of mantle convection with plate-like tectonics. Finally, we describe the numerical methods used to solve these equations.

Chapter 2 presents the mantle circulation inverse problem. We first introduce the general formalism of data assimilation and connect the mathematical objects that are defined to their meaning in the context of mantle circulation reconstruction. Then, we detail the data that is available to reconstruct mantle circulation. We finally review the reconstruction strategies that have so far been developed for mantle circulation.

Chapter 3 presents a proof of concept for the application of data assimilation methods to reconstruct mantle circulation from its surface tectonics history. We develop a suboptimal scheme derived from the Kalman filter, where surface velocities and seafloor age maps are used as data to assimilate. We test this algorithm on synthetic experiments using 2-D spherical annulus mantle evolutions. The results obtained show that it is possible, in principle, to reconstruct mantle circulation from the sole (and imperfect) knowledge of the history of surface tectonics.

In Chapter 4, we go one step further and apply the Ensemble Kalman Filter (EnKF) to this problem. The EnKF is a sequential Monte Carlo method particularly adapted to solve high dimensional data assimilation problems with a nonlinear dynamics. We test Introduction the EnKF using synthetic observations, on a 2D-spherical annulus model and compare it with the method developed in Chapter 3. The EnKF performs on average better and is more stable than the former method. Moreover, the efficiency and scalability of the code we developed during this work opens the way to the application of the EnKF to realistic data assimilation with plate reconstructions and a 3D spherical mantle model.

Chapter 5 discusses the results obtained during this thesis and the possible future directions.

Introduction Française

La Terre est une planète active, comme en témoignent l'occurence régulière de tremblements de terre, d'éruptions volcaniques et la présence de grandes régions fortement déformées. Cette activité est l'expression de surface d'un système dynamique constitué des deux couches les plus superficielles de la Terre: la croûte (d'épaisseur variable de 0 à 60 km) et le manteau (s'étendant de la base de la croûte à 2900 km en profondeur).

À ce jour, la communauté des sciences de la Terre utilise deux théories différentes pour décrire la dynamique du système manteau-croûte à une échelle globale: la théorie de la tectonique des plaques décrit les mouvements à la surface de la Terre alors que la théorie de la convection mantellique décrit sa dynamique interne.

La théorie de la tectonique des plaques est une théorie cinématique. Sa formulation originale [START_REF] Mckenzie | The north pacific: an example of tectonics on a sphere[END_REF][START_REF] Morgan | Rises, trenches, great faults, and crustal blocks[END_REF][START_REF] Le Pichon | Sea-floor spreading and continental drift[END_REF]) stipule que la surface de la Terre est faite d'un puzzle de blocs rigides, les plaques, qui se déplacent les unes par rapport aux autres. Le mouvement relatif de ces plaques est accomodé par la déformation le long des lignes séparant les plaques, appelées limites de plaques.

Cette théorie réussit à prédire les vitesses de surface et à expliquer la répartition des chaînes de montagnes, des tremblements de terre et du volcanisme (la figure 2 montre par exemple l'adéquation entre la sismicité et les limites des plaques). Elle fournit également un cadre global pour intégrer différents types de données (magnétiques, bathymétriques, géodésiques, paléontologiques, géologiques ...) et ainsi fournir des modèles des mouvements à la surface de la Terre à l'heure actuelle (modèles de mouvements de plaques géologiquement actuels, par exemple [START_REF] Demets | Geologically current plate motions[END_REF], et dans le passé (modèles classiques de reconstitution de la tectonique des plaques, par exemple [START_REF] Torsvik | Plate tectonics and net lithosphere rotation over the past 150my[END_REF][START_REF] Seton | Global continental and ocean basin reconstructions since 200ma[END_REF]Müller & Wessel, 2015).

FIGURE 2: Carte synthétique du monde prise de Müller & Wessel (2015). Cette carte montre: en blanc, les limites de plaques de Coffin et al. (1997), en points rouges la sismicité peu profonde (0 -70 km de profondeur) de Engdahl et al. (1998), et avec les zones illuminées les régions de limites de plaques diffuses, basées sur la figure 1 de [START_REF] Gordon | Diffuse oceanic plate boundaries: strain rates, vertically averaged rheology, and comparisons with narrow plate boundaries and stable plate interiors, The history and dynamics of global plate motions[END_REF]. Abréviations des plaques: B, Bornéo; AN, Antarctique; AR, Arabie; AU, Australie; CA, Caraïbes; CAP, Capricorne; CL, Caroline; CO, Cocos; UE, Eurasie; I, Indo-Chine; En Inde; JF, Juan de Fuca; NA, Amérique du Nord; NB, Nubia; NC, Chine du Nord; NZ, Nazca; OK, Okhotsk; PA, Pacifique; SA, Amérique du Sud; SC, Mer Scotia; SM, Somalie; Y, Yangtze; T, bassin de Tarim; PH, Philippine.

La théorie de la convection mantellique, au contraire, est une théorie physique. Elle est basée sur l'observation que le manteau, bien que solide, se comporte comme un fluide à des échelles de temps supérieures à quelques milliers d'années [START_REF] Nansen | The earth's crust : its surface-forms, and isostatic adjustment, Avhandlinger utgitt av Det norske Videnskapsakademi i Oslo[END_REF][START_REF] Haskell | The viscosity of the asthenosphere[END_REF]Griggs, 1939a), et utilise les principes et les méthodes de la dynamique des fluides pour étudier le manteau. La théorie de la convection mantellique affirme que le manteau se refroidit et que le mécanisme principal de transfert de chaleur est la convection: le matériel chaud et donc plus léger, à l'intérieur du manteau, s'élève tandis que le matériel froid et plus dense, à la surface, coule, créant ainsi des courants de convection. L'évolution de la théorie de la convection mantellique est étroitement liée à l'exploration de la dynamique de surface de la Terre. A l'origine, on a étudié la possibilité de la convection mantellique [START_REF] Holmes | Radioactivity and earth movements[END_REF][START_REF] Pekeris | Thermal convection in the interior of the earth[END_REF][START_REF] Hales | Convection currents in the earth[END_REF]Griggs, 1939b;[START_REF] Runcorn | Convection currents in the earth's mantle[END_REF] dans le but de fournir un mécanisme d'entraînement viable pour la théorie de la dérive des continents de Wegener (1924). La formulation de la théorie de la tectonique des plaques a fourni une description plus détaillée de la dynamique de la surface de la Terre et il a été bientôt reconnu que les plaques tectoniques correspondent à la couche limite thermique de surface dans la théorie de la convection mantellique [START_REF] Turcotte | Mantle convection and the new global tectonics[END_REF], ou autrement dit, que les plaques font partie de la convection mantellique. Depuis, deux des questions centrales des études de la convection mantellique ont été: Comment la tectonique des plaques découle-t-elle de la convection mantellique et quelles sont les interactions entre les plaques et le reste du manteau?

La description physique de la convection mantellique est basée sur la dérivation d'équations aux dérivées partielles régissant l'écoulement [START_REF] Turcotte | Mantle convection and the new global tectonics[END_REF]. Ces équations sont très non linéaires, donc le calcul d'une solution nécessite généralement le développement de méthodes numériques [START_REF] Mckenzie | Convection in the earth's mantle: towards a numerical simulation[END_REF]. Le développement de modèles de convection mantellique a suivi deux lignes de recherche complémentaires pour étudier la relation entre la tectonique de surface et la dynamique du manteau (voir Lowman, 2011, pour une revue).

La première consiste à utiliser des modèles de reconstitution de plaques, passés ou présents, comme conditions aux limites des modèles de convection mantellique. On obtient ainsi des estimations de la structure et de l'évolution de l'écoulement dans le manteau terrestre. Nous appellerons ces estimations les reconstitutions de la circulation mantellique. [START_REF] Hager | Kinematic models of large-scale flow in the earth's mantle[END_REF] ont été les pionniers de ce type d'études et ont montré que la présence de plaques à la surface affecte l'écoulement à grande Introduction Française échelle du manteau. Les reconstitutions de la circulation mantellique ont permis une meilleure compréhension de la dynamique et l'évolution des structures du manteau en profondeur [START_REF] Bunge | Time scales and heterogeneous structure in geodynamic earth models[END_REF][START_REF] Mcnamara | Thermochemical structures beneath africa and the pacific ocean[END_REF][START_REF] Davies | Reconciling dynamic and seismic models of earth's lower mantle: The dominant role of thermal heterogeneity[END_REF][START_REF] Bower | Lower mantle structure from paleogeographically constrained dynamic earth models[END_REF], de l'évolution des panaches mantelliques et de leur relation avec les points chauds [START_REF] Hassan | A rapid burst in hotspot motion through the interaction of tectonics and deep mantle flow[END_REF], des changements dans l'axe de rotation de la Terre [START_REF] Steinberger | Changes of the earth's rotation axis owing to advection of mantle density heterogeneities[END_REF], dans le niveau marin [START_REF] Moucha | Dynamic topography and long-term sea-level variations: There is no such thing as a stable continental platform[END_REF] ou la topographie dynamique [START_REF] Flament | A review of observations and models of dynamic topography[END_REF], par exemple.

La deuxième ligne de recherche consiste à implémenter la physique nécessaire pour générer des plaques dans les codes de convection, de sorte que la tectonique des plaques émerge naturellement dans les modèles de convection mantellique [START_REF] Bercovici | The generation of plate tectonics from mantle dynamics[END_REF]. Tackley (2000a) décrit ce but comme le "Saint Graal" des modélisateurs de la convection mantellique: on pourrait décrire physiquement la dynamique du manteau et de la surface, et, ce faisant, unifier les théories de la tectonique des plaques et de la convection mantellique. 16 ans plus tard, ce but n'est toujours pas atteint, mais les progrès dans ce sens ont conduit à l'émergence de modèles numériques de convection mantellique qui produisent une dynamique de surface semblable à celle de la Terre, au premier ordre [START_REF] Moresi | Mantle convection with a brittle lithosphere: thoughts on the global tectonic styles of the earth and venus[END_REF][START_REF] Moresi | Plate tectonics and convection in the earth's mantle: Toward a numerical simulation[END_REF][START_REF] Richards | Role of a lowviscosity zone in stabilizing plate tectonics: Implications for comparative terrestrial planetology[END_REF][START_REF] Stein | The effect of rheological parameters on plate behaviour in a self-consistent model of mantle convection[END_REF][START_REF] Bercovici | The generation of plate tectonics from mantle convection[END_REF]Tackley, 2000b;[START_REF] Van Heck | Planforms of self-consistently generated plates in 3d spherical geometry[END_REF]. Bien que les mécanismes de génération des plaques ne soient pas encore clairement compris, ces modèles sont suffisamment réalistes pour étudier la dynamique de formation des fonds océaniques [START_REF] Coltice | Dynamic causes of the relation between area and age of the ocean floor[END_REF][START_REF] Coltice | Convergence of tectonic reconstructions and mantle convection models for significant fluctuations in seafloor spreading[END_REF], le mouvement à long terme des continents [START_REF] Rolf | Statistical cyclicity of the supercontinent cycle[END_REF][START_REF] Yoshida | Preliminary three-dimensional model of mantle convection with deformable, mobile continental lithosphere[END_REF] ou encore la répartition des tailles des plaques et la dynamique de certains contextes tectoniques spécifiques [START_REF] Mallard | Subduction controls the distribution and fragmentation of earth's tectonic plates[END_REF].

Dans cette dissertation, nous proposons de tirer parti du pouvoir prédictif de tels modèles pour développer une nouvelle génération de reconstitutions de la circulation mantellique, dans laquelle les observations des mouvements de plaques passées et présentes ne sont pas imposées comme des conditions aux limites, mais assimilées dans les modèles de convection mantellique par une méthode inverse.

Le problème de la reconstitution de la circulation mantellique est en effet un problème inverse: les mouvements passés et présents à la surface de la Terre sont la conséquence de la structure et de l'évolution de l'écoulement dans le manteau. Par conséquent, nous cherchons à déterminer la cause (la structure et l'évolution de l'écoulement dans le manteau) à partir de ses conséquences (observations des mouvements de surface).

Du point de vue de la théorie des problèmes inverses (voir, par exemple [START_REF] Wunsch | Discrete inverse and state estimation problems: with geophysical fluid applications[END_REF] pour une introduction générale aux méthodes inverses), Le problème posé par la reconstitution de la circulation mantellique correspond à la question générale: Comment estimer l'évolution de l'état d'un système à partir d'une série chronologique de données et d'un modèle dynamique? Les méthodes de résolution de ce type de problème sont regroupées sous le terme générique de méthodes d'assimilation de données.

Les méthodes d'assimilation de données ont été initialement développées dans le contexte de la météorologie [START_REF] Daley | Atmospheric data analysis[END_REF], et ont depuis lors trouvé de nombreuses applications, en particulier dans l'étude des fluides géophysiques (l'atmosphère [START_REF] Bengtsson | Dynamic meteorology: data assimilation methods[END_REF][START_REF] Daley | Atmospheric data analysis[END_REF][START_REF] Kalnay | Atmospheric modeling, data assimilation, and predictability[END_REF], l'océan Chassignet & Verron, 2006 et le noyau externe liquide [START_REF] Fournier | An introduction to data assimilation and predictability in geomagnetism[END_REF]. Ces domaines de recherche fournissent des ressources étendues sur les méthodes permettant de résoudre les problèmes d'assimilation des données. Nous nous appuyons sur l'expérience et les progrès de l'assimilation de données dans ces domaines pour concevoir des algorithmes appliqués au problème de la reconstitution de la circulation mantellique.

Cette thèse est organisée comme suit:

Le chapitre 1 présente le problème direct de la circulation mantellique. Nous dérivons les équations aux dérivées partielles régissant la dynamique du manteau ainsi que les conditions aux limites et les conditions initiales. Nous discutons également les lois constitutives choisies pour décrire les propriétés mécaniques du manteau, en mettant l'accent sur la rhéologie du manteau, qui joue un rôle crucial dans la modélisation de la convection mantellique avec une tectonique de type plaque. Enfin, nous décrivons les méthodes numériques nécessaires pour résoudre ces équations.

Le chapitre 2 présente le problème inverse de la circulation mantellique. Nous introduisons d'abord le formalisme général de l'assimilation des données et relions les objets mathématiques définis à leur signification dans le contexte de la reconstitution de la circulation mantellique. Ensuite, nous détaillons les données disponibles pour reconstruire la circulation mantellique. Nous examinons enfin les stratégies de reconstitution qui ont été développées jusqu'ici pour la circulation mantellique.

Le chapitre 3 présente une preuve du concept pour l'application de méthodes d'assimilation de données pour reconstruire la circulation mantellique à partir de l'histoire de sa tectonique de surface. Nous développons un schéma sous-optimal dérivé du filtre de Kalman, dans lequel les vitesses de surface et les cartes d'âge du fond marin sont utilisées comme données à assimiler. Nous testons cet algorithme sur des expériences synthétiques, en utilisant des évolutions 2D-sphériques annulaires. Les résultats obtenus montrent qu'il est possible, en principe, de reconstituer la circulation mantellique à partir de la connaissance unique (et imparfaite) de l'histoire de la tectonique de surface.

Dans le chapitre 4, nous allons plus loin et appliquons le filtre de Kalman d'Ensemble (EnKF) à ce problème. L'EnKF est une méthode séquentielle de Monte Carlo particulièrement adaptée pour résoudre l'assimilation de données pour des problèmes de grande dimension avec une dynamique non linéaire. Nous testons l'EnKF en utilisant des observations synthétiques, sur un modèle 2D-sphérique annulaire et comparons les résultats avec la méthode développée dans le chapitre 3. L'EnKF fonctionne en moyenne mieux et est plus stable que l'ancienne méthode. De plus, l'efficacité et la scalabilité du code que nous avons développé au cours de ce travail ouvrent la voie à l'application de l'EnKF à une assimilation avec des données provenant des reconstitutions des plaques terrestres et un modèle de manteau sphérique en 3D.

Le chapitre 5 traite des résultats obtenus au cours de cette thèse et des orientations futures possibles.

CHAPTER 1

The Mantle Circulation Forward Problem

Solving the direct problem of mantle circulation consists of specifying how to compute an evolution of mantle circulation from a given initial state. To do so, we use fluid dynamics principles. After justifying the use of fluid dynamics to compute mantle circulation, we derive the partial differential equations governing the flow, then specify the initial and boundary conditions necessary to solve these equations. An additional section is dedicated to the rheology of the mantle, which plays a major role in its dynamics. Finally, we discuss the numerical approximation of such a problem.

Fluid Mechanics for Mantle Convection

The study of mantle convection focuses on the dynamics of the mantle at length scales greater than one kilometer. This enables the modelling of the mantle as a continuous medium, although mantle material is polycrystalline and affected by brittle deformation near the surface.

Moreover, as in many solids, the response of mantle material to stress depends on the frequency content of the forcing. At shorter timescales (typically < 100 years), the mantle behaves like an elastic body, able to propagate seismic shear waves for example. However when submitted to stress over a longer period of time, mantle material deforms permanently. The Maxwell time τ m gives a timescale for this transition from elastic to viscous behavior: τ m = µ 0 /G where µ 0 is the dynamic viscosity and G is the shear modulus (elasticity). For the mantle, the shear modulus is of the order of 10 11 Pa.

If we take a viscosity of 10 26 Pa.s for the lithosphere, and a viscosity of 10 21 Pa.s for the rest of the mantle, we obtain a Maxwell time of 30 Myr and 300 yr respectively.

So, at timescales of 1 Myr or greater, the bulk of the mantle can be modelled as a viscous fluid. In the following, we shall make the same hypothesis for the lithosphere, although elasticity should be taken into account [START_REF] Thielmann | Lithospheric stresses in rayleigh-bénard convection: effects of a free surface and a viscoelastic maxwell rheology[END_REF][START_REF] Moresi | A lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials[END_REF].

Geochemical and seismological evidence suggests that the mantle is chemically heterogeneous at all length scales (see for example [START_REF] Meibom | The statistical upper mantle assemblage[END_REF]; Gurnis (1986a) for geochemistry and [START_REF] Ricard | Is there seismic attenuation in the mantle?[END_REF] for seismology). Although chemical heterogeneities play an important role for mantle convection (see for example [START_REF] Tackley | Mantle geochemical geodynamics[END_REF], for an introduction to mantle geochemical geodynamics), we will consider only thermal heterogeneities in this development.

Conservation Equations

General derivations of equations governing fluid dynamics can be found in classical textbooks of fluid dynamics such as [START_REF] Batchelor | An introduction to fluid dynamics[END_REF]; [START_REF] Chorin | A mathematical introduction to fluid mechanics[END_REF]. Here, we present the series of approximations classically used to model mantle convection, as developed in more details in [START_REF] Turcotte | Mantle convection and the new global tectonics[END_REF], [START_REF] Schubert | Mantle convection in the Earth and planets[END_REF] or [START_REF] Ricard | Physics of mantle convection[END_REF], for example.

General Equations for the Conservation of Mass, Momentum and Energy

The conservation of mass implies that the variation in density ρ of an open volume is compensated by a flow of mass ρu (where u is the velocity) across its boundary. In differential form, it gives

∂ρ ∂t + ∇ • (ρu) = 0, (1.1)
with t the time. The variation of density with time is linked to sound wave propagation and is therefore negligible for the study of mantle convection.

Newton's second law (Newton, 1642(Newton, -1727) ) states that the rate of change of momentum of a particle of fluid is equal to the sum of surface and body forces applied to this particle of fluid. Combined with the equation of conservation of mass and the divergence theorem, this leads to the differential equation

ρ Du Dt = ∇ • σ + f , (1.2)
where σ is the stress tensor and f is the sum of body forces acting on the fluid per unit volume. The only body force acting on mantle material is the gravitational body force, so f = ρg with g the gravitational acceleration (see for example [START_REF] Ricard | Physics of mantle convection[END_REF], section 7.02.2.3.2 for a discussion on the effect of non-inertial forces on the mantle). The gravitational acceleration g depends on the repartition of mass in the mantle according

to g = -∇Ψ and ∇ 2 Ψ = 4πGρ, (1.3)
where Ψ is the gravitational potential and G is the gravitational constant.

The stress tensor is decomposed into an isotropic thermodynamic pressure P and a deviatoric velocity driven stress tensor τ , such that σ = -P I + τ , with I the identity matrix. As discussed previously, we suppose that mantle material behaves viscously. Although mantle anisotropy is established [START_REF] Karato | Seismic anisotropy in the deep mantle, boundary layers and the geometry of mantle convection[END_REF], and viscous anisotropy could affect the dynamics of the flow [START_REF] Pouilloux | Anisotropic rheology of a cubic medium and implications for geological materials[END_REF][START_REF] Mühlhaus | Emergent anisotropy and flow alignment in viscous rock[END_REF], we assume, for simplicity, that the mantle viscosity is isotropic. The viscosity is then described with a bulk viscosity µ v and a shear viscosity µ. This leads to a velocity driven stress tensor τ linked to velocity through

τ = µ v - 2 3 µ ∇ • u I + µ ∇u + [∇u] T . (1.4)
The balance of energy gives

ρ DU Dt = -∇ • q -P ∇ • u + ∇ • (τ • u) + ρH. (1.5)
where U is the internal energy per unit mass, H the rate of energy production (radiogenic heating) per unit mass and q the heat flux.

We apply Fourier's law (Fourier, 1768(Fourier, -1830) ) to link the heat flux q to temperature q = -k∇T with k the thermal conductivity.

By applying the first law of thermodynamics and supposing local thermodynamic equilibrium, we have

DU Dt = C p DT Dt - αT ρ DP Dt + P ∇ • u ρ (1.6)
where α is the thermal expansivity and C p is the heat capacity at constant pressure per unit volume.

We obtain the set of partial differential equations

∇ • (ρu) = 0, (1.7a) ρ Du Dt = -∇P + ∇ • τ + ρg, (1.7b) with g = -∇Ψ and ∇ 2 Ψ = 4πGρ, τ = µ v - 2 3 µ ∇ • u I + µ ∇u + [∇u] T , ρC p DT Dt = ∇ • (k∇T ) + αT DP Dt -∇ • (τ • u) + ρH. (1.7c)
For mantle convection, this system of equation is simplified either by using the anelastic liquid approximation or the Boussinesq approximation (Joseph Valentin Boussinesq, 1842-1929). The anelastic liquid approximation takes into account compressibility and the depth-dependence of density. Derivation of anelastic liquid equation for the mantle can be found in [START_REF] Jarvis | Convection in a compressible fluid with infinite prandtl number[END_REF]; [START_REF] Glatzmaier | Numerical simulations of mantle convection: Time-dependent, three-dimensional, compressible, spherical shell[END_REF]; [START_REF] Bercovici | Three-dimensional convection of an infinite-prandtl-number compressible fluid in a basally heated spherical shell[END_REF]; [START_REF] Alboussière | Reflections on dissipation associated with thermal convection[END_REF]. Here, we will derive the simpler Boussinesq equations.

The Boussinesq Equations

To further simplify the system of equations, all the variables are written as the sum of a reference static solution and a perturbation. In the Boussinesq approximation, the reference state is chosen hydrostatic with a constant temperature T 0 and a constant density ρ 0 . We write

P = P + P ′ , T = T 0 + T ′ and ρ = ρ 0 + ρ ′ .
(1.8)

with P the hydrostatic pressure and P ′ , T ′ and ρ ′ the perturbations for pressure, temperature and density. The gravity can also be decomposed into a one-dimensional (1D)

Chapter 1. The Mantle Circulation Forward Problem profile g and a perturbation g ′ : g = g + g ′ . Following Equation 1.3, we have:

g = -∇Ψ and ∇ 2 Ψ = 4πGρ 0 , (1.9) g ′ = -∇Ψ ′ and ∇ 2 Ψ ′ = 4πGρ ′ .
(1.10)

The hydrostatic pressure P is

∇P = ρ 0 g. (1.11)
We further suppose that the density perturbation depends only on temperature. In the hypothesis that ρ ′ ≪ ρ 0 , we obtain the equation of state

ρ ′ = -ρ 0 αT ′ .
(1.12) with α the constant expansivity.

Compared to the anelastic liquid approximation, the Boussinesq approximation does not take into account compressibility and its main effect on temperature: the presence of a vertical adiabatic temperature gradient. Although the adiabatic temperature difference across the mantle is not negligible, the nonadiabatic temperature difference is high enough to be responsible for most of mantle dynamics. For the energy conservation equation (Eq. 1.7c), it means that the terms αT DP Dt and ∇•(τ •u) are negligible in the Boussinesq approximation [START_REF] Bercovici | Three-dimensional convection of an infinite-prandtl-number compressible fluid in a basally heated spherical shell[END_REF][START_REF] Glatzmaier | Numerical simulations of mantle convection: Time-dependent, three-dimensional, compressible, spherical shell[END_REF][START_REF] Jarvis | Convection in a compressible fluid with infinite prandtl number[END_REF][START_REF] Tackley | Effects of strongly variable viscosity on three-dimensional compressible convection in planetary mantles[END_REF]. This leads to the Boussinesq equations

∇ • u = 0 (1.13a) ρ 0 Du Dt = -∇P ′ + ∇ • µ ∇u + [∇u] T -ρ 0 αT ′ g + g ′ ρ 0 , (1.13b) with g = -∇Ψ and ∇ 2 Ψ = 4πGρ 0 , g ′ = -∇Ψ ′ and ∇ 2 Ψ ′ = 4πGρ ′ , ρ 0 C p DT ′ Dt = ∇ • (k∇T ′ ) + ρ 0 H. (1.13c)
The term g ′ ρ 0 , related to self-gravitation was shown to only affect the computation of gravity, but not the geometry of the flow [START_REF] Zhong | A benchmark study on mantle convection in a 3-d spherical shell using citcoms[END_REF]: we will omit it in the following.

Further Simplifications and Nondimensionalization

We suppose g, C p , k and H constant. The gravity is written g = -ge r with e r the radial unit vector.

We use four characteristic scales for the nondimensionalization:

• the length D = 2900 km, corresponding to mantle thickness,

• the temperature difference ∆T = 1600 K, the temperature drop across the boundary layer (the lithosphere),

• the thermal diffusion time τ D = D 2 ρ 0 Cp k = 265 billion years, i.e. the time scale of diffusion across the whole mantle,

• the reference viscosity µ 0 = 10 23 P a.s -1 , which correspond to the viscosity at the mid-depth in the mantle, and perform the following change of variables in the Boussinesq equations

T ′ = T * ∆T , t = t * τ D , u = u * D/τ D , (1.14) µ = µ * µ 0 , P ′ = P * µ 0 /τ D and ∇ = ∇ * /D. (1.15)
This gives the nondimensional set of equations

∇ * • u * = 0, (1.16a) 1 Pr Du * Dt * = -∇ * P * + ∇ * • µ * ∇ * u * + [∇ * u * ] T + Rae r T * , (1.16b) DT * Dt * = ∇ * 2 T * + R h , (1.16c)
where Ra and Pr are the Rayleigh and Prandtl numbers (John William Strutt, 3rd Baron Rayleigh, 1842-1919;Ludwig Prandtl, 1875-1953), defined as

Ra = ρ 2 0 C p gα∆T D 3 µ 0 k and Pr = µ 0 C p k .
(1.17)

For the mantle, Ra is estimated to be of the order of 10 7 and Pr of the order 10 24 . The very high Prandtl number implies that the left-hand side of equation 1.16b is almost zero: the mantle does not have any inertia.

R h is the nondimensional rate of energy production

R h = ρ 0 HD 2 k∆T . (1.18)
Finally, we obtain the nondimensional equations of conservation of mass, momentum and energy that is used in this thesis to model mantle convection: (1.19b) (1.19c) To solve this set of equation, we need to specify the domain on which they are to be solved as well as the initial and boundary conditions. We also need to specify a constitutive equation for viscosity. These are the topics of the following two sections.

∇ • u = 0, (1.19a) ∇P -∇ • µ ∇u + [∇u] T = RaT e r ,
DT Dt = ∇ 2 T + R h ,

Initial and Boundary Conditions

The volume of the mantle, noted Ω, is at first order a spherical shell of inner radius 3480 km and outer radius 6371 km (if we include the crust). We denote its outer boundary with Γ a and its inner boundary with Γ b . The natural system to locate a point on such a geometry is spherical coordinates. We use the physical convention: the position vector x is defined by its colatitude φ, its longitude θ and its radius r and e φ , e θ , e r are the corresponding unit vectors.

The upper boundary Γ a is in contact with air and water (dynamic viscosity η = 10 -5

and 10 -3 Pa.s, respectively) and the lower boundary Γ b is in contact with the liquid core (dynamic viscosity of the order 10 -3 Pa.s). The drastic drop of viscosity along these boundaries implies free surface conditions, with a deformable surface. The topography of the boundaries are very small compared to the length scale of convection so it is generally neglected in global mantle convection models. Mechanical boundary conditions are then free-slip with rigid, perfectly spherical and impermeable bound-

aries Γ a and Γ b ∀t, ∀x ∈ Γ b , u(x, t) • e r = 0 and σ(x, t) • e r = 0 (1.20a) ∀t, ∀x ∈ Γ a , u(x, t) • e r = 0 and σ(x, t) • e r = 0 (1.20b)
However, regional studies have shown the importance of taking into account the deformation of the surface to obtain a surface dynamics closer to the Earth, especially in the context of subduction [START_REF] Gerya | Future directions in subduction modeling[END_REF]. Recently, a global mantle convection model with a deformable surface has been developed [START_REF] Crameri | A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the 'sticky air'method[END_REF][START_REF] Crameri | Spontaneous development of arcuate single-sided subduction in global 3-d mantle convection models with a free surface[END_REF], 2015) but its high numerical cost limits its application for now.

An alternative mechanical boundary condition would be free slip at Γ b and imposed velocities at Γ a ∀t, ∀x ∈ Γ b , u(x, t) • e r = 0 and σ(x, t)e r = 0, (1.21a)

∀t, ∀x ∈ Γ a , u(x, t) • e r = 0 and u H (x, t) = u a H (x, t), (1.21b)
with u H the horizontal velocity. This boundary condition is used to force the surface dynamics, and we will discuss it in the next chapter on mantle circulation reconstructions.

Both boundaries Γ a and Γ b are supposed isothermal with temperature T a and T b , re-

spectively ∀x ∈ Γ b , T (x) = T b , (1.22a) ∀x ∈ Γ a , T (x) = T a . (1.22b)
The only prognostic equation is the energy conservation Eq.1.19c so the initial condition is entirely defined by an initial temperature field

T i ∀x ∈ Ω, T (x, t = 0) = T i (x).
(1.23)

Rheology and Self-Consistent Plate Generation

Mantle material is polycrystalline. Its response to stress depends on various mechanisms of deformation at the crystal scale [START_REF] Karato | Rheological properties of minerals and rocks[END_REF][START_REF] Kohlstedt | Properties of rocks and minerals-constitutive equations, rheological behavior, and viscosity of rocks[END_REF]. The combination of all these mechanisms results in a complex and nonlinear rheology which is key to understanding mantle dynamics [START_REF] Turcotte | Mantle convection and the new global tectonics[END_REF][START_REF] Davies | Mantle convection[END_REF][START_REF] Bercovici | Mantle dynamics: An introduction and overview[END_REF].

In section 1.1, we supposed a relationship between the deviatoric stress tensor τ and the strain rate tensor ǫ of the form

τ = 2µ ǫ, (1.24)
with µ the viscosity and ǫ = ∇u + [∇u] T /2. This formulation can take into account a wide range of deformation mechanisms by considering a composite rheology, where the effective viscosity µ is

1 µ = 1 µ 1 + 1 µ 2 + ... (1.25)
with µ 1 , µ 2 ... the effective viscosities of each mechanism, potentially varying with temperature, pressure, strain rate, grain size and composition.

In this dissertation, we focus on the relationship between surface tectonics and mantle convection, so we limit the discussion on rheology to the question: what type of viscosity law generates Earth-like tectonics at the surface of convection models? One of the most important characteristics of mantle rheology is the dichotomy between a strong lithosphere and the rest of the mantle, which is softer. This is primarily due to the thermal activation of crystal deformation [START_REF] Gordon | Thermally activated processes in the earth: Creep and seismic attenuation[END_REF][START_REF] Weertman | High temperature creep of rock and mantle viscosity[END_REF], modelled by a viscosity following an Arrhenius law

µ ∝ exp E A + P V A T , (1.26)
with E A an activation energy and V A an activation volume. Both parameters set the amplitude of viscosity variations in the mantle and affect the style of convection. Systematic studies of convection with a temperature dependent rheology show three possible convective regimes, depending on the viscosity contrast between the surface and the bottom of the system [START_REF] Christensen | Convection with pressure-and temperature-dependent nonnewtonian rheology[END_REF][START_REF] Solomatov | Scaling of temperature-and stress-dependent viscosity convection[END_REF][START_REF] Ratcliff | Transitions in thermal convection with strongly variable viscosity[END_REF]. Low viscosity contrasts produce mobile-lid convection in which the flow resembles isoviscous convection. With a high viscosity contrast (greater than 10 4 ), the convection is limited to the lower part of the system, and the surface is frozen: this regime is called stagnant-lid convection. An intermediate regime exists, where the surface is sluggish.

Experimental studies, crystal deformation theory, and geophysical observations indicate an average viscosity profile with a strong lithosphere (µ > 10 25 Pa.s), and a softer underlying mantle with a viscosity ranging from 10 19 -10 21 Pa.s in its upper part to 10 22 -10 23 Pa.s in its lower part [START_REF] King | Reconciling laboratory and observational models of mantle rheology in geodynamic modelling[END_REF]. Consequently, the Earth's mantle convection regime should be stagnant-lid, if its viscosity was only temperature and pressure dependent. The presence of mobile plates at the surface of the Earth, and the high deformation observed at plate boundaries proves otherwise.

It is possible to obtain plate-like dynamics by artificially introducing weak zones in the top boundary layer of a stagnant-lid convective system [START_REF] Christensen | The interaction of a subducting lithospheric slab with a chemical or phase boundary[END_REF][START_REF] Gurnis | Controls of the structure of subducted slabs[END_REF][START_REF] Puster | Mantle convection experiments with evolving plates[END_REF][START_REF] Zhong | Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection[END_REF]. [START_REF] Zhong | Role of faults, nonlinear rheology, and viscosity structure in generating plates from instantaneous mantle flow models[END_REF] identify two diagnostics to evaluate how Earth-like the surface dynamics of a model is.

Plateness (first defined in 2-D by [START_REF] Weinstein | Thermal convection with non-newtonian plates[END_REF] quantifies the localization of deformation at the surface of a dynamic system. The plateness is maximum (generally set to one) when the surface is made of perfectly rigid blocks moving along infinitely small boundaries that concentrate surface deformation. Conversely, the plateness is low when the deformation is distributed uniformly at the surface. The plateness of the Earth is very high: [START_REF] Gordon | Diffuse oceanic plate boundaries: strain rates, vertically averaged rheology, and comparisons with narrow plate boundaries and stable plate interiors, The history and dynamics of global plate motions[END_REF] estimates that the deformation at the surface of the Earth is concentrated in less than 15% of its total area. The second diagnostic is the ratio of surface toroidal to poloidal motion. Any divergence free 3-D velocity field can be decomposed into a toroidal and a poloidal flow. At the surface, the poloidal flow creates divergence or convergence motions whereas the toroidal flow is associated with plate spin and strike-slip motion. On Earth, the ratio of surface toroidal to poloidal flow varies between 0.25 to 0.5 for the last 120 Myr [START_REF] Olson | On the equipartition of kinetic energy in plate tectonics[END_REF][START_REF] Lithgow-Bertelloni | Toroidal-poloidal partitioning of plate motions since 120 ma[END_REF]. Models with imposed weak zones are a valuable tool to study the evolution of preexisting plates and plate boundaries, but they cannot create new plate boundaries.

To generate weak zones at the surface of models self-consistently, it is necessary to take into account a nonlinear relationship between stress and strain rate. However, it is still not clear which relationship best describes lithosphere rheology, so various approaches coexist (see reviews of Tackley, 2000a;[START_REF] Lowman | Mantle convection models featuring plate tectonic behavior: An overview of methods and progress[END_REF][START_REF] Bercovici | The generation of plate tectonics from mantle convection[END_REF][START_REF] Bercovici | The generation of plate tectonics from mantle dynamics[END_REF]. Figure 1.1 summarizes the different relationships tested over the past 40 years.

The first models proposed a power-law relationship between stress and strain rate FIGURE 1.1: Curves of stress versus strain-rate for Newtonian, non-Newtonian power-law rheology (n=3), Bingham plastic, visco-plastic, stick-slip and pseudo-stick-slip rheologies. From [START_REF] Bercovici | The generation of plate tectonics from mantle dynamics[END_REF]. [START_REF] Parmentier | Studies of finite amplitude nonnewtonian thermal convection with application to convection in the earth's mantle[END_REF][START_REF] Christensen | Convection with pressure-and temperature-dependent nonnewtonian rheology[END_REF]. However, they produce very little plateness when the power-law index is 3, which is typically what is measured for mantle rocks [START_REF] Weertman | High temperature creep of rock and mantle viscosity[END_REF]Evans & Kohlstedt, 1995). An increase of the power-law index leads to better localization of the deformation in narrow weak zones [START_REF] Weinstein | Thermal convection with non-newtonian plates[END_REF], but 3-D models yield only moderate toroidal motion [START_REF] Christensen | 3-d convection with variable viscosity[END_REF][START_REF] Bercovici | A simple model of plate generation from mantle flow[END_REF][START_REF] Bercovici | A source-sink model of the generation of plate tectonics[END_REF][START_REF] Čadek | Comparison between newtonian and non-newtonian flow driven by internal loads[END_REF]. To enhance the localization of deformation, [START_REF] Bercovici | A simple model of plate generation from mantle flow[END_REF][START_REF] Bercovici | A source-sink model of the generation of plate tectonics[END_REF] proposed a power-law rheology with index -1 (pseudo-stick-slip in Figure 1.1), which would be the steady-state equivalent of a more complex, time dependent rheology (Tackley, 2000c). This type of rheology produces strain-rate weakening: at low strain-rates, the material behaves like a very viscous fluid but when the strain rate increases above a given threshold, the viscosity decreases so drastically that it causes a decrease in stress. 3-D tests produce encouraging plateness and toroidal to poloidal ratio (Tackley, 1998a). However, such a rheology is difficult to implement in 3-D mantle convection models because of its high nonlinearity (Tackley, 2000c). Viscoplastic rheology (see figure 1.1) is generally preferred. In this formulation, the material behaves like a Newtonian viscous material when the stress is low. When the stress reaches a given threshold, called the yield stress σ Y , the viscosity follows

µ = σ Y 2 ǫII , (1.27)
where ǫII = ǫij ǫij (in Einstein notation) is the second invariant of the strain rate tensor. This type of rheology aims at modelling the brittle and ductile deformation processes that take place in the lithosphere [START_REF] Moresi | Mantle convection with a brittle lithosphere: thoughts on the global tectonic styles of the earth and venus[END_REF][START_REF] Trompert | Mantle convection simulations with rheologies that generate plate-like behaviour[END_REF]. When combined with an asthenosphere of reduced viscosity, they produce plateness and toroidal to poloidal ratios comparable to that of the Earth (Tackley, 2000c;[START_REF] Richards | Role of a lowviscosity zone in stabilizing plate tectonics: Implications for comparative terrestrial planetology[END_REF][START_REF] Stein | The effect of rheological parameters on plate behaviour in a self-consistent model of mantle convection[END_REF]). The progresses of numerical methods (e.g. [START_REF] Yoshida | Application of the yin-yang grid to a thermal convection of a boussinesq fluid with infinite prandtl number in a three-dimensional spherical shell[END_REF][START_REF] Tackley | Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid[END_REF][START_REF] Burstedde | Large-scale adaptive mantle convection simulation[END_REF] and the increase of computational power led to a series of models in 3-D spherical shell geometry selfconsistently generating plate-like tectonics at their surface [START_REF] Richards | Role of a lowviscosity zone in stabilizing plate tectonics: Implications for comparative terrestrial planetology[END_REF][START_REF] Walzer | Mantle convection and evolution with growing continents[END_REF][START_REF] Van Heck | Planforms of self-consistently generated plates in 3d spherical geometry[END_REF][START_REF] Yoshida | Mantle convection with longest-wavelength thermal heterogeneity in a 3-d spherical model: Degree one or two?[END_REF][START_REF] Foley | Generation of plate-like behavior and mantle heterogeneity from a spherical, viscoplastic convection model[END_REF]. Figure 1.2 shows a snapshot of the surface and inner structure of a typical plate-like evolution produced by this type of models. Although these models display Earth-like plateness and toroidal:poloidal ratio, they have several major flaws. Firstly, the range of yield stresses (σ Y ) that produces plate-like tectonics is very limited [START_REF] Foley | Generation of plate-like behavior and mantle heterogeneity from a spherical, viscoplastic convection model[END_REF] and does not correspond to estimated values of the strength of the lithosphere [START_REF] Kohlstedt | Strength of the lithosphere: Constraints imposed by laboratory experiments[END_REF]. However, the introduction of continents [START_REF] Rolf | Focussing of stress by continents in 3d spherical mantle convection with self-consistent plate tectonics[END_REF], melting [START_REF] Lourenço | Melting-induced crustal production helps plate tectonics on earth-like planets[END_REF] and a free surface [START_REF] Crameri | The interaction between subduction-related mantle currents and surface topography[END_REF] tends to increase the yield stress at which plate tectonics appears on models. Secondly, during subduction, both sides dive into the mantle, whereas on Earth, subduction is exclusively singlesided [START_REF] Gerya | Why is terrestrial subduction onesided?[END_REF]. Free surface boundary conditions give encouraging results

for the development of single-sided subduction in such models (Crameri & Tackley, 2014), especially combined with a weak-crust layer [START_REF] Crameri | Parameters controlling dynamically self-consistent plate tectonics and single-sided subduction in global models of mantle convection[END_REF]. Thirdly, the strike-slip boundaries are very diffuse, contrary to the Earth.

Several further developments have been proposed to obtain a surface dynamics closer to the Earth, such as taking into account the viscosity variations with the concentration of fluids in rocks, grain size and the deformation history [START_REF] Bercovici | The generation of plate tectonics from mantle dynamics[END_REF]. However, for the time being, a composite rheology with a Newtonian component of the form of equation 1.26, combined with a plastic behavior (Eq.1.27) remains the best option to study global convection with Earth-like tectonics. Recent applications of such models include comparison of seafloor age distributions and the time scale of spreading fluctuations in models and on Earth [START_REF] Coltice | Dynamic causes of the relation between area and age of the ocean floor[END_REF][START_REF] Coltice | Convergence of tectonic reconstructions and mantle convection models for significant fluctuations in seafloor spreading[END_REF], studies on the supercontinent cycle [START_REF] Rolf | Statistical cyclicity of the supercontinent cycle[END_REF], and on the distribution of plate sizes [START_REF] Mallard | Subduction controls the distribution and fragmentation of earth's tectonic plates[END_REF].

Numerical Approximation

The final system of equations to solve is composed of three conservation equations

∇ • u = 0, (1.28) -∇P + ∇ • µ ∇u + [∇u] T = RaT e r , (1.29) DT Dt = ∇ 2 T + R h , (1.30)
a constitutive law for viscosity, which is a nonlinear function of temperature, pressure and stress

µ = 1 µ n + 1 µ y -1 , (1.31) with µ n ∝ exp E A + P V a T , (1.32)
and

µ y = σ Y 2 ǫII , (1.33)
and the initial and boundary conditions

∀x ∈ Ω, T (x, t = 0) = T i (x), (1.34) ∀x ∈ Γ a ∪ Γ b , u r (x) = 0 and σe r = 0, (1.35) ∀x ∈ Γ a , T (x, t) = T a and ∀x ∈ Γ b T (x, t) = T b . (1.36)
This set of equations (or a variation of it, discussed in the previous sections) is discretized in time using finite differences and in space using finite difference [START_REF] Kageyama | An overset grid in spherical geometry[END_REF], finite volume [START_REF] Harder | A finite-volume solution method for thermal convection and dynamo problems in spherical shells[END_REF][START_REF] Tackley | Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid[END_REF] or finite element [START_REF] Baumgardner | Three-dimensional treatment of convective flow in the earth's mantle[END_REF][START_REF] Zhong | A benchmark study on mantle convection in a 3-d spherical shell using citcoms[END_REF][START_REF] Burstedde | Large-scale adaptive mantle convection simulation[END_REF][START_REF] Kronbichler | High accuracy mantle convection simulation through modern numerical methods[END_REF] methods. Those grid-based methods are preferred to spectral methods, which, although in principle more accurate and more efficient, have so far not been able to treat efficiently large horizontal viscosity variations [START_REF] Balachandar | Viscous dissipation in threedimensional convection with temperature-dependent viscosity[END_REF][START_REF] Zhang | Intense local toroidal motion generated by variable viscosity compressible convection in 3-d spherical-shell[END_REF].

Regardless of the technique, the evolution of a solution is computed by repeating two stages:

1. Flow solution: at a given time step, the velocity and pressure fields are computed by solving Equations 1.28 and 1.29, for a given temperature field.

2. Time stepping: the temperature field at the next time step is computed by solving Equation 1.30, using the former temperature and velocity solution. In this dissertation, we use the finite volume code STAGYY [START_REF] Tackley | Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid[END_REF], which is an extension of STAG3D, continuously developed since 1992 [START_REF] Tackley | Effects of an endothermic phase transition at 670 km depth in a spherical model of convection in the earth's mantle[END_REF][START_REF] Tackley | Effects of strongly variable viscosity on three-dimensional compressible convection in planetary mantles[END_REF]Tackley, , 1998b[START_REF] Tackley | Strong heterogeneity caused by deep mantle layering[END_REF][START_REF] Tackley | Stag3d: A code for modeling thermo-chemical multiphase convection in earth's mantle[END_REF]. Recent developments of this code include the implementation of continents [START_REF] Rolf | Focussing of stress by continents in 3d spherical mantle convection with self-consistent plate tectonics[END_REF], free surface with the sticky air method [START_REF] Crameri | A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the 'sticky air'method[END_REF] and crust production by melting [START_REF] Lourenço | Melting-induced crustal production helps plate tectonics on earth-like planets[END_REF]. This code has been extensively tested with the mantle convection community benchmarks from Blankenbach et al., 1989 (2-D), Busse et al., 1994 (3D Cartesian) , Van Keken et al., 1997 (thermo-chemical convection), Stemmer et al., 2006 (spherical geometry), Tosi et al., 2015 (viscoplasticity) and showed good agreement with those [START_REF] Tackley | Effects of an endothermic phase transition at 670 km depth in a spherical model of convection in the earth's mantle[END_REF][START_REF] Tackley | Effects of strongly variable viscosity on three-dimensional compressible convection in planetary mantles[END_REF][START_REF] Tackley | Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations[END_REF][START_REF] Tackley | Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid[END_REF]. Additionally, Tackley (2000c) performed a convergence test for models with self-consistent generation of plates. STAGYY is parallelized by domain decomposition, using the MPI message-passing library, and shows good scalability [START_REF] Tackley | Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid[END_REF].

The capabilities of this code go far beyond the simple model setup that we use in this dissertation. We chose a relatively simple setup to focus on the development of methods of reconstruction of mantle circulation, but it was important to develop a mantle circulation reconstruction framework on a code which allows more complexities, for future applications.

Conclusion

In this chapter, we presented the solution for the direct problem of mantle circulation:

we showed how we compute an evolution from a mantle convection model featuring plate-like tectonics at its surface and an initial temperature field. However, the initial temperature field of the Earth's mantle is unknown. Knowledge of past and present mantle circulation come from indirect sources of information such as seismic tomography, for the 3D current structure of the mantle and plate tectonics reconstructions, for past mantle circulation. This makes the reconstruction of mantle circulation an inverse problem. The aim of the next chapter is to describe this inverse problem, review methods used to tackle it, and present data assimilation methods, that we apply in chapter 4 and 5.

Résumé du Chapitre

Résoudre le problème direct de la circulation dans le manteau consiste à spécifier comment calculer une évolution de la circulation dans le manteau à partir d'un état initial donné. Pour ce faire, nous utilisons les principes de dynamique des fluides.

Après avoir justifié l'utilisation de la dynamique des fluides pour calculer la circula-

CHAPTER 2 Data and the Mantle Circulation Inverse Problem

The aim of a mantle circulation reconstruction is to estimate the past and present 3-D flow in the mantle, using data on mantle circulation in combination with a dynamical model of mantle convection. Data on mantle circulation is available both at present, with the constraints provided by seismic tomography for example, and in the past, with constraints on the surface tectonics history, for example. Inverse problems dealing with asynchronous data fall in the category of data assimilation problems. This chapter aims at presenting the data assimilation problem for mantle circulation. We start by formulating the mantle circulation problem in the data assimilation framework (first section). Then we detail the data available on the history of mantle circulation (second section). The third section finally presents the applications to the reconstruction of past mantle flow.

A Data Assimilation Framework for Mantle Circulation Problems

Data assimilation methods have originally been developed as a way to estimate the state of the atmosphere at present, so as to obtain an initial state for numerical weather predictions (see for example [START_REF] Ghil | Meteorological data assimilation for oceanographers. part i: Description and theoretical framework[END_REF], for a short account of early data assimilation developments in weather forecasting). The improvement of both numerical weather models and the observational network, along with the increase in computational power, led to the evolution from subjective analysis of sparse observations, to objective analysis of a denser network of observations and finally to the assimilation of data not only spatially, but also through time, i.e. 4-D data assimilation [START_REF] Charney | Use of incomplete historical data to infer the present state of the atmosphere[END_REF][START_REF] Smagorinsky | The relative importance of variables in initial conditions for dynamical weather prediction1[END_REF][START_REF] Rutherford | Data assimilation by statistical interpolation of forecast error fields[END_REF]. Atmospheric data assimilation borrows many concepts from control and estimation theory (see for example [START_REF] Ghil | Applications of estimation theory to numerical weather prediction[END_REF][START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF]. The main difference is in the larger size of the problem and dataset used in atmospheric data assimilation, which pushed the development of specific techniques. The application of data assimilation methods has been progressively opening to various geophysical systems, following the development of observational networks and numerical models in the corresponding disciplines. Current applications include oceanography [START_REF] Chassignet | Ocean weather forecasting: an integrated view of oceanography[END_REF], geomagnetism [START_REF] Fournier | An introduction to data assimilation and predictability in geomagnetism[END_REF], oil and gas reservoir (Evensen, 2009a, Chapter 17), and glaciology [START_REF] Bonan | Data assimilation in glaciology[END_REF], for example. In parallel, the goals of data assimilation methods have diversified: the original aim was to estimate the current state of a system to forecast future evolutions, but other applications have emerged, such as reanalysis of past data (hindcasting), sensitivity analysis of models to different parameters and initial and boundary conditions, assessment of the numerical model, and parameter estimation. Concerning mantle convection, data assimilation techniques provide powerful tools to investigate fundamental questions such as:

• what is the predictive power of state of the art mantle convection models, or how well do we model mantle dynamics?

• what was the circulation of the mantle in the past?

• how sensitive are mantle circulation models to the variation of key parameters, such as rheological ones?

The integration of data in dynamical models is a longstanding problem in global mantle convection studies [START_REF] Hager | Kinematic models of large-scale flow in the earth's mantle[END_REF] and in geodynamics in general (see for example the introductory book by [START_REF] Ismail-Zadeh | Data-driven Numerical Modelling in Geodynamics: Methods and Applications[END_REF]. One of the ways to solve this problem has been to apply data assimilation methods to mantle circulation reconstructions. This approach was pioneered by [START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF] and [START_REF] Ismail-Zadeh | Inverse problem of thermal convection: numerical approach and application to mantle plume restoration[END_REF] who developed variational data assimilation methods for mantle circulation reconstructions. In these articles, they present the data assimilation problem of mantle circulation reconstruction in its continuous form (using a set of partial differential equations similar to the one described in Section 1.2 and associated continuous fields). Here, we choose a different approach and formulate the data assimilation problem directly on the discretized version of mantle flow equations. We use the standard data assimilation formulation advised by [START_REF] Ide | Unified notation for data assimilation: operational, sequential and variational[END_REF] and define the unknown (the state of the mantle), the sources of information (observations, background state and dynamical model) and how each of these components is related to the others.

The Unknown: the Evolution of the True State of the Mantle

The true state gathers all the variables that are necessary and sufficient to compute the 3-D flow of the Earth's mantle at a given time.

In the continuum approximation (see section 1.1), the evolution of mantle flow is described by a set of fields varying continuously in space and time. However, the computation of a dynamical evolution of those fields requires their discretization (see section 1.5). So, in practice, we aim at estimating the evolution of the discretized state of the mantle. In this case, the true state can take the form of a vector containing the values of the discretized fields at a given time step k. This vector is named the true state vector and noted x t k . Then, we write the evolution of the true state as a time series of true state vectors

x t 1 , x t 2 , ..., x t K , (2.1)
with indices {1, ..., K} referring to the numbers of the time steps.

The type and number of discretized fields included in the true state depend on the model chosen for mantle convection. For the model described in Chapter 1, the 3-D flow of the mantle at a given time is entirely defined from the sole knowledge of the temperature field. In this case, the true state vector contains the temperature values of the mantle at points defined on the numerical grid. For a model taking into account the variations of composition in the mantle, we would have to add the discretized compositional fields, for example.

Observed Data on Mantle Circulation

The set of observed data on mantle circulation is organized in the form of a series of observed data vectors

{y o 1 , y o 2 , ..., y o K } (2.2)
where y o k contains the data observed at the time step k ∈ {1, 2, ..., K}.

At each time step k, there is a function that associates a data vector to the set of continuous fields defining the state of the system. This function is approximated by a mapping H k between the discretized state space and the data space. H k is called the observation operator. For any state vector x k , the corresponding data vector y k is given by

y k = H k (x k ) + ǫ h k (2.3)
where ǫ h k represents the errors caused by the hypotheses used to define H k (e.g. discretization and description of the physics of the measurement). Moreover, the observed data are affected by measurement errors, so the relation between the true state and the observed data at time step k is

y o k = H k (x t k ) + ǫ o k + ǫ h k (2.4)
where ǫ o k is the error vector of observed data at time step k.

In mantle circulation problems, the term "observation" is to be understood in a broad sense. It is common to integrate raw observations into a model to obtain global estimates of variables that are directly linked to mantle circulation. We then use these global estimates as "observations" in the data assimilation procedure. We will discuss this point in the following sections.

The Dynamical Model

The dynamical model M predicts the evolution of the state vector from one time step to the next. An example of M for mantle convection has been derived in chapter 1.

As highlighted and discussed throughout chapter 1, the dynamical model M is a numerical approximation of an imperfect physical model describing the evolution of the continuous state of the system through time. From a state vector x k , we can write the state vector at time step k + 1 as

x k+1 = M(x k ) + η k , (2.5)
with η k the modelling errors due to discretization and simplifications in the physical modelling.

During this thesis, we used the numerical code STAGYY (see section 1.5) with a set of options and parameters described in the chapters 3 and 4 as the model M.

The Background State of the Mantle

The background state of the mantle corresponds to the estimation of the state of the mantle before the observed data have been integrated. This a priori state comes from our knowledge of the physics of mantle convection, or the integration of another set of observed data, for example. The background state is linked to the first true state vector through

x b = x t 1 + ǫ b , (2.6)
where ǫ b is the error vector of the background state vector.

In this work, we used the results of a very long evolution computed with the dynamical model M to infer the average temperature field in the mantle and evaluate the amplitude of the background error vector (see section 3.4.1 for more details).

Summary

We have described the components of the general problem of mantle circulation reconstructions:

• the unknown: a time series of true state vectors {x t 1 , x t 2 , ..., x t K }

• three sources of information, all affected by errors:

the background state x b the time series of observed data {y o 1 , y o 2 , ..., y o K } and their associated observation operators {H 1 , H 2 , ..., H K } -the dynamical model M.

In the next section, we focus on the description of the observed data currently used to reconstruct mantle circulation. Then, in Section 2.3 we describe the algorithms used to solve the mantle circulation inverse problem.

Data on Mantle Circulation

Currently, two types of "observed data" are used in dynamical mantle circulation reconstructions:

• a temperature field of the mantle at present, which is derived from seismic tomography,

• a time series of surface velocities, which is the synthesis of a wide range of observations on the history of the surface of the Earth.

In the following, we detail the procedure used to obtain these data, and discuss the possible errors that affect them.

Mantle Temperature Field at Present

If we consider purely thermal convection, the density at a given depth is solely dependent on temperature. For example, in the model developed in chapter 1, the density anomaly ρ ′ is linked to the temperature anomaly T ′ through ρ ′ = -ρ 0 αT ′ .

(2.7) with α the expansivity and ρ 0 the reference density. Hence, to obtain an estimate of the temperature field in the mantle at present, we need:

• a reference temperature profile (In the Boussinesq approximation, it is simply a constant temperature, see Section 1.2)

• a 3-D model of density anomaly This model of density anomaly is derived from wave speed anomaly models inferred from seismic tomography.

Seismic Tomography

Dziewonski & [START_REF] Dziewonski | Preliminary reference earth model[END_REF] computed the first widely accepted reference profile of the elastic properties of the Earth: the Preliminary Reference Earth Model (PREM).

In particular, the PREM provides estimates of the velocity of compressional and shear waves (P-and S-waves), and density in the Earth as a function of depth. With these estimates it is possible to predict the travel time of any seismic wave for a given earthquake and a given seismic station. In practice, the predicted travel times do not fit exactly the observed travel times. Such discrepancies are due to the lateral heterogeneity of Earth's material elastic properties, interpreted as variations of either temperature or composition: the S-wave velocity v s and the P-wave velocity v p are

v s = G ρ and v p = K + 4 3 G ρ (2.8)
where K , G and ρ are the bulk modulus, the shear modulus and the density. At a given depth, each of them varies with temperature and composition.

Seismic tomography refers to inverse methods that aim at mapping in 2-D or 3-D the lateral heterogeneity of the elastic properties of the Earth, using seismograms as data, and a forward model specifying how seismic waves propagate in the Earth. The three components of the technique, i.e. the inverse method, the forward model, and the type of data extracted from the seismograms vary depending on the target region of mapping, the type of physical properties that are to be determined, the restrictions on computational power, and the availability of data. A full review of seismic tomography techniques is beyond the scope of this dissertation, and we refer the reader to [START_REF] Thurber | 1.10 -theory and observations -seismic tomography and inverse methods[END_REF] for an introduction to the field and to [START_REF] Rawlinson | Seismic tomography: a window into deep earth[END_REF]; [START_REF] Trampert | Global imaging of the earth's deep interior: Seismic constraints on (an)isotropy, density and attenuation[END_REF]; [START_REF] Liu | Seismic imaging: from classical to adjoint tomography[END_REF] for recent reviews.

Seismic tomography studies of the whole mantle have focused on estimating 3-D models of S-wave speed anomalies (e.g. [START_REF] Gu | Models of the mantle shear velocity and discontinuities in the pattern of lateral heterogeneities[END_REF][START_REF] Panning | A three-dimensional radially anisotropic model of shear velocity in the whole mantle[END_REF][START_REF] Kustowski | Anisotropic shear-wave velocity structure of the earth's mantle: A global model[END_REF][START_REF] Ritsema | S40rts: a degree-40 shear-velocity model for the mantle from new rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements[END_REF][START_REF] Auer | Savani: A variable resolution whole-mantle model of anisotropic shear velocity variations based on multiple data sets[END_REF][START_REF] Chang | Joint inversion for global isotropic and radially anisotropic mantle structure including crustal thickness perturbations[END_REF], P-wave speed anomalies (e.g. [START_REF] Boschi | Whole earth tomography from delay times of p, pcp, and pkp phases: Lateral heterogeneities in the outer core or radial anisotropy in the mantle?[END_REF][START_REF] Li | A new global model for p wave speed variations in earth's mantle[END_REF], or both (e.g. Su &

The differences in tomographic models are attributed to the variations in datasets, but also to the non-uniqueness of the solution to the inverse problem, which is dealt with differently depending on studies, using a wide range of inverse methods and parametrization of the forward model. The seismological community has developed several techniques to assess the reliability of tomographic models (see [START_REF] Rawlinson | Chapter one-seismic tomography and the assessment of uncertainty[END_REF], for a review). A direct computation of the a posteriori error covariance matrix associated to the estimated S-and/or P-wave anomaly models would provide a much needed way to evaluate the information that could be extracted from those, and avoid geological and geodynamical misinterpretations (see [START_REF] Foulger | Caveats on tomographic images[END_REF], for a review on the pitfalls of seismic tomography interpretation for non-specialists). However, the size of the whole mantle seismic tomography problem makes the computation of this matrix very costly (and even intractable for some models). In practice, partial resolution tests are preferred (see examples in the previously cited works). Such tests do not quantify directly the absolute amplitude of errors on tomographic models. A project of 3-D Reference Earth Model for S-and P-wave speed, density, and anisotropy with the estimation of uncertainties is under construction [START_REF] Moulik | 3-dimensional Reference Earth Model Webpage[END_REF]. Such model would provide a sound basis for geodynamical interpretation.

Meanwhile, several studies have adopted a more pragmatic approach to assess the robustness of the structures illuminated by seismic tomography. The idea is to compare various tomographic models in order to extract their common features, supposing that if tomographic models developed using different parametrizations, observations and inverse methods have common structures, then it is likely that these structures are not artefacts but reflect true heterogeneities of Earth's material properties. For example, [START_REF] Lekic | Cluster analysis of global lower mantle tomography: A new class of structure and implications for chemical heterogeneity[END_REF] followed this approach and performed a cluster analysis of tomographic models of the lower mantle to infer its large scale structure. For the whole mantle, [START_REF] Becker | A comparison of tomographic and geodynamic mantle models[END_REF] propose tomography-based benchmarks for geodynamics by computing a weighted average of available P-wave and S-wave models (called Pmean and Smean). They regularly update these models taking into account new tomographic studies (Becker & Boschi, 2016).

Conversion of Wave Speed Anomaly Models to Temperature Fields

Equation 2.8 shows that S-and P-wave speeds depend on three material properties: the bulk modulus K, the shear modulus G and the density ρ. These material properties vary with temperature, pressure and composition. The first order relationship between a lateral relative density anomaly δρ/ρ = δ ln ρ and a lateral relative S-or P-wave speed anomaly δv s /v s = δ ln v s or δv p /v p = δ ln v p can be written as

δ ln ρ = R ρ/s δ ln v s and δ ln ρ = R ρ/p δ ln v p (2.9)
with R ρ/s and R ρ/p the partial derivatives ∂ρ/∂v s and ∂ρ/∂v p . At constant depth, δ ln ρ, δ ln v s , δ ln v p , R ρ/s and R ρ/p depend on temperature and composition, so the interpretation of wave speed anomaly models in terms of density variations of purely thermal origin is not rigorously exact [START_REF] Mattern | Lower mantle composition and temperature from mineral physics and thermodynamic modelling[END_REF][START_REF] Stixrude | Thermodynamics of mantle minerals-i. physical properties[END_REF], 2011). However, it is generally the assumption taken in mantle circulation reconstructions (e.g. [START_REF] Conrad | Seismic tomography, surface uplift, and the breakup of gondwanaland: Integrating mantle convection backwards in time[END_REF][START_REF] Steinberger | Changes of the earth's rotation axis owing to advection of mantle density heterogeneities[END_REF][START_REF] Spasojevic | Adjoint models of mantle convection with seismic, plate motion, and stratigraphic constraints: North america since the late cretaceous[END_REF], where the density anomaly of the mantle below the lithosphere is computed with either constant or depth-dependent R ρ/s or R ρ/p , estimated using mineral physics constraints (e.g. [START_REF] Karato | Origin of lateral variation of seismic wave velocities and density in the deep mantle[END_REF][START_REF] Cammarano | Inferring upper-mantle temperatures from seismic velocities[END_REF]. The S-wave velocity anomalies are generally preferred in these computations due to their better resolution of upper mantle heterogeneities.

Instantaneous models of mantle dynamics have also played a role in justifying such assumption and providing a scaling between wave speed anomalies and density anomalies (see [START_REF] Forte | Constraints on seismic models from other disciplines-constraints on 3-d seismic models from global geodynamic observables: Implications for the global mantle convective flow[END_REF], for a review). Instantaneous models of mantle circulation estimate the flow of the mantle at present from:

• data such as seismic tomography models, the surface kinematics, the geometry of the geoid, the value of gravity anomalies, topography...

• a physical model of instantaneous mantle flow [START_REF] Richards | Geoid anomalies in a dynamic earth[END_REF][START_REF] Ricard | Geoid heights and lithospheric stresses for a dynamic earth[END_REF][START_REF] Ricard | A geodynamic model of mantle density heterogeneity[END_REF], combined with models of the material properties of the mantle. [START_REF] Hager | Lower mantle heterogeneity, dynamic topography and the geoid[END_REF] used the first models of seismic tomography to infer a density field and compute the geoid from it. In doing so, they were able to relate to first order seismic and geodynamical observables. By complementing the tomography-derived density model with a model of subducted slabs, [START_REF] Hager | Long-wavelength variations in earth's geoid: physical models and dynamical implications[END_REF] were able to explain 90% of the current geoid. The improvement of both observations and models has produced a variety of instantaneous mantle circulation models which aim at reconciling seismic and geodynamical observations. For example, [START_REF] Simmons | Constraining mantle flow with seismic and geodynamic data: a joint approach[END_REF][START_REF] Simmons | Joint seismic, geodynamic and mineral physical constraints on three-dimensional mantle heterogeneity: Implications for the relative importance of thermal versus compositional heterogeneity[END_REF][START_REF] Simmons | Gypsum: A joint tomographic model of mantle density and seismic wave speeds[END_REF] perform a joint inversion of geodynamical and seismic observable to determine the optimum R ρ/s profile and infer the temperature and compositional contributions to the variations of densities. [START_REF] Stadler | The dynamics of plate tectonics and mantle flow: From local to global scales[END_REF], [START_REF] Ghosh | Effects of lateral viscosity variations on the geoid[END_REF] and [START_REF] Ghosh | Plate motions and stresses from global dynamic models[END_REF] combine tomography-derived density anomalies and models of slabs to determine the state of stress in the lithosphere and [START_REF] Wang | Joint modeling of lithosphere and mantle dynamics: Evaluation of constraints from global tomography models[END_REF] build on this work to evaluate the capacity of various tomography models to produce accurate geodynamic observables.

So far, dynamical mantle circulation models have only used estimated temperature fields from tomography, or from an integration of other variables at present from joint geodynamical and seismic inversions (e.g. [START_REF] Glišović | A new back-and-forth iterative method for time-reversed convection modeling: Implications for the cenozoic evolution of 3-d structure and dynamics of the mantle[END_REF]. However, the continuous developments in both instantaneous and dynamical mantle circulation models bring these areas of research ever closer and opens the way to a direct integration of present-day raw observations in a dynamical model. [START_REF] Mckenzie | The north pacific: an example of tectonics on a sphere[END_REF], [START_REF] Morgan | Rises, trenches, great faults, and crustal blocks[END_REF] and Le Pichon (1968). We then describe the current models for past surface motions, which are to a large extent still based on these principles and observations.

Surface Kinematics History of the Mantle

In the beginning of the 20th century, Alfred Wegener (1880Wegener ( -1930) ) proposed the continental drift hypothesis: the continents once formed a single landmass (called Pangaea)

and have been drifting apart since then. Throughout his career, he collected observations supporting his theory and documented them in the successive editions of his book, "Die Entstehung der Kontinente und Ozeane" (see for example [START_REF] Wegener | The origin of oceans and continents, Trans. JGA Skerl. Originally published as Die Entstehung der Kontinente und Ozeane[END_REF] the English translation of the third edition). Wegener bases his theory on the complementary shape of continents and the continuity of geological formations, paleoclimates and fossil repartition once continents are reunited.

Although the continental drift theory was rejected by a large part of the community [START_REF] Oreskes | The rejection of continental drift: Theory and method in American earth science[END_REF], it gained a few early supporters. For example, Du Toit (1937) collected further geological observations supporting the theory of continental drift. [START_REF] Holmes | Radioactivity and earth movements[END_REF][START_REF] Holmes | Principles of physical geology[END_REF] proposed mantle convection as a viable mechanism for continental drift.

Geodesist Felix Andries Vening-Meinesz pioneered the oceanic geodetic studies and interpreted the negative gravity anomalies at oceanic trenches (the deepest regions of the ocean floor) as the locations of downwelling convective currents (see [START_REF] Vening-Meinesz | Gravity expeditions at sea 1923-1938[END_REF], for a synthesis of his gravity expeditions between 1923 and 1938). Griggs (1939b) designed experimental models of continental drift driven by convection currents.

The 1960 -1970 decade saw the progressive formulation and acceptance of the plate tectonics theory, starting with the sea-floor spreading theory [START_REF] Dietz | Continent and ocean basin evolution by spreading of the sea floor[END_REF][START_REF] Hess | History of ocean basins[END_REF] and new evidence of large scale lateral motion of continents [START_REF] Runcorn | Climatic change through geological time in the light of the palaeomagnetic evidence for polar wandering and continental drift[END_REF] and finishing with the "new global tectonics" (i.e. plate tectonics) [START_REF] Morgan | Rises, trenches, great faults, and crustal blocks[END_REF][START_REF] Le Pichon | Sea-floor spreading and continental drift[END_REF]).

The intense marine geophysical exploration of the previous decades (started for military purposes during WWII and continued from then on) had produced an extensive dataset of the bathymetry, heat flow, seismic properties, gravity and magnetic anomalies of the seafloor, giving a global view of the oceanic domain and requiring new interpretations. Parallelly, the progresses in seismology made possible the detection and precise localization of ever smaller and distant earthquakes, as well as a better determination of earthquake properties such as the focal mechanism, and the associated slip (see e.g. [START_REF] Isacks | Seismology and the new global tectonics[END_REF], for a synthesis of data interpreted in the light of plate tectonics theory). These studies led to the characterization of 5 main features of the ocean floor:

• abyssal plains, which represent most of the ocean floor, show very localized and limited deformation, and display "zebra patterns" of magnetic anomalies, i.e. elongated stripes of alternating positive and negative magnetic anomalies [START_REF] Mason | Magnetic survey off the west coast of north america, 32 n. latitude to 42 n. latitude[END_REF]. These zebra pattern are oriented parallely to the mid-ocean ridges.

• mid-ocean ridges, which form the largest mountain ranges on Earth and show remarkable continuity [START_REF] Ewing | Some problems of antarctic submarine geology[END_REF][START_REF] Heezen | The floors of the oceans i. the north atlantic[END_REF], see also figure 2.2). They are associated with a high heat flow, volcanic activity, shallow earthquakes and a fracture zone resembling a rift at their center. Mid-oceanic ridges have been identified as spreading centers in Hess's and Dietz's theory of seafloor spreading.

• deep trenches, which are associated with mountain ranges on adjacent continents or island arcs, show intense seismic activity with deep earthquakes organized on a plane dipping towards the continents or the island arc [START_REF] Benioff | Seismic evidence for the fault origin of oceanic deeps[END_REF].

Deep trenches are identified as zones of convergence where the oceanic floor dive into the mantle in Hess's and Dietz's theory of seafloor spreading, later labeled as subduction zones.

• fracture zones, linear features of several hundred to several thousand kilometers, perpendicular to mid-ocean ridges. These fracture zones are associated with an offset of the linear features they traverse (zebra pattern of magnetic anomalies and mid-ocean ridges) and hence appear as zones of shearing. Wilson (1965b) and [START_REF] Sykes | Mechanism of earthquakes and nature of faulting on the mid-oceanic ridges[END_REF] showed these fractures correspond to transform faults which link divergent or convergent boundaries to form a consistent network of mobile belts dividing the surface of the Earth into rigid plates. Moreover, it was recognized that fracture zones are seismically active only on the sections that connect two other active boundaries, so that the rest of them are interpreted as fossil transform faults.

• volcanic islands and seamounts generally organized in linear chains. The radioactive dating of the volcanic rocks of pacific chains showed a regular age progression, with the younger islands being closer to the east pacific rise (i.e. the Pacific ocean ridge). [START_REF] Wilson | A possible origin of the hawaiian islands[END_REF]Wilson ( , 1965a) interpreted these chains as the result of successive volcanic eruptions due to localized and fixed rising mantle currents below the moving plates. These hotspots were later connected to deep mantle plumes by [START_REF] Morgan | Convection plumes in the lower mantle[END_REF].

Simultaneously, the study of continental rocks brought key elements that led to the acceptance of plate tectonics theory. These continental arguments are the result of the interaction of studies on the remanent magnetism of rocks, radioactive dating and progresses in the understanding of the origin of the Earth's magnetic field, with the emergence of the geodynamo theory (Elsasser, 1946a,b).

First, a collection of paleomagnetic data of continental rocks indicated an apparent variation of the location of the magnetic pole through geological times. This apparent polar wander appeared to be different depending on continents. [START_REF] Runcorn | Climatic change through geological time in the light of the palaeomagnetic evidence for polar wandering and continental drift[END_REF] showed that lateral motions of continents could reconcile the apparent polar wandering paths. Supposing that the magnetic pole has not moved over geological times, it is then possible to track the paleolatitude of continents with paleomagnetic data (the paleolongitude cannot be estimated because the Earth's magnetic field is radially symmetrical).

The second observation was the frequent reversed polarity measured in the remanent magnetism of rocks of various ages. [START_REF] Cox | Geomagnetic polarity epochs and pleistocene geochronometry[END_REF] used the newly developed K-Ar dating method on samples coming from 6 volcanic eruptions and were able to show the consistency of the polarity among rocks of the same age, which suggested the existence of global polarity reversals of the magnetic field. Additionally, they could date the reversal events, and obtained in this way the first geomagnetic polarity time scale, constantly updated since then [START_REF] Ogg | Chapter 5 -geomagnetic polarity time scale[END_REF]. [START_REF] Morley | An explanation of magnetic banding in ocean basins, The Oceanic Lithosphere[END_REF] and [START_REF] Vine | Magnetic anomalies over oceanic ridges[END_REF] connected these findings with the seafloor spreading theory, stating that if the ocean floor is produced at mid-ocean ridges and moving away from them with time, then the zebra patterns of magnetic anomalies on the seafloor correspond to the fossilization of the magnetic polarity through time and should be organized symmetrically around the mid-ocean ridges. This prediction was then confirmed by [START_REF] Heirtzler | Crustal structure of the mid-ocean ridges: 3. magnetic anomalies over the mid-atlantic ridge[END_REF] and [START_REF] Heirtzler | Magnetic anomalies over the reykjanes ridge[END_REF][START_REF] Heirtzler | Marine magnetic anomalies, geomagnetic field reversals, and motions of the ocean floor and continents[END_REF]. Moreover, the identification of specific anomalies from the center of the ridge outward, along with a magnetic polarity timescale provided a way to estimate the age of the ocean floor depending on the distance of the ridge, and thus to deduce the rate of seafloor spreading.

The final concept to formulate plate tectonics theory was found independently by [START_REF] Mckenzie | The north pacific: an example of tectonics on a sphere[END_REF] and [START_REF] Morgan | Rises, trenches, great faults, and crustal blocks[END_REF]: since the Earth's surface is made of a jigsaw of internally rigid plates in constant motion, then the plate motions can be described with principles of solid body rotation on a sphere. More precisely, we can use Euler's rotation theorem (Leonhard Euler 1707-1783): a solid body motion on a sphere can always be described as a rotation around an axis passing through the center of the sphere, and the composition of two rotations is also a rotation. Consequently, absolute plate motions and relative motions of plate pairs are rotations, and are uniquely defined by their respective angular velocities and rotation poles (called the Euler poles and corresponding to one of the two points defined by the intersection of the rotation axis with the surface of the Earth). [START_REF] Morgan | Rises, trenches, great faults, and crustal blocks[END_REF] then deduced that transform faults could only develop along "small circles", i.e. circles at the surface of the Earth of constant latitude with respect to the Euler pole (see figure 2.3). By the end of 1968, Earth scientists had built a theory which at first order explains the surface dynamics of the Earth, and allows the integration of local observations of different nature, coming from almost all the disciplines of Earth sciences, into a consistent framework. The following decades saw the emergence of space geodesy (see [START_REF] Blewitt | [END_REF] for a review), which made possible:

• the direct measurements of surface motions with techniques such as Very Long Baseline Interferometry, Satellite Laser Ranging, and the Global Positioning System. We can now measure surface motions and infer plate motions which could not be deduced with classical plate tectonics principles (e.g. [START_REF] Demets | Geologically current plate motions[END_REF].

Those techniques also led to the production of instantaneous maps of plate motions (e.g. [START_REF] Argus | The angular velocities of the plates and the velocity of earth's centre from space geodesy[END_REF][START_REF] Kogan | Current global plate kinematics from gps (1995-2007) with the plate-consistent reference frame[END_REF], and the measurement of plate deformation (e.g. [START_REF] Kreemer | A geodetic plate motion and global strain rate model[END_REF].

• a better localization of ships during measurement campaigns at sea. This resulted in an improved mapping of magnetic anomalies identifications (e.g. [START_REF] Seton | Community infrastructure and repository for marine magnetic identifications[END_REF], see also Figure 2.5).

• the production of high resolution maps of gravity anomalies (e.g. Sandwell &

Smith, 2009, see also Figure 2.6), which considerably increased the database of fracture zones, seamounts and other features of the ocean floor (e.g. [START_REF] Wessel | Global distribution of seamounts inferred from gridded geosat/ers-1 altimetry[END_REF][START_REF] Matthews | The tectonic fabric of the ocean basins[END_REF]Wessel et al., 2015). symmetrical spreading rates, it is possible to estimate earlier finite rotations between these two plates by studying the sole pacific plate data (e.g. [START_REF] Seton | Global continental and ocean basin reconstructions since 200ma[END_REF], but this method introduces further uncertainties in the reconstruction. When both sides of the ridge have been destroyed, the reconstructions rely on the interpretations of onshore geology and extrapolations of the direction and the rate of spreading at more recent times (e.g. [START_REF] Rowley | Extrapolating oceanic age distributions: lessons from the pacific region[END_REF]. In this case, the uncertainties on the reconstructions are high, but difficult to quantify. [START_REF] Torsvik | Plate tectonics and net lithosphere rotation over the past 150my[END_REF] defines the "world uncertainty" to quantify the increase of the global uncertainty of plate reconstructions through time, due to the absence of older ocean floor: the world uncertainty corresponds to the fraction of the Earth's lithosphere which has been subducted since a given time. For example, the oceanic floor older than 140 Myr has almost entirely been subducted, which brings the world uncertainty at 60% at that time. Besides, the time resolution of the finite rotation series is directly dependent on the frequency of polarity reversal of the magnetic field. One of the most important lack of data is the Cretaceous superchron [START_REF] Granot | Geomagnetic field variability during the cretaceous normal superchron[END_REF]: between 121 and 83 million years ago, the polarity of the Earth's magnetic field appears to have been stable, so a large portion of the ocean seafloor cannot be dated with precision using magnetic anomaly identifications, as is shown in Once the motion of divergent plate boundaries is determined, it is possible to infer the relative motions of other boundaries using geometrical constraints. To determine the relative motions of all plate pairs, we need to construct a plate circuit such that each plate is connected to the plate circuit by at least one divergent boundary where the finite rotations are computable, as Le Pichon (1968) had constructed at present. However, it is not possible to build such a circuit for the last 200 Myr. Some plates, such as the Philipine plate at present [START_REF] Bird | An updated digital model of plate boundaries[END_REF], have very limited portions of divergent boundaries and are thus almost isolated from the rest of the global circuit. Before 83 Myr, the Pacific plate cannot be linked to the global circuit either [START_REF] Seton | Global continental and ocean basin reconstructions since 200ma[END_REF].

Moreover, the uncertainties on finite rotations propagate through the circuit and two alternative circuits can lead to different results (e.g. [START_REF] Steinberger | Prediction of emperor-hawaii seamount locations from a revised model of global plate motion and mantle flow[END_REF]. Finally, some parts of the surface of the Earth have undergone deformation, for example South East Asian plates. In these places, plate tectonics geometrical rules do not apply (e.g. [START_REF] Holt | On the determination of selfconsistent strain rate fields within zones of distributed continental deformation[END_REF][START_REF] Gordon | Diffuse oceanic plate boundaries: strain rates, vertically averaged rheology, and comparisons with narrow plate boundaries and stable plate interiors, The history and dynamics of global plate motions[END_REF]. To produce a global model of plate reconstruction it is thus necessary to take into account other sources of information. These complementary sources of informations are of two types: additional observations on the relative motion of plates, and observations on the absolute motion of plates.

Observations on the relative motions of plates include for example paleontological, structural, stratigraphic and regional geology constraints (e.g. [START_REF] Hosseinpour | Tectonic evolution of western tethys from jurassic to present day: coupling geological and geophysical data with seismic tomography models[END_REF][START_REF] Barnett-Moore | Assessing discrepancies between previous plate kinematic models of mesozoic iberia and their constraints[END_REF]. The integration of such constraints into global reconstructions is not done in a systematic way. On the contrary, it generally requires choices and interpretations from the geoscientist who builds the reconstruction. The resulting uncertainties are then difficult to quantify.

Observations on the absolute plate motions primarily include data from island and seamount chains, and paleomagnetism. Figure 2.8 represents the linkage between the relative motion of a plate with its neighboring plates, and its absolute motion. All these motions are represented for a plate which is stationary for the last 80 Myr, for illustrative purposes. The figure also depicts the predicted apparent polar wander path (i.e. the time evolution of the direction of the paleopole as recorded by the remanent magnetism of the rocks belonging to this plate) and the geometry and age progression of an island and seamount chain. These predictions are made with the following assumptions: 1) the Earth's magnetic field has always been primarily dipolar, and its poles have always remained stationary 2) the seamounts and ocean islands are created by static and localized hotspots below the plate (first hypothesized by [START_REF] Wilson | A possible origin of the hawaiian islands[END_REF]. With these hypotheses, it is possible to infer the paleolatitude of the plate from APW paths (the paleolongitude cannot be inferred because the Earth magnetic field is radially symmetrical) and the absolute motion of the plate by considering several seamount and island chains and estimating finite rotation series, similarly to what is done for the relative plate motions from magnetic anomalies and fracture zones. The respective uncertainties on observations can also be estimated (see for example Torsvik et al., 2008 for APW uncertainty estimations, and[START_REF] O'neill | On the uncertainties in hot spot reconstructions and the significance of moving hot spot reference frames[END_REF][START_REF] Wessel | Toward a self-consistent, high-resolution absolute plate motion model for the pacific[END_REF] for uncertainty estimations on absolute plate motions from seamount and island chain motions).

In practice, the estimation of absolute plate motions is more complex. First, the study of APW paths reveals a significant component of whole mantle and crust rotation, called true polar wander (e.g. [START_REF] Goldreich | Some remarks on polar wandering[END_REF][START_REF] Besse | Apparent and true polar wander and the geometry of the geomagnetic field over the last 200 myr[END_REF]Evans, FIGURE 2.8: Diagram of the relationship between relative plate motions (shown on the plate boundaries), absolute plate motions and predicted observations of apparent polar wander and seamount and island chain geometry, if the magnetic pole of the Earth and the hotspots were static.

Taken from [START_REF] Torsvik | Global plate motion frames: toward a unified model[END_REF].

2003; [START_REF] Maloof | Combined paleomagnetic, isotopic, and stratigraphic evidence for true polar wander from the neoproterozoic akademikerbreen group, svalbard, norway[END_REF]Steinberger & Torsvik, 2008). The true polar wander is due to the redistribution of masses in the convecting mantle, which changes the axis of maximum non-hydrostatic moment of inertia and leads to a whole mantle and crust rotation so as to realign this axis to the spin axis of the Earth. Second, the global study of seamount and ocean island chains shows that the hotspots are slightly moving through time, due to convective currents in the mantle (e.g. [START_REF] Molnar | Relative motions of hotspots in the pacific, atlantic and indian oceans since late cretaceous time[END_REF][START_REF] Tarduno | Large-scale motion between pacific and atlantic hotspots[END_REF][START_REF] Tarduno | Paleomagnetic evidence for motion of the hawaiian hotspot during formation of the emperor seamounts[END_REF][START_REF] Divenere | Are the pacific and indo-atlantic hotspots fixed? testing the plate circuit through antarctica[END_REF][START_REF] Tarduno | The emperor seamounts: Southward motion of the hawaiian hotspot plume in earth's mantle[END_REF]. As a result, the APW paths and seamount and ocean island chains actually give indications on the relative motions of plates with respect to the magnetic pole of the Earth and the moving mantle, respectively. To integrate these observations together with relative plate motions in a global model of plate tectonics reconstruction, it is necessary to define a common reference frame and relate all the data to it.

Reference frames vary depending on studies. [START_REF] Torsvik | Global plate motion frames: toward a unified model[END_REF] review and compare four types of reference frames generally used in plate tectonics reconstructions: paleomagnetic, African fixed hotspot, African moving hotspot, and global moving hotspot.

The paleomagnetic reference frame uses APW paths corrected from their estimated true polar wander to determine the variations in paleolatitudes of continents. Since the longitude is not constrained by APW paths, the reference frame also considers one plate as fixed in longitude (generally, Africa is chosen due to its central place in the Pangaea, and its supposed small amplitude of motion). The paleomagnetic reference frame is the only reference frame available before 130 Myr. The three other reference frames are variations on hotspot-based reference frames using different approximations. The African fixed hotspot reference frame considers that the hotspots in the Atlantic and Indian oceans are totally fixed (e.g. [START_REF] Müller | Revised plate motions relative to the hotspots from combined atlantic and indian ocean hotspot tracks[END_REF], and anchors relative plate motions in a frame using data coming from those hotspots. On the contrary, the African moving hotspot, and global moving hotspot reference frames estimate the motions of hotspots, and propose corrections to construct a reference frame. The estimation of the motion of hotspots is done using mantle circulation models (e.g. [START_REF] O'neill | On the uncertainties in hot spot reconstructions and the significance of moving hot spot reference frames[END_REF][START_REF] Doubrovine | Absolute plate motions in a reference frame defined by moving hot spots in the pacific, atlantic, and indian oceans[END_REF], that we will discuss in the next section.

Recent plate tectonics reconstructions use a hybrid reference frame, as proposed by [START_REF] Torsvik | Global plate motion frames: toward a unified model[END_REF], where a moving hotspot reference frame is used for recent times (up to 100 Ma) and is linked to a paleomagnetic reference frame for earlier times (e.g. [START_REF] Torsvik | Plate tectonics and net lithosphere rotation over the past 150my[END_REF][START_REF] Seton | Global continental and ocean basin reconstructions since 200ma[END_REF][START_REF] Müller | Ocean basin evolution and global-scale plate reorganization events since pangea breakup[END_REF].

Additionally, mantle structures estimated from tomography models have been included

as absolute plate motion constraints in some plate tectonics reconstruction models. [START_REF] Van Der Meer | Towards absolute plate motions constrained by lower-mantle slab remnants[END_REF] identified paleoslabs (remnants of subducted lithosphere in the mantle) using seismic tomography, and opened the way to the introduction of such features into plate tectonic reconstructions (e.g. [START_REF] Williams | Absolute plate motions since 130 ma constrained by subduction zone kinematics[END_REF]. [START_REF] Torsvik | Deep mantle structure as a reference frame for movements in and on the earth[END_REF] used the location of deep mantle structures (Large Low Shear Velocity Provinces, LLSVP), that are thought to be responsible for the production of Large igneous provinces, to estimate the longitude of continents in the Paleozoic.

The final result of these plate tectonics reconstructions is a time-series of plate layouts associated with pairs of rotation poles and angular velocities, which describe the history of the motions at the surface of the Earth. Modelling techniques and software developments have recently allowed the production of plate tectonics reconstructions continuous in time (e.g. [START_REF] Gurnis | Plate tectonic reconstructions with continuously closing plates[END_REF][START_REF] Müller | Ocean basin evolution and global-scale plate reorganization events since pangea breakup[END_REF], which made their integration into mantle convection models both easier and more consistent (see next section). Given the amount of observations used for these reconstructions, and the variety of specific techniques used to integrate those observations, the quantification of uncertainties in the final result is difficult to estimate. To our knowledge, a rigorous quantification of errors in the global plate tectonics reconstructions has not been performed yet. The progress made in successive plate reconstruction models are generally evaluated by comparing them with previous ones on their ability to explain observations [START_REF] Müller | Ocean basin evolution and global-scale plate reorganization events since pangea breakup[END_REF] and the plausibility of the dynamics they predict [START_REF] Williams | Absolute plate motions since 130 ma constrained by subduction zone kinematics[END_REF].

Conclusions on Observed Data

In this section, we described the data currently used in mantle circulation models: the temperature field at present and the history of surface kinematics. These data are not raw observations, but the result of models that integrate a wide range of observations. This setup is not ideal, since the addition of successive layers of models and/or inversions make the link between raw observations and the final result (the estimated mantle circulation) complex and the uncertainties difficult to assess. In this sense, the ultimate goal of data assimilation for mantle convection is to develop a framework able to integrate directly raw observations. The construction of community databases gathering raw data in a unified format [START_REF] Seton | Community infrastructure and repository for marine magnetic identifications[END_REF]Wessel et al., 2015;[START_REF] Moulik | 3-dimensional Reference Earth Model Webpage[END_REF], the progresses made in each of the disciplines to model and understand mantle structure and dynamics, and the techniques developed to assimilate data into models all contribute to this aim.

Dynamical Mantle Circulation Reconstructions

Three alternative strategies have been proposed to reconstruct the history of mantle circulation: backward advection, semi-empirical sequential methods and variational data assimilation. In the present dissertation, we propose a fourth one: sequential data assimilation, that we develop in Chapters 3 and 4. Backward advection and empirical sequential methods have already produced numerous estimations of global mantle circulation evolutions in the last decades. The variational data assimilation methods are more recent, and have so far been applied to regional reconstructions [START_REF] Liu | Simultaneous inversion of mantle properties and initial conditions using an adjoint of mantle convection[END_REF][START_REF] Spasojevic | Adjoint models of mantle convection with seismic, plate motion, and stratigraphic constraints: North america since the late cretaceous[END_REF] and just recently to global mantle circulation models [START_REF] Glišović | Reconstructing the cenozoic evolution of the mantle: Implications for mantle plume dynamics under the pacific and indian plates[END_REF][START_REF] Horbach | The adjoint method in geodynamics: derivation from a general operator formulation and application to the initial condition problem in a high resolution mantle circulation model[END_REF][START_REF] Glišović | A new back-and-forth iterative method for time-reversed convection modeling: Implications for the cenozoic evolution of 3-d structure and dynamics of the mantle[END_REF].

Backward advection and empirical sequential methods are direct methods: they integrate data in mantle convection models as boundary or initial conditions, without taking into account the uncertainties on the data and/or the model. We present these methods in the first part of this section. On the contrary, variational and sequential data assimilation consider mantle circulation reconstruction as an inverse problem, and explore (at least partially) the mantle state space to determine the evolution that best fits the observed data, given uncertainties on both data and the convection model.

We describe these techniques in the second part of this section.

Direct Methods

By the end of the 90's, the increase in computational power and the advances in numerical modelling had led to the emergence of numerical models able to solve the governing equations of mantle convection in 3D spherical geometry, with radial viscosity variations and at a Rayleigh number 1 to 2 orders of magnitude lower than that of the Earth's mantle [START_REF] Tackley | Effects of an endothermic phase transition at 670 km depth in a spherical model of convection in the earth's mantle[END_REF][START_REF] Bunge | Effect of depthdependent viscosity on the planform of mantle convection[END_REF]. Simultaneously, the resolution of mantle tomographic models had improved, revealing the distribution and geometry of subducted slabs [START_REF] Su | Simultaneous inversion for 3-d variations in shear and bulk velocity in the mantle[END_REF][START_REF] Masters | A shear-velocity model of the mantle [and discussion[END_REF][START_REF] Li | Global mantle shear velocity model developed using nonlinear asymptotic coupling theory[END_REF][START_REF] Grand | High resolution global tomography: a snapshot of convection in the earth[END_REF][START_REF] Van Der Hilst | Evidence for deep mantle circulation from global tomography[END_REF]. Finally, semi-dynamical models of mantle density evolution derived from surface kinematic history had shown their ability to reproduce the current geoid and seismic wave speed anomalies, and to predict plate velocities in the past (e.g. [START_REF] Ricard | A geodynamic model of mantle density heterogeneity[END_REF][START_REF] Zhang | Geoid anomalies from cenozoic subduction in semi-dynamical flow models including a phase boundary[END_REF][START_REF] Deparis | Mass anomalies due to subducted slabs and simulations of plate motion since 200 my[END_REF]Lithgow-Bertelloni & Richards, 1998). The Backward advection [START_REF] Steinberger | Changes of the earth's rotation axis owing to advection of mantle density heterogeneities[END_REF] and semi-empirical sequential [START_REF] Bunge | Time scales and heterogeneous structure in geodynamic earth models[END_REF] methods emerged from these progresses, as a way to investigate the dynamics of the mantle and its interaction with surface tectonics.

Backward advection

Backward advection consists in reversing the time in mantle convection equations, to compute the evolution of mantle circulation from its present state (derived from seismic tomography) backward in time, using plate tectonics reconstructions as surface boundary conditions. As described in section 1.3, the conditions at the upper boundary Γ a of the model are then ∀t, ∀x ∈ Γ a , u(x, t) • e r = 0 and u H (x, t) = u a H (x, t), (2.10) with u the total velocity and u H the horizontal velocity. t, x and e r are the time, the position vector and the unit radial vector. u a H is the surface velocity field given by plate tectonics reconstructions. In this formulation, it is necessary to have plate tectonics reconstructions which are continuous in time.

For mantle convection, the only prognostic equation is the equation of energy conservation. In the Boussinesq approximation, the nondimensional form of this equation is

∂T ∂t + u • ∇T = ∇ 2 T + R h , (2.11)
where T and R h are the nondimensional temperature and rate of energy production.

The backward heat advection-diffusion equation should be:

∂T ∂t -u • ∇T = -∇ 2 T -R h . (2.12)
This problem is an ill-posed problem: the backward diffusion term -∇ 2 T will trans- (2.13)

The backward advection method should be valid for a timespan which is short compared to the diffusion time of temperature heterogeneities in the mantle. If we consider a temperature anomaly of characteristic size 100 km (5 times smaller than the current resolution of global tomography models, see e.g. [START_REF] Ritsema | S40rts: a degree-40 shear-velocity model for the mantle from new rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements[END_REF], we obtain a characteristic diffusion time of

τ d = d 2 κ = 300 Myr, (2.14)
where κ = 10 -6 m 2 .s -1 is the coefficient of thermal diffusivity. [START_REF] Conrad | Seismic tomography, surface uplift, and the breakup of gondwanaland: Integrating mantle convection backwards in time[END_REF] add Furthermore, uncertainties on the state of the mantle at present (derived from seismic tomography) and on plate tectonics reconstructions add to the intrinsic time limit of retrodictions for the backward advection method.

Semi-empirical Sequential Methods

The semi-empirical sequential methods estimate mantle circulation by integrating plate tectonics reconstructions chronologically into a mantle convection model. Plate tectonics reconstructions are either introduced as velocity boundary conditions (following Equation 2.10), as first described by [START_REF] Bunge | Time scales and heterogeneous structure in geodynamic earth models[END_REF], or with a more sophisticated method, by blending a convection solution with thermal and kinematic models of plates and slabs [START_REF] Bower | Assimilating lithosphere and slab history in 4-d earth models[END_REF]. Contrary to backward advection, the timespan of the reconstruction is not limited by the assimilation method, but the availability of plate tectonics reconstructions. This makes possible the computation of mantle circulation models for much longer timespans (so far, up to 450 Myr, as done by [START_REF] Zhang | A model for the evolution of the earth's mantle structure since the early paleozoic[END_REF]. Such long term evolutions allow the investigation of the effect of surface tectonics on the deep mantle structure (e.g. [START_REF] Mcnamara | Thermochemical structures beneath africa and the pacific ocean[END_REF]Schuberth et al., 2009b,a;[START_REF] Zhang | Heat fluxes at the earth's surface and core-mantle boundary since pangea formation and their implications for the geomagnetic superchrons[END_REF][START_REF] Davies | Reconciling dynamic and seismic models of earth's lower mantle: The dominant role of thermal heterogeneity[END_REF][START_REF] Bower | Lower mantle structure from paleogeographically constrained dynamic earth models[END_REF][START_REF] Bull | The effect of plate motion history on the longevity of deep mantle heterogeneities[END_REF][START_REF] Hassan | Provenance of plumes in global convection models[END_REF], and the investigation of alternative scenarios of absolute plate motions beyond 70 Myr (e.g. [START_REF] Shephard | Testing absolute plate reference frames and the implications for the generation of geodynamic mantle heterogeneity structure[END_REF][START_REF] Nerlich | Mantle convection models constrain central neo-tethys ocean reconstructions[END_REF]Zahirovic et al., in press).

One important issue of semi-empirical sequential methods is that the initial condition Both semi-empirical sequential and backward advection methods have shown their sensitivity to variations in mantle convection models (e.g. [START_REF] O'neill | On the uncertainties in hot spot reconstructions and the significance of moving hot spot reference frames[END_REF][START_REF] Bello | Assessing the role of slab rheology in coupled plate-mantle convection models[END_REF] and in observed data (e.g. [START_REF] Doubrovine | Absolute plate motions in a reference frame defined by moving hot spots in the pacific, atlantic, and indian oceans[END_REF][START_REF] Shephard | Testing absolute plate reference frames and the implications for the generation of geodynamic mantle heterogeneity structure[END_REF].

The limits of these methods motivate the development of methods to assimilate data on mantle circulation into mantle convection models while taking into account the uncertainties on both the model and the observed data. We describe these techniques in the next section.

Data Assimilation Methods

In Section 2.1, we described the problem of mantle circulation reconstruction as as a time-dependent state estimation problem, following the classical notations of data assimilation. Two approaches exist to find a solution to this problem: sequential data assimilation and 4-D variational data assimilation.

So far, data assimilation methods for mantle circulation have been developed either as simplified schemes directly derived from 4-D variational data assimilation (e.g. [START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF], or as techniques derived from optimal control theory, bearing similarities with 4-D variational data assimilation (e.g. Ismail-Zadeh et al., 2003b,a, 2007). In this dissertation, we present two applications of sequential data assimilation methods to the mantle circulation problem. To ease the comparison between techniques, we keep the general formulation of data assimilation derived from the discretized state and model, as in Section 2.1. We first describe the common objective of data assimilation techniques. Then we describe the principle of 4-D variational data assimilation and the applications to mantle circulation. Finally, we describe the general principle of sequential data assimilation.

In the following, we consider only the second order statistics of the errors affecting the states and data. Other formulations exist, taking into account a more thorough description of error statistics (see e.g. [START_REF] Bocquet | Beyond gaussian statistical modeling in geophysical data assimilation[END_REF], for a review). More detailed descriptions of data assimilation techniques can be found in classical textbooks of atmospheric sciences and numerical weather prediction (e.g. [START_REF] Bengtsson | Dynamic meteorology: data assimilation methods[END_REF][START_REF] Daley | Atmospheric data analysis[END_REF][START_REF] Kalnay | Atmospheric modeling, data assimilation, and predictability[END_REF], oceanography (e.g. [START_REF] Chassignet | Ocean weather forecasting: an integrated view of oceanography[END_REF] or data assimilation in general (e.g. Evensen, 2009a;[START_REF] Swinbank | Data assimilation for the earth system[END_REF][START_REF] Blayo | Advanced data assimilation for geosciences[END_REF].

Objective of Data Assimilation

As described in Section 2.1, the unknown is defined as a time series of true state vectors {x t 1 , x t 2 , ..., x t K } and three sources of information are available, all affected by errors:

• the background state x b , which is related to the first true state by

x b = x t 1 + ǫ b (2.15)
with ǫ b the error on the background state. The error is modelled as a random vector of 0 mean (the background state is supposed unbiased) and covariance matrix P b .

• the time series of observed data {y o 1 , y o 2 , ..., y o K } and their associated observation operators {H 1 , H 2 , ..., H K }. At a given time step k the observed data vector y o k is related to the true state x t k by

y o k = H k (x t k ) + ǫ o k + ǫ h k (2.16)
with ǫ o k the measurement errors and ǫ h k the errors linked to the observation operator. The sum of errors ǫ o,h k = ǫ o k + ǫ h k is modelled as a random vector of 0 mean (observations are supposed unbiased) and covariance matrix R k .

• the dynamical model M, which relate a given state x k-1 at time step k -1 to the state vector x k at time step k by

x k = M(x k-1 ) + η k (2.17)
with η k the errors linked to modelling approximations. The error is modelled as a random vector of 0 mean (the model is supposed unbiased) and covariance matrix Q k .

We cannot compute directly the error made by a given evolution {x 1 , x 2 , ..., x K } with respect to the true state evolution {x t 1 , x t 2 , ..., x t K }, since it is unknown. However, we can compute the distances between the estimated evolution and the constraints given by the different sources of information. More precisely, we compute the Mahalanobis distances [START_REF] Mahalanobis | On the generalized distance in statistic[END_REF]. The Mahalanobis distance ∆ a (x) is a measure of the distance of a vector x to the multivariate continuous random vector a of mean µ and covariance matrix Σ:

∆ a (x) = (x -µ) T (Σ) -1 (x -µ).
(2.18)

For example, when the covariance matrix Σ is the identity matrix, this means that the variance of each component of the random vector a is one and the components of a are uncorrelated. In this case, the Mahalanobis distance corresponds to the Euclidean distance between the vector x and the mean µ (i.e. the L2 norm of xµ).

For our problem, we compute the following Mahalanobis distances:

∆ ǫ b (x 1 -x b ) 2 = (x 1 -x b ) T P b -1 (x 1 -x b ) (2.19a) ∀k ∈ {1, .., K}, ∆ ǫ o,h k (H k (x k ) -y o k ) 2 = (H k (x k ) -y o k ) T (R k ) -1 (H k (x k ) -y o k ) (2.19b) ∀k ∈ {2, .., K}, ∆ η k (M(x k-1 ) -x k ) 2 = (M(x k-1 ) -x k ) T (Q) -1 (M(x k-1 ) -x k ). (2.19c)
The global aim of data assimilation methods is to find a state evolution which minimizes these distances. 4-D variational and sequential data assimilation methods differ by the strategies adopted to determine this evolution.

4-D Variational Data Assimilation

4-D variational data assimilation methods aim at providing a time series of the estimated evolution of the state of the system such that the objective function J defined by

J (x 1 , x 2 , ..., x K ) =(x 1 -x b ) T P b -1 (x 1 -x b ) + K k=1 (H k (x k ) -y o k ) T (R k ) -1 (H k (x k ) -y o k ) + K k=2 (x k -M(x k-1 )) T (Q k ) -1 (x k -M(x k-1 )) (2.20)
is minimized. The minimization of this function is called weak constraint four dimensional variational assimilation or weak constraint 4D-var, following the terminology of [START_REF] Sasaki | Some basic formalisms in numerical variational analysis[END_REF]. The term "weak constraint" comes from the fact that we take into account modelling errors.

J (x 1 , x 2 , ..., x K ) =(x 1 -x b ) T P b -1 (x 1 -x b ) + (H K (x K ) -y o K ) T (R K ) -1 (H K (x K ) -y o K ) + K k=2 (x k -M(x k-1 )) T (Q k ) -1 (x k -M(x k-1 )) (2.21)
The second difference is the mechanical boundary conditions at the surface of the model. Since the mantle convection model cannot produce plate-like tectonics at the surface self-consistently, plate tectonics reconstructions have to be imposed as surface boundary conditions (Equation 2.10) to obtain an "Earth-like" evolution. In this formulation, the mantle model M includes the plate tectonics reconstructions.

In this specific case, [START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF] derives Euler-Lagrange equations associated to the continuous objective function. In principle, the solution to the Euler-Lagrange equations provides the state evolution which minimizes the objective function. However, in practice, further assumptions are made before computing the solution to this problem. [START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF] suppose no error on the model (which means in their formulation also no error on plate tectonics reconstructions). If we neglect modelling errors, then the evolution of the state of the system can be determined solely from the first state:

∀k ∈ {2, ..., K}, x k = M(x k-1 ) (2.22)
Under this assumption, the objective function J depends only on the initial state x 1 and becomes

J (x 1 ) =(x 1 -x b ) T P b -1 (x 1 -x b ) + (H K (x K ) -y o K ) T (R K ) -1 (H K (x K ) -y o K ).
(2.23)

We consider further that the covariance matrices associated to the background state error P b and the observed data error R K are of the form:

P b = V b I and R K = V o I (2.24)
where I is the identity matrix and V b and V o are the variances associated to each component of the background state and observation error random vectors.

We obtain the cost function

J (x 1 ) = 1 V b (x 1 -x b ) T (x 1 -x b ) + 1 V o (H K (x K ) -y o K ) T (H K (x K ) -y o K ) .
(2.25)

In the continuous form, and considering that the states and data contain only temperature fields, this translates to

J c (T i ) = 1 V b Ω T i (a) -T b (a) 2 da + 1 V o Ω [T f (a) -T o (a)] 2 da. (2.26)
Where T i is the initial temperature field to be determined, T f the final temperature, computed from the evolution of T i using the forward model of mantle convection.

T b is the a priori initial temperature field, and T o is the temperature field at present, derived from seismic tomography. Ω represents the volume of the mantle and a is the position vector.

The function J c , (or a further simplification, taking V b → ∞ and V o = 1), along with a mantle convection model, form the basis of the development of current data assimilation methods applied to the mantle circulation inverse problem. In practice, iterative methods for the minimization of the function J c are developed in the continuous form, and the resulting equations are then discretized to solve the system.

The general organization of algorithms is:

1. Initialization with a discretized temperature field T

(1) i .

2. Iteration until convergence of the solution. At the nth iteration (a) Computation of the state evolution using the forward code of mantle convection, from a given discretized initial state T

(n) i

we obtain a final state

T (n) f .
(b) Computation of the correction δT (n) to add to the initial state to minimize the cost function. The new initial temperature field is then T

(n+1) i = T (n) i + δT (n)
Algorithms differ by the techniques used to estimate the correction δT (n) in stage 2b.

One solution, adopted by [START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF], and further investigated by Horbach et al. (2008) perform synthetic regional tests using this method and show that the result of the assimilation depends strongly on the chosen initial state in stage 1 of the algorithm.

To obtain convergence towards the reference evolution, it is necessary to compute the initial temperature field (stage 1 of the algorithm) from a simple backward advection, as done in section 2.3.1.

Finally, [START_REF] Ismail-Zadeh | Quasi-reversibility method for data assimilation in models of mantle dynamics[END_REF] propose an algorithm based on the regularization of the backward heat equation (Equation 2.12) for stage 2b. This method is called the quasi-reversibility method [START_REF] Lattés | The method of quasi-reversibility: applications to partial differential equations[END_REF]. With synthetic tests, Ismail-Zadeh et al. ( 2007) finds a lower accuracy of the quasi-reversibility method compared to adjoint methods. However, the quasi-reversibility method is computationally lighter (see Ismail-Zadeh et al., 2016, Chapter 7 for a comparison between backward advection, adjoint and quasi-reversibility methods). [START_REF] Glišović | Reconstructing the cenozoic evolution of the mantle: Implications for mantle plume dynamics under the pacific and indian plates[END_REF], 2016) recently applied the quasi-reversibility method to the reconstruction of mantle circulation for the last 65 Myr.

All the techniques described previously, when applied to global reconstructions of mantle circulation, impose plate tectonic reconstructions as mechanical boundary conditions, without taking into account the uncertainties on these reconstructions. As

shown for example by the synthetic experiments of [START_REF] Vynnytska | Restoring past mantle convection structure through fluid dynamic inverse theory: regularisation through surface velocity boundary conditions[END_REF], the imposition of surface constraints during the assimilation is decisive to obtain solutions that converge to the true evolution (also highlighted by [START_REF] Glišović | A new back-and-forth iterative method for time-reversed convection modeling: Implications for the cenozoic evolution of 3-d structure and dynamics of the mantle[END_REF].

Given the uncertainties of plate tectonics reconstructions as we go back in time (see section 2.2.2), it then appears important to be able to take them into account in the data assimilation algorithms developed.

In this dissertation, we focus on the development of data assimilation methods that consider the history of plate tectonics not as a boundary condition that is imposed at the surface of the model, but as observed data, which are assimilated in a consistent fashion in mantle circulation models. To do so, we use direct models of mantle convection that are able to produce plate-like tectonics at their surface self-consistently.

These models have viscosity laws that strongly vary with both temperature and stress, as is discussed in Section 1.4. This introduces strong nonlinearities in the flow equations, which complicates the development of 4-D variational data assimilation methods. Hence, we turn to sequential data assimilation methods, the principle of which we describe in the next section.

Sequential Data Assimilation

The exact formulation of the sequential data assimilation algorithms that were developed during this thesis are described in the next two chapters. Here, we simply highlight the principle of the method. Although sequential methods find their roots in estimation theory more than optimal control theory (see the following Chapters), we present here the different estimations as the minimization of objective functions, so as to ease the comparison with the previous 4-D variational data assimilation methods.

Sequential methods divide the problem into several consecutive stages:

1. The best guess for the initial state x 1 is computed by minimizing the objective function:

J (x 1 ) =(x 1 -x b ) T P b -1 (x 1 -x b ) + (H 1 (x 1 ) -y o 1 ) T (R 1 ) -1 (H 1 (x 1 ) -y o 1 ) (2.27)
the result of this minimization is called the first analyzed state x a 1 , which is linked to the first true state by x a 1 = x t 1 + ǫ a 1 , with ǫ a 1 the error on the first analyzed state.

ǫ a 1 is modelled as a random vector, of mean 0 and covariance matrix P a 1 (called the analyzed error covariance matrix). P a 1 is also evaluated during this stage by taking into account the errors on both the background state and the first observed data.

2. Two steps are repeated, until all observations have been taken into account. At timestep k ∈ {2, 3, ..., K}:

(a) the Forecast state x f k is computed using the model M:

x f k = M(x a k-1 ).
(2.28)

The forecast state is linked to the kth true state by x f k = x t k + ǫ f k , with ǫ f k the error on the kth forecast state. ǫ f k is modelled as a random vector, of mean 0 and covariance matrix P f k (called the forecast error covariance matrix). P f k is also evaluated during this stage, by taking into account the errors on the analyzed state and the modelling errors.

(b) the Analyzed state x a k is determined by minimizing the objective function

J (x k ) =(x k -x f k ) T P f k -1 (x k -x f k ) + (H 1 (x k ) -y o k ) T (R k ) -1 (H k (x k ) -y o k ).
(2.29)

x a k is linked to the first true state by x a k = x t k + ǫ a k , with ǫ a k the error on the kth analyzed state. ǫ a k is modelled as a random vector, of mean 0 and covariance matrix P a k (called the analyzed error covariance matrix). P a k is also evaluated during this stage, by taking into account errors on the kth forecast state and errors on the kth observed data.

In this scheme, the observations are integrated sequentially: the evolution of the estimated state is computed forward in time, and the observed data, as they become available, are assimilated into the estimated state evolution. The estimated state evolution should then converge towards the true evolution, as more and more data are assimilated.

The following chapters 3 and 4 propose two algorithms based on this principle to reconstruct mantle circulation. We performed synthetic experiments on both these algorithms to determine whether it is possible to reconstruct mantle circulation from the sole and imperfect knowledge of the surface tectonics history, and under which conditions a satisfactory reconstruction is possible. The first one is a suboptimal scheme derived from the classical Kalman Filter (see e.g. Todling & Cohn, 1994, for a review on suboptimal schemes based on Kalman filters). This scheme provided encouraging results and drove us to develop a more efficient sequential data assimilation algorithm for the mantle circulation problem. This is what we do in Chapter 4 with the implementation of an Ensemble Kalman Filter for this problem.

Résumé du chapitre

Le but d'une reconstitution de la circulation mantellique est d'estimer l'écoulement 3-D dans le manteau, au présent et dans le passé. Pour atteindre ce but, on combine des données sur la circulation mantellique avec un modèle dynamique de convection mantellique. On dispose de données sur la circulation mantellique au présent, avec par exemple les contraintes fournies par la tomographie sismique, et dans le passé, avec par exemple les contraintes sur l'histoire de la tectonique de surface. Les problèmes inverses s'intéressant à des données asynchrones appartiennent à la catégorie des problèmes d'assimilation de données. Ce chapitre vise à présenter le problème de l'assimilation de données pour la circulation mantellique. Nous commençons par formuler le problème de la circulation mantellique dans le cadre de l'assimilation de données. Nous décrivons ensuite les données actuellement utilisées dans les modèles de circulation mantellique: le champ de température actuel et l'histoire de la cinématique de surface. Finalement on présente les applications à la reconstruction de la circulation mantellique. Trois stratégies alternatives ont été proposées pour reconstruire l'histoire de la circulation mantellique: l'advection rétrograde, les méthodes séquentielles semi-empiriques et l'assimilation de données variationnelle. Dans la présente thèse, nous proposons une quatrième méthode: l'assimilation de données séquentielle, que nous développons dans chapitres 3 et 4. Ces quatre méthodes sont présentées dans la troisième section de ce chapitre.

CHAPTER 3

A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics

The following Chapter has been published in january 2016 as:

Introduction

Global tectonics is a surface expression of mantle convection [START_REF] Bercovici | The generation of plate tectonics from mantle convection[END_REF]: the motions of continents and seafloor are generated by forces acting within the mantle and the lithosphere. For example, [START_REF] Ricard | Mantle heterogeneities, geoid, and plate motion: A monte carlo inversion[END_REF] and [START_REF] Alisic | Multi-scale dynamics and rheology of mantle flow with plates[END_REF] obtained a consistent description of surface kinematics by converting the long wavelength heterogeneities of seismic velocity into buoyancy forces. Mantle convection studies also attest to this link, showing for example that a downwelling in a context of large-scale convection draws continents to aggregate [START_REF] Zhong | Role of ocean-continent contrast and continental keels on plate motion, net rotation of lithosphere, and the geoid[END_REF][START_REF] Rolf | Statistical cyclicity of the supercontinent cycle[END_REF].

Reconstructing the convective history of the Earth's mantle is a long-standing problem in geosciences [START_REF] Bunge | Time scales and heterogeneous structure in geodynamic earth models[END_REF]. Indeed, an accurate determination of mantle temperature and velocity fields evolution only for the last 200 My would already lead to a better understanding of ore deposits formation, water resources distribution, the evolution of the geodynamo on geological time scales or deep material properties. The geodynamics community has taken advantage of the link between mantle convection and tectonics to build a first approach to reconstruct the convection history of the past 100 My. For this approach, plate tectonics theory is used to describe surface kinematics.

By driving convection currents at the surface with plate kinematic reconstructions, it is possible to propose a scenario of the evolution of mantle heterogeneities. Results so obtained proved to be relatively consistent with long wavelength seismic tomography or/and the geoid [START_REF] Bunge | Mesozoic plate-motion history below the northeast pacific ocean from seismic images of the subducted farallon slab[END_REF][START_REF] Flament | A review of observations and models of dynamic topography[END_REF]. These 3D spherical convection models are called mantle circulation models [START_REF] Bunge | Mantle-circulation models with sequential data assimilation: inferring present-day mantle structure from platemotion histories[END_REF], as a reference to oceanography. They require very fine parametrization in order to generate internal structures consistent with both the physics of convection and geophysical observations [START_REF] Bower | Assimilating lithosphere and slab history in 4-d earth models[END_REF]. To go one step further, several groups used tomographic models as data, in addition to plate reconstructions [START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF][START_REF] Ismail-Zadeh | Quasi-reversibility method for data assimilation in models of mantle dynamics[END_REF][START_REF] Liu | Simultaneous inversion of mantle properties and initial conditions using an adjoint of mantle convection[END_REF]. A temperature field of the present-day Earth's mantle is generated from these tomographic models, and the past mantle circulation is retrieved following various methods such as backward advection, variational methods and quasi-reversibility methods (see Ismail-Zadeh & Tackley (2010) for a full report on mantle convection and surface tectonics these methods). There are two drawbacks in this methodology, both difficult to overcome. First, the convection models used so far do not naturally produce Earth-like tectonics: they generate either sluggish tectonics at the surface or a stagnant lid when larger viscosity contrasts are prescribed [START_REF] Solomatov | Scaling of temperature-and stress-dependent viscosity convection[END_REF]. Second, although tectonic data are the most important source of information (they are time-dependent, unlike seismic data), they are used only to drive the flow, not as data to match. To drive the flow, the kinematics needs to be known at every time-step of the calculation. This is not the case for plate tectonic reconstructions, even though [START_REF] Gurnis | Plate tectonic reconstructions with continuously closing plates[END_REF] designed a procedure to interpolate the plate tectonics geometry between times at which it is known.

Significant progress has been made in the past 10 years on modelling convection that produces more realistic surface tectonics. Convection with a pseudo-plastic rheology generates surface tectonics with a plate-like behaviour [START_REF] Moresi | Mantle convection with a brittle lithosphere: thoughts on the global tectonic styles of the earth and venus[END_REF][START_REF] Moresi | Plate tectonics and convection in the earth's mantle: Toward a numerical simulation[END_REF][START_REF] Stein | The effect of rheological parameters on plate behaviour in a self-consistent model of mantle convection[END_REF][START_REF] Bercovici | The generation of plate tectonics from mantle convection[END_REF]Tackley, 2000b;[START_REF] Van Heck | Planforms of self-consistently generated plates in 3d spherical geometry[END_REF][START_REF] Bercovici | Plate tectonics, damage and inheritance[END_REF]. Recent models display seafloor spreading and continental drift comparable to that of the Earth to first order: seafloor age distributions and the time scale of spreading fluctuations are consistent with what has been inferred for the Earth for the last 200 My [START_REF] Coltice | Dynamic causes of the relation between area and age of the ocean floor[END_REF][START_REF] Coltice | Convergence of tectonic reconstructions and mantle convection models for significant fluctuations in seafloor spreading[END_REF]. This opens the way to producing a mantle circulation model using convection models with plate-like behaviour, under the observational constraint provided by tectonic data. This constraint is not to be enforced through kinematic boundary conditions imposed at every time step. Instead, it should be taken into account in a statistically consistent fashion, whereby the prediction due to the model and the observations are combined in a way that respects the uncertainties affecting both.

We explore here this possibility using a controlled case, and present a proof of concept. We have developed a sequential data assimilation method which is a suboptimal scheme based on the Kalman filter (see Todling & Cohn (1994) for a review of suboptimal schemes based on the Kalman filter). Such data assimilation methodology provides the best linear unbiased estimate of the temperature field and surface velocities at times data are available, and a prediction of the evolution of the system state between those times. Using synthetic tests in spherical annulus geometry [START_REF] Hernlund | Modeling mantle convection in the spherical annulus[END_REF], we show that our data assimilation method is efficient for reconstructions over 200 My with data gaps of 5 -10 My and relative uncertainties on surface observations below 30%. This method also allows the refinement of surface tectonics, taking advantage of the information within the physics of convection.

In the following, we first describe the general approach to Kalman Filtering in Section 2. Section 3 provides details on the definition of the different operators and vectors forming the backbone of a sequential data assimilation framework: the physical model, the data and the state for mantle convection, as well as the observation operator. Section 4 presents our sequential data assimilation algorithm and the adjustments used for the reconstruction of mantle convection. Then, we present in Section 5 an evaluation of the method with a series of synthetic experiments, which leads to the discussion of section 6.

The Extended Kalman Filter

We introduce in this section the Extended Kalman Filter that we adapt and apply to the convection reconstruction problem. We use here classical notations, taken from [START_REF] Ide | Unified notation for data assimilation: operational, sequential and variational[END_REF]. Our goal is to estimate the evolution of the mantle temperature field, using information from a mantle convection model and tectonic data including global kinematics and seafloor ages. The following is a rather general description of the Extended Kalman Filter; similar developments can be found in [START_REF] Ghil | Data assimilation in meteorology and oceanography[END_REF]; [START_REF] Kalnay | Atmospheric modeling, data assimilation, and predictability[END_REF]; [START_REF] Wunsch | Discrete inverse and state estimation problems: with geophysical fluid applications[END_REF] for meteorology and oceanography and in Fournier et al.

(2010) for geomagnetism, for example.

Data assimilation aims at estimating the state of a dynamical system and its evolution for a given period of time, combining information from observations and a physical model. The evolution of the state of the system is described by the set of true state Chapter 3. A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics vectors, storing the chosen scalar fields that are discretized:

x t 1 , x t 2 , ..., x t nt , (3.1)
the subscripts {1, 2, ..., n t } being the discrete time steps.

Likewise, the observations are defined as a set of column vectors

y o t 1 , y o t 2 , ..., y o tn with {t 1 , t 2 , ..., t n } ⊂ {1, 2, ..., n t }. (3.2)
The true state and the observations are linked linearly by the observation operator H

∀t i ∈ {t 1 , t 2 , ..., t n }, y o t i = Hx t t i + ǫ o t i , (3.3)
where ǫ o t i is the observation error at the timestep t i . The precise value of ǫ o t i is unknown. However, it is possible to model it as a random vector and estimate its probability density function (PDF). We assume ǫ o t i follows a centred Gaussian statistics, and its associated covariance matrix

R t i = Covar(ǫ o t i ) = ǫ o t i (ǫ o t i ) T (3.4)
is called the observation error covariance matrix. The operator [] T means transpose.

The Extended Kalman Filter belongs to the Sequential Data Assimilation techniques.

It consists of integrating the observation vectors one after the other into a numerical model, when they become available. Fig. 3.1 gives an overview of this procedure. The initialization, at timestep t 1 , provides an a priori guess of the state of the system. Then two steps called analysis and forecast are performed sequentially, every time a new observation vector becomes available, until all the observations have been taken into account. Hence the name of the technique: sequential data assimilation. The analysis corrects the current estimate of the system state by considering the new observations.

The forecast provides an estimation of the evolution of the system until the next time new observations are available. The evolution of the estimated state should converge x a t 1 (blue dot at t 1 ), which is closer to the true state x t t 1 (red cross at t 1 ). The evolution of the state (green line) is then computed by the model M until t 2 , leading to a new forecast state x f t 2 (green dot). A new analyzed state x a t 2 is computed considering new observations. The sequence is repeated until all observations have been taken into account. The aim of data assimilation is to get as close as possible to the true state evolution (red line).

mantle convection and surface tectonics towards the true state as more observations are assimilated.

Initialization

The Extended Kalman Filter algorithm starts at the timestep t 1 , when the first observations are available. Before taking into account the first observations, an a priori state of the system has to be estimated. It is called the background state x b . The relationship between the background state and the true state at t 1 is

x b = x t t 1 + ǫ b , (3.5)
where ǫ b is the background state error. ǫ b is a random vector with ǫ b = 0 and

Covar(ǫ b ) = P b .
The background state vector is considered as the forecast state vector at timestep t 1 , written x f t 1 and its covariance matrix P b is the forecast error covariance matrix at t 1 , written P f t 1 .

Analysis and forecast sequence

At t i , a new estimate of the state of the system is calculated considering y o t i , R t i , x f t i and P f t i . This step is called analysis, the new estimate being the analyzed state x a t i . The relationship between the analyzed state and the true state is

x a t i = x t t i + ǫ a t i , (3.6)
where ǫ a t i is the analysis error at t i . ǫ a t i is modelled as a random vector following a Gaussian law, and its covariance P a t i is the analysis error covariance matrix. The analyzed state x a t i is chosen so that • x a t i is a linear combination of x f t i and y o t i ,

• Tr(P a n ) is minimised,

• ǫ a n = 0.

These conditions lead to the classical Best Linear Unbiased Estimate (see [START_REF] Ghil | Data assimilation in meteorology and oceanography[END_REF] for example)

x a t i = x f t i + K t i y o t i -Hx f t i , (3.7)
K t i = (HP f t i ) T H(HP f t i ) T + R -1
, (3.8)

P a t i = (I -K t i H)P f t i , (3.9)
where K t i is the Kalman gain at t i . If the probability density functions of the errors on the forecast state and the observed data are Gaussian and the observation operator is linear, then x a t i is not only the estimate of minimum variance but as well the most likely state: in this case, the analysis will be optimal.

The evolution of the estimated state of the system from t i to t i+1 , the next timestep observations are available, is computed using the direct numerical model M:

x f t i+1 = M x a t i + η t i .
(3.10)

Its associated covariance matrix is

P f t i+1 = MP a t i M T + Q t i , (3.11)
where M is the tangent linear model, i.e. the linearised version of the operator M and Q t i is the model error covariance matrix. Equations 3.10 and 3.11 describe the forecast step. The estimated state of the system at a timestep t i+1 x f t i+1 is called the forecast state and P f t i+1 is the forecast error covariance matrix at timestep t i+1 . P f t i+1 takes into account both the propagation of errors of the analyzed state at timestep t i , ǫ a t i and the model errors η t i . We assume a perfect model in this study, so that ∀t i , η t i = 0.

(3.12)

mantle convection and surface tectonics

Convection model, Mantle State Vector and Tectonic

Data

In this section, we describe in detail what M, H, x and y refer to in our mantle convection context.

Convection Model with Plate-Like Behaviour

The forward model M is our source of prior information. It solves the equations of conservation of mass, momentum and energy with classical simplifications for mantle convection: infinite Prandtl number and Boussinesq approximation. We further assume an isochemical mantle, and non-dimensionalize the equations to thermal diffusion scales (for a full development of the equations, see [START_REF] Ricard | Physics of mantle convection[END_REF] for example).

We obtain The dynamic viscosity µ varies with temperature and stress following the equation

∇ • u = 0, ( 3
µ = µ -1 T + µ -1 y -1 , (3.16)
µ T decreasing exponentially with temperature (according to Arrhenius law), and divided by β when reaching the solidus:

µ T = exp E A T + T 1 - E A 2T 1 if T < T s , (3.17) µ T = β -1 exp E A T + T 1 - E A 2T 1 if T > T s . (3.18)
T 1 is the temperature for which µ T = 1, E A is the activation energy and T s = T s 0 + ∇ r T s (r a -r) with r a the surface value of r. T s models the variation of solidus with depth and is tuned so that the viscosity drop is located at the base of the top boundary layer. This results in a weaker asthenosphere and favours plate-like behaviour [START_REF] Richards | Role of a lowviscosity zone in stabilizing plate tectonics: Implications for comparative terrestrial planetology[END_REF]Tackley, 2000b).

µ y is defined by

σ yield = σ Y + (r a -r).∇ r σ Y , µ y = σ yield 2 ǫ , (3.19)
with σ Y , ∇ r σ Y and ǫ being the yield stress at the surface, the depth-dependence of the yield stress and the second invariant of the strain rate tensor respectively. The strain rate tensor is linked to the velocity by

ǫ = 1 2 ∇u + (∇u) T . (3.20)
Solutions are computed using StagYY [START_REF] Tackley | Effects of an endothermic phase transition at 670 km depth in a spherical model of convection in the earth's mantle[END_REF], a finite-volume, multigrid convection code. We use a spherical annulus grid which provides results closer to the spherical grid than a cylindrical geometry [START_REF] Hernlund | Modeling mantle convection in the spherical annulus[END_REF]. The grid is refined in the radial direction near the upper boundary of the model. In the following, the longitudinal coordinate of a point is written φ m , with m ∈ {1, 2, ..., M } and its radial coordinate is written r n with n ∈ {1, 2, ..., N }, r varying from r b to r a . The value of the parameters used for this work are given in Table 3.1.

Since this work is a proof of concept, we chose a fairly simple model, with equation governing the flow relying on strong assumptions (incompressible, isochemical mantle). The Rayleigh number is 10 6 which is one or two orders of magnitude lower than that of the Earth, but high enough to ensure chaotic convection. The viscosity law selfconsistently generates plates at the surface, as it can be seen from the surface velocity in Fig. 3.10. Ignoring many more complexities does not mean they are not fundamental for reconstructing Earth mantle evolution, but we focus in this manuscript on the data assimilation methodology.

The State of the mantle

Given Equations 3.13 to 3.15 above, the field variables describing the state of the Earth's mantle x are velocity, temperature and pressure. Inertial forces are negligible for mantle convection, which means that the velocity and pressure fields can be deduced from the temperature field at any time given the viscosity law and other parameters, as

shown in the diagnostic Equations 3.13 and 3.14. Since it is possible to compute u and p from the sole knowledge of T , x should be restricted to T alone.

The Data: Surface Heat Flux and Surface velocities

As a first approach, the data y o we use are not direct measurements per se, but plate reconstruction models. For instance, [START_REF] Seton | Global continental and ocean basin reconstructions since 200ma[END_REF] or [START_REF] Shephard | The tectonic evolution of the arctic since pangea breakup: Integrating constraints from surface geology and geophysics with mantle structure[END_REF] proposed plate tectonics reconstructions for the last 200 My, using the continuously closing plates methodology [START_REF] Gurnis | Plate tectonic reconstructions with continuously closing plates[END_REF]) so that a reconstruction can be numerically computed for any time between a given t i and the subsequent t i+1 . These reconstructions integrate paleomagnetic, paleobiological and geological data to provide continuous maps of surface velocity and seafloor age as well as the position and geometry of continents. It is this type of data that is used today in convection reconstructions with imposed boundary conditions. One fundamental difference between these methods and our sequential assimilation method is that the latter naturally takes into account uncertainties in the reconstructions. A second difference is that we do not need the surface data to be known at all times.

Plate reconstructions provide estimates of the velocity at any location on the surface of the Earth in the approximation of the plate tectonics theory, as well as the age of the seafloor. In the model we use in this manuscript, surface heat flux is an excellent proxy for the age of the seafloor [START_REF] Coltice | Dynamic causes of the relation between area and age of the ocean floor[END_REF]. Consequently, we propose to consider surface heat flux and surface velocity as the data to assimilate. However, with more sophisticated models, small scale convection would require an explicit computation of the age of the seafloor.

At each time t i , the data vector for the present study will be y = [q s (φ 1 ), . . . , q s (φ N ), V s (φ 1 ), . . . , V s (φ N )] T , (3.21)

where q s (φ m ) and V s (φ m ) are the surface heat flux and tangential velocity values at longitude φ m , respectively. Hence, for the model parameters set of Table 3.1, the data vector contains 2N = 768 values.

Chapter 3. A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics

The Observation Operator and the Augmented State

As described by Equation 3.3, the link between the state of the system x and the data y has to be expressed in the form of a linear equation. This is straightforward for surface heat flux. At a given longitude φ m , the surface heat flux is approximated by a first order discretization of Fourier's law

∀φ m q(φ m ) = -k T (φ m , r N ) -T (φ m , r N -1 ) r N -r N -1 , (3.22)
where k is the thermal conductivity, which is 1 for our non-dimensional framework.

However, the link between the velocity field and the temperature field is highly nonlinear, because of our choice of rheology. This is why we consider an augmented state of the mantle, containing both temperature and surface velocity fields:

x =[T (φ 1 , r 1 ) , . . . , T (φ M , r 1 ), T (φ 1 , r 2 ), . . . , (3.23)
T (φ M , r N ), V (φ 1 ), . . . , V (φ M )],

(3.24)

where T (φ m , r n ) is the temperature value at longitude φ m and radius r n and V (φ m ) is the surface tangential velocity value at longitude φ m . For the model parameters set of Table 3.1, the state vector contains N M + N = 18816 values.

Sequential Data Assimilation Algorithm for Mantle Convection

We describe here the steps of our sequential algorithm, which are

• an initialization, which evaluates the background state x b and its error covariance matrix P b ,

• an analysis which estimates

x f t i , P f t i , y o t i , R o t i → x a t i ,
• a forecast, which computes

x a t i → x f t i+1 , P f t i+1 .

Initialization

The background state x b and its associated covariance matrix P b describe an estimate of the state of the system when no observation is yet available. x b components correspond to the average temperature field and surface velocities of the mantle. P b contains the two-point spatial correlation of the temperature field and surface velocity field, as described by [START_REF] Balachandar | Eigendecomposition of the two-point correlation tensor for optimal characterization of mantle convection[END_REF].

To estimate x b and P b , we use a set of K = 200 states of the system x k with k ∈ {1, ..., K}, computed with the mantle convection code. These states are extracted from a free run of the convection model (i.e. unconstrained by data). To ensure decorrelation between the states, the time between two snapshots is equivalent to 5 Lyapunov times.

The Lyapunov time is the e-folding time for the growth of an initial error in a dynamic system. In the calculations we describe below, we obtain a typical Lyapunov time of 140 My, similar to the Lyapunov times computed in 3D spherical convection models [START_REF] Bello | On the predictability limit of convection models of the earth's mantle[END_REF].

The background state is

x b = x = 1 K K k=1 x k (3.25)
and the error covariance matrix is

P b = Covar(x) = 1 K -1 K k=1 x k -x b x k -x b T .
(3.26)

The probability density functions of all the variables are also estimated, to compare them to Gaussian probability density functions. Indeed, the best linear unbiased estimate used during analysis is optimal only in the linear Gaussian case, as described in section 2.

As noted before, the initial setting of the model is spherically symmetric. As a consequence ∀(φ, r), T (φ, r) = T (0, r) , (3.27) x b , P b , y b and R b are used to standardize the state x and data y for the analysis step.

∀(φ 1 , φ 2 , r 1 , r 2 ), Cov(T (φ 1 , r 1 ), T (φ 2 , r 2 )) = Cov(T (0, r 1 ), T (φ 1 -φ 2 , r 2 )), (3.28)
∀(φ 1 , φ 2 , r 1 , r 2 ), Cov(T (φ 1 , r 1 ), T (φ 2 , r 2 )) = Cov(T (0, r 1 ), T (φ 2 -φ 1 , r 2 )). (3.29) Likewise ∀φ, V s (φ) = V s (0) , (3.30) ∀(φ 1 , φ 2 ), Cov(V s (φ 1 ), V s (φ 2 )) = Cov(V s (0), V s (φ 1 -φ 2 )), (3.31) ∀(φ 1 , φ 2 , r 1 ), Cov(T (φ 1 , r 1 ), V s (φ 2 )) = Cov(T (0, r 1 ), V s (φ 2 -φ 1 )), (3.32) Cov(T (φ 1 , r 1 ), V s (φ 2 )) = -Cov(T (0, r 1 ), V s (φ 1 -φ 2 )). ( 3
We define their corresponding standard score state x and standard score data y by Hence, the two points spatial temperature-velocity correlations can be mapped onto a spherical annulus. Fig. 3.5 represents this correlation map. Fig. 3.5 shows that surface velocities are weakly correlated with temperature, with a maximum correlation value of 0.4. Nevertheless, the correlations between surface velocity and temperature at the bottom are not negligible, and this is mainly due to the sinking of slabs to the bottom of the model.

x = S b -1 x -x b , ( 3
Although we ignore the non-Gaussianity of some distributions, second order statistics succeed in summing up the main features of the convection model dynamics.

Furthermore, the matrix P b is eigendecomposed leading to

P b = V ΛV T . (3.45)
where V contains the eigenvectors of P b and Λ its eigenvalues. The background correlation matrix P b is reduced to

P b r = V r Λ r V T r (3.46)
where Λ r contains only the 1928 largest eigenvalues, which account for 99.98% of the cumulative variance of P b . V r is composed of the corresponding eigenvectors. The reduction of the correlation matrix is equivalent to assuming that the correlations associated with the lowest eigenvalues and corresponding eigenvectors are not meaningful.

This technique, described by [START_REF] Cane | Mapping tropical pacific sea level: Data assimilation via a reduced state space kalman filter[END_REF], not only lightens the computational cost, but also corrects the state of the system only in those directions followed by the dynamical model.

We can now describe the first standard score forecast state as

x f t 1 = 0, (3.47)
P f t 1 = P b r = V r Λ r V T r .
(3.48)

Analysis

At every time t i , i ∈ {1, ..., n}, new observed data y o t i are available. These observations are assumed to be unbiased, which means ǫ o t i = 0. An error covariance matrix R t i is associated with these observations. The observation error covariance matrix R is assumed here to be the same for any time. R and y o t i are standardized with respect to y b and S b y

y o t i = S b y -1 y o t i -y b , (3.49) R = S b y -1 R S b y -1 . (3.50)
Independently, a prior estimate of the state is computed. It is the forecast state x f t i , assumed to be unbiased ( ǫ f t i = 0). Its associated error covariance matrix is P f t i . P f t i and x f t i are standardized with respect to x b and S b

x f t i = S b -1 x f t i -x b , (3.51) P f t i = S b -1 P f t i S b -1 . (3.52)
To further simplify this problem, we assume that

P f t i ≈ P b ≈ V r Λ r V T r . (3.53)
Under this assumption, the reduced state Kalman filter is given by

K r = ( H P f t i ) T H( H P f t i ) T + R -1 . (3.54)
The result of the analysis is finally rescaled using the forecast standard deviation

x a t i = x f t i + S f t i K r y o t i -H x f t i .
(3.55)

Forecast

The forecast state at time t i+1 , x f t i+1 is computed by the convection code STAGYY, taking

x a t i as starting state.

mantle convection and surface tectonics

The associated covariance matrix P f t i+1 is estimated by assuming

P f t i+1 =S f t i+1 V r Λ r V T r S f t i+1
(3.56)

S f t i+1 =αS f t i . (3.57)
The coefficient α is estimated by considering errors on data and using the fact that the error on forecast data and on observed data are uncorrelated:

α = Tr(Covar(Hx f t i+1 -y o t i+1 ) -R) Tr(Covar(Hx f t i -y o t i ) -R)
.

(3.58)

Synthetic Experiments

We test here our sequential assimilation algorithm on synthetic cases. The parameters chosen for the model 1 are described in Table 3.1, along with the parameters of two additional models: model 2, which has a higher bottom heat flux to surface heat flux ratio (37%), and model 3, which has a higher Rayleigh number (10 7 ). The time in the evolutions computed is rescaled as classically done, using the transit time of the mantle as the relevant scaling time (see Gurnis, 1986b, for example). The Earth's mantle transit time t E t estimates the average time spent by a particle to move from the surface to the core-mantle boundary. For the Earth, it is evaluated by considering the thickness of the mantle D E and the root mean square of surface velocities v E rms

t E t = D E v E rms . (3.59) v E
rms is obtained from plate tectonics reconstructions [START_REF] Shephard | The tectonic evolution of the arctic since pangea breakup: Integrating constraints from surface geology and geophysics with mantle structure[END_REF]Seton et al., 2012, for example). t M t , the transit time for the model can be computed by the same procedure. The rescaling of time for a computed evolution is

t = t * t E t t M t , (3.60)
where t * is the non-dimensional time. This rescaling is designed so that the time of the evolution can compare to that of the Earth, even if the vigour of convection differs.

Setup of the Experiments

We test our data assimilation method on synthetic experiments. They consist of two phases.

First, we synthesize a true state evolution x t 1 , x t 2 , ..., x t nt and corresponding observed data sets y o t 1 , y o t 2 , ..., y o tn . This is done by computing a convection evolution spanning 600 My and considering it as the true state evolution x t 1 , x t 2 , ..., x t nt . From this evolution, we extract sets of surface heat flow and surface velocities at regular intervals ∆t. These data are randomly noised to produce y o t 1 , y o t 2 , ..., y o tn . We set the average amplitude of the random noise as a fraction γ of q rms and V rms , the root mean squares of surface heat flux and velocity, respectively. The observed error covariance matrix is then and γ = 10%. The exact combinations of parameters tested are presented in Table 3.2.

R = γ 2             q 2 rms ... 0 0 ... 0 0 ... 0 0 ... 0 0 ... q 2 rms 0 ... 0 0 .. 0 V 2 rms ... 0 0 ... 0 0 ... V 2 rms             . ( 3
Table 3.2 defines as well the name of the different assimilations that we will use in the following.

Quality of the data assimilation estimate

We evaluate the quality of the data assimilation scheme on its ability to retrieve the true temperature fields and to match surface data. ). The first column displays the true temperature field for each experiment, the second column shows the analyzed temperature field, and the third column is the associated error multiplied by 10 and represented using the same color scale than the field itself, in order to ease visual inspection.

ǫ f T (φ m , r n ) = |T f (φ m , r n ) -T t (φ m , r n )|. ( 3

mantle convection and surface tectonics

Below the boundary layer, the locations of positive temperature anomalies are better retrieved for t10γ10 than for t50γ10. The general geometry of slabs at the bottom is representative of the true state, even though some details are missing. The geometry of upwellings is retrieved for t10γ10, but only their approximate location is well estimated for t50γ10. Some minor positive anomalies are missing at the bottom of the domain for both cases. The relative error shows that although the geometry of the temperature field is reasonably well recovered, the actual temperature values for downand upwellings are not as accurately estimated. This tendency is especially noticeable for t50γ10. Overall, the convection pattern is correctly estimated, but for t50γ10 slabs are blurred. The structures are also well recovered for t10γ10Q40 and t10γ10Ra7. We can see that the general shape of the hot upwellings is well estimated for both cases, except on one upwelling, on the bottom right corner of the t10γ10Ra7 experiment, where the hot upwelling with two branches is estimated as only one bigger upwelling. The structure of slabs are also reasonnably well estimated. The value of the errors for both t10γ10Q40 and t10γ10Ra7 are higher than those of t10γ10, on some localized regions.

We compute the spatial average error of the temperature field for each timestep i ∈ {1, 2, ..., n t }

ǫ f T (i) = M m=1 N n=1 (T f i (φ m , r n ) -T t i (φ m , r n )) 2 ∆V (φ m , r n ) M m=1 N n=1 T t i (φ m , r n ) 2 ∆V (φ m , r n ) , (3.63)
where ∆V (φ m , r n ) is the volume of the cell centred on (φ m , r n ). For each combination of ∆t and γ, ǫ f T (i) is averaged over 21 synthetic experiments differing only by their initial conditions for the model 1. The synthetic experiments t10γ10Q40 (higher bottom to surface heat flux ratio) and t10γ10Ra7 (higher Rayleigh number) showed respectively 20 and 19 successes out of 21. We identify a data assimilation experiment as failing if there is a time after which the global error ǫ f T is consistently increasing when analyses are performed. The three failure cases show a consistently increasing difference between observed data and forecast data, making them easy to identify even without the knowledge of the true state evolution. Picking another random noise of the 1.10 -1 9.10 -2 8.10 -2 7.10 -2 6.10 -2 5.10 -2

-2

Error FIGURE 3.7: Evolution of the average relative error on the forecast temperature field with γ = 10% and ∆t = 1, 5, 10, 15, 20 and 50 My. Note that the scale on the y-axis is logarithmic.

same magnitude for each of the three experiments that failed produced a successful data assimilation trial. Hence, peculiar random noise structure could prevent accurate retrieval by the methodology presented here. We compute the average evolution of errors for synthetic experiments t10γ10Q40 and t10γ10Ra7 taking into account only the successes, ie 20 experiments for t10γ10Q40 and 19 for t10γ10Ra7. Hence, the estimation of the error is not specific to a case that is easier to retrieve than others. On every curve, the analysis times are characterized by a sudden drop of the error (see Figs 3.7, 3.8 and 3.9). This is remarkable for the first analysis, at time 0, where the error is reduced by a factor of two for any set of parameters tested. These results show that taking into account only the second order statistics is efficient at correcting a forecast state and improving the estimate of the true state.

As observed in Figs 3.7, 3.8 and 3.9, the general evolution of errors with time follows three phases: an error reduction phase, which may be followed by an error stabilization phase around the lower value, and finally an error growth phase. The duration and existence of each phase depends on the chosen combination of ∆t, γ and the model If analysis intervals are equal to or longer than 15 My, growth of the error happens within 600 My of assimilation. In contrast, for t1γ10, t5γ10 and t10γ10, no error growth occurs, and errors are stable at around 5% after 200 My. Throughout the assimilation, errors for t1γ10 become greater than the errors for t5γ10. This is due to the method used for updating forecast errors in our scheme. Indeed, their estimation is based on the comparison between observed and forecast surface data. A 1 My interval between analyses does not allow the possible errors made internally to be propagated onto the surface. It follows that the error on the forecast is underestimated for short intervals of analysis, leading to results that are not as good as for longer intervals of analysis. The error evolutions for t10γ10Ra7 and t10γ10Q40 are shown in Fig. 3.9, along with the evolution of errors of t10γ10. The error for t10γ10Ra7 is on average smaller than the error for t10γ10. On the contrary, the error for t10γ10Q40 is on average bigger than the error of t10γ10. We note as well that, between two analyses, the errors are increasing much more for t10γ10Ra7 than for the other cases. Fig. 3.10 demonstrates that our data assimilation scheme provides a posteriori a better estimate of surface data than the a priori observation for extreme cases. Indeed, for t10γ50, after 100 My of data assimilation, the a posteriori retrieved surface heat flux and surface velocities are better estimates of the true values than the ones used as the observations.

We conducted here synthetic experiments, where the true state is known, and the result of data assimilation can be compared to this true state evolution. However, in a data assimilation using real observations, the true state of the system is unknown. In this case, it is necessary to use other diagnostics to evaluate the quality of the assimilation.

One diagnostic is based on the study of the innovation, which is the difference between the observed data y o i and the forecast data Hx f i ,

d i = y o i -Hx f i . (3.64)
The simplest evaluation is to consider the evolution of the cumulative mean innovation

d k = k i=1 d i (3.65)
which is supposed to converge to zero if the scheme is unbiased [START_REF] Talagrand | A posteriori validation of assimilation algorithms[END_REF].

Discussion

We chose a sequential method because of its relative ease of implementation compared to alternative methods such as variational data assimilation (see [START_REF] Talagrand | Assimilation of observations, an introduction[END_REF][START_REF] Evensen | Data assimilation[END_REF], for discussions on both methodologies). Moreover, sequential data assimilation has proved to be efficient in solving a wide range of geophysical problems [START_REF] Aubert | Inferring internal properties of earth's core dynamics and their evolution from surface observations and a numerical geodynamo model[END_REF]Hoteit & Pham, 2004, for example). We opted for a suboptimal scheme based on the Kalman filter instead of driving the flow with surface velocities so as to exploit the surface information a step further. The initial condition, or first guess, is given as a compromise between a 1D temperature profile and the inversion of the first data available. Then, the whole temperature field is updated mantle convection and surface tectonics whenever observations are available. This means our models are less sensitive to initial conditions than mantle circulation models in which plate velocities are imposed as boundary conditions.

The specificity of sequential data assimilation is that the estimated state of the mantle improves through time, as more data are assimilated into the model (Fig. 3.1). This raises an issue since the information at hand for the Earth tends to improve in accuracy as it gets closer to present-day: classical sequential data assimilation schemes are not able to extend the recent information back in time. This difficulty can be overcome by implementing Kalman smoothers [START_REF] Cohn | A fixed-lag kalman smoother for retrospective data assimilation[END_REF][START_REF] Cosme | Implementation of a reduced rank square-root smoother for high resolution ocean data assimilation[END_REF][START_REF] Nerger | On the influence of model nonlinearity and localization on ensemble kalman smoothing[END_REF] which use time correlations to update the previous states with new observations.

Applying these techniques would also enable us to add tomographic models as observations in mantle circulation models, and to propagate this information back in time.

However, one of the benefits of the method we developped would be to provide an image of the present state of the Earth mantle which is independant from seismic data, and would give the opportunity to work on the interpretation of tomographic models.

Our sequential assimilation technique gives encouraging results for our test cases: data assimilation gradually reduces the distance between the true models and the estimated ones for the first 300 My, even for highly noised data (up to 30% of noise). The analysis is based on a linear correction, the amplitude of which depends on the distance between the predicted data and the observed data. The linearity of the correction is an approximation that is efficient as soon as the prediction is close enough to the observation. When these two are too different, the analysis tends to over-correct the state, which leads to an unrealistic analyzed state. The algorithm is stable as long as the time between two observations is short enough: if it is 15 My and greater the estimated temperature field will eventually diverge from the true temperature field. If the time between two analyses is too long, the linear approximation for correcting the state field is not valid any more. The temperature field is over-corrected and this leads to analyzed temperatures inconsistent with the physical model.

We evaluate hereafter how the required ∆t and γ in our test cases compare to the corresponding values available for Earth's data: time between two plate reconstruction stages and uncertainty on plate kinematics and seafloor age maps. The seafloor spreading isochrons used for tectonic reconstructions are based on the identification of magnetic anomalies (Müller et al., 2008a). This gives a strong constraint on the motion of plates at least every 10 My, except for the Cretaceous superchron spanning between 83.5 to 125 Ma. Other geological arguments have been provided to infer plate evolutions for this period [START_REF] Torsvik | A new scheme for the opening of the south atlantic ocean and the dissection of an aptian salt basin[END_REF]. For more recent times (the last 20 My), more precise surveys of magnetic anomalies have led to regional reconstructions with a temporal resolution of around 1 My [START_REF] Merkouriev | Constraints on indian plate motion since 20 ma from dense russian magnetic data: Implications for indian plate dynamics[END_REF]. Moreover, Gurnis et al. ( 2012) developed a way to estimate kinematic states between two stages of reconstructions. With this method, having data to assimilate at least every 10 My is ensured, although interpolated solutions could miss peculiar changes in plate motions, and errors on the reconstructed states propagate into the interpolated states.

Regarding these errors, noise in plate tectonic reconstructions is difficult to estimate since the process involves taking into account various types of data and human syntheses. Hence, uncertainties on plate reconstructions increase as we go back in time, since less data are available and more interpretation is required. The accuracy of reconstructions also varies in space. For instance the evolution of the Atlantic for more than 100 My is well known, however, the eastern Pacific before 60 My is mostly unknown because the oceanic seafloor has subducted. Efforts to estimate errors on plate tectonic reconstructions have been directed to the estimation of errors on present maps of oceanic sea floor (Müller et al., 2008a), and on the noise reduction in finite rotations deduced from them [START_REF] Iaffaldano | Redback: Open-source software for efficient noise-reduction in plate kinematic reconstructions[END_REF]. Both of these works show that the error of age maps is smaller than 10% for recent times (< 50 My). However, errors back in time are more difficult to estimate since accuracy for certain plates decreases (like those in the eastern Pacific as discussed before for instance), and an error on the angular velocity or the position of the Euler pole would constrain the errors on velocities over the whole plate domain [START_REF] Molnar | A method for bounding uncertainties in combined plate reconstructions[END_REF]. These errors would have to be propagated to seafloor age maps.

Chapter 3. A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics

As successful as it may appear, the method presented here has several limitations.

First, our technique does not take into account non-linearities of the model physics, and the non-Gaussianity of the distribution of temperature, heat flux and velocity values. Indeed, we computed only second order statistics to estimate the link between temperatures and velocities. However, the method we present here already captures the essential features of mantle convection, and provides accurate convection reconstructions with the approximations made.

In this study, we assumed a perfect model (with η = 0 as stated in Equation 3.12).

However, this would not be the case for applying the technique with Earth data. Model equations assume simplifications of the physics, detailed in [START_REF] Ricard | Physics of mantle convection[END_REF]. Any parametrization we may use would have shortcomings. The most important, least known, ingredient is thought to be rheology, which is a strong limitation according to [START_REF] Worthen | Towards adjointbased inversion for rheological parameters in nonlinear viscous mantle flow[END_REF]. Applying our method to Earth data would help evaluate the forecasting power of the state-of-the-art convection models. Explicitly introducing a model error term in the data assimilation scheme could be a way to improve the reconstructions [START_REF] Evensen | Data assimilation[END_REF]. Also, we used heat flux as a proxy for seafloor age because there is no small scale convection in our test cases. For Earth data, monitoring age with tracers, giving the exact estimate of seafloor age would be required. The state vector and observation matrix would differ slightly, but the technique would still apply. Complexifying the direct model used in data assimilation (ie having a 3D model at higher Rayleigh number, with a more complex rheology and a multi component system) will as well increase dramatically the size of the data assimilation problem. In particular, the state covariance matrix will reach a size at which the direct computations carried out in this paper will be extremely heavy. One solution would be to take advantage of the structure and properties of the covariance matrix (symmetries, periodicity) to design faster algorithms. The method could also be modified to avoid the computation of the covariance matrix and calculate directly the smaller matrix H P f t i .

Conclusion

We have applied sequential data assimilation methodology to reconstruct mantle convection using surface data (velocity and heat flux as a proxy for seafloor age). This work is a first approach to test the efficiency of assimilating surface data to recover the evolution of convecting structures. Our technique belongs to the suboptimal schemes of sequential data assimilation based on the Kalman filter. We modified the optimal interpolation method to update the forecast error covariance matrix.

Our scheme of sequential data assimilation proves to be efficient in recovering the temperature field of a convective system with plate-like tectonics at its surface over several 100 My. The only observations used were surface heat fluxes and surface velocities. We tested the robustness of the method by conducting synthetic experiments in 2D spherical annulus geometry. We obtained accurate results, even for the location of plumes at the base of the models, for periods of at least 300 My provided the time between analyses is shorter than 15 My and the noise in observation is lower than 30%. These requirements are satisfactory, since current plate reconstruction models already provide estimates for velocities and ages at least every 10 My and with a fine accuracy for the past 60 My, and possibly beyond.

Application to the Earth would first involve the use of 3D-spherical sophisticated models, with high resolution, high convective vigour, continents and a more realistic rheology. In 3D the ratio of number of data points to unknowns would be the same as in 2D spherical annulus geometry by definition. However, increasing the resolution of our model would make this ratio deteriorate, which could be compensated by increasing the correlation between closer nodes. Moreover, although convection models can account for complex parametrizations, their limitations could introduce errors in the reconstructions.

These difficulties should not stop the exploitation of the method: if successful, the convection reconstruction would provide a new image of the mantle, alternative to seismic tomography, which would appear as an independent source of information and assess the quality of the estimate of the state of the mantle at present day; if unsuccessful, it would provide a quantitative evaluation of forward models and help us decide how to improve them.

acknowledgments

We like to thank the editor as well as the two anonymous reviewers for their comments that helped to improve the initial manuscript.

Résumé du chapitre

Grâce aux avancées réalisées au cours de la dernière décennie dans le domaine de la modélisation de convection mantellique, il est maintenant possible de modéliser au

CHAPTER 4 Ensemble Data Assimilation For Mantle Circulation

To be submitted in: [START_REF] Bunge | Time scales and heterogeneous structure in geodynamic earth models[END_REF][START_REF] Mcnamara | Thermochemical structures beneath africa and the pacific ocean[END_REF][START_REF] Bower | Lower mantle structure from paleogeographically constrained dynamic earth models[END_REF][START_REF] Davies | Reconciling dynamic and seismic models of earth's lower mantle: The dominant role of thermal heterogeneity[END_REF], to study the evolution of mantle plumes and their relationship to hotspots [START_REF] Hassan | A rapid burst in hotspot motion through the interaction of tectonics and deep mantle flow[END_REF], to infer changes in the Earth's rotation axis [START_REF] Steinberger | Changes of the earth's rotation axis owing to advection of mantle density heterogeneities[END_REF], sea-level [START_REF] Moucha | Dynamic topography and long-term sea-level variations: There is no such thing as a stable continental platform[END_REF] or dynamic topography [START_REF] Flament | A review of observations and models of dynamic topography[END_REF].

The geodynamics community has developed three alternative approaches to the problem of the reconstruction of mantle circulation. The first approach, backward advection, consists in starting at present by estimating the current density field of the mantle from seismic tomography models (see [START_REF] Conrad | Seismic tomography, surface uplift, and the breakup of gondwanaland: Integrating mantle convection backwards in time[END_REF], for a description of this method). This density field is then advected backward in time with plate tectonic reconstructions as imposed boundary condition [START_REF] Steinberger | Changes of the earth's rotation axis owing to advection of mantle density heterogeneities[END_REF]. This method has a limited numerical cost and exploits the two most instructive constraints on mantle circulation: plate tectonic reconstructions and seismic tomography. However, this technique neglects thermal diffusion, so it is not able to reconstruct past thermal structures that have completely diffused before present and it is limited to times and regions for which the effect of diffusion is thought to be small. This limits reconstructions to the last 50 to 75 Myr [START_REF] Conrad | Seismic tomography, surface uplift, and the breakup of gondwanaland: Integrating mantle convection backwards in time[END_REF] or even to shorter periods if we consider the uncertainties on tomographic models [START_REF] Bello | On the predictability limit of convection models of the earth's mantle[END_REF]. The second approach, the semi-empirical sequential method, estimates mantle circulation by integrating plate tectonics reconstructions chronologically into a mantle convection model. Plate tectonic reconstructions are either introduced as velocity boundary conditions, as first described by [START_REF] Bunge | Time scales and heterogeneous structure in geodynamic earth models[END_REF], or with a more sophisticated method, by blending a convection solution with thermal and kinematic models of plates and slabs [START_REF] Bower | Assimilating lithosphere and slab history in 4-d earth models[END_REF]. This approach allows the use of models of convection with chemical heterogeneities [START_REF] Mcnamara | Thermochemical structures beneath africa and the pacific ocean[END_REF]. Also, it is not anymore the reconstruction method that limits the timespan of the reconstruction, but the availability of plate tectonic reconstructions. This led to mantle circulation models integrating up to 450 Myr of plate reconstruction history [START_REF] Zhang | A model for the evolution of the earth's mantle structure since the early paleozoic[END_REF]. However, this method considers plate tectonic reconstructions as perfect estimates of surface tectonics: uncertainties affecting the reconstructions are not taken into account although they are substantial, especially as reconstructions go further in the past (for example, there is almost no information on the state of the ocean floor before 140 Myr, see e.g. [START_REF] Torsvik | Plate tectonics and net lithosphere rotation over the past 150my[END_REF]. This method also requires the choice of an arbitrary initial temperature field to compute the evolution. The third approach uses data assimilation methods to solve the mantle circulation problem. Data assimilation methods are inverse methods dealing with the specific problem of estimating the evolution of a dynamical system from asynchronous data and a physical model (Evensen, 2009a). The full inverse problem for mantle circulation, as stated by [START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF], would take into account model errors, numerical approximations, errors on plate reconstructions and on the estimation of the current tomography-derived temperature field to provide the best fit given all sources of information. However, solving the full inverse problem of mantle circulation is still a great challenge given the nonlinearities in mantle convection dynamics and the computational power required to compute a realistic forward mantle convection evolution alone [START_REF] Stadler | The dynamics of plate tectonics and mantle flow: From local to global scales[END_REF][START_REF] Burstedde | Large-scale adaptive mantle convection simulation[END_REF]. So far, variational data assimilation dominates over other methods to estimate mantle circulation [START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF][START_REF] Horbach | The adjoint method in geodynamics: derivation from a general operator formulation and application to the initial condition problem in a high resolution mantle circulation model[END_REF]Ghelichkhan & Bunge, 2016). To simplify the problem, they minimize the misfit between the final temperature field of the mantle circulation model and the one deduced from seismic tomography. These mantle circulation models impose plate tectonic reconstructions as boundary conditions, as in the first two approaches.

Here, we take a different view on data assimilation methods for mantle circulation models by focusing on how to take into account the uncertainties in plate tectonic reconstructions. For almost a decade, 3-D spherical mantle convection models have shown the capability to self-consistently produce plate-like tectonics at their surface [START_REF] Walzer | Mantle convection and evolution with growing continents[END_REF][START_REF] Van Heck | Planforms of self-consistently generated plates in 3d spherical geometry[END_REF][START_REF] Yoshida | Mantle convection with longest-wavelength thermal heterogeneity in a 3-d spherical model: Degree one or two?[END_REF][START_REF] Foley | Generation of plate-like behavior and mantle heterogeneity from a spherical, viscoplastic convection model[END_REF]. These models physically link surface tectonics comparable to that of the Earth to mantle convection processes [START_REF] Coltice | Dynamic causes of the relation between area and age of the ocean floor[END_REF][START_REF] Rolf | Statistical cyclicity of the supercontinent cycle[END_REF][START_REF] Mallard | Subduction controls the distribution and fragmentation of earth's tectonic plates[END_REF]. In [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF], we took advantage of this link to build a sequential data assimilation algorithm able to integrate plate reconstructions into a mantle convection code while taking into account the uncertainties on those plate tectonic reconstructions.

This technique assimilates a time series of surface observations chronologically, by repeating two stages, analysis and forecast, until all observations are taken into account.

Whenever an observation is available, the analysis evaluates the most likely state of the mantle at this time, considering a prior guess (supplied by the forecast) and the new observations at hand. For this evaluation, we used the classical best linear unbiased estimate [START_REF] Talagrand | Assimilation of observations, an introduction[END_REF]. Then, the forward model of mantle convection forecasts the evolution of the mantle until the next observation time. We tested this algorithm on synthetic experiments. It proved to be efficient in recovering mantle circulation given constraints on the amplitude of errors affecting observations and the timespan between observations. Here we extend this work by applying a more advanced sequential data assimilation method, the Ensemble Kalman Filter (EnKF, described in [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics[END_REF][START_REF] Burgers | Analysis scheme in the ensemble kalman filter[END_REF]. This method is particularly suited for high dimensional nonlinear dynamical models (Evensen, 2009b). Instead of estimating the most likely state of the mantle, the Ensemble Kalman filter provides at each time an approximation of the probability density function of the state of the system in the form of a finite ensemble of states. During the forecast stage, each member of the ensemble evolves independently. For the analysis, we use the second order statistics of the ensemble to correct each ensemble member with the new observations at hand. We evaluate this method with synthetic experiments in 2D-spherical annulus geometry [START_REF] Hernlund | Modeling mantle convection in the spherical annulus[END_REF] and compare it to the algorithm developed by [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF]. The EnKF provides more accurate estimations than the former method, and is even able to reconstruct evolutions that the former method could not. Moreover, the EnKF also estimates locally the error on the reconstruction. The optimal size of the ensemble for our test case is 300 members. Both covariance inflation and localization eliminate spurious correlations arising from the finite size of the Ensemble that is used to compute them.

This paper is organized as follows. In section 2, we present our simplifications on the general mantle circulation reconstruction problem and the correspondence with the notation in the EnKF algorithm. Then, in section 3, we detail the EnKF method and justify the variants chosen for the application to mantle circulation. Section 4 presents the results obtained on synthetic experiments and compares them to results obtained by the method described in [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF]. Section 5 is a discussion on the choice of the method and the challenges involved in the application of such a method to a realistic setting.

Presentation of the Problem

We aim at reconstructing mantle circulation for the last hundreds of millions of years by combining a mantle convection model with plate tectonic reconstructions, using an Ensemble Kalman Filter. To study the behavior of the Ensemble Kalman filter on such problem, we consider a simplified mantle convection model. This section describes the model used to compute a mantle evolution, the data set assimilated in this evolution, and finally the backbone of Ensemble Kalman Filtering.

Mantle Convection Model

At the timescales and lengthscales we are interested in (≥ 10 kyr, ≥ 1000 km), the mantle can be modelled as a continuous visco-plastic medium. To compute mantle circulation, we solve the equations of conservation of mass (Eq. 4.1 below), momentum 4.2. Presentation of the Problem 123 (Eq. 4.2 below) and energy (Eq. 4.8 below) for an isochemical mantle under the Boussinesq approximation. The system of equations is non-dimensionalized to the thermal diffusion time scale (see [START_REF] Ricard | Physics of mantle convection[END_REF]. Given the high Prandtl number of the mantle (of the order of 10 24 ), inertia is neglected. With these assumptions, the equations of conservation of mass and momentum become diagnostic equations of the form

∇ • u = 0, (4.1)
∇ • σ -∇p + RaT e r = 0, (4.2)
where σ, u, p, and T are the non-dimensional deviatoric stress, velocity, dynamic pressure, and temperature, respectively. The equations are written in spherical coordinates (r, θ, φ), using the physical convention with r the radius, θ the colatitude and φ the longitude. The associated unit vectors are (e r , e θ , e φ ).

Ra is the Rayleigh number, defined as

Ra = ρ 0 g 0 α 0 ∆T a 3 µ 0 κ 0 , (4.3)
with ρ 0 the density for T = 0, g 0 the gravitational acceleration, α 0 the thermal expansivity, ∆T the temperature drop, a the depth of the layer, κ 0 the thermal diffusivity, µ 0 the dynamic viscosity of the system. The Rayleigh number in our model is 10 6 . It is one or two orders of magnitude lower than that of the Earth, but high enough to ensure chaotic convection. The vertical velocities and shear-stress at the surface and the base of the model are set to zero.

The deformation response of mantle material to stress is implemented as a linear relationship linking the strain rate tensor ǫ to the deviatoric stress tensor σ as

σ = 2µ eff ǫ = µ eff ∇u + (∇u) T . (4.4)
The effective viscosity µ eff takes into account both a viscous Newtonian behavior with a viscosity µ n and a pseudo-plastic behavior with an equivalent "pseudo-plastic viscosity" µ y , µ eff = min(µ n , µ y ).

(4.5)

The Newtonian viscosity µ n follows an Arrhenius law

µ n = µ 0 exp E A T + T 1 (4.6)
with µ 0 = exp -E A 2T 1 , T 1 the temperature at which µ n = 1, and E A the nondimensional activation energy. We implement the decrease of viscosity in the asthenosphere by reducing by a factor of 10 the viscosity µ n when the temperature is above a solidus equation T s = T s 0 + ∇ r T s (r a -r) with r a the surface value of r. The implementation of a weak asthenosphere tends to favor plate-like behavior (Tackley, 2000b;[START_REF] Richards | Role of a lowviscosity zone in stabilizing plate tectonics: Implications for comparative terrestrial planetology[END_REF].

The pseudo-plastic part of the effective viscosity µ y is defined by

µ y = σ yield 2 ǫII , (4.7)
where ǫII is the second invariant of the strain rate tensor and

σ yield = σ Y + (r a -r)∇ r σ Y ,
with σ Y and ∇ r σ Y the yield stress at the surface and the depth-dependence of the yield stress, respectively.

The energy conservation equation is the only prognostic equation of the system

DT Dt = ∇ 2 T + R h . (4.8)
R h is the non-dimensional internal heating rate defined as

R h = ρ 0 D 2 H k 0 ∆T (4.9)
with H the dimensional heating rate and k 0 the thermal conductivity. We set isothermal top and bottom boundaries with temperatures T a and T b , respectively. The models presented here have 10% basal heating and 90% internal heating.

These equations are solved using the finite volume, multigrid parallel code STAGYY [START_REF] Tackley | Effects of an endothermic phase transition at 670 km depth in a spherical model of convection in the earth's mantle[END_REF], on a spherical annulus grid. This geometry provides results closer to the spherical grid than cylindrical geometry [START_REF] Hernlund | Modeling mantle convection in the spherical annulus[END_REF]. In the following, the longitudinal coordinate of a point is φ l , with l ∈ {1, 2, ..., L} and its radial coordinate is r m with m ∈ {1, 2, ..., M }, r varying from r b to r a .

Note that this paper focuses on the methodology of ensemble data assimilation for a convecting system similar to that of the Earth's mantle. Hence, we choose a rather simple model that can reproduce plate-like tectonics at the surface. We rely on simplifications such as 2D geometry, incompressible and isochemical mantle and a rheology which does not take into account the history of the material. Although some of the complexities we ignore may play a fundamental role in the reconstruction of the Earth's mantle evolution, we choose to focus in this manuscript on the data assimilation methodology. Moreover, we choose to keep the same parameters as the test case of [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF] in order to allow direct comparison between the methods. Table 4.1 lists the chosen parameter values.

To ease the comparison with Earth's mantle convection, we rescale the nondimensional time in the evolution, t, by the transit time of the convective system. By definition, the transit time of the Earth's mantle is t E t = a E /v E rms , with a E the thickness of the mantle and v E rms the root mean square of surface velocities of the Earth, as estimated by plate tectonic reconstructions [START_REF] Seton | Global continental and ocean basin reconstructions since 200ma[END_REF]. We compute the same value for the model

t m t = a/v m rms . The scaled time t s is then t s = t t E t t m t .

Observations of Mantle Circulation

The state of the Earth's surface is the time integrated expression of mantle circulation.

At a global scale, the main source of information for the last 100 Myr is the database of the localization and identification of magnetic anomalies on the seafloor, translated into maps of seafloor ages (Müller et al., 2008a;[START_REF] Seton | Community infrastructure and repository for marine magnetic identifications[END_REF]. This information is complemented with regional geological studies giving constraints on the timing and geometry of tectonic events as well as a synthesis of paleontological, structural geology, stratigraphical, magnetic anomalies, gravity data and seismic studies. In addition, paleomagnetic data provide constraints on the paleolatitude of continental blocks [START_REF] Besse | Apparent and true polar wander and the geometry of the geomagnetic field over the last 200 myr[END_REF].

Plate tectonic reconstructions use the geometric theory of plate tectonics to integrate all these observations. The result is a time series of maps of seafloor ages, plate layout and kinematics. The continuously closed plate algorithm [START_REF] Gurnis | Plate tectonic reconstructions with continuously closing plates[END_REF] produces plate tectonic reconstruction maps continuous in space and time [START_REF] Seton | Global continental and ocean basin reconstructions since 200ma[END_REF][START_REF] Müller | Ocean basin evolution and global-scale plate reorganization events since pangea breakup[END_REF].

Although we are aware that these plate tectonic reconstruction maps are in themselves models and not direct observations, we propose to develop an assimilation method that use them as data to assimilate in our mantle convection model. This solution is generally chosen in mantle circulation reconstructions [START_REF] Bunge | Mantle-circulation models with sequential data assimilation: inferring present-day mantle structure from platemotion histories[END_REF][START_REF] Zhang | A model for the evolution of the earth's mantle structure since the early paleozoic[END_REF][START_REF] Bower | Assimilating lithosphere and slab history in 4-d earth models[END_REF], because it provides continuous surface boundary conditions in space and time for the period of reconstruction. One advantage of the technique we develop is that it is possible to consider errors on the data that is assimilated, another is that the reconstructions do not need to be known at all times and at all points on the surface. Hence it is possible, in principle, to design a data assimilation scheme using direct observations. However, this would require further developments both on the database design and on the data assimilation algorithm. Sequential data assimilation methods for mantle circulation are still in their infancy, so we opt for a simpler structure of the data to be assimilated: a time series of maps of surface velocity and seafloor age, as given by plate tectonic reconstructions.

In this study, we limit ourselves to the test of data assimilation in synthetic experiments. In the model described in subsection 4.2.1, the absence of small scale convection at the base of the boundary layer makes the surface heat flux an excellent proxy for the age of the seafloor [START_REF] Coltice | Dynamic causes of the relation between area and age of the ocean floor[END_REF]. Consequently, we consider surface heat flux and surface velocity as the data to assimilate.

To our knowledge, the amplitude of the uncertainty on global plate tectonic reconstructions has not yet been assessed. For the synthetic tests we perform in section 4.4, we choose an arbitrary value of 10% of the root mean square value of heat flux and surface velocity, respectively. We further discuss this choice in section 4.5.

Ensemble Kalman Filtering Framework: Notations

Our aim is to assimilate a time series of observations (surface velocities and heat fluxes)

into a mantle convection model to estimate the evolution of the state of the mantle. We introduce here the general formulation of Ensemble Kalman Filtering and link them to our problem. We use the notation system recommended by [START_REF] Ide | Unified notation for data assimilation: operational, sequential and variational[END_REF].

The time series of data is defined as a set of column vectors {y o 1 , y o 2 , ..., y o K }, where the subscripts {1, 2, ..., K} refer to the times at which observations are available. As seen in the previous section, the data used for our experiments are surface velocity and surface heat flux. The data vector at time k is thus defined as

y o k = q o k (φ 1 ), q o k (φ 2 ), ..., q o k (φ L ), u o φk (φ 1 ), u o φk (φ 2 ), ..., u o φk (φ L ) T , (4.10)
where q o k (φ l ) and u o φk (φ l ) are the observed values of surface heat flux and surface horizontal velocity at the k-th timestep and longitude φ l , and (•) T means transpose. We model errors on observations by a random vector of zero mean and covariance matrix R k (we suppose unbiased observations). Although R k is a diagonal matrix of constant value and size in our experiments, it is not generally the case. Correlations between errors on observations could be specified in R k .

The evolution of the state of the system is estimated sequentially during the period where observations are available. At each timestep k ∈ {1, 2, ..., K}, we define two state vectors: the a priori state, or forecast state x f k and the analysis state x a k , which is the state corrected after having assimilated the observations y o k . The system of equations developed in section 4.2.1 shows that we can compute velocity, viscosity and pressure values at each grid point from the temperature field. Nevertheless, the relation between surface velocities and the temperature field is nonlinear. We choose to include the whole temperature field and the surface velocities to form an augmented state vector. This simplifies the computations thereafter. The state of the mantle at a timestep k ∈ [1, K] is defined as

x k = [T k (φ 1 , r 1 ), T k (φ 1 , r 2 ), ..., T k (φ L , r M ), u φk (φ 1 ), u φk (φ 2 ), ..., u φk (φ L )] T , (4.11)
where T k (φ l , r m ) and u φk (φ l ) are the values of temperature at the kth timestep, longitude φ l and radius r m and surface horizontal velocity at the kth timestep and longitude φ l .

The forecast and analyzed states are uncertain as well. Their uncertainties are represented by two random vectors of zero expectancy and covariance matrices P f k and P a k , respectively. We do not compute explicitely these covariance matrices. Instead, we Finally, we introduce the observation operator, which maps a given state vector x e kn (e being f or a) to the corresponding data y e kn . If the surface heat flux is approximated by a first order discretization of Fourier's law, then the observation operator is linear, and can be represented by the matrix H such that ∀k ∈ {1, 2, ..., K}, ∀n ∈ {1, 2, ..., N }, y e kn = Hx e kn .

(4.12) Table 4.2 summarizes the dimensions of the vectors and matrices for our problem.

Ensemble Kalman Filter with Localization and Inflation

The Ensemble Kalman Filter [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics[END_REF][START_REF] Burgers | Analysis scheme in the ensemble kalman filter[END_REF] is a sequential data assimilation algorithm using the same equations as the Kalman Filter for the analysis step, but Monte Carlo methods to forecast the error statistics on the state. We explain here how we adapt the Ensemble Kalman Filter to our problem and justify the choice of the starting ensemble.

To implement the EnKF, we used the software environment Parallel Data Assimilation Framework (PDAF, [START_REF] Nerger | The parallel data assimilation framework: experiences with kalman filtering, in Use of High Performance in Meteorology[END_REF][START_REF] Nerger | Software for ensemble-based data assimilation systems-implementation strategies and scalability[END_REF].

Initialization: First Analysis and Generation of the starting ensemble

As in [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF], we compute second order statistics from a series of 400 decorrelated snapshots of convection simulations. We obtain the first forecast state of average x f 1 and associated covariance matrix P f 1 . The background covariance matrix P f 1 is eigendecomposed and rank reduced into P f 1r = VΛV T , with Λ a diagonal matrix containing the n r = 1928 largest eigenvalues of P f 1 (which accounts for 99.98% of its cumulative variance) and V a matrix of the corresponding eigenvectors.

The first set of observations y o

1 is assimilated to obtain

x a 1 = x f 1 + VAV T H T R -1 (y o 1 -Hx f 1 ), (4.13)
P a 1 = VAV T , (4.14)
with

A = Λ -1 + V T H T R -1 HV -1 . (4.15)
We generate an ensemble of N initial states using the second order exact sampling method [START_REF] Hoteit | Filtres de Kalman réduits et efficaces pour l'assimilation de données en océanographie[END_REF][START_REF] Pham | Stochastic methods for sequential data assimilation in strongly nonlinear systems[END_REF]. First, A is eigendecomposed

A = V a Λ a V aT . (4.16)
The ensemble members are then computed following

X a 1 =       | | x a 11 ... x a 1N | |       =       | | x a 1 ... x a 1 | |       + √ N -1VV a Λ a1/2    Ω T N ×(N -1) 0 (nr-N )×N    , (4.17)
where Ω is a random matrix whose columns are vectors forming an orthonormal basis and each of them is orthogonal to 1 = [1, ..., 1] T . Ω is generated through the algorithm described in the appendix of [START_REF] Nerger | A unification of ensemble square root kalman filters[END_REF]. The matrix Ω is designed so that the sample mean of the starting ensemble is equal to x a 1 and its sample covariance matrix is equal to matrix P a 1 reduced to its N largest eigenvalues.

This method of generating the starting ensemble takes advantage of the extensive knowledge we have on the background statistics of the model. Several other methods have been tested to generate a starting ensemble, such as starting with random decorrelated snapshots of mantle convection obtained from a very long run, second order exact sampling from x f 1 and P f 1 , and several assimilations of the first observations y o 1 . None of these solutions were as efficient for our problem as the technique used here.

Forecast

Between timesteps k -1 and k, the forward numerical code STAGYY computes independently the evolution of each of the analyzed states

{x a k-1,n } n∈[1,N ] to produce a forecast ensemble {x f k,n } n∈[1,N ] .
The forecast state is the average of the ensemble

x f k = 1 N X f k 1, (4.18)
and the forecast error covariance matrix is given by the sample covariance matrix of the ensemble of forecast states

P f k = 1 N -1 X f k I - 1 N 11 T I - 1 N 11 T T X f T k . (4.19)
After several assimilation cycles, the finite size of the ensemble induces the underestimation of the error variance [START_REF] Van Leeuwen | Comment on "data assimilation using an ensemble kalman filter technique[END_REF], and can lead to filter divergence.

We observed this behavior in our case, and to stabilize the filter we apply covariance inflation, as suggested in [START_REF] Anderson | A monte carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts[END_REF] and [START_REF] Hamill | Distance-dependent filtering of background error covariance estimates in an ensemble kalman filter[END_REF].

We correct the forecast ensemble variance with an inflation factor γ according to

X f k ← 1 N X f k 11 T + X f k I - 1 N 11 T √ γ, (4.20)
where ← means that we replace the matrix on the left-hand side by the term on the right-hand side. γ is computed following the same principles as in the suboptimal Kalman Filter developed in [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF], i.e. by comparing the error on observations and the standard deviation of the innovation d k defined as

d k = y o k - 1 N HX f k 1. (4.21)
The inflation factor is

γ = V d -V o V f , (4.22)
with

V d = Tr d k d T k , (4.23) V o = Tr(R k ), (4.24)
V f = Tr HX f k I - 1 N 11 T I - 1 N 11 T T X f T k H T , (4.25)
where Tr(•) means the trace. The inflation factor is then truncated between a minimum value of 1 (to prevent further contraction of the ensemble spread) and a maximum value of γ + = 1.25 (to prevent overspread). Several values of maximum inflation factor have been tested, from γ + = 1.1 to γ + = 2, and showed little impact on the efficiency of the assimilation. A constant inflation factor was also tested, but the results with an adaptive inflation factor were substantially more accurate, especially for the first assimilation times.

Analysis

The analyzed state x a kn of the nth member of the ensemble is

x a kn = x f kn + K k y o kn -Hx f kn (4.26)
where K k is the Kalman Gain. y o kn is the observed data vector y o k to which a random perturbation of zero expectation and covariance matrix R k is added, as is recommended in [START_REF] Burgers | Analysis scheme in the ensemble kalman filter[END_REF].

The Kalman Gain is defined as

K k = (P f k • C)H T H(P f k • C)H T + R k -1 , (4.27)
where the matrix P f k is the sample covariance matrix of the ensemble of forecast states

{x f kn } n∈[1,N ]
. We use a limited ensemble size (maximum 768) to estimate P f k . Spurious correlations ensue, especially between distant points. To counteract this effect, we implement direct forecast error localization by Schur multiplying (symbol •) P f k by the localization matrix C, as introduced by [START_REF] Hamill | Distance-dependent filtering of background error covariance estimates in an ensemble kalman filter[END_REF] and [START_REF] Houtekamer | A sequential ensemble kalman filter for atmospheric data assimilation[END_REF]. The matrix C is itself the Schur product of a vertical localization matrix C v and a horizontal localization matrix C h . The value of C v (i, j) depends on the absolute radius difference of the i-th and the j -th components of the state vector and on vertical correlation length ℓ v . The value of C h (i, j) depends on the absolute angle difference of the i -th and the j -th components of the state vector and on vertical correlation length ℓ h . Both values follow a Gaspari-Cohn compactly supported fifth-order piecewise rational function (similar to a Gaussian but with a compact support, Eq. 4.10 of [START_REF] Gaspari | Construction of correlation functions in two and three dimensions[END_REF].

We also tested the domain localization strategy as described in [START_REF] Janjic | On domain localization in ensemble-based kalman filter algorithms[END_REF], since it is computationally more efficient and already implemented in PDAF. However, it led to a systematic failure of the assimilation. This is due to the nature of our problem: all the observations are located at the surface of the model and we aim at estimating the temperature field over the whole depth of the mantle. A vertical localization is as necessary as a horizontal localization: the localization has to be done in the state space and not only in the data space.

Implementation of the Ensemble Kalman Filter

We used the software environment PDAF [START_REF] Nerger | The parallel data assimilation framework: experiences with kalman filtering, in Use of High Performance in Meteorology[END_REF][START_REF] Nerger | Software for ensemble-based data assimilation systems-implementation strategies and scalability[END_REF] in combination with the mantle convection code STAGYY [START_REF] Tackley | Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid[END_REF] to develop an Ensemble Kalman Filter code for mantle convection . PDAF provides a set of core routines computing in parallel the analysis steps for a range of ensemble based data assimilation techniques. It provides as well a set of standard routines to adapt the parallelization of a preexisting parallel forward numerical model and integrate the data assimilation routines. The final product is a highly scalable ensemble data assimilation code running both forecasts and analyses in parallel.

We modified the STAGYY code following the procedure recommended by PDAF (see the online documentation wiki at [START_REF] Nerger | PDAF: the Parallel Data Assimilation Framework[END_REF]. We also made a few modifications in PDAF routines to allow for direct forecast error localization with the Ensemble Kalman filter. Additionally, we designed a basic observation database so as to load in a single step all the observations used in the data assimilation procedure.

A posteriori Evaluation of the Ensemble Kalman Filter Method

We test the data assimilation scheme on twin experiments using the model described in section 4.2.1. Throughout this section, we compare the results of the Ensemble Kalman Filter for mantle circulation reconstructions to the results computed using the method developed in [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF], hereafter referred to as method 1.

After describing the setup used for twin experiments, we test the robustness of the EnKF method and compare it to that of method 1. Then, we determine the range of data assimilation parameters which are suitable to conduct an ensemble data assimilation.

Finally, we assess the ability of the scheme to actually reconstruct specific geodynamic structures.

Twin Experiment Setup

Twin experiments are a way to assess the accuracy of a data assimilation procedure in a controlled environment, where the true evolution is perfectly known.

First, we compute a reference evolution using the forward numerical model, considered as the true evolution, from which we extract the set of true states {x t k } k∈ [1,K] . Here, the timespan of the evolution is 150 Myr and we sample true states every 10 Myr. From these states, we compute a time series of surface heat fluxes and surface velocities, following Equation 4.12. We noise these observations with a random Gaussian noise of standard deviation 10% of the root mean square of surface heat flux q rms and surface velocities v rms , to obtain the time series of observations to assimilate {y o k } k∈ [1,K] . It follows that the observation error covariance matrix R is diagonal and does not change with time.

Then, we perform ensemble data assimilation for the data set {y o k } k∈[1,K] , with the observation error covariance matrix R. We did not consider any model error in the filter we describe, so the parameters of the model used in the data assimilation realizations are the same as those of the reference model.

We present here tests with different assimilation parameters, varying the number of members N , the vertical correlation length ℓ v and the horizontal correlation angle ℓ h .

Table 4.3 details the range of parameters tested.

We compute four different evolutions to test the accuracy of the ensemble Kalman filter for different cases. Figure 4.3 shows the initial and final states of these evolutions, together with the result of global error evolution, and will be discussed in the next section.

Robustness of the Assimilation Algorithm

The evolution of the global error on the estimated temperature field over the time period {1, ..., K} is

ǫ f T (1), ǫ a T (1), ǫ f T (2), ..., ǫ f T (K), ǫ a T (K) (4.28)
where ǫ e T (k), e standing for a (analysis) or f (forecast) is with V(φ l , r m ) the volume of the grid cell at longitude φ l and radius r m , and T e k (φ l , r m ) the average temperature of the estimated ensemble (either forecast or analysis) at longitude φ l and radius r m .

ǫ e T (k) = L l=1 M m=1 T e k (φ l , r m ) -T t k (φ l , r m ) 2 V(φ l , r m ) L l=1 M m=1 T t k (φ l , r m ) 2 V(φ l ,
We test the EnKF on one evolution, with sizes of the ensemble N = 96, 288 and 768 and for each combination of the following data assimilation parameters values: vertical correlation length ℓ v = 0.3, 0.5, 0.7 and 1 and horizontal correlation angle ℓ h = π/10, π/8, π/6, π/4 and π/2 . We show in Figure 4.1, for each ensemble size, the maximum and minimum values of errors obtained for all these parameters, as a function of time.

To determine the best assimilation, we compute the average error after analysis

ǫ a T = 1 K K k=1 ǫ a T (k). (4.30)
For each ensemble size, the error evolution of the best assimilation (in the sense of minimum ǫ a T ) is also shown in Figure 4.1.

For any size of the ensemble, it is possible to find a set of parameters leading to a stabilization of the global error on the temperature field after a few analyses. The time after which the solution is stabilized corresponds to the transit time of the physical model (70 Myr). We can see that for N = 288 and N = 768, any combination of vertical and horizontal correlation lengths leads to errors lower than the first analysis.

Although the error is decreasing through time for any combination of data assimilation parameters, the difference between the maximum and the minimum errors obtained is greater than 1%, which is large given that the first analysis error is already below 8%. The best error evolutions for N = 288 and N = 768 are very similar, with a minimum error of 4.07% and 3.87% after 90 Myr, and an average global error after analysis of 5.01% and 4.85%, respectively. During the assimilation of a dataset, most of the computational time is dedicated to the forecast step, so the data assimilation with 768 members is 2.7 times longer than the assimilation with 288 members, on the account of the embarrassingly parallel nature of the forecast phase. Since we obtain very similar results for N = 288 and N = 768, we favor the assimilation with 288 members.

We compute the error on the estimated temperature from the true temperature field.

However, in a realistic case, the true temperature is not known, and the evaluation of the data assimilation algorithm is based on the study of the statistics of the innovation vector d k at analysis k We also tested the assimilation algorithm for 4 different evolutions, with the optimal parameters for an ensemble size of N = 288 members (ℓ v = 0.7 and ℓ h = π/10). The spread of the ensemble is an estimation of the error on the state. We compare the evolution of ǫ e T to the global standard deviation of the temperature field of the ensemble:

d k = y o k -Hx f k . ( 4
σ f T (1), σ a T (1), σ f T (2), ..., σ f T (K), σ a T (K) (4.33)
with σ e T (k) defined as

σ e T (k) = N n=1 L l=1 M m=1 T e kn (φ l , r m ) -T e k (φ l , r m ) 2 V(φ l , r m ) (N -1) L l=1 M m=1 T e k (φ l , r m ) 2 V(φ l , r m ) . (4.34)
We compute the error for an estimated evolution with the method 1 using Equation 4.29.

Although we ran the four evolutions with the same values of physical parameters (as described in Table 4.1), they show different configurations: Evolution A has a shorter wavelength of convection, with the persistence of 4 subductions, 3 ridges and 5 upwellings, the death of one ridge and creation of two. Evolutions B, C and D have longer wavelengths of convection, with two major downwellings, stable throughout the evolutions. In evolution B, one of these downwellings has a very large negative temperature anomaly at the bottom of the domain. In evolution C, the remnant of a subduction merges with a larger subduction into a single downwelling.

In the 4 cases, the errors on the estimated temperature field systematically decrease during the analysis step for the EnKF algorithm. The errors stay below the first analysis error for evolutions A, B and C, while they reach slightly higher values for evolution D. The error of the EnKF is always lower than that obtained with method 1 for the first 50 My. The average error is lower for the EnKF than for method 1 in 3 out of 4 cases. The average standard deviation of the ensemble (ensemble spread) is of the same order of magnitude as the true error. However, its evolution is not the same as the true error, with differences between both of more than 2% for some part of evolution C, for example. For evolutions C and D, the results of the two methods are comparable whereas the assimilation with EnKF performs better than method 1 for evolutions A and B.

For evolution B, method 1 fails to reconstruct accurately the evolution, with the error reaching values greater than 10% at the end of the assimilation. This case is further investigated on Figures 4.4

Effect of the data assimilation parameters on the quality of the estimation

As shown in Figure 4.1, the choice of N , ℓ v and ℓ h is critical to minimize errors in the assimilation, with errors on the estimated temperature field varying from 4% to more than 10% according to the choice of parameters. We investigate further the effect of these parameters by comparing the average global errors after analyses, ǫ a T , for different combinations of N , ℓ v and ℓ h . Figure 4.6 displays the values of ǫ a T for sizes of ensemble N = 96, 288 and 768 (Figures 4.6A, B and C respectively) with ℓ v varying between 0.3 and 1, and ℓ h between π/10 and π/2. As in Figure 4.1, we observe a dichotomy between assimilations with N = 96 members, with higher errors, and assimilations with N = 288 and 768, with lower errors.

For each size of ensemble N we identify the pair (ℓ v , ℓ h ) that leads to the assimilation with the lowest error ǫ a T min (N ). From this minimum value ǫ a T min (N ), we select all the pairs (ℓ v , ℓ h ) that lead to data assimilation with global errors less than ǫ a T min (N ) + 0.2%. As the size of the ensemble increases, the optimal lengths of correlations (ℓ v , ℓ h ) tend to increase. This is a classical effect [START_REF] Houtekamer | Data assimilation using an ensemble kalman filter technique[END_REF], observed in ensemble Kalman filters for various dynamical systems. As N increases, the amplitude of noise in the sample correlation matrix P f decreases, and small, yet real, correlations between distant points can be taken into account [START_REF] Hamill | Distance-dependent filtering of background error covariance estimates in an ensemble kalman filter[END_REF]. Between ensemble sizes of N = 96 and N = 288 the zone of optimal correlations is displaced towards the greater vertical correlation lengths. When we increase the size of the ensemble from N = 288 to N = 768, the zone of optimal correlations is displaced towards greater horizontal correlation angles. So the accurate estimation of correlations between points on the same vertical level needs less samples than between points on the same horizontal level. This is due to the specifics of mantle convection dynamics. The highly nonlinear rheology produces plates at the surface with values of velocity and temperature that may vary substantially (by one or two orders of magnitude) on short distances in the horizontal direction, especially because of pseudoplasticity. On the contrary, highly viscous cold downwellings establish a strong continuity in the vertical direction. For the ensemble size N = 288 and all the values of (ℓ v , ℓ h ), we additionally evaluate the average global ensemble spread (4.35) the average forecast error on data

σ a T = 1 K K k=1 σ a T (k),
ǫ f y = 1 K K k=1 d k y o k (4.36)
and the norm of the cumulative innovation for K analyses:

d K = 1 K K k=1 y o k -Hx f k . (4.37)
These three values are indicators of the accuracy of the assimilation and can be computed in the case of an assimilation with Earth data, unlike ǫ a T .

Figure 4.7 represents these results along with the true error ǫ a T . The ensemble of optimal data assimilation parameters is also outlined (ǫ a T < ǫ a T min (N ) + 0.2%). Overall, the average ensemble spread σ a T (Figure 4.7B) decreases when ℓ h and ℓ v increase, with a minimum for ℓ h = π/2 and ℓ v = 1. The higher the correlation lengths, the more covariances will be taken into account in the analysis, and the analyzed members will be closer to each others and σ a T lower. The average ensemble spread σ a T is of the same order of magnitude as the true error ǫ a T . Moreover, there is a local minimum of σ a T at ℓ v = 0.7 and ℓ h = π/10. These parameters correspond to the minimum true error ǫ a T .

The average forecast errors and the norm of the cumulative innovations display the same behavior: they decrease with increasing vertical and horizontal correlation lengths.

The longer the correlations lengths, the closer the forecast data are to the observations, and the less biased the assimilation. This means that a better fit to the observations does not necessarily imply a better fit to the true temperature field. In a realistic context, the result of the assimilation should be checked against independent data to evaluate its accuracy. In the case of the Earth's mantle, independent data could be for example the geoid or tomographic models.

Accuracy of the Reconstruction of Geodynamic Structures

We focus on three key flow structures: 1) downwelling slabs (subduction) 2) ridges, i.e. shallow structures resulting from divergent plates at the surface, 3) plumes, hot upwellings raising from the base of the model. Globally, the EnKF288 and EnKF96 solutions for the temperature field are smoother than the solution of method 1. We observe this difference especially in the asthenosphere, the part of the mantle below the top boundary layer. For method 1, the asthenosphere shows short wavelength temperature variations. These variations are absent from the true temperature field and are inconsistent with convection solutions with the chosen parameters. They stem from the amplification of the noise in the observations during the analysis. Moreover, the asthenosphere of the analyzed temperature field of method 1 is hotter than the true temperature.

Both EnKF96 and EnKF288 reconstruct successfully the ridges locations and structures, as testified by their error fields. On the contrary, method 1 fails to reconstruct the ridge on the top right of the domain. It also predicts a ridge that does not exist in the true state (in the top left quadrant). On the right of the domain, another ridge is associated with a vertical positive temperature anomaly underneath. This pattern is found regularly under ridges when applying method 1. This is due to the use of a constant forecast error covariance matrix, P f 0 for the analysis. This constant matrix does not take into account the specifics of the dynamics under a ridge, where the positive anomaly is generally shallow. We do not observe this detrimental effect in the EnKF assimilations, where we compute the forecast error covariance matrix P f k at each analysis time from the forecast ensemble.

All three assimilations reconstruct the subductions and predict accurately the bending direction of slabs at the base of the model. Method 1 tends to underestimate the amplitude of the negative temperature anomalies whereas both EnKF assimilations overestimate them. This is especially noteworthy for the bottom left subduction. Moreover, the estimated slabs are wider than the true slabs. However, we note two arguments in favor of the EnKF: first, the estimation of the slab improves when the size of the en- that the estimation in this part of the domain is less accurate.

Both EnKF288 and EnKF96 solutions do not show any plume at the base of the mantle. However, the ensemble spread shows a greater uncertainty on the places where plumes occur. Method 1 predicts the approximate location of all plumes, but their geometry is not accurate. Method 1 provides only one estimate of the temperature field.

In this evolution, the plumes are allowed to develop. EnKF96 and EnKF288 provide an ensemble of states. Each state develops plumes at different locations and their averages show only a slightly hotter anomaly over a wide area of possible location for the plumes, as we showed earlier in Figure 4.5 for another assimilation.

To illustrate how different flow structures are reconstructed, we represent on Figure 4.9 the evolution of temperature through time of the ensemble members of EnKF288 at points on the same vertical for a subduction, a plume, a ridge initiation and a stable ridge. Figure 4.8 shows the location of these geodynamical features on the true temperature field. We plot the temperature evolutions at the surface, mid-mantle and at the bottom of the domain. Note that the surface and bottom values of temperature actually correspond to the values of the first points below the surface and above the bottom of the domain, respectively.

At the surface, the temperature is corrected accurately at each analysis, with a difference between the true temperature and the analyzed temperature of less than 0.01. The correction associated with the analysis gradually decreases with depth due to both covariance localization and the dynamics of the system. After 70 My (i.e. one transit time), the true value of temperature falls within the range of values predicted by the ensemble for all geodynamical contexts and all depths.

For the subduction, the correction is first done on the surface, and then propagates gradually in depth. The reconstruction of mid-mantle temperature becomes accurate after 40 My, and at the bottom of the model after 70 My, which is the value of the transit time. At the surface, the spread of the ensemble decreases as more data are assimilated.

On the contrary, the spread of the ensemble remains steady for mid-mantle depths and at the bottom of the domain. For these depths, only the average temperature varies. At the surface for the plume, the spread of the ensemble is very low except for a peak at 40 My, which corresponds to an instability, corrected after one analysis. We note that this instability affects greatly method 1 since it leads to the false prediction of the ridge seen in Figure 4.8. At mid-mantle, the ensemble average is slowly converging to the true temperature. At the bottom, the estimated temperature is lower than the true temperature, although it slightly increases throughout the assimilation.

The ridge initiation shows how new observations affect the spread of the ensemble. At the surface, the spread of the ensemble remains low until 100 My, the time of initiation of the ridge. From then on, the estimated temperature increases and the ensemble members follow the cycle of increasing spread during forecast and dramatic decrease of spread during analysis. The temperature in the mid-mantle is estimated with a very good accuracy after 50 My. On the contrary, the assimilation does not predict the evolution of the temperature at the bottom of the domain, although the true temperature falls within the zone defined by the standard deviation of the ensemble after 50 My.

For the stable ridge, the spread of the ensemble at the surface is increasing during forecast and decreasing dramatically during the analysis. At mid-mantle, the estimated temperature becomes accurate after 100 My. At the bottom of the domain the temperature is underestimated although it follows the variations of the true temperature:

increase of temperature at the beginning of the assimilation and slight decrease at the end of the assimilation.

Discussion

We chose the Ensemble Kalman Filter method for its ease of implementation and flexibility to adapt to different forward numerical models. Indeed, as long as the nature of the state and observations does not change, the computation of the analysis step remains the same regardless of the convection code used. On the contrary, the alternative method, variational data assimilation, requires the development of an adjoint code that needs further development for each additional complexity added to the forward model (see Kalnay et al., 2007, for a comparison of EnKF and 4D variational methods). For the mantle circulation problem, this results in a series of derivation of the adjoint considering different approximations (Ismail-Zadeh et al., 2003a;[START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF]Ghelichkhan & Bunge, 2016;[START_REF] Worthen | Towards adjointbased inversion for rheological parameters in nonlinear viscous mantle flow[END_REF]. The ability of a data assimilation scheme to adapt to different numerical codes is a particularly important issue for mantle convection since models are in constant evolution, with current developments including the implementation of chemistry, nonlinear rheologies, elasticity, phase transition and compressibility (see e.g. [START_REF] Zhong | Numerical methods for mantle convection[END_REF], for a review of recent developments of mantle convection codes). In particular, this ease of implementation allows us to work on models producing self-consistently plate-like tectonics at their surface, and hence to obtain forecasts whose data can be ultimately compared with plate reconstructions.

The application of the Ensemble Kalman Filter to the mantle circulation problem is the continuation of the simpler sequential filter that we developed in an earlier work [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF]. The main difference between the two filters is that the EnKF evaluates the state covariance matrix with an ensemble of members. This ensemble approach allows the nonlinear evolution of errors during the forecast stage. This leads to a higher precision in the reconstruction, but also to a more robust scheme, able to reconstruct evolutions which could not be reconstructed with the former method. Moreover, the ensemble assimilation provides an estimate of the errors on the reconstruction at each point of the domain. The estimation of errors could be a valuable information for plate tectonic reconstructions for regions and times where data are scarce, by showing the possible alternative scenarios supported by the ensemble.

This gain in information and quality for reconstructions comes with a price. While we could perform the former assimilation method in one core hour, the method developed here requires several hundreds to several thousands of core hours. However, an efficient parallelization using the PDAF software [START_REF] Nerger | Software for ensemble-based data assimilation systems-implementation strategies and scalability[END_REF] in combination with the parallel code STAGYY produces a highly parallel ensemble filter, able to perform the assimilations on 768 cores in 20 min for an ensemble of 96 members and 3 hours for an ensemble of 768 members.

The important computational cost of the EnKF limited us in the number of assimilations we could test. After checking the stability of the assimilation results on four different evolutions, we chose to focus on studying the effect of the parameters of the Ensemble data assimilation: the size of the ensemble and the vertical and horizontal correlation lengths. The optimum size of the ensemble for our problem is of the order of 288 members. Indeed, almost tripling the number of members leads to a decrease of the average error of less than 0.2%, and on the contrary, dividing the size of the ensemble by 3 leads to an increase of the average error of more than 1%. Although these differences in errors appear to be small, they affect the quality of the reconstruction of thermal structures, as is illustrated in figure 4.8. The average errors on the temperature field for the estimates shown in figure 4.8 range between 4.8 and 5.8%, so the difference in errors are less than 1%. Locally, this translates into the presence (or absence) of geodynamic structures (like ridges and upwellings) which are artefacts. Covariance localization proved to be important to minimize the error in the reconstruction of mantle structure: for 288 members, the difference in the average error is of 0.85% between the optimal correlation length and the least favorable one. We also investigated the statistics of the cumulative innovation and of the average forecast error for different ensemble sizes and correlation lengths. In a realistic case, these are the only variables available to evaluate the quality of the assimilation. The variation of both cumulative innovation and average forecast error as a function of ensemble size show the same tendency as the average error on the temperature field. On the contrary, the correlation lengths minimizing the norm of the cumulative innovation and the average forecast error were different from the ones minimizing the error on the temperature field.

This shows the limits of these indicators to determine the optimal parameters for the assimilation. In a realistic case, rigorous a posteriori evaluation of a data assimilation result would require comparison of the prediction made with independent observations [START_REF] Talagrand | Errors. a posteriori diagnostics[END_REF]. For mantle circulation, seismic tomography, topography, true polar wander or the geoid could play this role.

By construction, sequential data assimilation methods do not propagate new information back in time. In the case of the reconstruction of mantle circulation, this is a clear disadvantage since the information on the Earth's surface tectonics tends to become more reliable as we get closer to present-day. Consequently, a natural extension of the present work would be to implement an Ensemble Kalman Smoother [START_REF] Evensen | An ensemble kalman smoother for nonlinear dynamics[END_REF][START_REF] Van Leeuwen | An ensemble smoother with error estimates[END_REF]. In the same way as the EnKF uses sample spatial correlations of the ensemble to update the state of the system with new observations, the Ensemble Kalman smoother uses sample time and space correlations of successive ensembles to update former states with the new observations. Evensen (2003) shows how an Ensemble Kalman Smoother can be implemented with a minimal computational cost alongside a preexisting EnKF. Moreover, [START_REF] Nerger | On the influence of model nonlinearity and localization on ensemble kalman smoothing[END_REF] shows that such algorithm is efficient for nonlinear models, and that in their test case, optimal localization parameters for the Ensemble Kalman Smoother coincide with optimal localization parameters for the EnKF.

As a first approach to test the EnKF for mantle circulation reconstructions, we chose a fairly simple convection model. As already discussed in Section 4.2.1, a more realistic mantle model would have, among other things, a 3D-spherical shell geometry and a higher Rayleigh number. This would substantially increase the size of the data assimilation problem. However, we followed the procedure as described in [START_REF] Nerger | Software for ensemble-based data assimilation systems-implementation strategies and scalability[END_REF] to implement the EnKF. This results in a highly scalable filter, enabling the computation of the EnKF assimilation in a reasonable time. An increase in the Rayleigh number also implies thiner boundary layers, slabs and plumes. This could translate into lower optimum correlation lengths for the EnKF. A more realistic model would additionally include a viscosity increase in the lower mantle [START_REF] Ricard | A geodynamic model of mantle density heterogeneity[END_REF], and the presence of continents. This would tend to lengthen the wavelength of convection in the lower mantle and therefore might ease the mantle circulation reconstruction (see for example Ricard, 2015, sections 7.02.6.3.2 and 7.02.6.7 for a discussion of both effects on mantle convection).

In the synthetic experiments of section 4.4, the convection model used to produce the series of data is the same as the forward model used during the assimilation. For an application with Earth data, this will not be the case. The equations solved in models of mantle convection still hold some shortcomings [START_REF] Ricard | Physics of mantle convection[END_REF]. Moreover, theories, observations and experiments do not yet fully constrain parameters, especially rheological ones [START_REF] King | Reconciling laboratory and observational models of mantle rheology in geodynamic modelling[END_REF], and variations in rheology affect the reconstructions of mantle circulation [START_REF] Bello | Assessing the role of slab rheology in coupled plate-mantle convection models[END_REF]. Hence it could be fundamental to take into account model errors. A first order solution is to increase the inflation parameter γ in Equation 4.20: this would overall increase the a priori uncertainty on the mantle estimation.

Performing experiments where the model used to compute the observation is different from the model used for the assimilation would provide us with more information on how to implement model errors. Another solution would be to consider the joint assimilation of the state and model parameters. Although it is in principle possible for the EnKF (Evensen, 2009b), it could be computationally not tractable. Indeed, the response of mantle dynamics to different rheological parametrization is highly nonlinear, and their inversion calls for the development of techniques focusing on rheology, such as adjoint based inversions of rheological parameters [START_REF] Worthen | Towards adjointbased inversion for rheological parameters in nonlinear viscous mantle flow[END_REF][START_REF] Ratnaswamy | Adjoint-based estimation of plate coupling in a non-linear mantle flow model: theory and examples[END_REF] or further applications of the recently developed pattern recognition techniques for mantle convection [START_REF] Atkins | Using pattern recognition to infer parameters governing mantle convection[END_REF].

The choice of the synthetic experiments assimilation window of 150 Myr is a compromise between having the possibility to compute assimilations for various cases and having an assimilation window covering most part of the timespan of plate tectonic reconstructions [START_REF] Seton | Global continental and ocean basin reconstructions since 200ma[END_REF][START_REF] Müller | Ocean basin evolution and global-scale plate reorganization events since pangea breakup[END_REF][START_REF] Torsvik | Plate tectonics and net lithosphere rotation over the past 150my[END_REF]. A real assimilation could take into account a longer timespan and therefore improve the assimilation results. However, the structure of the dataset used for the synthetic experiments is a very idealized version of the actual plate reconstruction models. We already discuss this issue in [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF]. In the following, we supplement and update this discussion in the light of research that has recently come to the fore.

First, we set a time series of data covering the whole surface of the domain and regularly available, every 10 Myr. Plate tectonic reconstructions data are more complex.

They are based on the estimation of finite relative rotations between individual plates, structured into a hierarchy describing global relative motions and anchored in an absolute reference frame. The span of each finite relative rotation is determined depending on the amount and quality of information available for a specific context and therefore varies depending on plate pairs and times. The average span of finite rotations of recent plate models is of the order of 10 Myr [START_REF] Torsvik | Plate tectonics and net lithosphere rotation over the past 150my[END_REF] to 5 Myr [START_REF] Müller | Ocean basin evolution and global-scale plate reorganization events since pangea breakup[END_REF], but varies over time with for example 1 Myr resolution for the last 20 Myr in some regions (Merkouriev & DeMets, 2014), or some gaps in the data such as during the cretaceous superchron from 121 to 83 Myr ago [START_REF] Granot | Geomagnetic field variability during the cretaceous normal superchron[END_REF]. The continuously closed plate algorithm [START_REF] Gurnis | Plate tectonic reconstructions with continuously closing plates[END_REF] produces plate tectonic reconstruction maps continuous in space and time which allows the creation of a series of global plate reconstructions at regular intervals. Nonetheless, creating such a regularized time series of reconstruction might miss tectonic events. Instead, we could adapt the frequency of analyses to the varying plate reconstruction resolution. Additional synthetic experiments with a time-series of data whose frequency evolves through time are necessary to explore the limits of such method.

Second, the observations were perturbed independently with a Gaussian noise of 10% of the respective root mean square value of surface heat flux and surface velocities.

The estimation of uncertainties on absolute plate motion models involves estimation of both uncertainties in relative plate motion and on the absolute reference frame (Müller & Wessel, 2015). The main source of information on the motion of plates comes from the map of seafloor magnetic anomalies. [START_REF] Hellinger | The uncertainties of finite rotations in plate tectonics[END_REF] developed a method to compute relative motion of plates and associated uncertainties inferred from magnetic anomaly identifications. Recently, [START_REF] Seton | Community infrastructure and repository for marine magnetic identifications[END_REF] built an open source community database. It gathers seafloor magnetic anomaly identifications, and estimation with [START_REF] Hellinger | The uncertainties of finite rotations in plate tectonics[END_REF] method of plates relative motion and associated uncertainties. This database could be used in the future as a basis to automatically produce global plate motion histories and assess their uncertainties. To our knowledge, this has not been done so far at a global scale. On a regional scale and for recent time (5 to 20 Myr),

Conclusion

We applied the Ensemble Kalman Filter algorithm to the reconstruction of mantle circulation through time. We chose a formulation with covariance inflation and localization to minimize the effect of sampling errors in the estimation of the forecast error covariance matrix. Synthetic "twin" experiments with different evolutions and for different parameters allowed us to assess the efficiency of the algorithm and to determine the optimal parameters for the assimilation.

This work builds on the developments of a first approach to sequential data assimilation for mantle circulation made in [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF]. The EnKF is more robust and on average more accurate than the former method. Additionally, the Ensemble Kalman Filter provides not only an estimate of mantle circulation, but also detailed maps of uncertainties on this estimation.

We evaluate the accuracy of the EnKF as a function of three main parameters: the size of the ensemble, and two covariance localization parameters, namely the vertical correlation length and horizontal correlation angle. We find that the optimal size of the ensemble is of the order of 300 members. For this ensemble size, the optimal vertical correlation length corresponds to two thirds of the domain thickness, and the optimal horizontal correlation angle is of π/10 (around 2000 km). These values should be reevaluated as the dynamical model becomes more realistic.

The EnKF was implemented using the parallel data assimilation framework PDAF in a preexisting mantle convection code, STAGYY. The resulting code is highly scalable, which means that the application of the EnKF to realistic data assimilation with plate reconstructions and a 3D spherical mantle model is within reach in a foreseeable future.

Résumé du chapitre

Les récentes avancées en modélisation de la convection mantellique ont conduit à l'apparition d'une nouvelle génération de codes de convection, capables de générer de façon auto-cohérente des tectoniques de type plaque à leur surface. Ces modèles relient physiquement la dynamique du manteau à la tectonique de surface. Associés à des reconstructions de la tectonique des plaques, ils ont le potentiel de produire une nouvelle génération de modèles de circulation mantellique qui utilisent des méthodes d'assimilation de données et tiennent compte des incertitudes sur les reconstructions de la tectonique des plaques. Nous avons récemment fourni une preuve de ce concept en appliquant un filtre de Kalman sous-optimal à la reconstruction de la circulation mantellique [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF] CHAPTER 5

Discussion and Conclusion

During this thesis, we successively proposed two algorithms to reconstruct mantle circulation by combining observations on surface tectonics with a mantle convection model, using sequential data assimilation methods.

We first developed a suboptimal sequential data assimilation scheme based on the Kalman filter (Chapter 3). This scheme is efficient in recovering the temperature field of a convective system with plate-like tectonics at its surface over several 100 Myr. We tested the robustness of the method by conducting synthetic experiments in 2D spherical annulus geometry. We obtained accurate results, even for the location of plumes at the base of the models, for periods of at least 300 Myr, provided the time between analyses is shorter than 15 Myr and the noise in observation is lower than 30%.

These encouraging results motivated us to further explore the data assimilation problem of reconstructing mantle circulation from surface tectonics (Chapter 4). The Ensemble Kalman Filter seemed to be a method of choice given its ability to deal with high dimensional problems and nonlinear models. It indeed proved to be more robust and on average more accurate than the former method, while additionally providing detailed maps of uncertainties on the estimated evolution.

The synthetic experiments we performed in both studies are by no means exhaustive, and many more tests remain to be done to thoroughly evaluate data assimilation methods in the context of mantle circulation reconstruction. Among them, we identify the need for:

• further twin experiments using a more realistic setting:

from the dynamical model side, this means, among other things, tests in 3D spherical geometry, the introduction of continents, a higher Rayleigh number, a viscosity jump at mid-mantle, but also taking into account the chemical aspects of mantle convection (see for example Ricard, 2015, section 6 for an overview of the specificities of mantle convection).

from the observed data side, this means using an observation error covariance matrix which takes into account the correlations between errors on adjacent points, and the correlation between errors on observed surface velocities and observed surface heat fluxes. For models developing small-scale convection, It will also be necessary to use directly seafloor age maps instead of surface heat fluxes as observed data.

• the design of "sister experiments", to assess the effect of errors in mantle convection models. The principle of these sister experiments would be to compute a reference evolution with a given set of model parameters, and then to apply the data assimilation algorithm using a different set of parameters. This type of experiment will inform us on the impact an error in the forward model will have in the estimation of mantle circulation. This information could be of prior importance to succeed in the application of these data assimilation algorithms to Earth data.

Of course, the ultimate test for a data assimilation scheme is its application to Earth data. In chapter 4, we developed a data assimilation code for mantle circulation reconstruction by combining a parallel ensemble data assimilation framework (PDAF [START_REF] Nerger | The parallel data assimilation framework: experiences with kalman filtering, in Use of High Performance in Meteorology[END_REF][START_REF] Nerger | Software for ensemble-based data assimilation systems-implementation strategies and scalability[END_REF] with the state-of-the-art parallel mantle convection code STAGYY (see [START_REF] Tackley | Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid[END_REF]. The high scalability of the resulting code, and the modelling capabilities of STAGYY (see Section 1.5) make the application of data assimilation to real data technically possible. The computation of a mantle circulation reconstruction will further require:

• a time series of the layout of plates and associated plate motions. This time series is easily accessible thanks to the open-access policy of teams working on plate tectonics reconstructions (e.g. [START_REF] Müller | Ocean basin evolution and global-scale plate reorganization events since pangea breakup[END_REF].

• a comprehensive assessment of the errors on plate tectonics reconstructions, taking into account the errors in both relative and absolute plate motions, and consistent with the formalism used in data assimilation. This task can only be achieved by a tight collaboration with teams working on plate tectonics reconstructions.

As highlighted in section 2.2.2, tools have been developed to estimate the uncertainties on relative and absolute plate motions (see e.g. Wessel & Müller, 2015 for a recent review on these techniques). Hence, a first order estimation of uncertainties in plate tectonics reconstructions is possible, even though, as we go back in time, this estimation might not reflect true uncertainties on plate tectonics reconstructions (see discussion at the end of Section 2.2.2).

• choices in the parametrization for the mantle convection model. This choice of parameters will be a compromise between obtaining an internal dynamics and a surface tectonics as close to the Earth as possible, and the limits in computational power. In section 1.4, we argued for the use of a composite Newtonian and pseudoplastic rheology to obtain a mantle convection model self-consistently producing plate-like tectonics at its surface. This type of models has been investigated in pioneering studies using 3D spherical geometry (e.g. [START_REF] Richards | Role of a lowviscosity zone in stabilizing plate tectonics: Implications for comparative terrestrial planetology[END_REF][START_REF] Walzer | Mantle convection and evolution with growing continents[END_REF][START_REF] Van Heck | Planforms of self-consistently generated plates in 3d spherical geometry[END_REF][START_REF] Yoshida | Mantle convection with longest-wavelength thermal heterogeneity in a 3-d spherical model: Degree one or two?[END_REF][START_REF] Foley | Generation of plate-like behavior and mantle heterogeneity from a spherical, viscoplastic convection model[END_REF]. Further work has also involved the statistical comparison of the surface dynamics of this type of models with plate tectonics reconstructions [START_REF] Coltice | Dynamic causes of the relation between area and age of the ocean floor[END_REF][START_REF] Coltice | Convergence of tectonic reconstructions and mantle convection models for significant fluctuations in seafloor spreading[END_REF][START_REF] Rolf | Statistical cyclicity of the supercontinent cycle[END_REF][START_REF] Mallard | Subduction controls the distribution and fragmentation of earth's tectonic plates[END_REF]. Finally, they have recently been used to reconstruct mantle circulation using the semi-empirical sequential method [START_REF] Bello | Assessing the role of slab rheology in coupled plate-mantle convection models[END_REF]Coltice, in prep). These studies provide a solid basis to design a mantle convection model for data assimilation applications.

To evaluate the result of such mantle circulation reconstructions, we will need to test it against independent data. We propose two types of tests:

1. comparison of mantle circulation reconstructions against data independent from plate tectonics reconstructions. At present, these include seismic tomography (see Section 2.2.1), dynamic topography and the geoid (see e.g. [START_REF] Steinberger | Topography caused by mantle density variations: observationbased estimates and models derived from tomography and lithosphere thickness[END_REF], and measurement of heat fluxes at the surface [START_REF] Jaupart | Heat flow and thermal structure of the lithosphere[END_REF]. In the past, this could include sea level changes (Müller et al., 2008b) and evidences of vertical motions [START_REF] Flament | A review of observations and models of dynamic topography[END_REF].

2. Another possible test is to perform data assimilation over a restricted timespan to obtain an estimation of the state of the mantle at some time in the past. Then, we can compute a free evolution (i.e. unconstrained by data) of this estimated state up to the present day. By comparing the free evolution with observations, we can assess the predicting power of mantle circulation reconstructions.

Section 2.3 gives a partial overview of the possible applications of mantle circulation reconstructions: studies on the effect of surface tectonics on the deep mantle structure (e.g. [START_REF] Mcnamara | Thermochemical structures beneath africa and the pacific ocean[END_REF]Schuberth et al., 2009b,a;[START_REF] Zhang | Heat fluxes at the earth's surface and core-mantle boundary since pangea formation and their implications for the geomagnetic superchrons[END_REF][START_REF] Davies | Reconciling dynamic and seismic models of earth's lower mantle: The dominant role of thermal heterogeneity[END_REF][START_REF] Bower | Lower mantle structure from paleogeographically constrained dynamic earth models[END_REF][START_REF] Bull | The effect of plate motion history on the longevity of deep mantle heterogeneities[END_REF][START_REF] Hassan | Provenance of plumes in global convection models[END_REF], and investigation of alternative scenarios of absolute plate motions (e.g. [START_REF] O'neill | On the uncertainties in hot spot reconstructions and the significance of moving hot spot reference frames[END_REF][START_REF] Doubrovine | Absolute plate motions in a reference frame defined by moving hot spots in the pacific, atlantic, and indian oceans[END_REF][START_REF] Shephard | Testing absolute plate reference frames and the implications for the generation of geodynamic mantle heterogeneity structure[END_REF][START_REF] Nerlich | Mantle convection models constrain central neo-tethys ocean reconstructions[END_REF]Zahirovic et al., in press). Other applications include studies of changes in the Earth's rotation axis [START_REF] Steinberger | Changes of the earth's rotation axis owing to advection of mantle density heterogeneities[END_REF], the sea-level [START_REF] Moucha | Dynamic topography and long-term sea-level variations: There is no such thing as a stable continental platform[END_REF], dynamic topography [START_REF] Flament | A review of observations and models of dynamic topography[END_REF], and the effect of mantle circulation on core dynamics [START_REF] Olson | Core evolution driven by mantle global circulation[END_REF].

Using data assimilation allows the integration of tectonic data into convection models in a more consistent way, since the convection model we use generates naturally plate-like tectonics at its surface, and plate tectonics reconstructions are not enforced through kinematic boundary conditions, but taken into account in a statistically consistent fashion.

Compared to other disciplines, data assimilation for mantle convection is still in its infancy. Knowledge, experience and technical advances of data assimilation in other fields, especially atmospheric sciences and oceanography, provide an extensive amount of resources to design algorithms for the mantle circulation reconstruction problem. A direct extension of our work would be for example to develop an ensemble Kalman smoother algorithm [START_REF] Evensen | An ensemble kalman smoother for nonlinear dynamics[END_REF][START_REF] Van Leeuwen | An ensemble smoother with error estimates[END_REF] to better estimate past mantle circulation using data on more recent states of the mantle. New methods are also being developed to take better into account the nonlinearities in dynamical models and the resulting nongaussianity in the statistics of errors [START_REF] Bocquet | Beyond gaussian statistical modeling in geophysical data assimilation[END_REF]. Such solutions should be considered to build and improve data assimilation algorithms for mantle circulation.

We started this dissertation by remarking the coexistence of two different theories that set the framework for understanding the dynamics of the Earth's mantle and crust system: plate tectonics and mantle convection. The successive data assimilation developments that we proposed are an attempt to combine the products of these two theories: tonics theory provided such global framework and proved to be accurate at first order in its description of surface dynamics. However, plate tectonics is, in essence, a kinematic theory: it describes the motions at the surface of the Earth, but not their cause.

Moreover, it is an idealized model of surface motions: plates are not perfectly rigid (e.g. [START_REF] Gordon | The plate tectonic approximation: Plate nonrigidity, diffuse plate boundaries, and global plate reconstructions[END_REF], and plate boundaries are not infinitely small but can be on the contrary diffuse, in the oceans [START_REF] Gordon | Diffuse oceanic plate boundaries: strain rates, vertically averaged rheology, and comparisons with narrow plate boundaries and stable plate interiors, The history and dynamics of global plate motions[END_REF][START_REF] Wiens | A diffuse plate boundary model for indian ocean tectonics[END_REF] as well as in the continents, as the presence of massive mountain ranges attest. Mantle convection theory provides, on the contrary, a global physical framework which takes into account those diffuse and continuous deformations, and describes both the solid Earth inner dynamics and surface tectonics. Many challenges remain before we can directly integrate raw observations in mantle convection models to produce a physically consistent reconstruction of both the Earth's mantle and surface tectonics. Undoubtedly, data assimilation methods will be one of the pieces of this puzzle. Dans la section 1.4, nous avons plaidé en faveur de l'utilisation d'une rhéologie composite newtonienne et pseudoplastique pour obtenir un modèle de convection mantellique produisant de manière auto-cohérente une tectonique de type plaque à sa surface. Ce type de modèles a été étudié dans des études pionnières utilisant la géométrie sphérique 3D (par exemple [START_REF] Richards | Role of a lowviscosity zone in stabilizing plate tectonics: Implications for comparative terrestrial planetology[END_REF][START_REF] Walzer | Mantle convection and evolution with growing continents[END_REF][START_REF] Van Heck | Planforms of self-consistently generated plates in 3d spherical geometry[END_REF][START_REF] Yoshida | Mantle convection with longest-wavelength thermal heterogeneity in a 3-d spherical model: Degree one or two?[END_REF][START_REF] Foley | Generation of plate-like behavior and mantle heterogeneity from a spherical, viscoplastic convection model[END_REF].

Discussion et Conclusion

Des travaux ultérieurs ont également porté sur la comparaison statistique de la dynamique de surface de ce type de modèles avec des reconstitutions de tectonique des plaques [START_REF] Coltice | Dynamic causes of the relation between area and age of the ocean floor[END_REF][START_REF] Coltice | Convergence of tectonic reconstructions and mantle convection models for significant fluctuations in seafloor spreading[END_REF][START_REF] Rolf | Statistical cyclicity of the supercontinent cycle[END_REF][START_REF] Mallard | Subduction controls the distribution and fragmentation of earth's tectonic plates[END_REF]. Enfin, ils ont été récemment utilisés pour reconstruire la circulation mantellique en utilisant la méthode séquentielle semi-empirique [START_REF] Bello | Assessing the role of slab rheology in coupled plate-mantle convection models[END_REF]Coltice, in prep). Ces études fournissent une base solide pour concevoir un modèle de convection mantellique adapté à l'application d'algorithmes d'assimilation de données.

Pour évaluer les reconstitutions de la circulation mantellique obtenues par assimilation de données, nous devrons les comparer avec des données indépendantes. Nous proposons deux types de tests:

1. la comparaison des reconstitutions de la circulation mantellique avec des données indépendantes des reconstitutions de la tectonique des plaques. Pour la circulation mantellique actuelle, on peut envisager par exemple les modèles de tomographie sismique (voir section 2.2.1), la topographie dynamique et le géoïde (voir, par exemple Steinberger, 2016), et les mesures de flux de chaleur à la surface [START_REF] Jaupart | Heat flow and thermal structure of the lithosphere[END_REF]. Dans le passé, les données indépendantes pourraient inclure les changements du niveau marin (Müller et al., 2008b) et les indices de mouvements verticaux [START_REF] Flament | A review of observations and models of dynamic topography[END_REF].

2. Un autre test possible consiste à effectuer l'assimilation de données sur un intervalle de temps restreint pour obtenir une estimation de l'état du manteau à un moment donné dans le passé. Ensuite, nous pouvons calculer une évolution libre (c'est-à-dire non contrainte par les données) de cet état estimé jusqu'à ce jour.

En comparant l'évolution libre avec les observations, nous pourrons évaluer le pouvoir prédictif des reconstitutions de la circulation mantellique.

La section 2.3 donne un aperçu partiel des applications possibles des reconstitutions de la circulation mantellique: étude de l'effet de la tectonique de surface sur la structure du manteau inférieur (e.g. [START_REF] Mcnamara | Thermochemical structures beneath africa and the pacific ocean[END_REF]Schuberth et al., 2009b,a;[START_REF] Zhang | Heat fluxes at the earth's surface and core-mantle boundary since pangea formation and their implications for the geomagnetic superchrons[END_REF][START_REF] Davies | Reconciling dynamic and seismic models of earth's lower mantle: The dominant role of thermal heterogeneity[END_REF][START_REF] Bower | Lower mantle structure from paleogeographically constrained dynamic earth models[END_REF][START_REF] Bull | The effect of plate motion history on the longevity of deep mantle heterogeneities[END_REF][START_REF] Hassan | Provenance of plumes in global convection models[END_REF], et étude de scénarios alternatifs de mouvements absolus des plaques (par exemple [START_REF] O'neill | On the uncertainties in hot spot reconstructions and the significance of moving hot spot reference frames[END_REF][START_REF] Doubrovine | Absolute plate motions in a reference frame defined by moving hot spots in the pacific, atlantic, and indian oceans[END_REF][START_REF] Shephard | Testing absolute plate reference frames and the implications for the generation of geodynamic mantle heterogeneity structure[END_REF][START_REF] Nerlich | Mantle convection models constrain central neo-tethys ocean reconstructions[END_REF]Zahirovic et al., in press). D'autres applications incluent les études des changements dans l'axe de rotation de la terre [START_REF] Steinberger | Changes of the earth's rotation axis owing to advection of mantle density heterogeneities[END_REF], le niveau marin [START_REF] Moucha | Dynamic topography and long-term sea-level variations: There is no such thing as a stable continental platform[END_REF], la topographie dynamique [START_REF] Flament | A review of observations and models of dynamic topography[END_REF] et l'effet de la circulation mantellique sur la dynamique du noyau [START_REF] Olson | Core evolution driven by mantle global circulation[END_REF].

L'utilisation de l'assimilation de données permet l'intégration de données tectoniques dans les modèles de convection d'une manière plus cohérente, puisque le modèle de convection que nous utilisons génère naturellement une tectonique de type plaque à sa surface et que les reconstitutions de la tectonique des plaques ne sont pas imposées par des conditions limites cinématiques, mais intégrées d'une manière statistiquement cohérente.

Comparée à d'autres disciplines, l'assimilation de données pour la convection mantellique est encore à ses débuts. Les connaissances, l'expérience et les progrès techniques de l'assimilation de données dans d'autres domaines, en particulier les sciences de l'atmosphère et l'océanographie, fournissent une quantité importante de ressources pour concevoir des algorithmes adaptés au problème de la reconstitution de la circulation mantellique. Une extension directe de notre travail serait par exemple de développer un lisseur de Kalman [START_REF] Evensen | An ensemble kalman smoother for nonlinear dynamics[END_REF][START_REF] Van Leeuwen | An ensemble smoother with error estimates[END_REF] pour mieux estimer la circulation mantellique passée en utilisant des données sur les états les plus récents du manteau. De nouvelles méthodes sont également développées pour mieux prendre en compte les non-linéarités dans les modèles dynamiques et la nongaussianité qui en résulte dans la statistique des erreurs [START_REF] Bocquet | Beyond gaussian statistical modeling in geophysical data assimilation[END_REF].

Ces solutions devraient être envisagées pour construire et améliorer les algorithmes ne sont pas infiniment petites, mais peuvent être au contraire diffuses, dans les océans [START_REF] Gordon | Diffuse oceanic plate boundaries: strain rates, vertically averaged rheology, and comparisons with narrow plate boundaries and stable plate interiors, The history and dynamics of global plate motions[END_REF][START_REF] Wiens | A diffuse plate boundary model for indian ocean tectonics[END_REF] ainsi que sur les continents, comme en témoignent la présence de massifs montagneux. La théorie de la convection mantellique fournit au contraire un cadre physique global qui tient compte de ces déformations diffuses et continues et décrit à la fois la dynamique interne de la Terre solide et la tectonique de surface. De nombreux défis restent avant que nous puissions intégrer directement les observations brutes dans les modèles de convection mantellique pour produire une reconstitution physiquement cohérente de la circulation dans le manteau terrestre et de la tectonique de surface. Sans aucun doute, les méthodes d'assimilation de données seront l'une des pièces de ce puzzle.

FIGURE 1

 1 FIGURE 1.2: Left: surface velocity and viscosity fields snapshots of a free (i.e. unconstrained by data) mantle convection evolution with composite viscoplastic and temperature-pressure Newtonian rheology. Right: corresponding temperature field of the interior of the system. From Bello et al. (2015)

(

  Ismail-Zadeh & Tackley, 2010 provide a detailed description of numerical methods for geodynamic modelling, and Zhong et al., 2015 review the recent advances of mantle convection numerical modelling).

FIGURE 2

 2 FIGURE 2.2: "World Ocean Floor Panorama" by[START_REF] Heezen | World ocean floor panorama[END_REF], two pioneers of ocean exploration, who produced the first comprehensive maps of the ocean floor.

FIGURE 2 . 3 :

 23 FIGURE 2.3: Diagram of two plates sharing a divergent plate boundary.The motion of plate 2 relative to plate 1 is described by the Euler pole A and an angular velocity. Transform faults develop along small circles of constant latitude with respect to the Euler pole A.[START_REF] Morgan | Rises, trenches, great faults, and crustal blocks[END_REF]

FIGURE 2 . 4 :

 24 FIGURE 2.4: The first consistent global model of plate motions as published by Le Pichon (1968).

FIGURE 2

 2 FIGURE 2.5: Map taken from Seton et al. (2014), which shows magnetic anomalies identifications (colored by age), recently collected to build a community database.

FIGURE 2 . 6 :

 26 FIGURE 2.6: Map of gravity anomalies (Sandwell & Smith, 2009) and a model of plates at present. The name of major plates are: AFR: Africa, ANT: Antarctica, ARA: Arabia, AUS: Australia, C: Cocos, CAP: Capricorn, CAR: Caribbean, EUR: Eurasia, IND: India, NAM: North America, NAZ: Nazca, PAC: Pacific, PH: Philippine, SAM: South America, SOM: Somalia. Taken from Seton et al. (2012).

FIGURE 2 . 7 :

 27 FIGURE 2.7: Diagram of the principle of the Hellinger method for relative plate motion reconstructions. Plates A and B are separated by a ridge (double lines), and active transform faults (single line). a) localization of magnetic anomaly identifications of the same age on both sides of the plate boundaries (pink and blue points) and of fossil transform faults (dashed lines) b) definition of isochrons and identification of conjugate great circle segments c) position of magnetic anomaly identifications after the finite rotation, and inferred past plate boundary. Taken from Wessel & Müller (2015)

Figure 2

 2 Figure 2.5. Müller et al. (2008b) produced continuous maps of seafloor age uncertainty based on the distance from the nearest magnetic anomaly.

  form a slight positive temperature anomaly at present into a very large and localized temperature anomaly in the past. It follows that the final solution (the temperature field at some time in the past) does not depend continuously on the initial condition (the temperature field at present): the problem is unstable, and following Hadamard (1902) definition, ill-posed (see Ismail-Zadeh et al., 2016, section 1.1 for a simple 1-D example of this behavior). In practice, the diffusion term is dropped. Then, the backward rate of energy production -R h also needs to be dropped, for global energy balance reasons. The backward conservation energy equation becomes a simple backward advection equation ∂T ∂t = u • ∇T.

  to the backward advection algorithm an analytic solution describing the thermal boundary layers (where diffusion is the principal mechanism of heat transfer) to model more accurately mantle circulation. They test the reliability of the backward advection method by first computing the state of the mantle in the past by backward advection, and then computing a forward evolution from this initial state. Finally, they compare the predicted state of the mantle at present with the state of the mantle derived from seismic tomography. With this setup, they estimate that backward advection is reliable to reconstruct mantle circulation for the last 50 to 75 Myr. O'Neill et al. (2005) push the reconstructions of mantle circulation to 120 Myr, using a different technique for the boundary layers. From their mantle circulation models, they estimate the motions of hotspots through time, in order to build a moving hotspot reference frame (see previous section). They take into account uncertainties in the plate tectonics reconstructions and convection models by successively running backward advection evolutions and modifying the parameters of the model and plate tectonics reconstructions, to converge to a result consistent with both the physics and the data. Doubrovine et al. (2012) also use backward advection to estimate a moving hotspot reference frame for the last 130 Myr. However, they stop the backward advection at 70 Myr, and estimate the motion of hotspots between 70 and 130 Myr using the velocity field at 70 Myr. Williams et al. (2015) compare the moving hotspot reference frame of O'Neill et al. (2005); Doubrovine et al. (2012) to the absolute reference frame recently derived from the identification of paleoslabs in tomography models[START_REF] Van Der Meer | Towards absolute plate motions constrained by lower-mantle slab remnants[END_REF].[START_REF] Williams | Absolute plate motions since 130 ma constrained by subduction zone kinematics[END_REF] find discrepancies between reference frames for times older than 70 Myr, further confirming the results of[START_REF] Conrad | Seismic tomography, surface uplift, and the breakup of gondwanaland: Integrating mantle convection backwards in time[END_REF]: backward advection models are probably not reliable beyond 70 Myr.

  FIGURE 2.9: Maps of the estimated nondimensional temperature at present at 2811 km depth for three mantle circulation models computed from three different initial states, all other conditions being equal (thermochemical mantle model, 119 Myr of assimilation of plate tectonics reconstructions). Initial states are: a) 1-D temperature profile b) 60 Myr of imposed boundary conditions corresponding to the velocity and layout of plates 119 Myr ago c) same as b) but with 120 Myr. Taken from McNamara & Zhong (2005).

  observation operator H K corresponds to the interpolation of the temperature field from the grid defined in mantle convection models to the different discretization used in seismic tomography models. The equivalent discrete objective function becomes:

(

  2014),[START_REF] Vynnytska | Restoring past mantle convection structure through fluid dynamic inverse theory: regularisation through surface velocity boundary conditions[END_REF] andGhelichkhan & Bunge (2016) is to use the adjoint method to obtain a first order estimation of the gradient of the objective function ∇J c around a given initial temperature field. This method requires the derivation of a set of partial differential adjoint equations associated with the set of partial differential equations defining the forward model of mantle convection. The stage 2b corresponds to the backward computation of the discretized adjoint equations. The advantage of this method is that the partial differential adjoint equations are very similar to the forward differential equations, which allows the use of the same code to perform stage 2a and stage 2b. However, this also means that the computational time required for stage 2a is equivalent as the one needed for stage 2b. Moreover, the computation of stage 2b requires to store the whole set of temperature and velocity evolutions computed in stage 2a.[START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF] performs a preliminary test in 3 dimensions for a 100 Myr evolution and obtain convergence of the objective function value after 100 iterations.[START_REF] Horbach | The adjoint method in geodynamics: derivation from a general operator formulation and application to the initial condition problem in a high resolution mantle circulation model[END_REF] apply this method to the reconstruction of mantle circulation for the past 40 Myr. They test the sensitivity of the method to the initial temperature field chosen in stage 1 of the algorithm, and show that the solutions, after several iterations, all converge to the same evolution.[START_REF] Vynnytska | Restoring past mantle convection structure through fluid dynamic inverse theory: regularisation through surface velocity boundary conditions[END_REF] further test this variational algorithm using synthetic experiments in 2D cylindrical geometry. First, they compute a reference evolution using the forward model, with a free-slip boundary condition. The surface velocities and final temperature fields are extracted from the reference evolution. Then the variational algorithm is applied to estimate an evolution from these two constraints: surface velocities are used as boundary conditions, and the final temperature field is the data to fit. The background state is a spherically symmetric temperature field. They limit the timespan of the evolution to 50 Myr. Within this period, the evolution inferred from the variational data assimilation algorithm converges towards the reference solution. They also perform synthetic experiments where the surface boundary conditions are not imposed during the assimilation. In this case, the assimilation algorithm produces an evolution whose final temperature field converges towards the reference final temperature field, but the whole estimated evolution diverges from the reference evolution, showing the presence of local minima of the objective function in Equation2.26 if velocities are not imposed at the surface of the mantle convection model. More generally, these results show the importance of surface velocity constraints to estimate the evolution of the temperature field in the past.Ismail-Zadeh et al. (2003a[START_REF] Steinberger | Prediction of emperor-hawaii seamount locations from a revised model of global plate motion and mantle flow[END_REF]) also use the adjoint method but derive the adjoint equations considering a model restricted to the equations of conservation of energy (and associated boundary conditions), supposing no effect of the temperature variation on the velocity field. In step 2a, they consider the complete mantle model with the equations of conservation of mass, momentum and energy and the associated boundary conditions, arguing that the iteration of step 2a and step 2b will then lead to an estimation of the evolution consistent with the physics of mantle convection. The computational cost of this technique is lighter than the former method and this formulation also allows the implementation of temperature-dependent viscosity.Liu & Gurnis

  .,[START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF]. A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics, Geophysical Journal International, 204(1),[200][201][202][203][204][205][206][207][208][209][210][211][212][213][214]

  .13) ∇ • µ ∇u + (∇u) T -∇p + RaT e r = 0, (3.14)DT Dt = ∇ 2 T + H,(3.15)where u, p, T and t are the non-dimensional velocity, dynamic pressure, temperature and time, respectively. We work in spherical coordinates (r, θ, φ) of unit vectors (e r , e θ , e φ ). Ra is the Rayleigh number and H is the non-dimensional internal heating rate. The models presented here have 10% basal heating and 90% internal heating.The temperatures at the top and bottom boundaries are set to T a and T b . The surface and the base of the model are shear-stress free.

FIGURE 3 . 2 :

 32 FIGURE 3.2: Probability density functions of background surface velocity (left) and background surface heat flux (right) compared to their Gaussian approximations (in light gray)

  FIGURE 3.4: Maps of the two-point spatial correlation of temperature. The arrow indicates the location of the anchor point. The maps represents the correlation between the temperature at anchor points and the temperature points on the rest of the domain. Considering the symmetries of correlations described on equations 3.28 and 3.29, only the correlations ∀φ ′ ∈ [0, π], r ′ ∈ [0, 1], Cov(T (0, r), T (φ ′ , r ′ )) for r = r b (A), r = (r a + r b )/2 (B) and r = r a (C) are shown.

Fig. 3 .

 3 Fig. 3.6 shows examples of the final forecast state (second column) for evolutions with different parameters, after 300 My of data assimilation: t10γ10 (31 observation times, model 1), t50γ10 (7 observation times, model 1), t10γ10Q40 (31 observation times, model 2) and t10γ10Ra7 (31 observation times, model 3). The two first cases are done using the same model parameters, so we display the data assimilation results of the same evolution for better comparison. The true temperature fields for each case are displayed on the first column. The local error on the third column of Fig. 3.6 is the absolute value of the error at each coordinate (φ m , r n )

FIGURE 3 . 6 :

 36 FIGURE 3.6: Comparison of the results of data assimilation after 300 My for four cases. First row: data assimilation with a ∆t = 10 My (experiment t10γ10), second row: data assimilation with a ∆t = 50 My (experiment t50γ10), third row: data assimilation with a ∆t = 10 My (experiment t10γ10Q40) and fourth row: data assimilation with a ∆t = 10 My (experiment t10γ10Ra7). The first column displays the true temperature field for each experiment, the second column shows the analyzed temperature field, and the third column is the associated error multiplied by 10 and represented using the same color scale than the field itself, in order to ease visual inspection.

  FIGURE 3.8: Evolution of the average relative error on the forecast temperature field with ∆t = 10 My and γ = 5, 10, 30 and 50%. Note that the scale on the y-axis is logarithmic.

  FIGURE 3.9: Evolution of the average relative error on the forecast temperature field with ∆t = 10 My and γ = 10% for the model 1, model 2 and model 3. Note that the scale on the y-axis is logarithmic.

Fig. 3 .

 3 Fig. 3.8 shows the evolutions of errors for data assimilations with analyses every 10 My and different errors on observations: t10γ5, t10γ10, t10γ30 and t10γ50. For levels of noise in observation ranging from 5 to 30%, the error reduction phase lasts more than 300 My. Even when observations are noised as much as 50%, data assimilation keeps errors to values lower than 10% on the estimate of the temperature field over the first 100 My. Fig.3.8 shows that low levels of noise in observations (5 to 10%) give results of similar quality.

  FIGURE 3.10: Comparison between true, observed and analyzed surface heat flux (left) and surface velocities (right) after 100 My of assimilation with the parameters t10γ50

  FIGURE 3.11: Evolution of the cumulative mean innovation as a function of the number of analyses. On the left: for model 1 with γ = 10% and ∆t = 1, 5, 10, 15, 20 and 50 My; on the right: for Model 1 with ∆t = 10 My and γ = 5, 10, 30 and 50%

  premier ordre la dynamique interne et la tectonique de surface du système croûtemanteau. Nous montrons ici que les données tectoniques (comme la cinématique de surface et la répartition de l'âge du fond marin) et les modèles de convection du manteau avec comportement en plaques peuvent en principe être combinés pour reconstruire la circulation mantellique. Nous présentons une méthode d 'assimilation de données séquentielle, basée sur un schéma sous-optimal dérivé du filtre de Kalman, où les vitesses de surface et les cartes d'âges des fonds océaniques ne sont pas utilisées comme conditions limites pour l'écoulement mais comme données à assimiler.Deux étapes (une prévision suivie d'une analyse) sont répétées séquentiellement pour intégrer des données observées à différents moments. Chaque fois que des observations sont disponibles, l'étape d'analyse détermine l'état le plus probable du manteau à ce moment, en tenant compte d'une estimation préalable (fournie par la prévision) et des nouvelles observations disponibles, en utilisant la meilleure estimation linéaire non biaisée. Entre deux temps d'observation, l'évolution du manteau est régie par le modèle direct de convection mantellique. Pour évaluer son efficacité, cette méthode est appliquée à des cas synthétiques en 2-D (géométrie d'anneau sphérique). Nous comparons les évolutions de référence aux estimations obtenues par assimilation de données. Deux paramètres contrôlent le comportement du schéma: le temps entre deux analyses, et l'amplitude du bruit dans les observations synthétiques. Notre technique se révèle efficace pour reconstruire des évolutions de champ de température, à condition que le temps entre deux analyses soit de 10 Myr. Si l'amplitude de l'erreur a priori sur les observations est importante (30%), notre méthode fournit une meilleure estimation de la tectonique de surface que les observations, en tirant profit de l'information contenue dans la physique de la convection.

  .,Fournier, A., & Coltice, N., in prep. Ensemble data assimilation for mantle circulation, Nonlinear Processes in Geophysics4.1 IntroductionMantle circulation models are estimates of mantle flow history. They combine two sources of information: observations on the dynamics or 3D structure of the Earth's mantle and a numerical model of mantle convection. In their effort to reconcile both observations and our physical understanding of mantle dynamics, they serve a wide variety of purposes and disciplines.[START_REF] Hager | Kinematic models of large-scale flow in the earth's mantle[END_REF] originally built instantaneous mantle circulation models to understand the effect of plates on large-scale mantle flow. Since then, they have been used, among other applications, to understand the dynamics and evolution of the deep earth mantle structures

  FIGURE 4.1: Time evolution of the errors on the estimated temperature field obtained from data asssimilations with the same 150 Myr dataset, but different assimilation parameters. The size of the ensemble is A) N = 96, B) N = 288 and C) N = 768. The assimilations are computed for any combination of data assimilation parameters: γ + = 1.25, ℓ v = 0.3, 0.5, 0.7 and 1 and ℓ h = π/10, π/8, π/6, π/4 and π/2. The black line represent the evolution of the error for the best assimilation A)N = 96, ℓ z = 0.5, ℓ h = π/6 and γ + = 1.25 B)N = 288, ℓ z = 0.7, ℓ h = π/10 and γ + = 1.25 and C)N = 768, ℓ z = 0.5, ℓ h = π/4 and γ + = 1.25. The grey area is delimited by the maximum and minimum values of errors at each time, for all data assimilations.

  .31) At each analysis time, we compute the Euclidean norm of the instantaneous innovation d k and the Euclidean norm of the cumulative mean innovation d k d k = d k and d computing these norms, we normalize the part of the innovation corresponding to surface heat flux and velocities by their respective root mean square values. The norm of the instantaneous innovation d k measures the distance between the forecast data and the observation, and therefore gives indications on the success or failure of the assimilation. The evolution of the cumulative mean of the innovation d k allows us to check the consistency of the data assimilation algorithm at first order. Indeed, the derivation of the EnKF equations assumes that the error on observations y o and the error on the forecast data Hx f are unbiased. Such hypotheses imply that the statistically expected value of d is zero, which means that the norm of the cumulative innovation should converge to zero as the number of analyses increases.

Figure 4 .

 4 Figure 4.2 shows the evolution of d k and d k as a function of the number of analyses for data assimilations with different sizes of ensemble and their respective optimum vertical and horizontal correlation lengths. The norm of the instantaneous innovation (Figure 4.2B) first decreases during the first 8 analyses, i.e. 70 My, and then oscillates for the rest of the assimilation. We observe the same behavior for the evolution of the error on temperature ǫ e T (k): the norm of the instantaneous innovation for ensemble sizes of N = 288 and 768 are very similar, and lower than that obtained with N = 96.

Figure 4 .

 4 Figure 4.2A shows the cumulative innovation constantly decreasing throughout the assimilation, with comparable values for N = 288 and N = 768, and slightly higher values for N = 96.

  FIGURE 4.2: Evolution of A) the cumulative mean innovation and B) the norm of the instantaneous innovation, as a function of number of analyses and for different ensemble sizes. For each size of the ensemble, the evolutions correspond to the best combinations of correlation length parameters: N = 96, ℓ z = 0.5, ℓ h = π/6 and γ + = 1.25 ; N = 288, ℓ z = 0.7, ℓ h = π/10 and γ + = 1.25 and N = 768, ℓ z = 0.5, ℓ h = π/4 and γ + = 1.25.

  and 4.5. Figure4.4 compares the true temperature field evolution with the analyzed temperature field of method 1 and of the Ensemble Kalman filter with N = 288, ℓ v = 0.7 and ℓ h = π/10. The sudden increase in the error of the estimated temperature field for method 1 seen on Figure4.3B happens after around 80 Myr of assimilation, when the direction of bending at the bottom of the domain changes for the downwelling on the left side (see Figure4.4, second row). The analyzed temperature field of Method 1 does not predict this change of direction (see Figure4.4, first row), while the analyzed temperature field of the ENKF predicts it (see Figure4.4, third row). Method 1 computes only the evolution of the best estimate of the system. The computation of only one estimate ignores that, in this case, a slight perturbation of the estimated state could lead to a totally different dynamics. On the contrary, the EnKF method computes the evolution of an ensemble of perturbed solutions and thus takes into account the nonlinearity of the solution, at least for the forecast stage.

  FIGURE 4.4: Comparison of temperature field evolutions for evolution B. The first row depicts the evolution of the analyzed temperature field with method 1. The second row is the true evolution of the temperature field. The third row is the evolution of the analyzed temperature field with ensemble Kalman filter, N = 288, ℓ v = 0.7 and ℓ h = π/10.

  FIGURE 4.6: Value of the mean analyzed error for the same evolution and different vertical and horizontal correlation lengths. A) for 96 ensemble members, B) for 288 ensemble members, C) 768 ensemble members. The dashed lines delimit the zones for which errors are less than ǫ a T min (N ) + 0.2%.

FIGURE 4 . 7 :

 47 FIGURE 4.7: Value of A) mean analyzed error, B) mean ensemble spread, C) average forecast error on data, D) norm of cumulative innovation after K analyses for N = 288, and different vertical and horizontal correlation lengths. The dashed line delimits the zone for which errors are less than ǫ a T min (288) + 0.2%.

Figure 4 .

 4 Figure 4.8 shows the final state of the assimilation after 150 My for the evolution A of Figure 4.3. We selected 3 assimilations: EnKF96, an Ensemble Kalman Filter with N = 96, ℓ v = 0.5 and ℓ h = π/6 (first row), EnKF288 an EnKF with N = 288, ℓ v = 0.7 and ℓ h =π/10 (second row) and the assimilation with method 1 (third row). We do not show the Ensemble Kalman Filter with 768 members since the resulting temperature field is almost indistinguishable from that of EnKF288. The first column represents the true temperature field, which is the same for all assimilations. The second column is the analyzed temperature field, i.e. the average of the temperature fields of the analyzed

  FIGURE 4.8: Comparison of estimated states after 150 Myr for the evolution A of Figure 4.3. First row: Ensemble Kalman Filter with N = 96, ℓ v = 0.5 and ℓ h = π/6; second row: N = 288, ℓ v = 0.7 and ℓ h = π/10; third row for method 1 (third row). The first column represents the true temperature field at 150 Myr, the second column is the analyzed temperature field, the third is the absolute error on temperature value and the fourth is the estimated error on the analyzed field (spread of the ensemble). On the true temperature field of EnKF288 we framed the location of the subduction (A), plume (B), ridge initiation (C) and stable ridge (D) studied in Figure 4.9.

  FIGURE 4.9: Evolution of the temperature value during assimilation at the surface (first row), in the mid mantle (second row) and at the bottom of the domain (third row), for a ridge initiation (first column), a subduction (second column), an upwelling (third column) and a stable ridge (fourth column).

acknowledgements

  The research leading to these results has received funding from the European Research Council within the framework of the SP2-Ideas Program ERC-2013-CoG, under ERC grant agreement 617588. Calculations were performed using HPC resources from GENCI-IDRIS (grant 2014-047243). The contribution of Alexandre Fournier is IPGP contribution number XXX.

  plate tectonics reconstructions on one hand and mantle convection codes on the other hand. The continuous progress of mantle convection codes in modelling Earth's inner dynamics and surface tectonics suggest a more ambitious, long-term goal: the raw observation on surface tectonics (the mapping of magnetic anomalies on the ocean floor, the observation and dating of tectonic features, paleomagnetic data...) could be integrated directly into mantle convection models, without using the theoretical framework of plate tectonics.The revolution of Earth sciences which took place in the 60's (see Section 2.2.2), with the emergence of plate tectonics theory, was not merely a change in paradigm. It required the Earth science community to gather and connect observations coming from various disciplines to build a unifying global theory of Earth's surface motions. In doing so, it set a new standard for the functioning of solid Earth Science. It brought the Earth science community to consider geological processes as a multidisciplinary and multiscale problem: to understand the global history of the solid Earth, we need to navigate between micro, macro, local, regional and global scales and reconcile observations of different nature, coming from various disciplines, in a consistent framework. Plate tec-

  Au cours de cette thèse, nous avons successivement proposé deux algorithmes pour reconstruire la circulation mantellique en combinant les observations sur la tectonique de surface avec un modèle de convection mantellique, en utilisant des méthodes d'assimilation de données séquentielle.Nous avons d'abord développé un schéma d'assimilation séquentielle sous-optimal basé sur le filtre de Kalman (Chapitre 3). Ce schéma est efficace pour reconstituer le champ de température d'un système convectif avec une tectonique de type plaque à sa surface sur plusieurs 100 Myr. Nous avons testé la robustesse de la méthode en effectuant des expériences synthétiques en géométrie sphérique 2D. Nous avons obtenu des reconstitutions fidèles, même pour l'emplacement des panaches à la base des modèles, après des périodes d'assimilation d'au moins 300 Myr, à condition que le temps entre les analyses soit inférieur à 15 Myr et que l'amplitude du bruit dans les observations soit inférieure à 30%.Ces résultats encourageants nous ont poussés à explorer davantage le problème de l'assimilation des données pour la reconstitution de la circulation mantellique à partir de la tectonique de surface (Chapitre 4). Le filtre de Kalman d'Ensemble semblait être une méthode de choix compte tenu de sa capacité à traiter des problèmes aux dimensions élevées et des modèles non linéaires. Il s'est en effet révélé plus robuste et en moyenne plus précis que l'ancienne méthode, tout en fournissant en plus des cartes détaillées des incertitudes sur l'évolution estimée.Les expériences de synthèse que nous avons réalisées dans les deux études ne sont nullement exhaustives et de nombreux autres tests restent à faire pour évaluer soigneusement les méthodes d'assimilation des données dans le contexte de la reconstitution de la circulation mantellique. Parmi eux, nous identifions la nécessité:• d'expériences jumelles supplémentaires utilisant des paramètres plus réalistes:du côté du modèle dynamique, cela signifie, entre autres, des tests en géométrie sphérique 3D, l'introduction de continents, un nombre de Rayleigh plus élevé, un saut de viscosité au milieu du manteau, mais également prendre en compte les aspects chimiques de la convection mantellique (voir par exemple Ricard, 2015, section 6 pour un aperçu des spécificités de la convection mantellique).du côté des données observées, cela signifie utiliser une matrice de covariance des erreurs sur les données observées qui prend en compte les corrélations entre les erreurs sur des points adjacents et la corrélation entre les erreurs sur les vitesses de surface observées et les flux de chaleur de surface observés. Pour les modèles développant une convection à petite échelle, Il sera également nécessaire d'utiliser directement les cartes de l'âge des fonds marins au lieu des flux de chaleur de surface comme données observées.• la conception d'«expériences soeurs», pour évaluer l'effet des erreurs dans les modèles de convection mantellique. Le principe de ces expériences soeurs serait de calculer une évolution de référence avec un ensemble de paramètres physiques donnés, puis d'appliquer l'algorithme d'assimilation de données en utilisant un ensemble différent de paramètres.Bien sûr, le test ultime pour un schéma d'assimilation de données est son application aux données réelles. Dans le chapitre 4, nous avons développé un code d'assimilation de données pour la reconstitution de la circulation mantellique en combinant un code d'assimilation de données d'ensemble parallélisé(PDAF Nerger et al., 2005;[START_REF] Nerger | Software for ensemble-based data assimilation systems-implementation strategies and scalability[END_REF] avec le code parallèle de convection STAGYY(voir Tackley, 2008). La grande scalabilité du code obtenu et les capacités de modélisation de STAGYY (voir section 1.5) rendent techniquement possible l'application de l'assimilation de données à des données terrestres. Le calcul d'une reconstitution de la circulation mantellique exige en outre:• une série chronologique de la disposition des plaques et de leurs mouvements.Cette série chronologique est facilement accessible grâce à la politique d'accès libre des équipes travaillant sur les reconstitutions de tectonique des plaques (par exemple[START_REF] Müller | Ocean basin evolution and global-scale plate reorganization events since pangea breakup[END_REF].• une évaluation complète des erreurs sur les reconstitutions de la tectonique des plaques, en tenant compte des erreurs dans les mouvements de plaques relatifs et absolus, et en cohérence avec le formalisme utilisé dans l'assimilation des données. Cette tâche ne peut être réalisée que par une étroite collaboration avec les équipes travaillant sur les reconstitutions de la tectonique des plaques. Comme on l'a souligné dans la section 2.2.2, des outils ont été développés pour estimer les incertitudes sur les mouvements de plaques relatifs et absolus (voir, par exemple, Wessel & Müller, 2015 pour une revue récente de ces techniques). Par conséquent, une estimation au premier ordre des erreurs sur les reconstitutions de la tectonique des plaques est possible, même si cette estimation pourrait ne pas refléter les incertitudes véritables sur les reconstitutions, en particulier pour les temps très anciens (voir discussion à la fin de la section 2.2.2).• des choix dans la paramétrisation pour le modèle de convection mantellique. Ce choix de paramètres constituera un compromis entre d'une part l'obtention d'une dynamique interne et d'une tectonique de surface aussi proche que possible de la Terre et d'autre part les limites en termes de puissance de calcul disponible.

  d'assimilation de données pour la circulation mantellique.Nous avons commencé cette thèse en soulignant la coexistence de deux théories différentes qui sont utilisées pour comprendre la dynamique du système croûte-manteau terrestre: la tectonique des plaques et la convection mantellique. Les développements successifs d'algorithmes d'assimilation de données que nous avons proposés sont une tentative de combiner les produits de ces deux théories: les reconstitutions de la tectonique des plaques d'une part et les codes de convection mantellique d'autre part. Le développement continu des codes de convection mantellique dans la modélisation de la dynamique interne de la Terre et de la tectonique de surface suggère un objectif à long terme plus ambitieux: les observations brutes sur la tectonique de surface (cartes des anomalies magnétiques des fonds océaniques, données paléomagnétiques...) pourraient être intégrées directement dans les modèles de convection mantellique, sans utiliser le cadre théorique de la tectonique des plaques.La révolution des sciences de la Terre qui a eu lieu dans les années 60 (voir la section 2.2.2), avec l'émergence de la théorie de la tectonique des plaques, n'était pas simplement un changement de paradigme. Il a fallu que la communauté des sciences de la Terre rassemble et relie les observations provenant de diverses disciplines pour construire une théorie globale expliquant les mouvements de la surface terrestre. Ce faisant, elle a établi une nouvelle norme pour le fonctionnement de sciences de la Terre Solide. Elle a amené la communauté des sciences de la terre à considérer les processus géologiques comme un problème multidisciplinaire et multi-échelles: pour comprendre l'histoire globale de la Terre solide, nous devons naviguer entre les échelles microscopique, macroscopique, locale, régionale et globale et concilier des observations de nature différente, provenant de diverses disciplines, dans un cadre cohérent. La théorie de la tectonique des plaques a fourni un tel cadre global et s'est révélée juste dans sa description de la dynamique de surface, au premier ordre. Cependant, la tectonique des plaques est une théorie cinématique: elle décrit les mouvements à la surface de la Terre, mais pas leur cause. En outre, il s'agit d'un modèle idéalisé des mouvements de surface: les plaques ne sont pas rigoureusement rigides, et les limites des plaques

  However, the mantle model they consider differs from the model used in this dissertation on two major points.First, in these models, the viscosity depends neither on temperature, nor on the state of stress, but only on depth. This has important consequences for the development of data assimilation techniques. As discussed in Section 1.4, the temperature and stress dependence of the viscosity is necessary to obtain mantle convection models which display plate-like tectonics at their surface. If we use a mantle model which does not have this type of viscosity law, then it is not possible to integrate the data derived from

[START_REF] Bunge | Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography[END_REF]

writes the equivalent of the objective function of Equation

2

.20 in the continuous case, using the classical conservation equations governing mantle flow, as described in Chapter 1. plate reconstructions as is done in classical data assimilation schemes, i.e. by defining an operator H that links a state computed from the model to surface velocities. In this case, the only data which is usable as proper "observed data", as defined previously, is the temperature field at present, estimated from seismic tomography. The time series of observed data {y o 1 , y o 2 , ..., y o K } shrinks to {y o K }. The vector y o K contains the values of the discretized temperature field derived from seismic tomography. The associated

TABLE 3

 3

		.1: Values of the parameters of the forward models
	Symbol Meaning	model 1 model 2 model 3
	Ra	Rayleigh number	10 6	10 6	10 7
	H	Non-dimensional internal heating rate	20.5	13.5	20.5
	M	Resolution of the grid in the longitude direction	384	384	512
	N	Resolution of the grid in the radial direction	48	48	64
	r a	Radius of the top of the domain	2.2	2.2	2.2
	r b	Radius of the bottom of the domain	1.2	1.2	1.2
	T a	Temperature at the top of the domain	0	0	0
	T b	Temperature at the bottom of the domain	0.9	1.3	1
	E A	Activation Energy	23.03	23.03	23.03
	T 1	Temperature at which µ T = 1	1	1	1
	β	Factor of viscosity reduction for partial melting	10	10	10
	T s 0	Solidus Temperature at r = r a	0.6	0.6	0.6
	∇ r T s	Radial gradient of the solidus temperature	2	2	2
	σ Y	Yield Stress	1.10 4	1.10 4	4.10 4
	∇ r σ Y	Radial gradient of the yield stress	2.10 5	2.10 5	2.10 5

  The research leading to these results has received funding from the European Research Council within the framework of the

SP2-Ideas Program ERC-2013-CoG, under ERC grant agreement 617588. Calculations were performed on LGLTPE Seisglob high-performance computing cluster and using HPC resources from GENCI-IDRIS (grant 2014-047243). Figures were generated using matplotlib library Hunter (2007). The contribution of Alexandre Fournier is IPGP contribution number 3675.

TABLE 4 .

 4 1: Values of the parameters of the forward model

	Symbol Meaning

TABLE 4 .

 4 2: Notations and dimensions of data assimilation variables compute two ensembles of N states {x f kn } n∈[1,N ] and {x a kn } n∈[1,N ] , such that their average equals x f k and x a k , respectively, and their respective sample covariance matrices approximate P f k and P a k . The ensemble of states {x f kn } n∈[1,N ] and {x a kn } n∈[1,N ] are stored in the matrices X f k and X a k , where the nth column is the state of the nth ensemble member x f kn and x a kn , respectively.

	Symbol Meaning

TABLE 4 .

 4 3: Notations and range of values tested for data assimilation parameters

	Symbol Meaning	value
	N	number of ensemble members	96 to 768
	K	number of observation times	16
	γ +	maximum inflation factor	1.25
	ℓ v	vertical correlation length	0.3 to 1
	ℓ h	horizontal correlation angle	π/10 to π/2

  in finite rotation data and produce time series of high resolution plate relative motions. More recently,[START_REF] Iaffaldano | Rapid plate motion variations through geological time: Observations serving geodynamic interpretation[END_REF] applied this technique to the relative motion of the pacific plate with North America for the last 75 Myr. The uncertainties on relative plates velocities ranges from 5 to 40% of the root mean square surface velocity. As we go further back in time, the quantification of relative plate motion uncertainties becomes hazardous: most of the seafloor created before 150 Myr has been destroyed by subduction. These plate tectonic reconstructions involve interpretation of different types of data, with a limited spatial coverage and relies heavily on human expertise. For these epochs, maintaining very high uncertainties on the regions where few data supports the reconstructions would be a solution.

[START_REF] Iaffaldano | Reconstructing plate-motion changes in the presence of finite-rotations noise[END_REF]

applied the trans-dimensional hierarchical Bayesian method to reduce noise

Abstract

With the progress of mantle convection modelling over the last decade, it now becomes possible to solve for the dynamics of the interior flow and the surface tectonics to first order. We show here that tectonic data (like surface kinematics and seafloor age distribution) and mantle convection models with plate-like behaviour can in principle be combined to reconstruct mantle convection. We present a sequential data assimilation method, based on suboptimal schemes derived from the Kalman filter, where surface velocities and seafloor age maps are not used as boundary conditions for the flow, but as data to assimilate. Two stages (a forecast followed by an analysis) are repeated sequentially to take into account data observed at different times. Whenever observations are available, an analysis infers the most probable state of the mantle at this time, considering a prior guess (supplied by the forecast) and the new observations at hand, using the classical best linear unbiased estimate. Between two observation times, the evolution of the mantle is governed by the forward model of mantle convection. This method is applied to synthetic 2-D spherical annulus mantle cases to evaluate its efficiency. We compare the reference evolutions to the estimations obtained by data assimilation. Two parameters control the behaviour of the scheme: the time between two analyses, and the amplitude of noise in the synthetic observations. Our technique proves to be efficient in retrieving temperature field evolutions provided the time between two analyses is 10 Myr. If the amplitude of the a priori error on the observations is large (30 per cent), our method provides a better estimate of surface tectonics than the observations, taking advantage of the information within the physics of convection.

Abstract

Recent advances in mantle convection modelling led to the release of a new generation of convection codes, able to generate self-consistently plate-like tectonics at their surface. Those models physically link mantle dynamics to surface tectonics. Combined with plate tectonic reconstructions, they have the potential to produce a new generation of mantle circulation models that use data assimilation methods and where uncertainties on plate tectonic reconstructions are taken into account. We recently provided a proof of this concept by applying a suboptimal Kalman Filter to the reconstruction of mantle circulation [START_REF] Bocher | A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics[END_REF]. Here, we propose to go one step further and apply the Ensemble Kalman Filter (EnKF) to this problem. The EnKF is a sequential Monte Carlo method particularly adapted to solve high dimensional data assimilation problems with nonlinear dynamics. We tested the EnKF using synthetic observations consisting of surface velocity and heat flow measurements, on a 2D-spherical annulus model and compared it with the method developed previously. The EnKF performs on average better and is more stable than the former method. Less than 300 ensemble members are sufficient to reconstruct an evolution. We use covariance adaptive inflation and localization to correct for sampling errors. We show that the EnKF results are robust over a wide range of covariance localization parameters. The reconstruction is associated with an estimation of the error, and provides valuable information on where the reconstruction is to be trusted or not.