# Accrétion et différenciation de la Terre primitive :
## Évolution thermochimique jusqu'à aujourd'hui.
## Mémoire présenté en vue de l'obtention de l'Habilitation à Diriger des Recherches

Julien Monteux

► To cite this version:

Julien Monteux. Accrétion et différenciation de la Terre primitive : Évolution thermochimique jusqu'à
aujourd'hui. Mémoire présenté en vue de l'obtention de l'Habilitation à Diriger des Recherches.
Géophysique [physics.geo-ph]. Université Clermont Auvergne, 2024. tel-04420188

HAL Id: tel-04420188

https://hal.science/tel-04420188

Submitted on 26 Jan 2024

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
Julien Monteux

Accrétion et différenciation de la Terre primitive :
Évolution thermochimique jusqu'à aujourd'hui
<table><tr><td>10 Annexes</td><td>565</td></tr><tr><td>10.1 Annexe 1 : Curriculum Vitae</td><td>566</td></tr><tr><td>10.2 Annexe 2 : Valorisation des travaux</td><td>571</td></tr><tr><td>10.3 Annexe 3 : Publications dans d'autres thématiques</td><td>582</td></tr></table>
Remerciements

Je souhaite remercier ici...

Les membres du jury qui ont accepté de juger les travaux résumés dans ce manuscrit, Anne Davaille, Chloé Michaut, Philippe Cardin, Denis Andrault, Maud Boyet et Gael Choblet.

Les personnes qui m'ont aidées dans le parcours scientifique qui m'a amené aujourd'hui à Clermont-Ferrand en particulier Riad Hassani, Jean-Luc Got, Nicolas Coltice, Yanick Ricard, Mark Jellinek, Catherine Johnson, Gabriel Tobie, Gaël Choblet, Laurent Husson.

Le Laboratoire Magmas et Volcans pour son accueil chaleureux.

Ma famille et mes amis.

Pauline, Martin, Valentine, Manon.
Contexte

La géodynamique occupe une place centrale dans la compréhension de la structure et de
l'évolution thermique interne des planètes telluriques du système solaire puisqu'elle permet
d'interpréter et de relier les différents types d'observations (sismique, géochimique, minéra-
logiques ...). L'activité de recherche que je mène est dédiée à la modélisation numérique de
phénomènes géologiques et planétologiques.

Mon parcours scientifique m'a permis de réaliser mon stage de recherche de M2R dans l'équipe "Géophysique des Volcans" du Laboratoire de Géophysique Interne et de Tectono-physique à Chambéry (encadrants : J.L. Got et R. Hassani), puis ma thèse au Laboratoire de Sciences de la Terre de Lyon (encadrants : Y. Ricard et N. Coltice), suivi d'un stage post-doctoral au département Earth and Ocean Sciences de Vancouver (collaborateurs : M. Jellinek et C. Johnson). Ce parcours m'a mené ensuite au Laboratoire de Planétologie et de Géody-namique de Nantes où j'ai obtenu un financement A.N.R. Retour Post-Doctorant afin de dé-velopper le projet nommé ACCRETIS afin de contraindre la différenciation des satellites de glace (collaborateurs : G. Tobie, G. Choblet, B. Langlais, H. Amit). Aujourd'hui, je continue de développer des modèles numériques afin de caractériser l'évolution thermo-chimique de la Terre primitive au Laboratoire Magmas et Volcans de Clermont-Ferrand (collaborateurs : A. Bouhifd, M. Boyet, D. Andrault, G. Manthilake, S. Demouchy, M. Guitreau).
# 1 Introduction

La Terre est une formidable machine thermique qui tend à effacer les traces des évènements les plus anciens. Cependant, les données géochimiques, la présence de la Lune et certains observables géophysiques gardent encore aujourd'hui la mémoire d'évènements cruciaux pour l'évolution de la Terre qui se sont déroulés il y a plus de 4 milliards d'années. Notre connaissance de l'histoire de l'intérieur de la Terre a bénéficié de l'avènement de l'exploration planétaire, des outils de calcul numérique et d'analyse isotopique mais aussi des avancées importantes en pétrologie expérimentale.

Le bombardement météoritique à l'origine de la formation de la Terre et des autres objets du système solaire a eu lieu il y a environ 4.5 Ga. Ce phénomène d'accrétion s'est déroulé en peu de temps, quelques millions d'années après la formation des premiers solides (i.e. Calcium-Aluminium-rich Inclusions, C.A.I.) (Kleine and Rudge, 2011). Les impacts météoritiques ont ainsi joué un rôle crucial dans les étapes primitives d'évolution de la Terre gouvernant notamment son budget énergétique et chimique ainsi que les modalités de sa différenciation interne et de l'initiation de son champ magnétique. L'accumulation de matériel a formé des corps de plus en plus volumineux et ainsi de plus en plus rares qui ont eux même impacté des corps plus gros. Au bout de quelques millions d'années, les collisions ont impliqué des objets de plusieurs centaines de kilomètres de rayon (Agnor et al., 1999; Kokubo and Ida, 2000).

Au cours de l'accrétion de la Terre, les contributions des grands impacts entre planétismaux et la désintégration des nucléides radioactifs à courte durée de vie ont induit des tempéritures internes très élevées. Moins de 100 millions d'années plus tard, l'impact majeur qui a formé la Lune a fourni une énorme quantité d'énergie, faisant peut-être fondre entièrement la Terre. Ensuite, le refroidissement rapide de l'océan magmatique global a entraîné la cristallisation progressive du manteau et la formation d'une importante couche limite thermique à la base du manteau inférieur. Cette couche inférieure (souvent appelée couche "D") a toujours un impact important sur le bilan thermique de la Terre, car elle est étroitement liée au flux de chaleur à la limite noyau-manteau (CMB). Aujourd'hui, 4,5 milliards d'années (Ga) après sa formation, la Terre est toujours habitable, non seulement en raison des conditions clémentes qui règnent à sa surface, mais aussi grâce à son activité interne. La persistance du champ géomagnétique, en particulier, est essentielle pour protéger la vie des radiations solaires.

Il existe un consensus sur l'importance des premières étapes de formation des planètes et lunes du Système Solaire dans l'acquisition de leurs caractéristiques actuelles en termes de composition chimique, de dynamique interne et de formation de leur atmosphère (Kleine et al., 2002; Halliday, 2004). Cependant, les modalités d'évolutions primitives et les processus mis en jeu restent largement débattus. Pour répondre aux grandes questions et défis majeurs des Sciences de la Terre et de l'Univers, la planétologie comparée et les simulations numériques constituent des outils performants pour comprendre et reconstituer les phénomènes qui ont pu se dérouler il y a plus de 4.5 milliards d'années. Les planètes telluriques se sont formées
rapidement, dans les 100 premiers millions d'années après la formation des premiers solides du Système Solaire (Kleine et al., 2009), et les impacts géants étaient probablement fréquents durant les dernières étapes de l'accrétion (Canup and Agnor, 2000; Jacobson and Morbidelli, 2014). Le développement de méthodes d'analyses géochimiques a permis de contraindre la chronologie de la différenciation des planètes telluriques. Les grands événements de différenciation (formation du noyau, cristallisation de l'océan magmatique) se sont produits dans les premiers 100 Ma (Kleine et al., 2002; Boyet and Carlson, 2005).

L'intérieur de la Terre a ainsi subi une succession d'évènements majeurs depuis la formation de la Terre jusqu'à ce que nous en connaissons aujourd'hui notamment grâce à l'apport de la tomographie sismique. Ces évènements tels que les impacts géants, la séparation manteau/noyau ou la solidification du manteau primitif ont eu des conséquences non seulement thermiques avec des redistributions importantes de la chaleur primitive, dynamiques avec des transitions majeures au niveau des régimes d'écoulements mais aussi chimiques avec la mise en place de ségrégation plus ou moins complètes au sein du réservoir mantellique. Ces évènements ont aussi pu conditionner les modalités de démarrage du champ magnétique terrestre. Si le manteau constitue environ 85% du volume de la Terre, il n'en demeure pas moins mystérieux et son évolution primitive reste encore aujourd'hui largement débattue dans la communauté scientifique. Aujourd'hui l'essentiel du manteau est solide et son évolution thermo-chimique semble plutôt lente. Cependant il n'en a pas toujours été ainsi.

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Chronologie des grandes étapes de l'évolution interne de la Terre.</figcaption></figure>

Au vu de la grande diversité des échelles caractéristiques de temps, de températures, de vitesses et de tailles, les simulations numériques sont des outils particulièrement adaptés pour
atteindre cet objectif et ainsi mieux comprendre :

— Par quels processus les planètes ont-elles acquis leur structure interne ?

— Quel est le rôle des impacts dans l'évolution thermique des intérieurs planétaires ?

— Quelle a été l'influence des paramètres d'accrétion sur la composition des grands réservoirs à travers la formation du noyau et la cristallisation d'océans magmatiques ?

Pendant mon parcours scientifique, j'ai étudié le rôle des impacts météoritiques et des grandes phases d'accrétion sur l'évolution de différents objets du Système Solaire : planètes telluriques, lunes, satellites de glace. Pour cela j'ai développé des outils de modélisation (numériques et analogiques) dans diverses géométries (1D, 2D, 3D, cartésien, cylindrique, sphérique). Ces modèles m'ont permis de mettre en évidence le rôle central des étapes primitives de formation et notamment des impacts géants dans l'acquisition des caractéristiques principales des objets du Système Solaire. Dans ce document, je propose de retracer l'évolution interne de la Terre depuis les dernières étapes de son accrétion jusqu'à la mise en place d'un champ magnétique, d'une convection mantellique et d'une tectonique des plaques (Figure 1). Nous verrons notamment le rôle central qu'ont joué les impacts météoritiques dans l'acquisition de ces caractéristiques et notamment sur le budget thermique primitif de la Terre, des planètes telluriques plus généralement mais aussi des satellites de glace.
## 2 Le rôle des impacts dans l'évolution thermique précoce des planètes

Les impacts météoritiques ont joué un rôle crucial dans les étapes primitives d'évolution de la Terre gouvernant notamment sa différenciation interne jusqu'à la formation de son atmosphère (Abe et al., 2000; Monteux et al., 2009; Ricard et al., 2009; Lebrun et al., 2013). Au cours de la dernière étape d'accrétion, l'essentiel de la masse des objets telluriques s'accumule et les protoplanètes deviennent des planètes. Cette étape a fait intervenir des phénomènes d'impact extrêmement énergétiques entre des embryons planétaires potentiellement différenciés. Ces impacts ont provoqué des phénomènes de fusion des roches terrestres à grandes échelle et la formation d'océans magmatiques au sein desquels le fer et les silicates se sont séparés pour former le noyau et le manteau (Tonks and Melosh, 1993). Malgré l'importance qu'ont eu les impacts météoritiques pour façonner la structure thermo-chimique interne de la Terre telle que nous la connaissons aujourd'hui, ces phénomènes sont encore peu contraints tant qualitativement que quantitativement et de nombreuses questions crucales sont encore en suspens. Par exemple il est actuellement difficile de savoir comment est dissipée l'énergie cinétique de l'impacteur et quel est l'état thermique de la protoplanète après un impact. Les impacts météoritiques étant des phénomènes hautement complexes et énergétiques, mieux les comprendre pourrait permettre de mieux anticiper les différentes étapes de fusion et de recristallisation de la planète au cours de son accrétion et l'état thermique final auquel abouti la planète après des impacts répétés.

Un impact météoritique s'il se produit à des vitesses supérieures aux vitesses des ondes élastiques crée une onde de choc dans le matériel impacté ce qui provoque un changement de l'état thermodynamique du matériel impacteur et impacté. Le passage de l'onde de choc provoque une compression brutale qui se traduit par une addition irréversible d'énergie provoquant ainsi une hausse de la pression et de l'entropie (Figure 2). Après le passage de l'onde de choc, le matériel se décomprime jusqu'à pression ambiante de façon adiabatique. En connaissant les équations d'état et les caractéristiques thermodynamiques des matériaux mis en jeu (capacité calorifique, chaleur latente de fusion ...), on peut connaitre l'état thermodynamique (pression, température, taux de fusion ...) du matériel post-impacté.

### 2.1 Modélisation des impacts

Les impacts géants peuvent entraîner une augmentation significative de la pression et de la température sous le point d'impact et ainsi des déformations majeures au sein du corps impacté. En plus de la taille et de la vitesse des impacteurs, d'autres paramètres peuvent jouer un rôle important sur comment l'énergie cinétique de l'impacteur est transmise au corps impacté : la rhéologie, la porosité, la composition chimique des corps impliqués ou encore l'angle d'impact. Pour caractériser les effets de ces différents paramètres, il est désormais possible d'utiliser des modèles numériques pour simuler des impacts entre deux corps. Les hydrocodes (ou codes
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – La formation et la propagation des ondes de choc et de raréfaction par analogie avec la collision de deux trains. P, T, E et V sont respectivement la pression, la température, l'énergie interne et le volume dans la matière non comprimée (état 0) et comprimée (état 1). $u_p$ est la vitesse des particules derrière l'onde de choc et $U_s$ la vitesse du choc (D'après Langenhorst and Deutsch (2012)).</figcaption></figure>

d'impacts) permettent de simuler numériquement des impacts météoritiques et le réajustement thermo-mécanique post-impact (ex. Jutzi and Asphaug (2011); Ćuk and Stewart (2012)).

Les simulations de processus dynamiques nécessitent la résolution des équations de conservation (masse, quantité de mouvement et énergie) pour un matériau donné. En plus de résoudre ces équations, les codes qui permettent de simuler les processus d'impact sont couplés à un modèle constitutif qui combine deux parties : la réponse volumétrique du matériau à l'onde de choc représentée par l'équation d'état du matériau (EOS) (eg. Tillotson, ANEOS...) , et la réponse aux déformations déviatoriques représentée par un modèle de résistance (Johnson-Cook, Ivanov, Collins...). La précision des prédictions du modèle dépend de la précision avec laquelle les modèles de matériaux (EOS + résistance) reproduisent le comportement du matériau, et de la connaissance des propriétés incluses dans ces modèles (Pierazzo et al., 2008; Stickle et al., 2020). Bien que les équations de conservation soient cohérentes entre les différents modèles numériques d'impact, la façon dont elles sont résolues peut varier. Il existe deux méthodes principales pour résoudre les équations du mouvement : une solution eulérienne et une solution lagrangienne. La solution eulérienne (ex. code CTH) traite les variables du domaine (ex. densité, température, pression) à partir d'un maillage fixe tandis que la méthode lagrangienne (ex. code Smooth Particle Hydrodynamic) traite les variables du domaine via un maillage qui se déplace dans l'espace et se déforme en même temps que le matériau. Une approche hybride
appelée Arbitrary Lagrangian Eulerian (ALE) permet d'utiliser un mélange de schémas lagrangiens et eulériens dans lesquels la solution d'un maillage lagrangien déformé est reportée sur un maillage eulérien fixé dans l'espace (Benson, 1992).

Ces codes d'impact (ou hydrocodes) se sont récemment avérés être des outils d'un grand intérêt pour la compréhension des phénomènes d'impacts géants et ont notamment permis de confirmer qu'un impact géant entre la Terre et un objet avec un rayon compris entre 3400 et 6370 km (nommé Théia) avait pu générer un champ de débris suffisant autour de la Terre pour former la Lune (Canup, 2004; Ćuk and Stewart, 2012). Ces codes permettent de modéliser l'impacteur et le corps impacté sous forme d'assemblage de plusieurs millions de blocs plurikilométriques et intègrent des équations d'état de plusieurs matériaux géologiques (ex. fer, basalte, dunite...) permettant ainsi de calculer précisément l'évolution Pression/Température de chaque matériau. Parmi eux, le code iSALE (Collins et al., 2008) est un hydrocode open-source (bénéficiant d'une version 2D et 3D) largement utilisé par la communauté scientifique de la planétologie. Il a été développé pour modéliser les écoulements de fluides à différentes vitesses. iSALE est ainsi particulièrement pertinent pour mesurer les conséquences internes de larges impacts météoritiques (Davison et al., 2010).

## 2.2 Augmentation de pression au cours d'un impact

Malgré la fréquence élevée des impacts à faibles vitesses (<10 km/s) pendant l'accrétion des corps telluriques du Système Solaire (Agnor et al., 1999), les caractéristiques de la propagation des ondes de choc lors de tels impacts étaient jusqu'à récemment peu contraintes (Pierazzo et al., 1997). Dans une étude développée en collaboration avec J. Arkani-Hamed, nous avons utilisé des simulations hydrocodes iSALE-2D afin de calculer l'évolution de la pression liée à l'impact et les vitesses de propagation de l'onde de choc sur un corps tellurique de la taille de Mars et pour des vitesses d'impact comprises entre 4 et 10 km/s (Figure 3) (Monteux and Arkani-Hamed, 2016). Nous avons considéré dans cette étude de gros impacts avec des diamètres compris entre 100 et 400 km comparables à ceux ayant généré des cratères d'impacts martiens comme Hellas ou Argyre. Nos résultats montrent que l'intensité de l'onde de choc se répartit en 3 zones sphériques dans le manteau autour du point d'impact et que cette intensité diminue avec l'éloignement au point d'impact. Nous avons caractérisé la répartition de la pression post-impact par des lois d'échelles ce qui représente une première étape primordiale pour déterminer l'intensité du chauffage lié à un impact sur une protoplanète en croissance.

Alors que les impacts majeurs lors de l'accrétion tardive d'une planète de type Mars se produisent sur un corps différencié, les caractéristiques de la propagation de l'onde de choc sont mal connues à l'intérieur de ces objets avec un manteau rocheux et un noyau liquide. Dans une deuxième étude (Monteux and Arkani-Hamed, 2019), nous avons utilisé d'autres simulations
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Évolution temporelle du champ de composition et de la pression lors du passage de l'onde de choc générée par un impact avec une vitesse $V_{imp} = 10$ km/s entre un impacteur avec un diamètre $D_{imp} = 100$ km et un corps de la taille de Mars (D'après Monteux and Arkani-Hamed (2016)).</figcaption></figure>

d'hydrocodes iSALE-2D pour calculer la pression de choc dans un corps différencié de type Mars pour des vitesses d'impact allant de 5 à 20 km/s, des rayons d'impacteur allant de 50 à 200 km et différentes rhéologies. Pour mieux représenter la distribution de la pression de choc en fonction de la distance du site d'impact à la surface, nous avons proposé deux régions distinctes dans le manteau : une région de champ proche qui s'étend jusqu'à 7-15 fois le rayon du projectile dans la cible, où la pression de choc maximale décroît de manière exponentielle avec l'augmentation de la distance du site d'impact, et une région de champ lointain où la pression décroît fortement avec la distance en suivant une loi de puissance. À la frontière entre le noyau et le manteau (CMB), la pression de choc maximale augmente du côté du manteau vers le côté du noyau. L'onde de choc réfractée se propage à l'intérieur du noyau où la pression de choc diminue en suivant une seconde loi de puissance.
## 2.3 Chauffage résultant d'un impact

L'énergie disponible lors d'un impact est l'énergie cinétique de l'impacteur. Ainsi elle est fonction de deux paramètres : la vitesse et la masse de l'impacteur. Lors de la collision, 10 à 20% est utilisée pour détruire l'impacteur, 40 à 50% de l'énergie d'impact permet de déformer le corps impacté et 30 à 40% est convertie en chaleur (O'Keefe and Ahrens, 1977). Au cours de la collision, une partie de l'énergie cinétique est donc dissipée sous forme de chaleur et permet une augmentation locale de la température en profondeur, en dessous du point d'impact dans une région sphérique à hémisphérique (Figure 4). L'évolution temporelle de la taille des impacteurs a probablement joué un rôle important dans l'évolution thermique primitive des planètes puisque la taille des impacteurs conditionne le taux d'accrétion mais aussi le volume de la région chauffée suite à un impact (Pierazzo et al., 1997).

Il est possible d'estimer au premier ordre l'augmentation de température en dessous du point d'impact en considérant que :

1. le volume chauffé est 3 fois le volume de l'impacteur

2. la vitesse de l'impacteur est au moins égale à la vitesse d'échappement de la planète impactée

3. une fraction $\gamma$ de l'énergie cinétique de l'impacteur est dissipée sous forme de chaleur ($\gamma \approx 30\%$) le reste étant utilisé pour déformer et fragmenter l'impacteur et sa cible.

On aboutit ainsi à l'estimation suivante :

$$ \Delta T \approx \frac{4\pi\gamma G\rho R^2}{9C_p} \qquad (1) $$

avec $G$ la constante gravitationnelle. Pour un corps de la taille de la Lune l'augmentation de température est d'environ 100 K alors que pour un corps de la taille de la Terre, cette augmentation de température est de l'ordre de 2000 K (Monteux et al., 2007).

Les vitesses d'impact sont généralement comprises entre 1 et 10 fois la vitesse d'échappement de la planète impactée (vitesse nécessaire pour échapper à l'attraction gravitationnelle du corps en croissance). Si les vitesses d'impact sont essentiellement fonction de la taille du corps impacté, l'apport d'énergie après impact augmentera donc au cours de l'accrétion (Figure 5). Les parties les plus externes des planètes en croissance seront ainsi plus susceptibles de fondre en premier (Kaula, 1979). Aussi lors des dernières étapes d'accrétion de grand volumes de manteau primitif sont susceptibles d'être chauffés par impact et de fondre si la température locale dépasse le solidus du manteau au cours de la croissance de la planète. Les énergies d'impact sont ainsi de plus en plus susceptibles de faire fondre le fer et les silicates du manteau. Cette énergie peut s'ajouter à la chaleur de désintégration radioactive des éléments à courtes périodes tels que $^{26}$Al ou $^{60}$Fe et entraîner des phénomènes de fusion significative au cours desquels se forment des océans magmatiques.
<figure><img src="image_5.png" /></figure>

FIGURE 4 – Augmentation de la température après un impact au sein d'une planète différenciée de la taille de Mars. La première ligne : $R_{imp} = 50\text{km}$ et $V_{imp} = 10\text{ km/s}$, et la deuxième ligne : $R_{imp} = 150\text{ km}$ et $V_{imp} = 10\text{ km/s}$ (gauche : modèle constitutif avec résistance, droite : modèle hydrodynamique sans résistance). Les quatre figures représentent la combinaison des lois d'échelle obtenues à partir de codes iSALE avec le modèle d'augmentation de température induite par l'impact de Watters et al. (2009). La ligne blanche en pointillés représente le CMB (D'après Monteux and Arkani-Hamed (2019)).

Grâce aux outils de modélisation disponibles (voir section 2.1), il est désormais possible d'estimer plus précisément la répartition du chauffage post-impact à l'intérieur du corps impacté. Dans l'étude évoquée précédemment (Monteux and Arkani-Hamed, 2019), à partir des résultats obtenus grâce aux simulations de l'hydrocode iSALE, nous avons dérivé des lois d'échelles qui illustrent l'influence de la rhéologie de la cible, de la taille de l'impacteur et de la vitesse de l'impact sur la propagation de l'onde de choc et la répartition de la pression au
cours d'un impact. En combinant les lois d'échelles obtenues avec un formalisme qui relie augmentation de température et pression de choc comme celui proposé par Watters et al. (2009), nous avons pu déterminer le chauffage induit par des impacts de diamètres supérieurs à 100 km au sein d'une planète différenciée. Nous avons ainsi pu comparer les champs de température post-impact obtenus avec l'hydrocode iSALE et les champs de température prédits par nos lois d'échelle. Ces résultats sont illustrés dans la Figure 4. Cette figure montre notamment l'importance de la taille de l'impacteur qui conditionne fortement le volume chauffé. Cette figure illustre aussi que la géométrie de la zone chauffée localisée en dessous du point d'impact n'est ni parfaitement hémisphérique (Tonks and Melosh, 1992) ni sphérique (Croft, 1982; Pierazzo et al., 1997).

<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Représentation de l'évolution topographique et thermique d'une planète en croissance à partir d'un objet tellurique de rayon 500 km (gauche) jusqu'à l'accrétion d'un planète de 3000 km de rayon. Après chaque impact la protoplanète grossit et une partie de l'énergie cinétique est dissipée sous forme de chaleur en profondeur (centre). À la fin de l'accrétion, les vitesses d'impact étant plus élevées, l'intensité du chauffage par impact est plus grande dans les parties externes de la planète (droite).</figcaption></figure>

Une fois que l'anomalie thermique est formée en dessous du point d'impact, un réajustement mécanique à plus long terme va avoir lieu afin que le matériel chaud et moins dense puisse se refroidir. Au cours de mes travaux de thèse (Monteux et al., 2007), je me suis intéressé à la dynamique de relaxation de l'anomalie thermique engendrée par un impact météoritique sans fusion. J'ai développé un code numérique permettant de résoudre les équations de Navier-Stokes en différences finies, en géométrie axisymétrique cylindrique et en formulation en fonction de courant (Monteux et al., 2007). Nos résultats montrent que si l'impacteur est trop petit, l'anomalie thermique va seulement se refroidir par diffusion. Cependant, lorsque l'impacteur dépasse une taille critique, l'anomalie thermique subit d'abord une phase advective
isotherme puis refroidit par diffusion. L'épaisseur de l'anomalie thermique à la fin de l'étape advective ne dépend que des paramètres rhéologiques du corps impacté.

## 2.4 Impacts géants et fusion à grande échelle

Au cours des dernières étapes de la croissance des planètes telluriques, les collisions entre protoplanètes ne sont pas des évènements exceptionnels (Agnor et al., 1999). Les impacts géants sont maintenant connus pour être inévitables durant les étapes tardives de l'accrétion des planètes et ont pu avoir des conséquences majeures sur l'évolution tardive des planètes telluriques. Nous savons que certaines caractéristiques du Système Solaire peuvent être expliquées par les impacts géants. Les impacts géants sont ainsi invoqués pour expliquer la faible épaisseur du manteau de Mercure par rapport à la taille de son noyau (Benz et al., 1988), la dichotomie martienne (Marinova et al., 2008) ou encore les moments angulaires actuels des planètes (Kurosaki and Inutsuka, 2019).

La formation de la Lune est probablement due à un ou plusieurs impacts géants durant l'histoire primitive de la Terre. La faible taille du noyau métallique lunaire par rapport à l'épaisseur de son manteau et le fait que la Lune et la Terre présentent des similitudes de composition isotopique notamment en oxygène (Wiechert et al., 2001) supportent l'idée que, tôt dans l'histoire de la Terre, à un moment où la Terre était déjà différenciée (70-100 Ma après la formation des C.A.I.), la Terre a été impactée par un objet avec un rayon de plusieurs milliers de kilomètres (Théia). Le champ de débris suivant cet impact, essentiellement formé de matériel terrestre mantellique et donc appauvri en fer, s'est ensuite ré-accrété autour de la Terre permettant la formation de la Lune en orbite autour de la Terre. Les noyaux de l'impacteur et de la Terre se sont mélangés pour former le noyau terrestre (Hartmann and Davis, 1975). Les études du moment angulaire du système Terre-Lune menées par Cameron and Ward (1976) ont ensuite permis d'arriver à la même conclusion en donnant un ordre de grandeur de la taille de l'impacteur : entre 1/3 et 1/2 rayon terrestre. L'impact ayant formé le système Terre-Lune aurait eu lieu avec une vitesse de l'ordre de quelques km.s⁻¹. Cette théorie a inspiré de nombreux travaux de modélisation de cet impact notamment à l'aide des codes S.P.H. (Smoothed Particle Hydrodynamics) développés par Canup and Asphaug (2001). La similarité de composition isotopique en oxygène (Wiechert et al., 2001) et en tungstène (Touboul et al., 2007) laisse entrevoir deux possibilités. Soit la Lune s'est formée essentiellement à partir de matériel terrestre, soit les manteaux lunaire et terrestre se sont équilibrés à la suite de l'impact géant. Or, dans les modèles dynamiques (Canup, 2004), la majorité du matériel lunaire proviendrait de l'impacteur.

Pour expliquer la composition de la Lune et sa similarité isotopique avec la Terre, Lock et al. (2018) ont proposé un modèle faisant intervenir un impact géant, plus énergétique que dans le modèle traditionnel. Cet impact entraînerait la Terre dans un état de vaporisation et de
rotation rapide qui s'étend sur des dizaines de milliers de kilomètres. De tels états planétaires sont appelés "synestias". Lorsque la synestia se refroidit, la matière se condense et forme la Lune. Lock et al. (2018) ont constaté que la Lune se forme à l'intérieur de la synestia, entourée de vapeur de composition terrestre à des pressions de plusieurs dizaines de bars. La Lune orbite dans la synestia suffisamment longtemps pour s'équilibrer chimiquement avec la Terre vaporisée ce qui favoriserait les similitudes chimiques entre la Terre et la Lune.

Outre la formation de la Lune, cet impact géant a probablement eu des conséquences ma-
jeures sur l'état thermique du manteau terrestre. Le timing de l'impact (moins de 100 Ma après
la formation des premiers solides) suppose qu'il a eu lieu sur une planète très chaude déjà si-
gnificativement affectée par les précédents impacts, le chauffage radioactif et la formation du
noyau (dont je discute dans la section suivante). Selon les modèles de Nakajima and Stevenson
(2015), l'impact qui a formé le système Terre-Lune a entraîné la fusion complète du manteau
terrestre et des processus de mélanges plus ou moins importants selon la taille de l'impacteur
et la vitesse de rotation de la Terre avant impact.

L'idée que l'impact à l'origine du système Terre-Lune ait fondu et mélangé efficacement
le manteau primitif est difficile à justifier géo-chimiquement notamment si on s'intéresse au
système Hf-W. Les roches terrestres présentent en effet des rapport variés en ¹⁸²W/¹⁸⁴W no-
tamment dans les roches anciennes (Lee and Halliday, 1996; Willbold et al., 2011; Touboul
et al., 2012). Cette variété pourrait indiquer que le manteau primitif était hétérogène et que
la signature de cette hétérogénéité a été préservée pendant plusieurs milliard d'années. Ainsi,
moins de 1000 ans après l'impact géant (ce qui correspond au temps d'accrétion de la Lune), le
manteau se résumerait à un océan magmatique de 2900 km de profondeur mais conservant po-
tentiellement des hétérogénéités chimiques (Nakajima and Stevenson, 2015). La préservation à
long terme de telles hétérogénéités indiquerait que la vigueur convective du manteau au cours
des temps géologiques n'a pas été suffisante pour détruire tous les réservoirs mantelliques qui
se sont formés au début de l'histoire de la Terre (Jackson et al., 2010). Des modèles dyna-
miques du manteau terrestre suggèrent que des parties vierges du manteau pourraient échapper
à la différenciation au cours de l'âge de la Terre grâce à un isolement dans des tourbillons
convectifs (Brandenburg et al., 2008). La question du mélange qui suit l'impact géant (1000
ans après l'impact) doit cependant être envisagée grâce à des modèles de convection adaptés à
la caractérisation de processus turbulents (Solomatov, 2000)(voir Section 6).

2.5 Autres sources de chaleur

En plus du chauffage par impact, d'autres sources de chaleur ont pu participer au budget
thermique primitif de la Terre. Ainsi, le chauffage par les éléments radioactifs à courte période
ou la dissipation thermique liée à la formation du noyau terrestre ont pu jouer un rôle essentiel
dans l'évolution des jeunes planètes.
### 2.5.1 Chauffage par les éléments radioactifs à courte période

La chaleur dégagée par la désintégration des éléments radioactifs à courtes périodes a probablement joué un rôle important sur les températures des planétésimaux primitifs. En effet, la désintégration radioactive de certains éléments s'accompagne d'un dégagement d'énergie sous forme de rayonnement (γ) ou de l'émission de particules α ou β dont la diffusion dans les matériaux environnants dégage d'importantes quantités de chaleur. Les éléments isotopiques (parents) à courtes demi-vies ont maintenant disparu mais un excès d'éléments (fils) issus de cette désintégration peut se retrouver dans les météorites intégrant des C.A.I.. L'étude de ces radioactivités maintenant "éteintes" dans les météorites donne des informations importantes sur l'environnement initial de formation des premiers corps et des étoiles (Meyer and Clayton, 2000).

<table><caption>TABLE 1 – Couples d'isotopes radioactifs ayant contribué au bilan thermique du manteau terrestre (selon Lebrun et al. (2013) et McDonough et al. (2020)).</caption><thead><tr><td>Élément père</td><td>τ<sub>1/2</sub>(Myr)</td><td>Élément fils</td><td>Production de chaleur par unité de masse (W/kg)</td></tr></thead><tbody><tr><td><sup>26</sup>Al</td><td>0.7</td><td><sup>26</sup>Mg</td><td>1.5 × 10<sup>-7</sup></td></tr><tr><td><sup>238</sup>U</td><td>4468.3</td><td><sup>206</sup>Pb</td><td>9.17 × 10<sup>-5</sup></td></tr><tr><td><sup>40</sup>K</td><td>1248</td><td><sup>40</sup>Ar</td><td>2.97 × 10<sup>-5</sup></td></tr><tr><td><sup>232</sup>Th</td><td>14010</td><td><sup>208</sup>Pb</td><td>2.56 × 10<sup>-5</sup></td></tr><tr><td><sup>60</sup>Fe</td><td>1.5</td><td><sup>60</sup>Ni</td><td>0.3 – 1.2 × 10<sup>-8</sup></td></tr></tbody></table>

Si l'on considère *H* le taux de production interne de chaleur par désintégration isotopique d'un planétésimal de rayon *R*. Ce chauffage prend le pas sur la diffusion thermique lorsque $\rho C_p H$ est supérieur à $kT/R^2$ avec *k* la conductivité thermique du planétésimal, *T* sa température, $C_p$ sa capacité calorifique et $\rho$ sa densité. $H \sim T/t$ donc lorsque $\rho C_p R^2/(kt) \gg 1$, le planétésimal se réchauffe plus vite qu'il ne se refroidit par diffusion avec *t* l'âge du planétésimal. En 1 Ma et pour des objets avec des rayons inférieurs à 5 km, le refroidissement par diffusion était prépondérant. En 1 Ma, selon Kokubo and Ida (2000), les planétésimaux les plus volumineux auraient atteint des rayons de ~ 10 – 100 km. La chaleur associée à la désintégration était donc globalement retenue. La température maximale atteinte au centre d'un planétésimal de rayon *R* = 30 km aurait pu dépasser 1270 K (Schramm et al., 1970; Hutcheon and Hutchison, 1989). Les plus gros objets auraient ainsi pu subir une éventuelle différence associée à des phénomènes de fusion via la désintégration d'isotopes radioactifs à courte période puis évacuer leur chaleur par diffusion. Au moment où les protoplanètes commencent à atteindre des tailles importantes (~ 1000 km), elles ont évacué par diffusion une partie importante de leur chaleur issue de la désintégration radioactive des éléments à courte période.

Deux de ces radioactivités éteintes ont pu participer à l'état thermique initial des planétésimaux (Tableau 1) : la désintégration du <sup>26</sup>Al (principalement localisé dans la phase silicatée des planétésimaux) et la désintégration du <sup>60</sup>Fe (principalement localisé dans la phase métallique des planétésimaux). Des excès importants de <sup>26</sup>Mg ont été trouvés confinés dans des
C.A.I. dans des météorites carbonnées et dans chondrites ordinaires mais aussi dans du matériel météoristique non réfractaire (Hutcheon and Hutchison, 1989). Cette présence de $^{26}$Mg en excès traduit la présence passée de $^{26}$Al. Le $^{26}$Al était ainsi présent en quantité significative pendant les premiers millions d'années du système solaire.

<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Évolution temporelle du chauffage radioactif pour les éléments chimiques qui ont le plus contribué au bilan thermique de la Terre depuis sa formation il y a plus de 4.5 Ga (à partir de McDonough et al. (2020)).</figcaption></figure>

Il est possible d'estimer la puissance dissipée par ce mode de chauffage en fonction du temps (Figure 6). Cette puissance est fonction de la composition chimique considérée en conditions initiales. McDonough et al. (2020) estiment qu'au cours des 20 premiers millions d'années après la formation des C.A.I., le $^{60}$Fe et le $^{26}$Al ont été les éléments qui ont le plus contribué au chauffage radioactif dissipant une puissance allant jusqu'à $10^6$ TW. Aujourd'hui, le chauffage radioactif est principalement fourni par la désintégration de $^{238}$U, $^{40}$K et $^{232}$Th qui permettent de produire environ 20 TW de chauffage en profondeur principalement dans la lithosphère.
### 2.5.2 Dissipation thermique lors de la formation du noyau

Une autre source majeure de chaleur ayant fortement influencé l'état thermique primitif du manteau terrestre provient de la formation de son noyau. En effet si l'on imagine que le fer du noyau était initialement réparti uniformément au sein de la Terre après l'accrétion, la séparation fer/silicates et la chute du fer vers le noyau provoquent un transfert de masse important du fait que le fer du noyau est plus dense que les silicates du manteau. Ce transfert de masse diminue l'énergie potentielle du fer dans le manteau qui se convertit en chaleur à l'intérieur de la jeune Terre. Flasar and Birch (1973) estiment que cette conversion d'énergie potentielle en chaleur libérerait suffisamment d'énergie pour provoquer une augmentation de la température à l'intérieur de la Terre de plus de 2000 K.

La répartition de cette chaleur à l'intérieur de la Terre dépend fortement de la façon dont le fer va migrer vers son centre (Ke and Solomatov, 2009). Si le fer plonge en suivant de long chenaux, la chaleur ira principalement dans le noyau. Cependant si le fer plonge sous la forme de grosses gouttes ou de larges diapirs, l'énergie potentielle se dissipera principalement dans le manteau. Dépendamment des mécanismes de ségrégation, cela donnera lieu à des profils thermiques différents avec soit un noyau bien plus chaud que le manteau, soit un manteau très chaud voir fondu. Les modalités de différenciation et notamment les volumes caractéristiques de fer (gouttelettes centimétriques ou diapir pluri-kilométriques) mis en jeu lors de la séparation fer/silicates conditionnent donc l'état thermique primitif du manteau terrestre. Samuel et al. (2010) ont aussi montré qu'un grand nombre de petits diapirs de fer plongeant dans un océan magmatique favorise l'échange de chaleur fer/silicates et conduit à un noyau relativement chaud comme sur la Terre. À contrario, un petit nombre de grands diapirs a l'effet inverse et pourrait correspondre à une planète semblable à Mars. La séparation fer/silicates par diapirisme tend à laisser le manteau profond nettement plus chaud que sa partie supérieure, ce qui pourrait provoquer la fusion d'une fraction significative au dessus de la frontière manteau/noyau. Nous reviendrons plus en détail sur ce processus dans la section 3.

### 2.5.3 Chauffage par dissipation de marées

L'impact à l'origine de la formation du système Terre-Lune a probablement induit la formation d'un océan de magma terrestre profond et généré un disque de débris autour de la proto-Terre qui s'est rapidement ré-accrété pour former la Lune. À cette époque, la Lune était plus proche de la Terre qu'aujourd'hui, ce qui a pu avoir des conséquences importantes sur l'interaction entre les deux objets. En particulier, les interactions de marée qui ont suivi la formation de la Lune ont contrôlé l'évolution de la rotation et ont potentiellement influencé le refroidissement de la Terre et de la Lune primitives (Chen and Nimmo, 2016).

L'évolution de l'orbite lunaire dépend fortement de la manière dont l'énergie des marées a été dissipée sur la Terre primitive et est donc étroitement liée à son évolution thermique. La friction des marées associée aux mouvements de cisaillement dissipe une énorme quantité
d'énergie. Sur la Terre moderne, la dissipation des marées est dominée par les processus de dissipation océanique (Ray et al., 2001). Actuellement, la dissipation associée aux marées représente 3,5 TW, mais il y a 4,5 Ga, la puissance produite par la friction des marées pouvait être supérieure de plus de 2 ordres de grandeur (Farhat et al., 2022). Au cours de l'Hadéen, la fraction fondue du manteau silicaté de la Terre devait avoir un impact important sur les processus de dissipation des marées et, ainsi, sur l'ensemble de l'évolution du système Terre-Lune. En retour, l'énergie dissipée dans la Terre hadéenne a pu affecter sa structure thermique, sa dynamique convective et étendre la fraction interne en matériel fondu, conduisant à un régime thermique induit par les marées comparable à ce qui est observé sur Io, une des lunes de Jupiter (Kervazo et al., 2021). Même après la première étape de cristallisation de l'océan magmatique, une quantité importante de matière fondue peut être préservée dans le manteau hadéen (Ray et al., 2001) ce qui peut contrôler à la fois la production de chaleur par friction tidale et le transfert de chaleur par extraction de la matière fondue. Cette période de dissipation intense des marées (dépassant plusieurs centaines de TW) a pu durer pendant les premières dizaines de millions d'années suivant l'impact formant la Lune, ce qui pourrait avoir eu une influence significative sur la dynamique de refroidissement de la Terre partiellement fondue et sur la séquence de cristallisation.

Si l'état orbital et rotationnel récent du système Terre-Lune est bien documenté (Williams, 2000), les interactions du système Terre-Lune à ses débuts n'ont pas été étudiées et, en particulier, l'influence de la Lune primitive sur l'évolution thermique et la fusion de la Terre primitive reste encore inexplorée. Nous reviendrons plus en détails sur ce processus dans la section 6.4.

## 2.6 Perspectives : Le rôle des impacts dans l'évolution primitive de l'atmosphère terrestre (Thèse de T. Engels)

Les impacts météoritiques ont aussi participé à la formation de l'atmosphère terrestre :

1. en érodant l'atmosphère préexistante : Le long de sa trajectoire dans l'atmosphère, l'impacteur forme un canal de gaz chaud moins dense que l'atmosphère environnante. Ce gaz chaud migre à travers ce canal et est éventuellement éjecté dans l'espace (Shuvalov, 2009).

2. en apportant du matériel au corps impacté (Raymond et al. (2006), Engels et al. (in revision for Icarus)) : Ce matériel peut être soit dégazé lors de son impact à la surface du corps en croissance soit être dégazé ultérieurement après son dépôt sous forme solide à la surface du corps impacté.

3. en dégazant du matériel profond : Un impact entraîne la formation d'une atmosphère transitoire à haute température et riche en silicate, mais aussi la fusion de matériel rocheux dans lequel du gaz est emprisonné (Figure 7) (Marchi et al., 2016). Lors de son refroidissement rapide, l'atmosphère transitoire se condense, tandis que la masse fondue générée par l'impact se propage à la surface. Lorsque la roche fondue se propage, elle libère des gaz (CO, CO₂, H₂, SO₂, CH₄, H₂O, H₂S). Des études en laboratoire ont
permis de quantifier la quantité de volatils dégazés sur des matériaux impactés comme de la calcite ou de la serpentine (Lange and Ahrens, 1982, 1986). Ces études montrent notamment que la serpentine commence à dégazer son H₂O à des pressions de choc d'environ 20 GPa ou que la calcite commence à dégazer son CO₂ à 10–15 GPa. Au contact de la surface, un panache de gaz peut se développer au dessus du point d'impact où les molécules de gaz sont aussi chauffées et peuvent être accélérées vers le haut et échapper à la gravité du corps impacté.

La compétition de ces 3 phénomènes peut ainsi entraîner l'augmentation de l'épaisseur de l'atmosphère ou son érosion. L'efficacité du dégazage ou de l'érosion par impact est fonction de la vitesse d'impact ($v_{imp}$), de la taille de l'impacteur ($R_{imp}$), de la taille du corps impacté ($R_p$), de la composition chimique et des propriétés rhéologiques (porosité, viscosité) des matériaux mis en jeu.

<figure><img src="image_8.png" /><figcaption>FIGURE 7 – Représentation schématique d'un épisode de fusion et du dégazage provoqué par impact sur la Terre primitive (À partir de Marchi et al. (2016)).</figcaption></figure>

Le dégazage des éléments volatils est encore peu contraint et attribué à l'action combinée des hautes températures, à une amplification de la diffusion ou au développement de micro-fractures (Zhang, 2014). Comme la combinaison de ces mécanismes est peu connue, les lois d'échelles permettant d'obtenir la fraction de matériel dégazé par impact géant à partir d'expériences en laboratoires sur de petits volumes sont difficiles à obtenir. De plus, il existe encore de grandes incertitudes sur l'histoire accrétive des planètes et de la Terre notamment. Grâce à la modélisation numérique on connaît désormais mieux les conséquences internes des gros impacts météoritiques (Monteux and Arkani-Hamed, 2016, 2019). Les lois d'échelles obtenues à partir de modélisations hydrocodes montrent que la pression en dessous du point d'impact peut atteindre plusieurs dizaines de GPa sur des volumes comparables à celui de l'impacteur. Ensuite l'augmentation de pression subie lors de l'impact diminue rapidement en s'éloignant du point d'impact. Les paramètres d'accrétion ont ainsi pu jouer un rôle majeur sur l'évolution de l'épaisseur de cette atmosphère (Ahrens, 1993; Ahrens et al., 2004). Cependant, à l'heure actuelle, aucun modèle quantitatif ne permet d'intégrer l'effet du dégazage du manteau, des gaz initiaux et des composés volatils produits par l'impact lors de l'accrétion des planètes (Zhang, 2014).

Si globalement la taille des impacteurs va augmenter au cours de l'histoire primitive du
Système Solaire, cette taille n’est pas unique et obéit à une distribution où les gros corps sont plus rares que les petits. L’accrétion de petits planétésimaux va favoriser la formation d’une atmosphère alors que les impacts géants vont entrainer un dégazage massif mais surtout une érosion massive de l’atmosphère préexistante voir éventuellement sa disparition. Un impacteur de la taille de la Lune serait ainsi susceptible d’ôter 50% de l’atmosphère terrestre (Ahrens, 1993; Ahrens et al., 2004). Les objectifs de la thèse de Tristan Engels sont de contraindre l’évolution primitive de l’atmosphère terrestre lors des dernières étapes de son accrétion (i.e. après l’impact géant à l’origine du système Terre-Lune) et en particulier de :

1. Contraindre l'influence des paramètres d'accrétion ($v_{imp}$, $R_{imp}$, taux d'accrétion) sur l'évolution des atmosphères planétaires. La taille des impacteurs va gouverner non seulement la quantité de matériel potentiellement dégazé par impact mais aussi la fraction d'atmosphère érodée. La répartition temporelle des tailles d'impacteurs et des leurs vitesses a donc surement joué un rôle sur l'épaisseur finale de l'atmosphère primitive à la fin des grandes étapes d'accrétion. L'évolution de la taille du corps en croissance et ainsi son taux d'accrétion vont aussi contrôler sa capacité à retenir son atmosphère.

2. Quantifier l'importance de l'accrétion par rapport au dégazage profond. Au cours de son histoire primitive, la Terre a pu perdre son atmosphère de manière répétée. Les derniers impacts géants qui auraient chassé tout ou partie de l'atmosphère pourraient être le point de départ de la formation de l'atmosphère terrestre par dégazage du manteau ou dégazage par impact (Zhang, 2014). En apportant des contraintes sur les caractéristiques (épaisseur et composition) de l'atmosphère primitive et en comparant avec les caractéristiques de l'atmosphère actuelle, il est possible de contraindre l'importance du dégazage généré lors de l'accrétion par rapport au dégazage profond lié à des phénomènes plus tardifs comme le volcanisme ou le bombardement tardif « Late Heavy Bombardment » (LHB).

Nous testerons si un dégazage sélectif est un modèle viable pour expliquer le spectre des éléments volatils (C, H, N et S) (Marty, 2012) ou si d'autres processus doivent être pris en compte (formation du noyau, la nature des briques élémentaires ou bien l'érosion de la surface terrestre via les impacts). Il est à noter que ce projet est très complémentaire des expériences menées au LMV pour déterminer le partage des volatils entre le noyau et le manteau (ex. Clesi et al. (2016, 2018); Boujibar et al. (2014)). Nous développons actuellement avec le Laboratoire de Météorologie Physique de Clermont-Ferrand (L. Deguillaume, A. Bianco, L. Pailler) un modèle numérique dans lequel il sera possible de choisir les caractéristiques de l'accrétion de la Terre en termes de taille d'impacteurs, de vitesse d'impacts et de taux d'accrétion et de connaître les caractéristiques de l'atmosphère primitive générée par une telle accrétion. Les informations obtenues par ces modèles seront l'épaisseur mais aussi la composition de cette atmosphère en éléments tels que C, H, N, O et S.

Dans ce cadre, l'expertise du LaMP permettra d'aider à la compréhension des mécanismes de réactivité chimique avec l'objectif d'estimer l'évolution temporelle de l'atmosphère primitive générée par l'accrétion. Les composés chimiques (CO, $CO_2$, $CH_4$, $H_2S$, $SO_2$, ...) émis sont
des précurseurs à la formation de nouvelles espèces chimiques après réactions chimiques. Ces
processus chimiques de transformations dépendent des conditions atmosphériques (pression,
température,...) et du flux actinique solaire initiateur potentiel de réactions photochimiques.
Des mécanismes chimiques devront être construits sur la base des précurseurs chimiques en
évaluant les paramètres cinétiques importants (constantes de réactivité et leur dépendance à
la température et à la densité de l'air, par exemple). Les travaux de recherche sur la réactivité
chimique des panaches volcaniques pourront également servir à la création de ces mécanismes
car les conditions de températures et les précurseurs sont similaires. Des mécanismes d'oxyda-
tion pourraient également aboutir à la formation de composés moins volatils et à la formation
potentielle de particules solides (nucléation). Une fois ces transformations mieux appréhen-
dées, des estimations de temps de vie de ces composés chimiques pourront être évaluées afin
d'estimer l'évolution temporelle de l'atmosphère primitive. Les codes numériques de chimie
atmosphérique développés par le LaMP sont en cours d'adaptation par T. Engels afin de simu-
ler l'évolution chimique de l'atmosphère primitive.

Références personnelles

Monteux, J., Coltice, N., Dubuffet, F., Ricard, Y., 2007, Thermo-mechanical adjustment after impacts during planetary growth, Geophys. Res. Lett., 34, 24201-24205. LIEN HAL

Monteux, J., Arkani-Hamed, J., 2016b, Scaling laws of impact induced shock pressure in planetary mantle, Icarus, 264, 246–256. LIEN HAL

Monteux, J., Arkani-Hamed J., 2019, Shock wave propagation in layered planetary interiors : Revisited. Icarus, 331, 238-256. LIEN HAL

Engels T., Monteux, J., Boyet M., Bouhifd M.A., Large impacts and their contribution to the water budget of the Early Moon, in revision for Icarus.
<figure><img src="image_9.png" /></figure>

# Thermo-mechanical adjustment after impacts during planetary growth

Julien Monteux,<fnref n="1" /> Nicolas Coltice,<fnref n="1" /> Fabien Dubuffet,<fnref n="1" /> and Yanick Ricard<fnref n="1" />

Received 29 August 2007; accepted 7 November 2007; published 18 December 2007.

[1] The thermal evolution of planets during their growth is strongly influenced by impact heating. The temperature increase after a collision is mostly located next to the shock. For Moon to Mars size planets where impact melting is limited, the long term thermo-mechanical readjustment is driven by spreading and cooling of the heated zone. To determine the time and length scales of the adjustment, we developed a numerical model in axisymmetric cylindrical geometry with variable viscosity. We show that if the impactor is larger than a critical size, the spherical heated zone isothermally flattens until its thickness reaches a value for which motionless thermal diffusion becomes more effective. The thickness at the end of advection depends only on the physical properties of the impacted body. The obtained timescales for the adjustment are comparable to the duration of planetary accretion and depend mostly on the physical properties of the impacted body. Citation: Monteux, J., N. Coltice, F. Dubuffet, and Y. Ricard (2007), Thermo-mechanical adjustment after impacts during planetary growth, Geophys. Res. Lett., 34, L24201, doi:10.1029/2007GL031635.

## 1. Introduction

[2] Impacts have strongly influenced the evolution of planets: a collision of the Earth with a Mars-sized body is at the origin of the formation of the Moon [Hartmann and Davis, 1975] and the impact by a kilometer-sized body could be responsible for the mass extinction at the K-T boundary [Alvarez et al., 1980]. It is during accretion that impacts played the most significant role, depositing and burying heat into growing planetary bodies.

[3] When the impact velocity becomes larger than the elastic velocities, a shock wave develops. The shock pressure, increasing with the size of the impacted body, is nearly uniform in a spherical region next to the impact (the isobaric core), and strongly decays away from it [Croft, 1982]. Following the adiabatic pressure release, the peak pressure being independent of impactor size, a temperature increase of several hundred degrees remains on Moon to Mars size bodies [Senshu et al., 2002] (see equation (2)). Hence, the hotter temperatures are located close to the surface during planetary growth [Kaula, 1979] and large impacts have caused extensive melting and formation of magma oceans on Earth [Tonks and Melosh, 1993].

[4] The thermal anomaly caused by an impact generates a buoyant thermal anomaly that ultimately drives an isostatic adjustment. If the impact velocity is larger than 7.5 km.s⁻<fnref n="1" />, a significant volume of the isobaric core is molten [O'Keefe

and Ahrens, 1977] hence the adjustment is controlled by two-phase flow and probably hydrofracturation [Solomatov, 2000]. For smaller planets or planetesimals, melting is nearly absent therefore the thermo-mechanical adjustment is dominated by the slow viscous deformation and thermal diffusion of the hot thermal anomaly.

[5] In this study, we investigate the thermal relaxation and viscous deformation after the shock of an impactor on a small planet or planetesimal in order to derive scalings for the relevant length and time scales of the thermo-mechanical adjustment.

## 2. Model Description

### 2.1. Thermal State After an Impact

[6] Energy balance calculations and shock simulations suggest that the radius of the isobaric core $R_{ic}$ is comparable or slightly larger than that of the impactor $R_{imp}$ and we use $R_{ic} = 3^{1/3} R_{imp}$ [Senshu et al., 2002; Pierazzo et al., 1997]. Away from the isobaric core, the shock wave propagates and the peak pressure decays with the square of the distance $r$ from the center of the isobaric core [Pierazzo et al., 1997]. Just after the adiabatic pressure release, the thermal perturbation corresponds to an isothermal sphere of radius $R_{ic}$ and temperature $T_0 + \Delta T$ that decays when $r > R_{ic}$ as

$$T(r) = T_0 + \Delta T \left(\frac{R_{ic}}{r}\right)^m, \quad (1)$$

with $m \sim 4.4$ as proposed by Senshu et al. [2002].

[7] The energy dissipated as heat following the shock is a fraction of the kinetic energy of the impactor. The impactor velocity $v_{imp}$ should be comparable to the escape velocity $v_{imp} = \sqrt{2gR}$, where $g = 4/3 G\pi\rho R$, $\rho$ and $R$ are the gravity, density and radius of the impacted growing planet [Kokubo and Ida, 1996]. Assuming $\rho \sim \rho_{imp} \sim \rho_{ic}$, the temperature increase $\Delta T$ is

$$\Delta T = \frac{4\pi}{9} \frac{\gamma}{f(m)} \frac{\rho GR^2}{C_p}, \quad (2)$$

where $C_p$ is the heat capacity of the impacted body and $G$ is the gravitational constant. The efficiency of kinetic to thermal energy conversion $\gamma$ is close to 0.3 according to physical and numerical models [O'Keefe and Ahrens, 1977]. The function $f(m)$ represents the volume effectively heated normalized by the volume of the isobaric core (i.e., $f(m) = 1$ if only the isobaric core is heated). Assuming $R_{ic} \ll R$ and integrating equation (1) leads to $f(m) \sim 2.7$ and 37% of the impact heating is released within the isobaric core. The temperature increase does not depend on the size of the impactor but on the square of the radius of the impacted body.

<footnote n="1">¹Laboratoire de Sciences de la Terre, UMR 5570, CNRS, Université de Lyon, Université Claude Bernard Lyon 1, Ecole Normale Supérieure de Lyon, Villeurbanne, France.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 1.** Close up view of the non-dimensional temperature fields and streamlines computed for a temperature dependent viscosity at non-dimensionalized times (top) $t^* = 0$, (middle) $t^* = 3$ and (bottom) $t^* = 10$ (we use $\lambda = 10^{-2}$ and $Ra_{ic} = 10^2$ (see 2.2)).</figcaption></figure>

[8] Immediately after the shock, a fraction of the isobaric core is removed during crater excavation. However, it is only for small impactors (less than 5 km of radius) and large planets (more than 3000 km of radius) that a significant fraction of the heated zone is redistributed [Maxwell, 1977]. Because our study is focused on large impacts (>10km) on small planets (<3000km), we can neglect the excavated heat or consider it as an additional uncertainty on the parameters $\gamma$ or $f(m)$ (see equation (2)). Following excavation, reduction

of the crater depth to at most few kilometers takes place through isostatic rebound. The total amount of vertical motion that occurs during the isostatic process (a few km), and its duration (typically a few $10^3$ y) are negligible compared to the length scales and timescales of thermal readjustment. After the crater formation, most of the initial energy still remains below the crater [Turtle et al., 2003]. Modelling these processes of mass and energy redistribution is beyond the scope of this paper since we are interested in the long term consequences of shock heating as in work by Reese et al. [2002].

[9] The proposed thermal state following an impact sketched in Figure 1 (top) is that of a cold body of homogeneous temperature $T_0$ perturbed by the impact of a sphere of radius $R_{ic}$ (the depth of the crater being negligible compared to the heated zone).

## 2.2. Thermo-Mechanical Model

[10] The governing non-dimensional equations for the extremely viscous flow of a cooling hot drop are

$$ -\nabla P^* + \vec{\nabla} \cdot \left( \frac{\eta(T^*)}{\eta_0} \vec{\nabla} \vec{v}^* + \left[ \frac{\eta(T^*)}{\eta_0} \vec{\nabla} \vec{v}^* \right]^T \right) + T^* \vec{e}_z = 0, \quad (3) $$

$$ \frac{\partial T^*}{\partial t^*} = \frac{\nabla^2 T^*}{Ra_{ic}} - \vec{v}^* \cdot \nabla T^*, \quad (4) $$

$$ \vec{\nabla} \cdot \vec{v}^* = 0, \quad (5) $$

where distances, temperature and velocity are normalized by $R_{ic}$, $\Delta T$ and the characteristic Stokes velocity $v_s$ of the isobaric core

$$ v_s = \frac{\alpha \rho g \Delta T R_{ic}^2}{\eta_0}, \quad (6) $$

where $\eta_0$ is the viscosity far from the impact and $\alpha$ the thermal expansivity of the impacted body. $Ra_{ic}$ is the Rayleigh number based on the isobaric core radius:

$$ Ra_{ic} = \frac{\alpha \rho g \Delta T R_{ic}^3}{\kappa \eta_0}. \quad (7) $$

We define a Rayleigh number based on the size of the isobaric core $R_{ic}$ since in all our experiments, the radius of the planet R remains much larger than $R_{ic}$ and thus does not affect the dynamics except through the gravity and the temperature increase (see equation (2)).

[11] For planets of Moon to Mars size the gravity and the temperature increase are not very large (e.g. $g \simeq 3$ m s$^{-2}$, $\Delta T \simeq 300$ K). We also consider impactors with radius small compared to the planet radius (e.g. $R_{ic} \simeq 300$ km). In this case, $Ra_{ic}$ should remain moderate (say lower than $10^5$) assuming that the coldest material of the growing planetesimal have a viscosity $\eta_0$ comparable to that of the present day Earth, say around $10^{21}$ Pa s (e.g. $Ra_{ic} \simeq 4700$ for $\kappa = 10^{-6}$ m$^2$s$^{-1}$, $\alpha = 5 \times 10^{-5}$ K$^{-1}$, $C_p = 1200$ J K$^{-1}$ kg$^{-1}$, $\rho = 3870$ kg m$^{-3}$). The viscosity is temperature-dependent
<figure><img src="image_11.png" /><figcaption>**Figure 2.** (top) Thickness $a^*/2$ (black solid line) and radial extent $R^*$ (red solid line) and (bottom) maximal temperature at $r=0$ as functions of time $t^*$ for $Ra_{ic} = 10^2$. Power-law fits following a diffusive solution are depicted by dashed lines. The equilibrium between diffusive (gray field) and advective (white field) stages is obtained at $t^* = t_{adv}^*$.</figcaption></figure>

$\eta(T^*) = \eta_0\lambda^{T^*}$, $\lambda$ being the viscosity ratio (lower than 1) between the hottest ($T^* = 1$) and the coldest ($T^* = 0$) material. This viscosity decreases sharply with temperature and its expression is simpler to implement than the usual Arrhenius law.

[12] We developed and benchmarked against analytical solutions a finite difference code to solve equations (3) and (5) in axisymmetric cylindrical geometry. We use a stream function formulation for the Stokes equation with a direct implicit inversion method [Schubert et al., 2001]. The heat equation (4) is solved by an Alternating Direction Implicit (ADI) scheme. All the equations are discretized with a second order accuracy in space and time. Advection terms are calculated with a centered differences scheme. This scheme is known to produce numerical oscillations for grid Peclet numbers larger than $\simeq 20$ (the grid Peclet number is $Ra\Delta x^2 v^*/\Delta t$ where $\Delta x$ is the grid dimension and $\Delta t$ the time step) [Dubuffet et al., 2000] but in our simulations, the maximum Peclet number remains of order unity. We used at least 129 $\times$ 129 grid points and the dimensions of the calculation domain varies between $10R_{ic} \times 10R_{ic}$ and $100R_{ic} \times 100R_{ic}$. Boundary conditions are free-slip, isothermal at the surface and insulating on other walls. In order to keep a good accuracy in a large enough domain and with strong viscosity variations we run the simulations up to

$Ra_{ic} = 10^4$ only, from which we fit scaling laws that can be extrapolated to higher Rayleigh numbers.

[13] The geometrical evolution of the post-impact thermal anomaly as a function of time is monitored by its non-dimensional radial extent $R^*(t^*)$, its thickness $a^*(t^*)$ and its maximum temperature $T_{max}^*(t^*)$. $a^*(t^*)$ is the depth where the second derivative of the vertical temperature profile at $r^* = 0$ is zero. Along this profile the maximum temperature value $T_{max}^*(t^*)$ is reached at $z^* = z_{max}^*$. $R^*(t^*)$ is the distance where the second derivative of the horizontal temperature profile at $z^* = z_{max}^*$ is zero.

### 3. Results

[14] For large enough Rayleigh numbers, the thermal relaxation consists in an early advective stage corresponding to an isothermal flattening of the hot drop, followed by a later stage of diffusive cooling. For $Ra_{ic} \le 4.9$, cooling is motionless.

#### 3.1. Advective Stage and Diffusive Stage

[15] Figure 2 shows a first stage in the thermal relaxation corresponding to isothermal spreading of the buoyant hot region below the surface. This phenomenon of viscous gravity currents has been widely studied [Bercovici and Lin, 1996; Koch and Koch, 1995; Huppert, 1982; Koch and Manga, 1996]. The evolution of the shape is comparable to these works even though they were either designed to study mantle plumes fed by a deeper conduit [Bercovici and Lin, 1996] or compositional plumes [Koch and Manga, 1996]. Moreover, the analytical results and scaling laws given by Koch and Koch [1995] have been mostly derived in a regime where $R^* \gg a^*$ which is not really the case in our calculations.

[16] During the advective stage (see Figure 1 (middle)), the aspect ratio of the drop is increasing while the temperature and the volume of the thermal anomaly remain nearly constant, i.e.

$$ \frac{a^*}{2} R^{*2} \sim 1. \qquad (8) $$

The second stage of thermal relaxation is dominated by diffusion. After the hot drop stops flattening, heat is diffused in all directions and more efficiently through the top isothermal cold surface. As a consequence, $R^*$ and $a^*$ increase with time as seen in Figure 2. The evolution of $a^*(t^*)$ is rapidly consistent with a purely diffusive model: $a^*(t^*) \sim \sqrt{2\kappa t^*}$. The lateral extent, $R^*(t^*)$, evolves more slowly but reaches a similar diffusive behavior after a long time. The temperature decreases rapidly with the power of $-1.8$ (see Figure 1 (bottom) and Figure 2).

#### 3.2. Time and Length Scales

[17] The transition from the advective to the diffusive stage happens when the diffusion velocity, $\kappa/a$, overcomes the advection velocity which is of order $\alpha\rho g\Delta T a^2 / \eta_0$. This simple balance implies that

$$ \frac{a_{\min}^*}{2} = c_1 Ra_{ic}^{-1/3}, \qquad (9) $$

where $c_1$ is a constant. The volume of the hot anomaly being constant, the radius of the thermal anomaly at the end
<figure><img src="image_12.png" /><figcaption>**Figure 3.** (top) Minimum thickness $a_{min}^*/2$, (middle) radial extent $R^*$ and (bottom) time of equilibrium $t_{adv}^*$ as functions of $Ra_{ic}$ for different viscosity contrasts (black circles for a uniform viscosity, brown squares and black triangles for $\lambda = 10^{-1}$ and $\lambda = 10^{-2}$). The dashed lines correspond to the predictions of equations (9),(10), and (12) (we use $c_1 = 1.7$, $c_2 = 0.2$).</figcaption></figure>

of the advective stage, $R_{adv}^*$ is easily obtained by combining equations (8) and (9):

$$R_{adv}^* = c_1^{-1/2} Ra_{ic}^{1/6}. \qquad (10)$$

[18] The time $t_{adv}$ at the end of the advection stage corresponds to the time needed to advect the bottom of the thermal anomaly from its initial depth $2R_{ic}$ to its final depth $a_{min}$. As the vertical velocity $\partial a/\partial t$ is of order of $-\alpha \rho g \Delta T a^2 / \eta_0$, we get

$$\frac{\partial a}{\partial t} = -c_2 \frac{\alpha \rho g \Delta T a^2}{\eta_0}, \qquad (11)$$

where $c_2$ is a geometrical factor.

[19] Integration of equation (11) from $a(0) = 2R_{ic}$ to $a(t_{adv}) = a_{min}$ using equation (9) implies that the end of the advection phase occurs at

$$t_{adv}^* = \frac{1}{c_2} \left( \frac{1}{2c_1} Ra_{ic}^{1/3} - \frac{1}{2} \right). \qquad (12)$$

[20] These scalings of equations (9), (10), and (12) are confirmed by fitting the results of the numerical experiments shown in Figure 3 with $c_1 \sim 1.7$ and $c_2 \sim 0.2$.

[21] Of course, the transition between advective and diffusive stages only occurs when the initial size of the isobaric core is larger than the minimum thickness given by equation (9). This threshold $R_{ic}^t$ obtained when $a_{min}^* = 2$ corresponds to the threshold Rayleigh number

$$Ra_{ic}^t = c_1^3 = 4.9. \qquad (13)$$

For $Ra_{ic} < Ra_{ic}^t$, $a_{min}^* = 2R_{ic}$ and $t_{adv}^*$ is not defined. Below $Ra_{ic}^t$ the heat is diffused out without advection.

[22] The previous scalings obtained for a uniform viscosity are also valid for large viscosity contrasts. Our simulations depicted in Figure 3 (squares for $\lambda = 10^{-1}$ and triangles for $\lambda = 10^{-2}$) show that large viscosity contrasts enable the drop an easier spreading below the surface in agreement with Koch and Koch [1995]. As the resistance to internal shearing decreases with $\lambda$, horizontal velocity contrasts are more important for low viscosities. As a result, the thickness decreases by about 10%, the radial extent increases by a similar amount and the advection time decreases by a factor $\sim 2$. The temperature dependence of the viscosity affects our results by a minor amount because the readjustment is mostly controlled by the viscosity far from the isobaric core.

[23] The scaling laws with physical dimensions can be easily expressed. Using equation (2) and assuming that the planet density remains uniform so that $g = 4/3\pi G\rho R$, the minimal thickness of the thermal anomaly and the time to reach this thickness are

$$a_{min} = 2b_1 \frac{L^2}{R} \qquad (14)$$

and

$$t_{adv} = b_2 \frac{L^2}{\kappa} \left(\frac{L}{R}\right)^2 \left(1 - b_1 \frac{L^2}{RR_{ic}}\right). \qquad (15)$$

In these expressions, $b_1$ and $b_2$ are dimensionless constants,

$$b_1 = \frac{3}{2} c_1 \left( \frac{f(m)}{2\gamma\pi^2} \right)^{1/3} \sim 1.96, \qquad b_2 = \frac{b_1^2}{2c_1^3 c_2} \sim 1.96, \qquad (16)$$

and the properties of the impacted planet appear through a characteristic length

$$L = \left( \frac{C_p \kappa \eta_0}{\alpha \rho^3 G^2} \right)^{1/6} \sim 212 \text{ km.} \qquad (17)$$

## 4. Discussion and Conclusion

[24] We developed a thermo-mechanical model for the long term relaxation after an impact. In a first stage, the heated region spreads below the surface until diffusive cooling becomes more effective. The transition between the advective and diffusive stages is described by a thickness $a_{min}$ and timescale $t_{adv}$ for which we proposed scalings
<figure><img src="image_13.png" /><figcaption>**Figure 4.** (top) Temperature increase $\Delta T$, (middle) thickness $a_{\min}/2$ and (bottom) advection time $t_{adv}$ as functions of the planetary radius. For too small impactors (right side label), advection does not occur. The advection time is plotted for different impactor radii.</figcaption></figure>

laws. Hence we can predict geometrical and time evolution of the thermal anomaly caused by a meteorological impact as functions of rheological parameters of the impacted planetesimal and impactor. All our results are summarized in Figure 4.

[25] The temperature increase (Figure 4, top) and the thickness of the thermal anomaly after advection (Figure 4, middle) do not depend on the initial size of the impactor but only on the properties of the impacted body (and therefore its radius assuming known its other properties, see equations (2) and (14)). As the volume of the isobaric core is proportional to that of the impactor the minimum thickness of the thermal anomaly corresponds also to the minimum radius of the impactor that can trigger advection (Figure 4, middle). For a Mars size planet ($R = 3400$ km), impacts increase the temperature by 390 K (in relative agreement with Senshu et al. [2002]) and post-impact advection only occurs for impactors with radius larger than 18 km. After advection, the final thickness of the thermal anomaly is 52 km. For smaller impactors, only heat diffusion occurs.

[26] The duration of advection depends on the impactor size (Figure 4, bottom). As shown in equation (15), the time of advection is lower than a threshold value obtained for an infinitely large impactor (which of course would disrupt the planet). For a Mars size planet impacted by bodies with $1/10$ to $1/100$ smaller radii, advection ends up after, 10 Myr,

5 Myr, respectively. After this advective stage, heat is slowly removed by diffusion in ~20 Myr.

[27] These timescales are of the same order as those for accretion and differentiation [Yin et al., 2002]. Hence, until impact melting is efficient, heat brought by impacts is stored within the mantle even taking into account of the deformation of the heated region. The scalings proposed here could be used to compute more accurate one dimensional thermal evolution models of growing planets.

[28] **Acknowledgments.** The authors thank J. Schmalzl, H. Samuel and an anonymous reviewer for their useful comments about the manuscript. They also thank J. Melosh and M. Monnereau for helpful discussions. This project was funded by the A.N.R. ETHER.

## References

Alvarez, L. W., W. Alvarez, F. Asaro, and H. V. Michel (1980), Extraterrestrial cause for the Cretaceous Tertiary extinction, *Science*, **208**, 1095–1108.

Bercovici, D., and J. Lin (1996), A gravity current model of cooling mantle plume heads with temperature-dependent buoyancy and viscosity, *J. Geophys. Res.*, **101**, 3291–3310.

Croft, S. K. (1982), A first-order estimate of shock heating and vaporization in oceanic impacts, in *Geological Implications of Impacts of Large Asteroids and Comets on Earth*, edited by T. L. Silver and P. H. Schultz, *Spec. Pap. Geol. Soc. Am.*, **90**, 143–152.

Dubuffet, F., M. Rabinowicz, and M. Monnereau (2000), Multiple scales in mantle convection, *Earth Planet. Sci. Lett.*, **178**, 351–366.

Hartmann, W. K., and D. R. Davis (1975), Satellite-sized planetesimals and lunar origin, *Icarus*, **24**, 504–514.

Huppert, H. E. (1982), The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface, *J. Fluid Mech.*, **121**, 43–58.

Kaula, W. M. (1979), Thermal evolution of Earth and Moon growing by planetesimal impacts, *J. Geophys. Res.*, **84**, 999–1008.

Koch, D. M., and D. L. Koch (1995), Numerical and theoretical solutions for a drop spreading below a free fluid surface, *J. Fluid Mech.*, **287**, 251–278.

Koch, D. M., and M. Manga (1996), Neutrally buoyant diapirs: A model for Venus coronae, *Geophys. Res. Lett.*, **23**, 225–228.

Kokubo, E., and S. Ida (1996), On runaway growth of planetesimals, *Icarus*, **123**, 180–191.

Maxwell, D. E. (1977), Simple Z model for cratering, ejection, and the overturned flap, in *Impact and Explosion Cratering: Planetary and Terrestrial Implications*, edited by D. J. Roddy, R. O. Pepin, and R. B. Merrill, pp. 1003–1008, Elsevier, New York.

O'Keefe, J. D., and T. J. Ahrens (1977), Impact-induced energy partitioning, melting, and vaporization on terrestrial planets, *Lunar Planet. Sci. Conf.*, **8th**, 3357–3374.

Pierazzo, E., A. M. Vickery, and H. J. Melosh (1997), A reevaluation of impact melt production, *Icarus*, **408**, 423.

Reese, C. C., V. S. Solomatov, and J. R. Baumgardner (2002), Survival of impact-induced thermal anomalies in the Martian mantle, *J. Geophys. Res.*, **107(E10)**, 5082, doi:10.1029/2000JE001474.

Schubert, G., D. L. Turcotte, and P. Olson (2001), *Mantle Convection in the Earth and Planets*, Cambridge Univ. Press, New York.

Senshu, H., K. Kuramoto, and T. Matsui (2002), Thermal evolution of a growing Mars, *J. Geophys. Res.*, **107(E12)**, 5118, doi:10.1029/2001JE001819.

Solomatov, V. S. (2000), *Fluid Dynamics of a Terrestrial Magma Ocean, Origin of the Earth and Moon*, edited by R. M. Canup et al., pp. 323–338, Univ. of Ariz. Press, Tucson.

Tonks, W. B., and H. J. Melosh (1993), Magma ocean formation due to giant impacts, *J. Geophys. Res.*, **98**, 5319–5333.

Turtle, E. P., E. Pierazzo, and D. P. O'Brien (2003), Numerical modeling of impact heating and cooling of the Vredefort impact structure, *Meteorit. Planet. Sci.*, **38**, 293–303.

Yin, Q., S. B. Jacobsen, K. Yamashita, J. Blichert-Toft, P. Télouk, and F. Albarède (2002), A short timescale for terrestrial planet formation from Hf-W chronometry of meteorites, *Nature*, **418**, 949–952.

<footnote>N. Coltice, F. Dubuffet, J. Monteux, and Y. Ricard, Laboratoire de Sciences de la Terre, UMR 5570, CNRS, Université de Lyon, Université Claude Bernard Lyon 1, 2 rue R. Dubois, F-69622, Villeurbanne cedex, France. (nicolas.coltice@univ-lyon1.fr; fabien.dubuffet@univ-lyon1.fr; julien.monteux@univ-lyon1.fr; ricard@ens-lyon.fr)</footnote>
<figure><img src="image_14.png" /></figure>

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

# Scaling laws of impact induced shock pressure and particle velocity in planetary mantle

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

CrossMark

J. Monteux<sup>a,*</sup>, J. Arkani-Hamed<sup>b,c</sup>

## ARTICLE INFO

**Article history:**

Received 13 March 2015

Revised 24 September 2015

Accepted 24 September 2015

Available online 3 October 2015

**Keywords:**

Accretion

Impact processes

Thermal histories

Terrestrial planets

Interiors

## ABSTRACT

While major impacting bodies during accretion of a Mars type planet have very low velocities (<10 km/s), the characteristics of the shockwave propagation and, hence, the derived scaling laws are poorly known for these low velocity impacts. Here, we use iSALE-2D hydrocode simulations to calculate shock pressure and particle velocity in a Mars type body for impact velocities ranging from 4 to 10 km/s. Large impactors of 100–400 km in diameter, comparable to those impacted on Mars and created giant impact basins, are examined. To better represent the power law distribution of shock pressure and particle velocity as functions of distance from the impact site at the surface, we propose three distinct regions in the mantle: a near field regime, which extends to 1–3 times the projectile radius into the target, where the peak shock pressure and particle velocity decay very slowly with increasing distance, a mid field region, which extends to ~4.5 times the impactor radius, where the pressure and particle velocity decay exponentially but moderately, and a more distant far field region where the pressure and particle velocity decay strongly with distance. These scaling laws are useful to determine impact heating of a growing protoplanet by numerous accreting bodies.

© 2015 Elsevier Inc. All rights reserved.

## 1. Introduction

Small planets are formed by accreting a huge number of planetesimals, a few km to a few tens of km in size, in the solar nebula (e.g. Wetherill and Stewart, 1989; Matsui, 1993; Chambers and Wetherill, 1998; Kokubo and Ida, 1995, 1996, 1998, 2000; Wetherill and Inaba, 2000; Rafikov, 2003; Chambers, 2004; Raymond et al., 2006). An accreting body may generate shock wave if the impact-induced pressure in the target exceeds the elastic Hugoniot pressure, ~3 GPa, implying that collision of a planetesimal with a growing planetary embryo can generate shock waves when the embryo's radius exceeds 150 km, assuming that impact occurs at the escape velocity of the embryo and taking the mean density of the embryo and projectile to be 3000 kg/m³. Hundreds of thousands of collisions must have occurred during the formation of small planets such as Mercury and Mars when they were orbiting the Sun inside a dense population of planetesimal. Such was also the case during the formation of embryos that later were accreted to produce Venus and Earth. Terrestrial planets have also experienced large high velocity impacts after their formation. Over 20 giant impact basins on Mars with diameters larger than

1000 km (Frey, 2008), the Caloris basin on Mercury with a 1550 km diameter, and the South Pole Aitken basin on Moon with a 2400 km diameter are likely created during catastrophic bombardment period at around 4 Ga. The overlapping Rheasilvia and Veneneia basins on 4-Vesta are probably created by projectiles with an impact velocity of about 5 km/s within the last 1–2 Gyr (Keil et al., 1997; Schenk et al., 2012).

The shock wave produced by an impact when the embryo is undifferentiated and completely solid propagates as a spherical wave centered at the impact site until it reaches the surface of the embryo in the opposite side. Each impact increases the temperature of the embryo within a region near the impact site. Because impacts during accretion occur from different directions, the mean temperature in the upper parts of the embryo increases almost globally. On the other hand, the shock wave produced by a large impact during the heavy bombardment period must have increased the temperature in the mantle and the core of the planets directly beneath the impact site, enhancing mantle convection (e.g. Watters et al., 2009; Roberts and Arkani-Hamed, 2012, 2014), modifying the CMB heat flux which could in turn favor a hemispheric dynamo on Mars (Monteux et al., 2015), or crippling the core dynamo (e.g. Arkani-Hamed and Olson, 2010a).

<footnote><sup>a</sup> Laboratoire Magmas et Volcanes, Clermont-Ferrand, France</footnote>

<footnote><sup>b</sup> Department of Physics, University of Toronto, Toronto, Canada</footnote>

<footnote><sup>c</sup> Department of Earth and Planetary Sciences, McGill University, Montreal, Canada</footnote>

<footnote>* Corresponding author.</footnote>
The impact-induced shock pressure inside a planet has been investigated by numerically solving the shock dynamic equations using hydrocode simulations (e.g. Pierazzo et al., 1997; Wünnemann and Ivanov, 2003; Wünnemann et al., 2006; Barr and Citron, 2011; Kraus et al., 2011; Ivanov et al., 2010; Bierhaus et al., 2012) or finite difference techniques (e.g. Ahrens and O'Keefe, 1987; Mitani, 2003). However, these numerical solutions demand considerable computer capacity and time and are not practical for investigating the huge number of impacts that occur during the growth of a planet. Hence, the scaling laws derived from field experiments (e.g. Perret and Bass, 1975; Melosh, 1989) or especially from hydrocode simulations (Pierazzo et al., 1997) are of great interest when considering the full accretionary history of a planetary objects (e.g. Senshu et al., 2002; Monteux et al., 2014) or when measuring the influence of a single large impact on the long-term thermal evolution of deep planetary interiors (e.g. Monteux et al., 2007, 2009, 2013; Ricard et al., 2009; Roberts et al., 2009; Arkani-Hamed and Olson, 2010a; Arkani-Hamed and Ghods, 2011). Although the scaling laws provide approximate estimates of the shock pressure distribution, their simplicity and the small differences between their results and those obtained by the hydrocode simulations of the shock dynamic equations (that are likely within the numerical errors that could have been introduced due to the uncertainty of the physical parameters used in the hydrocode models) make them a powerful tool that can be combined with other geophysical approaches such as dynamo models (e.g. Monteux et al., 2015) or convection models (e.g. Watters et al., 2009; Roberts and Arkani-Hamed, 2012, 2014).

During the decompression of shocked material much of the internal energy of the shock state is converted into heat leading to a temperature increase below the impact site. The present study focuses on deriving scaling laws of shock pressure and particle velocity distributions in silicate mantle of a planet on the basis of several hydrocode simulations. The scaling laws of Pierazzo et al. (1997) were derived using impact velocities of 10–80 km/s, hence may not be viable at low impact velocities. For example, at an impact velocity of 5 km/s, comparable to the escape velocity of Mars, the shock pressure scaling law provides an unrealistic shock pressure that increases with depth. Here we model shock pressure and particle velocity distributions in the mantle using hydrocode simulations for impact velocities of 4–10 km/s and projectile diameters ranging from 100 to 400 km, as an attempt to extend Pierazzo et al.'s (1997) scaling laws to low impact velocities and reasonable impactor radii occurring during the formation of terrestrial planets. Hence, on the basis of our scaling laws it is possible to estimate the temperature increase as a function of depth below the impact site for impact velocities compatible with the accretionary conditions of terrestrial protoplanets. These scaling laws can easily be implemented in a multi-impact approach (e.g. Senshu et al., 2002; Monteux et al., 2014) to monitor the temperature evolution inside a growing protoplanet whereas it is not yet possible to adopt hydrocode simulations for that purpose.

The hydrocode models we have calculated are described in the first section, while the second section presents the scaling laws derived from the hydrocode models. The concluding remarks are relegated to the third section.

## 2. Hydrocode models of shock pressure distribution

The huge number of impacts during accretion makes it impractical to consider oblique impacts. Not only it requires formidable computer time, but more importantly because of the lack of information about the impact direction, i.e. the impact angle relative to vertical and azimuth relative to north. Therefore, we consider only head-on collisions (vertical impact) to model the thermo-mechanical

evolution during an impact between a differentiated Mars size body and a large impactor. We use the iSALE-2D axisymmetric hydrocode, which is a multi-rheology, multi-material hydrocode, specifically developed to model impact crater formation on a planetary scale (Collins et al., 2004; Davison et al., 2010). To limit computation time, we use a 2 km spatial resolution (i.e. more than 25 cells per projectile radius, cppr) and a maximum time step of 0.05 s which is sufficient to describe the shockwave propagation through the entire mantle. The minimum post impact monitoring time is set to the time needed by the shockwave to reach the core–mantle boundary (≈5 min for the impact velocities studied here).

We investigate the shock pressure and particle velocity distributions inside a Mars size model planet for impact velocities $V_{imp}$ of 4–10 km/s and impactors of 100–400 km in diameter. Such impactors are capable of creating basins of 1000–2500 km in diameter according to Schmidt and Housen (1987) and Holsapple (1993) scaling relationships between the impactor diameter and the resulting basin diameter. These basins are comparable with the giant impact basins of Mars created during the heavy bombardment period at around 4 Ga (Frey, 2008).

In our models, the impactor was simplified to a spherical body of radius $R_{imp}$ with uniform composition while the target was simplified to a two layers spherical body of radius $R$ and an iron core radius of $R_{core}$. The silicate mantle has a thickness of $\delta_m$ (see Table 1). We adopt physical properties of silicates (dunite or peridotite) for both the mantle and the impactor to monitor the shock pressure and the particle velocity in a Mars type body. We approximate the thermodynamic response of both the iron and silicate material using the ANEOS equation of state (Thompson and Lauson, 1972; Benz et al., 1989). To make our models as simple as possible we do not consider here the effects of porosity, thermal softening or low density weakening. However, as a first step towards more realistic models, we investigate the influence of acoustic fluidization and damage. All these effects can be accounted for in iSALE-2D and we will consider each effect in a separate study in near future.

<table><caption>**Table 1**
Typical parameter values for numerical hydrocode models.</caption><tbody><tr><td>Target radius</td><td>R</td><td>3400 km</td></tr><tr><td>Target core radius</td><td>R<sub>core</sub></td><td>1700 km</td></tr><tr><td>Silicate mantle thickness</td><td>δ<sub>m</sub></td><td>1700 km</td></tr><tr><td>Impactor radius</td><td>R<sub>imp</sub></td><td>50–200 km</td></tr><tr><td>Impact velocity</td><td>V<sub>imp</sub></td><td>4–10 km/s</td></tr><tr><td colspan="3"><strong>Mantle properties (silicates)</strong></td></tr><tr><td>Initial density</td><td>ρ<sub>m</sub></td><td>3314 kg/m<sup>3</sup></td></tr><tr><td>Equation of state type</td><td></td><td>ANEOS</td></tr><tr><td>Poisson</td><td></td><td>0.25</td></tr><tr><td>Strength Model<br>(iSALE parameters)</td><td></td><td>Rock<br>(Y<sub>i0</sub> = 10 MPa, μ<sub>i</sub> = 1.2,<br>Y<sub>im</sub> = 3.5 GPa)</td></tr><tr><td>Acoustic Fluidization Model<br>(iSALE parameters)</td><td></td><td>Block<br>(t<sub>off</sub> = 16 s, c<sub>vib</sub> = 0.1 m/s,<br>vib<sub>max</sub> = 200 m/s)</td></tr><tr><td>Damage Model<br>(iSALE parameters)</td><td></td><td>Ivanov<br>(ε<sub>fb</sub> = 10<sup>-4</sup>, B = 10<sup>-11</sup>,<br>ρ<sub>c</sub> = 3 × 10<sup>8</sup> Pa)</td></tr><tr><td>Thermal softening and porosity<br>models</td><td></td><td>None</td></tr><tr><td colspan="3"><strong>Core properties (iron)</strong></td></tr><tr><td>Initial density</td><td>ρ<sub>c</sub></td><td>7840 kg/m<sup>3</sup></td></tr><tr><td>Equation of state type</td><td></td><td>ANEOS</td></tr><tr><td>Poisson</td><td></td><td>0.3</td></tr><tr><td>Strength Model<br>(iSALE parameters)</td><td></td><td>Von Mises<br>(Y<sub>0</sub> = 100 MPa)</td></tr><tr><td>Acoustic Fluidization Model<br>(iSALE parameters)</td><td></td><td>Block<br>(t<sub>off</sub> = 16 s, c<sub>vib</sub> = 0.1 m/s,<br>vib<sub>max</sub> = 200 m/s)</td></tr><tr><td>Damage, thermal softening and<br>porosity models</td><td></td><td>None</td></tr></tbody></table>
<figure><img src="image_17.png" /><figcaption>**Fig. 1.** Pre-impact temperature (left) and lithostatic pressure (right) as a function of depth. The dashed lines illustrate the core–mantle boundary.</figcaption></figure>

The early temperature profile of a Mars size body is difficult to constrain because it depends on its accretionary history, on the amount of radiogenic heating and on the mechanisms that led to its core formation (e.g. Senshu et al., 2002). The uncertainties on the relative importance of these processes as well as the diversity of the processes involved in the core formation lead to a wide range of plausible early thermal states after the full differentiation of Mars. Since we do not consider here the thermal softening during the impact, we assumed the same radially dependent preimpact martian temperature field in all our simulations. Fig. 1 shows the pre-impact temperature profile as a function of depth. As the pre-impact pressure is governed by the material repartition and as we consider a differentiated Mars, the pre-impact pressure is more straightforward. Fig. 1 also shows the pre-impact hydrostatic pressure used in our models as a function of depth considering a 1700 km thick silicate mantle surrounding a 1700 km radius iron core. We emphasize that the peak pressure shown in our study does not include this hydrostatic pressure. However, the hydrostatic pressure is taken into consideration in calculating the hydrocode models. The peak pressure presents the shock induced pressure increase and is expected to depend on the physical properties of the target, but not on the size of the target as long as the size is large enough to allow shock wave propagates freely without interference with reflected waves. A scaling law should reflect the properties of the shock wave propagation in a uniform media. Fig. 2 shows the typical time evolution of the compositional and pressure fields after a 100 km diameter impact with $V_{\text{imp}} = 10$ km/s. Immediately after the impact, the shockwave propagates downward from the impact site. The shock front reaches the core–mantle boundary in less than 5 min while the transient crater is still opening at the surface. It is worth mentioning that the main goal is to derive a scaling law which is useful for numerous impacts during accretion where no information is available about the impact direction, i.e. the impact angle relative to vertical and azimuth relative to north. Moreover, the pressure reduction near the surface due to interference of the direct and reflected waves can easily be accommodated following the procedure by Melosh (1989), which is adapted to spherical surface by Arkani-Hamed (2005), when applying the scaling law to a particular accretion scenario.

In Fig. 3, we monitor the peak pressure as a function of the distance from the impact site *d* normalized by the impactor radius $R_{\text{imp}}$ along the symmetry axis for the case illustrated in Fig. 2. In our iSALE models, the impact-induced pressure fields (as well as temperature and velocity fields) are extracted from a cell-centered Eulerian grid data. To validate our models, we have tested different spatial resolutions expressed here in terms of cells per projectile radius (cppr). In Fig. 3a, we represent the peak pressure

<figure><img src="image_18.png" /><figcaption>**Fig. 2.** A close up view of the material repartition (left column) and total pressure (right column) as functions of time (from top to bottom) in the model planet (for $V_{\text{imp}} = 10$ km/s and $D_{\text{imp}} = 100$ km). In this model, the grid resolution is 2 km in all directions. The silicate mantle and the impactor are made of dunite.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 3. Peak pressure decrease as a function of depth normalized by the radius of the impactor for the impact velocity of 10 km/s. The silicate mantle as well as the impactor are made of dunite. (a) Influence of the spatial resolution. Here we only consider the case with R<sub>imp</sub> = 50 km. The results from our hydrocode models are shown by colored curves with a spatial resolution ranging from 5 to 50 cppr and (b) comparison of our results with R<sub>imp</sub> = 50 km (red curve, 25 cppr) and R<sub>imp</sub> = 10 km (green curve, 10 cppr) with the results from a similar model of Pierazzo et al. (1997) (black squares). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

decrease as a function of depth for cppr values ranging from 5 to 50, showing convergence for cppr values larger than 25. As illustrated in Fig. 3a, the difference between the results with 25 and 50 cppr is less than 10%. This resolution study is in agreement with Pierazzo et al. (2008) who have shown that the iSALE models converge for resolutions of 20 cppr or higher, although a resolution of 10 cppr still provides reasonable results (a resolution of 20 cppr appears to underestimate peak shock pressures by at most 10%). The resolution is 25 cppr or higher in our models.

In Pierazzo et al. (1997), the impactor radius ranged between 0.2 and 10 km. In Fig. 3b, we compare our results obtained with R<sub>imp</sub> = 10 km, R<sub>imp</sub> = 50 km and the results obtained by Pierazzo et al. (1997). Fig. 3b shows that even with a radius of 50 km, both our results and the results from Pierazzo et al. (1997) are in good agreement, confirming that the impactor size has minor effects on the peak pressure evolution with depth. The small differences between our results with different impactor radii (discussed further in more details) are plausibly the direct consequence of using increasing cppr values with increasing impactor radii. Consequently, we will use the normalized distance in all of our figures, as equations of motion should be invariant under rescaling of distance (Melosh, 1989; Pierazzo et al., 1997). However, to monitor the peak pressure evolution with R<sub>imp</sub> = 10 km and to maintain a reasonable computational time, we have used only 10 cppr. This figure confirms that below 10 cppr, the spatial accuracy is insufficient and our results diverge from the results obtained by Pierazzo et al. (1997).

### 3. Shock pressure and particle velocity scaling laws at low impact velocities

A given hydrocode simulation may take on the order of 48 h to determine a 2D shock pressure and particle velocity distributions in the mantle of our model planet. The impact velocity is about 4 km/s for a protoplanet with a radius of 2860 km and mean density of 3500 kg/m³, assuming that impacts occur at the escape velocity of the protoplanet. Mars is more likely a runaway planetary embryo formed by accreting small planetesimals and medium size neighboring planetary embryos. This indicates that the accreting bodies had velocities higher than 4 km/s when Mars

was growing from 2860 km radius to its present radius of about 3400 km. Taking the mean radius of the impacting bodies to be 100 km, which is larger than that of a typical planetesimal, more than 15,000 bodies must have accreted at impact velocities higher than 4 km/s. Calculating the impact induced shock pressure and particle velocity inside the growing Mars would be formidable if hydrocode simulation is adopted for each impact. Because the impact-induced shock pressure P and particle velocity V<sub>p</sub> inside a planetary mantle decrease monotonically with distance from the impact side, simple exponential functions have been proposed to estimate peak pressure and particle velocity in the mantle of an impacted body. Solving the shock dynamic equations by a finite difference technique for silicate target and projectile, Ahrens and O'Keefe (1987) showed that pressure distribution in the target displays three regimes: an impedance match regime, Regime I, which extends to 1–3 times the projectile radius into the target where the peak shock pressure is determined by the planar impedance match pressure (McQueen et al., 1970); a shock pressure decay regime, Regime II, where the pressure decays exponentially as

$$ 
\begin{gathered}
P = P_o(d/R_{imp})^n, \quad \text{for } d > 1-3 \text{ times } R_{imp}, \text{ and} \\ n = -1.25 - 0.625 \log(V_{imp}) \quad (1)
\end{gathered}
$$

and the elastic regime, Regime III, where the shock pressure is reduced below the strength of target, the Hugoniot elastic limit, and the shock wave is reduced to an elastic wave with pressure decaying as d<sup>-3</sup>. In Eq. (1) d is the distance from the impact site at the surface, R<sub>imp</sub> is the projectile radius, and V<sub>imp</sub> is the impact velocity in km/s. The peak pressure measurements in the nuclear explosions (Perret and Bass, 1975) led Melosh (1989) to propose a scaling law for the particle velocity distribution assuming conservation of momentum of the material behind the shock front which is taken to be a shell of constant thickness. Using several different target materials, and adopting hydrocode simulations Pierazzo et al. (1997) showed that the shock pressure P and particle velocity V<sub>p</sub> decrease slowly in the near field zone, but rapidly in the deeper region,

$$ P = P_{ic}(d_{ic}/d)^n, \quad n = -1.84 + 2.61 \log(V_{imp}), \quad d > d_{ic} \quad (2a) $$

$$ V_p = V_{pic}(d_{ic}/d)^m \quad m = -0.31 + 1.15 \log(V_{imp}), \quad d > d_{ic} \quad (2b) $$
<figure><img src="image_20.png" /><figcaption>Fig. 4. Shock pressure versus normalized distance from the impact site at the surface produced by a 100 km diameter impactor with an impact velocity of 10 km/s. The dashed curve represents the hydrocode model, while the straight lines are fitted to three different parts of the hydrocode model.</figcaption></figure>

The authors coined an isobaric zone of shock pressure $P_{ic}$ and particle velocity $V_{pic}$ for the near field of radius $d_{ic}$, about 1.5 $R_{imp}$. Eqs. (2a) and (2b) were derived by averaging results from many different materials and impactor sizes. The impact velocities adopted were 10–80 km/s and the projectile diameter ranged from 0.4 to 20 km.

In a log-log plot the peak shock pressure and the corresponding particle velocity are linear functions of distance from the impact site,

$$ \text{Log } P = a + n \text{Log}(d/R_{\text{imp}}) \quad (3a) $$

$$ \text{Log } V_p = c + m \text{Log}(d/R_{\text{imp}}) \quad (3b) $$

where a is the logarithm of pressure and c is the logarithm of particle velocity both at $R_{imp}$ and n and m are decay factors. All parameters are impact velocity dependent:

$$ a = \alpha + \beta \text{Log}(V_{\text{imp}}) \quad (4a) $$

$$ c = \gamma + \Omega \text{Log}(V_{\text{imp}}) \quad (4b) $$

and

$$ n = \lambda + \delta \text{Log}(V_{\text{imp}}) \quad (5a) $$

$$ m = \eta + \zeta \text{Log}(V_{\text{imp}}) \quad (5b) $$

Fig. 4 shows the peak shock pressure inside the mantle of the model planet we obtained by hydrocode simulation and using a projectile of 100 km diameter at 10 km/s impact velocity. There is actually no isobaric region, rather the peak pressure decays slowly in the near field zone but much rapidly in the deeper parts, similar to the results by Ahrens and O’Keefe (1987). Although the regression lines fitted to the near field and far field are good representatives of the shock pressure distribution, they intersect at a much higher pressure than that of the hydrocode model and overestimate the pressure by as much as 30% in a large region located between the near field and the far field. Therefore, to better approximate the pressure distribution we fit the pressure curve by three lines, representing the near field, mid field and far field regions, which render a much better fit as seen in Fig. 4.

Fig. 5a shows the peak shock pressure versus distance from the impact site for impact velocities of 4–10 km/s and an impactor of 100 km in diameter. Because of different phenomenon such as excavation, melting, vaporization and intermixing between the target and projectile material, the shock front in near field is more difficult to characterize by our numerical models even for high cppr. Also, as the shock wave is not yet detached from the impactor, it cannot be treated as a single shock wave. Hence, the near field–mid field boundary and the scaling laws for the near field are less accurate than for the two other fields especially for large impact velocities ($V_{imp} > 7$ km/s). In Fig. 5, the larger dots show the intersections of the linear regression lines. For example a dot that separates near field from mid field is the intersection of the regression lines fitted to the near field and midfield. The regression lines are determined from fitting to the hydrocode data. Visually, we first divide the hydrocode data of a given model into three separate sections with almost linear trends, and then fit the regression lines to those three trends. Fig. 5a shows the hydrocode data, small dots, and the regression lines, demonstrating the tight fitting of the lines to the data. The large dots in the figure show the intersection of the regression lines of the adjacent regions. For example the dot that shows the near-field mid-field boundary is the intersection of the regression lines fitted to the near field and mid field hydrocode data. Note that a big dot does not necessarily coincide with the exact hydrocode result, a small dot. As the near field–mid field boundary is relatively less resolved, Fig. 5a shows a scatter of the dots separating those two fields and a slight slope change from $V_{imp} > 7$ km/s. In the average the near field–mid field boundary is at ~2.3 $R_{imp}$ (~115 km) from the impact site. Fig. 5a shows that the depth to the mid field–far field boundary is almost independent of the pressure. It is at about 4.5 $R_{imp}$ (~225 km) from the impact site.

<figure><img src="image_21.png" /><figcaption>Fig. 5. (a) Shock pressure versus normalized distance from the impact site at the surface for an impactor of 100 km diameter and impact velocities ranging from 4 to 10 km/s. The numbers on the curves are the impact velocities. The hydrocode results are presented by dots, while the regression lines to the near field, mid field and far field regions are straight lines. The larger dots show the intersections of the linear regression lines. For example a dot that separates near field from mid field is the intersection of the regression lines fitted to the near field and midfield data and (b) shows the corresponding particle velocity.</figcaption></figure>
We propose three scaling laws for shock pressure, and three for particle velocity:

$$ \text{Log } P_i = a_i + n_i \text{Log}(d_i/R_{\text{imp}}), \quad i = 1, 2, \text{ and } 3 \qquad (6a) $$

$$ \text{Log } V_{\text{pi}} = c_i + m_i \text{Log}(d_i/R_{\text{imp}}), \quad i = 1, 2, \text{ and } 3 \qquad (6b) $$

For the near field $d_i < 2.3 R_{\text{imp}}$, mid field $2.3 R_{\text{imp}} < d_i < 4.5 R_{\text{imp}}$, and far field $d_i > 4.5 R_{\text{imp}}$. Table 2 lists the values of the parameters in the above equations as well as the misfits obtained from our regressions (smaller than 0.001 in all the regressions calculated here). Fig. 5a shows the three regression lines fitted to the near field, mid field, and far field of each model for shock pressure, and Fig. 5b displays those for particle velocity. The misfits from Table 2 are based on the fixed end points, shown as dots in Fig. 5.

The close agreement between Pierazzo et al. (1997) model for impact velocity of 10 km/s derived by averaging results of impactors with diameters 0.4–20 km, and our result for the same impact velocity but an impactor of 100 km in diameter indicates

that the shock pressure distribution is less sensitive to impactor size in a log-log plot of pressure versus distance normalized to the impactor radius. To further investigate this point, we calculate models with impactor sizes of 100–400 km in diameter. Fig. 6a shows the hydrocode results for impact velocity of 10 km/s using different impactor size. The result obtained for $R_{\text{imp}} = 10$ km is not included in the figure because of its low cppr value (see Fig. 3b). The curves have almost the same slope in the far field, and small deviations in the mid field and near field.

Pierazzo et al. (1997) used several different rock types for the solid mantle and concluded that their scaling laws are less sensitive to the rock types. Bearing in mind that dunite and peridotite are probably the most representative rocks for solid mantle, we run a hydrocode model using an impactor of 100 km in diameter, impact velocity of 10 km/s, and peridotite as a representative mantle and impactor rock. Fig. 6b compares the results for the dunite and peridotite models. They are indeed very similar, especially in the far field region, where the exponential factor n

Table 2

Parameters of the peak shock pressure distribution and the corresponding particle velocity in the mantle of the Mars size model planet. The pressure is expressed as:

$$ \text{Log}(P) = a + n \text{Log}(d/R_{\text{imp}}) $$

and the particle velocity as:

$$ \text{Log}(V_p) = c + m \text{Log}(d/R_{\text{imp}}) $$

where the pressure P is in GPa, the particle velocity $V_p$ is in km/s, d is the distance from the impact site at the surface, and $R_{\text{imp}}$ is the impactor radius. a and c are the logarithm of pressure and particle velocity at the distance $R_{\text{imp}}$ from the impact site, and n and m are the decay exponents of pressure and particle velocity with distance from the impact site. a, c, n and m are impact velocity dependent:

$$ a = \alpha + \beta \text{Log}(V_{\text{imp}}) $$

$$ n = \lambda + \delta \text{Log}(V_{\text{imp}}) $$

$$ c = \gamma + \Omega \text{Log}(V_{\text{imp}}) $$

$$ m = \eta + \zeta \text{Log}(V_{\text{imp}}) $$

A misfit value is obtained by calculating the standard deviation of a line fitted to the hydrocode data within a given region: $\epsilon = 1/N\sqrt{\sum_{i=1}^{N}(Y_{\text{data}} - Y_{\text{regression}})^2}$, where N is the total number of points, $Y_{\text{data}}$ is the hydrocode result and $Y_{\text{regression}}$ denotes the value obtained by the linear regression. The zero misfit implies that the regression line is fitted to only 2 points, hence an exact fitting.

<table><thead><tr><th>V<sub>imp</sub> (km/s)</th><th>a</th><th>n</th><th>misfit</th><th>c</th><th>m</th><th>misfit</th></tr></thead><tbody><tr><td colspan="7"><strong>Near field</strong></td></tr><tr><td>4</td><td>1.1717</td><td>-0.4530</td><td>1.074E-07</td><td>0.1276</td><td>-1.1132</td><td>3.071E-08</td></tr><tr><td>5</td><td>1.3963</td><td>-0.6296</td><td>4.616E-03</td><td>0.2901</td><td>-1.0437</td><td>6.6837E-04</td></tr><tr><td>6</td><td>1.5137</td><td>-0.4713</td><td>8.429E-08</td><td>0.3722</td><td>-0.8573</td><td>0.000</td></tr><tr><td>7</td><td>1.6527</td><td>-0.3237</td><td>5.960E-08</td><td>0.3315</td><td>-0.4286</td><td>1.490E-08</td></tr><tr><td>8</td><td>1.8093</td><td>-0.3302</td><td>5.960E-08</td><td>0.5736</td><td>-0.6567</td><td>0.000</td></tr><tr><td>9</td><td>1.8853</td><td>-0.1228</td><td>8.429E-08</td><td>0.6817</td><td>-0.2156</td><td>4.214E-08</td></tr><tr><td>10</td><td>1.9072</td><td>-0.1364</td><td>1.332E-07</td><td>0.7354</td><td>-0.2622</td><td>8.411E-03</td></tr><tr><td colspan="7">α = 0.040, β = 1.914, λ = -1.214, δ = 1.058<br>γ = -0.795, Ω = 1.502, η = -2.602, ζ = 2.368</td></tr><tr><td colspan="7"><strong>Mid field</strong></td></tr><tr><td>4</td><td>1.3714</td><td>-0.9459</td><td>2.0576E-03</td><td>0.0917</td><td>-1.0211</td><td>3.8861E-04</td></tr><tr><td>5</td><td>1.5978</td><td>-0.9995</td><td>6.5599E-04</td><td>0.2911</td><td>-1.0402</td><td>7.3517E-04</td></tr><tr><td>6</td><td>1.6367</td><td>-0.8038</td><td>4.7266E-03</td><td>0.4143</td><td>-0.9942</td><td>8.0339E-04</td></tr><tr><td>7</td><td>1.8821</td><td>-0.9792</td><td>2.8677E-03</td><td>0.5141</td><td>-0.9563</td><td>1.8048E-03</td></tr><tr><td>8</td><td>1.9735</td><td>-0.8576</td><td>3.5979E-03</td><td>0.6750</td><td>-0.9888</td><td>6.5566E-05</td></tr><tr><td>9</td><td>2.0060</td><td>-0.7072</td><td>6.2120E-03</td><td>0.8611</td><td>-1.1259</td><td>2.0130E-03</td></tr><tr><td>10</td><td>2.0224</td><td>-0.6059</td><td>8.2497E-03</td><td>0.9317</td><td>-1.1190</td><td>1.3516E-03</td></tr><tr><td colspan="7">α = 0.346, β = 1.736, λ = -1.469, δ = 0.768<br>γ = -1.206, Ω = 2.114, η = -0.864, ζ = -0.208</td></tr><tr><td colspan="7"><strong>Far field</strong></td></tr><tr><td>4</td><td>1.5136</td><td>-1.1453</td><td>6.4139E-04</td><td>0.2397</td><td>-1.2158</td><td>1.1549E-03</td></tr><tr><td>5</td><td>1.7356</td><td>-1.1640</td><td>5.0752E-04</td><td>0.4635</td><td>-1.2389</td><td>1.0468E-03</td></tr><tr><td>6</td><td>1.9107</td><td>-1.1864</td><td>6.9751E-04</td><td>0.6248</td><td>-1.2531</td><td>1.2862E-03</td></tr><tr><td>7</td><td>2.0602</td><td>-1.2182</td><td>6.0862E-04</td><td>0.7628</td><td>-1.2783</td><td>1.3663E-03</td></tr><tr><td>8</td><td>2.2186</td><td>-1.2816</td><td>3.3932E-04</td><td>0.8936</td><td>-1.3220</td><td>1.0471E-03</td></tr><tr><td>9</td><td>2.4057</td><td>-1.3818</td><td>1.1730E-03</td><td>1.0620</td><td>-1.4091</td><td>8.7868E-04</td></tr><tr><td>10</td><td>2.5440</td><td>-1.4492</td><td>1.5957E-03</td><td>1.1887</td><td>-1.4687</td><td>1.3772E-03</td></tr><tr><td colspan="7">α = -0.056, β = 2.558, λ = -0.647, δ = -0.744<br>γ = -1.177, Ω = 2.333, η = -0.818, ζ = -0.600</td></tr></tbody></table>
<figure><img src="image_22.png" /><figcaption>**Fig. 6.** (a) Shows the hydrocode results for impactors of 50–400 km diameter and impact velocity of 10 km/s, (b) compares the hydrocode results using dunite and peridotite as mantle rock types, for an impactor of 100 km diameter and impact velocity of 10 km/s and (c) illustrates the shock pressure as a function of $d/R_{\text{imp}}$ with (red curve) and without (black curve) acoustic fluidization. The green curve represents the results considering an Ivanov damage model. (For comparison, the black squares represent the results from Pierazzo et al. (1997), which has not considered acoustic fluidization). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

in Eq. (3a) is −1.449 for dunite and −1.440 for peridotite. The major differences between the two models arise from the near field zone, hence propagates down to the deeper regions.

Among the other parameters and phenomenon that may influ-
ence the shockwave propagation (porosity, thermal softening...),
accounting for the acoustic fluidization is required to accurately
simulate the formation of a complex crater (e.g. Melosh, 1979;
Bray et al., 2008; Potter et al., 2012). Indeed, acoustic fluidization
is invoked to explain the collapse of complex craters by modifying
the frictional strength of the damaged target. Fig. 6c compares the
results obtained for R<sub>imp</sub> = 50 km and V<sub>imp</sub> = 10 km/s considering
acoustic fluidization and an Ivanov damage model which pre-
scribes damage as a function of plastic strain. This figure shows
that for the near and mid fields, the results are similar. Pierazzo
et al. (1997) did not include acoustic fluidization or damage in their
models. Hence, our results without acoustic fluidization or damage
model and the results from Pierazzo are in good agreement
(Fig. 6c). As soon as the far field is reached, acoustic fluidization
and damage tend to reduce the intensity of the shock pressure. This
indicates that building more sophisticated models will be neces-
sary in the near future. As the impact heating is mainly localized
in the near- and mid-fields, including a damage model or acoustic
fluidization should only weakly affect the thermal evolution of a
growing protoplanet. However, it is worth mentioning that we
are not concerned with the shape of the crater produced by a large
impact, rather the main goal of our study is to extend the scaling
laws of Pierazzo et al. (1997) to lower impact velocities which
are more compatible with accretionary conditions.

Fig. 7 shows the impact velocity dependence of *a*, *n*, *c*, and *m*. Also included in Fig. 7b is the model by Ahrens and O'Keefe (1987) which was derived using impact velocities of 5 km/s and higher. Pierazzo et al. (1997) used impact velocities higher than those considered in the present study, except for their 10 km/s model. Hence, their results are shown in Fig. 7b by only one point, asterisk, at the impact velocity of 10 km/s.

The shock pressure along a non-vertical profile is not supposed
to be the same as the one along a vertical profile, largely because of
the pre-impact lithostatic pressure. As emphasized by Pierazzo
et al. (1997, 2008) the shock front in deeper regions appears rela-
tively symmetric around the impact point. Of course, it is not real-
istic in the case of an oblique impact (not studied here) and for the
shallowest angles where the surface significantly affects the shock
pressure decay. We have monitored the effect of the shockwave
propagation angle $\theta$ with values varying between 90° (vertical pro-
file) and 27° (Fig. 8). Similarly to Pierazzo et al. (1997), we did not
find a significant angle dependence on our results especially when
$\theta$ is ranging between 90° and 45°. For smaller values of $\theta$, the sur-
face effects appear to modify the shockwave propagation by reduc-
ing its intensity (Fig. 8a). Except in the mid field, where the *n* value
decreases from −0.6 to −1.31, and in the far field, where the *a*
value decreases from 2.54 to 1.93, the coefficients *a* and *n* from
our scaling laws do not change significantly with the angle
(Fig. 8b). This is particularly true in the near field where most of
the impact heating occurs.

Scaling laws have been used by many investigators (e.g. Senshu
et al., 2002; Tonks and Melosh, 1992, 1993; Watters et al., 2009;
<figure><img src="image_23.png" /><figcaption>**Fig. 7.** Dependence of regression parameters *a*, *n*, *c*, and *m* from Eqs. (4) and (5) on the impact velocity. Dots are based on hydrocode models and lines are regression fits, see Table 2.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 8.** The shockwave propagation as a function of the angle from the horizontal for the classical case ($V_{\text{imp}}$ = 10 km/s and $R_{\text{imp}}$ = 50 km). (a) Shock pressure decrease as a function of the distance to the impact site for propagation angles ranging between 27° and 90° and (b) dependence of regression parameters *a* and *n* (from Eqs. (4) and (5)) on the angle from horizontal.</figcaption></figure>

Roberts et al., 2009; Arkani-Hamed and Olson, 2010a,b; Arkani-Hamed and Ghods, 2011), mainly because they require a much smaller computer and much less computer time and the difference between a hydrocode model and a corresponding scaling model is minute. For example, Fig. 9 shows the 2D distribution of the peak shock pressure determined for an impactor of 100 km in diameter and an impact velocity of 10 km/s calculated using our scaling laws

in near field, mid field, and far field, and the parameter values from Table 2. The grid spacing is 2 km in radial direction and 0.03° in the colatitude direction. The entire computer time in a PC, CPU: 2.4 GHz, was only 16 s, which also calculated the 2D distribution of shock-induced temperature increase using Watters et al.'s (2009) foundering shock heating model. The computer time is substantially shorter than 48 h taken by our corresponding hydrocode
<figure><img src="image_25.png" /><figcaption>**Fig. 9.** 2D distribution of the peak shock pressure (left panel), and corresponding shock induced temperature increase (central panel) determined from our scaling laws for an impactor of 100 km in diameter and an impact velocity of 10 km/s calculated using the scaling laws in near field, mid field, and far field and the impact-induced temperature increase model from Watters et al. (2009) (see Eqs. (7)–(9)). The parameter values are based on regression lines from Fig. 7 and Table 2. The grid spacing is 2 km in radial direction and 0.03° in the colatitude direction. Only direct shock wave is considered. The black region in the mantle does not receive direct waves. The thin line extending from the impact site at the surface is the tangent line to the core. It shows that the part of the core receiving direct shock wave is not more than 25%. The color is saturated at higher values to better illustrate the regions with lower values. The right panel shows a close-up view of the corresponding temperature increase obtained from iSALE more than 30 min after the impact and illustrates the modifications that occur during the late stages of the impact process (excavation, central peak formation, slumping...). However, the intensity of impact heating is in agreement with our scaling laws. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

model using a CPU: 2.9 GHz laptop. This shows that it is feasible to determine impact heating during the accretion of a terrestrial planet using scaling laws, whereas it is almost impossible to adopt hydrocode simulations for that purpose.

During the decompression of shocked material much of the internal energy of the shock state is converted into heat (O'Keefe and Ahrens, 1977). Using thermodynamic relations, the waste heat used to heat up the impacted material can be estimated (Gault and Heitowit, 1963; Watters et al., 2009) and the corresponding temperature increase ΔT can be obtained. Hence, on the basis of our scaling laws it is possible to estimate the temperature increase as a function of depth below the impact site for impact velocities compatible with the accretionary conditions of terrestrial planets. These scaling laws can easily be implemented in a multi-impact approach (e.g. Senshu et al., 2002; Monteux et al., 2014) to monitor the temperature evolution inside a growing protoplanet whereas it is not yet possible to adopt hydrocode simulations for that purpose. For example, included in Fig. 9 is the impact induced temperature increase corresponding to the shock pressure shown in the figure. The temperature increase is determined on the basis of foundering model of Watters et al. (2009) using constant values for the acoustic velocity C (6600 m/s) and the parameter S (0.86) in their expressions:

$$ \Delta T(P) = \frac{P}{2\rho_0 S} (1 - f^{-1}) - (C/S)^2 [f - \ln f - 1] \quad (7) $$

$$ f(P) = -\frac{P}{\beta} \left( 1 - \sqrt{\frac{2P}{\beta} + 1} \right)^{-1} \quad (8) $$

$$ \beta = \frac{c^2 \rho_0}{2S} \quad (9) $$

with P the shock-increased pressure and ρ₀ the density prior to shock compression (see Table 1 for values).

Due to small size the impactor is not capable of increasing the lower mantle temperature of the model planet significantly, and only minor impact heating of the core has occurred. The thermal evolution model has to be combined to a topographical evolution model to account for the growth of the protoplanet as in Monteux et al. (2014). In these models, the impact angle (considered here as vertical) probably plays a key role because it influences both the morphology of the impact heating and the shape of the post-impact topography. A more elaborated scaling laws built upon 3D hydrocode models will be developed for that purpose in the near future.

It is worth emphasizing that our scaling laws, like those of others (Ahrens and O'Keefe, 1987; Pierazzo et al., 1997; Mitani, 2003), are derived from a few hydrocode models. Fig. 10 shows the profiles of the pressure along the axis of symmetry for comparison. The differences between the hydrocode model and the scaling law are small for the most part, but the exact scaling law differs by ~10 GPa for d/R<sub>imp</sub> = 2–3. This difference arises from the difficulty of correctly describe the near field as previously mentioned. Note that the interpolated model is in much better agreement with the hydrocode model.

A linear relationship has been proposed between the shock wave velocity V<sub>s</sub> and particle velocity V<sub>p</sub> on the basis of laboratory measurements (McQueen et al., 1967; Trunin, 2001)

$$ V_s = C + SV_p \quad (10) $$

where C is the acoustic velocity and S is a constant parameter. We estimate the acoustic velocity in the mantle of the model planet on the basis of our hydrocode models (Fig. 5a and b) using Eq. (10) and the Hugoniot equation

$$ P = \rho V_p V_s \quad (11) $$

where ρ (=3320 kg/m³) is the pre-shock density. Fig. 11 shows the variations of C with depth for models with impact velocities of 4–10 km/s and an impactor of 100 km diameter, using S = 1.2 which
<figure><img src="image_26.png" /><figcaption>**Fig. 10.** Shock pressure versus distance from the impact site normalized to the impactor radius. The hydrocode model is the output of the hydrocode simulation. In the exact model we have used $\alpha$ and $n$, while $\alpha$, $\beta$, $\lambda$ and $\delta$ are used in the interpolated model, see Table 2.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. 11.** The acoustic velocity determined using Eq. (10) for an impactor of 100 km diameter, and impact velocities ranging from 4 to 10 km/s. The S parameter in the equation is set to 1.2. The red dashed line is a line fit to C in the deeper region. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

is within the values proposed by the authors for dunite (e.g. Trunin, 2001). Adopting $S = 0.86$ (McQueen et al., 1967) does not change the results significantly, especially in the far field, where the acoustic velocity is less sensitive to particle velocity and linearly increases with depth. However C shows particle velocity dependence in the mid field and near field.

## 4. Conclusions

We have modeled the shock pressure and particle velocity distributions in the mantle of a Mars size planet using hydrocode simulations (iSALE-2D) for impact velocities of 4–10 km/s and protjectile diameters ranging from 100 to 400 km. We have extended Pierazzo et al.'s (1997) scaling laws to low impact velocities and also considered large impactor radii occurring during the formation of terrestrial planets. We propose three distinct regions in the mantle: a near field region, which extends to 1–3 times the projectile radius into the target, where the peak shock pressure and particle velocity decay very slowly with increasing distance,

a mid field region, which extends to ~4.5 times the impactor radius, where the pressure and particle velocity decay exponentially but moderately, and a more distant far field region where the pressure and particle velocity decay strongly with distance. The mid field–far field boundary is well constrained, whereas that of the near field–mid field is a relatively broad transition zone for the impact velocities examined.

## Acknowledgments

This research was supported by Agence Nationale de la Recherche (Oxydeep decision No. ANR-13-BS06-0008) to JM, and by Natural Sciences and Engineering Research Council (NSERC) of Canada to JAH. We gratefully acknowledge the developers of iSALE (www.isale-code.de), particularly the help we have received from Gareth S. Collins. We are also grateful to the two reviewers for very helpful suggestions.

## References

* Ahrens, T.J., O'Keefe, J.D., 1987. Impact on the Earth, ocean and atmosphere. Int. J. Impact Eng. 5, 13–32, 122.

* Arkani-Hamed, J., 2005. Magnetic crust of Mars. J. Geophys. Res. 110, E08005. http://dx.doi.org/10.1029/2004JE002397.

* Arkani-Hamed, J., Ghods, A., 2011. Could giant impacts cripple core dynamos of small terrestrial planets? Icarus 212, 920–934.

* Arkani-Hamed, J., Olson, P., 2010a. Giant impact stratification of the martian core. Geophys. Res. Lett. 37, L02201. http://dx.doi.org/10.1029/2009GL04141.

* Arkani-Hamed, J., Olson, P., 2010b. Giant impacts, core stratification, and failure of the martian dynamo. J. Geophys. Res. 115, E07012.

* Barr, A., Citron, R., 2011. Scaling of melt production in hypervelocity impacts from high-resolution numerical simulations. Icarus 211, 913–916.

* Benz, W., Cameron, A.G.W., Melosh, H.J., 1989. The origin of the Moon and the single-impact hypothesis III. Icarus 81, 113–131.

* Bierhaus, M., Wünnemann, K., Elbeshausen, D., 2012. Numerical modeling of basin-forming impacts: Implications for the heat budget of planetary interiors. Lunar Planet. Sci. 43. Abstract #2174.

* Bray, V.J. et al., 2008. The effect of target properties on crater morphology: Comparison of central peak craters on the Moon and Ganymede. Meteorit. Planet. Sci. 43, 1979–1992.

* Chambers, J.E., 2004. Planetary accretion in the inner Solar System. Earth Planet. Sci. Lett. 223, 241–252.

* Chambers, J.E., Wetherill, G.W., 1998. Making the terrestrial planets: N-body integrations of planetary embryos in three dimensions. Icarus 136, 304–327.

* Collins, G.S., Melosh, H.J., Ivanov, B.A., 2004. Modeling damage and deformation in impact simulations. Meteorit. Planet. Sci. 39 (2), 217–231. http://dx.doi.org/10.1111/j.1945-5100.2004.tb00337.x.

* Davison, T.M., Collins, G.S., Ciesla, F.J., 2010. Numerical modeling of heating in porous planetesimal collisions. Icarus 208, 468–481. http://dx.doi.org/10.1016/j.icarus.2010.01.034.

* Frey, H., 2008. Ages of very large impact basins on Mars: Implications for the late heavy bombardment in the inner Solar System. Geophys. Res. Lett. 35, L13203. http://dx.doi.org/10.1029/2008GL033515.

* Gault, D.E., Heitowit, E., 1963. The partition of energy for hypervelocity impact craters formed in rock. Paper Presented at 6th Hypervelocity Impact Symposium. Firestone Tire and Rubber Co., Cleveland (30 Apr. to 2 May).

* Holsapple, K.A., 1993. The scaling of impact processes in planetary sciences. Annu. Rev. Earth Planet. Sci. 21, 333–373.

* Ivanov, B.A., Melosh, H.J., Pierazzo, E., 2010. Basin-forming impacts: Reconnaissance modeling. Geol. Soc. Am. Spec. Pap. 465, 29–49.

* Keil, K. et al., 1997. Constraints on the role of impact heating and melting in asteroids. Meteorit. Planet. Sci. 32, 349–363.

* Kokubo, E., Ida, S., 1995. Orbital evolution of protoplanets embedded in a swarm of planetesimals. Icarus 114, 247–257.

* Kokubo, E., Ida, S., 1996. On runaway growth of planetesimals. Icarus 123, 180–191.

* Kokubo, E., Ida, S., 1998. Oligarchic growth of protoplanets. Icarus 131, 171–178.

* Kokubo, E., Ida, S., 2000. Formation of protoplanets from planetesimals in the solar nebula. Icarus 143, 15–27.

* Kraus, R.G., Senft, L.E., Stewart, S.T., 2011. Impacts onto H₂O ice: Scaling laws for melting, vaporization, excavation, and final crater size. Icarus 214, 724–738.

* Matsui, T., 1993. Early evolution of terrestrial planets: Accretion, atmospheric formation, and thermal history. In: Oya, H. (Ed.), Premitive Solar Nebula and Origin of Planets. Terra Scientific Publication Company, Tokyo, pp. 545–559.

* McQueen, R.G., Marsh, S.P., Fritz, J.N., 1967. Hugoniot equation of state of twelve rocks. J. Geophys. Res. 72, 4999–5036.

* McQueen, R.G. et al., 1970. The equation of state of solids from shock wave studies. In: Kinslow (Ed.), High Velocity Impact Phenomena. Academic, New York, pp. 294–491.

* Melosh, H.J., 1979. Acoustic fluidization: A new geologic process? J. Geophys. Res. 84, 7513–7520.
Melosh, H.J., 1989. Impact Cratering: A Geologic Process. Oxford Univ. Press, New York.

Mitani, N.K., 2003. Numerical simulations of shock attenuation in solids and reevaluation of scaling law. J. Geophys. Res. 108 (E1), 5003. http://dx.doi.org/10.1029/2000JE001472.

Monteux, J. et al., 2007. Thermo-mechanical adjustment after impacts during planetary growth. Geophys. Res. Lett. 34, L24201. http://dx.doi.org/10.1029/2007GL031635.

Monteux, J. et al., 2009. A model of metal-silicate separation on growing planets. Earth Planet. Sci. Lett. 287, 353–362. http://dx.doi.org/10.1016/j.epsl.2009.08.020.

Monteux, J., Jellinek, M., Johnson, C.J., 2013. Dynamics of core merging after a martian mega-impact. Icarus 226 (1), 20–32. http://dx.doi.org/10.1016/j.icarus.2013.05.008.

Monteux, J. et al., 2014. Can large icy Moons accrete undifferentiated? Icarus 237, 377–387. http://dx.doi.org/10.1016/j.icarus.2014.04.041 (15 July 2014).

Monteux, J. et al., 2015. Giant impacts, heterogeneous mantle heating and a past hemispheric dynamo on Mars. Phys. Earth Planet. Int. 240, 114–124. http://dx.doi.org/10.1016/j.pepi.2014.12.005.

O'Keefe, J.D., Ahrens, T.J., 1977. Impact-induced energy partitioning, melting, and vaporization on terrestrial planets. Lunar Planet. Sci. 8, 3357–3374.

Perret, W.R., Bass, R.C., 1975. Free-Field Ground Motion Induced by Underground Explosions. Sandia Report SAND74-0252.

Pierazzo, E., Vickery, A.M., Melosh, H.J., 1997. A re evaluation of impact melt production. Icarus 127, 408–423.

Pierazzo, E. et al., 2008. Validation of numerical codes for impact and explosion cratering: Impacts on strengthless and metal targets. Meteorit. Planet. Sci. 43 (12), 1917–1938.

Potter, R. et al., 2012. Estimating transient crater size using the crustal annular bulge: Insights from numerical modeling of lunar basin-scale impacts. Geophys. Res. Lett. 39, L18203.

Rafikov, R.R., 2003. The growth of planetary embryos: Orderly, runaway, or oligarchic. Astron. J. 125, 942–961.

Raymond, S.N., Quinn, T., Lunine, J.I., 2006. High-resolution simulations of the final assembly of Earth-like planets, Terrestrial accretion and dynamics. Icarus 183, 265–282.

Ricard, Y., Srámek, O., Dubuffet, F., 2009. A multi-phase model of runaway core-mantle segregation in planetary embryos. Earth Planet. Sci. Lett. 284, 144–150.

Roberts, J.H., Arkani-Hamed, J., 2012. Impact-induced mantle dynamics on Mars. Icarus 218, 278–289.

Roberts, J.H., Arkani-Hamed, J., 2014. Impact heating and coupled core cooling and mantle dynamics on Mars. J. Geophys. Res. Planets 119, 729–744. http://dx.doi.org/10.1002/2013JE004603.

Roberts, J.H., Lillis, R.J., Manga, M., 2009. Giant impacts on early Mars and the cessation of the martian dynamo. J. Geophys. Res. 114, E04009.

Schenk, P. et al., 2012. The geologically recent impact basins at Vesta's south pole. Science 336, 694–697.

Schmidt, R.M., Housen, K.R., 1987. Some recent advances in the scaling of impact an explosion cratering. Int. J. Impact Eng. 5, 543–560.

Senshu, H., Kuramoto, K., Matsui, T., 2002. Thermal evolution of a growing Mars. J. Geophys. Res. 107 (E12), 5118. http://dx.doi.org/10.1029/2001JE001819.

Thompson, S., Lauson, H., 1972. Improvements in the CHART D Radiation-Hydrodynamic Code III: Revised Analytic Equations of State. Sandia National Laboratory Report SC-RR-71 0714. 113 p.

Tonks, W.B., Melosh, H.J., 1992. Core formation by giant impacts. Icarus 100, 326–346.

Tonks, W.B., Melosh, H.J., 1993. Magma ocean formation due to giant impacts. J. Geophys. Res. 98, 5319–5333.

Trunin, R.F., 2001. Experimental Data on Shock Compression and Adiabatic Expansion of Condensed Matter. RFNC-VNIIEF, Sarov, p. 446. ISBN: 5-85165-624-7.

Watters, W.A., Zuber, M.T., Hager, B.H., 2009. Thermal perturbations caused by large impacts and consequences for mantle convection. J. Geophys. Res. 114, E02001.

Wetherill, G.W., Inaba, S., 2000. Planetary accumulation with a continuous supply of planetesimals. Space Sci. Rev. 92, 311–320.

Wetherill, G.W., Stewart, G.R., 1989. Accumulation of a swarm of small planetesimals. Icarus 77, 330–357.

Wünnemann, K., Ivanov, B.A., 2003. Numerical modelling of impact crater depth-diameter dependence in an acoustically fluidized target. Planet. Space Sci. 51, 831–845.

Wünnemann, K., Collins, G.S., Melosh, H.J., 2006. A strain-based porosity model for use in hydrocode simulations of impacts and implications for transient-crater growth in porous targets. Icarus 180, 514–552.
<figure><img src="image_28.png" /></figure>

# Shock wave propagation in layered planetary interiors: Revisited

J. Monteux<sup>a,*</sup>, J. Arkani-Hamed<sup>b,c</sup>

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

## ABSTRACT

While major impacts during late accretion of a Mars type planet occur on a differentiated body, the characteristics of the shockwave propagation are poorly known within these layered objects. Here, we use iSALE-2D hydrocode simulations to calculate shock pressure in a differentiated Mars type body for impact velocities ranging from 5 to 20 km/s, impactor radii ranging from 50 to 200 km, and different rheologies. To better represent the distribution of shock pressure as a function of distance from the impact site at the surface, we propose two distinct regions in the mantle: a near field region that extends to 7-15 times the projectile radius into the target, where the peak shock pressure decays exponentially with increasing the distance from the impact site, and a far field region where the pressure decays strongly with the distance following a power law. At the core-mantle boundary, the peak shock pressure increases from the mantle side to the core side. The refracted shockwave travels within the core where the shock pressure decreases following a second power law. In this study, we fit the output obtained from iSALE hydrocode simulations to determine scaling laws that illustrate the influence of the distance from the impact site, the ray angle, the target rheology, the impactor size and the impact velocity. Finally we combine these shock-pressure scaling laws with the formalism proposed by Watters et al. [2009] to determine the impact heating induced by large impacts within a differentiated Mars.

## 1. Introduction

Terrestrial planets are formed by accreting a huge number of planetesimals in the solar nebula [e.g., Chambers and Wetherill, 1998; Kokubo and Ida, 2000; Rafikov, 2003; Chambers, 2004; Raymond et al., 2006]. The ultimate goal is to investigate the thermal evolution of a growing terrestrial planet while taking into account the impact heating, heating by short-lived and long-lived radioactive elements, and by the core formation process. Because partial melting and core formation occur when a protoplanet's radius exceeds ~2700 km [e.g., Senshu et al., 2002], about 47% of Mars' and 90% of the Earth's and Venus' masses are expected to be accreted in the presence of a liquid core. To estimate the thermal state of a growing planet it is necessary to determine the impact heating of its solid mantle and liquid core by each impactor.

The impact-induced shock pressure and temperature increase inside a planet has been investigated using elaborated numerical models developed to better understand shock physics [e.g., Ahrens and O'Keefe, 1987; Pierazzo et al., 1997; Mitani, 2003; Wünnemann et al., 2006; Wünnemann et al., 2008; Ivanov et al., 2010; Barr and Citron, 2011; Kraus et al., 2011; Bierhaus et al., 2012; Leinhardt and Stewart, 2009]. However, the numerical solutions demand considerable computer capacity and time and are not practical for investigating a huge number of impacts that occur during the growth of a planet. For example, the

formation of a planet like Mars requires about $3 \times 10^5$ planetesimals of ~100 km in diameter. On the other hand, the scaling laws derived on the basis of numerical solutions of the shock dynamic equations [e.g., Ahrens and O'Keefe, 1987; Pierazzo et al., 1997; Mitani, 2003; Monteux and Arkani-Hamed, 2016] require remarkably less computer power and time, and have been used by many investigators [e.g., Tonks and Melosh, 1992, 1993; Watters et al., 2009; Roberts et al., 2009; Arkani-Hamed and Olson, 2010a, 2010b; Ghods and Arkani-Hamed, 2011], mainly because of their simplicity and partly because the difference between their results and those obtained by the numerical simulations of the shock dynamic equations is likely within the tolerance of the errors introduced due to the uncertainty of pertinent physical parameters. We note that the existing scaling laws are applicable for uniform mantle models of terrestrial planets and cannot be used for planets consisting of a solid silicate mantle and a liquid iron core with drastically different physical properties. Extra formulations are required to relate the shock pressure and particle velocity at the base of the mantle to those at the top of the core. Arkani-Hamed and Ivanov [2014] derived the required boundary conditions at the core-mantle boundary (CMB) of a Mars type planet for a vertical impact on the basis of Hugoniot equations, the equality of the reflected and the refracted shock pressures, the continuity of the vertical component of particle velocity, and the Snell's law relating the incident angle to the refracted angle of a shock ray at the boundary. This resulted in a set of coupled 8

<footnote><sup>a</sup> Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, F-63000 Clermont-Ferrand, France</footnote>

<footnote><sup>b</sup> Department of Physics, University of Toronto, Toronto, Canada</footnote>

<footnote><sup>c</sup> Department of Earth and Planetary Sciences, McGill University, Montreal, Canada</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: j.monteux@opgc.univ-bpclermont.fr (J. Monteux).</footnote>
trigonometric equations to be solved simultaneously for each shock ray. At each point of the core mantle boundary, with spacing of 1° colatitude, the equations were solved iteratively until the difference between the two successive solutions were within a set error limit. Although such a simulation is much faster than the corresponding hydrocode simulation, it still takes considerable time to determine impact heating in a growing planet that requires several hundred thousand impacts during the accretion. Recently, Monteux and Arkani-Hamed [2016] derived scaling laws in the mantle of a Mars size planet on the basis of systematic hydrocode modeling, using impact velocities ranging from 4 to 10 km/s for impactor sizes ranging from 50 to 200 km in radius, and two different mantle rocks of dunite and peridotite overlying the iron core. The ANEOS type equations of state were used for the dunitic mantle rocks and the iron core.

In the present study, we adopt a new technique based on hydrocode models by Monteux and Arkani-Hamed [2016] while deriving a new scaling law for the iron core to assess shock pressure and the shock-related temperature increase in the entire Mars type planet. A given shock front propagates from the impact site at the surface of the planet down to the CMB. It then partly reflects back to the mantle, and partly enters the core and continues propagating until it impinges the CMB in the antipodal hemisphere, where it partly enters the antipodal mantle and partly reflects and propagates in the antipodal core. The second section characterizes the shockwave propagation within a differentiated Mars size planet for different sets of parameters (impactor radius, impact velocity, and mantle rheology). The third section compares the results from the scaling laws with those obtained using hydrocode simulations. In the fourth section, impact-heated models determined on the basis of the scaling laws to demonstrate the feasibility of the models derived by the scaling laws. Discussions and concluding remarks are relegated to the fifth, final section.

## 2. Shockwave propagation inside a Mars-type planet

Fig. 1 shows the 2D axi-symmetric total pressure (shock + litho-static) distribution inside a Mars size body of radius 3400 km with a core radius of 1700 km, created by an impactor of 200 km in radius and at an impact velocity of 10 km/s obtained by hydrocode simulation. The shock pressure produced near the surface decays as the shock wave propagates downward until the wave impinges the CMB, where it partitions into two waves, a reflected wave that propagates outward in the mantle, and a refracted wave that enters the core and propagates downward. The reflected wave arriving at a given point in the mantle travels much longer, hence is always weaker than the direct wave that propagates directly from the impact site to that point. The peak pressure at the point is produced by the direct wave. The effects of the two reflected waves, one at the surface and the other at the CMB are not considered on the temperature of the mantle in this study (see below). An abrupt pressure jump is created at the top of the core by the refracted wave. The pressure then continues decreasing while the wave propagates in the core.

A total of 11 hydrocode models were simulated by Monteux and Arkani-Hamed [2016] that included impact velocities of 4 to 10 km/s at increments of 1 km/s for an impactor of 50 km in radius, and the impactor radii of 50 to 200 km at increments of 50 km for an impact velocity of 10 km/s. Two mantle rock types of dunite and peridotite were considered for a Mars-type planet. The authors considered shock wave propagation in the mantle of the impacted hemisphere. Here we extend the models allowing the wave enter the core, hence providing a means to estimate the characteristics of shock pressure and the related temperature increase throughout the mantle and the core.

### 2.1. The model set up

At the end of the planetary accretion, large impacts between planetsimals and protoplanets are frequent. The amount of energy

dissipated during this regime of planetary evolution is large enough to significantly melt the growing planets and efficiently separate the metal phase from the silicate phase, resulting in the core formation. We consider a large impact between an undifferentiated bolide and a differentiated Mars size planet using numerical hydrocode simulations.

### 2.1.1. Physical model

We investigate the shock pressure inside a Mars size model planet for impact velocities $V_{\text{imp}}$ of 5 to 20 km/s and impactors radii $R_{\text{imp}}$ of 50 to 200 km. Such impactors are likely to create large impact basins with sizes comparable to the giant impact basins of Mars that were formed between 4.3 and 4 Ga [Frey, 2008; Robbins et al., 2013]. The impactors are large enough to have important consequences on their target. However the impactors might be too small to have experienced a complete melting and full metal/silicate separation before the impact unless they formed during the early Solar System evolution [Ricard et al., 2017] Impacts involving differentiated impactors will lead to the fragmentation of the impactor's core before sinking through the target's mantle [Kendall and Melosh, 2016]. In our models, an impactor is simplified by a spherical body of radius $R_{\text{imp}}$ with uniform dunitic composition. Since we do not consider a realistic impactor with a metallic core, we underestimate the actual impactor mass and, as a consequence, the kinetic energy available for impact heating. In our models, the impactor and the impacted mantle have the same properties (composition, strength, rheology, and equation of state).

The ultimate goal is to determine the impact heating during the accretion of terrestrial planets such as Mars and the Earth. The huge number of impacts during accretion makes it impractical to consider oblique impacts for two main reasons: it requires formidable computer time, but more importantly because of the lack of information about the impact direction, i.e. the impact angle relative to vertical and azimuth relative to north. The peak shock pressure produced along the vertical direction beneath the impact site by an oblique impact appears to depend approximately on the sine of the impact angle relative to horizon (see Fig. 3 of Pierazzo and Melosh, 2000a, 2000b). This implies a shock pressure reduction by a factor of ~0.71 for the most probable impact angle of 45° [Shoemaker, 1962] compared to that of a vertical impact, assuming that all other parameters are equal. For example an oblique impact velocity of about 15 km/s produces almost the same shock pressures as that produced by a vertical impact velocity of about 11 km/s. The difference between these two velocities is well within the variations of the impact velocities of the N-body accretion models of terrestrial planets [e.g., Agnor et al., 1999; Chambers, 2013]. Studies of more realistic accretion scenarios, which take into account the oblique impacts and allow escape of material, conclude that the final state of an accreted body is less sensitive to the details of the collision [Agnor and Asphaug, 2004; Kokubo and Genda, 2010]. Therefore, we consider only head-on collisions (vertical impacts) to model the thermo-mechanical evolution during an impact between a differentiated Mars size body and a large uniform impactor.

As emphasized by Monteux and Arkani-Hamed [2016], the rheology of the impacted body plays a key role on the propagation of the shockwave and on the subsequent impact heating. However, the rheology of protoplanets is difficult to constrain precisely and is strongly governed by its thermal evolution. Depending on its accretion rate, on its composition in radiogenic elements and on the occurrence of large to giant impacts, growing Mars-sized protoplanets likely consist of a partially to fully molten mantle overlying a liquid metallic core [Kaula, 1979; Solomatov, 2000; Senshu et al., 2002; Ricard et al., 2009]. In the following study we consider that the rheology of the impacted protoplanets ranges between a “hydrodynamic model” with no-strength where both the dunitic mantle and the metallic core are inviscid materials, and a “model with strength” where the metallic core is an inviscid fluid while the dunitic mantle is treated using a pressure-, temperature-, and damage- dependent strength model.

Damage describes the transition from intact to fractured strength.
<figure><img src="image_30.png" /></figure>

**Fig. 1.** Time evolution of the total pressure (shock + lithostatic) inside a differentiated Mars size body during the first 600 s after the impact, caused by a 200 km radius undifferentiated impactor at 10 km/s impact velocity. In our models, the Martian core is liquid. The lithostatic pressure is adopted from Monteux and Arkani-Hamed [2016] and the pre-impact temperature profile is illustrated in Fig. 4 (cold profile). The four left panels represent 4 specific times (from top to bottom): t = 0 s, t = 50 s, t = 100 s and t = 225 s with a focus on the impacted mantle. The four right panels represent 4 times (from top to bottom): t = 225 s, t = 350 s, t = 475 s and t = 600 s with a focus on the impacted core. Note the change in the pressure scales between the left and right panels.

Under a constant pressure, damage begins to accumulate when the material begins to flow plastically. Damage is maximum when the material is completely fractured. Among the damage models available in iSALE, we adopt the Ivanov damage model [Ivanov et al., 1997], which prescribes damage as a function of plastic strain. In this model the plastic strain at failure is an increasing function of the total pressure P (lithostatic pressure plus the shock-induced pressure). This model allows rapid brittle failure in tensile and low-pressure modes, and semi-brittle failure at high pressures. The damage model from Ivanov et al. [1997] is simpler than the model developed by Collins and Melosh [2014] where the plastic strain at failure is a function of the brittle-ductile and brittle-plastic transition pressures, both these parameters being poorly constrained in the early protoplanets context.

The thermal softening is an efficient process for large-scale events [Potter et al., 2012]. Depending on the temperature and more specifically on the melt fraction, the impacted material can behave as a solid material or as fluid. The transition from solid-like to fluid-like behavior occurs at a critical melt fraction of ≈40% for early mantle material (Solomatov, 2000). As it depends on the temperature, the thermal softening process also strongly depends on the pre-impact temperature that is unfortunately poorly constrained in the Martian case. In our iSALE models, we use the Ohnaka thermal softening model where the shear strength of rocks depends on temperature, and reduces to zero at the melting temperature.

In the “model with strength”, we approximate the thermodynamic response of both the iron and dunitic material using the ANEOS equation of state [Thompson and Lauson, 1972; Benz et al., 1989]. To make our models as simple as possible we do not consider here the effects of porosity or low density weakening which are likely to be negligible in the deep interior of a planet. However, as a first step toward more

realistic models, we investigate the influence of acoustic fluidization, thermal softening and damage. The simulation parameters are listed in Table 1.

### 2.1.2. Numerical model

We use the iSALE-2D axisymmetric hydrocode, which is a multi-rheology, multi-material hydrocode, specifically developed to model impact crater formation on a planetary scale [Amsden et al., 1980; Collins and Melosh, 2014, Davison et al., 2010]. To enable a comparison of our models, for all our simulations, we use a 2 km grid resolution computational domain. The horizontal and vertical number of cells per projectile radius (cppr) is maintained constant within the impacted mantle (=1700) and within the impacted core (=850) (see Table 2). The number of cppr varies only within the impactor when varying the impactor radius (cppr = 25 for R<sub>imp</sub> = 50 km and cppr = 100 for R<sub>imp</sub> = 100 km). Such spatial resolutions should lead to an underestimation of the peak shock pressure ranging between 5% (close to the impact site) and 10% (far from the impact site) according to resolution studies [e.g., Pierazzo et al., 2008; Wünnemann et al., 2008; Monteux and Arkani-Hamed, 2016].

iSALE requires the length of each time step (dt) and the maximum time step (dt<sub>max</sub>). iSALE limits the time step by the speed of sound according to the Courant–Friedrichs–Lewy (CFL) condition while dt<sub>max</sub> limits a timestep if it is smaller than the CFL-limited timestep. In our models, we use a dt<sub>max</sub> value of 0.05 s. The minimum post impact monitoring time is set to the time needed by the shockwave to travel from the impact site through the entire core of the impacted planet.

In our iSALE models, we investigate the effects of the ray angle $\beta$ relative to the symmetric axis (Fig. 2). To compare the shock wave propagation along each path illustrated in Fig. 2, we developed a
<table><caption>**Table 1**

Typical physical parameter values for numerical hydrocode models.</caption><tr><td>Target radius</td><td>R</td><td>3400 km (Rivoldini et al., 2011)</td></tr><tr><td>Target core radius</td><td>R<sub>core</sub></td><td>1700 km (Rivoldini et al., 2011)</td></tr><tr><td>Impactor radius</td><td>R<sub>imp</sub></td><td>50–200 km</td></tr><tr><td>Impact velocity</td><td>V<sub>imp</sub></td><td>5–20 km/s</td></tr><tr><td>Mantle properties (Dunite)</td><td></td><td>(Benz et al., 1989)<br/>Model with strength</td></tr><tr><td>Initial density</td><td>ρ<sub>m</sub></td><td>3314 kg/m<sup>3</sup></td></tr><tr><td>Heat capacity</td><td>c<sub>p</sub></td><td>1200 J/kg/m<sup>3</sup></td></tr><tr><td>Equation of state type</td><td></td><td>ANEOS for dunite (also stands for hydrodynamic models)</td></tr><tr><td>Poisson ratio</td><td></td><td>0.25</td></tr><tr><td>Strength Model<br/>(iSALE parameters)</td><td></td><td>Rock<br/>(Y<sub>i0</sub> = 10 MPa,<br/>μ<sub>i</sub> = 1.2, Y<sub>lm</sub> = 3.5 GPa)</td></tr><tr><td>Acoustic Fluidization Model<br/>(iSALE parameters)</td><td></td><td>Block<br/>(t<sub>off</sub> = 16 s, c<sub>vib</sub> = 0.1 m/s,<br/>vib<sub>max</sub> = 200 m/s)</td></tr><tr><td>Damage Model<br/>(iSALE parameters)</td><td></td><td>[Ivanov et al., 1997]<br/>(ε<sub>fb</sub> = 10<sup>-4</sup>, B = 10<sup>-11</sup>,<br/>P<sub>c</sub> = 3 × 10<sup>8</sup> Pa)</td></tr><tr><td>Thermal softening model<br/>(iSALE parameters)</td><td>ξ</td><td>1.2</td></tr><tr><td></td><td>T<sub>m</sub></td><td>1375 K [Davison, 2010]</td></tr><tr><td>Low density weakening and<br/>porosity models</td><td></td><td>None</td></tr><tr><td>Core properties (Iron)</td><td></td><td>(Thompson, 1990)<br/>All models</td></tr><tr><td>Initial density</td><td>ρ<sub>c</sub></td><td>7840 kg/m<sup>3</sup></td></tr><tr><td>Heat capacity</td><td>c<sub>p</sub></td><td>800 J/kg/m<sup>3</sup></td></tr><tr><td>Equation of state type</td><td></td><td>ANEOS for iron</td></tr><tr><td>Damage, low density weakening<br/>thermal softening and<br/>porosity models</td><td></td><td>None</td></tr></table>

<table><caption>**Table 2**

Typical numerical parameter values for numerical hydrocode models.</caption><tr><td>Horizontal cells</td><td>2000</td></tr><tr><td>Vertical cells</td><td>4000</td></tr><tr><td>Grid spacing</td><td>2 km</td></tr><tr><td>Cells per projectile radius (CPPR)</td><td></td></tr><tr><td>Impactor</td><td>25 (for R<sub>imp</sub> = 50 km) to 100 (for R<sub>imp</sub> = 200 km)</td></tr><tr><td>Target mantle</td><td>1700</td></tr><tr><td>Target core</td><td>850</td></tr><tr><td>Maximum time step (dt<sub>max</sub>)</td><td>0.05 s</td></tr></table>

<figure><img src="image_31.png" /><figcaption>**Fig. 2.** The rays along which the shock pressure is monitored. The different rays illustrated here with thin solid lines correspond to different values of β ranging between 0° (vertical) and 80°. The black thick line at D<sub>Z1</sub> represents the boundary between zones 1 and 2 (in this figure, the boundary is placed arbitrarily). Between the impact site and D<sub>Z1</sub>, the pressure decreases exponentially while beyond D<sub>Z1</sub>, the pressure decrease follows a power law (see Section 2.3). The numbers 30, 45, and 60 show the β values in degrees.</figcaption></figure>

numerical method to extract the shock pressure and temperature at a fixed position from all our simulations. As a consequence, instead of using tracers, the impact-induced pressure and temperature fields are extracted from a cell-centered Eulerian grid points [Monteux and Arkani-Hamed, 2016]. Hence, the distance D corresponds to the straight-line distance between the impact site and the grid point where the peak pressure is located. It approximates the distance along a shock ray (see the Discussion section). Hence, we refer to the straight line as the ray from now on. We also investigate the effects of the angle β of the straight line, i.e., the so called ray, relative to the symmetric axis (Fig. 2). We emphasize that a ray does not enter the core when the ray angle β is larger than 30° in our model planet.

## 2.2. The reference models

Fig. 3 (left) shows the peak shock pressure distribution along the axis of symmetry (β = 0°) in a model with R<sub>imp</sub> = 50 km and V<sub>imp</sub> = 10 km/s. From Fig. 3, we can identify 2 zones of particular interest: Zone 1 where the peak shock pressure P<sub>Z1</sub> decreases exponentially with the distance from the impact site, and Zone 2 where the peak shock pressure P<sub>Z2</sub> still decreases significantly though less rapidly with the distance. As we will see later, these 2 zones are particularly distinguishable for small values of R<sub>imp</sub>. Considering V<sub>imp</sub> = 10 km/s is the average asteroid impact velocity on Mars [Bottke et al., 1994], we will define the case with R<sub>imp</sub> = 50 km and V<sub>imp</sub> = 10 km/s as our reference model. Fig. 3 (right) illustrates that considering a differentiated impactor (with a dunicut mantle and a liquid iron core) slightly increases the peak shock pressure magnitude during the propagation of the shockwave within the impacted planet because of the different mass repartition in the impactor. However, considering a differentiated impactor affects neither the pressure decay in 2 zones nor the position of the limit between the two zones.

The shock wave generated by an impactor propagates toward the center of the planet, while the shock pressure decays as the wave deposits energy in the form of irreversible work into the target [Ahrens and O'Keeffe, 1987]. Hence, the ability of the impacted material to deform when experiencing large shock pressure is of the first importance for the characterization of the peak shock pressure evolution. Fig. 3 illustrates the influence of the mantle rheology on the shock wave propagation. The peak shock pressure experienced within the mantle decreases more rapidly in the model with strength than in the hydrodynamic model. In our models, the acoustic fluidization and the damage play minor roles on the shockwave propagation compared to that of the strength, because the lithostatic pressure is too high to allow fracturing of rocks.

If strength is present, shock waves attenuate faster because the rarefaction waves travel faster in material with strength and, therefore, cause a more rapid decay behavior [Curran et al., 1977; Bierhaus et al., 2013]. We note that the speed of sound in a strength-less material is given by $c = \sqrt{K/\rho}$ where K is the bulk modulus and ρ the density. If the material has some strength, the longitudinal wave speed is determined by $c_L = \sqrt{(K + 4/3\mu)/\rho}$ where μ is the shear modulus. Because $c_L > c$, the speed of the first, direct wave is larger when the strength is accounted for than in the hydrodynamic case, hence the shock waves attenuate faster in the models with strength [Curran et al., 1977; Bierhaus et al., 2013]. Ultimately, if the speed of the rarefaction wave depends on strength, then differences in strength should result in different attenuation behaviors. However this point is beyond the scope of our study. Outside Zone 1, the peak shock pressure continues to decrease with distance from the impact site and both models show a similar behavior: a power-law decrease followed by a pressure jump at the CMB. The more moderate decay in Zone 2 might be the consequence of a solid-state phase transformation that is considered in the ANEOS for dunite [Collins and Melosh, 2014]. This process does not occur in the hydrodynamic models, because the shock pressure does not drop below the critical pressure.
<figure><img src="image_32.png" /><figcaption>Fig. 3. Left: Shock pressure along the axis of symmetry (β = 0°) as a function of the distance from the impact site. In this (reference) case, R<sub>imp</sub> = 50 km and V<sub>imp</sub> = 10 km/s. The green, blue and orange dashed lines represent the fits obtained using Eqs. (1a), (1b) and (1c) respectively. Black stars represent the results from Pierazzo et al. [1997] for V<sub>imp</sub> = 10 km/s and R<sub>imp</sub> ranging between 0.4 and 10 km. Right: Same case as left panel (Model with strength) considering a dunitic (black solid line) impactor and a differentiated impactor with a 24 km core radius (red solid line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Included in Fig. 3 are the results from Pierazzo et al. [1997] for V<sub>imp</sub> = 10 km/s and R<sub>imp</sub> ranging between 0.4 and 10 km, where the strength was not included. Our results for the hydrodynamic model, and the results from Pierazzo are in good agreement (Fig. 3). The small differences between our results and the results from Pierazzo et al. [1997] are plausibly the direct consequence of using different cppr values (for R<sub>imp</sub> = 50 km, cppr = 25 while cppr = 20 for Pierazzo et al., 1997). Fig. 3 illustrates that accounting for strength significantly reduces the intensity of the shock pressure. This indicates that building more sophisticated models is necessary to estimate the shock pressure evolution after a large impact, as suggested by Monteux and Arkani-Hamed [2016].

For large impacts, the thermal softening is an efficient process that may influence the shockwave propagation. Indeed, accounting for thermal softening changes the strength and, thus, the speed of the rarefaction wave, which may result in a different attenuation behavior. This process strongly depends on the pre-shocked temperature profile that is still poorly constrained during the early evolution of terrestrial planets. To illustrate the influence of the thermal softening, we have run two models with different pre-impact spherically symmetric temperature distributions, characterized by a cold boundary layer at the surface overlying a convective-type, adiabatic temperature within the mantle. The difference between these two models is the surface temperature that is 293 K for the cold case and 500 K for the hot case. Fig. 4

shows that the thermal softening (as damage and acoustic fluidization) plays a minor role on the shockwave attenuation after a large impact. Accounting for thermal softening leads to a peak shock pressure decrease in Zone 1 and to an increase of the depth of Zone 1-Zone 2 transition. An initially hotter mantle only slightly enhances these 2 effects.

## 2.3. Scaling Laws

We have slightly modified the procedure by Monteux and Arkani-Hamed [2016]. By combining their near field and midfield, we reduce the number of parameters and make it easy to adopt the scaling laws in the impact heating of a terrestrial proto-plant during its accretion. We determine new scaling laws of shock pressure P for deriving universal scaling laws. From Fig. 3, we have identified 2 zones of particular interest: Zone 1 (associated to a peak shock pressure P<sub>Z1</sub>) and Zone 2 (associated to a peak shock pressure P<sub>Z2</sub>). This dichotomy is particularly pronounced in the models with strength.

We suggest two distinct scaling laws along the axis of symmetry to characterize the peak shock pressure in the model planet:

$$P_{Z1} = A_0 \exp(-A_1 D/R_{imp}) \quad (1a)$$

$$P_{Z2} = A_2 (D/R_{imp})^{-A_3} \quad (1b)$$

<figure><img src="image_33.png" /><figcaption>Fig. 4. Left: Pre-impact temperature profiles and solidus used for models including thermal softening. The cyan line represents the “cold” profile with a pre-impact surface temperature of 273 K while the red line represents the “hot” profile with a pre-impact surface temperature of 500 K. Right: Shock pressure along the axis of symmetry (β = 0°) as a function of the distance from the impact site. The black line represents the reference case (same as Fig. 3 for model with strength, R<sub>imp</sub> = 50 km and V<sub>imp</sub> = 10 km/s). The cyan and red lines represent the shock pressure decay for the models accounting for thermal softening considering the “cold” and the “hot” pre-impact temperature profile respectively illustrated on the left figure. The green, blue and orange dashed lines represent the fits obtained using Eq. (1a), (1b) and (1c) respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
where D is the distance from the impact site, R<sub>imp</sub> is the impactor radius. A<sub>0</sub> is the maximum shock pressure and A<sub>1</sub>, A<sub>2</sub>, and A<sub>3</sub> are constant obtained by fitting the results from the hydrocode models. All these parameters are impact velocity and impactor radius dependent. They also depend on the ray angle β. Zone 1 corresponds to a hemisphere centered at the impact site with radius D<sub>Z1</sub>. Zone 2 corresponds to the mantle shell beyond D<sub>Z1</sub>. Zone 2 lies between D<sub>Z1</sub> and the CMB for the ray angle β < 30°. The zone is limited by the surface of the impacted planet for β > 30° (see Fig. 2).

The scaling law used to fit the peak shock pressure decay within Zone 1 is different from the classical power-law that usually uses other empirical or semi-empirical relations in impact studies [e.g. Ahrens et al., 1977]. However it has two advantages: 1) with this formalism the peak shock pressure does not go to infinity as D/R<sub>imp</sub> goes to 0 compared to other power law forms [e.g., Ruedas, 2017], and 2) the peak shock pressure decays relatively slowly for D/R<sub>imp</sub> < 3–4, which includes the “isobaric core” with slowly decaying shock pressure and the “pressure decay regime” regions proposed by Pierazzo et al. [1997]. On the other hand, the Zone 2 represents the far field region in the mantle where the shock pressure decays monotonically with distance from the impact site.

As the shock wave propagates in the mantle it creates a high pressure behind the shock front, which is determined on the basis of the Hugoniot equations. When the shockwave crosses the CMB, the peak shock pressure jumps from P<sub>m</sub> (immediately above the CMB) to P<sub>c</sub> (immediately below the CMB) because of the abrupt changes in the physical parameters, as previously identified [e.g., Ivanov et al., 2010; Arkani-Hamed and Ivanov, 2014]. After this sudden increase, the peak shock pressure continues to decrease as the wave propagates inside the core, which is characterized by a separate scaling law:

$$P_{\text{core}} = A_4 (D/R_{\text{imp}})^{-A_5} \qquad (1c)$$

where A<sub>4</sub> and A<sub>5</sub> are constants. Fitting the peak shock pressure from our numerical models along the ray paths shown in Fig. 2 enables us to fully characterize the pressure increase experienced within the martian mantle and core. We note that a vertical impact on a spherically symmetric planet adopted in this study artificially enhances the shock pressure along the axis of symmetry in the core, hence Eq. (1c) overestimate the pressure near the axis of symmetry, while such an effect does not occur in the mantle (see the Discussion section below). Eqs. (1a), (1b) and (1c) will be used to fit the results from both the models with strength and the hydrodynamic models. The values of the different coefficients resulting from the corresponding fits are provided in the following sections. We would like to emphasize here that these scaling laws are only empirical and can be implemented in large-scale mantle models [e.g. Monteux and Arkani-Hamed, 2014] or in accretionary models [e.g. Monteux et al., 2014]. These scaling laws are determined by least square fitting to the results from the iSALE hydrocode models.

## 2.4. Shock wave propagation in the mantle: effects of the ray angle

During an impact, the maximum shock pressure is experienced at the impact site on the surface. Then the impact energy dissipates as the shockwave propagates within the model planet. Fig. 5 shows that the maximum shock pressure A<sub>0</sub>, in either the models with strength or the hydrodynamic models, increases with the impact velocity and the impactor radius. The peak shock pressure is nearly constant within an almost hemispheric region of a radius 1–2 times the impactor radius, centered at the impact site (see Fig. 3), which is coined as “isobaric core” by Pierazzo et al. [1997]. We note that A<sub>0</sub> is the shock pressure at the impact site and it remains the same for all of the ray angles considered. It is used in our scaling law (Eq. (1a)) to obtain the A<sub>1</sub> values from our models (Fig. 6). Fig. 5 shows the maximum shock pressure predicted from the planar impact model [Melosh, 1989]. As detailed by Osinski and Pierazzo [2012], the planar impact approximation yields

an upper limit of the pressure expected during an impact event. Hence, it is not surprising that our A<sub>0</sub> values are slightly smaller than the results from the planar impact approximation.

Previous studies [Pierazzo et al., 1997; Monteux and Arkani-Hamed, 2016] monitored the influence of the ray angle β of the shock wave, and found no significant angle dependence in their results for β ranging between 0° (vertical) and 45°. Here we extend the investigation to β values up to 80°. A<sub>1</sub> is the exponential factor (Eq. (1a)), characterizing the ability of the upper mantle to attenuate the shockwave as it propagates. For large A<sub>1</sub> values, the peak shock pressure decreases rapidly with distance. Hence, higher A<sub>1</sub> values are expected for the models with strength than for the hydrodynamic models. In Fig. 6a and b, we present the A<sub>1</sub> values obtained by fitting the results from iSALE models to Eq. (1a) along different ray angles illustrated in Fig. 2. A<sub>1</sub> increases as a function of β, indicating that the peak shock pressure decreases more efficiently with distance at large ray angles. This behavior is common to both models with strength and hydrodynamics models. However, the increase is not significant for β values smaller than 45°, as also mentioned by Pierazzo et al. [1997]. We note that the rays impinging the core mantle boundary have β values smaller than 30°, indicating that propagation of shock waves inside the core is almost independent of the β values. Fig. 6a shows that A<sub>1</sub> is weakly affected by the impactor size for the range of impactors used in our study (i.e. 50 < Rimp < 200 km). For larger impactors (not studied here), one can expect an influence of the impacted protoplanet curvature on the shockwave propagation. However, Fig. 6b shows that A<sub>1</sub> is more sensitive to the impact velocity. Pierazzo et al. [1997] showed that the power-law exponent used to fit the angle dependency (for β ranging between 0° and 45°) was increasing with the velocity. From our results (for β ranging between 0° and 80°) it is difficult to extract such a tendency.

For the large ray angles, the path followed by the shockwave is shallower (Fig. 2). The shock wave propagating through shallower regions will cross material where the pressure and the uncompressed target density are smaller in average compared to shock wave propagating along the axis of symmetry. As a consequence, along small ray angle paths, the deformation is more difficult than along large ray angle paths. In addition, in shallower regions a significant part of the impact energy is dissipated during processes such as the excavation or the post-impact rebound. Finally, the interference of the direct pressure wave with the rarefaction wave that is reflected at the surface reduces the amplitude of the shock pressure [Melosh, 1989; Arkani-Hamed, 2005; Louzada and Stewart, 2009]. The consequence of these factors is that the isobaric surfaces are slightly elliptical with vertical major axis rather than spherical.

As previously mentioned, the shock wave decay rate with distance should be more rapid for the model with strength compared to the hydrodynamic models (as illustrated in Fig. 3). Fig. 6a and b confirm this point by showing that A<sub>1</sub> values are larger for models with strength than for hydrodynamic models for all β values. The correlation coefficients obtained by fitting our hydrocode data with the scaling laws from Eqs. (1a), (1b) and (1c) show that our choice is relevant. The correlation coefficients are close to 0.9 for all the range of impactor radii and impact velocities used in our models (except for V<sub>imp</sub> = 20 km/s). The correlation coefficients are even higher for β values smaller than 65°. The correlation coefficients are also larger in the models with strength than in the hydrodynamic models because the dichotomy between Zone 1 and Zone 2 is more pronounced in the models with strength than in the hydrodynamic models (see Fig. 3).

In our study, the distance D from the impact site is normalized to the impactor radius R<sub>imp</sub> to eliminate the dependence on the size of the projectile. Fig. 7 illustrates peak shock pressure profiles along a ray angle of β = 36° meaning that the ray does not cross the CMB (see Fig. 2). Fig. 7 shows that the boundary between Zone 1 and Zone 2 depends on both the impactor size and the impact velocity, hence the total pressure (shock plus lithostatic) through the Equation of State adopted in the hydrocode models. This suggests a relation between the
<figure><img src="image_34.png" /><figcaption>Fig. 5. The maximum peak shock pressure $A_0$ (Eq. (1a)) as a function of the impactor radius $R_{imp}$ (left), and the impact velocity $V_{imp}$ (right). The dotted lines represent the result from the planar impact approximation [Melosh, 1989] considering that both the impactor and the target are made of dunite (with $\rho_0 = 3314.9$, $C_0 = 6290$ m/s and $S = 0.86$ for both iSALE models and planar impact approximation).</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. 6. a: $A_1$ as a function of the ray angle $\beta$ (Eq. (1a): $P_{Z1} = A_0 \exp(-A_1 D/R_{imp})$) for different impactor radii, but the same impact velocity ($V_{imp} = 10$ km/s). $A_1$ values are obtained by fitting the peak shock pressure in Zone 1 (close to the impact site) from iSALE models. The correlation coefficient ($r$), between the values obtained by the iSALE model and by the scaling law is mentioned in the legend box.</figcaption></figure>

<caption>b: $A_1$ as a function of the ray angle $\beta$ (Eq. (1a): $P_{Z1} = A_0 \exp(-A_1 D/R_{imp})$) for different impact velocities ($R_{imp} = 50$ km). $A_1$ values are obtained by fitting the peak shock pressure in Zone 1 (close to the impact site) from iSALE models. The correlation coefficients ($r$), between the values obtained by the iSALE model and by the scaling law is mentioned in the legend box.</caption>
<figure><img src="image_36.png" /><figcaption>Fig. 6. (continued)</figcaption></figure>

observed transition and a phase change in the Equation of State from our iSALE models. The Zone 1 - Zone 2 boundary occurs at a constant distance close to 10 times the radius of the impactor for the range of impactor sizes studied (Fig. 7, left). In the cases shown in Fig. 7, $\beta$ is close to the critical value of 30° where the ray crosses the CMB (Fig. 2). As R<sub>imp</sub> increases, the minimal value D/R<sub>imp</sub> for which the shockwave crosses the CMB decreases from 34 (for R<sub>imp</sub> = 50 km, black curve) to 8.5 (for R<sub>imp</sub> = 200 km, blue curve). At this singular point (Fig. 7, left),

a small jump of the peak pressure occurs because of the reflection of the shock waves with $\beta$ < 30° at the CMB. The ratio of the isobaric core radius to the projectile radius weakly depends on the impact velocity, as also concluded by Pierazzo et al. [1997]. For the range of impact velocities considered, the Zone 1 - Zone 2 boundary ranges between 7 and 15 times the radius of the impactor (Fig. 7, right). The small peak shock pressure jumps identified in Fig. 7 (left) for large impactors does not occur here because of the small, 50 km, impactor radius.

<figure><img src="image_37.png" /><figcaption>Fig. 7. Shock pressure along the ray angle $\beta = 36^\circ$ as a function of the distance from the impact site for different impactor sizes with an impact velocity of 10 km/s (left) and different impact velocities for an impactor radius of 50 km (right). On the left panel, the vertical dotted line represents $D_{Z1} = 10 R_{imp}$. On the right panel, the four vertical dotted lines represent different values of $D_{Z1}$ corresponding to the different values of $V_{imp}$.</figcaption></figure>
Once the shockwaves enter into Zone 2, the peak pressure continues to decrease. Fig. 8a and b show the dependency of $A_2$ on the ray angle (see Eq. (1b)). Because peak shock pressure decreases more rapidly with distance in the models with strength than in the hydrodynamic models (see Fig. 3), $A_2$ values are smaller in the models with strength than in the hydrodynamic models. Also, $A_2$ only slightly decreases with the ray angle $\beta$ in both cases. For all impact velocities and impactor radii studied here, the $A_2$ values obtained from the models with strength are nearly constant up to the values of $\beta$ close to 65°. For $\beta > 65°$, the data is more scattered (see Discussion at the end of the section). For the large impactor, $R_{imp} = 200$ km, the evolution of $A_2$ as a function of $\beta$ is significantly different. This is probably due to the fact that, in that case, Zone 1 overlaps the CMB.

Fig. 9a and b show the relationship between $A_3$ (see Eq. (1b)) and the ray angle $\beta$. As for $A_1$, $A_3$ is smaller in the models with strength than in the hydrodynamic models. This feature is related to the higher shear strength in the models with strength, which consumes a portion of the impact energy to plastically deform the impacted material. Again, the more moderate decay in Zone 2 (i.e. smaller $A_3$ values) in models with strength compared to hydrodynamic models might be the consequence of a solid-state phase transformation that is considered in the ANEOS for dunite [Collins and Melosh, 2014]. Apart from this difference, both models behave similarly. Fig. 9a and b show that $A_3$ slightly increases

with $\beta$ up to $\beta \approx 65°$ indicating that the peak shock pressure decreases more efficiently with distance at large ray angles. From the results illustrated in Fig. 9b, the $A_3$ coefficient does not significantly depend on the impact velocity. However, Fig. 9a shows that the tendencies described earlier are completely different for $R_{imp} \ge 150$ km. Indeed, for large impactors, the size of Zone 2 is significantly reduced and the impacted hemisphere mainly consists of Zone 1 and the metallic core within the ray angle $< 30°$.

To illustrate the dependencies of $A_1$, $A_2$ and $A_3$ on $\beta$, we express them as $A_1 = A_6 \exp(A_7 \beta)$ where i denotes 1, 2, or 3 (see legend boxes in the corresponding figures). Figs. 6, 8 and 9 show small values for $A_7$, typically ranging between $10^{-2}$ and $10^{-3}$. This confirms that the shock wave propagation weakly depends on the ray angle that is ranging between 0° (vertical) and 45°, as already concluded by previous studies [Pierazzo et al., 1997; Monteux and Arkani-Hamed, 2016]. This can be extended for $\beta$ up to 60° on the basis of Figs. 6, 8 and 9. For the hydrodynamic models the trend is similar although the iSALE data are more dispersed. Figs. 8 and 9 also show that for $R_{imp}$ values larger than 150 km, the trend is less clear than for small impactors. As already mentioned, we attribute this feature to the overlap between Zone 1 and the CMB, which leads to the progressive disappearance of Zone 2 as $R_{imp}$ increases. These differences are confirmed by the correlation coefficients that are larger for small impactors and are generally higher for

<figure><img src="image_38.png" /><figcaption>Fig. 8. a: $A_2$ as a function of the ray angle $\beta$ (Eq. (1b)): $P_{Z2} = A_2 (D/R_{imp})_{3}^{A}$ for different impactor radii, but an impact velocity of 10 km/s. $A_2$ values are obtained by fitting the shock pressure in Zone 2 (beyond $D_{Z1}$) from iSALE models. In these figures, $D_{Z1} = 10$. The correlation coefficient (r) for each scaling law is mentioned in the legend box (note that fits are restricted to $\beta \le 65°$).</figcaption></figure>

<caption>b: $A_2$ as a function of the ray angle $\beta$ (Eq. (1b)): $P_{Z2} = A_2 (D/R_{imp})_{3}^{A}$ for different impactor velocities, but an impactor radius of 50 km. $A_2$ values are obtained by fitting the shock pressure in Zone 2 (beyond $D_{Z1}$) from iSALE models. In these figures, $D_{Z1} = 7, 10, 13$ and $15$ for $V_{imp} = 5, 10, 15$ and $20$ km/s respectively. The correlation coefficient (r) for each scaling law is mentioned in the legend box (note that fits are restricted to $\beta \le 65°$).</caption>
<figure><img src="image_39.png" /><figcaption>Fig. 8. (continued)</figcaption></figure>

the models with strength than for the hydrodynamic models where the iSALE data is more dispersed.

Figs. 8 and 9 also exhibit a large variability for β values larger than 65°. For these large β values the shock waves propagate through shallower regions (see Fig. 2) where the total pressure is much smaller than that in the deep mantle. Moreover, in the shallowest region the largest initial impact is followed by multiple smaller impacts caused by ejecta fall back. Finally, the near surface region is also affected by the wave reflected at the surface that interferes with the direct wave and complicates the determination of the peak shock pressure [Melosh, 1989; Arkani-Hamed, 2005; Louzada and Stewart, 2009]. These 3 effects may explain the variability of the A₁, A₂, A₃ constants for large β values.

## 2.5. Effects of the CMB

The distinctly different physical properties of the mantle and the core have two major effects on the shock wave as it passes from the mantle to the core. The shock pressure increases abruptly as seen in Fig. 3, and a shock ray bends, obeying the Snell's law. Here we discuss each of these two effects.

Because the core radius is half the size of the impacted planet, the shock waves encounter the CMB only for β values smaller than 30 degrees (Fig. 2). Fig. 10 shows the ratio of the shock pressure immediately above the CMB, Pₘ, to that immediately below the CMB, P_c, in the impacted hemisphere. Both are calculated along a straight line determined by a ray angle and at the grid points of the hydrocode models

immediately above and below the CMB, respectively. We note that Pₘ and P_c are actually the pressures slightly off the shock ray, hence are approximations to the actual pressures on the ray. A given shock ray bends toward the axis of symmetry as it crosses the CMB and the refraction angle determined through the Snell's law is smaller than the incident angle (see below). The bend is more pronounced at higher ray angles.

To connect the peak shock pressure between the mantle and the core we fit the ratio P_c/Pₘ as a function of the ray angle β with a quadratic expression for the models with strength:

$$P_c/P_m = b_1 + b_2\beta + b_3\beta^2 \quad (2a)$$

and a linear expression for the hydrodynamic models:

$$P_c/P_m = b_4 + b_5\beta \quad (2b)$$

where b₁, b₂, b₃, b₄ and b₅ are constant. The values of these constants for both models and as a function of Rₘ and Vₘ are listed in Table 3. The corresponding results are illustrated in Fig. 10a and b. We emphasize here that the Eqs. (2a) and (2b) used to fit the pressure jump as a function of the ray angle β are chosen empirically according to our results illustrated in Fig. 10 and b.

For all the impact velocities and impactor radii used in our models, the behavior is quite similar, although the shockwave that arrives at the CMB has much higher amplitude in the hydrodynamic case than in the models with strength (Fig. 3). Fig. 10 shows that the pressure jump at the CMB is more pronounced in the models with strength than in the
hydrodynamic models, regardless of the amplitude of the arriving shockwave. Because the pressure jump at the CMB is likely related to the amplitude of the arriving shockwave for a given model, our results emphasize the importance of the impacted material properties on the amplitude of the pressure jump.

For large impactors ($R_{imp} > 150$ km) and for $\beta$ values smaller than 30°, we have shown that Zone 1 directly interferes with the impacted planet's core thus leading to a disappearance of Zone 2. Above this critical value of $\beta > 30°$, even if the direct shockwave does not cross the core (Fig. 2), one can expect that the transition from Zone 1 to Zone 2 would still be identified in our models. However, Fig. 8a (bottom panels) and Fig. 9a (bottom panels) also show a different behavior for large impactors for $\beta > 30°$. We attribute this phenomenon to the interaction of the direct shockwave with the core which produce a variety of diffracted and reflected shockwaves that affect the peak shock pressure far from the impact site and for $\beta > 30°$. As the shockwave behavior is completely different, our scaling laws for $A_2$, $A_3$ and $P_c/P_m$ proposed in Eqs. (1b), (2a) and (2b) are less viable. The tendency for $A_2$ and $A_3$ is difficult to identify from our results (Figs. 8a and 9a (bottom panels)) and this complexity probably illustrate the complexity of the shockwave propagation far from the impact site.

The incident angle of a shock ray at the CMB increases from 0° to 90° as the ray angle $\beta$ increases from 0° to 30°. The refraction angle of a

given shock ray is determined on the basis of the Snell's law,

$$U_m/\sin \gamma_m = U_c/\sin \gamma_c \quad (3)$$

where $U_m$ and $U_c$ are the shock wave velocities in the mantle and the core, $\gamma_m$ is the incident angle, and $\gamma_c$ is the refraction angle. $U_m$ and $U_c$ are determined from our numerical models (see the next section).

## 2.6. Shock wave propagation in the core

Fig. 11 shows the dependency of the coefficients $A_4$ and $A_5$ in Eq. (1c) on $V_{imp}$ and $R_{imp}$ determined along the axis of symmetry using the hydrocode models. $A_4$ shows a power-law behavior while $A_5$ shows almost linear dependency, especially those associated with the models with strength.

For a given ray angle $\beta$ we calculate the incident and the refraction angles at CMB, hence the corresponding shock ray path in the core, assuming linear. The peak shock pressure decay along this ray path is determined from Eq. (1c), where $D$ is the distance from the impact site to a point in the core measured along the ray path assuming that $A_4$ and $A_5$ are independent of the ray angle. Fig. 12 shows the distance of the peak shock pressure from the impact site in the impacted hemisphere of the mantle and the entire core versus the time since the impact. Table 4 lists the corresponding shock wave velocity in the mantle and the core.

<figure><img src="image_40.png" /><figcaption>Fig. 9. a: $A_3$ as a function of the ray angle $\beta$ (Eq. (1b): $P_{Z2} = A_2 (D/R_{imp})^{\gamma_2}$) for different impactor radii. $A_3$ values are obtained by fitting the shock pressure in Zone 2 (beyond $D_{Z1}$) from iSALE models. In these figures, $D_{Z1} = 10$. The correlation coefficient ($r$) for each scaling law is mentioned in the legend box (note that fits are restricted to $\beta \le 65°$).</figcaption></figure>

<caption>b: $A_3$ as a function of the ray angle $\beta$ (Eq. (1b): $P_{Z2} = A_2 (D/R_{imp})^{\gamma_2}$) for different impact velocities. $A_3$ values are obtained by fitting the shock pressure in Zone 2 (beyond $D_{Z1}$) from iSALE models. In these figures, $D_{Z1} = 7, 10, 13$ and $15$ for $V_{imp} = 5, 10, 15$ and $20$ km/s respectively. The correlation coefficient ($r$) for each scaling law is mentioned in the legend box (note that fits are restricted to $\beta \le 65°$).</caption>
<figure><img src="image_41.png" /><figcaption>Fig. 9. (continued)</figcaption></figure>

The shock wave velocities vary slightly among the models. More importantly, the velocity in the core is consistently lower than that in the mantle, indicating that the refracted ray angle is always smaller than the incident ray angle.

### 3. Comparison of the scaling law models with the iSALE models

Before implementing our scaling laws to the accretion of a planet it is required to compare the scaling law models with the iSALE models. Here we compare the hydrocode models and the models with strength to the corresponding models determined using our scaling laws for the impact velocity of 10 km/s and the impactor radius of 50 km (Fig. 13 a and b). In the scaling law models the refracted angles of the shock rays at the CMB are determined by the Snell's law. Hence, the shock pressure is not estimated in the equatorial region of the core where no refracted wave propagates. This results in the distinct lack of shock pressure in the equatorial part of the core and in the lower mantle beneath the shock ray that is tangent to the core mantle boundary. Several other features are not modeled by the scaling laws because we do not consider the reflection at the CMB and at the surface or the reflection and refraction at the antipodal CMB, which are included in solving the shock dynamic equation in the hydrocode models.

Fig. 13 (right panel) shows the relative difference between the two models. For the impacted hemisphere, the mean errors are 13% within the mantle and 14% within the core for the model with strength, and 10% within the mantle and 12% within the core for the hydrodynamic model. The fact that iSALE model combines both direct and indirect waves, while our model shows only the direct wave could contribute to

the error between our scaling law parameterization and the iSALE data. The error pattern is similar in both the hydrodynamic and strength models. As mentioned above, the direct shock pressure within the core for $β$ values larger than 30° is not determined by the scaling laws. The shock pressure of the hydrocode model in this region must be due to the secondary shock waves that are reflected at the surface and then refracted at the CMB. The scaling law models overestimate the peak shock pressure close to the impact site (Zone 1) where the error can locally reach 50%. However, despite the large peak shock pressure range in this region, our scaling laws provide a first order approximation to the hydrocode models.

An impact generates two compressional shockwaves, one propagates downward in the planet and the other propagates upward in the penetrating impactor. The latter reaches the top of the impactor, and reflects back and propagates downward as a rarefaction wave. The turbulence induced by the compressional wave in a given location inside the planet is diminished by the trailing rarefaction wave. This limits the thickness of the entire shock front with appreciable turbulence. The velocity of the rarefaction wave is higher than that of the original down going compressional wave, because the rarefaction wave propagates inside the penetrating impactor in the early times and then inside the planet that is already set in motion by the initially down going compressional wave. At a shock wave velocity of about 8 km/s, the passage of the entire shock front at a given location takes less than a minute for the largest impactor of 200 km considered in this study, and much shorter for other impactors. A given location in the mantle experiences more than one shock wave, such as the shock wave reflected at the surface and that reflected at the CMB. However, these shock
waves have low amplitudes because they travel longer distances than the direct wave and they arrive at much later times relative to the direct pressure wave. They rarely interfere with the direct wave except near the surface and near the CMB. In the case of the core, parts of the antipodal core is shocked twice, first by the direct wave that is refracted at the CMB of the impacted hemisphere beneath the impact site and second by the shock wave reflected at the CMB of the antipodal hemisphere. The first wave travels a shorter distance than the second one, hence passes a given location sooner than the arrival of the reflected wave. As in the case of the mantle, there is negligible interaction between these two waves, except for very close to the antipodal CMB. In general, a given location experiences multiple pressure increases due to the multiple reflected and refracted waves that travel through the planet's interior. However, the peak pressure is due to the direct shockwave.

## 4. Impact heating

During the decompression, a fraction of the energy is converted into heat [Bjorkman and Holsapple, 1987]. The temperature increase has been related to the shock pressure [e.g., Watters et al., 2009] as well as to the particle velocity [e.g., Gault and Heitowit, 1963]. The shock pressure in a given location diminishes within a time scale that is much

shorter than the cooling timescale of the impact induced thermal anomaly. Although secondary shock waves lead to local pressure increase, the corresponding temperature increase is negligible compared to that of the direct shock wave. Here we only consider the direct shock heating in our estimation of the impact induced temperature increase.

The impact heating induced by the shockwave propagation is determined by implementing our scaling laws within the theoretical model proposed by Watters et al. [2009]

$$ \Delta T(P_{\delta}) = \left( \frac{P_{\delta}}{2\rho_{0}S} (1 - f^{-1}) - \left(\frac{C}{S}\right)^{2} [f - \ln f - 1] \right) / c_{p} \quad (4) $$

$$ f(P_{\delta}) = -\frac{P_{\delta}}{\beta} \left( 1 - \sqrt{\frac{2P_{\delta}}{\beta} + 1} \right)^{-1} \quad (5) $$

$$ \beta = \frac{C^2 \rho_0}{2S} \quad (6) $$

where $P_{\delta}$ is the shock-induced pressure, $\rho_0$ is the density prior to the shock compression and $c_p$ is the heat capacity of the impacted material (see values in Table 1). The shock pressure implemented in Eqs. (4)–(6), is the shock-induced pressure predicted by our scaling laws (Eqs. (1a)–(1c) and (2a)–(2b)). Fig. 14 shows the temperature increase obtained using our scaling laws for both the hydrodynamic model and the

<figure><img src="image_42.png" /><figcaption>Fig. 10. a: Shock pressure jump $P_c/P_m$ at the core mantle boundary as a function of the ray angle $\beta$ for different impactor sizes, but a single impact velocity of 10 km/s. Data obtained from iSale models are shown with symbols while the analytical expression from Eq. (1c) is illustrated with dashed lines. For values of $\beta$ larger than 30°, the rays do not cross the CMB.</figcaption></figure>

b: Shock pressure jump $P_c/P_m$ at the core mantle boundary as a function of the ray angle $\beta$ for different impact velocities, but a single impactor radius of 50 km. Data obtained from iSale models are shown with symbols while the analytical expression from Eq. (1c) is illustrated with dashed lines. For values of $\beta$ larger than 30°, the rays do not cross the CMB.