# Thermal evolution of the Moon

Matthieu Laneuville

► To cite this version:

Matthieu Laneuville. Thermal evolution of the Moon. Sciences of the Universe [physics]. Institut de Physique du Globe de Paris, 2013. English. NNT: 2013GLOB1201. tel-04021007

HAL Id: tel-04021007

https://hal.science/tel-04021007

Submitted on 9 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Institut de Physique du Globe de Paris

Ecole doctorale des Sciences de la Terre

Thèse de Doctorat

pour le titre de

Docteur en Science

de l'Institut de Physique du Globe de Paris

Specialité : GEOPHYSIQUE

Défendu par

Matthieu LANEUVILLE

Évolution thermique de la Lune

Équipe PLANÉTOLOGIE ET SCIENCES SPATIALES,
défendu le 5 Décembre, 2013.

Jury :

Directeur: Mark WIECZOREK - IPG (Paris)

Co-directeur: Doris BREUER - DLR (Berlin)

Rapporteur : Tilman SPOHN - DLR (Berlin)

Rapporteur : Stephane LABROSSE - ENS (Lyon)

Examinateur : Julien AUBERT - IPG (Paris)

Examinateur : Claude JAUPART - IPG (Paris)

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

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

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

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

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>
Remerciements

Sans ordre particulier, je voudrais remercier Mark Wieczorek et Doris Breuer pour leur suivi continu ces 3 dernières années. J'aimerais aussi remercier les équipes de sciences planétaires de l'IPG et du DLR pour leur soutien et les nombreuses discussions passionnantes. J'aimerais particulièrement remercier Ina et Nicola pour leur patience à m'enseigner les trucs et astuces de l'utilisation de GAIA. Des commentaires de Julien Aubert, Chloé Michaut, Guillaume Morard et Francis Nimmo ont beaucoup aidé à l'amélioration de ce manuscrit. J'aimerais aussi bien sûr remercier mes parents et amis pour avoir supporté mes bavardages lunaires incessants ces dernières années!

Ce travail a été financé par une bourse de la région Ile-de-France. Les simulations numériques ont été faites sur le cluster TITANE installé au Centre de Calcul Recherche et Technologie, JUROPA installé au Jülich Supercomputing Center et HLRN, le centre pour calcul intensif de l'Allemagne du nord.

In no particular order, I would like to thank Mark Wieczorek and Doris Breuer for their continuous guidance during these 3 years. I also would like to thank both IPGP and DLR's planetary sciences teams for the support and many enlightening discussions. I especially would like to thank Ina and Nicola for their invaluable help using GAIA. Comments by Julien Aubert, Chloé Michaut, Guillaume Morard and Francis Nimmo on this manuscript were also greatly appreciated. I of course also would like to thank my parents and friends for supporting my constant jabbering about the Moon these last few years!

This work was funded by the Region Ile-de-France. Numerical simulations were made on the TITANE cluster from the Centre de Calcul Recherche et Technologie, JUROPA at the Jülich Supercomputing Center and HLRN, north Germany's supercomputing center.
<figure><img src="image_8.png" /></figure>
# Contents

<table><tr><td>0</td><td>Résumé de la problématique et résultats principaux</td><td>1</td></tr><tr><td>1</td><td>The Moon – Context and overview</td><td>11</td></tr><tr><td>1.1</td><td>Comparative planetology</td><td>11</td></tr><tr><td>1.2</td><td>The Moon, what do we know?</td><td>13</td></tr><tr><td>1.2.1</td><td>Missions overview</td><td>13</td></tr><tr><td>1.2.2</td><td>Formation scenarios</td><td>15</td></tr><tr><td>1.2.3</td><td>The magma ocean hypothesis</td><td>15</td></tr><tr><td>1.2.4</td><td>Volcanic activity</td><td>16</td></tr><tr><td>1.2.5</td><td>Present day state</td><td>17</td></tr><tr><td>1.3</td><td>Summary and outstanding questions</td><td>18</td></tr><tr><td>I</td><td>Lunar asymmetry and mantle evolution</td><td>21</td></tr><tr><td>2</td><td>The lunar asymmetry</td><td>23</td></tr><tr><td>2.1</td><td>A brief history</td><td>23</td></tr><tr><td>2.2</td><td>Previous models</td><td>25</td></tr><tr><td>2.3</td><td>Motivation</td><td>26</td></tr><tr><td>3</td><td>Lunar thermo-chemical evolution</td><td>27</td></tr><tr><td>3.1</td><td>Convective thermochemical evolution modeling</td><td>27</td></tr><tr><td>3.2</td><td>Initial conditions</td><td>29</td></tr><tr><td>3.2.1</td><td>Heat sources</td><td>29</td></tr><tr><td>3.2.2</td><td>Initial temperature profile</td><td>33</td></tr><tr><td>4</td><td>Thermo-chemical evolution results</td><td>35</td></tr><tr><td>4.1</td><td>Thermal evolution</td><td>35</td></tr><tr><td>4.2</td><td>Magmatism</td><td>41</td></tr><tr><td>4.3</td><td>Predicted gravity field</td><td>45</td></tr><tr><td>4.4</td><td>Dynamo generation</td><td>50</td></tr><tr><td>4.5</td><td>Discussion</td><td>54</td></tr><tr><td>4.6</td><td>Conclusions & perspectives</td><td>55</td></tr><tr><td>II</td><td>Core evolution and paleomagnetism</td><td>57</td></tr><tr><td>5</td><td>The lunar magnetic history</td><td>59</td></tr><tr><td>5.1</td><td>(Paleo)Magnetic field measurements</td><td>59</td></tr><tr><td>5.2</td><td>Past and current models of magnetic history</td><td>62</td></tr><tr><td>5.3</td><td>Existence of a core</td><td>64</td></tr></table>
<table><tbody><tr><td>5.4</td><td>Our approach</td><td>66</td></tr><tr><td>6</td><td>Evolution models of the lunar dynamo</td><td>69</td></tr><tr><td>6.1</td><td>Lunar thermal evolution</td><td>69</td></tr><tr><td>6.2</td><td>Energy budget of the core</td><td>70</td></tr><tr><td>6.3</td><td>Entropy budget</td><td>73</td></tr><tr><td>6.4</td><td>Magnetic field scaling</td><td>75</td></tr><tr><td>6.5</td><td>Model parameters</td><td>76</td></tr><tr><td>7</td><td>Lunar core evolution</td><td>81</td></tr><tr><td>7.1</td><td>Foreword on 2D/3D geometry</td><td>82</td></tr><tr><td>7.2</td><td>Core-mantle evolution coupling</td><td>82</td></tr><tr><td>7.2.1</td><td>Influence on the core</td><td>82</td></tr><tr><td>7.2.2</td><td>Influence on the mantle</td><td>83</td></tr><tr><td>7.3</td><td>Core evolution : Typical case</td><td>84</td></tr><tr><td>7.4</td><td>Parameter space exploration</td><td>85</td></tr><tr><td>7.4.1</td><td>Influence of core composition</td><td>85</td></tr><tr><td>7.4.2</td><td>Influence of initial core thermal state</td><td>86</td></tr><tr><td>7.4.3</td><td>Other influences</td><td>87</td></tr><tr><td>7.5</td><td>Discussion</td><td>92</td></tr><tr><td>7.6</td><td>Conclusion</td><td>93</td></tr><tr><td>III</td><td>Summary and perspectives</td><td>95</td></tr><tr><td>8</td><td>The 21st century Moon</td><td>97</td></tr><tr><td>8.1</td><td>Summary</td><td>97</td></tr><tr><td>8.1.1</td><td>Mantle thermal state</td><td>98</td></tr><tr><td>8.1.2</td><td>Core properties</td><td>98</td></tr><tr><td>8.2</td><td>Perspectives</td><td>99</td></tr><tr><td>IV</td><td>Appendices and bibliography</td><td>101</td></tr><tr><td>A</td><td>Gravity field calculation</td><td>103</td></tr><tr><td>A.1</td><td>Gravity field</td><td>103</td></tr><tr><td>A.2</td><td>Calculation of static and dynamic topography</td><td>104</td></tr><tr><td>B</td><td>Core budget calculations</td><td>107</td></tr><tr><td>B.1</td><td>Entropy budget</td><td>107</td></tr><tr><td>B.2</td><td>Simplified equations</td><td>108</td></tr><tr><td></td><td>List of notations</td><td>111</td></tr><tr><td></td><td>Bibliography</td><td>115</td></tr></tbody></table>
Résumé de la problématique et
résultats principaux

Des avancées colossales concernant la compréhension de la structure interne de la Terre, de sa composition ainsi que des processus géologiques sous-jacents ont été faites durant le siècle passé. Néanmoins, il existe des aspects qui resteront toujours difficiles voire impossibles à étudier du point de vue terrestre, notamment à cause de la tectonique des plaques qui assure le recyclage de la croûte terrestre. De part son activité géologique limitée et sa quasi-absence d'érosion, la surface lunaire garde la trace des 4.5 milliards d'années de son évolution, ce qui permet de mieux comprendre non seulement son évolution initiale mais aussi les processus affectant le Système Solaire peu de temps après sa formation. La compréhension de l'histoire lunaire est donc une clé pour comprendre l'évolution de la Terre.

La Lune a longtemps été considérée comme un objet n'ayant eu aucune activité géologique majeure depuis sa formation. Ce fût donc une surprise quand, à partir de la fin des années 1960, les missions Apollo et Luna ont apporté la preuve que des épisodes de volcanisme intense ont eu lieu sur la face visible ainsi que de la présence d'un champ magnétique à la surface. Un élément frappant est que le volcanisme semble avoir été principalement présent sur la face visible, où environ un tiers de la surface est recouverte de laves, contre moins d'un pour-cent de la face cachée (Figure 1). Ces éléments ont entrainé un regain d'intérêt pour la Lune, notamment car, étant donné sa petite taille, l'existence d'une activité géologique durable a de fortes chances d'être liée directement à sa formation.

Il existe 4 différents scénarios pour expliquer la formation de la Lune: formation par fission, capture, co-accrétion ou impact géant. La formation suite à l'impact d'un corps de la taille de Mars sur la proto-Terre est le modèle privilégié aujourd'hui mais reste le sujet de nombreuses discussions. Les rapports isotopiques de certains éléments sont en effet identiques entre la Terre et la Lune, ce qui suggère que la Lune a été formée avec la même matière que la Terre. De plus, le moment cinétique actuel du système Terre-Lune pose de grosses contraintes sur la géométrie d'un tel impact. Celle-ci est encore débattue, mais une constante dans ces modèles est que la Lune s'est accretée à partir du disque de débris en 100 à 1000 ans, ce qui est assez rapide pour faire fondre une couche épaisse à la surface. Ce processus, connu sous le nom d'océan de magma, a une très grande influence sur l'évolution thermo-chimique de la Lune.

En effet, la cristallisation fractionnée d'un océan de magma entraîne sa stratification. La cristallisation débute au fond de l'océan et le liquide restant s'enrichit en éléments incompatibles. Les éléments lourds ont tendance à sédimenter alors que les éléments légers vont à la surface. Une particularité de l'océan de magma sur la Lune
<figure><img src="image_9.png" /></figure>

Figure 1: Mosaïque d'images de Lunar Reconnaissance Orbiter (WAC) présentant la face visible (gauche) et cachée (droite) de la lune. NASA/GSFC/Arizona State University.

est que le plagioclase, qui est plus léger que le liquide environnant quand il cristallise, est stable à une profondeur assez importante pour lui permettre d'être extrait du liquide efficacement et être transporté à la surface. C'est la raison pour laquelle la surface lunaire est formée d'anorthosite quasiment pure, une roche riche en anorthite, pole calcique du plagioclase. Sous cette croûte doivent donc se trouver les éléments les plus incompatibles de l'océan de magma. Cet aspect est d'autant plus intéressant que les éléments radioactifs, principal moteur de l'évolution thermique des planètes, font partie de ces éléments incompatibles.

Cependant, le modèle prédit un enrichissement global d'éléments radioactifs sous la croûte, alors qu'en surface, ils apparaissent concentrés dans une région située sur la face visible appelée le Procellarum KREEP Terrane (qui tire son nom de sa concentration en potassium (K), terres rares (Rare Earth Elements) et phosphate (P)). De plus sous la condition d'une cristallisation fractionnée, la stratification engendrée possède un profil de densité instable et le devenir des différentes couches juste après la solidification de cet océan magmatique est incertain. Enfin, il n'est pas entièrement clair si l'asymétrie du volcanisme est dû à un processus de transport asymétrique du magma produit dans le manteau, ou si la production magmatique elle-même était asymétrique. Depuis les résultats de la mission Lunar Prospector en 1998, la distribution en éléments radioactifs à la surface est connue globalement et est étroitement corrélée au volcanisme extrusif, comme montré dans la Figure 2.

Un lien semble donc exister entre ces deux observations, mais leur relation exacte est incertaine. La première partie de ce travail de thèse a donc été de modéliser l'évolution thermo-chimique du manteau lunaire en testant différentes hypothèses concernant la distribution en éléments radioactifs et le profil de température initial
<figure><img src="image_10.png" /></figure>

Figure 2: (haut) Concentration de thorium à la surface par Lunar Prospector Lawrence et al. (2003). Le contour du Procellarum KREEP Terrane est définie par la ligne blanche et correspond au contour de 4 ppm de thorium. Les étoiles représentent les sites d'atterrisage des missions Apollo, et les cercles ceux des missions Luna. Les étoiles jaunes sont l'emplacement des deux mesures de flux de chaleur des missions Apollo. Les ellipses sur la face cachée correspondent à l'anomalie interne et limite structurale externe du basin South Pole-Aitken. La face visible est à gauche, la face cachée à droite et la projection est Lambert azimuthale à surface constante (image modifié depuis Mimoun et al. (2011)). (bas) Carte des ages des mare basaltiques de Hiesinger et al. (2003) affichée sur la carte des mare basaltiques de l'USGS. Ces cartes suggèrent que le volcanisme et la distribution des sources de chaleur dans la croûte sont liés.

pour tester si la compréhension que nous avons aujourd'hui de la distribution des éléments radioactifs est cohérente avec l'évolution du manteau. Une description plus approfondie des connaissances que nous avons sur la Lune ainsi que des questions ouvertes est faite dans le Chapitre 1. Je vais ici présenter maintenant les principaux résultats de la première partie de notre étude.
## Evolution asymétrique du manteau

Les observations à distance ne permettent pas de savoir si l'enrichissement en éléments radioactifs est un phénomène de surface, ou si l'enrichissement se poursuit en profondeur sur des dizaines de kilomètres. Un indice qui laisse à penser que l'enrichissement n'est pas que superficiel mais concerne probablement toute l'épaisseur de la croûte est que ces éléments sont présents principalement dans les roches qui forment la croûte, et non dans les laves. Il est possible d'estimer le volume de cet enrichissement correspondant à une cristallisation symétrique de l'océan de magma, puis de supposer que ce même volume se trouvait sur la face visible au début de nos simulations. Plusieurs processus ont été avancés pour expliquer cette localisation, et pour cette étude nous supposons que l'un de ces processus fonctionne, sans pour autant le modéliser. Une description complète de ces problématiques est faite lors de l'introduction du Chapitre 2.

Nous avons utilisé un programme de modélisation 3D de la convection manuelle pour simuler l'évolution du manteau lunaire avec différentes hypothèses concernant son état initial. Nous considérons deux distributions extrêmes pour comprendre leurs implications: (1) les sources de chaleur sont localisées sous la croûte (et ont donc la même densité que le manteau) et (2) les sources de chaleur sont mélangées à la croûte (et ont donc la même densité que la croûte). Un schéma de notre modèle est visible dans la Figure 3. Dans les deux cas, l'enrichissement est uniquement dans une région circulaire de 40 degrés de rayon centrée sur le Procolarum KREEP Terrane. La concentration en éléments radioactifs dans le manteau a été choisi pour que la composition lunaire totale soit chondritique, ou encore enrichie d'un facteur 2 en éléments radioactifs par rapport à la composition chondritique, ce qui correspond aux deux modèles principaux généralement considérés dans la littérature. Le modèle numérique que nous utilisons ainsi que les hypothèses sur son état initial est décrit en détail dans le Chapitre 3.

Le résultat principal de ces simulations est que, lorsque l'enrichissement est situé sous la croûte, la majorité du volcanisme a lieu sur la face visible. Il n'est donc pas nécessaire de trouver un mécanisme de transport de magma asymétrique pour expliquer l'asymétrie dans les volumes d'extrusions. Le temps nécessaire à faire fondre le manteau supérieur correspond aussi au délai d'emplacement des coulées de lave et nous sommes donc confiant qu'un tel modèle permet d'expliquer les observations concernant le volcanisme. Ces simulations ont aussi permis de se rendre compte que l'anomalie thermique induite par la concentration en éléments radioactifs restait présente jusqu'à aujourd'hui. Cela a une influence sur la densité du manteau, et donc sur les mesures gravitaires. L'évolution de la température du manteau lunaire en fonction du temps d'après nos simulations est représentée en la Figure 4.

Nos modèles d'évolution thermo-chimique nous permettent de prédire la distribution 3D de la densité en fonction du temps et nous avons donc pu modéliser les différentes contributions chimiques et thermiques résultant de la dynamique 3D sur l'anomalie gravitaire prédite aujourd'hui. Une première contribution est liée à la dilatation thermique, qui change la densité du manteau mais aussi le rayon de la
<figure><img src="image_11.png" /></figure>

Figure 3: Diagramme schématique du modèle de convection thermo-chimique (pas à l'échelle). L'équivalent de 10 km de basaltes 'KREEP' est placé sous 40 km de croûte (bleu), dans les 20 km les plus profonds de la croûte (haché), ou redistribué dans toute l'épaisseur de la croûte (orange). L'étendue latérale de la région est 80° de diamètre, et la résolution radiale est 20 km. La température du noyau $T_c$ varie avec le temps, mais la cristallisation potentielle d'une graine n'est pas pris en compte.

Lune. La deuxième est liée à l'aspect chimique des simulations avec une réduction de la densité dans le manteau par fusion partielle, et l'emplacement de laves denses à la surface. En fonction des hypothèses faites sur l'état de compensation de la surface, différents résultats sont trouvés, mais une conséquence importante reste: une anomalie négative est prédite au centre du Procellarum KREEP Terrane aujourd'hui (cf. Figure 5). Ce résultat a des implications immédiates pour l'interprétation des données de sondage électromagnétique, la modélisation des déformations de marées ainsi que l'interprétation des données sismiques Apollo. Les détails des calculs que nous avons menés et les discussions concernant leurs applications sont décrits dans le Chapitre 4 ainsi que dans l'annexe A.

## Paléo-magnétisme et thermodynamique du noyau

Le deuxième aspect troublant de l'évolution lunaire est l'existence d'un champ magnétique rémanent à la surface de plus d'une dizaine de microteslas et qui semble avoir été enregistré il y a 3 à 4 milliards d'années. A cause de la faible taille de la Lune, il est difficile d'expliquer la présence d'un champ si important et de longue durée dans le passé (pour comparaison, le champ terrestre est de l'ordre de 50 mi-
<figure><img src="image_12.png" /><figcaption>Figure 4: Coupe du manteau lunaire montrant l'évolution de la température pendant 4.5 milliards d'années dans le cas où un profil de température est intermédiaire et où les sources de chaleur sont situées sous la croûte (modèle 'T-0LB'). Le temps avant le présent est indiqué sur chacune des graphiques, en milliards d'années. Le cercle noir est le noyau et le blanc correspond aux régions partiellement fondues. Les lignes de courant sont représentés par les lignes pointillées.</figcaption></figure>

croteslas). Une synthèse des mesures paléomagnétiques est disponible dans la Figure 6. L'existence d'une dynamo thermique dans le noyau lunaire a été proposé, mais elle ne peut expliquer la durée du champ observé. D'autres processus, plus exotiques, ont été proposés depuis, notamment la possibilité d'un champ magnétique induit par la rotation différentielle entre manteau et noyau. Des variations de ce modèle existent, mais le processus même de génération d'une dynamo par friction
<figure><img src="image_13.png" /></figure>

Figure 5: Champ gravitaire de la Lune. Contribution (a) de l'expansivité thermique, (b) du changement de densité du manteau dû à la fusion partielle, (c) de la topographie induite, et (d) des coulées volcaniques non compensées avec une épaisseur maximale de 5 km. Pour ces contributions, le degré 0 de l'expansion en harmoniques sphériques a été supprimé, mais le degré 1 est conservé. L'anomalie gravitaire prédite (e) dans le cas d'une lithosphère flexible et les basaltes compensés, et (f) dans le cas d'une lithosphère rigide avec des basaltes non compensés et pas de topographie induite. Pour les deux dernières images, le degré 1 de l'expansion en harmonique sphérique (qui correspond à un déplacement du centre de masse) a été supprimé. Toutes les images utilisent une projection de Mollweide à surface constante centré sur le Procellarum KREEP Terrane. Le cercle noir correspond à l'étendue de cette région. Les résultats sont ceux de notre modèle préféré (modèle 'T-0LB').
à l'interface noyau-manteau n'ayant pas été étudié en profondeur, il est incertain si
les lois d'échelle habituelles sont applicables à l'estimation du champ induit.

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

Figure 6: Synthèse des études paléomagnétiques anciennes (rouge) et récentes (vert). Les flèches vertes signalent que les valeurs données sont les valeurs du champ magnétique minimal. La zone orange est la zone initialement considérée comme correspondant à l'existence d'un champ magnétique actif. La ligne pointillée bleu est une value typique du champ magnétique terrestre.

Dans cette partie de la thèse nous avons étudié la possibilité de la présence d'une dynamo chimique dans le noyau lunaire liée à la croissance de la graine. Ce processus permet la génération d'un champ magnétique pendant plus longtemps qu'une dynamo thermique standard, car même quand le flux de chaleur est inférieur au flux critique pour la dynamo thermique, la flottabilité engendrée à l'interface graine-noyau due au chauffage par libération de chaleur latente et à la concentration d'éléments légers dans le noyau externe est suffisante pour entraîner une convection globale dans le noyau. Une description détaillée de la problématique est disponible dans l'introduction du Chapitre 5.

Pour modéliser ce phénomène, nous avons couplé les simulations effectuées dans la première partie de la thèse à un modèle unidimensionel du noyau. Pour contraindre l'évolution du noyau, nous avons établi son budget en énergie et entropie. Le budget d'énergie associé à une équation d'état définissant le comportement du fluide constituant le noyau externe en fonction de la température nous a permis de calculer la taille de la graine en fonction du temps. Une fois celle-ci obtenu, il est possible d'étudier les différentes sources et puits d'entropie pour estimer le taux de dissipation maximal d'énergie magnétique et donc l'existence d'un champ magnétique
stable sur des temps géologiques. La description complète du modèle est disponible
au Chapitre 6 et dans l'annexe ??.

Nos simulations montrent que la croissance de la graine coïncide avec la durée estimée du paléo-magnétisme. Cependant, une fois ce processus engagé il devient difficile de le stopper, et nos simulations prédissent la présence d'un champ magnétique aujourd'hui. Un exemple d'évolution typique se trouve en Figure 7. Ce désaccord évident avec les données peut être expliqué de plusieurs façons. Une possibilité est que l'enrichissement continuel du noyau externe en éléments légers le fasse passer dans le régime de "pluie de fer" – c'est à dire que la graine cristallise à l'interface noyau-manteau à la place de l'interface graine-noyau. Un tel processus est moins efficace pour produire un écoulement de grande échelle dans le noyau puisque la libération de chaleur latente et l'enrichissement en éléments légers s'effectuent principalement à l'interface noyau-manteau, et n'induisent alors plus de convection et il est probable que la dynamo cesse à ce point.

<figure><img src="image_15.png" /><figcaption>Figure 7: Résultats obtenus pour un cas avec 6wt.% de souffre et ΔT = 300 K; modèle dT300_X6 dans la Table 7.1. (a) Flux de chaleur à l'interface noyau-manteau, (b) température à cette même interface, (c) taille de la graine et (d) champ magnétique estimé à la surface en fonction du temps.</figcaption></figure>

L'autre possibilité est liée aux lois d'échelles utilisées. Dès que la graine at-
teint une fraction non négligeable du noyau en taille, la couche disponible pour
la convection dans le noyau externe est très faible et il est possible que le champ
magnétique généré par l'écoulement dans cette géométrie, si écoulement il y a, ne
corresponde pas aux paramètres habituels et donc que nos prédictions deviennent
alors fausses. En conclusion nous avons montré que, soit le noyau lunaire est par-
tiellement cristallisé, auquel cas la génération d'un champ magnétique est inévitable
et une solution doit être trouvé pour expliquer son absence aujourd'hui, soit le noyau
lunaire est encore totalement liquide, ce qui implique une température relativement
élevée et/ou une concentration en éléments légers importante. Dans les deux cas,
des mesures complémentaires contraignant notamment la taille de la graine, cou-
plées à nos prédictions permettront une meilleure compréhension de l'évolution du
noyau lunaire. Les résultats sont présentés avec plus de détails et leurs implications
sont discutées dans le Chapitre 7.
# The Moon – Context and overview

## Contents

<table><tr><td>1.1 Comparative planetology</td><td>11</td></tr><tr><td>1.2 The Moon, what do we know?</td><td>13</td></tr><tr><td>1.2.1 Missions overview</td><td>13</td></tr><tr><td>1.2.2 Formation scenarios</td><td>15</td></tr><tr><td>1.2.3 The magma ocean hypothesis</td><td>15</td></tr><tr><td>1.2.4 Volcanic activity</td><td>16</td></tr><tr><td>1.2.5 Present day state</td><td>17</td></tr><tr><td>1.3 Summary and outstanding questions</td><td>18</td></tr></table>

## 1.1 Comparative planetology

Modern planetary sciences date back to the works of Copernic, Kepler, Galileo and Newton during the 16th and 17th centuries. They understood the relative place of the Earth and Sun in our Solar System as well as the different forces necessary to explain the observed planetary motion. But the advent of the space age also revolutionised our understanding of the Solar System. The first major space age success can be defined as the soviet mission Luna 3 in 1959, which provided the first images of the lunar farside. In the following decades, every major body of the Solar System has been visited providing a wealth of new information. Even more recently, the discovery of more than 700 confirmed exoplanets provides context for our Solar System and leads to a renewed interest in comparative planetology.

The current view of our Solar System is now that of a chaotic system with complex resonances and orbital variations, to be compared with the old view of an Earth-centric system with a clock-like perfection to it. The early Solar System formation may have included giant planets migration (Lin et al., 1996) and an impact cataclysm in the Earth-Moon region known as the late heavy bombardment (e.g., Wetherill, 1975). Delivery of water to the Earth is also believed to have occurred through impact with planetary embryos formed in the outer parts of the Solar System (Morbidelli et al., 2000). Orbit resonances such as those of Io and Europa around Jupiter play a major role in their evolution (Hussmann and Spohn, 2004), while the effect of tides on early lunar evolution is still under investigation
(e.g., *Meyer et al.*, 2010). Thus the Solar System always has to be seen as a whole and complex system in order to have a chance of catching subtle influences.

The idea behind comparative planetology is that despite their similarities, every
planet is different and allows us to test the limits of our theories. The theory of plate
tectonics for example, was developed on Earth to explain the observed continental
drift. It is a very efficient way to cool planetary mantles, and its major influence on
atmospheric composition over geologic timescale links it tightly with habitability. It
has therefore been used since on Mars to try and understand conditions under which
an early magnetic field could be produced (*Nimmo and Stevenson*, 2000; *Breuer*,
2003). Before being mapped by Magellan, it was thought that Venus could have
had plate tectonics to explain what looked like oceanic spreading centers (*Head
and Crumpler*, 1989). It was later shown that plate tectonics were actually absent
and this has been attributed to the relative dryness of the lithosphere (*Nimmo and
McKenzie*, 1998). Since the dawn of the exoplanets era, the likelihood of plate tec-
tonics on super-Earth has also been investigated in order to assess the potential for
habitability in distant solar systems (e.g., *Korenaga*, 2010). Therefore even if the
physics behind the onset of plate tectonics is still under active research, it is the
knowledge we have from multiple planets that forces us to setup a consistent frame-
work valid throughout a large range of conditions (e.g., *Tackley*, 2000; *Bercovici*,
2003).

After accretion, most planetary bodies differentiate, i.e., the interior is segregated
by density, forming different layers known as core, mantle, crust and sometimes an
atmosphere. This process depends mostly on the size and bulk content of the planet,
while gravity and magnetic field influence the stability of the atmosphere. Once this
separation is made, exchanges are possible only through specific phenomena. Crust
is formed by mantle melting and can be recycled back in the mantle mostly through
plate tectonics. Outgassing of volatiles is often associated with crustal processes thus
providing a way to exchange material between mantle and atmosphere. Estimates
made on one of this reservoir can often be coupled with the others in order to test
its viability.

The diversity of planetary environments is clearly observable when considering
coupling between these different reservoirs. The thermal regime of planetary man-
tles and the composition of their potential atmosphere are tightly linked. The range
of these interactions was studied recently on several terrestrial bodies. A positive
feedback between surface temperature and mantle outgassing was shown to be im-
portant on Venus (*Phillips et al.*, 2001), while mantle outgassing is a possible source
of methane for Titan's atmosphere (*Tobie et al.*, 2006). On Mars, (*Jakosky and
Phillips*, 2001) discussed the complex interactions between the different planetary
reservoirs, highlighting their effect on habitability.

The wide range of (or lack of) planetary dynamos is also a prominent example of
the interest of comparative planetology. The terrestrial planets have very different
magnetic histories (*Stevenson et al.*, 1983). The Earth and Mercury possess a present
day magnetic field, while the Moon and Mars present only remanent crustal magnetic
fields and Venus lacks a global detectable field. Understanding the mantle and core
conditions required to explain this range of behavior forces us to provide a more consistent picture of solar system formation and subsequent evolution (Stevenson, 2003).

In addition, some bodies have a better record of older times due to a limited amount of erosion. The Moon is one of these. The very old lunar surface is a record of what other planetary bodies might have looked like 4.5 Ga ago. Most of the lunar crust is primordial whereas the Earth or Venus are comparatively much younger. There are processes that are hard to understand on the Earth because surface erosion limit our direct investigation of the past. The cratering record for example is clearly visible on the lunar surface, but only major or recent features remain on Earth.

The formation of the Moon being tightly linked to the formation of the Earth, a better understanding of the Moon's evolution will also provide clues about Earth's evolution. And understanding Earth's early evolution better, is a key to understanding the origin of life. Especially, any dramatic event that occurred on the Moon, such as the putative late heavy bombardment, is bound to have influenced the Earth as well. Study of the Earth and Moon as a pair thus provides synergies and insight both on lunar and earth science.

In this project we therefore focused on our Moon because of its tight link with the evolution of the Earth and the renewed international interest of the last decade. This interest and the wealth of data from older missions makes lunar science less speculative than other planets. I will now present what the current understanding of lunar history is, and what outstanding questions remain.

## 1.2 The Moon, what do we know?

The most important lunar nearside features, the maria, take their names from early astronomers who thought they were actual seas. Their nature remained debated until the Apollo/Luna era, where the picture of a cold, undifferentiated body was definitely discarded. Figure (1.1) shows these surface features with a near and farside mosaic from the Lunar Reconnaissance Orbiter mission's wide-angle camera. The next few subsections provide a very quick history of major lunar science topics since the Apollo/Luna, with emphasis on the latest developments.

### 1.2.1 Missions overview

Since the end of the Apollo/Luna era with Luna 24 in 1976, only a few lunar missions were conducted. The first object to orbit the Moon again was the Japanese Hiten spacecraft in 1990, which was only a technology demonstration probe and demonstrated the aerobreaking technique. The 1990s was also the decade of Clementine (1994) and Lunar Prospector (1998), both orbiters part of the NASA program. Clementine's main achievement was complete mapping of the surface composition and altimetry. Radar data also suggested the presence of ice around the pole for the first time. Lunar Prospector allowed global mapping of the crust's major heat
<figure><img src="image_16.png" /><figcaption>Figure 1.1: Lunar Reconnaissance Orbiter WAC mosaic of the lunar near- and farside. NASA/GSFC/Arizona State University</figcaption></figure>

sources and several other elements which lead to the definition of the nearside region now known as the Procellarum KREEP Terrane due to its high enrichment in potassium, rare earth elements and phosphorous.

After the SMART-1 mission from ESA in 2003 designed to demonstrate new technologic developments, a few missions from different space agencies were launched, proving the renewed international interest towards the Moon. At one month interval both SELENE and Chang'e 1 were launched in 2007, from the Japanese and Chinese space agency, respectively. Chandrayaan 1 was launched a year later by the Indian space agency. All three missions had a similar profile with 1 to 2 years orbit concluded by a controlled impact. The instruments ranged from spectrometers to altimeters, radars, magnetometers which all contributed to improvements of several global fields. The Lunar Reconnaissance orbiter, launched by NASA in 2009, is mostly famous for its high resolution imaging of the surface and investigation of the lunar surface water content.

Since 2010, three new missions have already been successfully launched to the Moon. Chang'e 2, by the Chinese space agency, orbited the Moon for almost a year before continuing its mission at Earth L2 Lagrange point and to the asteroid 4179 Toutatis. The Gravity Recovery And Interior Laboratory (GRAIL) mission by NASA was launched in 2011 with the objective of determining precise crustal properties through high degree gravity field modeling, and was a success. Finally, the Lunar Atmosphere and Dust Environment Explorer (LADEE) was launched on September 6 this year and is on its way to the Moon. Follow-up missions by the Chinese, Indian, Japanese and possibly Russian space agencies are in development but the current budget situation may make some of the plans hard to achieve.
### 1.2.2 Formation scenarios

The four types of lunar origin models are fission, capture, co-accretion and giant impact. To differentiate between these models, three main constraints are used. First, the angular momentum of the current Earth-Moon system has to be explained, second, the isotopic resemblance between the bulk silicate Earth and the Moon suggests a common origin, and finally the relative depletion of the Moon in iron suggests an anomalous history.

The fission model assumes that early Earth was spinning fast enough so that surface material became unstable and was ejected into orbit (Darwin, 1880). The capture model assumes that the Moon was formed elsewhere in the Solar System and was later captured into a stable orbit around the Earth (Wood and Mitler, 1974). In the co-accretion model, the Moon accretes from a circum-terrestrial disk that was present during the complete Earth accretion process. The giant impact scenario, assumes that this disk formed from the giant impact of the Earth with a large projectile (e.g., Stevenson, 1987).

Giant impact is the current accepted model, but some discrepancies with observations still have to be cleared. For example, a recent work by van Westrenen and de Meijer (2012) stresses that the bulk lunar chemistry, in which most isotopic system are identical to that of the Earth, is inconsistent with dynamical model predictions, which find that the Moon will re-accrete mostly from the impactor's mantle and thus has no reason to have similar isotopic composition than the Earth. Though recent dynamical models, using a previously neglected Sun-Moon resonance, now achieve better post-impact predictions (Canup, 2012; Cuk and Stewart, 2012). Regardless of the exact scenario, it takes about 100 to 1000 years to re-accrete the Moon, which is fast enough to produce a thick layer of molten material on top of the Moon (e.g., Mizutani et al., 1972). This very particular initial condition, the so-called lunar magma ocean, will have a tremendous influence on the initial conditions of thermal evolution models.

### 1.2.3 The magma ocean hypothesis

Historically, the concept of the lunar magma ocean was introduced before the giant impact scenario was the preferred formation hypothesis. As early as 1970, it was recognised that the lunar crust was mostly anorthositic and therefore probably the product of global differentiation (Warren, 1985). There is little doubt that cooling a lunar magma ocean would crystallise buoyant plagioclase that would rise to the surface and form the observed ancient anorthositic crust. This is typical of lunar lithology because it requires a large body of differentiating magma that is not frequent on Earth. An example of such a formation on Earth is the Sept Iles layered intrusion in Quebec, Canada (Namur et al., 2011). In addition, a plagioclase crust is less likely on large bodies because it would crystallise too close to the surface to efficiently separate. Since the early proposal, other arguments in favor of the lunar magma ocean appeared, the main ones being maybe the strong homogeneity
of refractory elements ratio among lunar samples and the complementary europium
anomaly between mare basalts and ferroan anorthosites. This strongly suggests that
they originate from the same parent magma.

However the depth of this magma ocean remains uncertain today. Numerous
studies have been made since the 1970s and the possible range now appears to be
200 to 1000 km (e.g., *Solomon and Longhi, 1977*; *Elkins-Tanton et al., 2011*). The
two main constraints on that size are the thickness of the anorthositic crust, which
provide an estimate of the bulk aluminum content and the total contraction and
expansion during the subsequent thermal evolution. Such radius changes being very
small on the Moon, it limits the amount of cooling and warming that can occur
during history. Consequent cooling of the magma ocean (contraction) and warming
of the primordial lower mantle (expansion) have often been used to justify the small
radius changes and would favor magma oceans less than about 600 km deep (*Kirk
and Stevenson, 1989*).

But obviously these constraint are only as good as the models that are used to
test them. Assuming fractional or equilibrium crystallisation of the magma ocean
changes the composition profile, the bulk initial composition also has a large in-
fluence on what chemical phases can crystallise and therefore, on the composition
profile. In addition, few actual crystallisation experiments have been made (*Elardo
et al.*, 2011; *Tronche and van Westrenen*, 2011), and thus constraints on these as-
sumptions are scarce. However, it is known that fractional crystallisation of a magma
ocean will produce a gravitationally unstable layering. The latest layers to crystallise
below the crust are also the densest and are enriched in heat sources. The behavior
of that layering is crucial for the understanding of the subsequent thermochemical
evolution. Thus coupling such results to dynamical modeling should shed light on
the lunar bulk composition and its initial thermochemical structure.

1.2.4 Volcanic activity

The first global view of lunar magmatism dates back to 1959 and the first images
we obtained from the farside by the soviet mission Luna 3. It appeared that surface
volcanism was very heterogeneously distributed on the Moon, with about one third
of the nearside covered by maria against less than a percent on the far side (see
*Head, 1976*, for a review). Dating of lunar samples, in combination with crater
counting techniques, showed that mare basalts on the nearside range in age from
about 4 to 1 Ga (*Hiesinger et al.*, 2010) and those on the far side from 3 to 2.5 Ga
(*Haruyama et al.*, 2009).

The origin of this volcanic asymmetry has long been disputed and there are
now strong arguments in favor of a genetic link between the heat source enrichment
also observed in that region and volcanism *Jolliff et al.* (2000). However the role
of that enrichment remains unclear. *Zhong et al.* (2000) argue that it is merely a
by-product of the volcanic phase 4 Ga ago and that a large part of the KREEP layer
was initially on top of the core-mantle boundary. On the contrary, *Wieczorek and
Phillips* (2000) argue that it is the strong heat sources enrichment in the Procellarum
KREEP Terrane that is the cause of the 4 Ga volcanic phase. Arguments have been recently put forward for both types of models *Zhang et al.* (2013); *Laneuville et al.* (2013) and it is still unclear what scenario best fits observations.

A large dataset on surface composition has been gathered in the last decade and an important result is that the far side crust, which is mainly anorthositic, appears to be made of rocks that crystallised from less evolved magma than the nearside crust *Ohtake et al.* (2012). This suggest that the lunar dichotomy is directly linked to magma ocean crystallisation and was thus already emplaced at the beginning of lunar thermal models.

1.2.5 Present day state

The current view about the lunar interior comes from the compilation of a few decades worth of data. Initial studies from Apollo seismic data suggested a discontinuity at about 560 km depth *Nakamura et al.* (1982). This has been confirmed (e.g., *Khan et al.*, 2000) but as a result of the large uncertainties in the arrival times of the deep moonquakes, the transition could not be unambiguously identified. Conceptually, it could well be the limit between the upper mantle, which was part of the magma ocean and an undifferentiated part of the mantle below. But other explanations exist, and crystallisation of the magma ocean itself could produce such a discontinuity (*Wieczorek et al.*, 2006). Another result from seismology is the evidence for a region of partial melt in the lower mantle. The lack of shear waves arrivals at some stations, added to the very effective dissipation required in the lower mantle to explain global dissipation rates both suggest, though not unambiguously, the presence of partial melt below a depth of about 1000 km *Nakamura et al.* (1973).

Recent re-interpretation of Apollo seismic data allowed a better constraint of the deep lunar interior. Using different approaches, *Weber et al.* (2011) and *Garcia et al.* (2011) both provided an estimate of the core size. *Weber et al.* (2011) used a new processing method to conclude that the core is about 330 km in radius and is about 60% liquid by volume (i.e., possess an inner core of about 240 km in radius). *Garcia et al.* (2011) on the other hand found that the core is 380 ± 40 km and does not observe an inner core. The discrepancies between the two conclusions in such recent papers clearly show that new data is required if we wish to decipher the structure of the lunar interior and thus, its early evolution.

Despite the discrepancies on the deep interior, new structural data has been obtained in the past few years. The GRAIL mission has been used to construct a very high spherical harmonic degree gravity field of the Moon (Zuber et al., 2013). This huge improvement with respect to previous data from Clementine, Lunar Prospector and SELENE helped resolve small surface features and thus improve crustal thickness models. For example, *Wieczorek et al.* (2013) found that the crust is considerably less dense than previously thought, probably due to an average porosity of 12%. Another crucial result is that the crust is now thought to be between 34 and 43 km on average, against close to 50 km before GRAIL. This reduces the amount of aluminum in the crust and thus that the Moon is not enriched in refractory elements
with respect to the Earth, which has strong implications on formation scenarios and thermal evolution models.

Finally, several unambiguous evidence for water on the lunar surface have been found in the past years. The idea of cold traps in permanently shadowed regions is not new (e.g., *Watson et al.*, 1961), but its discovery is relatively recent. Clementine bistatic radar experiment provided the first hints (*Nozette et al.*, 1996), which have since been confirmed multiple times by radar measurements (*Lawrence*, 2011) and spectrometry of the ejecta cloud of LCROSS (*Colaprete et al.*, 2010). However these correspond to volatiles from extra-lunar origin. Another line of observations have shown that the Moon also has a non-negligible water content in its mantle. This was observed first in Apollo samples, using new experimental techniques (*Saal et al.*, 2008) and later using remote observation of the surface (*Klima et al.*, 2013). All this suggest that the Moon may not be as dry as previously thought, and this should be integrated in lunar mantle simulations as the influence on viscosity will be non-negligible.

## 1.3 Summary and outstanding questions

Tremendous progress has been made in the last 50 years following the Apollo/Luna era on our comprehension of the Moon. From a cold, primordial body from unknown origin the Moon has now become an initially molten body formed by re-accretion of the debris from a giant impact involving the Earth. The many missions of the last decades and number of ground breaking articles in all aspects of lunar science is proof that our nearest neighbour still has a lot to teach us. They are naturally all linked and can be traced back to the formation scenario and early evolution of the Earth-Moon system. The main questions are:

### What is the bulk composition of the Moon?

Linked to the problematics of crustal structure and composition and the initial magma ocean is the question of the bulk composition of the Moon. Long thought to be enriched in refractory elements, new developments now tend to argue for a chondritic Moon. This has a large influence on the formation scenario and early history of the Moon.

### What is the origin of the lunar asymmetry?

As will be presented in the next chapter, the Moon possess a global asymmetry both in structure and composition. The very high correlation between heat sources and volcanism on about one third of the lunar surface has been the source of many discussions since it was globally determined by the mission Lunar Prospector in 1998. As gamma-ray spectrometry samples only the very surface of the planet, the total enrichment and its links to volcanism are still debated.
**What are the characteristics of the putative magma ocean? How does it
influence subsequent evolution?**

Although its existence is widely accepted, the characteristics of the putative lunar magma ocean are highly debated. The predictions made by models of the magma ocean evolution are crucial for the subsequent thermochemical evolution of the Moon, yet no consensus has been reached to date, with thicknesses ranging from a few hundred kilometres to whole Moon melting.

**How large is the lunar core? What is its composition and present thermal
state?**

As previously seen, large discrepancies still exist between core size estimates and the potential existence of an inner core. The amount and type of alloying element in the core has yet to be determined and could be constrained by using data from seismology and the mantle thermal state.

**How can lunar paleomagnetism be explained?**

As I will discuss in the second part of this thesis, the history of lunar magnetism is still unclear. The duration and strength of the active magnetic field era is not well constrained, and to date, no model has been able to match all observations.

**What is the deep structure of the lunar mantle, especially, what is the
cause of the low velocity zone?**

Analysis of Apollo missions' seismic data demonstrated the existence of a zone of low velocity at the base of the mantle, which is currently understood as a zone of partial melting. However that constraint is still vague, as its global extent is unsure.

This manuscript deals with a part of the story. It is organised in two main
parts, the first one deals with the mantle evolution and the influence of heat sources
distribution on present day predictions. The second part focuses on the core evolu-
tion in order to explain paleomagnetic measurements. Each part is organised in an
introductory section focused on that theme, then the model is described and results
are presented and discussed.
<figure><img src="image_17.png" /></figure>
# Part I

## Lunar asymmetry and mantle evolution
<figure><img src="image_18.png" /></figure>
# The lunar asymmetry

## Contents

<table><tr><td>2.1 A brief history</td><td>23</td></tr><tr><td>2.2 Previous models</td><td>25</td></tr><tr><td>2.3 Motivation</td><td>26</td></tr></table>

## 2.1 A brief history

The asymmetric distribution of lunar volcanism became apparent when the farside of the Moon was imaged for the first time by the soviet mission Luna 3 in 1959. It was surprising to realize that almost no basaltic lava flows were present on the farside and that most of the lunar maria were located on the nearside. Whereas about one third of the nearside hemisphere of the Moon has been resurfaced by basaltic lava flows, only about 1% of the farside is covered, principally within a few large impact basins such as South Pole Aitken, Moscovience, and Apollo. Dating of the lunar samples, in combination with crater counting techniques, showed that mare basalts on the nearside range in age from about 4 to 1 Ga (Hiesinger et al., 2003) and those on the farside from about 3 to 2.5 Ga (Haruyama et al., 2009). Taking into account the likelihood that older volcanism is now buried either below the maria or beneath the ejecta blankets of large ancient impact basins, the total duration of volcanism is greater than 3 billion years (Antonenko et al., 1995).

Analyses of Apollo γ-ray measurements and the thorium content of lunar samples suggested that heat producing elements were localized on the nearside of the Moon (Metzger et al., 1977; Haskin, 1998), but it was not until the Lunar Prospector mission in 1998 that the global distribution of heat producing elements was mapped (Lawrence et al., 1998). As shown in Figure (2.1), a high concentration of incompatible elements, in particular thorium and uranium, was observed on the nearside in a region strongly correlated with the maria (Lawrence et al., 1998), suggesting a genetic link between the two. Since thorium and uranium are usually the main heat sources that affect a planet's thermal evolution, their asymmetric distribution most likely lead to an asymmetric magmatic and geologic evolution (Wieczorek and Phillips, 2000).

Jolliff et al. (2000) defined the region of high thorium concentration on the nearside as the Procellarum KREEP Terrane (PKT), where the acronym KREEP stands
<figure><img src="image_19.png" /></figure>

Figure 2.1: (top) Lunar Prospector surface thorium concentrations of *Lawrence et al.* (2003). The PKT is outlined in white and is defined by the 4 ppm thorium contour. Stars are the Apollo landing sites, and the circles are the Luna sample return sites. The yellow stars are the locations of the two Apollo heat flow measurements. The farside ellipses are the inner compositional anomaly and outer structural rim of the South Pole-Aitken basin. The nearside is on the left, the farside on the right and the projection is Lambert azimuthal equal area (image modified from *Mimoun et al.* (2011)). (bottom) Mare basalt ages map of the nearside using data from *Hiesinger et al.* (2003) and the mare basalt map from the USGS. These maps suggest that volcanism and crustal heat sources are genetically related.

for potassium, Rare Earth Elements and Phosphorus, which occupies about 17% of
the lunar surface, or one third of the nearside hemisphere. Non-mare materials exca-
vated from beneath mare basalts show high thorium concentrations as well, implying
that KREEP is not solely enriched in the near-surface lava flows, but also in the un-
derlying crust (Jolliff *et al.*, 2000). Furthermore, the PKT contains both highlands
and volcanic flows, both of which can have high thorium concentrations. The total
extent of the concentration of heat sources in the crust is somewhat uncertain, as
it is not simple to determine the composition of the deep crust. Nevertheless, the
lack of KREEP signatures in the material ejected from large impact basins, such as Crisium, Humboldtianum, Moscovience, and Orientale, imply that KREEP in the crust is localized to the nearside hemisphere (see, e.g. *Warren*, 2001). Furthermore, a recent study on the viscoelastic deformation of lunar impact basins also suggests that the farside crust is depleted in radioactive elements with respect to the crust within the PKT (*Kamata et al.*, 2013).

It is widely accepted that the Moon formed hot and that part of its outer layers was molten to form what is known as the lunar magma ocean. As the magma ocean crystallized, the remaining liquid became progressively enriched in incompatible elements. When plagioclase began to crystallize, it was lighter than the surrounding liquid and rose to form the crust. The last part of the magma ocean to crystallize formed a highly evolved layer between the crust and mantle composed of materials enriched in KREEP, and also in iron-rich minerals such as ilmenite. The KREEP layer was once thought to form a global layer that was about 2 km thick (*Warren and Wasson*, 1979), but some process appears to have concentrated this material within the present-day PKT, forming an equivalent layer that might be about 10 km thick.

## 2.2 Previous models

Three different classes of models have been proposed to explain this宽scale seg-
regation. *Wasson and Warren* (1980) first suggested that the crystallization of the
magma ocean was asymmetric. A locally thicker crust on the farside would give rise
to a thinner underlying magma ocean, and thus concentrate KREEP-rich materials
on the nearside, where the crust is thinnest. This model requires an initially thicker
crust on the farside, possibly caused by a giant impact on the nearside (*Neumann
et al.*, 1996; *Byrne*, 2008), or a global convection pattern within the magma ocean
that could have transported crust preferentially to the farside (*Loper*, 2002).

Second, *Zhong et al.* (2000) assumed that a global layer of dense, late-stage ilmenite cumulates quickly sank to the core-mantle boundary, carrying along with it a large fraction of the KREEP layer. This layer then became hot and unstable and rose again on the nearside as a degree-1 upwelling, explaining both the present-day distribution of heat sources and the timing of volcanism. A more recent study by *Qin et al.* (2012) showed that the correlation between deep moonquakes and mare basalt could be the consequence of this ilmenite cumulate layer, if it was enriched in water. A related model, proposed by *Parmentier et al.* (2002), showed that the downwelling of a mixed ilmenite cumulate layer itself could follow a degree-1 pattern and concentrate KREEP beneath the PKT (though, see also *Elkins-Tanton et al.* (2002) for comments about that model).

Finally, impacts are often invoked to explain KREEP localization. *Ghods and Arkani-Hamed* (2007) showed that impacts would generate subsurface thermal anomalies that might redistributute KREEP-rich materials located beneath the crust. Depending on their size, the impacts would either mix this layer into the
underlying mantle, or simply concentrate KREEP at the edge of the basin, thus enhancing volcanic activity there. Another view involving impact is from *Jutzi and Asphaug* (2012) who claimed that a low-velocity impact from a 1200 km diameter companion to our moon might have displaced the global KREEP-rich layer to the nearside.

## 2.3 Motivation

Regardless of the origin of the PKT, its bare existence is expected to have completely changed the subsequent post-magma ocean thermal evolution of the Moon. The influence of the PKT on lunar evolution has already been investigated by *Wieczorek and Phillips* (2000) and *Hess and Parmentier* (2001), but these studies used highly simplified thermal models. *Wieczorek and Phillips* (2000) developed an axially symmetric 3D conduction model which showed that partial melting of the underlying mantle is an inevitable outcome of a thick KREEP layer on the nearside hemisphere, and that volcanism should span most of lunar history. This result was confirmed by *Hess and Parmentier* (2001) in a 1D thermal conduction study but they also noted that the wide, partially molten, region caused by that layer could form an impenetrable barrier to the eruption of mare basalts. They concluded that the hypothesis of a thickened KREEP layer below the PKT imposed strong constraints on the concentration of heat sources in the PKT and crustal thickness to remain consistent with both geological and petrological observations. Recently, *Grimm* (2013) re-analysed the results of *Wieczorek and Phillips* (2000) and showed that such models also predict large gravity or topography anomalies and electrical conductivity signatures that may be inconsistent with observations.

In this study, we use 3D thermochemical convection models to determine the consequences of KREEP localization in the PKT and compare these results with observations. Predictions are made that help determine if our current view of the PKT fits with the available data. In addition to previous models, we also calculate core temperatures and obtain core-mantle boundary heat flow estimates that can be used to test the origin and timing of a lunar dynamo. We also calculate surface gravity anomalies that can be compared with orbital observations. In Chapter 3, we present the details of the thermochemical convection model we use and describe the model's initial conditions. In Chapter 4, the results are first tested against previous models and observables, and then new consequences are presented.
# Lunar thermo-chemical evolution

## Contents

<table><tr><td>3.1</td><td>Convective thermochemical evolution modeling</td><td>27</td></tr><tr><td>3.2</td><td>Initial conditions</td><td>29</td></tr><tr><td>3.2.1</td><td>Heat sources</td><td>29</td></tr><tr><td>3.2.2</td><td>Initial temperature profile</td><td>33</td></tr></table>

The thermal evolution of the Moon is studied using a thermochemical convection model for a fluid with a temperature-dependent viscosity in a spherical shell. We consider both core cooling and time-dependent radioactive decay as heat sources. In this section, we present the equations to be solved followed by a description of the initial conditions that we used in our simulations.

## 3.1 Convective thermochemical evolution modeling

When subjected to stress, planetary mantles behave as fluids on geological timescales. Their evolution is described by the general Navier-Stokes equations supplemented by suitable assumptions. In particular, for a fluid dominated by diffusion creep (i.e. Newtonian), with an infinite Prandtl number and within the Boussinesq approximation, the non-dimensional equations for the conservation of mass, momentum and energy are written as follows

$$ \nabla \cdot \vec{u} = 0, \tag{3.1} $$

$$ -\nabla p + \nabla \cdot [\eta(\nabla \vec{u} + (\nabla \vec{u})^T)] + (T Ra_T + C Ra_C) \vec{e}_r = 0, \tag{3.2} $$

$$ \left( \frac{\partial}{\partial t} + \vec{u} \cdot \nabla \right) T = \nabla^2 T + Q - L_m \frac{\partial F}{\partial t}, \tag{3.3} $$

where $\vec{e}_r$ is the radial vector, $\vec{u}$ is the velocity, $p$ the dynamic pressure, $\eta$ the viscosity, $Ra_T$ and $Ra_C$ the thermal and compositional Rayleigh number, respectively, $t$ is time, $T$ the temperature, $Q$ the internal heat production rate, $L_m$ the latent heat of melting, $F$ the partial melt fraction and $C$ the depletion field, which tracks the chemical properties of the material and which corresponds to the cumulated melt fraction. We use the assumption of Newtonian flow (as has been used in other studies such as Konrad and Spohn (1997), Spohn et al. (2001), Stegman et al. (2003), Ziethe et al. (2009)), but we note, as previously noted by Christensen (1984), that the effect of non-Newtonian rheology could be important, especially at low stresses.
The diffusion creep viscosity $\eta$, normalized by $\eta_0$ can be written in non-dimensional form as follows (Roberts and Zhong, 2006)

$$ \eta(T) = \exp \left( \frac{E}{T + T_{surf}} - \frac{E}{T_0 + T_{surf}} \right), \qquad (3.4) $$

where $E$ is the activation energy, $T_{surf}$ the surface temperature and $T_0$ the reference temperature, at which $\eta = \eta_0$. Assuming dry olivine dominates mantle rheology, we use a reference viscosity of $10^{21}$ Pa s at 1600 K (Karato and Wu, 1993). Every quantity in the above four equations is non-dimensional and the nondimensionalization factors can be found in Table 3.1. We also assume that the viscosity dependence on composition and melt fraction is negligible to first order. For computational purposes, the viscosity is limited to a maximum value of $\eta = \eta_{max}$ where $\eta_{max}$ is taken such that the stagnant lid regime is reached, i.e. the viscosity is high enough so that the crust is no longer mobile. The Rayleigh numbers are given by

$$ Ra_T = \frac{\alpha_0 g \rho_0 \Delta T D^3}{\kappa_0 \eta_0}, \qquad (3.5) $$

$$ Ra_C = \frac{g \Delta \rho D^3}{\kappa_0 \eta_0}, \qquad (3.6) $$

where $\alpha_0$ is the reference thermal expansivity, $g$ the surface gravity acceleration, $\rho_0$ the reference density, $\Delta T$ the initial temperature drop across the mantle, $D$ the mantle thickness, $\kappa_0$ the reference thermal diffusivity, and $\Delta \rho$ is the density change upon 30% mantle depletion, corresponding to the change from peridotite to harzburgite (Scott and Stevenson, 1989). The two Rayleigh numbers arise in the adimensionalization process from the fact that the density is a function of temperature through thermal expansion and of the local depletion $C$, i.e.

$$ \rho = \rho_0 \left[ 1 - \alpha(T - T_0) - \frac{\Delta \rho C}{\rho_0 C_0} \right]. \qquad (3.7) $$

where $C_0$ is the depletion required to obtain harzburgite. The depletion field models the current chemical state of the material and is calculated using the following transport equation:

$$ \frac{\partial C}{\partial t} + \vec{u} \cdot \nabla C = \frac{\partial F}{\partial t}. \qquad (3.8) $$

As for the temperature-induced density anomalies, the Boussinesq approximation requires density variations due to depletion to be small with respect to the background density. We note also that gravity is assumed constant in the mantle, which will tend to overestimate buoyancy effects at the core-mantle boundary. Consumption of latent heat through melting is taken into account by a sink in the energy equation (Ita and King, 1994), where we assume a linear consumption of latent heat with melt fraction (see Eq. 3.3). Although melt transport is not directly taken into account, melt is assumed to leave the system instantaneously, thus no latent heat is released when the system cools down. The melt is assumed to rise vertically, and the total amount of melt generated below each element on the surface is tracked as
a function of time. For simplicity, we do not fractionate heat producing elements into the melt. For low degrees of melting (de Smet, 1999), the melt fraction can be written as

$$F = \frac{T - T_{sol}}{T_{liq} - T_{sol}}. \qquad (3.9)$$

This convective thermochemical modeling approach has been widely used Konrad and Spohn (1997); Spohn et al. (2001); Ziethe et al. (2009) and discussions about mantle composition are postponed to section 4. We assume here a peridotite composition with solidus and liquidus of KLB-1 peridotite from Hirschmann (2000):

$$T_{sol} = 1409 + 134.2P - 6.581P^2 + 0.1054P^3, \qquad (3.10)$$

$$T_{liq} = 2035 + 57.46P - 3.4872P^2 + 0.0769P^3, \qquad (3.11)$$

where P is the hydrostatic pressure in GPa. The effect of an insulating crust is also taken into account by assigning a lower thermal conductivity to the crust. Here, we assume a constant thickness of 40 km Khan and Mosegaard (2002); Lognonné et al. (2003); Wieczorek et al. (2013) and a thermal conductivity of $k_c = 1.5$ W m⁻¹ K⁻¹ in order to account for the presence of low-conductivity megaregolith. This value has been estimated from the case where ~5 km of megaregolith with $k_r = 0.3$ W m⁻¹ K⁻¹ is on top of 35 km of crust with $k_c = 3$ W m⁻¹ K⁻¹. All model parameters and their non-dimensionalization are described in Table 3.1.

We solve the set of equations (3.1)-(3.3) in a 3D spherical geometry using the finite-volume convection code Gaia, which has already been validated in other publications (e.g., Hüttig and Stemmer, 2008) and which can accurately handle localized viscosity variations of up to several orders of magnitude. Equation (3.8) is solved using a tracer approach using 20 particles per cell (e.g. Keken et al., 1997), and whose implementation in Gaia has been extensively tested by Plesa et al. (2013). This tracer approach essentially eliminates numerical diffusion and allows us to treat accurately sharp composition contrasts and to resolve sub-grid structures. Our simulations were carried out using a 20 km radial and a 60 km lateral resolution. The boundary conditions are free-slip both at the surface and core-mantle boundary. The surface temperature is fixed at 250 K while the core-mantle boundary temperature is decreasing as the core cools. The core in our model acts as an isothermal heat bath, whose cooling depends on its density and heat capacity (Table 3.1). Given the high viscosity of the surface, its velocity is negligible, and the assumption of a no-slip surface has no effect on the results. A schematic of the model geometry can be found in Fig (3.1). For comparison with our convection results, we also ran a set of simulations that are purely conductive.

# 3.2 Initial conditions

## 3.2.1 Heat sources

Whether the Moon possesses a composition in refractory elements similar to the Earth or is enriched by a factor of about two is a long lasting debate (Taylor, 1982;
<table><caption>Table 3.1: Model parameters</caption><thead><tr><th>Symbol</th><th>Description</th><th>Value</th><th>Scaling</th></tr></thead><tbody><tr><td>R<sub>p</sub></td><td>Planet radius</td><td>1740 km</td><td>R<sub>p</sub> - R<sub>c</sub></td></tr><tr><td>R<sub>c</sub></td><td>Core radius</td><td>390 km</td><td>R<sub>p</sub> - R<sub>c</sub></td></tr><tr><td>D<sub>c</sub></td><td>Crustal thickness</td><td>40 km</td><td>R<sub>p</sub> - R<sub>c</sub></td></tr><tr><td>D<sub>K</sub></td><td>KREEP layer thickness</td><td>10 km</td><td>R<sub>p</sub> - R<sub>c</sub></td></tr><tr><td>T<sub>surf</sub></td><td>Surface temperature</td><td>250 K</td><td>ΔT</td></tr><tr><td>T<sub>0</sub></td><td>Reference temperature</td><td>1600 K</td><td>ΔT</td></tr><tr><td>η<sub>0</sub></td><td>Reference viscosity</td><td>10<sup>21</sup> Pa s</td><td>η<sub>0</sub></td></tr><tr><td>η<sub>max</sub></td><td>Maximum viscosity</td><td>10<sup>28</sup> Pa s</td><td>η<sub>0</sub></td></tr><tr><td>E</td><td>Activation energy</td><td>3 10<sup>5</sup> J mol<sup>-1</sup></td><td>RΔT</td></tr><tr><td>L</td><td>Latent heat of melting</td><td>6 10<sup>5</sup> J kg<sup>-1</sup></td><td>c<sub>p</sub>ΔT</td></tr><tr><td>R</td><td>Universal gas constant</td><td>8.314 J mol<sup>-1</sup> K<sup>-1</sup></td><td>-</td></tr><tr><td>c<sub>p</sub></td><td>Mantle specific heat capacity</td><td>1000 J kg<sup>-1</sup> K<sup>-1</sup></td><td>c<sub>p</sub></td></tr><tr><td>k<sub>c</sub></td><td>Crust thermal conductivity</td><td>1.5 W m<sup>-1</sup> K<sup>-1</sup></td><td>k<sub>m</sub></td></tr><tr><td>k<sub>m</sub></td><td>Mantle thermal conductivity</td><td>3 W m<sup>-1</sup> K<sup>-1</sup></td><td>k<sub>m</sub></td></tr><tr><td>κ<sub>0</sub></td><td>Reference thermal diffusivity</td><td>10<sup>-6</sup> m<sup>2</sup> s<sup>-1</sup></td><td>k<sub>m</sub>/(ρ<sub>0</sub>c<sub>p</sub>)</td></tr><tr><td>ρ<sub>0</sub></td><td>Reference density</td><td>3400 kg m<sup>-3</sup></td><td>ρ<sub>0</sub></td></tr><tr><td>α<sub>0</sub></td><td>Thermal expansivity</td><td>2 10<sup>-5</sup> K<sup>-1</sup></td><td>α<sub>0</sub></td></tr><tr><td>c<sub>p,core</sub></td><td>Core specific heat capacity</td><td>800 J kg<sup>-1</sup> K<sup>-1</sup></td><td>c<sub>p</sub></td></tr><tr><td>ρ<sub>core</sub></td><td>Core density</td><td>7400 kg m<sup>-3</sup></td><td>-</td></tr><tr><td>k<sub>core</sub></td><td>Core thermal conductivity</td><td>25-50 W m<sup>-1</sup> K<sup>-1</sup></td><td>-</td></tr><tr><td>α<sub>core</sub></td><td>Core thermal expansivity</td><td>10<sup>-4</sup> K<sup>-1</sup></td><td>-</td></tr><tr><td>ΔT</td><td>Temperature drop across the mantle</td><td>1750 K</td><td>-</td></tr><tr><td>Δρ</td><td>Peridotite to harzburgite density difference</td><td>60 kg m<sup>-3</sup></td><td>-</td></tr><tr><td>C<sub>0</sub></td><td>Peridotite to harzburgite melt fraction</td><td>0.3</td><td>-</td></tr><tr><td>Q<sub>m</sub></td><td>Radiogenic heating</td><td>cf. Table 3.2</td><td>k<sub>m</sub>/(D<sup>2</sup>ρ<sub>0</sub>)</td></tr><tr><td>g</td><td>Surface gravity acceleration</td><td>1.62 m s<sup>-2</sup></td><td>g<sub>0</sub></td></tr><tr><td>Ra</td><td>Reference thermal Rayleigh number</td><td>2.24 10<sup>5</sup></td><td>-</td></tr><tr><td>Ra<sub>C</sub></td><td>Reference chemical Rayleigh number</td><td>2.4 10<sup>5</sup></td><td>-</td></tr></tbody></table>

Warren, 2001; Taylor et al., 2006). This is critical for thermal evolution models of the Moon, as the bulk abundance of uranium is estimated to range from 17 ppb for primitive Earth values to 34 ppb in the enriched case.

The abundance of heat producing elements in the lunar mantle can be estimated from the composition of primary mantle melts and mineral/melt partition coefficients. As the most primitive primary melts, the picritic glasses, have a large range of compositions, different authors have advocated different values for the bulk mantle composition. We consider two representative numbers: 25 ppb thorium from Warren and Wasson (1979), and 40 ppb thorium from Jolliff et al. (2000). We estimate the present day mantle abundance of U and K using an average Th/U ratio of 3.7, and a K/U ratio of 2500 (Taylor, 1982).

Our model contains a layer of KREEP-rich material that is concentrated in the Procellarum KREEP Terrane (PKT). We note simply that if a global 2 km KREEP-layer (as was once envisioned during the Apollo era (Warren and Wasson, 1979)) was segregated beneath the PKT, such a layer would be about 10 km thick. We thus
assume that the KREEP layer is 10 km thick, as in Wieczorek and Phillips (2000), but acknowledge that its true thickness is somewhat uncertain. We choose the uranium concentration of this layer to be 3.4 ppm, which corresponds to the average uranium content of Apollo 15 KREEP basalts (Korotev, 2000). Though KREEP basalts are rare in the Apollo sample collection, this composition is representative of a large portion of the materials that were excavated by the Imbrium impact basin, which lies within the PKT. Because the radial resolution of our model is 20 km, when the KREEP layer is thinner than this, we place the equivalent amount of heat producing elements in a 20 km thick layer.

We model the KREEP layer as a cylindrical cap of either 40 or 80° diameter and the equivalent of 10 km KREEP basalt is placed either below the crust, at the bottom of the crust, or redistributed over the entire crust (see Figure 3.1). When the KREEP layer is redistributed over the entire crust, the resulting concentration is obtained from mass conservation. The crustal regions that are not enriched in heat sources have a uranium concentration of 0.14 ppm (Jolliff et al., 2000), which is consistent with a viscoelastic study of lunar impact basins deformations that suggests that the farside crust is more than 10 times depleted in heat sources with respect to the PKT region (Kamata et al., 2013).

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

Figure 3.1: Schematic diagram of the thermochemical convection model setup (not to scale). The equivalent of 10 km of KREEP basalt is placed below a 40 km thick crust (blue), in the lower 20 km of the crust (cross hatch), or redistributed over the entire crust (orange). The lateral extent of the PKT is 80° in diameter and the radial resolution is 20 km. The temperature of the core $T_c$ evolves with time, but does not include the energetics of core crystallization.
The KREEP layer is not fixed to its initial position and is free to flow with time. When the KREEP layer is placed within the mantle, we assume that its density is equal to that of the surrounding mantle, even though the Mg-rich composition of KREEP basalt implies a density that is lower than typical mantle materials. Since our KREEP layer is only one grid element thick, we have run higher resolution simulations in 2D with a 5 km radial resolution to better quantify how this material is displaced with time. Our results show that only a very small portion of the KREEP layer at its edge gets entrained into the mantle flow, and we thus expect that our poor resolution of the KREEP layer will not have any significant influence on the results.

Our assumed distribution of heat sources leads to extreme cases with bulk uranium contents between 19.6 ppb for a 40° diameter PKT and a mantle with the lowest concentration of uranium, and 28.7 ppb for a 80° diameter PKT and the highest concentration of uranium in the mantle. The bulk compositions lie well within the Earth-like and enriched scenarios discussed above. Our preferred model has a 80° diameter PKT, a mantle with the most depleted concentrations, and a KREEP layer emplaced below the crust. For this model the bulk uranium concentration is 19.6 ppb, which is similar to the value for the bulk silicate Earth, and which is consistent with recent bulk Moon composition estimates derived from the GRAIL mission (Wieczorek et al., 2013). This corresponds to having about one third of the Moon's heat sources in the PKT region. A summary of the heat-source concentrations can be found in Table 3.2.

<table><caption>Table 3.2: Uranium concentration for the different model layers. For all cases, K/U = 2500 and Th/U = 3.7 (Taylor, 1982). In this table, ‘KREEP’ refers to the 10-km layer below the crust, ‘enriched crust’ refers to the case where the KREEP layer is redistributed over the entire crust, and ‘crust’ refers to the unenriched part of the crust. ‘0’ and ‘1’ refer to the mantle heat source content, ‘0’ being the lowest and ‘1’ the highest. ‘L’ and ‘S’ refer to the PKT diameter (ie. large and small). ‘B’, ‘D’ and ‘W’ correspond to the KREEP emplacement geometry: below the crust, distributed at the bottom of the crust, or within the entire crust, respectively.</caption><thead><tr><th>Model</th><th>PKT<br>diameter</th><th>Mantle<br>ppb</th><th>KREEP<br>ppm</th><th>Enriched crust<br>ppm</th><th>Crust<br>ppm</th><th>Bulk silicate<br>ppb</th></tr></thead><tbody><tr><td>0LB</td><td>80°</td><td>6.8</td><td>3.4</td><td>-</td><td>0.14</td><td>25.1</td></tr><tr><td>1LB</td><td>80°</td><td>10.8</td><td>3.4</td><td>-</td><td>0.14</td><td>28.7</td></tr><tr><td>0SB</td><td>40°</td><td>6.8</td><td>3.4</td><td>-</td><td>0.14</td><td>19.7</td></tr><tr><td>1SB</td><td>40°</td><td>10.8</td><td>3.4</td><td>-</td><td>0.14</td><td>23.4</td></tr><tr><td>0LD</td><td>80°</td><td>6.8</td><td>3.4</td><td>-</td><td>0.14</td><td>24.5</td></tr><tr><td>1LD</td><td>80°</td><td>10.8</td><td>3.4</td><td>-</td><td>0.14</td><td>28.2</td></tr><tr><td>0SD</td><td>40°</td><td>6.8</td><td>3.4</td><td>-</td><td>0.14</td><td>19.6</td></tr><tr><td>1SD</td><td>40°</td><td>10.8</td><td>3.4</td><td>-</td><td>0.14</td><td>23.3</td></tr><tr><td>0LW</td><td>80°</td><td>6.8</td><td>-</td><td>0.82</td><td>0.14</td><td>25.1</td></tr><tr><td>1LW</td><td>80°</td><td>10.8</td><td>-</td><td>0.82</td><td>0.14</td><td>28.7</td></tr><tr><td>0SW</td><td>40°</td><td>6.8</td><td>-</td><td>0.82</td><td>0.14</td><td>19.7</td></tr><tr><td>1SW</td><td>40°</td><td>10.8</td><td>-</td><td>0.82</td><td>0.14</td><td>23.4</td></tr></tbody></table>
### 3.2.2 Initial temperature profile

To cover the range of possible Moon formation and initial differentiation scenarios we considered a range of initial temperature profiles. Every profile starts initially from a surface temperature of 250 K at 4.5 Ga and ends at a value of 2000 K at the core-mantle boundary in order to account for an excess temperature in the core due to differentiation of several hundred degrees (Konrad and Spohn, 1997). The core is then allowed to cool as an isothermal heat bath. The choice of the temperature increase with depth in the mantle is then linked to an assumed magma ocean crystallization scenario. The dynamics of such a system are complex, and we therefore consider two extreme cases. The ‘cold’ case assumes an adiabatic gradient within the whole mantle (as in Wieczorek and Phillips, 2000) corresponding to a well-mixed, initially convecting interior. We note that, as gravity is assumed constant as a function of depth, we slightly overestimate the adiabatic gradient. The temperature profile of the ‘hot’ case follows the mantle solidus for the first 700 km, below which the mantle is adiabatic. For this scenario, the region at the solidus represents the upper portion of the solidified magma ocean that did not convectively readjust. In both cases, the initial temperature profile is linear in the crust, reaching the mantle solidus at the crust-mantle interface. Our preferred case ends-up having an intermediate profile, following the solidus down to 350 km depth (Figure 3.2.2). We do not consider stable post-magma ocean overturn temperature profiles here (such as in Hess and Parmentier, 1995; Elkins-Tanton et al., 2011), in an attempt to keep the setting of our model as simple as possible. We note simply that if there were any density gradients with depth, that this would limit the amount of convection, and that the thermal evolution would approach a model that was purely conductive (as in Wieczorek and Phillips, 2000).
<figure><img src="image_21.png" /><figcaption>Figure 3.2: Initial temperature profiles. For all models the temperature at the crust-mantle interface is set to the solidus of the mantle. The cold case starts with an adiabatic gradient in the mantle, the intermediate case follows the mantle solidus down to 350 km and is adiabatic below, and the hot case follows the mantle adiabat to 700 km and is adiabatic below. Our preferred model corresponds to an intermediate case following the solidus down to 350 km depth (see text for details). Solidus and liquidus are taken from *Hirschmann* (2000).</figcaption></figure>
# Thermo-chemical evolution results

## Contents

<table><tr><td>4.1</td><td>Thermal evolution</td><td>35</td></tr><tr><td>4.2</td><td>Magmatism</td><td>41</td></tr><tr><td>4.3</td><td>Predicted gravity field</td><td>45</td></tr><tr><td>4.4</td><td>Dynamo generation</td><td>50</td></tr><tr><td>4.5</td><td>Discussion</td><td>54</td></tr><tr><td>4.6</td><td>Conclusions & perspectives</td><td>55</td></tr></table>

Before describing our simulation results, we describe the conditions that we consider constitute a successful model. First, a successful model must give rise to volcanism primarily within the PKT, and only a small amount exterior to this region and on the farside. Second, mare volcanism must last for a significant amount of time, from > 4 Ga to about 1 Ga. Third, the total volume of extrusive lava flows has been estimated to lie between $2 \cdot 10^6$ and $7 \cdot 10^6$ km³ (e.g., Wieczorek et al., 2006; Shearer et al., 2006). If the ratio of intrusive to extrusive magmatism is 5:1 (White et al., 2006), then the total volume of generated magma should lie somewhere between $1.2 \cdot 10^7$ and $4.2 \cdot 10^7$ km³. Given that the intrusive to extrusive ratio is not known with any certainty, and that the thicknesses of the mare are also uncertain, this estimate should be considered uncertain by a factor of about 5. Fourth, based on the ages of mare basalts in Figure 2.1, the youngest lavas should erupt in the center of the PKT. Finally, the surface heat flow must be compatible with the measurements made at the Apollo 15 and 17 landing sites.

## 4.1 Thermal evolution

Our nominal model (model ‘0LB’, see Table 4.1) started with an intermediate initial temperature profile, with the KREEP layer located below the 40 km thick crust and has Earth-like bulk abundances of refractory elements About one third of the heat source budget concentrated in a small region on the nearside hemisphere, and as shown in Figure 4.1 as a series of temperature slices, this has a dramatic influence on the thermal history of the Moon. Other examples of thermal evolutions can be found in Figures 4.5 to 4.7 at the end of the section. Contrary to thermal evolution models with a symmetric distribution of heat sources, where the global cooling rate dictates the mantle behavior, the PKT region is the driver in our case as it heats
<table><caption>Table 4.1: Models summary including the main melting characteristics. ‘I’, ‘C’ and ‘H’ correspond to intermediate, cold and hot initial temperature profile respectively. The rest of the nomenclature is the same as in Table 2.</caption><thead><tr><th rowspan="2">Model</th><th colspan="3">Nearside</th><th colspan="3">Farside</th></tr><tr><th>Duration</th><th>Max depth</th><th>Total vol.</th><th>Duration</th><th>Max depth</th><th>Total vol.</th></tr><tr><th></th><th>Ga</th><th>km</th><th>km<sup>3</sup></th><th>Ga</th><th>km</th><th>km<sup>3</sup></th></tr></thead><tbody><tr><td>I-0LB</td><td>4.5 - 0.2</td><td>660</td><td>3.7e8</td><td>3.8 - 3.1</td><td>680</td><td>3.5e7</td></tr><tr><td>I-0LD</td><td>3.8 - 2.5</td><td>700</td><td>5.2e7</td><td>3.8 - 3.1</td><td>680</td><td>2.5e7</td></tr><tr><td>I-0LW</td><td>3.8 - 1.6</td><td>700</td><td>7.5e7</td><td>3.8 - 3.0</td><td>680</td><td>2.5e7</td></tr><tr><td>H-0LB</td><td>4.5 - 0.0</td><td>980</td><td>6.7e8</td><td>4.1 - 1.6</td><td>990</td><td>4.3e8</td></tr><tr><td>H-0LD</td><td>4.3 - 0.6</td><td>880</td><td>7.0e8</td><td>4.3 - 1.6</td><td>860</td><td>6.6e8</td></tr><tr><td>H-1LW</td><td>4.3 - 0.0</td><td>1150</td><td>1.4e9</td><td>4.3 - 0.0</td><td>1160</td><td>1.2e9</td></tr><tr><td>C-0LD</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>C-0LB</td><td>4.5 - 0.4</td><td>540</td><td>1.6e8</td><td>-</td><td>-</td><td>-</td></tr><tr><td>C-1LW</td><td>2.6 - 0.5</td><td>800</td><td>2.0e8</td><td>1.3 - 1.1</td><td>770</td><td>2.7e4</td></tr></tbody></table>

the underlying mantle. This is not a new result (e.g., Wieczorek and Phillips, 2000), but contrary to previous conductive models that studied the PKT, our simulations show that heating of the underlying mantle helps develop a stronger convection in the nearside mantle than in the farside. As shown in Figure (4.2), this leads to a more efficient cooling than would occur in a purely conductive model. For a description of symmetrical lunar thermal evolutions models, the reader is referred to Spohn et al. (2001) and Ziethe et al. (2009).

Figure (4.3) shows how the present day heat flux varies as a function of distance from the center of the Procellarum KREEP Terrane. In the center of the PKT, we find a maximum value of about 25 mW/m² and a background value outside of this terrane of about 10 mW/m². These values are similar to the measurements made at the Apollo 15 and 17 landing sites, i.e. 21±3 mW/m² and 14±2, respectively (Langseth et al., 1976), especially when considering that the precise locations of these landing sites with respect to the edge of the PKT is somewhat uncertain. As a result of lateral variations in megaregolith thickness (Warren and Rasmussen, 1987), the Apollo heat flow estimates should perhaps be corrected downward to 18 and 12 mW/m². As shown in Figure (4.4a), the heat flux in the PKT is predicted to be about a factor of two greater than in the surrounding highlands during most of the Moon's thermal evolution. The present day heat flux is found to be rather insensitive to the various initial conditions that were considered (Figure 4.3).

The high concentration of heat sources in the PKT leads to nearly instantaneous partial melting of the underlying mantle. For this reference model, after about 600 million years, melting is also observed on the farside, but lasts for only a few 100 million years. The time at which farside melting occurs is related to the time to initiate global convection, and hence is directly related to the Rayleigh number (e.g. Schubert et al., 1969). On the nearside, the region where melt occurs increases with
<figure><img src="image_22.png" /></figure>

Figure 4.1: Temperature cross sections of the lunar mantle for a complete thermal evolution for the case with an intermediate initial temperature profile and the KREEP layer located below the crust (model ‘T-0LB’). Numbers correspond to time before present in Ga. The black circle is the lunar core and white corresponds to regions that are partially molten. The streamlines are shown as dashed lines.

depth over time, reaching 600 km after 2 billion years. Afterwards, this region starts to shrink, with the last magmas being confined to depths of about 500 km a few 100 million years before present. The PKT-induced mantle circulation on the nearside
<figure><img src="image_23.png" /><figcaption>Figure 4.2: Present day temperature profiles beneath the center of the PKT and beneath the center of the farside (solid lines) for our preferred case (model T-0LB) compared with the results of a purely conductive simulation (dotted).</figcaption></figure>

produces two slowly convecting cells that bring down cold materials approximately 80° away from the center of the PKT at about 0.25 cm/yr. Overall, mantle mixing is inefficient, with velocities below the lithosphere reaching 0.2 cm/yr on average on the nearside and 0.15 cm/yr on the farside 3.5 Ga ago, decreasing to about 0.05 cm/yr at the present day. This is an order of magnitude slower than on the Earth with velocities estimated from sinking slabs (e.g., van der Meer et al., 2009) and corresponds to about one complete overturn of the lunar mantle in 4.5 billion years. As a result, even when KREEP is emplaced below the crust, with the same density as the mantle, it is dynamically stable over the whole evolution of the Moon.

We estimate the average change in radius of the Moon on its near and farside hemispheres due to the changes in density that occur from thermal expansion and compositional depletion. The estimated change of radius from mapping of contrac- tional and compressional features takes only into account changes since the end of basin formation at about 3.8 Ga. As shown in Figure (4.4b), up to 7 km of expan- sion occurs on the nearside during the first 600 million years, but surface expressions from this epoch are not expected to be seen today. Since 3.8 Ga, about 1 km of expansion occured on the nearside until about 3 Ga, at which point contraction started to dominate. A present-day total contraction of about 2 km is expected to be seen on the nearside (Fig 4.4b). In contrast, on the farside, contraction is found to be continuous throughout lunar history and produces a present-day contraction
<figure><img src="image_24.png" /></figure>

Figure 4.3: Average surface heat flow as a function of distance from the center of the PKT for several simulations. The solid lines correspond to cases where KREEP is emplaced below the crust and the dotted line when it is redistributed within the entire crust. ‘C-’, ‘I-’ and ‘H-’ correspond to a cold, intermediate and hot initial temperature profile. The nomenclature follows Table 3. The green dot and diamond are the surface heat flow estimates at the approximate locations of the Apollo 15 and 17 landing sites, respectively (Langseth et al., 1976) and the gray boxes include the correction due to heat flow focussing from Warren and Rasmussen (1987).

of about 3 km. The presence of an early period of expansion on the nearside is in agreement with what appear to be giant magmatic dikes, as observed in the GRAIL data by Andrews-Hanna et al. (2013). However, as his model requires global expansion to account for similar features on the farside hemisphere, this is not compatible with our model results which show that the farside hemisphere was always contracting. Finally, our predicted contraction on the nearside is in agreement with the global value of about 1 km estimated by Watters et al. (2010), but continuous contraction of the farside leads to a much larger change in radius than is implied by the observed contractional features. Nevertheless, these authors state that their estimate is a lower limit and may be uncertain by a factor of about 2.

We ran a series of simulations varying the heat source content and distribution, the size of the PKT and the initial temperature profile (see Table 4.1 for a summary and Figures A1 to A3 for examples of thermal evolutions). Three major consequences always arise from the initial hypothesis of localized heat sources. First,
<figure><img src="image_25.png" /><figcaption>Figure 4.4: (a) Surface heat flow in the center of the PKT and on the farside as a function of time before present. (b) Radius change in the center of the PKT and farside due to secular cooling as a function of time before present. The vertical line denotes the approximate end of the putative late heavy bombardment. (c) Cumulate melt volume averaged over the near and farside (solid lines). Total mare basalt volume estimate (thin black) is from Shearer et al. (2006), and the maximum melt volume (dashed line) takes into account a 5:1 intrusive to extrusive volcanism ratio (White et al., 2006). (d) Average melt production rate on the near and farside. These figures were obtained for our preferred case (model T-0LB).</figcaption></figure>
melt is localized on the nearside while the farside remains mostly undisturbed by the PKT. The amount of farside melting is highly dependent on the mantle's initial temperature profile. Second, a thermal anomaly is preserved until the present day in the mantle underlying the PKT. And third, the asymmetric heating from the PKT has an influence down to the core-mantle boundary. The consequence of these results are described in more detail in the following subsections, with an emphasis placed on understanding the general attributes of models with an enhancement of heat production in the Procellarum KREEP Terrane, and not models that fit a specific observation exactly.

## 4.2 Magmatism

The main influence on the Moon's magmatic history is the assumed initial temperature profile. A higher initial temperature leads to more melting, as expected, but also to earlier and deeper melting. This is a very strong constraint, as models with a cold initial temperature profile fail to produce any significant amount of melt on the farside. On the other hand, the hot initial start tends to produce too much melting (on the order of $10^9$ km$^3$ on the farside alone). As discussed later, the assumed solidus for the mantle has a large influence on the total volume of magma that is generated as well, but this (and the mantle composition) is poorly constrained.

The exact emplacement and configuration of the KREEP layer does not have much influence on the farside evolution, but it does have a strong influence on the nearside magmatic history. When the KREEP layer is emplaced below the crust, melting starts immediately and stops a few hundred million years ago for all models, independently of the initial temperature profile. The influence of the initial temperature profile simply reflects the maximum depth of melting. While the hot model produces melt down to about 1000 km, the intermediate and cold cases reach between 600 and 500 km depth, respectively. When the KREEP layer is initially evenly redistributed in the entire crust, it has a smaller effect on the mantle as heat is lost more easily to space. In this case, the melting duration and depth are highly dependent on the initial temperature profile: a hot initial temperature profile leads to melting at greater depths, but also to a longer duration than for a cold initial start. The case where the KREEP layer is emplaced in the lower portion of the crust is very similar to the case where it is redistributed in the entire crust: the exact melting time and volume is subject to a trade-off between initial temperature and bulk heat sources content. Nevertheless, none of the models that we explored with a KREEP layer either redistributed within the entire crust or emplaced at the bottom of the crust succeeded in matching the estimated timing and amount of volcanism on the nearside (see Table 4.1 for a summary and Figure A4 for a comparison). Under the assumptions of our modeling, this suggests that KREEP was most probably emplaced below the crust.

The model that best fits the estimates of near and farside volcanism corresponds to an intermediate temperature profile with the KREEP layer emplaced directly
<figure><img src="image_26.png" /><figcaption>Figure 4.5: Temperature cross sections of the lunar mantle for a complete thermal evolution for the case with an intermediate initial temperature profile and the KREEP layer redistributed within crust (model ‘T-0LW’). Numbers correspond to time before present in Ga. The black circle is the lunar core and white corresponds to regions that are partially molten. The streamlines are shown as dashed lines.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 4.6: Temperature cross sections of the lunar mantle for a complete thermal evolution for the case with a cold initial temperature profile and the KREEP layer located below the crust (model ‘A-0LB’). Numbers correspond to time before present in Ga. The black circle is the lunar core and white corresponds to regions that are partially molten. The streamlines are shown as dashed lines.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Figure 4.7: Temperature cross sections of the lunar mantle for a complete thermal evolution for the case with a hot initial temperature profile and the KREEP layer located below the crust (model ‘D-0LB’). Numbers correspond to time before present in Ga. The black circle is the lunar core and white corresponds to regions that are partially molten. The streamlines are shown as dashed lines.</figcaption></figure>
below the crust (model ‘T-0LB’ in Table 4.1). As shown in Figure (4.4), melting then occurs on the nearside from 4.5 to 0.2 Ga ago, reaches about 600 km depth below the PKT and generates about $5 \cdot 10^8$ km$^3$ of magma with a peak in magma production between 4 and 3.5 Ga ago. Consistent with observations (Figure 2.1), the youngest lavas are found to erupt in the center of the PKT. Farside melting is characterized by a main volcanic phase starting 3.7 Ga ago, that generates $2 \cdot 10^7$ km$^3$ of magma (see Fig. 4.4 for details). For this particular model, about 10 times more melt is produced on the nearside than on the farside, consistent with current estimates (Wieczorek et al., 2001; Morota et al., 2009).

If all the generated magma was assumed to erupt onto the surface, the thickness of the lava flows would be on average about 30 km on the nearside and 10 km on the farside. However, if we assume an intrusion to extrusion ratio of five, these thicknesses would be reduced to 5 and 1.7 km, respectively. The estimated thickness of mare basalts vary from a few hundred meters to several kilometers in large basins (Williams and Zuber, 1998). The averages that we find for our nominal model are therefore probably a factor of 2 to 5 too large. Nevertheless, we note that our simulations will tend to overestimate the amount of melting, as we will discuss later on, because the increase in temperature of the solidus with mantle depletion is not considered in our models, nor is the fractionating of heat producing elements into the magma considered. Considering these uncertainties, we find the model we show in Fig. (4.4) to be largely compatible with the observations. A comparison of melt volumes on the near and farside for different initial conditions can be found in Figure 4.8.

The bulk composition of the Moon is a long standing debate. It was first argued to be enriched in refractory elements in comparison to the Earth (Taylor, 1982). With revised estimates of the surface heat flow, Warren and Rasmussen (1987) claimed that the bulk composition could be similar to that of the primitive Earth. We find that to stay within observational constraints on the volume of mare basalts and the timing of mare volcanism, the lower bulk heat source contents are preferable, which is consistent with the latest view of an Earth-like composition for the Moon (Wieczorek et al., 2013).

## 4.3 Predicted gravity field

For acceptable models, mantle melting should be almost over today as no present day volcanism is observed. Nevertheless, even though melting ceases before the present, we always observe a temperature anomaly below the PKT, and this has direct consequences on the Moon’s global gravity field. We calculate the predicted radial gravity anomalies, $g$, using the spherical harmonics expansion

$$g(r, \theta, \phi) = \frac{GM}{r^2} \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \left(\frac{R_0}{r}\right)^l (l+1)C_{lm}Y_{lm}(\theta, \phi), \quad (4.1)$$

where $r$, $\theta$ and $\phi$ are the coordinates of observation, $G$ is the gravitational
<figure><img src="image_29.png" /><figcaption>Figure 4.8: Cumulate melt volume on the nearside (left) and farside (right). The color indicates the initial temperature profile (black is intermediate, blue is cold and red is hot). Solid lines correspond to cases with the KREEP layer located below the crust whereas the dashed line corresponds to KREEP initially redistributed within the crust. The nomenclature follows Table 4.1.</figcaption></figure>

constant, $M$ the mass of the Moon, $R_0$ the reference radius of the spherical harmonic coefficients $C_{lm}$, and $Y_{lm}$ the spherical harmonic functions of degree $l$ and order $m$ (e.g., Wieczorek, 2007). The radial gravity anomalies are positive when directed downward. The spherical harmonic coefficients $C_{lm}$ have contributions from several sources, such as buoyancy generated by temperature and compositional variations, thermal uplift of the surface, and the extruded mare basaltic lava flows.

For the temperature and compositional contributions, the spherical harmonic coefficients of the gravitational field are obtained by integrating the contributions of a large number of thin spherical shells:

$$C_{lm} = \frac{4\pi}{M(2l+1)} \sum_{i=1}^{N} r'_{i}^{2} \left(\frac{r'_{i}}{R_{0}}\right)^{l} \rho_{lm}(r'_{i}) \Delta R. \quad (4.2)$$

where $\rho_{lm}(r)$ are the spherical harmonics coefficients of the density distribution at radius $r$. For all the following calculations, we removed the degree-0 coefficient, which correspond to a globally constant gravitational contribution. As shown in Figure (4.9a), the contribution to the gravity field from temperature induced lateral variations in density is maximal in the center of the PKT and has a value of about -500 mGal using a thermal expansion coefficient of $2 \cdot 10^{-5} K^{-1}$. As can be seen in the temperature-density relation of Eq. (3.7) and the linearity between density and gravity in Eq. (A.8), the magnitude of this gravity anomaly is linearly related to the thermal expansivity $\alpha$. If we were to have used $\alpha = 3 \cdot 10^{-5} K^{-1}$, as in Grimm (2013), we would obtain a maximum value for the gravity anomaly of -800 mGal, consistent with his findings. In addition to this effect of temperature on the Moon's gravity field, partial melting of a peridotitic mantle leads to the formation of harzburgite residue with a lower density, which will influence the observed field
as well. This density contribution is calculated using the chemical depletion part of Eq. 7. This contribution was neglected in the work of Grimm (2013), and we find that this compositional effect contributes an additional negative gravity anomaly of about -600 mGal (Fig 4.9b).

In addition to the temperature and compositional contributions, a positive con-
tribution from surface uplift will also occur. Under the assumption of a strengthless
lithosphere, mantle flow driven by thermochemical density variations will produce
surface deformation. Although this is not taken into account in our convection
model explicitly, it can be estimated a posteriori using a well-established formalism
(e.g. *Richards and Hager*, 1984). The total predicted topography is found to be
about 4.8 km (7.2 km when using $\alpha = 3 \cdot 10^{-5}$ K$^{-1}$), in close agreement with the
prior estimate of Grimm (2013). Using a modification of Eq. 13 (e.g., *Wieczorek,
2007*), we calculate in Figure (4.9c) the gravity contribution of this surface relief to
be +700 mGals. However, the presence of an elastic lithosphere could have a signif-
icant effect in reducing this dynamic topography (*Zhong, 2002; Golle et al., 2012*).
If the lithosphere were perfectly rigid, there would be zero uplift, and hence there
would be no positive gravity anomaly associated with the surface. The strengthless
and rigid lithosphere scenarios represent two extreme scenarios that will bracket
all intermediate cases. We note here that our case with a rigid lithosphere should
be treated with caution: in incompressible models, considering the effect of thermal
expansion without including surface uplift in essence removes mass from the system.

In this work, we do not perform the exact flexure calculation as there are several
complications that would make this calculation highly uncertain. For example, a
significant part of the load acting on the lithosphere comes from degree-0 and degree-
1 terms, but these are not taken into account in standard spherical shell flexural
models, such as those based on Turcotte et al. (1981). Furthermore, the degree-0
term would act to place the lithosphere into extension or compression, which could
give rise to fractures that would invalidate the spherical shell flexure formalism, and
which would also set up in-plane forces that are not generally considered. Further-
more, the elastic thickness of the lithosphere is expected to vary laterally, which
again would invalidate the assumptions of most simple elastic shell flexure models
(though see Beuthe, 2008). As an example, if we take the 750 K isotherm as rep-
resenting the elastic thickness of the lithosphere (Watts and Zhong, 2000), at 3.5
billion years ago during the main phase of mare volcanism, the elastic thickness is
predicted to be less than 30 km in the PKT, and 80 km for the farside highlands,
consistent with estimates from Crosby and McKenzie (2005). Finally, any flexure
calculation would need to consider the load caused by the mare basalts and intrusive
magmas, which are not well constrained.

The final gravitational signature to consider is that due to the eruption of dense lavas onto the surface, and the formation of dense intrusions in the crust. The magnitude of this anomaly, however, will depend sensitively on their assumed compensation state. If the lithosphere attained a state of isostatic equilibrium, perhaps because of the temperatures associated with melting in the underlying mantle, the final gravity anomaly would be nearly zero (maximum of about 30 mGal for 5
<figure><img src="image_30.png" /><figcaption>Figure 4.9: Gravity field of the Moon. Contributions from (a) thermal expansion, (b) compositional depletion in the mantle, (c) uplift of the surface for a strengthless lithosphere, and (d) uncompensated mare basalts with a maximum thickness of 5 km. For these contributions the degree-0 term has been removed, but the degree-1 term is included. Predicted gravity anomaly (e) for a strengthless lithosphere with isostatically compensated mare basalts, and (f) for a completely rigid lithosphere with uncompensated mare basalts and zero thermal uplift. For the last two figures, the degree-1 gravity field (which corresponds to change in the planets center of mass) has been removed. All figures are presented in a Mollweide equal-area projection centered on the PKT. The black circle corresponds to the extent of the PKT. The results are calculated for our preferred case (model T-0LB).</figcaption></figure>
km-thick maria). However, if the lavas were perfectly uncompensated, the gravity anomaly would be about +140 mGals for each km of uncompensated lava. If the mare basalts were about 5 km thick, as observed in the previous section, then this gravity anomaly could reach up to +700 mGals (Figure 4.9d). Considering the intrusive components would generate an even larger anomaly. Nevertheless, given that the mare do not possess distinct gravity signatures (with the exception of mascon basins, see Zuber et al. (2013)), it is probable that the loads on the lithosphere resulting from intrusive and extrusive lavas are at least partially compensated. We note in addition that downward attenuation of the anomaly in Eq. 13 implies that the positive contribution from dense material taken from depth will be larger than its negative counterpart.

In this paragraph, we present the total predicted gravity field corresponding to two extreme cases. We first assume that the lithosphere did not possess significant strength at the time of mare volcanism (i.e., the basalts are assumed to be compensated), and plot the total gravity anomaly in Figure (4.9e). For our final gravity map, we remove the degree-1 gravity signature, which corresponds to a slight change in the center of mass of the planet of about 1 km towards the nearside. The combined signal is seen to be small, with an anomaly of about -200 mGals within the PKT. If only about 1-2 km of basalts (or their intrusive equivalent) were partially uncompensated, the total gravity anomaly would be close to zero, consistent with the observations (the observed gravity anomalies in the PKT range from about -100 to 200 mGals, excluding the mascon basins). Second, if we assume that the lithosphere were completely rigid, the lack of thermal uplift in the PKT, combined with 5 km of uncompensated basalts, would generate an anomaly of -400 mGal, as shown in Figure (4.9f). If we were to include the gravity contribution from the uncompensated magmatic intrusions at the crust-mantle interface, the gravity anomaly could be as high as +3000 mGal, which tends to argue against the scenario of a rigid lithosphere. In conclusion, under reasonable assumptions, it is possible to find models that do not predict a large gravity anomaly to be present within the PKT.

In contrast to the gravity field, our thermal model does predict surface topography that might be inconsistent with the observations. For example, the model presented in Figure (4.9e), where the lithosphere possesses no strength in the PKT, dynamic topography generated by thermal expansion in the underlying mantle predicts about 5 km of uplift (see also Grimm, 2013). Since this model also generates a shift in the center of mass of about 1 km toward the nearside, this topographic uplift in the PKT should be reduced by the same amount. Still, given that the PKT is located in a region where the surface elevations are about 2 km below the average, this is a potentially discrepant model result. If, on the other hand, we assume that the lithosphere is perfectly rigid, there will be no topographic uplift in the PKT. This model predicts a change in the center of mass of about 2 km toward the farside hemisphere, which would result in an uplift of 2 km on the nearside, which is also inconsistent with the observations.

The magnitude of our modeled contributions depends both on the depth and amplitude of the subsurface thermal anomaly as well as on the lithospheric thickness.
Models where the KREEP layer is distributed within the crust lose more heat to space and therefore tend to produce a smaller thermal anomaly and surface uplift than models where the KREEP layer is emplaced below the crust. The model ‘T-0LW’ for example predicts about 2.5 km uplift and surface gravity anomalies on the order of ± 40 mGal. Similarly, a lower initial temperature model exhibits a smaller anomaly at the present day. Model ‘A-0LB’ for example predicts gravity anomalies that are about 60% smaller than model ‘T-0LB’.

The most recent spacecraft observations show no significant gravity anomaly associated with Oceanus Procellurum, which suggests that it is probably regionally compensated (Zuber et al., 2013). The observed topography in the PKT however, is negative, not positive as might be expected from thermal uplift. This can be seen in Figure (4.10). We offer two possible explanations for this apparent conundrum. First, if the low elevations of the nearside crust are a result of crustal thickness variations (Wieczorek et al., 2013), this could easily mask the topographic signature of thermal uplift. Using a simple Airy isostasy model, a reduction in the crustal thickness by about 20 km would reduce the surface topography by 5 km. Second, since the high-titanium mare basalts erupted only on the nearside hemisphere, it is possible that the nearside mantle is denser than the farside because of its higher titanium content. This would give rise to a shift in the center of mass towards the nearside, which would reduce the nearside surface topography by the same amount.

Our observations are in general agreement with those of Grimm (2013), but in contrast to that author’s viewpoint, we argue that the small gravity anomalies and/or positive topography are not necessarily in contradiction with the hypothesis that about 10-km of KREEP-rich materials are present within, or just below, the crust of the Procellarum KREEP Terrane. Though minor differences are found with his model regarding the amplitude of the gravity and topography anomalies, these are largely because his study did not consider the degree-1 gravity signature of the planet (his calculations were done in a cartesian geometry), and also because his study did not consider gravity signals resulting from the compositional depletion of the mantle.

## 4.4 Dynamo generation

Figure (4.11) shows the average heat flow out of the core as a function of time. Similar to symmetric thermal evolution models of the Moon (Konrad and Spohn, 1997; Stegman et al., 2003), the core heat flux drops quickly below the adiabat threshold (within about 200 million years), below which a purely thermal dynamo is no longer possible. The adiabatic heat flow depends on the core properties and the largest uncertainty comes from thermal conductivity. We use the same range of plausible core adiabatic heat flux as Stegman et al. (2003), which are between 5 and 9 mW/m². The effect of the PKT on the core heat flow is small. The small peak at about 4 Ga is due to the onset of mantle convection and the observed difference between near and farside is due to the thermal perturbation from the
<figure><img src="image_31.png" /><figcaption>Figure 4.10: (A) Free-air and (B) Bouguer gravity anomaly maps from GRAIL model GL0420A. Maps are in Molleweide projection centered on 270E and show the nearside on the right and farside on the left. (Zuber et al., 2013).</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 4.11: Averaged CMB heat flow for the near and farside hemispheres as a function of time, compared with the average for a purely conductive case. These results are for our preferred case (model T-0LB).</figcaption></figure>

PKT reaching the core. The exact timing of the peak near 4 Ga is dependent on several parameters that control the onset of convection, such as viscosity, and the initial conditions of the thermal model. Nevertheless, even though the effect of the PKT on the average heat flow is small, its spatial distribution on the core-mantle boundary is asymmetrical, as shown in Figure (4.12). Directly below the PKT, the core is heated while much more heat is extracted from the surroundings.

A long-wavelength pattern could perhaps change the core flow patterns and influence the strength and geometry of a dynamo generated magnetic field. Takahashi and Tsunakawa (2009) have shown that a degree-1 asymmetrical heat flow pattern could lead to a dipolar magnetic field with the dipole axis circulating about a great circle about every 100 years. This magnetic pole variability would certainly influence rock magnetization. Crustal rocks would not be magnetized by a directionally stable field, complicating paleomagnetic analyses and interpretations of crustal magnetic anomalies.

We next predict the surface magnetic field at the surface using our core-mantle boudary heat flux and the scaling law for magnetic field strength from Christensen et al. (2008)

$$B = f \sqrt{0.63} \sqrt{2\mu_0\rho}^{1/6} (Fq_0)^{1/3} \left(\frac{R_c}{R_p}\right)^3. \qquad (4.3)$$

In this equation, B is the average magnetic field strength on the surface of the Moon,
<figure><img src="image_33.png" /><figcaption>Figure 4.12: Present day core-mantle boundary heat flow centered on the PKT. These results are for our preferred case (model T-0LB).</figcaption></figure>

$\mu_0$ is the permeability of free space, $\rho$ is the average core density, $F$ an efficiency factor, $q_0$ is the superadiabatic energy flux out of the core, $R_c$ and $R_p$ are the core and lunar radius, respectively, and $f$ is a prefactor expressing the efficiency of ohmic dissipation and the fraction of the field that is in the dipole term exterior to the core. Following Le Bars et al. (2012), we use $\rho = 7500 \text{ kg m}^{-3}$ and $f = 0.13$. $F$ is calculated from Christensen et al. (2008) and is 0.03 when the heat flux does not depend on radius within the core and when there is no inner core present. We note that this parameter is only weakly dependent upon these assumptions. Depending on the assumed initial temperature difference between the mantle and core, the thermal dynamo can last up to about 100 million years and produce a surface magnetic field of about $0.5 \ \mu\text{T}$. Given that $F$ and $f$ are both uncertain, the actual field could perhaps be up to about 10 times larger. Paleomagnetic studies of the oldest lunar samples imply surface fields that are at least $1 \ \mu\text{T}$ 4.2 Ga ago (Garrick-Bethell et al., 2009), and perhaps higher than $12 \ \mu\text{T}$ between 3.7 and 3.5 Ga (Shea et al., 2012; Suavet et al., 2013). Even though this thermal dynamo is shortlived and can not explain all of the paleomagnetic data, it could perhaps explain the magnetization of crustal materials early in lunar history, such as those associated with the South Pole-Aitken impact basin (Wieczorek et al., 2012). Other models that might explain the younger magnetic signatures include mechanical stirring at the core-mantle boundary from either precession (Dwyer et al., 2012) or changes in the lunar rotation rate from large impact events (Le Bars et al., 2012).

Even if the temporal evolution of the core heat flux is similar among our various simulations and cannot explain a long-lasting lunar dynamo, chemical buoyancy induced by inner core growth might help power the lunar dynamo for a much longer
period of time. A recent seismic study suggests that the Moon might possess a solid inner core (Weber et al., 2011), but this remains subject to debate (Garcia et al., 2011). The liquidus of iron alloys at lunar core pressures is highly dependent on the light element content and can vary between 2050 K for pure iron to 1300 K at the Fe-FeS eutectic. Given that our present day core-mantle boundary temperature is about 1700 K, it is therefore possible that a portion of the core could have crystallized. Using the phase diagrams from Buono and Walker (2011), we find that solid inner core crystallization would have occurred if the sulfur content was less than about 6 wt%. To account for a 240 km inner core as predicted by Weber et al. (2011), our core temperatures imply a sulfur content of about 3 wt%.

## 4.5 Discussion

As demonstrated in the previous section, mantle convection has a non-negligible influence on the thermal history of the Moon. From Fig. (4.2), we see that the lower mantle is cooler today than it would be for a purely conductive case, and this implies a larger heat flow out of the core between 4 and 3.5 billion years ago (as can also be seen in Fig. 4.11). Similarly, Fig. (4.4) shows that the onset of convection leads to a peak in the production of magma at some time near 4 Ga that the conductive models do not have, due to the upwelling of hot material from depth. Convective models therefore produce more volcanism, even though present day temperatures are slightly lower than the conductive cases. However, the effects of convection on the gravity observations of Fig. (4.9) are small as the largest temperature differences occur in the lower mantle, which have only a small effect on the surface gravity.

In addition, viscosity has an important effect on mantle dynamics. A lower reference viscosity, appropriate for more wet rheologies, would lead to a more vigorous convection, with a faster onset of convection and more up and downwellings. But as initiation of convection is also dependent on the initial temperature profile, by varying the initial conditions, it is always possible to find scenarios where the initial peak of melting occurs between about 4 and 3 Ga ago. A lower reference viscosity would lead to somewhat lower gravity anomalies as the mantle would be more homogeneous at the present time.

Although many of the general features of our thermal model in section 3 are robust to changes in the model parameters (such as the timing and duration of mare volcanism, the spatial distribution of mare volcanism, and dynamo generation), we note that the total volume of generated magma is highly sensitive to the assumed mantle solidus. First, the depletion of the mantle by melting would make it more refractory, causing its solidus temperature to increase. This effect was not considered in our model, nor is it in other lunar evolution models (e.g., Zhong et al., 2000; Hess and Parmentier, 2001; Spohn et al., 2001; Stegman et al., 2003; Ziethe et al., 2009), as the composition of the lunar mantle is not well known. Lacking better constraints, the entire mantle was set to a uniform peridotitic composition.
Nevertheless, we acknowledge that magma ocean crystallization could have led to a mantle that was zoned in composition. The first magma ocean cumulates would be rich in magnesium, but as crystallization progressed, the cumulates would become increasingly iron rich. To complicate matters, such a sequential cumulate pile is expected to be gravitationally unstable, and a mantle overturn could have brought buoyant nearly pure magnesian-rich dunite to the upper mantle (Hess and Parmen-dier, 1995; Elkins-Tanton et al., 2011). If this were to occur, the solidus of the upper mantle would likely be higher than that used in our study, and this could dramatically decrease the amount of magma generated by our model.

Heat source partitioning upon melting was also not considered and would help reduce the total amount of melt as a large part of the heat sources would be extracted from the mantle after reaching small melt fractions. This aspect of the melting process has been recognized to be important (Shearer et al., 2006), but was neglected in this study for the sake of simplicity. Kirk and Stevenson (1989) found that partitioning of heat sources in the melt reduces the total amount of melt produced in a given model. In addition, a stable mantle stratification, resulting from the gravitational overturn of the crystallized magma ocean cumulates would dramatically decrease convective motions within the mantle. If the density gradient with depth is great enough, the subsequent thermal evolution would be almost entirely due to the process of heat conduction and thus no decompressional melting would occur.

Our models are in general able to match qualitatively the duration of lunar volcanism, the asymmetry in nearside and farside magma production, the location of the youngest eruptions, and the present day heat flow measured at the Apollo 15 and 17 landing sites. Nevertheless, we acknowledge that there are two model predictions that are difficult to reconcile. First, our favored models generate about 10 time more magma than is predicted to exist. Second, our models predict several kilometers of dynamic topography within the PKT that is not observed. Regardless, these shortcomings are probably not critical, as there were several model parameters that could not be explored in this study. For example, changing the thickness of the KREEP layer or its distribution within the crust and mantle would help obtaining a better fit of these constraints. In addition, the mantle solidus could be modified to account for a post mantle overturn scenario. Finally, crustal thickness variations, which might be expected to occur during the formation of the PKT were not considered and would directly influence the predicted surface uplift.

## 4.6 Conclusions & perspectives

A 3D thermochemical convection code was used to study the influence of an enrichment in heat sources in the Procellarum KREEP Terrane on the nearside of the Moon. Our simulations show that the enhancement in heat production in this region partially melts the underlying mantle to depths of about 600 km over a duration of about 3 to 3.5 billion years. In contrast, the lower heat production on the
farside, combined with the initiation of mantle covection, give rise to farside magma
production rates that are 10 times smaller than the nearside, and which lasts for
only about 0.5 billion years. These results are consistent with the observed distri-
bution of the mare basalts, and well as the duration of mare volcanism as inferred
from lunar samples and crater counting techniques.

Our results imply that a large temperature anomaly is preserved at the present
day below the PKT. This result has immediate implications for the interpretation
of electromagnetic sounding data (Grimm, 2013), the modeling of tidal deformation
Love numbers (Zhong et al., 2012), the concentration of tidal stresses, as well as the
interpretation of the Apollo seismic data. This thermal anomaly predicts uplift of
the surface in the PKT, and when combined with the compositional depletion of the
mantle and mare basalts, a small gravity anomaly as well. If the lithosphere were
weak, the different gravity contributions would compensate to an almost zero gravity
anomaly with a surface uplift of a few kilometers. For a rigid lithosphere, no uplift
is generated, but the negative gravity anomaly would be larger. For the strengthless
lithosphere interpretation, this discrepancy could be resolved if the crust in the PKT
were 10-20 km thinner than the surroundings, or if the nearside mantle were denser
than average. Both possibilities are plausible given the unique manner by which
this province formed.

Finally, using the predicted core temperatures and core heat flux from our model,
we show that the Moon could have powered a dynamo with a surface field strength
of about 1 µT for only the first 200 million years of its evolution. This could perhaps
explain some of the oldest crustal magnetic anomalies, but not the paleomagnetic
results that require the dynamo to have lasted until about 3.56 billion years ago
(Suavet et al., 2013). Nevertheless, if the sulfur content of the lunar core is less
than 6 wt.%, our thermal history implies that some portion of the core could would
have crystallized to form a solid inner core. Not only is this consistent with recent
seismic data, but it also offers the possibility of driving a dynamo based on core
crystallization for a much longer period of time.
# Part II

## Core evolution and paleomagnetism
<figure><img src="image_34.png" /></figure>
# The lunar magnetic history

## Contents

<table><tr><td>5.1 (Paleo) Magnetic field measurements</td><td>59</td></tr><tr><td>5.2 Past and current models of magnetic history</td><td>62</td></tr><tr><td>5.3 Existence of a core</td><td>64</td></tr><tr><td>5.4 Our approach</td><td>66</td></tr></table>

## 5.1 (Paleo) Magnetic field measurements

Before the Apollo era, the Moon was thought to be a cold remnant from Solar System formation and as such, was not expected to possess any sign of interior activity. It was especially interesting for the Apollo/Luna missions to find signs of both melting and a surface magnetic field. Magnetometers were brought to the lunar surface by four Apollo missions to measure both spatial and temporal variations. Measurement of the temporal variations over several months permits the separation of the external and internal contributions. Measurements at the surface showed that the field varies strongly over short distances, with values going from 43 to 103 nT over a kilometer next to Apollo 14 landing site. The presence of these variations coupled with orbital measurements suggest the presence of strong localized sources, and its presence presented a challenge due to the assumed absence of activity of the Moon (e.g. *Dyal et al.*, 1970).

The better coverage and resolution associated with the more recent data from Lunar Prospector (e.g., *Purucker*, 2008) better constrains the origin of that field. Figure (5.1) shows a global projection of the surface magnetic field centered on the near and farside. The field distribution is uneven and not uniformly correlated with surface features. Specific anomalies occur in the center of basins and could originate from melt sheet cooling in the presence of a magnetic field. Other anomalies are not linked to any morphological expression; some anomalies are antipodal to basin impacts while others are not. This suggests a complex history for the evolution of the lunar magnetic field.

The lunar magnetic field being of remanent origin, the main question is what caused that remanent magnetisation. The two possibilities for the Moon are thermoremanent magnetisation (TRM), where a sample cooled below its Curie temperature in an external magnetic field, and shock-induced remanent magnetisation
<figure><img src="image_35.png" /><figcaption>Figure 5.1: Surface magnetic field strength from the sequential Lunar Prospector model from *Purucker and Nicholas* (2010) evaluated at 30 km above the surface. The figure is adapted from *Wieczorek et al.* (2012). The projection is Lambert azimuthal equal-area projection.</figcaption></figure>

(SRM), where magnetisation is induced by the passage of a shock wave after an impact occurs. Samples with traces of shock can possess often have a more complex history and therefore paleomagnetic studies usually use unshocked material in order to have a better understanding of its thermal history and therefore better constraints on its magnetic age.

Most lunar paleomagnetic studies are carried out either with the saturation Isothermal Remanent Magnetisation (sIRM), the Thellier-Thellier method or the an-hysteretic remanent magnetisation method (ARM) (Collinson, 1993). The Thellier-Thellier is an absolute method based on the independence and additivity principle of partial thermoremanent magnetisation (pTRM). By heating up a sample by small steps and monitoring the remaining TRM, one can clean magnetic overprints and recover the intensity of the stable paleomagnetic field (Dunlop, 2011). The problem with this method applied on lunar samples is that phase changes can occur and therefore destroy some of the information. The sIRM and ARM method are both isothermal and therefore preserves the sample. In the sIRM method, the paleointensity is recovered by assuming a poorly constrained linear relation between the NRM and the saturation magnetisation obtained by applying a very large magnetic field in the lab. In the second method, the paleomagnetic field is recovered by ap-
plying an alternating field of declining amplitude. The problem with these methods
is that the interpretation of the results are based on several assumptions and are
thus subject to discussion (Lawrence et al., 2008; Tikoo et al., 2012, e.g.,).

Initial paleomagnetic studies of this remanent field suggested that it was induced
by a global field 3.9 Ga ago that lasted a few hundred million years (e.g., *Cisowski
et al.*, 1983) with some recent samples (< 1.5 Ga) also presenting a much weaker field.
Figure (5.2) shows a summary of early paleomagnetic studies (in red). In addition,
the field strength, up to about 50 µT, similar to Earth observations, is very high for
such a small body. However, since lunar rocks are poor magnetic recorders and their
history are often uncertain, these estimates should be used with caution (Lawrence
et al., 2008). Since then, new estimates have been conducted on carefully selected
samples. They have extended the duration of the global magnetic field up to about
one billion years, while revising its strength to lower values (Garrick-Bethell et al.,
2009; Shea et al., 2012; Suavet et al., 2013). New values are pictured in Figure (5.2)
in green.

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

Figure 5.2: Ancient (red) and recent (green) paleomagnetic studies. The green arrows stress the fact that quoted values are minimal magnetizing fields. The orange region shows the span of the original estimated magnetic era. The dashed blue line represents a typical value of Earth's magnetic field.

The current view suggests that a field of about at least several µT was present
at the surface of the Moon between 4.2 and 3.2 Ga, with a peak of more than 10
µT 3.7 Ga ago (Garrick-Bethell et al., 2009; Shea et al., 2012; Tikoo et al., 2012).
It is interesting to note that samples older than 4.2 billion years have not been
analysed and therefore an even older dynamo cannot be formally discarded. On
the other hand, *Tikoo et al.* (2012) showed that experimental uncertainties could have been the cause of the observed magnetic field on recent (< 1.5 Ga) samples and thus it is probable that the lunar dynamo did shut down between about 3 Ga and 1.5 Ga ago. Though several types of models have been proposed to try and explain these observations, the strength and timing of the lunar magnetic field has not convincingly been explained to date.

In summary, even if paleomagnetic studies have made a lot of progress in the last decade, it should not be unexpected that our view of lunar magnetism could change in the next decades. What is accepted today, is that there was a global field of *at least* a few micro teslas between *at least* 4.2 and 3.2 Ga and that no active field is present today. Several processes have been proposed to explain specific observations but no consensus exists to date. In the next section, I present some of these models and what they manage to explain.

## 5.2 Past and current models of magnetic history

Models of the origin of the lunar magnetic field can be separated in three classes: those of internal origin, external origin and transient, impact associated fields. The external field models argue that either the interplanetary field or the geomagnetic one could be the source of lunar magnetism. Internal field models include a core dynamo, but also more exotic cases such as bulk magnetisation of the interior or localised dynamo field from two lava pools of different temperature acting as a thermoelectric generator. Finally impact associated field models include both cases where the field is generated by the impact process and cases where the impact acts as a field magnifier only. These models are described in some details in *Collinson* (1993) and we describe a few of them in slightly more details here.

A typical way of explaining long lasting planetary magnetic fields is to invoke a core dynamo (*Stevenson et al.*, 1983). Planetary magnetic fields can be generated by convection in a liquid conducting core. In the case of a thermal dynamo, convection is triggered by standard Rayleigh-Benard instabilities, i.e. when more heat is extracted from the core than can be conducted along the adiabat. The strength of the magnetic field is governed by the vigor of convection. However, the Moon is a small body and mantle thermal evolution studies which looked at this approach found that heat extraction from the core is not large enough to produce a magnetic field for more than a few hundred million years, between 4.5 and about 4.2 Ga (*Konrad and Spohn*, 1997; *Laneuville et al.*, 2013). A more exotic model by *Stegman et al.* (2003) investigated the influence of an initially stratified mantle that could arise after the crystallised magma ocean gravitationally readjusted, but arrived to the same conclusion that a thermal dynamo could not sustain a magnetic field for more than a few 100 million years.

Another approach investigates dynamos driven by differential rotation between mantle and core induced by either precession of the lunar spin axis (*Dwyer et al.*, 2012) or impacts (*Le Bars et al.*, 2012). In both cases, differential rotation induces
large-scale flow in the core, which could have powered a lunar dynamo (e.g., *Le Bars et al.*, 2010). We note here that it is unclear if the efficiency and magnetic scaling of such dynamos is similar to standard thermo/chemical ones. The precession induced magnetic field is predicted to last for more than a billion years with intensities of about 1 µT, whereas the impact scenario predicts magnetic anomalies associated with large impact basins. Magnetic anomalies are indeed observed in the interior of some basins and could be due to the melt sheet cooling in the presence of a dynamo. The timescale of such impact-induced dynamo is about $10^3$ years, consistent with the melt sheet cooling timescale. Both could have therefore explained different features, but the predicted surface magnetic field is still an order of magnitude too small compared to paleomagnetic studies, which is the case with all lunar paleofields to date.

An alternative hypothesis is that the origin of this field was transient and linked to impact processes. In such a scenario, the partially ionized vapor cloud ejecta will interact strongly with any ambient field (either from internal or external source) due to its high electrical conductivity, leading to transient field amplification. The expansion of that vapor cloud will exclude the external magnetic field from the impact zone and compress it as it converges to the antipodal point. This effect could explain the strength of the observed field for anomalies located at the antipodes of impact craters. The idea dates back to *Hide* (1972) which suggested that the magnetic effects of high-velocity impacts cannot be ruled out. Several studies have since analysed the idea and *Hood and Artemieva* (2008) showed that the ejecta cloud released after impact could indeed converge to the antipode and that the time of maximum antipodal field amplification corresponds to the arrival of impact shock stresses able to produce a stable shock remanent magnetisation. It thus appears that this process can explain some particular features, but could be hard to reconcile with the observed global nature of the field, especially thermo-remanent anomalies.

Another way to explain the strength of the magnetic anomalies is to use highly magnetic extralunar materials. *Wieczorek et al.* (2012) showed that the ejecta blanket induced by an oblique impact corresponding to the South Pole-Aitken basin could well have distributed material in a large part of the farside. In such a scenario, the deposit of a few kilometers of high magnetic susceptibility material in an ambient field of about 10 µT (without any assumption about the origin of that field) would be enough to explain the observed 10 nT field at 30 km altitude. However this model may have more problems explaining nearside anomalies. Extralunar chondrites from impacts have a magnetic susceptibility about 100 times larger than mare basalts, and thus a much thinner layer is required to explain the observed remanent field. Note here that no single process is required to explain all the observation and a simple combination of a few processes could be the solution.
## 5.3 Existence of a core

A typical way to explain a global magnetic field is by dynamo action in a fluid core. Even if the existence of a core was not believed initially, several lines of argument now exist to suggest its existence.

### Moment of inertia

The moment of inertia is one way to constrain the mass distribution inside a planet. It does not uniquely determine the distribution but describes to what extent mass is concentrated towards the center of the body. The lunar moment of inertia is very close to that of a homogeneous sphere and thus its core should be small. A wide range of interior models have been designed from cases where the Moon has no core at all to cases with a 460 km core (Hood, 1986). Considered alone, this measurement does not require nor exclude the presence of a core (e.g., Ferrari et al., 1980). However, when used with other datasets, such as seismology, the moment of inertia can provide additional constraints.

### Induced magnetic dipole

Another historical method to determine the size of the lunar core is through measurement of the dipole induced as the Moon traverses the geomagnetic tail lobes. The almost spatially uniform external field diffuses into the lunar mantle and after a few hours induces current at the surface of a highly conducting core. The induced dipole moment is then directly relatable to the core radius and conductivity (Goldstein and Phillips, 1976). Using that method on the Apollo dataset, Russell et al. (1982) found that a highly conducting core of about 440 km was present. Later, Hood et al. (1999) found that Lunar Prospector measurements were consistent with a core between 250 and 540 km. The uncertainty associated with such measurements is large and these estimates are significant only when coupled to other constraints.

### Laser ranging

Study of laser ranges between Earth and lunar retroreflectors on the Moon can be used to determine the time-varying rotation of the Moon. This can then be used to infer the elastic properties of the lunar interior. Dissipation can originate either from solid body tides raised by the Earth and Sun, or from an at least partially liquid core with a rotation distinct from the solid body. These two contributions can be separated (Williams et al., 2001) and studies show that the observed signal is best explained by a sum of these two contributions and thus, that the lunar core is at least partially liquid. Williams et al. (2001) predicted a maximum radius for the core of 352 km for a pure iron case and 374 km for a Fe-FeS eutectic composition.
## Seismology

In theory, it is also possible to constrain the size of the core from seismologic measurements. This was made possible by the Apollo network, but due to the limited amount of data, it took more than 40 years to obtain the first direct estimates. Initial studies used an indirect approach, computing a compositional model for the crust and mantle and inferring core properties with an additional constraint from the moment of inertia. The range of core sizes is about 200 to 530 km, depending on the assumed composition (see Wieczorek et al., 2006, for a review). The first direct detection by Weber et al. (2011) predicted a core radius of about 330 km and an inner core radius of 240 km. At the same time, Garcia et al. (2011) estimated a core of about 380 km, while not looking for the existence of an inner core.

## Geochemical constraints

During differentiation, the siderophile elements are extracted into the core according to their partitioning coefficient. Thus their post-differentiation mantle composition, assuming the initial bulk composition is known, provides constraints on the amount of iron that differentiated and consequently, on the size of the core. Righter (2002) showed that depending on the Moon-forming giant impact conditions, a core of up to 2wt.% (~370 km radius) is consistent with siderophile elements concentrations.

## Thermochemical constraints

Inner core size, composition and thermal state are coupled; knowing two of them can provide constraints on the third one. Thermal evolution models predict present day core-mantle boundary temperatures in the range 1300 to 1750 K (see Shearer et al., 2006, for a review). The existence of a liquid outer core today implies that the core is not pure iron, because it would then be totally solid at such temperatures. For example, to match their 240 km inner core radius estimate, Weber et al. (2011) constrained the core sulfur content to less than about 6 wt.%.

In addition, coupling seismic data with other observables such as the moment of inertia allows for the inversion of the lunar composition and thermal state (e.g. Khan et al., 2007). In these inversions, the estimated lunar core density is lower than that of pure iron and predictions then can be made on light element content. The average core density predicted by either Garcia et al. (2011) or Khan et al. (2007) (5.2 and 5.7 g/cm³, respectively) corresponds to a core with about 10 or more wt.% S. However, as acknowledged by the authors, when taking the error bars into account, composition constraints from these density estimates should be used with care.

Similar debates and similar techniques exist on Earth, with the additional complexity of high pressures. No consensus has been found on what is the alloying element in the core, and the view that there is probably more than one slowly becomes more accepted (e.g. Hillgren et al., 2000). The Moon has probably been
formed after a giant impact with the proto-Earth and this adds another layer of
uncertainty. However, two points are surely accepted today: the Moon possess a
300 to 400 km core and that core is at least partially molten. Figure (5.3) shows
the possible interior structure proposed by Wieczorek et al. (2006).

<figure><img src="image_37.png" /><figcaption>Figure 5.3: Proposed lunar interior structure from Wieczorek et al. (2006).</figcaption></figure>

**5.4 Our approach**

In this part, we study the influence of inner core growth on the lunar dynamo and propose a scenario which explains the duration of the Moon's internal magnetic field. Such an approach is standard in the case of the Earth (e.g., Aubert et al., 2009) or asteroids (Nimmo, 2009), but has yet to be considered on the Moon. As proposed for other planetary bodies, depending on the light element content, the relative slopes of the adiabat and liquidus can change and crystallisation can start at the core-mantle boundary instead of starting at the center of the body, resulting in the so-called iron-snow regime (Hauck et al., 2006; Williams, 2009). Williams (2009) showed that for a wide range of initial sulfur content, crystallisation should start at the center of the Moon. We therefore only consider cases with initial bottom-up crystallisation.

In that case, compositional buoyancy due to release of light elements at the
inner core boundary helps sustain convection in the outer core against dissipation,
even when the heat extracted by the mantle is smaller than the core’s adiabat.
However, a transition to the iron-snow regime could still occur later in history due
to sulfur enrichment in the outer core from the crystallisation process. This effect
is not considered here, but will be discussed in the results section. Our recent lunar
thermal evolution study (Laneuville et al., 2013) suggested that core crystallization
could indeed occur on the Moon, depending on the initial light element content and
we therefore investigate this process into more details in the next chapters.

We will start by presenting our thermodynamic approach in Chapter 6, then
the influence of different alloying elements and thermodynamic properties will be
described in Chapter 7, along some discussion about other effects.
<figure><img src="image_38.png" /></figure>
# Evolution models of the lunar dynamo

## Contents

<table><tr><td>6.1</td><td>Lunar thermal evolution</td><td>69</td></tr><tr><td>6.2</td><td>Energy budget of the core</td><td>70</td></tr><tr><td>6.3</td><td>Entropy budget</td><td>73</td></tr><tr><td>6.4</td><td>Magnetic field scaling</td><td>75</td></tr><tr><td>6.5</td><td>Model parameters</td><td>76</td></tr></table>

In order to estimate the surface magnetic field, the power available to sustain the dynamo has to be estimated (Christensen and Aubert, 2006; Aubert et al., 2009). This power is directly linked to the sum of the thermal and chemical buoyancy fluxes within the core. As core and mantle evolution are coupled, we need to make assumptions on the mantle evolution. They are described in Section 6.1. The temporal evolution of the inner core can then be obtained through the core energy budget and relevant phase diagrams, described in Section 6.2. We use the entropy budget to determine the part of the power available to dynamo action in Section 6.3. Finally, the magnetic field scaling law is described in Section 6.4 and a discussion on the model's parameters can be found in Section 6.5.

## 6.1 Lunar thermal evolution

Laneuville et al. (2013) proposed thermochemical evolution models for the Moon considering the asymmetric distribution of heat sources at its surface (Jolliff et al., 2000). We use here the same assumptions and initial conditions and couple the mantle to our core evolution model. The influence of this feedback between the inner core growth and mantle has a non-negligible influence on the core evolution and should therefore not be underestimated in future studies. For example for a simulation predicting an inner core 50% of the core radius at present, neglecting latent heat release would have underestimated core-mantle boundary heat flow by a factor of about 2 during most of the thermal history. In addition, the assumption made in Laneuville et al. (2013) of localized heat sources does not play a major role on inner core formation, as long as they are still concentrated in the upper layers of the Moon. Details about these effects will be discussed in Chapter 7.
<figure><img src="image_39.png" /><figcaption>Figure 6.1: Preferred thermochemical evolution from Chapter 4. (left) Thermal state after 1 billion years evolution. (right) Core-mantle boundary heat flow as a function of time (without including inner core growth). The dashed line is the limit of a heat-flow equal to zero. The core adiabatic heat flux is expected to be about 5 mw/m².</figcaption></figure>

The only parameter that we will vary that was not studied in *Laneuville et al.* (2013) is the initial temperature difference between core and mantle. Assuming that the gravitational energy difference upon core differentiation from the silicate mantle heats the core preferentially, an initial temperature excess of about 700 K would be observed (e.g., *Solomon*, 1979). Heat partitioning between iron and silicates is not well constrained, we thus consider extreme cases to assess the effect. Following *Konrad and Spohn* (1997), we take a lower value of 200 K and keep the same high value as in *Laneuville et al.* (2013), i.e. 400 K. Another way to estimate this value is to assume that the entire Moon were molten and that the silicate mantle crystallized quickly. The temperature difference between the core and mantle would then be equal to the temperature difference between liquidus and solidus at the core mantle boundary (about 200 K). Figure (6.1) shows characteristics of our chosen mantle model from Chapter 4 (i.e., without including inner core growth).

## 6.2 Energy budget of the core

The core energy budget is dominated by secular cooling in the absence of core crystallisation. When the inner core starts to freeze, the latent heat term becomes non-negligible and a compositional buoyancy term is also introduced. We assume that no radioactive heat sources are present in the core. The energy budget can then be written as

$$Q_{CMB} = Q_{S} + Q_{L} + Q_{g}, \qquad (6.1)$$

where $Q_{CMB}$ is the average heat flux out of the core, $Q_{S}$, $Q_{L}$ and $Q_{g}$ are the part of the energy budget due to secular cooling of the core, latent heat release upon
inner core growth and compositional buoyancy due to light elements concentration in the remaining fluid, respectively. Following Nimmo (2007), the exact expression of these contributions is

$$Q_S = \int \rho c_p \frac{dT_c}{dt} dV, \qquad (6.2)$$

$$Q_L = 4\pi r_i^2 L_c \rho^{ic} \frac{dr_i}{dt}, \qquad (6.3)$$

$$Q_g = \int \rho \psi \alpha_c \frac{D\chi}{Dt} dV, \qquad (6.4)$$

where $\rho$ is the average core density, $\rho^{ic}$ the density of the solid inner core, $c_p$ the specific heat, $T_c$ the core-mantle boundary temperature, $t$ the time, $V$ the core volume, $r_i$ the inner core radius, $L_c$ the latent heat of crystallisation, $T_i$ the temperature at the inner core boundary, $\psi$ the gravitational potential, $\alpha_c$ is a dimensionless coefficient which indicates the sensitivity of the core density to the presence of light element $\chi$ and $D\chi/Dt$ is the total derivative of $\chi$ (including advection).

$$\alpha_c = \frac{\Delta\rho_c}{\rho^{ic}\Delta\chi} \qquad (6.5)$$

where $\Delta\rho_c$ is the density difference across the inner core boundary due to the compositional difference $\Delta\chi$ and $\rho^{ic}$ is the inner core density.

We now make a few assumptions to simplify equations (6.2)-(6.4). For an isentropic outer core with homogeneous composition, the temperature profile is adiabatic and can be written as in (Labrosse, 2003)

$$T_a(r, T_c) = T_c \exp \left[ -\frac{r^2 - r_c^2}{D^2} \right], \qquad (6.6)$$

where $r$ is the radius, $r_c$ is the core radius and $D$ the scale height

$$D = \left( \frac{3c_p}{2\pi\alpha\rho^{oc}G} \right)^{1/2}, \qquad (6.7)$$

where $G$ is the gravitational constant and $\alpha$ the thermal expansivity. The outer core density $\rho^{oc}$ is a linear function of light element content $\chi$, in wt.%

$$\rho^{oc}(\chi) = \rho_0^{oc} - a\chi, \qquad (6.8)$$

where $a$ is an alloy-dependent constant of proportionality. Assuming an adiabatic profile for the outer core temperature and an isotherm in the inner core implies that the complete core state is determined by one temperature only, which we choose to be the core-mantle boundary temperature. The inner core growth is therefore directly proportional to the core-mantle boundary temperature evolution and can be obtained through the equality of melting and adiabatic temperature at the inner-core boundary. After the core cooled by $\delta T_c$, the adiabat and liquidus temperature at $r = r_i + \delta r_i$ have to be equal. This can be written as

$$T_m(r_i + \delta r_i) = T_a(r_i + \delta r_i, T_c - \delta T_c). \qquad (6.9)$$
After a first order Taylor expansion of both sides of the equation and assuming hydrostatic equilibrium (i.e. dP=ρg(z)dz), we obtain

$$ \frac{dr_i}{dt} = \frac{1}{(dT_m/dP - dT_a/dP)\rho^{ic}g} \frac{T_i}{T_c} \frac{dT_c}{dt}, \qquad (6.10) $$

where $dT_m/dP$ and $dT_a/dP$ are the change of melting temperature with pressure and adiabatic gradient, respectively. Assuming $T_i \sim T_c$ (using Eq. (6.6), this leads to an error of less than 1%) and $dT_a/dr = 2rT/D^2$ (consistent with the order 1 Taylor expansion), Equation (6.10) can be further simplified to

$$ \frac{dr_i}{dt} = \frac{D^2}{2T_c f r_c (\Delta - 1)} \frac{dT_c}{dt}, \qquad (6.11) $$

where $\Delta = (dT_m/dP)/(dT_a/dP)$ and $f = r_i/r_c$. Similarly, the rate of light element release is written

$$ \frac{D\chi}{Dt} = \frac{4\pi r_i^2 \rho^{ic} \chi}{M_{oc}} \frac{dr_i}{dt}. \qquad (6.12) $$

The mass fraction of light element $\chi$ is assumed homogeneous in the outer core, but evolves with $r_i$ according to mass balance between inner and outer core

$$ \chi(r_i) = \frac{\chi_0}{1 - f^3}. \qquad (6.13) $$

where $\chi_0$ is the initial light element mass fraction in the core. Note this is exact only in the case where pure Fe is crystallised, for other cases, a complete mass balance has to be computed. However an Fe-FeS core is the case where the enrichment effect is the strongest because no sulfur is partitioned in the solid inner core. Still the effect is relatively small, as an inner core size of about 80% of $R_c$ is required before doubling the initial sulfur concentration in the outer core. Following Nimmo (2009), we truncate the Taylor expansion of (6.6) to the second order in ($r_c/D$) and write the secular term as

$$ Q_S = M_c c_p \left( 1 + \frac{2r_c^2}{5D^2} \right) \frac{dT_c}{dt}. \qquad (6.14) $$

The latent heat term is then

$$ Q_L = \frac{3}{2} M_c \frac{f L_H D^2}{T_c r_c^2} \frac{1}{\Delta - 1} \frac{dT_c}{dt}. \qquad (6.15) $$

Finally, the compositional buoyancy contribution can be written as

$$ Q_g = 3\pi G\rho M_c F \frac{\Delta\rho_c}{\rho^{ic}} \frac{1}{\Delta - 1} \frac{D^2}{T_c} \frac{dT_c}{dt}, \qquad (6.16) $$

where $F = f(1/5 + 2f^5/15 - f^2/3)/(1 - f^3)$, $\Delta\rho/\rho^{ic}$ is the density change upon inner core solidification with $\Delta\rho_c = \rho^{oc} - \rho^{ic}$. The ratio ($r_c/D$) is about 0.1 thus a second order truncation of Equation (6.6) creates an error of less than 0.01 K. For a derivation of these terms, see Appendix B. To derive this equation, we assumed
that the density jump at the inner-core boundary is just a small perturbation of the gravity potential and we neglected it. The inner core density is a linear function of light element content

$$ \rho^{ic}(\chi) = \rho_0^{ic} - b\chi, \qquad (6.17) $$

where b is an alloy-dependent constant and can be equal to zero (in the Fe-FeS case for example). As can be seen from Eq. (6.14)-(6.16), each term on the right-hand side of (6.1) is proportional to $dT_c/dt$. We can therefore use $Q_{CMB}$ to compute $dT_c/dt$ from which $dr_i/dt$ follows. Figure (6.2) is an example of the different contributions to the energy budget for a typical case. It is clear that, as soon as inner core crystallization occurs (after about 500 Ma in that case), release of latent heat accounts for the most part of the heat extracted from the core and a strong reduction in cooling rate is to be expected.

<figure><img src="image_40.png" /><figcaption>Figure 6.2: Different contributions to the energy budget for a typical case (case dT300_X6 in Table 7.1). $Q_s$, $Q_g$ and $Q_L$ are defined in the text and are the secular, compositional buoyancy and latent heat term, respectively. In that case, core crystallisation occurs after about 500 million years.</figcaption></figure>

## 6.3 Entropy budget

Ohmic dissipation generated by dynamo activity did not appear in the energy budget because buoyancy forces generate kinetic and magnetic energy which are directly converted to heat by viscous and ohmic dissipation. It is just a conversion of magnetic energy to heat within the core and does not enter the global energy budget.

However entropy is a measure of the available energy in a closed system and, as dissipation is a non-reversible process, it reduces the available energy with time (i.e. it is a source of entropy). Investigation of the entropy budget can therefore be used to estimate the amount of dissipation possible by dynamo action. It is written as

$$ E_{\phi} + E_{k} = E_{S} + E_{g} + E_{L}, \qquad (6.18) $$
where $E_\phi$ is the dissipated entropy, $E_S$ the entropy production due to secular cooling, $E_g$ the entropy production due to light element release at the inner core boundary, $E_L$ the entropy production due to latent heat release at the inner core boundary and $E_k$ is the entropy sink due to conduction along the adiabat. Note here that other entropy producing sources, such as friction at the core mantle boundary between the liquid core and the mantle, could largely increase the available energy to sustain the dynamo. However the additivity of the processes is unclear and this is out of the scope of this study. For a stable dynamo to occur, enough entropy has to be produced to balance ohmic dissipation, i.e., $E_\phi$ has to be greater than zero. The derivation of these terms can be found in Appendix B and can be written as

$$E_S = \int \rho c_p \left( \frac{1}{T_c} - \frac{1}{T} \right) \frac{dT_c}{dt} dV, \qquad (6.19)$$

$$E_L = Q_L \left( \frac{1}{T_c} - \frac{1}{T_i} \right), \qquad (6.20)$$

$$E_g = \frac{Q_g}{T_c}, \qquad (6.21)$$

$$E_k = \int k \left( \frac{\nabla T}{T} \right)^2 dV. \qquad (6.22)$$

Before inner core nucleation, $E_g$ and $E_L$ are both zero, but as soon as it starts, they dominate the budget. This is due to the thermodynamic efficiency linked to entropy production. Entropy is produced or consumed at a given temperature, which will then define their thermodynamic efficiency. Thus buoyancy sources that are distributed in the entire outer core (such as secular cooling) have a smaller average operating temperature than buoyancy sources released at the inner core boundary only. This can be seen in the efficiency factors in Equations (6.19)-(6.21). Entropy production linked to inner core growth is therefore more efficient than that linked with secular cooling, even though their contribution to the energy budget is modest. Using the same assumptions as for the energy budget, these terms can be written as

$$E_S = M_c c_p \frac{2r_c^2}{5T_c D^2} \frac{dT_c}{dt}, \qquad (6.23)$$

$$E_g = \frac{Q_g}{T_c}, \qquad (6.24)$$

$$E_L = \frac{3}{2} M_c \frac{f(1-f^2)L_H}{T_c^2} \frac{1}{\Delta-1} \frac{dT_c}{dt}, \qquad (6.25)$$

$$E_k = M_c \frac{12k}{5\rho D^4} r_c^2, \qquad (6.26)$$

where $k$ is the core thermal conductivity. Assuming most of the entropy is produced at the inner core boundary, the volumetric power available to drive the dynamo is the power that balances $E_\phi$, as written below:

$$p = \frac{E_\phi T_i}{V_{oc}}, \qquad (6.27)$$
where $T_i$ is the temperature at the inner core boundary and $V_{oc}$ the volume of the outer core. Figure (6.3) shows the different contributions to the entropy budget for a typical case. Even if its contribution to the energy budget was small, the entropy production due to release of light elements plays a large role in the entropy budget due to its intrinsic thermodynamic efficiency.

<figure><img src="image_41.png" /><figcaption>Figure 6.3: Different contributions to the entropy budget for a typical case (case dT300_X6 in Table 7.1). $E_s$, $E_g$ and $E_L$ are defined in the text and are the secular, gravity and latent heat term, respectively. $E_k$ is the entropy required to maintain the adiabatic gradient. In that case, core crystallisation occurs after about 500 million years.</figcaption></figure>

## 6.4 Magnetic field scaling

The magnetic Reynolds number is another constraint on the possibility of a dynamo within a planetary core. It is defined as $Re_m = vL/\eta$, where $v$ is the characteristic velocity, $L$ the characteristic length (taken to be the thickness of the convective shell) and $\eta$ the magnetic diffusivity. Following Nimmo (2009), we estimate the velocity induced by both thermal and chemical buoyancy by

$$v \approx 0.85 r_c^{3/5} \Omega^{-1/5} \left( \frac{4\pi G \Delta\rho}{3} \frac{dr_i}{dt} \right)^{2/5}, \quad (6.28)$$

where $\Omega$ is the rotation rate and the rest have already been defined previously. Using the present day rotation rate $\Omega = 2\pi/24 \text{ h}^{-1}$ in order to obtain a lower bound on $v$ and $\eta = 2 \text{ m}^2 \text{ s}^{-1}$, $v \approx 1 \text{ mm/s}$. Thus $Re_m \approx 200$, which is above the critical value for dynamo activity of 40 (Christensen and Aubert, 2006). This confirms that, if the available energy for dissipation is positive, as defined in the previous section, a dynamo is likely to occur.
Once we obtain the volumetric energy available for the dynamo, we can use a magnetic field scaling law to estimate the surface magnetic field. Christensen et al. (2008) derived a scaling law to determine the magnetic field strength of planets and stars, whereas Aubert et al. (2009) concentrated on a scaling law to study the palaeo-evolution of the geodynamo. Both scaling are based on similar arguments but using a different framework. We chose to use the scaling law presented in Aubert et al. (2009) because it the energy available for magnetic dissipation is calculated independently. It is defined as

$$B = f_{eff}(\rho\mu_0^3 p^2 L^2)^{1/6} \left(\frac{r_c}{r_m}\right)^3, \qquad (6.29)$$

where B is the average magnetic field strength at the surface of the Moon, μ the permeability of space, p the volumetric power available for the dynamo as computed in the previous section, L is the outer core thickness, $r_m$ the radius of the Moon and $f_{eff}$ is an efficiency factor including the scaling law prefactor from Aubert et al. (2009) which is about 1 and the ratio of magnetic field inside the shell to the dipole field at the CMB. f is about 0.1 because only about 1/2 of the energy is in the poloidal part and only a fraction of that energy is in the dipole term (the rest being in higher multipole terms). The $r_c/r_m$ product is just a scaling factor to obtain the surface dipole field. Note that the most important variable to determine the strength of the magnetic dipole field is the ratio $r_c/r_m$. Higher multipoles are neglected in this study.

## 6.5 Model parameters

Several core properties can be varied to assess the effect of different alloying elements. Here we focus on the Fe-FeS system, for which we have the most data. In addition we consider the effect of Fe-C and Fe-Si, but with looser constraints. The interest in using these 3 systems is that they provide a wide range of behavior. In the Fe-FeS system, pure Fe is crystallising and therefore the effect of light elements enrichment in the outer core is maximal. In the Fe-Si system, the inner core crystallises with the same concentration as the liquid and therefore no enrichment is observed. The Fe-C system is an intermediate case. We use data collected in Buono and Walker (2011), Chabot et al. (2008) and Yang and Secco (1999) on the Fe-FeS, Fe-C and Fe-Si system, respectively. We note here that interesting behavior can occur if a more complex case is investigated, such as a ternary system, but the lack of data do not allow us to consider it.
Buono and Walker (2011) propose an order 4 parameterisation of the melting temperature of sulfur both in pressure and composition

$$ T_m(P, \chi) = (-2.4724P^4 + 28.025P^3 - 9.1404P^2 + 581.71P + 3394.8)\chi^4 \\ + (1.7978P^4 - 6.7881P^3 - 197.69P^2 - 271.69P - 8219.5)\chi^3 \\ + (-0.1702P^4 - 9.3959P^3 + 163.53P^2 - 319.35P + 5698.6)\chi^2 \\ + (-0.2308P^4 + 7.1P^3 - 64.118P^2 + 105.98P - 1621.9)\chi \\ + (0.2302P^4 - 5.3688P^3 + 38.124P^2 - 46.681P + 1813.8). \quad (6.30) $$

Due to the scarcity of available data, we stick to 1st order for the Fe-Si and Fe-C systems and model the melting temperature as

$$ T_m(r) = (1 - \alpha_l \chi)(T_{m0} + T_{m1}P(r)), \quad (6.31) $$

where $T_{mi}$ are constants, $\chi$ is the mass fraction of the light element, $\alpha_l$ is an alloy-dependent constant and $P$ is the pressure.

Fitting Equation (6.31) to data in Chabot et al. (2008); Yang and Secco (1999), we find that the pressure dependence is $T_{m0} = 1800$ K and $T_{m1} = 36 \ 10^{-9}$ K/Pa. In the case of Fe-Si, the liquidus is indistinguishable from the solidus and we find $\alpha_l = 1.5$. In the Fe-C system, the solidus is different and we find $\alpha_l = 12$ for the solidus and $\alpha_l = 5.6$ for the liquidus, while the pressure dependence remains the same as for Fe-Si. Figure (6.4) shows the phase diagrams we consider. The red curves are the one required by our models and that we linearized. In all cases, we consider the Fe-rich side of the eutectic only, which is at about 27 wt.% in the case of Fe-FeS.

The densities of the solid and liquid phases are also linearized as follows:

$$ \begin{aligned} \rho^{ic}(\chi) &= \rho_0^{ic} - b\chi, \\ \rho^{oc}(\chi) &= \rho_0^{oc} - a\chi. \end{aligned} $$

Using data in Terasaki et al. (2010), we assume that the density of the liquid phase does not really depend on which is the light element and a linear regression of his data provides $a = 4$. The data point presented in Terasaki et al. (2010) for sulfur is inconsistent with that assumption, but different measurements have shown that it may be due to an experimental issue (e.g., Nishida et al., 2008). The density of pure liquid Fe at 4 GPa is taken from the equation of state presented in Komabayashi and Fei (2010). We find $\rho_0^{oc} = 7$ g/cm$^3$. The same reference is used for solid Fe at 4 GPa and we obtain $\rho_0^{ic} = 7.5$ g/cm$^3$. Fe-C and Fe-Si solid densities are obtained using a simple molecular weight ratio which give $b = 5.4$ and $b = 2.2$, respectively. The linearization found for the Fe-Si system is checked against data from Zhang and Guyot (1999) and is found to be consistent. All these parameters are summarized in Table 6.1, and the other model parameters can be found in Table 6.2.

Using these parameters, one can obtain the predicted temperature at the inner-core boundary as a function of inner core radius for a given initial alloy content.
<figure><img src="image_42.png" /><figcaption>Figure 6.4: Phase diagram scheme for different alloying elements. Red curves are the ones that are linearized and will be used in the models. Left is Fe-C and right is Fe-Si. The Fe-FeS liquidi are from *Buono and Walker* (2011).</figcaption></figure>

Figure (6.5) shows that, due to the cubic dependence of volume on radius, outer core enrichment in the alloying element is negligible until the inner core radius reaches about 60% of the outer core. After that the light element content rises quickly and the outer core can transition to a region of iron snow or up to the eutectic point, where both Fe and FeS crystallise. This should be kept in mind while discussing the results.

The effect of the kind of alloying element that is used can be seen in Figure (6.6). In that figure, we chose initial alloy contents such that the inner core starts to crystallise at the same temperature in all three cases (i.e., for $T_c = 1870$ K),
<table><caption>Table 6.1: Coefficients of the melting and core density linearization. (l) and (s) correspond to the liquidus curve and solidus parameterization, respectively.</caption><thead><tr><th>Coefficient</th><th>Fe-FeS</th><th>Fe-C</th><th>Fe-Si</th><th>Units</th></tr></thead><tbody><tr><td>a</td><td>4</td><td>4</td><td>4</td><td>-</td></tr><tr><td>b</td><td>0</td><td>5.4</td><td>2.2</td><td>-</td></tr><tr><td>T<sub>m0</sub></td><td>Eqn (6.30)</td><td>1800</td><td>1800</td><td>K</td></tr><tr><td>T<sub>m1</sub></td><td>Eqn (6.30)</td><td>36</td><td>36</td><td>K/GPa</td></tr><tr><td>α<sub>l</sub></td><td>Eqn (6.30)</td><td>5.6(l) 12(s)</td><td>1.5</td><td>-</td></tr></tbody></table>

<table><caption>Table 6.2: Alloy-independent model parameters.</caption><thead><tr><th>Variable</th><th>Value</th><th>Units</th><th>Reference</th></tr></thead><tbody><tr><td>R<sub>c</sub></td><td>390</td><td>km</td><td>-</td></tr><tr><td>c<sub>p</sub></td><td>835</td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>Nimmo (2009)</td></tr><tr><td>k</td><td>50</td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>de Koker et al. (2012)</td></tr><tr><td>L<sub>H</sub></td><td>300</td><td>kJ kg<sup>-1</sup></td><td>Stevenson et al. (1983)</td></tr><tr><td>α</td><td>9.2 10<sup>-5</sup></td><td>K<sup>-1</sup></td><td>Williams (2009)</td></tr><tr><td>f<sub>eff</sub></td><td>0.1</td><td>-</td><td></td></tr></tbody></table>

<figure><img src="image_43.png" /><figcaption>Figure 6.5: Inner-core boundary temperature as a function of inner core radius, assuming thermal equilibrium. Pictured is the Fe-FeS system with (left) $\chi_0 = 4$ wt.% and (right) $\chi_0 = 6$ wt.%.</figcaption></figure>

in order to compare the different behaviors. The steeper liquidus decrease as a function of inner core radius is observed in the case of sulfur, because the inner core is pure iron and thus the effect of light element enrichment in the outer core is the highest. On the other extreme, in the Fe-Si system, the core crystallises with the same composition as the outer core, and the liquidus slope is then due only to the pressure dependence. The Fe-C system is an intermediate case. It is important to note again that, for large inner core sizes, the outer core becomes very rapidly enriched in light element and could then transition to a non-standard state including top-down crystallisation or crystallisation at the eutectic.
<figure><img src="image_44.png" /><figcaption>Figure 6.6: Comparison of the effect of the alloy type on (left) inner-core boundary liquidus temperature as a function of inner core radius and (right) outer core enrichment as a function of inner core radius. The Si, S and C initial content is 3.64, 3.46 and 0.976 wt. %, respectively.</figcaption></figure>
# Lunar core evolution

## Contents

<table><tr><td>7.1</td><td>Foreword on 2D/3D geometry</td><td>82</td></tr><tr><td>7.2</td><td>Core-mantle evolution coupling</td><td>82</td></tr><tr><td>7.2.1</td><td>Influence on the core</td><td>82</td></tr><tr><td>7.2.2</td><td>Influence on the mantle</td><td>83</td></tr><tr><td>7.3</td><td>Core evolution : Typical case</td><td>84</td></tr><tr><td>7.4</td><td>Parameter space exploration</td><td>85</td></tr><tr><td>7.4.1</td><td>Influence of core composition</td><td>85</td></tr><tr><td>7.4.2</td><td>Influence of initial core thermal state</td><td>86</td></tr><tr><td>7.4.3</td><td>Other influences</td><td>87</td></tr><tr><td>7.5</td><td>Discussion</td><td>92</td></tr><tr><td>7.6</td><td>Conclusion</td><td>93</td></tr></table>

This chapter will present our main findings on lunar core evolution. The results were obtained both in a 2D cylindrical and 3D spherical geometry due to numerical reasons and we will discuss that influence in the first section. The second section will highlight the importance of core-mantle coupling both on core and mantle evolution. Then a typical case will be presented in order to understand the different processes at work. Finally, a systematic study of the influence of several parameters will be made. A summary of the 3D runs that were carried out and their main characteristics can be found in Table 7.1.

<table><caption>Table 7.1: Summary of the 3D simulations carried out for this section. The list of the 2D simulations is omitted as the qualitative results are irrelevant.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">ΔT</th><th rowspan="2">χ<sub>0</sub></th><th colspan="3">Final r<sub>i</sub></th></tr><tr><th>r<sub>c</sub></th><th>μT</th><th>Chem. dynamo onset</th></tr><tr><th></th><th>K</th><th>wt.%</th><th></th><th></th><th></th></tr></thead><tbody><tr><td>dT200_X8</td><td>200</td><td>8</td><td>0.31</td><td>0.07</td><td>-</td></tr><tr><td>dT300_X4</td><td>300</td><td>4</td><td>0.51</td><td>0.16</td><td>4.3</td></tr><tr><td>dT300_X6</td><td>300</td><td>6</td><td>0.42</td><td>0.12</td><td>4.0</td></tr><tr><td>dT400_X2</td><td>400</td><td>2</td><td>0.61</td><td>0.20</td><td>4.3</td></tr><tr><td>dT400_X4</td><td>400</td><td>4</td><td>0.50</td><td>0.17</td><td>3.9</td></tr></tbody></table>
## 7.1 Foreword on 2D/3D geometry

Thermochemical convection simulations in 3D are very demanding in computational resources. It was therefore impossible to conduct all tests in 3D and some were thus prepared in a 2D cylindrical geometry. To keep core cooling consistent between the two geometries, it is necessary to scale the inner core radius in 2D so that the ratio of the upper to lower boundary (i.e., surface and core-mantle boundary) stays the same (Hernlund and Tackley, 2008). This can be written as $r_{oc}/r_{ic} = (r_{os}/r_{is})^2$ where $r_{oc}$ and $r_{ic}$ refer to the outer and inner radius of the cylinder and $r_{os}$ and $r_{is}$ to the outer and inner radius of the sphere. Keeping the thickness of the mantle constant ($r_{oc} - r_{ic} = r_{os} - r_{is}$) assures uniqueness of the solution.

For the Moon, $r_{os} = 1740$ km and $r_{is} = 330$ km. Thus in 2D the core has to have a radius of $r_{ic} = 52$ km. This is particularly small and causes both numerical and theoretical issues. Indeed even if the core-mantle boundary heat flow is equivalent to the 3D case with such a radius, the implied curvature completely changes the nature of the flow. We therefore avoided scaling altogether when using 2D, knowing that the qualitative results will be incorrect, but still allowing us to test the effect of several parameters in a practical amount of time.

## 7.2 Core-mantle evolution coupling

Most models of planetary mantle thermochemical evolution use the core as a heat reservoir that can be cooled with time without considering inner core growth, whatever the temperature is. However, in systems where crystallisation of the inner core should occur, this simplification has an important effect both on the core and mantle.

### 7.2.1 Influence on the core

After a period of secular cooling, where the core is cooled to its liquidus temperature and during which a thermal dynamo can be sustained, crystallisation occurs. This is an exothermic reaction and therefore slows down core cooling. For our benchmark case in Figure (7.1), this happened after about 100 Ma.

After the inner core starts to grow (see Fig. 7.1b), a large part of the heat extracted from the core is the latent heat of crystallisation. This can be seen in Figure (7.1d) where almost half of the extracted heat is latent heat. That means that the core is cooling much less efficiently and its inner part is growing very slowly. During inner core growth, we assume thermal equilibrium, i.e., the core temperature at the inner core boundary has to be the liquidus temperature at that point. From Figure (6.5), we can see that for a core with 4 wt.% S (left panel), the drop in liquidus temperature after about 35% crystallisation is less than 10 K, and thus the core is not expected to cool more than that amount in the process.

Figure (7.1a) shows the core temperature evolution in the case where core crystallisation is and isn't considered (i.e., with and without feedback). When feedback
<figure><img src="image_45.png" /></figure>

Figure 7.1: (a) Core temperature, (b) inner core size, (c) core-mantle boundary heat flow and (d) contributions to the core-mantle boundary heat flow for a case with and without feedback between between inner core growth and mantle evolution. These results were obtained by non-scaled 2D simulations and should therefore only be read qualitatively.

is considered, both the latent heat and compositional buoyancy term of Equation (6.1) are used. It is clear that the case with core crystallisation cools less than the case without feedback. The difference of almost 100 K in this figure has a direct influence on the core mantle boundary heat flow, because in the meantime the mantle is still cooling. Figure (7.1c) shows that the heat flow starts to increase about 4 Ga ago with respect to the simulation which doesn't include feedback from core crystallisation, due to the larger temperature difference between core and mantle.

### 7.2.2 Influence on the mantle

When positive, the core-mantle boundary heat flow was a energy sink from the core point of view. From the mantle point of view, it is now a heat source. Therefore as the temperature difference between core and mantle increases, the importance of this heat source increases. This can be seen in Figure (7.1c). But even if there is a factor of about 2 difference between the cases with and without feedback from core crystallisation, the effect on the mantle remains small. Figure (7.2) shows a
comparison of the average temperature profile at the present day. A difference of 50 K is observed in the lower mantle. The effect on melt production is shown in Figure (7.2) and is seen to be relatively small. However the lower mantle thermal state is different and is thus dependent on the core evolution. In the case of the Moon, this could explain the presence of partially molten material at the core mantle boundary.

<figure><img src="image_46.png" /><figcaption>Figure 7.2: These results were obtained by non-scaled 2D simulations and should therefore only be read qualitatively. (a) Lower mantle temperature profile today and (b) melt production as a function of time.</figcaption></figure>

## 7.3 Core evolution : Typical case

The previous section showed qualitative results in the case of a 2D simulation. We now present our 3D benchmark case in Figure (7.3), corresponding to a simulation with 6 wt.% S initially, and a core-mantle temperature difference of 300 K.

As described earlier, the onset of inner core crystallisation almost stops core cooling. This occurs 4 Ga ago in this model. Figure (7.3d) presents the predicted magnetic field. A thermal dynamo is present during the first 100 million years and then no magnetic field is predicted for about 500 million years. This corresponds to a period during which the heat flow out of the core is not large enough to sustain a thermal dynamo and the growth of an inner core has not started yet to sustain the dynamo chemically.

The chemical dynamo is present from about 3.8 to 2.4 Ga ago, marginally consistent with paleomagnetic observations (which estimate the presence of a field between 4.1 and 3.1 Ga ago). This is the model that best fit the timing constraints we could get from 3D simulations as we could not do a full parameter study due to the heavy computational requirements. However a large discrepancy still exist with regards to the field's magnitude (about 0.1 µT versus about 10 µT expected). We will show in the next section the influence of several parameters on these results using 2D simulations.

Occulting the issue of the field strength, which is common to all lunar paleo-
magnetic models, a major problem of this benchmark run is that a new period of
chemical dynamo era is predicted to occur between 2 Ga ago and the present day.
This happens because the core is no longer cooling efficiently as its temperature
is dictated by the outer core liquidus and thus, as the mantle is cooling, the core-
mantle boundary heat flow increases. This in turn leads to a faster inner core growth
rate. We will discuss later how that problem could be avoided.

<figure><img src="image_47.png" /><figcaption>Figure 7.3: Results obtained in 3D for a case with 6 wt.% S and $\Delta T = 300$ K; model dT300_X6 in Table 7.1 (a) Core-mantle boundary heat flow, (b) core temperature, (c) inner core size and (d) estimated surface magnetic field as a function of time. The expected core adiabatic heat flow is about 5 mW/m$^2$.</figcaption></figure>

7.4 Parameter space exploration

7.4.1 Influence of core composition

The influence of the core concentration in light elements is twofold. It has an effect
on the liquidus temperature and thus on the timing of crystallisation and final size
of the inner core and an effect on the density drop between solid inner core and
liquid outer core, which will influence the dissipation rate.

The most important aspect is the effect on liquidus temperature. Part of that
effect is due to the initial concentration and part is due to the crystallisation behavior
of the particular alloy. The higher the initial light element concentration, the lower
the liquidus temperature will be. The onset of inner core crystallisation will thus be
delayed in such cases. This can be seen in Figure 7.4. Another direct consequence is
that the core will be on average warmer when less light element is present and thus
the core-mantle boundary heat flow will be larger (for a given mantle temperature).
This will lead to slightly larger growth rates for the inner core and consequently,
to larger dissipation rates. This qualitative behavior is independent of the alloying
system we consider.

However specificities between the different systems arise when considering the
enrichment of the outer core in light elements during inner core crystallisation. A
rapid increase in light element in the outer core (such as what happens in the Fe-FeS
system) leads to a rapid decrease of the outer core liquidus, and therefore allows for
a faster core cooling whereas in a case where the outer core is not enriched in light
elements (in the Fe-FeSi system, for example) the outer core liquidus only decreases
due to the pressure effect and the core cools only slightly. The Fe-FeS system is
therefore the most promising case initially (as it is the one where the chemical
differentiation is the strongest). As a consequence, we focused on that system and
tried to understand the different influences on the results better, keeping in mind
that another alloying element could slightly change the conclusions.

The second influence of core concentration in light elements is on the density
drop between inner and outer core during inner core growth. This is directly linked
to the magnitude of the buoyancy forces at the inner core boundary and thus,
on the dissipation rate. Again we expect the Fe-FeS system to have the largest
influence here as the light element is entirely segregated in the outer core. The higher
dissipation rate leads to a longer dynamo activity and a higher estimated surface
magnetic field. However the difference, on the order of 0.1 µT, is not large enough to
explain the 10 µT order-of-magnitude difference with paleomagnetic measurements.

7.4.2 Influence of initial core thermal state

The initial temperature difference between core and mantle influences mostly the
early lunar evolution. Figure (7.5) shows that a smaller temperature difference
results in an earlier growth of the inner core but that the present day value remains
unchanged. Most of the magnetic field predictions are also unchanged. The most
important consequence concerns the earliest stages of the evolution. Figure (7.5d)
shows that the case with the lowest temperature difference does not present a gap
in the magnetic field history whereas the case with 400 K has a 400 million years
gap between the end of the thermal dynamo and the onset of the chemical one.
This is due to the fact that in the small temperature difference case, the chemical
dynamo starts earlier. This effect can become interesting when more data will be
available from that time period to discriminate between cases with paleomagnetism
as ancient as 4.5 Ga and cases where the oldest recorded are 4.2 Ga.
<figure><img src="image_48.png" /><figcaption>Figure 7.4: Comparison between a case with 4 and 2 wt.% S in 3D. (a) Core-mantle boundary heat flow, (b) core temperature, (c) inner core size and (d) estimated surface magnetic field as a function of time. These are runs dT400_X2 and dT400_X4 from Table 7.1.</figcaption></figure>

### 7.4.3 Other influences

All the influences presented here have been tested using a 2D cylindrical version of the code for numerical reasons.

#### Core thermal expansivity

The thermal expansivity of the core has a major influence on the adiabatic gradient through the scale height factor D. The larger thermal expansivity is, the smaller D is, which in turns leads to a steeper gradient. This translates into lower values of $\Delta$, the ratio of liquidus to adiabatic gradient.

These parameters are non-trivially present in Equations (6.2)-(6.4) and the equivalent entropy terms. However from Figure (7.6) it appears that a low value of the thermal expansivity leads to a higher growth rate and thus to higher dissipation rates. In particular Figure (7.6) shows that everything else being equal, a case with a lower expansivity produced enough dissipation to sustain a magnetic field for about 1 billion years, where a higher value could not.

For this comparison we used the low value in Williams (2009) as the nominal
<figure><img src="image_49.png" /></figure>

Figure 7.5: Comparison between a case with initial core mantle temperature difference of 400 and 200 K. (a) Core-mantle boundary heat flow, (b) core temperature, (c) inner core size and (d) estimated surface magnetic field as a function of time. These are runs dT400_X4 and dT300_X4 from Table 7.1.

case ($9.2 \ 10^{-5} \ \text{K}^{-1}$), and Figure 7.6 was made with $\pm$ 10% variations around the nominal case (ie. $8.3$ and $10.1 \ 10^{-5} \ \text{K}^{-1}$). The range of reasonable values for thermal expansivity of pure Fe at lunar pressure is $8.2$ to $13.2 \ 10^{-5} \ \text{K}^{-1}$ in extreme cases (Williams, 2009). We therefore used one of the lowest possible value for thermal expansivity and the $\pm$ 10% range considered in this test appears to be non-negligible.

## Core latent heat of crystallisation

The latent heat of crystallisation directly influences the amount of heat produced while crystallising the inner core. As a consequence, the higher the value, the warmer the core will stay and thus, the smaller the inner core will be. A direct consequence is that less dissipation occurs in such a case. This can be seen in Figure (7.7). However the effect is very small considering the other uncertainties.

We used here a high and low value of $750$ and $250 \text{ kJ kg}^{-1}$, respectively. As discussed in *Stevenson et al.* (1983), the latent heat of crystallisation of pure Fe is $250 \text{ kJ kg}^{-1}$ at ambient pressure and its value increases with the melting temperature. The melting point increases by about 10% by 4 GPa thus a value close to 300 kJ
<figure><img src="image_50.png" /><figcaption>Figure 7.6: Influence of thermal expansivity on core processes. (a) Inner core size and (b) dissipated power in the core. These results were obtained by non-scaled 2D simulations and should therefore only be read qualitatively.</figcaption></figure>

$kg^{-1}$ may be considered nominal on the Moon. The high value was chosen from Nimmo (2009). Within the range of accepted values, variations in latent heat can therefore be neglected.

<figure><img src="image_51.png" /><figcaption>Figure 7.7: Influence of latent heat of crystallisation on core processes. (a) Inner core size and (b) dissipated power in the core. These results were obtained by non-scaled 2D simulations and should therefore only be read qualitatively.</figcaption></figure>

Core heat capacity

The heat capacity has a subtler influence. It as an effect on the adiabatic gradient, like the thermal expansivity and a direct effect on the efficiency of secular cooling as well. A higher value of the heat capacity leads to a smaller gradient and to less efficient secular cooling. Figure (7.8) shows that the effect on the gradient seems to be the most important as the case with a larger heat capacity produces a larger core. It has a strong effect on dissipation as well, similar to that of thermal expansivity.
The larger the heat capacity is, the stronger dissipation will be, because the entropy sink linked to conduction will be smaller.

In that comparison we used a nominal value of 835 J kg⁻¹ K⁻¹ *Williams* (2009) and arbitrarily considered a low and high case with values ± 10% different (ie. 750 and 920 J kg⁻¹ K⁻¹, respectively). However, even if the effect of this 10% is large, the value of heat capacity is fairly well known for a few decades already, and is not expected to vary much around 835 J kg⁻¹ K⁻¹ (*Desai*, 1986).

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

Figure 7.8: Influence of heat capacity on core processes. (a) Inner core size and (b) dissipated power in the core. These results were obtained by non-scaled 2D simulations and should therefore only be read qualitatively.

## Core thermal conductivity

The core thermal conductivity does not play any role in the energy budget. However, a large value increases the amount of heat that can be carried on the adiabat, which means that a larger heat flow is necessary to sustain the dynamo. Figure (7.9) shows that the effect of conductivity on dissipation is indeed straightforward, with the case having a low conductivity constantly higher than the case with a high conductivity.

Our nominal case is 50 W m⁻¹ K⁻¹ (corresponding to the low case from this Figure) and the high bound case has a value of 90 W m⁻¹ K⁻¹. We based our assumptions on a recent work by *de Koker et al.* (2012) which showed that at 2000 K, conductivity of pure Fe should be about 60 W m⁻¹ K⁻¹, with that value increasing with temperature and pressure. Even if the observed effect is fairly strong, such a large value is not expected so the influence of our result on that parameter will be small.

## Mantle heat sources distribution

In this section, we tested the influence of the heat sources distribution on the core evolution. Figure (7.10) shows that the case with heat sources localised in one region of the upper mantle triggers crystallisation earlier due to a faster onset of
<figure><img src="image_53.png" /></figure>

Figure 7.9: Influence of thermal conductivity on core processes. (a) Inner core size and (b) dissipated power in the core. These results were obtained by non-scaled 2D simulations and should therefore only be read qualitatively.

convection. However there is no major influence on the long run, therefore as long as these heat sources are localised in the upper mantle, their exact distribution plays a minor role. However scenarios where heat sources are localised in the lower mantle do exist.

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

Figure 7.10: Influence of bulk heat sources distribution on core processes. (a) Inner core size and (b) dissipated power in the core. These results were obtained by non-scaled 2D simulations and should therefore only be read qualitatively.

Such a scenario was presented by Hess and Parmentier (1995). We did not fully model such a case but we note that this was studied recently by Zhang et al. (2013) who claimed that such a model can explain the duration of thermal dynamo activity. We used their inner core growth rate and core-mantle boundary heat flow in our model to predict the influence of a chemical dynamo. While their thermal dynamo is indeed consistent with observations, their inner core growth rate suggest the existence of a chemical contribution as well up to the present day. These predictions can be found in Figure (7.11).
<figure><img src="image_55.png" /><figcaption>Figure 7.11: (left) Inner core growth rate from Zhang et al. (2013) and (right) surface magnetic field estimate from our model, using the inner core growth rate and core-mantle boundary heat flow from their study as input.</figcaption></figure>

## 7.5 Discussion

The first important point is that this study assumes bottom-up crystallization scenarios only. Williams (2009) showed that the Moon lies at the boundary between the two domains and an initial sulfur content higher than about 8 to 10 wt.% could lead to the iron snow regime (Hauck et al., 2006). In addition, a transition between the two regimes can occur as the outer core is enriched in light element during evolution. The generation of a chemical dynamo from iron-snow has never been studied in detail. Some authors have argued that remixing is efficient enough to produce the magnetic field on Ganymede for example Hauck et al. (2006), but in any case the magnetic scaling of such a flow is unknown. If it turns out that remixing is not efficient, or if the convecting zone is too thin, this could help shut down the long lasting dynamo observed in Figures (7.3)-(7.5). We therefore consider this possibility out of the scope of this study, but the possibility of such exotic scenarios should be kept in mind.

Concerning our results, as seen in Figures (7.3)-(7.5), it is not straightforward to obtain a 3D case which completely fits paleomagnetic data. A case with 6 wt.% S and initial core mantle temperature difference of 300 K is our best-fitting case, but it still predicts a period of dynamo action up to the present day. From our simulations we believe that, for the particular mantle assumptions that we used, a late stage increase in core mantle boundary heat flow and thus a renewal of dynamo action will always be present if one also wants to predict a global magnetic field around 4 Ga ago (but cases with no chemical dynamo and no inner core at all are also possible). We suggest that this could be explained if the outer core transitioned to an iron snow regime at some point, or if the convective layer became too small for convection to occur. A better fit can probably also be obtained by varying several of the parameters presented in the previous section. In any case, this should be the target for future studies.
For consistency, we used a similar mantle setup as *Laneuville et al.* (2013)'s preferred model. As can be seen in Figure (7.10), the exact heat source distribution within the crust plays little role in the core evolution and we are therefore confident that this is not an important factor for the lunar dynamo. However, other effects should be studied in the future: (1) an increased amount of heat sources in the mantle, (2) a wetter rheology (ie. lower reference viscosity to account for a non-negligible mantle water content) and (3) different initial mantle states.

More heat sources in the mantle would slow down mantle cooling and thus reduce the possibility of a long lasting dynamo. The effect of a lower viscosity should be to cool the mantle more efficiently, and thus to allow smaller initial core mantle temperature differences to still trigger dynamo action. Finally different initial mantle states (both temperature and density profiles) are important to consider as well, as their influence is harder to predict. For example *Zhang et al.* (2013) studied a case with a dense layer enriched in heat sources on top of the core-mantle boundary with important conclusions about paleomagnetism.

## 7.6 Conclusion

It is hard to reconcile a chemically induced dynamo with the lack of active magnetic field today. We have shown that two families of model exist: the first class correctly predicts a magnetic field era centered at about 3.5 Ga ago, but also predicts a magnetic field up to the present day. The second class avoids the present day magnetic field problem by not growing an inner core at all. In that case another scenario has to be used to explain the long lasting paleomagnetic era, and seismic evidence for the existence of an inner core also poses a problem.

If we trust the prediction about the inner core, then the first class of models has to be considered. They present two discrepancies with the observations. The first is that the predicted magnetic field lasts too long. Several arguments exist which could explain the difference, the most promising being the fact that outer core enrichment in light elements could produce a change in its dynamical behavior and shut down the dynamo earlier than currently predicted. This should be the target of future studies. Another possible explanation is a failure of the magnetic field scaling law in a thin shell.

The second discrepancy is that predicted field strengths are an order of magnitude too small when compared to paleomagnetic studies. One possible explanation is that several processes occur at the same time. For example if precession induced flow exists at the same time as the chemical dynamo, a larger magnetic field could be produced. The fluid dynamics of such a system has to be studied to assess the additivity of these processes. Another way to reconcile predictions with the measured paleointensity of the magnetic field is to use the magnetic field focussing effect proposed by *Hide* (1972). However this is only speculative and more modeling effort should be done in that area.
<figure><img src="image_56.png" /></figure>
# Part III

## Summary and perspectives
<figure><img src="image_57.png" /></figure>