core is probably fully formed and because of the potential onset of mantle convection, the thermal consequences of which are beyond the scope of this paper (Behounkova and Choblet, 2009).

For all the $R-R_{Fe}$ range studied here, conditions studied, we find that $Q \gg Q_A$ after impact, which means that core convection can occur (see Eq. 15). For example, once a metallic diapir with $R_{Fe}=300$ km has reached the center of a planet with radius $R=3500$ km, it starts to cool and over the 50 Myr period studied, $Q/Q_A$ decreases from 32.5 to 20. For $R=2000$ km, $Q/Q_A$ decreases from 15 to 10 over the same time interval (in the calculation of $Q_A$ we consider that $k_{Fe}=40$ W m$^{-1}$ K$^{-1}$). Assuming that the thickness of the thermal boundary layer around the protocore scales with the size of the protocore, this corresponds to a temperature difference of $\approx 300$ to 2700 K between the metallic phase and the overlying mantle material (see Fig. 2). However, as discussed in Section 3, this condition alone does not guarantee dynamo action, and so we next examine the additional requirements for a dynamo.

In the special case of the impacted planet being Earth-like from a magnetic perspective, convection within the protocore has to be sufficiently vigorous to supply energy at a rate that compensates for the loss due to ohmic decay. We monitor the ratio of the heat flow to the ohmic dissipation as a function of time. Assuming a conductive heat transfer across the protocore mantle boundary, $k = k_{Si}$ and $\delta \sim R_{Fe}$, from Eq. 13, Eq. 17 and Eq. 22 we obtain

$$ \left(\frac{Q_{cond}}{\Phi}\right) \sim \frac{3\mu k \Delta T_0}{\nu \bar{B}^2} \left(1 + a(\lambda)\frac{\Delta\theta}{\Delta T_0}\right), \qquad (24) $$

which is a conservative bound on whether a dynamo will occur (we consider here that $r_0/R \approx 1$). $a(\lambda)$ and $\Delta\theta/\Delta T_0$ are independent of $R$. However, from Eq. 5, $\Delta T_0 \propto R^2$. Hence, $(Q_{cond}/\Phi)$ varies with $R^2$. Fig. 3 shows $Q/\Phi$ as a function of the planet radius for different iron diapir sizes and different times (1 Myr and 50 Myr), after the diapir reaches the center of the impacted planet. The dashed line is the theoretical prediction from Eq. 24 for a uniform viscosity case (i.e., $a(\lambda) = 15\%$). Numerical results are consistent with Eq. 24 at very early time (just after the metallic diapir has reached the center of the planet). At larger times, the protocore cools down and $Q/\Phi$ decreases. The temperature of the metallic phase after sinking and hence the initial heat flow across the protocore is a sensitive function of the viscosity contrast between the metallic phase and the surrounding mantle (see Section 3.2). The theoretical value of the initial conductive heat flow when the viscosity of

<figure><img src="image_158.png" /><figcaption>Fig. 3. Temporal evolution (symbols from black to white) of the ratio $Q/\Phi$ as a function of the planetary radius. For each planetary radius, we investigated several diapir sizes at times $t=1$ and 50 Myr for a uniform viscosity. The length of the rectangles represents the variation of $Q/\Phi$ for the range of diapir sizes. The width of each rectangle is arbitrary. Theoretical predictions for $(Q_{cond}/\Phi)$ at $t=0$ Myr from Eq. 24 for uniform viscosity (black dashed line, $a(\lambda)=15\%$) and temperature dependent viscosity (black dotted line, $a(\lambda)=0\%$) are shown. Planetary and moon radii of small solar system objects are indicated on the top x-axis.</figcaption></figure>

the hot metallic fraction is much smaller than the viscosity of the relatively cold undifferentiated material is given by Eq. 13 and Eq. 22 with $a(\lambda)=0\%$. The black dotted line in Fig. 3 shows the corresponding theoretical prediction for $(Q_{cond}/\Phi)$ immediately after sinking as a function of $R$ for $a(\lambda)=0\%$. In this case, the theoretical value for $(Q_{cond}/\Phi)$ is different by a factor of 2.8 compared with the uniform viscosity case, consistent with expectations from Rybczynski–Hadamard effects.

In all sets of planetary and metallic diapir radii studied here $Q \gg Q_A$. Consequently, the condition for the presence of a dynamo from the combination of Eq. 16 and Eq. 18 simplifies to:

$$ \frac{Q}{\Phi} > \frac{3 C_{p,Fe}}{0.8\pi G \alpha_{Fe} \rho_{Fe} R_{Fe}^2} \qquad (25) $$

Fig. 4 shows $Q/\Phi$ as a function of the protocore radius resulting from a single impact for different planet sizes in uniform viscosity cases. For $t=1$ Myr after the end of the sinking (black filled symbols), and for a given planetary radius, different trends appear as a function of $R_{Fe}$. The ratio $Q/\Phi$ is initially larger for larger bodies because $Q \propto \Delta T_0 \propto R^2$. However, as mentioned in Section 2.2 when $R_{Fe} < 70$ km, no viscous heating occurs during the sinking (Monteux et al., 2009) as cooling by diffusion overwhelms viscous heating. As $R_{Fe}$ increases from 70 km radius, viscous heating becomes more and more efficient and leads to an increase of the initial heat flux and as a consequence of $Q/\Phi$. When $R_{Fe}$ reaches a critical value ($\approx 300$ km), the simplification $r_0/R=1$ is less and less valid since the diapir is initially buried deeper. Hence, for a given $R$, as $R_{Fe}$ increases $r_0/R$ decreases so the initial temperature of the protocore, the heat flux and $Q/\Phi$ also decrease. With subsequent cooling, the reduction in core temperature over a given time period $\Delta t$, say, scales as $Q\Delta t / (\bar{\rho}C_pR_{Fe}^3)$. Because heat flow will decline with the core–mantle temperature difference, and the incremental change in core temperature $\propto (1/R_{Fe}^3)$, small protocores cool faster than big ones. Hence, for a given $R$, and for $t=50$ Myr after the end of the sinking (white filled symbols), $Q/\Phi$ increases as a function of $R_{Fe}$.

In Fig. 4, the black dashed line represents the condition required to get a dynamo from Eq. 25. Results show that a dynamo is difficult to initiate with a small protocore (i.e. small impactors). Once the protocore

<figure><img src="image_159.png" /><figcaption>Fig. 4. $Q/\Phi$ 1 Myr after the diapir reaches the center of the impacted planet as a function of the metallic diapir radius (black symbols). The protocore radius differentiated after a single impact corresponds to an impactor radius indicated on the top x-axis. The symbol shape is related to the radius of the impacted planet. The black dashed line represents the threshold expression from Eq. 25 in the isoviscous case ($a(\lambda)=15\%$). Above this dashed line (white domain) a dynamo is plausible, while below it no dynamo is generated. The dotted line corresponds to the same threshold but assuming a protocore less viscous than the surrounding mantle ($a(\lambda)=0\%$). If $Q/\Phi$ after 50 Myr is within the dynamo domain, a second white filled symbol represents $Q/\Phi$ at this time.</figcaption></figure>
radius is larger than 400 km (i.e. $R_{\text{imp}} > 500$ km), however, an early dynamo can be initiated in the first 50 Myr in a planet with a radius greater than 2500 km and in the absence of a molten layer surrounding the protocore. As the planet size increases, the initial protocore temperature increases and the duration of the early dynamo is longer (over 50 Myr). As the protocore size increases, its characteristic cooling time increases and the duration of the early dynamo is also longer. However, a large viscosity contrast between hot iron and cold undifferentiated material will decrease the initial heat flow by a factor of 2.8 (see Fig. 3), and increase the dynamo threshold from Eq. 25. The dotted line in Fig. 4 illustrates this theoretical effect. Fig. 4 shows that a dynamo can still be initiated after an impact, but an increase of the impactor and impacted planet radii are needed relative to the isoviscous case. Hence an impact with $R_{\text{imp}} > 625$ km (i.e. $R_{\text{Fe}} > 500$ km) on a Mars size body (i.e. $R \sim 3500 \text{km}$) can initiate an early magnetic field. The duration of this early dynamo increases with the size of the impactor.

Dynamo generation requires a magnetic Reynolds number larger than a critical value (see Eq. 21). We monitored the mean heat flux $q$ as a function of time and calculated the corresponding magnetic Reynolds numbers. Fig. 5 shows the results for the two extreme planetary radii (i.e. $R=2000$ km and $R=4000$ km) and the range of protocore sizes studied here. The two lines (solid and dashed) correspond to $q = q_A = Q_A/(4\pi R^2)$ with $R=2000$ km and $R=4000$ km respectively. Three regimes emerge. Below, the solid line $q < q_A$ and no convection occurs. Above the solid line and below $Re_m^{\text{crit}}$, convection occurs but no dynamo is generated. Above these criteria, the thermally-driven convective flow within the protocore can generate a dynamo. As mentioned above, the magnetic Reynolds number will be a factor of about 1.4 too large in the isoviscous compared with the temperature-dependent viscosity case. The dotted line represents the values for $Re_m$ expected from the temperature-dependent viscosity model (i.e. $Re_m$ from isoviscous models divided by a factor 1.4). This figure illustrates three results. First, temperature-dependent viscosity has a weak influence on the $Re_m$ values obtained in this study. Second, as the heat flux out of the protocore obtained from our numerical models is ~20 larger than the adiabatic heat flux $q_A$ (see Eq. 15 and discussion above), then $Re_m(q) \gg Re_m(q_A)$. Third, $Re_m$ easily exceeds

the critical value needed to generate a dynamo in the whole range of planet and protocore radii studied here.

## 5. Discussion

These results provide an important link between the accretion and differentiation histories of the terrestrial planets and their potential for early magnetic fields. They suggest that the presence or absence of an early magnetic field on a planet could be related to its accretionary history.

### 5.1. Application to Solar System objects

Our results indicate that for bodies in the range of $R=2000–4000$ km and protocore radii in the range of 300–650 km, the core heat flow exceeds the adiabatic heat flow, and $Re_m$ easily exceeds $Re_m^{\text{crit}}$ for 50 Myr after the iron diapir reaches the center of the planet and in the absence of additional mantle melting at the protocore–mantle boundary. These conditions suggest that for the range of $R-R_{\text{Fe}}$ conditions studied, dynamo action is possible. A more restrictive criterion is that the heat flow must also exceed the ohmic dissipation under Earth-like conditions. In this case, and assuming no mantle melting, we find that a large impact (a 625 km radius impactor, or $R_{\text{Fe}}=500$ km) onto a Mars-size (i.e., $R=3500$ km or greater) undifferentiated planet is needed to initiate a dynamo that would persist at least for ~50 Myr (Fig. 4).

As noted earlier the ohmic dissipation criterion has several underlying assumptions about the dynamo action. In particular, we note that the Earth-like value of the magnetic field in the core $\bar{B}$ may not be applicable to other bodies. For example, Mercury, the only other inner solar system body to have a present-day dynamo, has a magnetic field at its surface that is 100 times weaker than Earth's surface field. While the origin of this weak field is not well understood, it is at least reasonable to propose that Mercury's core field could be weaker than Earth's. A factor of two decrease in the core field results in a factor of four increase in $Q/\Phi$, allowing dynamo action for smaller impactors and impacted bodies.

Independent of its mode of formation, such a scenario is more difficult to envision for the Moon for two reasons. First, the Moon is smaller ($R=1700$ km), and in our models, the heat flow across the core mantle boundary does not reach the critical value required to overcome ohmic dissipation. Second, its volumetric metallic fraction is only 1–2% meaning that the impactor size needed to differentiate a large volume of metal is unrealistically large. However, our models do not consider any contribution of the central body to the impact velocity. $v_{\text{esc}}$ for the Earth is 4–5 times bigger than $v_{\text{esc}}$ for the Moon. Considering that the impact velocity on the Moon is equal to the escape velocity on the Earth can increase the initial temperature of the Moon's protocore and the protocore–mantle boundary heat flow by a factor $\approx 20$. Moreover, our choice of $\gamma$ provides a lower bound on $\Delta T_0$. If $V_\infty$ is not zero, heating within and outside the core will increase. Finally, radiogenic heating (Nimmo, 2009; Williams and Nimmo, 2004), latent heat release and compositional convection within the core (Labrosse et al., 1997) can also lead to a significant increase of the heat flow across the core–mantle boundary. The solidification of a 100 km radius inner core can release $\sim 10^{10}$ W which is 1000 times larger than $\Phi$ for a 200 km core radius planet. Under the assumption that a heat flow high enough to overcome Earth-like ohmic dissipation is required, these sources of heat and/or compositional buoyancy could explain an early dynamo on the Moon or on small radii objects.

### 5.2. The effect of a molten lower mantle layer on dynamo duration

If the rocks at the protocore–mantle boundary are at sufficiently high temperature and pressure to be close to their melting temperature, heat

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

**Fig. 5.** Magnetic Reynolds number of the protocore as a function of its size from our numerical models. The protocore radius differentiated after a single impact corresponds to an impactor radius indicated on the top x-axis. The symbol shape is related to the radius of the impacted planet and its color (black or white) to the time after the diapir reaches the center of the impacted planet. The dotted line represents the numerical results expected from a temperature-dependent viscosity model (i.e., $Re_m$ from isoviscous models divided by a factor 1.4). The theoretical Reynolds number values obtained from Eq. 21 are represented with a solid line (with $\Delta T_0$ corresponding to $R=4000$ km) and a dashed line (with $\Delta T_0$ corresponding to $R=2000$ km). The gray domain illustrates $Re_m^{\text{crit}}$ above which a dynamo is plausible.
transfer from the protocore can cause the progressive growth of a mantle melt layer (see Eq. 5 and Fig. 2). Assuming that the undifferentiated mantle surrounding the protocore is homogenous, that gravity decreases linearly with depth and that $(R_{Fe}/R)^2 \ll 1$, the pressure at the core mantle boundary is:

$$P_{CMB} \approx \frac{2}{3} \pi G \rho_0^2 R^2 \qquad (26)$$

For the impacted planet radii range in our models, $P_{CMB}$ ranges from 10 to 50 GPa. The formation of a molten layer then depends on the initial temperature of the protocore which is governed by the impacted planet radius R (see Eq. 5). Melting occurs when the initial protocore temperature is larger than the solidus of the surrounding material. In the isoviscous case, Fig. 2 shows that the protocore temperature overcome the solidus temperature of the surrounding mantle only when R>3500 km with a chondritic composition (Agee, 1997) or when R>4500 km with a dry peridotite composition (Zerr and Boehler, 1994). When the viscosity of the diapir is smaller than the surrounding undifferentiated solid mantle, the initial temperature of the protocore is smaller and it is more difficult to form a molten layer unless the pre-impact temperature $T_0$ at the center is increased by ≈500 K. As mentioned in Section 2.1, the pre-impact temperature is a function of the accretion rate and of the early radiogenic heating. In our models the pre-impact thermal state is homogenous and $T_0$ = 1250 K. Hence, in our single impact scenario, the formation of a molten layer at the top of the protocore without additional silicate differentiation would require a pre-impact temperature that increases with depth.

The conductive heat flux applied to Eq. 25 is not appropriate if heat transfer from the core melts the lower mantle. Vigorous convection within a growing melt layer will enhance the heat transfer from the protocore. Thus, mantle melting can increase the likelihood of dynamo action over the conditions identified in Fig. 4 but, in return, can decrease the duration of an early dynamo as heat will be removed more quickly. Hence, a key issue is to know whether convection in the melt layer will cause the core to freeze. In the case of a molten layer, the duration of a dynamo is limited by the time needed to remove the excess heat (the specific heat plus the latent heat of crystallization)

<figure><img src="image_161.png" /><figcaption>Fig. 6. Time to remove the excess heat from the protocore for 3 different protocore-mantle boundary heat transfers as a function of the planetary radius. The dotted line represents the case of a conductive boundary layer (Eq. 22, $\delta = R_{Fe}$ and Eq. 27) and the dashed line represents the case of a time dependent conductive boundary layer (Eq. 22, $\delta = 2\sqrt{kt}$ and Eq. 27). The black line shows the case where a molten layer surrounds the protocore (Eq. 23, $\eta_{melt} = 100$ Pa.s and Eq. 27). This time also depends on the size of the protocore. For each case, 2 lines delimit a domain from $R_{Fe} = 150$ km (lower line) to 500 km (upper line).</figcaption></figure>

from the protocore with a heat flow Q. In that case, the duration of the dynamo is:

$$t_{\text{dynamo}} \approx \frac{4\pi R_{\text{Fe}}^3 \rho_{\text{Fe}} [C_{p,\text{Fe}}(T_{\text{core}} - T_{\text{melt,Fe}}) + L_{\text{Fe}}]}{3Q} \qquad (27)$$

where $T_{melt,Fe}$ and $L_{Fe}$ are the melting temperature and the latent heat of the metallic phase (see Table 1 for typical values). Fig. 6 shows the time needed to remove the heat excess as a function of the planetary radius for 3 different heat transfer types: solid state conduction, time-dependent conduction and convection within a molten layer. Our results show that this time can vary with several orders of magnitude. In the case of a conductive layer, the time to remove the heat excess is very long (≈10⁹ yr) and the duration of an early dynamo is essentially governed by the time during which vigorous convection occurs within the protocore (≈10⁶ – 10⁷ yr). However, in the case where a molten convecting layer forms above the protocore, core freezing occurs in 10³ – 10⁴ yr, potentially limiting the duration of the dynamo to <10 ky. As underlined in Fig. 2, the extent of melting depends on the composition of the surrounding rocks and on the easiness to cross the solidus. Hence, the composition of the early mantle is also a key parameter in the duration of the early dynamo.

### 5.3. Future investigations

This dynamo initiation model is a first step toward a more general model of core formation. From the results of our models, several questions need to be addressed:

1. How do hyper-velocity impacts affect the heat repartition after an impact and how does this in turn enhance an impact-induced early dynamo? The fraction of kinetic energy ($\gamma$) that is retained as heat below the impact site and the spatial heat repartition after an impact ($h_m$) are key parameters to characterize the early inner thermal state of growing planets. Better constraints on these parameters would greatly assist any models that investigate the thermal consequences of impacts early in a planet's history.

2. How will multiple diapirs from multiple impacts affect the early protocore thermal state? The sinking of multiple metallic diapirs will strongly affect the thermal state of both the core and the mantle and may play a key role on the temperature contrast between these two reservoirs (Golabek et al., 2009). This prospect is currently beyond the scope of this study because of the geometry used in our numerical models (axisymmetric spherical). However, our study gives a parametrization of the inner thermal consequences after a single large metallic diapir sinking. The scalings proposed here could be used to compute more accurate one dimensional thermal evolution models of growing planets.

3. What is the influence of the timescale between two successive impacts? If a second impact occurs shortly after the first (i.e., within, say, three thermal diffusion timescales $RFe^2/k$), then the thermal effect of the first impact will still be present within the mantle and the second diapir will sink into a warmer and less viscous mantle. The second diapir will descend more quickly. Moreover, the viscosity contrast between the second diapir and the hotter mantle will decrease and viscous heating will be relatively enhanced within the second metallic diapir (see Section 2.2).

4. In the framework of a late giant impact between two differentiated bodies, what is the effect of the core merging on the dynamo activity? The thermal effect of giant impact on dynamos has been recently studied (Reese and Solomatov, 2010; Roberts et al., 2009) but the dynamics of merging between the preexisting core and the impactor's core and the thermal consequences have also to be considered in future models. Within the framework of our model, a giant impact might kill a preexisting dynamo depending on size of the impactor as the thermal effect of both the impact and the
viscous heating during the sinking may reduce the total heat flow across the CMB.

5. How would a second large impact or merging phenomenon change the Hf-W systematics? Depending on the efficiency of equilibration between the impactor's core and the impacted body's mantle, this merging phenomenon may be a major feature of the current W isotopes composition of planets and influence the Hf-W chronology of accretion (Jacobsen, 2005).

**6. Conclusions**

The dynamics of terrestrial planetary accretion governs their differentiation histories and their potential for early magnetic fields. We have shown that under Earth-like conditions for ohmic dissipation in the core impactors (with $R_{imp} > 625$ km) leading to major differentiation events during early stages of planetary formations can initiate early dynamos. If the ohmic dissipation criterion is relaxed, allowing dynamo action simply if the critical magnetic Reynolds number in the core is exceeded, we find that dynamo action is possible for a larger range of impactors and planets (notably smaller impactors and smaller impacted bodies). Such early dynamos can persist for several kyr to several Myr depending on the size of both the impactor and the impacted body and depending on the heat transfer regime across the protocore. The duration of the early dynamo is also a strong function of the rheology of the early mantle of the impacted body. This rheology determines the sinking velocity of the metallic diapir, the amount of viscous heating during the sinking and the ease to remove the heat from the protocore in the case of a molten layer.

**Acknowledgments**

The authors thank Laboratoire de Sciences de la Terre de Lyon for its computing time and R. Pawlowicz for useful discussions. We also thank three anonymous reviewers for thoughtful and constructive comments. This project was funded by a Lavoisier fellowship and the Canadian Institute for Advanced Research. A.M. Jellinek and C.L. Johnson acknowledge support from NSERC and C.L. Johnson from NASA's Planetary Geology and Geophysics Program.

**References**

Acuña, M., et al., 1999. Global distribution of crustal magnetization discovered by the mars global surveyor MAG/ER experiment. Science 284, 790–793.

Agee, C.B., 1997. Melting temperatures of the Allende meteorite: implications for a Hadean magma ocean. Phys. Earth Planet. Int. 100, 41–47.

Behouinkova, M., Choblet, G., 2009. Onset of convection in a basally heated spherical shell, application to planets. Phys. Earth Planet. Int. 176, 157–173.

Breuer, D., Spohn, T., 2003. Early plate tectonics versus single-plate tectonics on Mars: evidence from magnetic field history and crust evolution. J. Geophys. Res. (Planets) 108, 5072–5085.

Buffett, B.A., 2002. Estimates of heat flow in the deep mantle based on the power requirements for the geodynamo. Geophys. Res. Lett. 29 (12), 1–4.

Christensen, U.R., 2010. Dynamo scaling laws and applications to the planets. Space Sci. Rev. 152, 565–590.

Christensen, U.R., Aubert, J., 2006. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields. Geophys. J. Int. 166, 97–114.

Christensen, U.R., Wicht, J., 2008. Models of magnetic field generation in partly stable planetary cores: applications to Mercury and Saturn. Icarus 196, 16–34.

Cisowski, S.M., Collinson, D.W., Runcorn, S.K., Stephenson, A., Fuller, M., 1983. A review of lunar paleointensity data and implications for the origin of lunar magnetism. J. Geophys. Res. 88, 691–704.

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

Douglas, J., 1955. On the numerical integration of $\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = \frac{\partial u}{\partial t}$ by implicit methods. J. Soc. Ind. Appl. Math. 3, 42–65.

Elkins-Tanton, L.T., Zaranek, S.E., Parmentier, E.M., Hess, P.C., 2005. Early magnetic field and magmatic activity on Mars from magma ocean cumulate overturn. Earth Planet Sci. Lett. 236, 1–12.

Fei, Y., Bertka, C.M., Finger, L.W., 1997. High-pressure iron-sulfur compound, Fe3S2, and melting relations in the FeFeS system. Science 275, 1621–1623.

Garrick-Bethell, I., Weiss, B.P., Shuster, D.L., Buz, J., 2009. Early lunar magnetism. Science 323, 356–359.

Golabek, G.J., Schmeling, H., Tackley, P.J., 2008. Earth's core formation aided by flow channelling instabilities induced by iron diapirs. Earth Planet Sci. Lett. 271, 24–33.

Golabek, G.J., Gerya, T.V., Kaus, B.J.P., Ziethe, R., Tackley, P.J., 2009. Rheological controls on the terrestrial core formation mechanism. Geochem. Geophys. Geosyst. 10, Q11007. doi:10.1029/2009GC002552.

Hadamard, J., 1911. Mouvement permanent lent d'une sphère liquide et visqueuse dans un liquide visqueux. C. R. Acad. Sci. 152, 1735–1738.

Honda, R., Mizutani, H., Yamamoto, T., 1993. Numerical simulation of Earth's core formation. J. Geophys. Res. 98, 2075–2090.

Hood, L.L., Richmond, N.C., Pierazzo, E., Rochette, P., 2003. Distribution of crustal magnetic fields on Mars: shock effects of basin-forming impacts. Geophys. Res. Lett. 30 (6), 1–4.

Jacobsen, S.B., 2005. The Hf-W isotopic system and the origin of the Earth and Moon. Ann. Rev. of Earth and Planet. Sci. 33, 531–570.

Johnson, C.L., Phillips, R.J., 2005. Evolution of the Tharsis region of Mars: insights from magnetic field observations. Earth and Planet Sci. Lett. 230, 241–254.

Kaula, W.M., 1979. Thermal evolution of earth and moon growing by planetesimal impacts. J. Geophys. Res. 84, 999–1008.

Ke, Y., Solomatov, V.S., 2009. Coupled core-mantle thermal evolution of early Mars. J. Geophys. Res. (Planets) 114 (13), 1–12.

Kleine, T., Münker, C., Mezger, K., Palme, H., 2002. Rapid accretion and early core formation on asteroids and the terrestrial planets from Hf-W chronometry. Nature 418, 952–955.

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

Kuang, W., Bloxham, J., 1997. An Earth-like numerical dynamo model. Nature 389, 371–374.

Labrosse, S., Poirier, J., Le Mouel, J., 1997. On cooling of the earth's core. Phys. of the Earth and Planet. Int. 99, 1–17.

Laney, C.B., 1998. Computational Gasdynamics. Cambridge University Press, Cambridge.

Lawrence, K., Johnson, C., Tauxe, L., Gee, J., 2008. Lunar paleointensity measurements: implications for lunar magnetic evolution. Phys. Earth Planet. Int. 168, 71–87.

Lillis, R.J., Frey, H.V., Manga, M., 2008. Rapid decrease in Martian crustal magnetization in the Noachian era: implications for the dynamo and climate of early Mars. Geoph. Res. Lett. 35, 14,203–14,209.

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

Monteux, J., Ricard, Y., Coltice, N., Dubuffet, F., Ulvrova, M., 2009. A model of metal-silicate separation on growing planets. Earth and Planet Sci. Lett. 287, 353–362.

Nimmo, F., 2009. Energetics of asteroid dynamos and the role of compositional convection. Geoph. Res. Lett. 36, 10,201–10,207.

Nimmo, F., Stevenson, D.J., 2000. Influence of early plate tectonics on the thermal evolution and magnetic field of Mars. J. Geophys. Res. 105, 11,969–11,980.

O'Keefe, J.D., Ahrens, T.J., 1977. In: Merrill, R.B. (Ed.), Impact-induced Energy Partitioning, Melting, and Vaporization on Terrestrial Planets: Lun. Planet. Sci. Conf., 8, pp. 3357–3374.

Peaceman, D.W., Rachford, H.H., 1955. The numerical solution of parabolic and elliptic differential equations. J. Soc. Ind. Appl. Math. 3, 28–41.

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

Ratcliff, J.T., Tackley, P.J., Schubert, G., Zebib, A., 1997. Transitions in thermal convection with strongly variable viscosity. Phys. Earth Planet. Int. 102, 201–212.

Reese, C.C., Solomatov, V.S., 2006. Fluid dynamics of local martian magma oceans. Icarus 184, 102–120.

Reese, C.C., Solomatov, V.S., 2010. Early martian dynamo generation due to giant impacts. Icarus 207, 82–97.

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., Lillis, R.J., Manga, M., 2009. Giant impacts on early Mars and the cessation of the martian dynamo. J. Geophys. Res. (Planets) 114 (13), 4009–4019.

Roe, P.L., 1986. Characteristic-based schemes for the Euler equations. Annu. Rev. Fluid Mech. 18, 337–365.

Rubie, D.C., Melosh, H.J., Reid, J.E., Liebske, C., Righter, K., 2003. Mechanisms of metal-silicate equilibration in the terrestrial magma ocean. Earth Planet. Sci. Lett. 205, 239–255.

Rybczynski, W., 1911. über die fortschreitende bewegung einer flüssigen kugel in einen medium. Bull. Acad. Sci. Cracowie 1, 40–46.

Safronov, V.S., 1978. The heating of the earth during its formation. Icarus 33, 3–12. doi:10.1016/0019-1035(78)90019-2.

Samuel, H., Tackley, P.J., 2008. Dynamics of core formation and equilibration by negative diapirism. Geochem. Geophys. Geosyst. 9, 6011–6026.

Samuel, H., Tackley, P.J., Evonuk, M., 2010. Heat partitioning in terrestrial planets during core formation by negative diapirism. Earth and Planet Sci. Lett. 290, 13–19.

Schubert, G., Russell, C.T., Moore, W.B., 2000. Geophysics: timing of the martian dynamo. Nature 408, 666–667.

Schubert, G., Turcotte, D.L., Olson, P., 2001. Mantle Convection in the Earth and Planets. Cambridge University Press.

Senshu, H., Kuramoto, K., Matsui, T., 2002. Thermal evolution of a growing Mars. J. Geophys. Res. 107, 1–13.

Shannon, M.C., Agee, C.B., 1996. High pressure constraints on percolative core formation. Geophys. Res. Lett. 23, 2717–2720.

Solomatov, V.S., 2000. Fluid Dynamics of a Terrestrial Magma Ocean. In: Canup, R.M., Righter, K., 69 collaborating authors (Eds.), Origin of the earth and moon. University of Arizona Press, Tucson, pp. 323–338.
Srámek, O., Ricard, Y., Dubuffet, F., 2010. A multiphase model of core formation. G. J. Int. 181, 198–220.

Stegman, D.R., Jellinek, A.M., Zatman, S.A., Baumgardner, J.R., Richards, M.A., 2003. An early lunar core dynamo driven by thermochemical mantle convection. Nature 421, 143–146.

Stevenson, D.J., 1990. Fluid dynamics of core formation. In: Newsom, H.E., Jones, J.H. (Eds.), Origin of the Earth. Oxford Univ, New York. 231–249 pp.

Stevenson, D.J., 2003a. Planetary magnetic fields. Earth Planet Sci. Lett. 208, 1–11.

Stevenson, D.J., 2003b. Planetary science: mission to Earth's core – a modest proposal. Nature 423, 239–240.

Stevenson, D.J., Spohn, T., Schubert, G., 1983. Magnetism and thermal evolution of the terrestrial planets. Icarus 54, 466–489.

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

Touboul, M., Kleine, T., Bourdon, B., Palme, H., Wieler, R., 2007. Late formation and prolonged differentiation of the Moon inferred from W isotopes in lunar metals. Nature 450, 1206–1209.

Turner, J.S., 1973. Buoyancy Effects in Fluids. Cambridge University Press, Cambridge.

Weiss, B.P., Berdahl, J.S., Elkins-Tanton, L., Stanley, S., Lima, E.A., Carporzen, L., 2008. Magnetism on the angrite parent body and the early differentiation of planetesimals. Science 322, 713–716.

Williams, J., Nimmo, F., 2004. Thermal evolution of the Martian core: implications for an early dynamo. Geol. 32.

Wood, B.J., Walter, M.J., Wade, J., 2006. Accretion of the Earth and segregation of its core. Nature 441, 825–833.

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

Yoshino, T., Walter, M.J., Katsura, T., 2003. Core formation in planetesimals triggered by permeable flow. Nature 422, 154–157.

Zerr, A., Boehler, R., 1994. Constraints on the melting temperature of the lower mantle from high-pressure experiments on MgO and magnetite. Nature 371, 506–508. doi:10.1038/371506a0.

Ziethe, R., Spohn, T., 2007. Two-dimensional stokes flow around a heated cylinder: a possible application for diapirs in the mantle. J. Geophys. Res. 112, 1–13.
# Can a sinking metallic diapir generate a dynamo?

Julien Monteux,<fnref n="1" /> Nathanaël Schaeffer,<fnref n="2" /> Hagay Amit,<fnref n="1" /> and Philippe Cardin<fnref n="2" />

Received 5 March 2012; revised 21 August 2012; accepted 21 August 2012; published 10 October 2012.

[1] Metallic diapirs may have strongly contributed to core formations during the first million years of planetary evolutions. The aim of this study is to determine whether the dynamics induced by the diapir sinking can drive a dynamo and to characterize the required conditions on the size of the diapir, the mantle viscosity and the planetary latitude at which the diapir sinks. We impose a classical Hadamard flow solution for the motion at the interface between a spherical sinking diapir and a viscous mantle on dynamical simulations that account for rotational and inertial effects in order to model the flow within the diapir. The flows are confined to a velocity layer with a thickness that decreases with increasing rotation rate. These 3D flows are then used as input for kinematic dynamo simulations to determine the critical magnetic Reynolds number for dynamo onset. Our results demonstrate that the flow pattern inside a diapir sinking into a rotating planet can generate a magnetic field. Large diapirs (R > 10 km) sinking in a mantle with a viscosity ranging from 10⁹ to 10¹⁴ Pa.s provide plausible conditions for a dynamo. Equatorial sinking diapirs are confined to a thicker velocity layer and are thus possibly more favorable for dynamo generation than polar sinking diapirs. In addition equatorial sinking diapirs produce stronger saturated magnetic fields. However, for the range of parameters studied here, estimation of the intensity of diapir-driven magnetic fields suggests that they could not have contributed to the lunar or Martian crustal paleomagnetic fields.

**Citation:** Monteux, J., N. Schaeffer, H. Amit, and P. Cardin (2012), Can a sinking metallic diapir generate a dynamo?, *J. Geophys. Res.*, 117, E10005, doi:10.1029/2012JE004075.

## 1. Introduction

[2] Episodic migration of large volumes of iron may have strongly contributed to core formations during early planetary evolution [Stevenson, 1981]. Complete separation of the metallic phase from the silicate chondritic material may occur within only several million years on terrestrial planets [Kleine et al., 2002; Yin et al., 2002; Touboul et al., 2007]. Such a fast process necessarily involves melting of the metallic phase which requires a large amount of energy. This energy can be provided by radiogenic heating [Yoshino et al., 2003], impact processes [Tonks and Melosh, 1993] or gravitational energy release during metal-silicate separation [Śrámek et al., 2010]. Once molten, the dense metallic phase can sink toward the center of the planet via various phenomena such as percolation through a solid matrix [Shannon and Agee, 1996], settling of metallic droplets within a convection magma ocean [Rubie et al., 2003; Höink et al., 2006],

two-phase flow dynamics [Ricard et al., 2009; Šrámek et al., 2010] or Rayleigh-Taylor instabilities [Stevenson, 1981; Honda et al., 1993]. These various processes may lead to a wide range of characteristic length scales [Karato and Murthy, 1997]. Percolation and metal rainfall occur at small length scales (~10⁻<fnref n="2" /> m) and can ultimately form a dense metallic layer at the bottom of a magma ocean. Then, depending on the thickness of this accumulated metallic layer, Rayleigh-Taylor instabilities can mobilize up to 100-km sized diapirs [Honda et al., 1993].

[3] Large scale metallic volumes can also be separated from silicates by differentiation events following a large impact on an undifferentiated planet [Tonks and Melosh, 1993; Monteux et al., 2009]. Indeed, after a large impact a fraction of the kinetic energy is transferred at depth to thermal energy via shock waves propagation and may separate locally the metallic phase from the silicate phase. Then, the dense metallic phase overcomes viscous heating during its sinking toward the center of the planet. If the separated metallic volume is hot enough, vigorous convection occurs within the metallic phase once it has reached the center of the planet, and a dynamo can be generated [Monteux et al., 2011a]. Ultimately, large iron cores of planetesimals could even have directly sunk into the pre-existing core just after an impact [Benz et al., 1987; Monteux et al., 2011b]. In this study, we term metallic diapir every volume of iron that sinks toward the center of the growing planet independently of its size or sinking process.

Corresponding author: J. Monteux, Laboratoire de Planétologie et de Géodynamique, UMR 6112, CNRS, Université de Nantes, Nantes Atlantiques Universités, 2 rue de la Houssinière, FR-44000 Nantes, France. (julien.monteux@univ-nantes.fr)

©2012. American Geophysical Union. All Rights Reserved.
0148-0227/12/2012JE004075

<footnote n="1">¹Laboratoire de Planétologie et de Géodynamique, UMR 6112, CNRS, Université de Nantes, Nantes Atlantiques Universités, Nantes, France.</footnote>

<footnote n="2">²ISTerre, Université de Grenoble 1, CNRS, Grenoble, France.</footnote>
[4] Large metallic diapirs lead to a large deviatoric stress that deforms the surrounding material and accommodates the sinking of the diapir toward the center of the planet. Depending on the rheology of the surrounding mantle, negative metallic diapirism may lead to viscous deformation [Monteux et al., 2009], non-linear rheology [Samuel and Tackley, 2008], elasto-plastic deformation [Gerya and Yuen, 2007] or even to fracturing if the ultimate strength of solid rocks is exceeded [Davies, 1982; Stevenson, 2003]. During the early times of planetary formation extinct radioactivities ($^{26}$Al and $^{60}$Fe) and large impact heating have played a major role in the thermal budget of the planets. This early heating was retained in the internal part of the planet, leading to a hot proto-mantle and favoring sinking via viscous deformation. Hence, large metallic diapirs potentially sank toward the center of the planet in a Hadamard-regime [Hadamard, 1911] with a timescale that is a function of the viscosity of the protomantle [Karato and Murthy, 1997]. This Hadamard flow constrains the internal dynamics of the diapir.

[5] Can a sinking metallic diapir generate a dynamo? Dynamo theory requires three main conditions for the generation of magnetic field: a fluid with large electrical conductivity, large velocities and some favorable flow pattern. In large-scale dynamos operating in planetary outer cores, the energy source generating large velocities is typically thermochemical convection driven by secular cooling and light element release due to a freezing inner core [Olson, 2007]. However, non-convective energy sources may also drive dynamos, for example precession [Tilgner, 2005], tidal distortion of the core-mantle boundary (CMB) [Le Bars et al., 2011] or heterogeneous mechanical boundary conditions [Guervilly and Cardin, 2010].

[6] Kinematic dynamos provide vital information about the efficiency of magnetic field generation for different flow morphologies [Gubbins et al., 2000a, 2000b]. Very simple laminar flows often lack significant helicity which is an important ingredient for a dynamo [Moffatt, 1978]. They are therefore inefficient in generating a magnetic field, which requires the magnetic Reynolds number to exceed large critical values [Gailitis, 1970; Dudley and James, 1989; Moss, 2008]. Even if the helicity is not indispensable in the dynamo process, these kinematic dynamo models have shown that it is a favorable factor. Self-consistent 3D numerical dynamo models show a hemispheric anti-symmetric pattern of helicity associated with the $\alpha^2$ dynamo mechanism [Olson et al., 1999]. In these models helicity is generated by the action of the Coriolis force that yields positive/negative radial vorticity correlated with downwelling in the northern/southern hemisphere respectively [Olson et al., 2002; Amit et al., 2007].

[7] We use kinematic dynamo simulations to determine the critical magnetic Reynolds number above which the flow strength (for a given flow pattern) is sufficient to amplify a magnetic seed field. To resolve the flow inside the diapir, we impose the simple flow pattern of a viscous bubble as a mechanical outer boundary condition on a dynamical simulation. Because the diapir is in a rotating planetary reference frame and its low viscosity may accommodate inertial effects, we take into account Coriolis and inertial forces. These rotational and inertial effects may help dynamo action by introducing helicity to the dynamics within the interior of the sinking diapir.

[8] The aim of this study is to determine whether such sinking dynamics can drive a dynamo, and to constrain the required conditions for the size of the diapir and the mantle viscosity of the planet. In addition, we also investigate the importance of the planetary latitude at which the diapir sinks on the conditions for obtaining a dynamo. In section 2, we show that there exists a considerable range of parameters in which large metallic diapirs are stable. The theory and method are described in section 3. We present the results of our models in section 4. In section 5 we discuss the diapir driven magnetic field growth rate and intensity from an initial interplanetary magnetic seed field. In section 6 we discuss the feasibility of obtaining diapir driven dynamos in geophysical conditions. Conclusions and possible planetary implications are highlighted in section 7.

## 2. Stability of Large Metallic Diapirs

[9] In the laminar flow regime, the sinking velocity of a metallic diapir obeys a Hadamard velocity [Hadamard, 1911; Batchelor, 1967; Monteux et al., 2009]

$$U = \frac{2}{9} \left( \frac{\eta + \eta_s}{\eta + \frac{2}{3}\eta_s} \right) \frac{\Delta\rho g R^2}{\eta_s} \approx \frac{1}{3} \frac{\Delta\rho g R^2}{\eta_s} \quad (1)$$

where $\eta$ is the viscosity inside the diapir, $\eta_s$ is the viscosity of the surrounding mantle (with $\eta_s \gg \eta$), $\Delta\rho$ is the density difference between the metallic diapir and the silicate mantle, $g$ is the gravitational acceleration at the depth of the diapir and $R$ is the radius of the metallic diapir. For early partially molten planetary mantles, $\Delta\rho \sim 5000 \text{ kg/m}^3$ and $\eta_s \sim 10^7 - 10^{13} \text{ Pa.s.}$ [Karato and Murthy, 1997]. Actually, the gravity and the sinking velocity decrease with depth. For simplicity however, we fix the gravity in the partially molten zone to a constant value of $g \sim 5 \text{ m/s}^2$ (intermediate value between the Moon's and the Earth's surface gravities). The sinking velocity from equation (1) is related to the Stokes velocity by $U_{Stokes} = 2U/3$ which is the velocity of a solid sphere sinking in a viscous medium (i.e. $\eta_s \ll \eta$).

[10] Whether a large metallic diapir will break depends on the stability of the surrounding flow, which is depicted by the Reynolds number $Re_s$ based on the sinking velocity $U$ and the fluid properties of the surrounding mantle:

$$Re_s = \frac{\rho_s U R}{\eta_s} = \frac{\rho_s \Delta\rho g R^3}{3\eta_s^2} \quad (2)$$

where $\rho_s \sim 5000 \text{ kg/m}^3$ is the estimated density of the mantle. For $Re_s \gtrsim 10^3$, viscous forces are negligible compared to inertia and diapir break-up occurs. When $1 \lesssim Re_s \lesssim 10^3$, the diapir may or may not break depending on the importance of inertia over surface tension forces. For $Re_s \lesssim 1$, inertia is negligible and no break-up occurs [Samuel, 2012]. In Figure 1, we represent the stability domain for large metallic diapirs as a function of diapir size $R$ and mantle viscosity $\eta_s$. The conservative upper bound $Re_s = 10$ plotted in Figure 1 provides a limit on candidate diapirs for dynamo action. Figure 1 shows that even very large diapirs can sink without break-up if the molten mantle viscosity is sufficiently large.

[11] Because the mantle viscosity is strongly dependent on the temperature, the thermal evolution of the protoplanet is a
<figure><img src="image_162.png" /><figcaption>**Figure 1.** Candidate diapirs for possible dynamo action (shaded region) must have $Re_s \approx 10$ to be stable and $Rm \gtrsim 100$ for magnetic induction to be possible. The $Ro = 1$ contour shows that larger $\eta_s$ gives more dominant rotational effects in the fluid dynamics inside the diapir (above this line $Ro > 1$). For definitions of $Re_s$, $Rm$ and $Ro$ see equations (2), (3) and (7).</figcaption></figure>

key parameter to define the plausible metal diapir sinking regimes. In the early stages of planetary formation and before the dissipation of potential energy associated with the full core formation, the growing planets were schematically divided in three domains: a completely molten shallow magma ocean, a partially molten deep magma ocean and a solid core [Karato and Murthy, 1997].

[12] The most superficial layer of a growing planet was probably completely molten forming a several hundred kilometers thick magma ocean. In this magma ocean, the viscosity ranges $10^{-4}-10^2$ Pa.s. In the shallow magma ocean and considering $R > 1$ km, $Re_s \gg 10^5$ and turbulent effects lead to emulsification and to rapid break-up of the metallic diapir until it reaches a cm size [Rubie et al., 2003; Samuel, 2012]. However, Dahl and Stevenson [2010] have shown that diapirs with radii larger than 10 km can survive complete erosion and preserve most of their initial volume as they reach the bottom of a 1000 km thick magma ocean. In the turbulent magma ocean, the metallic material rapidly segregates from silicates and forms a layer above the rheological transition separating the fully molten magma ocean to the partially molten magma ocean [Höink et al., 2005]. From the iron layer cumulated at the bottom of the molten magma ocean, Rayleigh-Taylor instabilities can mobilize large volume of iron by diapirism [Stevenson, 1981; Honda et al., 1993].

[13] From the bottom of the outermost magma ocean, temperature and melting fraction decrease with depth [Safronov, 1978; Kaula, 1979; Senshu et al., 2002] and the viscosity rapidly increases to $10^7-10^{13}$ Pa.s [Karato and Murthy, 1997] (see Table 1). Assuming that $U$ is a Stokes velocity, for a 100 km radius diapir the typical value of $Re_s$ within this partially molten reservoir ranges between $10^9$ and $10^{-3}$, while for a 1 km radius diapir $Re_s$ ranges between $10^3$ and $10^{-9}$. Hence,

in this zone and with a favorable diapir size and viscosity regime, diapirs can potentially settle in a laminar flow regime (see Figure 1). Within the solid deep core, the viscosity was probably larger than $10^{18}$ Pa.s thus also leading to a laminar Stokes flow regime.

[14] Severe heatings associated with large impacts are phenomena that often occur especially at the end of planetary accretion. In a local volume deep below the impact site, the temperature dramatically increases up to the solidus temperature of silicates which results in massive melting. If the molten volume is not differentiated before the impact, a local separation occurs in a nearly spherical volume [Tonks and Melosh, 1992; Monteux et al., 2009]. The impactor's core material can eventually be collected in the impact-produced magma pool where the local $Re_s \gg 1$ [Deguen et al., 2011]. Then, the sinking of this large collected volume is governed by the rheology of the surrounding material and may occur in a laminar Stokes flow regime depending on the viscosity of the mantle. We restrict our study to the laminar flow regime

<table><caption>**Table 1.** Reference Values for Some Diapirs of Radius $R$ Sinking in a Mantle of Viscosity $\eta_s$, Which Are Possible Candidates for Driving a Dynamo</caption><thead><tr><th>R (km)</th><th>η<sub>s</sub> (Pa.s)</th><th>U (m/s)</th><th>Re<sub>s</sub></th><th>Rm</th><th>E</th><th>Ro</th></tr></thead><tbody><tr><td>0.1</td><td>10<sup>7</sup></td><td>8</td><td>0.4</td><td>400</td><td>10<sup>-6</sup></td><td>10<sup>3</sup></td></tr><tr><td>1</td><td>10<sup>7</sup></td><td>800</td><td>400</td><td>4 × 10<sup>5</sup></td><td>10<sup>-8</sup></td><td>10<sup>4</sup></td></tr><tr><td>1</td><td>10<sup>9</sup></td><td>8</td><td>0.04</td><td>4 × 10<sup>3</sup></td><td>10<sup>-8</sup></td><td>100</td></tr><tr><td>10</td><td>10<sup>9</sup></td><td>800</td><td>40</td><td>4 × 10<sup>6</sup></td><td>10<sup>-10</sup></td><td>10<sup>3</sup></td></tr><tr><td>10</td><td>10<sup>11</sup></td><td>8</td><td>≤1</td><td>4 × 10<sup>4</sup></td><td>10<sup>-10</sup></td><td>10</td></tr><tr><td>10</td><td>10<sup>13</sup></td><td>0.08</td><td>≤1</td><td>400</td><td>10<sup>-10</sup></td><td>0.1</td></tr><tr><td>100</td><td>10<sup>11</sup></td><td>800</td><td>4</td><td>4 × 10<sup>7</sup></td><td>10<sup>-12</sup></td><td>100</td></tr><tr><td>100</td><td>10<sup>13</sup></td><td>8</td><td>≤1</td><td>4 × 10<sup>5</sup></td><td>10<sup>-12</sup></td><td>1</td></tr></tbody></table>
<figure><img src="image_163.png" /><figcaption>**Figure 2.** Streamlines of a sinking metallic sphere (red) in an undifferentiated planet (blue). The streamlines are in the moving frame of the sinking diapir. The ratio between the diapir radius to the planetary radius is 0.125. The three images show the flows for three different ratios between the mantle viscosity $\eta_s$ and the diapir viscosity $\eta$. These results were obtained using the thermo-chemical convection code developed by Monteux et al. [2009] in spherical axisymmetric geometry.</figcaption></figure>

where the diapir does not break. Hence, we particularly focus on the dynamics of a diapir sinking in the deep and partially molten magma ocean.

### 3. Theory and Method

[15] The growth or decay of a magnetic seed field $B_i$ depends on the diapir magnetic Reynolds number $Rm$ which represents the ratio of advection of magnetic field by the flow to magnetic diffusion

$$Rm = \frac{UR}{\lambda} = \frac{\Delta \rho g R^3}{3 \eta_s \lambda} \quad (3)$$

with $\lambda \sim 1 \text{ m}^2/\text{s}$ being the magnetic diffusivity of liquid iron inside the diapir [Olson, 2007]. Numerical dynamo simulations of rotating convection in a spherical shell show that $Rm \gtrsim 100$ is a necessary condition for magnetic field generation [Christensen and Aubert, 2006]. Figure 1 shows the limit $Rm = 100$ as a function of diapir radius $R$ and mantle viscosity $\eta_s$. Combining the conditions on $Re_s$ and $Rm$, we find a considerable range of diapirs (with $R \gtrsim 100 \text{ m}$ and the full estimated range of $\eta_s$) that are both stable against break-up as well as sufficiently fast for significant magnetic induction. It is not clear however whether the particular flow pattern inside the diapir is capable of generating a magnetic field. The goal of this study is to verify the feasibility of dynamo action by a diapir sinking flow.

[16] The laminar flow pattern within a sinking diapir in a non-rotating system is well known from theoretical [Hadamard, 1911; Batchelor, 1967] and experimental [Spells, 1952] studies. This circulation pattern consists of a purely poloidal one-roll axisymmetric velocity field [Batchelor, 1967]:

$$U_r(r, \theta) = U(R^2 - r^2)\cos \theta \quad (4)$$

$$U_{\theta}(r, \theta) = -U(R^2 - 2r^2)\sin \theta \quad (5)$$

where $r$ and $\theta$ are the radial and co-latitude spherical coordinates in the reference frame of the diapir. During its sinking, a spherical diapir might deform depending on the rheological properties of both the metallic and the silicate

phases [Monteux et al., 2009; Ulvrová et al., 2011]. For simplicity we assume that the diapir maintains a spherical shape during its sinking. Figure 2 shows the internal circulation obtained using the thermo-chemical convection code developed by Monteux et al. [2009]. The flow pattern is in good agreement with both the theoretical and the experimental results and weakly depends on the viscosity ratio between the sinking metallic diapir and the surrounding mantle.

[17] When the planet is rapidly rotating and the viscosity of the diapir is sufficiently small, the Coriolis and inertial forces must be taken into account. Two non-dimensional numbers characterize the dynamical system: the ratio of viscous to Coriolis forces is measured by the Ekman number

$$E = \frac{\nu}{\Omega R^2} \quad (6)$$

where $\nu = \eta/\rho$ is the kinematic viscosity of the molten iron and $\Omega$ is the rotation rate of the planet. The ratio of inertial to Coriolis forces is characterized by the Rossby number

$$Ro = \frac{U}{\Omega R} \quad (7)$$

where $U$ is the amplitude of the sinking diapir velocity field (equation (1)). The rotation rate $\Omega$ has probably decreased during the early stages of planetary accretion from a couple of hours [Agnor et al., 1999]. For simplicity, we consider in our study a constant $\Omega$ equal to the current Earth value (see Table 2).

[18] We solve numerically the non-dimensional Navier-Stokes equation including Coriolis and inertial forces in a spherical container representing the sinking diapir:

$$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} + 2\mathbf{e}_z \times \mathbf{u} = -\nabla \pi + E \nabla^2 \mathbf{u} \quad (8)$$

where $\mathbf{u}$ is the velocity field (in units of $R\Omega$), $t$ is time (in units of $\Omega^{-1}$), $\mathbf{e}_z$ is the unit vector in the direction of the planetary rotation axis and $\pi$ the reduced pressure. The Hadamard velocity field (equations (4) and (5)) scaled by $Ro$ is imposed at the spherical outer boundary and drives the flow within the diapir. When the hydrodynamic diapir Reynolds number $Re = Ro/E \ll 10^3$ (not to be confused with $Re_s$),

<table><caption>**Table 2.** Estimated Parameters and Non-dimensional Numbers for Present Earth's Outer Core, a Candidate Sinking Diapir, and Input Non-dimensional Numbers Used in Our Models<sup>a</sup></caption><thead><tr><th></th><th>Earth</th><th>Diapir</th><th>Models</th></tr></thead><tbody><tr><td>ν</td><td>10<sup>-6</sup></td><td>10<sup>-6</sup></td><td>-</td></tr><tr><td>λ</td><td>1</td><td>1</td><td>-</td></tr><tr><td>R</td><td>3 × 10<sup>6</sup></td><td>10<sup>2</sup>–10<sup>5</sup></td><td>-</td></tr><tr><td>Ω</td><td>7 × 10<sup>-5</sup></td><td>7 × 10<sup>-5</sup></td><td>-</td></tr><tr><td>U</td><td>5 × 10<sup>-4</sup></td><td>10<sup>-5</sup>–10<sup>3</sup></td><td>-</td></tr><tr><td>E</td><td>10<sup>-15</sup></td><td>10<sup>-12</sup>–10<sup>-6</sup></td><td>10<sup>-5</sup>–10<sup>-2</sup></td></tr><tr><td>Ro</td><td>10<sup>-6</sup></td><td>10<sup>-3</sup>–10<sup>4</sup></td><td>10<sup>-3</sup>–1</td></tr><tr><td>Rm</td><td>10<sup>3</sup></td><td>10<sup>2</sup>–10<sup>8</sup></td><td>10<sup>3</sup>–10<sup>5</sup></td></tr><tr><td>Pm</td><td>10<sup>-6</sup></td><td>10<sup>-6</sup></td><td>10<sup>-1</sup>–10<sup>5</sup></td></tr></tbody></table>

<footnote><sup>a</sup>All parameters are given in SI units. All values for Earth's core are from Olson [2007]. Diffusivities and rotation rate for diapirs are assumed to be as for the Earth's core. Note that in calculating the non-dimensional number ranges for the diapir not all combinations of $R$ and $U$ are admissible, since $U$ depends on $R$ (equation (1)). For the candidate sinking diapir, see Figure 1.</footnote>
<figure><img src="image_164.png" /><figcaption>**Figure 3.** Schematic illustration of the geometries of the two end-member sinking diapirs. In polar sinking (top) the local gravity is parallel to the rotation vector, whereas in equatorial sinking (right) the local gravity is perpendicular to the rotation vector.</figcaption></figure>

the expected flow pattern equations (4) and (5) is recovered in the whole volume to a very good accuracy. Details of the numerical method can be found in Appendix A.

[19] If the resulting three-dimensional flow **u** is stationary, it is then used as input to a kinematic dynamo code, which solves the induction equation

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B}) + \frac{E}{Pm} \nabla^2 \mathbf{B} \quad (9) $$

where the magnetic field **B** is in units of $\sqrt{\mu_0 \rho_{Fe} R \Omega}$, and $Pm = \nu/\lambda$ is the magnetic Prandtl number. $\mu_0$ is the magnetic permeability of free space ($\mu_0 = 4\pi \times 10^{-7} \text{ H.m}^{-1}$) and $\rho_{Fe}$ is the diapir density ($\rho_{Fe} = 10^4 \text{ kg.m}^{-3}$). Note that $Rm$ (equation (3)) can be rewritten as $Rm = RoPm/E$. A growing magnetic field indicates the presence of a dynamo.

[20] If the Navier-Stokes simulation leads to a time-dependent flow (as it is the case for high $Ro$ and $Re$), we resort to a much slower fully self-consistent dynamo code, which solves the induction equation (9) coupled with the Navier-Stokes equation (8) with the addition of the Lorentz force ($\nabla \times \mathbf{B}) \times \mathbf{B}$. We denote the critical magnetic Reynolds number for the onset of a dynamo by $Rm_c$. To characterize the efficiency of the flow pattern in generating a magnetic field, we also examine the helicity $H = \mathbf{u} \cdot (\nabla \times \mathbf{u})$.

[21] The relevant parameters and non-dimensional numbers are listed in Table 2, including their typical values for Earth's outer core, the possible ranges for the sinking diapir, as well as values used in our numerical models. The most significant variability and uncertainty appears in the size of the diapir $R$ and the viscosity of the surrounding mantle $\eta_s$. It is therefore worth-while writing the control parameters of the models in terms of these two quantities. For example,

substituting equation (1) into equation (7) allows rewriting the Rossby number as

$$ Ro = \frac{1}{3} \frac{\Delta \rho g R}{\eta_s \Omega} \quad (10) $$

Using fixed values for $\Delta\rho$, $g$, $\Omega$, $\nu$ and $\lambda$, the input control parameters $E$, $Ro$ and $Rm$ are given in terms of the diapir size $R$ and the mantle viscosity $\eta_s$ in Tables 1 and 2. For a successful dynamo ($Rm > Rm_c$) to be considered geophysical, it is required to have $E$, $Ro$ and $Rm$ within the (rather broad) estimated range for candidate diapirs (see Figure 1), but also the magnetic Prandtl number must have the value of molten iron, that is $Pm \sim 10^{-6}$ [Olson, 2007]. We limit our dynamical models to the rotational regime of $Ro \le 1$ but we also investigate the ability of the simple laminar Hadamard flow without rotation and inertia to generate a magnetic field.

[22] Two end-member scenarios are examined, one in which the diapir falls on the geographic pole of the planet, the other in which the diapir falls on the equator of the planet. In the first the Hadamard flow forcing is parallel to the planetary rotation axis, in the second the Hadamard flow forcing is perpendicular to the planetary rotation axis (see Figure 3). The consequences of these two end-member scenarios for the dynamo are compared.

# 4. Results

## 4.1. Fluid Flow Inside the Diapir

[23] We begin by reporting the velocity fields within the diapir obtained by our dynamical model. Figure 4a shows a solution for a large Ekman number of $E = 10^{-2}$ and a Rossby number of $Ro = 10^{-1}$ for a polar sinking diapir. The meridional flow consists of a single poloidal roll anti-symmetric with respect to the equator and axisymmetric with respect to the rotation axis, as in the classical Hadamard flow solution without rotation and inertia (equations (4) and (5)). However, in addition the solution contains an azimuthal axisymmetric toroidal flow component anti-symmetric with respect to the equator, which is due to the Coriolis force. Compared to Figure 2, the streamlines are distorted toward the outer part of the sphere in response to the rotational effect. Note that the radial flow is much weaker than the tangential flow. In the rotational regime where $Ro \ll 1$ the flow pattern is weakly dependent on $Ro$. In this regime, the main difference between flow models of two different $Ro$ values is that the magnitude scales linearly with $Ro$. For $Ro \sim 1$ and hence large Reynolds numbers $Re = Ro/E$, the flow is no longer stationary.

[24] For smaller Ekman numbers (Figures 4b and 4c), the flow pattern is qualitatively similar but is now confined to a thinner layer below the boundary of the diapir. Figure 5 shows the thickness of the velocity layer $h$ normalized by the diapir's radius as a function of $E^{-1}$ for $Ro = 10^{-2}$. The decrease in $h$ with decreasing $E$ obeys a power law

$$ \frac{h}{R} = AE^x \quad (11) $$

For a polar sinking diapir the prefactor is $A_p \sim 4$ and the power is $x_p \sim 0.5$ (Figure 5, red). The 1/2 power is suggestive of an Ekman boundary layer control [e.g. Pedlosky, 1987], as may be expected in a system governed by rotational effects. The same results are found for all $Ro < 1$ values (not shown).
<figure><img src="image_165.png" /><figcaption>**Figure 4.** Flow components and helicity for three different *E* values and *Ro* = 10<sup>-1</sup> in the case of a polar sinking diapir. The streamlines with black solid contours (representing counter-clockwise circulation) are superimposed on the *u*<sub>φ</sub> subplots. The sinking geometry is illustrated at the top of the figure where *u*<sub>r</sub> is represented. The mean flow is axisymmetric.</figcaption></figure>
<figure><img src="image_166.png" /><figcaption>**Figure 5.** The non-dimensional thickness of the velocity layer $h/R$ as a function of $E^{-1}$ for $Ro = 10^{-2}$ in log-log scale. The depth is defined where the azimuthal velocity drops to 10% of its maximum value at colatitude 60° (for polar sinking diapir, red) or at co-latitude 90° (for equatorial sinking diapir, green).</figcaption></figure>

[25] We repeat the flow analysis for equatorial sinking diapirs. Note that the visualization is in the planetary reference frame, so the north pole points in the direction of the rotation rate (and not in the diapir sinking latitude). Here the flow is symmetric with respect to the equator (Figure 6). In this configuration, the rotation and sinking axes are perpendicular, so that the axisymmetry of the flow is lost. Indeed, in the equatorial plane, the stationary flow is dominated by an $m = 1$ spiral (see Figure 7), and the effect of the Coriolis force tends to make the flow invariant along the rotation axis (Figure 6). The analysis of the thickness of the velocity layer for equatorial sinking diapir leads to the pre-factor $A_e \sim 2$ and a power of $x_e \sim 0.25$ (equation (11) and Figure 5, green), reminiscent of a Stewartson internal layer control [Stewartson, 1957; Schaeffer and Cardin, 2005].

[26] Because helicity favors dynamo action, we examine the pattern of helicity and its dependence on $E$ and $Ro$. For a polar sinking diapir, the helicity is hemispheric anti-symmetric negative/positive in the northern/southern hemisphere respectively (Figure 4). This helicity pattern is governed by the interaction between the latitudinal flow and the radial shear of the azimuthal flow. As with the flow pattern, the thickness of the layer of significant helicity decreases with decreasing $E$. For equatorial sinking diapirs the helicity is also hemispheric anti-symmetric (Figure 6). Here the hemispheric anti-symmetric helicity pattern is obtained by the interaction between the azimuthal flow and the radial shear of the latitudinal flow.

## 4.2. Magnetic Field Generation

[27] First, we used the analytic Hadamard flow solution without rotation and inertia (equations (4) and (5)). We found no dynamo. Next, we accounted for the effect of rotation: for each flow solution obtained by different combinations of $E$ and $Ro$, we search for the critical magnetic Reynolds number $Rm_c$ above which dynamo action takes place. Figure 8 shows

the results of such a parametric study for a polar sinking diapir. The smallest $Rm_c$ values (most favorable flow morphologies for dynamo action) appear for $E = 10^{-3}$. Weaker rotation yields weaker helicity and subsequent field generation, whereas stronger rotation confines the flow to a thinner layer and thus dynamo action is once again less efficient. In the rotational regime where $Ro \ll 1$, $Rm_c$ does not depend on $Ro$. This is expected because the flow pattern is almost independent of $Ro$ for $Ro \ll 1$. On approach to $Ro \sim 1$, the flow pattern changes and some non-trivial dependence of $Rm_c$ on $Ro$ emerges.

[28] The results of the parametric study for the equatorial sinking diapirs (Figure 9) are similar to those for the polar sinking diapirs. Because the flow of equatorial sinking diapirs is not axisymmetric, one must compute the velocity and magnetic fields with many coupled azimuthal modes, leading to much higher computation times (10 to 30 times longer) than in the polar sinking dynamos. Hence, we decided to focus on $E \le 10^{-3}$ and $0.1 \le Ro < 1$. The smallest $Rm_c$ for equatorial sinking diapir was found for $E = 10^{-4}$. From Figures 8 and 9 we can see that $Rm$ has to be at least larger than 3000 to envision a diapir driven dynamo generation which restricts the candidate diapir domain (Figures 1 and 10).

[29] Figure 11 shows images of the magnetic field for a diapir sinking from the equator for two sets of parameters (($E = 10^{-5}$, $Ro = 0.3$), top; ($E = 10^{-4}$, $Ro = 0.1$), bottom). The magnetic field is characterized by small scale features as the flow within the diapir becomes more turbulent (Figure 11, bottom). Note that in the equatorial sinking case, even at small Ekman number, the magnetic field is strong in a significant volume fraction, despite the confinement of the flow to a narrower layer.

## 5. Magnetic Seed Field Growth Rate and Saturation

[30] We have demonstrated that a sinking diapir can generate a dynamo. Next we estimate the intensity of the magnetic field, in particular whether the sinking time is long enough for the magnetic field to reach saturation before arriving at the bottom of the partially molten mantle. This requires the computation of the non-linear dynamo problem, i.e. the simultaneous solution of the coupled momentum and induction equation. As initial condition, we consider a uniform interplanetary magnetic seed field $B_i = 3 \times 10^{-9}$ T [Riedler et al., 1989]. We monitor the intensity of the dipolar component of the magnetic field at the diapir surface as a function of time from a successful dynamo model of equatorial sinking diapir. We have computed several such dynamos, all of them exhibiting a similar behavior: strong small scale magnetic field within the diapir (Figure 11), with relatively weak dipolar field at its surface. In Figure 12, the time is non-dimensionalized by the time a diapir with radius $R$ needs to travel a distance $R$. The magnetic field is saturated before a 10 km diapir reaches the bottom of a 200 km depth partially molten mantle. For polar diapirs (not shown), the saturated field is slightly weaker and the growth rate is slower compared with the equatorial case.

[31] We also monitor the corresponding magnetic field intensity at the surface of the planet as a function of time (Figure 12, red line). The intensity of the magnetic field
<figure><img src="image_167.png" /><figcaption>**Figure 6.** As in Figure 4 for equatorial sinking diapirs. As in Figure 4, the vertical direction points to the rotation axis of the planet.</figcaption></figure>
<figure><img src="image_168.png" /><figcaption>**Figure 7.** Illustration of the geometry and flow components in the equatorial plane of the planet for $E = 10^{-4}$ and $Ro = 10^{-1}$ for an equatorial sinking diapir.</figcaption></figure>

<figure><img src="image_169.png" /><figcaption>**Figure 8.** Regime diagram of kinematic dynamos for polar sinking diapirs. Black circles are no dynamos, red circles are dynamos. The values of $Pm$ for the dynamo cases in log-scale are represented by the size of the red circles (see legend).</figcaption></figure>
<figure><img src="image_170.png" /><figcaption>Figure 9. As in Figure 8 for equatorial sinking diapirs.</figcaption></figure>

decreases with the distance *d* from the surface of the diapir:

$$B \sim B_0 \frac{d^{\ell-2}}{R} \qquad (12)$$

where *ℓ* is the spherical harmonic degree. We consider the dipolar component of the magnetic field (i.e. *ℓ* = 1) which is the strongest contribution far from the generating region. The magnetic field intensity generated by the diapir deep in the mantle never overcomes the initial magnetic seed field intensity. This means that within the range of parameters of our study (Ro ≤ 1, Ek ≥ 10⁻⁵ and

Pm ≥ 10⁻¹), a dynamo generated by a sinking metallic diapir is unlikely to leave a footprint on the surface paleomagnetic field of a growing planet.

## 6. Discussion

[32] We studied dynamos generated by 3D circulation driven by Hadamard flow on the outer boundary under the influence of the Coriolis force (Ro ≲ 1, E ≪ 1). This is an important constraint, which restricts the candidate diapirs from our models to the region below the Ro ≈ 3 line of Figure 10, where rotational effects are important. In the highly rotational regime of Ro ≪ 1 the flow morphology as

<figure><img src="image_171.png" /><figcaption>Figure 10. A refined diagram of candidate diapirs for dynamo action. We constrain the candidate diapirs from the results of our models to the dark shaded zone. Successful dynamos are represented with green symbols (squares for polar and circles for equatorial sinking diapirs). The blue solid line for Rm = 3000 represents the minimum value for a successful dynamo from our models. The red dashed line for Ro = 1 represents the limit of the rotational regime investigated in our study. We also represent the Ro = 0.1 value. Larger but moderate Ro values marked by the Ro = 3 line (not studied here) may possibly lead to dynamo action. We extend the domain to ηs = 10¹⁴ Pa.s where dynamos were observed in our models.</figcaption></figure>
<figure><img src="image_172.png" /><figcaption>**Figure 11.** Example of dynamo magnetic fields for an equatorial sinking diapir. (left) The φ-component of the magnetic field in an arbitrary meridional cross-section and (right) the Z-component of the magnetic field in the equatorial plane. (top) $E = 10^{-4}$, $Ro = 0.1$, $Pm = 20$, $Re = 10^3$ and (bottom) $E = 10^{-5}$, $Ro = 0.3$, $Pm = 0.5$, $Re = 3 \times 10^4$.</figcaption></figure>

well as $Rm_c$ do not depend on $Ro$. The behavior of the system for $Ro \sim 1$ is complicated to investigate due to computational limitations. In these systems turbulence becomes important, the flow is unsteady and temporal fluctuations dominate the dynamics. This regime can eventually lead to more complex flows and favor dynamo generation, although it is unlikely [Peyrot et al., 2007; Ponty and Plunian, 2011].

[33] Can sinking diapirs produce dynamos in nature? In our numerical simulations, all physical parameters are in the geophysical range except for the liquid iron viscosity, which influences both Ekman and magnetic Prandtl number. Unfortunately, simulations with lower $E$ are computationally very costly. In Figures 8 and 9 the magnetic Prandtl number $Pm$ is represented by the size of the circles in the dynamo cases, so that larger circles represent smaller $Pm$, closer to

the geophysical value. Both these figures show that we have obtained dynamos at smaller $Pm$ values when decreasing $E$. In the regime which is accessible for our simulations, the lowest $Pm$ values are on the order of 0.1, five orders of magnitude too large. Such computational problems of obtaining Earth-like $E$ and $Pm$ values are ubiquitous to all numerical dynamo studies [Christensen and Aubert, 2006; Christensen and Wicht, 2007].

[34] From Figures 8 and 9, we obtain the critical magnetic Reynolds number $Rm_c$ for each combination of $E$ and $Ro$ values. This global $Rm_c$ scales with $R^3$ (equations (1) and (3)). Lowering the Ekman number confines the flow to a narrow velocity layer at the outer part of the diapir, which results in larger $Rm_c$. This motivates defining a local critical
<figure><img src="image_173.png" /><figcaption>**Figure 12.** Intensity of the dipolar magnetic field generated by an equatorial sinking diapir as a function of time with $Ro = 1$, $E = 10^{-5}$, $Rm = 10^5$, $Pm = 1$. The intensity of the magnetic field is non-dimensiona- lized by $\sqrt{\mu_0\rho_{Fe}R\Omega\rho_{Fe}} = 10^4$ kg.m$^{-3}$ and $R = 10$ km. The blue line represents the evolution of the dipolar magnetic field at the surface of the diapir, and the red line at the surface of the planet.</figcaption></figure>

magnetic Reynolds number $Rm_c^L$ based on the thickness of the velocity layer $h$:

$$Rm_c^L = \frac{h}{R}Rm_c \qquad (13)$$

[35] Using the best fits to $h$ from Figure 5, we plot $Rm_c^L$ as a function of $E^{-1}$ in Figure 13. It is difficult to extrapolate a relation between $Rm_c^L$ and $E$. Based on this figure, $Rm_c^L$ appears to have an asymptotic behavior with decreasing $E$, with polar sinking dynamos seeming to be favorable. How- ever, the thickness $h$ of that layer is much larger in the equatorial sinking case ($h \sim E^{1/4}$) than in the polar sinking case ($h \sim E^{1/2}$). It thus appears that dynamo generation from a sinking diapir could be favored by equatorial sinking dia- pirs for realistic Ekman numbers. Even in the very pessi- mistic case, where $Rm_c \sim 10^5$ (two orders of magnitude larger than our optimal model) for realistic Ekman numbers, giant diapirs (radius from 30 to 100 km) can still produce strong magnetic fields during their sinking. In addition, the turbulent fluctuations may limit the decrease of $h$ with decreasing $E$ (Figure 5), leading to the decrease of $Rm_c$ at low $E$ (equation (13)). However, lowering the Ekman num- ber will lead to smaller scale magnetic field generation (as suggested by Figure 11), leading to a possibly lower mag- netic field at the surface of the planet.

## 7. Conclusions and Perspectives

[36] Our results show that the flow pattern driven by a diapir sinking through a partially molten mantle within a

rotating planet can generate a magnetic field. This dynamo generation seems more favorable for a diapir sinking from the equator than from the planet's rotational pole. Large diapirs ($R > 10$ km) sinking in a mantle with a viscosity ranging from $10^9$ to $10^{14}$ Pa.s (see Table 1) provide plausible conditions to generate a transient dynamo that may have occurred in the early history of terrestrial planets or imme- diately after a giant impact and the subsequent core merging process. However, the magnetic field generated by a single diapir in our models is too weak to contribute to the paleo- magnetic field recorded on Mars or on the Moon.

[37] Several questions arise from our results:

[38] 1. What is the influence of the intensity of the mag- netic seed field on the feasibility of a diapir driven dynamo? We have considered in this study the very weak homoge- neous interplanetary magnetic seed field. However, mag- netic fields are expected to be generated in planetary cores during the early evolution of planets and moons [Monteux et al., 2011a]. The presence of a stronger magnetic field during early planetary evolutions may change the flow inside the diapir, leading to faster and stronger magnetic field amplification.

[39] 2. How will multiple diapirs affect the early planetary magnetic fields? Large diapirs are not isolated events during planetary differentiations. The presence of multiple diapirs may influence the sinking dynamics of each diapir [Manga and Stone, 1993] and hence the dynamo generation.

[40] 3. Can the magnetic field of the sinking diapirs help to start the planetary dynamo? If the convective planetary dynamo is subcritical [Christensen et al., 2001], i.e. the dynamo cannot start without a strong initial field, the core
<figure><img src="image_174.png" /><figcaption>**Figure 13.** Local magnetic Reynolds number $Rm_c^L$ as a function of $E^{-1}$. $Rm_c^L$ is obtained by combining the results for $Rm_c$ from Figures 8 and 9 (choosing the most favorable $Ro$ value for each value of $E$) and the fits for $h$ from Figure 5.</figcaption></figure>

merging process could bring the required magnetic field intensity inside the planetary core to trigger a dynamo.

## Appendix A: Numerical Method

[41] Our three-dimensional spherical code uses a second order finite differences scheme in the radial direction and pseudo-spectral spherical harmonic expansion in the tangential direction using the high performance SHTNs library (N. Schaeffer, Efficient spherical harmonic transforms aimed at pseudo-spectral numerical simulations, preprint, 2012, arXiv:1202.6522). The time stepping uses a semi-implicit Crank-Nicholson scheme for the diffusive terms, while the non-linear terms are handled by an Adams-Bashforth scheme (second order in time). To calculate the 3D flow, the momentum equation is time stepped in the spherical fluid shell with imposed velocity field at the spherical outer boundary. For kinematic dynamo simulations the induction equation is solved with an insulator outside the outer boundary and an imposed stationary 3D velocity field. The full dynamo problem can also be solved with coupled induction and momentum equations.

[42] We want to emphasize that there is no solid inner core: our code supports fluid that fills the whole sphere. The method to overcome numerical instability near $r = 0$ is based on the fact that the minimum length-scale resolved by a spherical harmonic expansion truncated at degree $\ell_{max}$ is about $\pi r / \ell_{max}$. Hence when $r \to 0$, we resolve smaller and smaller angular length-scale, which is at best useless, and can lead to numerical instabilities. We overcome this limitation with a spherical harmonic truncation $\ell_{tr}$ that depends on $r$. Specifically, we use

$$ \ell_{tr}(r) = 1 + (\ell_{max} - 1) \frac{r}{r_{max}} \quad (A1) $$

Although the finite difference scheme has a theoretical error that scales like $1/r$, the numerical solutions of the velocity field and magnetic field are perfectly smooth near $r = 0$, with flow and magnetic field that are allowed to cross $r = 0$. We use $\ell_{max}$ ranging from 20 to 90 and the number of radial grid points ranges from 100 to 300 depending on the strength of inertial effects (measured by the hydrodynamic Reynolds number $Re = Ro/E$).

[43] In the case of a diapir falling on the planetary pole, the boundary conditions are symmetric by rotation around the rotation axis of the planet. This implies that the stationary flow is axisymmetric (when there are no instabilities). When computing the kinematic dynamo problem with an axisymmetric flow, the various azimuthal wave numbers $m$ of the magnetic field are independent. Moreover, it is known from the Cowling theorem that the magnetic field must include non-axisymmetric ($m > 1$) terms [Cowling, 1934]. Hence we search for magnetic field with $m = 1, m = 2$ and $m = 3$ separately. We find that most often the first growing magnetic field mode is the one with $m = 2$.

[44] In the case of a diapir falling on the planetary equator, the boundary conditions are not symmetric by rotation around the rotation axis of the planet. However, they are purely $m = 1$, and we find that the stationary flow, without inertial effects (small $Ro/E$) is mainly a spiralling $m = 1$ flow. When inertial effects are not negligible, all azimuthal wave numbers are present in the flow, leading to longer computations, where we usually set $m_{max} = \ell_{max}/2$, which seems to be a reasonable truncation scheme when looking at the resulting spectra. The kinematic dynamo computation also requires all azimuthal wave numbers to be computed at once because they are coupled by the $m > 0$ flow. We use the same truncation as for the flow, and for some cases we verified that larger truncations did not change the dynamo onset threshold significantly.
**[45]** Acknowledgments. The authors thank B. Langlais and T. Seon for useful discussions. The authors also thank the anonymous reviewers for constructive comments. J.M. is funded by Agence Nationale de Recherche (Accretis decision ANR-10-PDOC-001-01). Part of the numerical simulations were run at the Service Commun de Calcul Intensif de l'Observatoire de Grenoble (SCCI).

## References

Agnor, C. B., R. M. Canup, and H. F. Levison (1999), On the character and consequences of large impacts in the late stage of terrestrial planet formation, *Icarus*, **142**, 219–237, doi:10.1006/icar.1999.6201.

Amit, H., P. Olson, and U. Christensen (2007), Tests of core flow imaging methods with numerical dynamos, *Geophys. J. Int.*, **168**, 27–39.

Batchelor, G. (1967), *An Introduction to Fluid Dynamics*, Cambridge Univ. Press, Cambridge, U. K.

Benz, W., W. L. Slattery, and A. G. W. Cameron (1987), The origin of the Moon and the single-impact hypothesis, *Icarus*, **71**, 30–45, doi:10.1016/0091-1035(87)90160-6.

Christensen, U., and J. Aubert (2006), Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields, *Geophys. J. Int.*, **166**, 97–114.

Christensen, U., and J. Wicht (2007), Numerical dynamo simulations, in *Core Dynamics, Treatise Geophys.*, vol. 8, edited by P. Olson, pp. 245–282, Elsevier, Oxford, U. K.

Christensen, U. R., et al. (2001), A numerical dynamo benchmark, *Physics Earth Planet. Inter.*, **128**(1–4), 25–34, doi:10.1016/S0031-9201(01)00275-8.

Cowling, T. (1934), The magnetic field of sun spots, *Mon. Not. R. Astron. Soc.*, **94**, 39–48.

Dahl, T. W., and D. J. Stevenson (2010), Turbulent mixing of metal and silicate during planet accretion and interpretation of the Hf-W chronometer, *Earth and Planet. Sci. Lett.*, **295**, 177–186, doi:10.1016/j.epsl.2010.03.038.

Davies, G. F. (1982), Ultimate strength of solids and formation of planetary cores, *Geophys. Res. Lett.*, **9**, 1267–1270.

Deguen, R., P. Olson, and P. Cardin (2011), Experiments on turbulent metal-silicate mixing in a magma ocean, *Earth Planet. Sci. Lett.*, **310**, 303–313, doi:10.1016/j.epsl.2011.08.041.

Dudley, M., and R. James (1989), Time-dependent kinematic dynamos with stationary flows, *Proc. R. Soc. London, Ser. A*, **425**, 407–429.

Gailitis, A. (1970), Self-excitation of a magnetic field by a pair of ring vortices [in Russian], *Mag. Gidrodin.*, **6**, 19–22. [*Magnetohydrodynamics*, Engl. Transl., **6**, 14–17].

Gerya, T. V., and D. A. Yuen (2007), Robust characteristics method for modelling multiphase visco-elastic thermo-mechanical problems, *Phys. Earth Planet. Inter.*, **163**, 83–105.

Gubbins, D., C. Barber, S. Gibbons, and J. Love (2000a), Kinematic dynamo action in a sphere: I. Effects of differential rotation and meridional circulation on solutions with axial dipole symmetry, *Proc. R. Soc. London, Ser. A*, **456**, 1333–1353.

Gubbins, D., C. Barber, S. Gibbons, and J. Love (2000b), Kinematic dynamo action in a sphere: II. Symmetry selection, *Proc. R. Soc. London, Ser. A*, **456**, 1669–1683.

Guervilly, C., and P. Cardin (2010), Numerical simulations of dynamos generated in spherical Couette flows, *Geophys. Astrophys. Fluid Dyn.*, **104**, 221–248, doi:10.1080/03091920903550955.

Hadamard, J. (1911), Mouvement permanent lent d'une sphère liquide et visquose dans un liquide visqueux, *C. R. Acad. Sci.*, **152**, 1735–1738.

Höink, T., J. Schmalzl, and U. Hansen (2005), Formation of compositional structures by sedimentation in vigorous convection, *Phys. Earth Planet. Inter.*, **153**, 11–20.

Höink, T., J. Schmalzl, and U. Hansen (2006), Dynamics of metal-silicate separation in a terrestrial magma ocean, *Geochem. Geophys. Geosyst.*, **7**, Q09008, doi:10.1029/2006GC001268.

Honda, R., H. Mizutani, and T. Yamamoto (1993), Numerical simulation of Earth's core formation, *J. Geophys. Res.*, **98**, 2075–2089.

Karato, S.-I., and V. Murthy (1997), Core formation and chemical equilibrium in the Earth—I. Physical considerations, *Phys. Earth Planet. Inter.*, **100**, 61–79.

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

Kleine, T., C. Münker, K. Mezger, and H. Palme (2002), Rapid accretion and early core formation on asteroids and the terrestrial planets from Hf-W chronometry, *Nature*, **418**, 952–955.

Le Bars, M., M. Wieczorek, O. Karatekin, D. Cébron, and M. Laneuville (2011), An impact-driven dynamo for the early Moon, *Nature*, **479**, 215–218.

Manga, M., and H. A. Stone (1993), Buoyancy-driven interactions between two deformable viscous drops, *J. Fluid Mech.*, **256**, 647–683, doi:10.1017/S0022112093002915.

Moffatt, H. (1978), *Magnetic Field Generation in Electrically Conducting Fluids*, Cambridge Univ. Press, Cambridge, U. K.

Monteux, J., Y. Ricard, N. Coltice, F. Dubuffet, and M. Ulvrova (2009), A model of metal-silicate separation on growing planets, *Earth Planet. Sci. Lett.*, **287**, 353–362.

Monteux, J., A. M. Jellinek, and C. L. Johnson (2011a), Why might planets and moons have early dynamos?, *Earth. Planet. Sci. Lett.*, **310**, 349–359, doi:10.1016/j.epsl.2011.08.014.

Monteux, J., M. Jellinek, and C. L. Johnson (2011b), Core merging after the Martian giant impact, *Lunar Planet. Sci. Conf., XLII*, Abstract 1665.

Moss, D. (2008), Simple laminar dynamos: From two rolls to one, *Geophys. Astrophys. Fluid Dyn.*, **102**, 195–203.

Olson, P. (2007), Overview on core dynamics, in *Core Dynamics, Treatise Geophys.*, vol. 8, edited by P. Olson, pp. 1–30, Elsevier, Oxford, U. K.

Olson, P., U. Christensen, and G. Glatzmaier (1999), Numerical modeling of the geodynamo: Mechanisms of field generation and equilibration, *J. Geophys. Res.*, **104**, 10,383–10,404.

Olson, P., I. Sumita, and J. Aurnou (2002), Diffusive magnetic images of upwelling patterns in the core, *J. Geophys. Res.*, **107**(B12), 2348, doi:10.1029/2001JB000384.

Pedlosky, J. (1987), *Geophysical Fluid Dynamics*, Springer, New York.

Peyrot, M., F. Plunian, and C. Normand (2007), Parametric instability of the helical dynamo, *Phys. Fluids*, **19**, 054109, doi:10.1063/1.2734118.

Ponty, Y., and F. Plunian (2011), Transition from large-scale to small-scale dynamo, *Phys. Rev. Lett.*, **106**, 154502, doi:10.1103/PhysRevLett.106.154502.

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

Riedler, W., K. Schwingenschuh, D. Moehlmann, V. N. Oraevskii, E. Eroshenko, and J. Slavin (1989), Magnetic fields near Mars—First results, *Nature*, **341**, 604–607, doi:10.1038/341604a0.

Rubie, D., H. Melosh, J. Reid, C. Liebske, and K. Righter (2003), Mechanisms of metal-silicate equilibration in the terrestrial magma ocean, *Earth Planet. Sci. Lett.*, **205**, 239–255.

Safronov, V. S. (1978), The heating of the Earth during its formation, *Icarus*, **33**, 3–12, doi:10.1016/0019-1035(78)90019-2.

Samuel, H. (2012), A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans, *Earth Planet. Sci. Lett.*, **313**, 105–114, doi:10.1016/j.epsl.2011.11.001.

Samuel, H., and P. Tackley (2008), Dynamics of core formation and equilibration by negative diapirism, *Geochem. Geophys. Geosyst.*, **9**, Q06011, doi:10.1029/2007GC001896.

Schaeffer, N., and P. Cardin (2005), Quasigeostrophic model of the instabilities of the stewartson layer in flat and depth-varying containers, *Phys. Fluids*, **17**(10), 104111, doi:10.1063/1.2073547.

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

Shannon, M. C., and C. B. Agee (1996), High pressure constraints on percolative core formation, *Geophys. Res. Lett.*, **23**, 2717–2720.

Spells, K. (1952), A study of circulation patterns within liquid drops moving through a liquid, *Proc. Phys. Soc. B*, **65**, 541–546.

Šrámek, O., Y. Ricard, and F. Dubuffet (2010), A multiphase model of core formation, *Geophys. J. Int.*, **181**, 198–220, doi:10.1111/j.1365-246X.2010.04528.x.

Stevenson, D. (1981), Models of the Earth's core, *Science*, **214**, 611–619.

Stevenson, D. J. (2003), Planetary science: Mission to Earth's core— A modest proposal, *Nature*, **423**(6937), 239–240, doi:10.1038/423239a.

Stewartson, K. (1957), On almost rigid rotations, *J. Fluid Mech.*, **3**(01), 17–26, doi:10.1017/S0022112057000452.

Tilgner, A. (2005), Precession driven dynamos, *Phys. Fluids*, **17**(3), 034104, doi:10.1063/1.1852576.

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

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

Touboul, M., T. Kleine, B. Bourdon, H. Palme, and R. Wieler (2007), Late formation and prolonged differentiation of the Moon inferred from W isotopes in lunar metals, *Nature*, **450**, 1206–1209.

Ulvrová, M., N. Coltice, Y. Ricard, S. Labrosse, F. Dubuffet, J. Velímský, and O. Šrámek (2011), Compositional and thermal equilibration of particles, drops, and diapirs in geophysical flows, *Geochem. Geophys. Geosyst.*, **12**, Q10014, doi:10.1029/2011GC003757.

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.

Yoshino, T., M. J. Walter, and T. Katsura (2003), Core formation in planetesimals triggered by permeable flow, *Nature*, **422**, 154–157.
<figure><img src="image_175.png" /></figure>

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

# Dynamics of core merging after a mega-impact with applications to Mars' early dynamo

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

CrossMark

J. Monteux<sup>a,b,*</sup>, A.M. Jellinek<sup>b</sup>, C.L. Johnson<sup>b,c</sup>

## ARTICLE INFO

**Article history:**

Received 3 October 2012

Revised 23 April 2013

Accepted 9 May 2013

Available online 18 May 2013

**Keywords:**

Mars, Interior

Magnetic fields

Impact processes

Planetary dynamics

## ABSTRACT

A giant impact occurring within the first 500 Myr of martian history may have been responsible for the dichotomy between the northern lowlands and the southern highlands and may have influenced the initiation or cessation of early and short-lived core dynamo. We hypothesize that a significant volume of metallic iron from a differentiated impactor merged with a preexisting martian core. We investigate the dynamics and thermal effects of this core merging, assuming that the impactor's core sank as a single metallic diapir through a solid mantle. We explore the consequences of this process for dynamo action and for Mars' magnetic field history. For large impacts (with radii larger than 100 km) and plausible mantle viscosities, merging is expected to occur in less than 1 Myr. Depending on the temperature-dependence of the mantle viscosity, viscous dissipation within the diapir may be very large. Where thermal mixing of the hot diapir into a preexisting core is complete, merging can increase the temperature gradient to the surrounding mantle and consequently drive a dynamo until this additional heat is transferred to the mantle, which takes on the order of 100 Myr. If merging leads to strong thermal stratification in the core, however, dynamo action may be inhibited.

© 2013 Elsevier Inc. All rights reserved.

## 1. Introduction

Mega impacts (with an impactor/target radius ratio between one tenth and 1) probably played an important role in the late history of terrestrial planetary accretion. An impact between the pro-to-Earth and a Mars-size protoplanet can, for example, explain the Earth-Moon system (Hartmann and Davis, 1975) and a large impact removing part of the silicate mantle is a hypothesis for the high iron/silicate ratio on Mercury (Smith, 1979; Benz et al., 1988). Among the hypotheses for the origin of the martian dichotomy including an endogenic origin (Elkins-Tanton et al., 2003; Roberts and Zhong, 2006; Citron and Zhong, 2012) or a plate-technology feature (Sleep, 1994), an exogenic origin by a mega impact that displaced crustal material from the northern to the southern hemisphere seems to be the most plausible candidate (Wilhelms and Squyres, 1984; Nimmo et al., 2008; Andrews-Hanna et al., 2008; Marinova et al., 2008, 2011). Models for the martian impact suggest that the impactor was 800–1300 km in radius and hit the planet with a speed comparable to or larger than the martian escape velocity (i.e. $v_{imp} > 5 \text{ km s}^{-1}$) within the first 500 Myr of martian history (Frey, 2006). Such a dramatic event can generate a

debris disk around Mars that could have re-accreted and formed the martian moons (Rosenblatt, 2011).

In addition to a history of large impacts, Earth, Mercury and Mars have, or have had, an internally generated magnetic field. Mercury and the Earth have active core dynamos, whereas widespread crustal magnetism strongly suggest that Mars had an early internally-generated magnetic field (Acuña et al., 1999; Hood et al., 2003; Lillis et al., 2008a) that ceased by around 4.0 Ga (Acuña et al., 1999; Johnson and Phillips, 2005; Lillis et al., 2008b). The timing of the initiation of the martian dynamo is difficult to constrain and strongly depends on the differentiation processes that occurred during the first million years of martian history (Monteux et al., 2011). The cause of the cessation of the martian dynamo is also still currently debated. Recent models explore the effects of large impacts on the dynamo generation process, and in particular on the cessation of dynamo action as a result of a reduction of the core-mantle boundary (CMB) heat-flux (Roberts et al., 2009; Watters et al., 2009; Roberts and Arkani-Hamed, 2012). Shock heating within the core can also increase the CMB heat flow and create a thermal stratification that prevents heat loss from the inner part of the core, inhibiting core convection (Arkani-Hamed and Olson, 2010). Alternatively, other models show that the thermal anomaly induced by a large impact and the formation of a hot molten iron layer from the impactor's core at the CMB can favor dynamo generation (Reese and Solomatov, 2010).

<footnote><sup>a</sup> Laboratoire de Planétologie et de Géodynamique de Nantes, UMR CNRS, Université de Nantes, France</footnote>

<footnote><sup>b</sup> Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, Canada</footnote>

<footnote><sup>c</sup> Planetary Science Institute, Tucson, AZ, United States</footnote>

<footnote>* Corresponding author at: Laboratoire de Planétologie et de Géodynamique de Nantes, UMR CNRS, Université de Nantes, France.</footnote>

<footnote>E-mail address: julien.monteux@univ-nantes.fr (J. Monteux).</footnote>
At the time of the proposed giant impact, both Mars and the impactor were probably differentiated (Yoshino et al., 2003). Models suggest that although some material was ejected far from Mars, the majority of the mass of the impactor's core was retained within the planet and merged with the pre-impact martian core (Canup, 2004; Ćuk and Stewart, 2012). The aim of this study is to characterize the dynamics of core merging as a result of a diapiric descent of molten iron (Monteux et al., 2009) on a Mars-size planet. In addition, we investigate how the processes of impact and core merging might influence the thermal regime of Mars' core and, in turn, magnetic field generation (cf., Monteux et al., 2011).

## 2. Thermo-chemical state before the martian mega-impact

### 2.1. Pre-impact interior of Mars

The initial structure and thermal state of a growing planet is determined by the characteristics of its accretion from chondritic material (Safronov, 1978; Kaula, 1979; Agee, 1997). During accretion, heating driven by a combination of the dissipation of impact energy and the decay of short lived radionuclides such as $^{26}$Al and/ or $^{60}$Fe (Yoshino et al., 2003; Monteux et al., 2007) increases the mean internal temperature and gives rise to a radial temperature gradient that depends on the accretion rate relative to the rate of radiative cooling to space (Kaula, 1979; Senshu et al., 2002). If the growth rate is very high in comparison to surface cooling, this heating can ultimately cause partial or complete melting of the chondritic material (Yoshino et al., 2003) and lead to extensive metal/silicate separation (Tonks and Melosh, 1992; Senshu et al., 2002; Monteux et al., 2009).

Hf/W chronology suggests that core formation happened during the first 10–30 myr of Mars' history (Lee and Halliday, 1997; Nimmo and Kleine, 2007). Such a rapid process involves extensive melting potentially enhanced by radiogenic heating as a result of the decay of short-lived radionuclides (Yoshino et al., 2003), impact heating (Tonks and Melosh, 1992; Senshu et al., 2002; Monteux et al., 2009) and gravitational energy conversion during metal/silicate separation (Stevenson, 1989; Ricard et al., 2009). Metal/silicate separation can occur via a wide range of phenomena such as percolation (Shannon and Agee, 1996), the sedimentation of metallic rain through a magma ocean (Rubie et al., 2003; Höink et al., 2005) or a large diapir sinking through a solid mantle after an impact (Tonks and Melosh, 1992; Monteux et al., 2009). Whatever the mechanism, Mars' internal structure characterized by a ~1700 km diameter Fe-core was mostly established within ~10 Myr of the planet's formation (Yoder et al., 2003) (cf., Fig. 1a).

The gravitational heat released during martian core formation was partitioned between the planet's core and mantle. The fraction of gravitational heat taken up by the metal or the silicate fraction depends strongly on the rheology of the planet and on the segregation mechanisms (Samuel and Tackley, 2008; Monteux et al., 2009; Ke and Solomatov, 2009). The combined processes leading to core formation yield a wide range of possible early thermal states, depending on the nature and timescale of the core formation process and the heat transfer properties of Mars' early mantle. In particular, the core could initially have had a temperature close to the deep mantle temperature if thermal equilibration was efficient. Alternatively, it could have been hotter than the mantle if the gravitational potential energy released during core formation was largely retained within the core itself, a situation which would lead to potentially strong cooling to the mantle (Fig. 1a).

### 2.2. Interior structure of the impactor

In this study, we consider impactors with a radius in the range 200–800 km. The lower bound for our range is motivated by the

impactor size needed to create large impact basins such as Hellas or Utopia, and the upper bound is motivated by the minimum impactor radius needed in exogenic models for the dichotomy boundary (Marinova et al., 2008). Assuming that both the impactor and the target body had chondritic compositions, their volumetric metal fractions, $f_0$, should be similar (we consider that the impactor has the same metal content as Mars and we use $f_0$ = 12.5% (Stevenson, 2001)). Hence, for 200–1300 km diameter impactors, an additional volume of core material with a radius between 100 km and 700 km merges with the preexisting core (Fig. 1).

## 3. Thermo-chemical state after a mega-impact

### 3.1. Mantle heating and melting

Kinetic energy of the impactor is dissipated as a result of the irreversible work done by shock waves in damaging crustal rocks (Senshu et al., 2002; Monteux et al., 2011) as well as heating and melting the target material. This dissipation process is a complex mechanism that is still poorly constrained specially for giant impact events. In our models, we consider that post-impact heating and melting mostly occurs within a spherical region with a volume $V_{ic}$ (and a radius $R_{ic}$) that is typically taken to be 3 times larger than the volume $V_{imp}$ of the impactor itself (O'Keefe and Ahrens, 1977; Croft, 1982; Pierazzo et al., 1997). The energy available to heat and melt the target planet is $\Delta E = \gamma m_{imp} v_{imp}^2 / 2$, where $\gamma$ is the fraction of the kinetic energy of the impactor ultimately dissipated to heat up the mantle (O'Keefe and Ahrens, 1977), $m_{imp}$ is the impactor mass and $v_{imp}$ is the impact velocity. The energy needed to melt a silicate volume $V_{ic}$ is $\Delta E_{m,si} = \rho_{si}V_{ic}L_{si}$, where $\rho_{si}$ and $L_{si}$ are the density and latent heat of the silicate material. The energy needed to melt the impactor core is $\Delta E_{m,Fe} = \rho_{Fe}f_0V_{imp}L_{Fe}$ with $L_{Fe}$ and $\rho_{Fe}$ are the density and latent heat of the impactor's core.

In the heated region the temperature increases uniformly from an initial value $T_0$ by an amount $\Delta T_0$ (Fig. 1c). The excess temperature $\Delta T_0$ decreases rapidly and smoothly with distance $r$ from the boundary of the isothermal anomaly as approximately $\Delta T_0(R_{ic}/r)^m$ (Fig. 1c). Following Senshu et al. (2002), and fitting the decay of peak pressure with distance away from the edge of the isobaric core $m \approx 4.4$ (Monteux et al., 2007). The energy needed to increase the temperature inside and outside the isobaric core is $\Delta E_{th} = h_m \rho_{si} V_{ic} C_{p,si} \Delta T_0$, where $C_{p,si}$ is the specific heat of the silicate material and $h_m$ is a geometric parameter representing the amount of heat that is used to increase the temperature inside and outside the isobaric core relative to the amount of heat used to increase the temperature by $\Delta T_0$ within the isobaric core (Senshu et al., 2002; Monteux et al., 2011) and

$$h_m = 1 + \frac{3(2m - 5)}{2(m - 3)(m - 2)} \approx 2.7 \quad (1)$$

Hence, from the following energy balance:

$$\Delta E = \Delta E_{th} + \Delta E_{m,si} + \Delta E_{m,Fe} \quad (2)$$

and using $m_{imp} = \rho_0 V_{imp} = \rho_0 V_{ic}/3$, we obtain:

$$\Delta T_0 = \frac{1/6\gamma\rho_0 v_{imp}^2 - L_{si}\rho_{si} - f_0/3L_{Fe}\rho_{Fe}}{h_m\rho_{si}C_{p,si}} \quad (3)$$

Assuming that the impact velocity is equal to the escape velocity of the planet we obtain a minimum estimate for the kinetic energy where $v_{imp} = v_{esc} = \sqrt{2gR}$ with $g = 4/3\pi G\rho_0 R$, G the gravitational constant and R the radius of the target planet. After some algebra:

$$\Delta T_0 = \frac{4/9\pi\gamma\rho_0^2GR^2 - L_{si}\rho_{si} - \frac{f_0}{3}L_{Fe}\rho_{Fe}}{h_m\rho_{si}C_{p,si}}, \quad (4)$$
<figure><img src="image_178.png" /><figcaption>**Fig. 1.** Schematic representation of the thermo-chemical evolution following a mega-impact on a Mars-size body. After the impact (a) and the subsequent shockwave propagation (b), the impactor's core, radius $R_{Fe}$, is buried within the martian mantle and a large thermal anomaly remains (c). The post-impact temperature increase is constant within a volume delimited by the dashed line and rapidly decreases away from it. The dense metallic material from the impactor rapidly sinks towards the pre-impact martian core. During the sinking, potential energy is converted into heat in the mantle via viscous dissipation. Two sinking regimes occur (c and d) depending on the distance $h$ between the impactor and the martian core (see Section 4). Finally the intruded metallic material merges with the preexisting core, spreads at its top (e) and eventually mixes with it later (f). The hatched volumes represent the metallic material from the martian core and the impactor. Gray shading represents temperature of the silicate (unhatched) or metallic (hatched) material and shows schematically the evolution of the thermal anomaly in the mantle during the sinking and merging phases (see also Fig. 5, left column).</figcaption></figure>

For a Mars-sized planet with a mantle temperature close to its solidus the impact heating is always sufficiently large to completely melt the impacted mantle volume $V_{ic}$ (i.e., $\Delta T_0 > 0$). For smaller or much colder planets, impact heating may only partially melt the impacted region, depending on the slope of the mantle solidus in temperature-composition space (Ernst et al., 2010; Roberts and Barnouin, 2012).

Noting that $f_0$ is O ($10^{-1}$), the term representing the energy needed to melt the impactor's core has a minor influence compared to the mantle melting term in Eq. (4). Hence, Eq. (4) can be simplified to give:

$$ \Delta T_0 \approx \frac{4/9\pi\gamma\rho_0^2GR^2 - L_{Si}\rho_{Si}}{h_m\rho_{Si}C_{p,si}}, \qquad (5) $$

Indeed using the set of parameters in Table 1 and making the assumption that the martian mantle before the impact is solid and close to its solidus, $\Delta T_0 \sim 400$ K. Including the impactor's core melting term in Eq. (4) would only reduce $\Delta T_0$ by $\sim 10$ K.

## 3.2. The fate of the impactor's material

During the mega-impact, a significant fraction of the impactor's silicate mantle is ejected far from the impact site to form a disk of

debris, but most of the impactor's core is retained within the target planet (Canup, 2004; Čuk and Stewart, 2012) (cf., Fig. 1b and c). For simplicity we assume that the full volume of the impactor's core is retained within the martian mantle after the giant impact. This assumption is reasonable for impact velocities close to the escape velocity of the target planet (Asphaug, 2010). We assume a solid mantle prior to impact, and that the impactor core sinks as a single, spherical diapir.

## 4. Dynamic models of diapir sinking

Once buried below the surface, the dense metallic core of the impactor sinks towards the center of the target planet (Fig. 1). Dissipation of the work done by buoyancy forces driving motion of this diapir occurs by a combination of viscous, viscoplastic (Samuel and Tackley, 2008), elasto-plastic (Gerya and Yuen, 2007) and fracturing (Davies, 1982; Stevenson, 2003) processes and causes heating in the mantle and/or the diapir, depending on the deformation regime. Where the effective mantle viscosity is very large in comparison to that of the diapir dissipation is concentrated in the mantle. As viscosity variations decline in response to heat transfer from the diapir to the mantle, however, dissipation will increasingly oc-
<table><caption>**Table 1**

Typical parameter values for numerical models.</caption><tr><td>Mars radius</td><td>R</td><td>3400 km</td></tr><tr><td>Mars core radius</td><td>R<sub>c</sub></td><td>1700 km</td></tr><tr><td>Thickness of the core thermal boundary layer</td><td>δ</td><td>0–0.4R<sub>c</sub></td></tr><tr><td>Impactor radius</td><td>R<sub>imp</sub></td><td>200–800 km</td></tr><tr><td>Impactor core radius (=diapir size)</td><td>R<sub>Fe</sub></td><td>100–400 km</td></tr><tr><td>Initial gap thickness</td><td>h<sub>0</sub></td><td>580 km</td></tr><tr><td>Average density of the planet</td><td>ρ<sub>0</sub></td><td>4060 kg m<sup>-3</sup></td></tr><tr><td>Iron density</td><td>ρ<sub>Fe</sub></td><td>8000 kg m<sup>-3</sup></td></tr><tr><td>Silicate density</td><td>ρ<sub>Si</sub></td><td>3500 kg m<sup>-3</sup></td></tr><tr><td>Density difference (=ρ<sub>Fe</sub> − ρ<sub>Si</sub>)</td><td>Δρ<sub>0</sub></td><td>4500 kg m<sup>-3</sup></td></tr><tr><td>Average coefficient of thermal expansion</td><td>α</td><td>4.5 × 10<sup>-5</sup> K<sup>-1</sup></td></tr><tr><td>Iron coefficient of thermal expansion</td><td>α<sub>Fe</sub></td><td>1.5 × 10<sup>-5</sup> K<sup>-1</sup></td></tr><tr><td>Silicate coefficient of thermal expansion</td><td>α<sub>Si</sub></td><td>5 × 10<sup>-5</sup> K<sup>-1</sup></td></tr><tr><td>Iron heat capacity</td><td>C<sub>p,Fe</sub></td><td>800 J K<sup>-1</sup> kg<sup>-1</sup></td></tr><tr><td>Silicate heat capacity</td><td>C<sub>p,Si</sub></td><td>1000 J K<sup>-1</sup> kg<sup>-1</sup></td></tr><tr><td>Average specific heat of the target body</td><td>ρ<sub>0</sub>C<sub>p</sub></td><td>4 × 10<sup>6</sup> J K<sup>-1</sup> m<sup>-3</sup></td></tr><tr><td>Buoyancy ratio</td><td>Γ</td><td>94</td></tr><tr><td>Iron latent heat</td><td>L<sub>Fe</sub></td><td>2.7 × 10<sup>5</sup> J kg<sup>-1</sup></td></tr><tr><td>Silicates latent heat</td><td>L<sub>Si</sub></td><td>4 × 10<sup>5</sup> J kg<sup>-1</sup></td></tr><tr><td>Pre-impact mantle temperature</td><td>T<sub>0</sub></td><td>1600 K</td></tr><tr><td>Pre-impact core temperature</td><td>T<sub>c</sub></td><td>2000 K</td></tr><tr><td>Heat diffusivity</td><td>χ</td><td>10<sup>-6</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Thermal conductivity</td><td>k = k<sub>Si</sub> = k<sub>Fe</sub></td><td>4 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Metal content</td><td>f<sub>0</sub></td><td>12.5%</td></tr><tr><td>Reference viscosity</td><td>η<sub>0</sub></td><td>10<sup>20</sup>–10<sup>22</sup> Pa s</td></tr><tr><td>Viscosity ratio between the hot metallic diapir and the surrounding mantle</td><td>λ<sub>0</sub></td><td>10<sup>-2</sup>–1</td></tr><tr><td>Impact energy conversion coefficient</td><td>γ</td><td>0.3</td></tr><tr><td>Volume effectively heated by impact over the impactor volume</td><td>h<sub>m</sub></td><td>2.7</td></tr><tr><td>Surface gravity</td><td>g<sub>0</sub></td><td>3.85 m s<sup>-2</sup></td></tr><tr><td>Core surface gravity</td><td>g<sub>c</sub></td><td>=g<sub>0</sub></td></tr><tr><td>Gravitational constant</td><td>G</td><td>6.67 × 10<sup>-11</sup> m<sup>3</sup> kg<sup>-1</sup> s<sup>-2</sup></td></tr><tr><td>Average magnetic field strength</td><td>B̄</td><td>2.5 mT</td></tr><tr><td>Magnetic diffusivity</td><td>ν</td><td>2 m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Magnetic permeability</td><td>μ</td><td>4π × 10<sup>-7</sup> H m<sup>-1</sup></td></tr><tr><td>Stokes regime geometrical constant</td><td>a<sub>0</sub></td><td>4/15–1/3</td></tr><tr><td>"Early time" regime geometrical constant</td><td>a<sub>1</sub></td><td>1/9</td></tr><tr><td>"Late time" regime parameter</td><td>a<sub>2</sub></td><td>1.3 × 10<sup>-9</sup> m s<sup>1/3</sup></td></tr></table>

cur in both the mantle and the diapir (Section 4.4.3) (Samuel and Tackley, 2008; Monteux et al., 2009).

### 4.1. Physical model: Descent of a single diapir

To investigate the dynamics of the metallic diapir, we adapt the numerical finite volume model in spherical axisymmetric geometry of Monteux et al. (2009). We assume that the mantle is deforming in a diffusion creep limit. Conservation of energy applied to a planet of radius R then leads to

$$ \frac{DT}{Dt} = \frac{\nabla^2 T}{Ra_\chi} + Di \left( \frac{\eta}{\eta_0} \Gamma \Omega - v_r \left( T + \frac{T_0}{\Delta T_0} \right) r \right), \quad (6) $$

where T, t and r are dimensionless temperature, time and radius. $v_r$ is the dimensionless radial velocity. $Ra_\chi = \frac{\Delta\rho_0 g_0 R^3}{\eta_0 \eta_0}$ is the compositional Rayleigh number, $Di = \frac{\alpha \rho_0 g_0 R}{\rho_0 C_p}$ is the dissipation number and $\Gamma = \Delta\rho_0/(\rho_0 \alpha \Delta T_0)$ is the buoyancy ratio (Table 1). Within the target planetary mantle gravity is set equal to a constant: $g_0 = \frac{4}{3} G \pi \rho_0 R$. Within the planetary core, however, gravity is proportional to r, where $g(r) = g_0 \frac{r}{R_c}$. Last, $\Omega$ is the dimensionless dissipation function and expresses the conversion of potential energy into heat through viscous dissipation:

$$ \Omega = 2\varepsilon : \varepsilon. \quad (7) $$

where $\varepsilon$ is the dimensionless shear strain rate tensor. Prior to impact, we assume a homogenous temperature, $T_0$, in the martian mantle. The viscosity is $\eta = \eta_0 \lambda_0^T$, where $\lambda_0$ is the ratio (<1) of the viscosity of the hottest material (iron or silicates) to that of the colder, surrounding solid silicate mantle far from the impact site

where $T = T_0$ (see value in Table 1). This viscosity decreases sharply with temperature and its expression is simpler to implement than the usual Arrhenius law (Ratcliff et al., 1997; Ziethe and Spohn, 2007). We neglect the compositional dependence of the viscosity as this effect will be much smaller than the temperature-dependence. $\eta_0$ is the reference viscosity of mantle material far from the impact site (see value in Table 1). The viscosity contrast between molten iron and mantle under martian conditions may exceed 10–20 orders of magnitude, depending on the mantle temperature. Such viscosity variations are impossible to resolve with our numerical method. However, experiments show that the dynamic influence of these large viscosity variations on diapir sinking are similar to viscosity variations of order $10^2$ (Jellinek et al., 2003; Thayalan et al., 2006), the effects of which we can explore parametrically.

The other dimensionless governing equations are continuity

$$ \mathbf{\nabla} \cdot \mathbf{v} = 0, \quad (8) $$

and momentum conservation, assuming infinite Prandtl number

$$ -\nabla P + \mathbf{\nabla} \cdot \left( \frac{\eta}{\eta_0} [\mathbf{\nabla} \mathbf{v} + [\mathbf{\nabla} \mathbf{v}]^T] \right) + \left( \frac{T}{\bar{T}} - f \right) \mathbf{r} \mathbf{e}_r = 0, \quad (9) $$

where $\mathbf{v}$ and $P$ are the non-dimensional velocity and pressure and $\mathbf{e}_r$ is the radial unit vector. Coupling core and mantle convection models is a difficult problem, given the different timescales and material properties for each layer. We assume here that Eq. (9) is also valid in the core. The core and the mantle are both treated as highly viscous materials. In reality, the viscosity of the core is much smaller than the mantle viscosity and the assumption that Prandtl number is
infinite is no longer valid. We will return to this issue in the conclusion.

The buoyancy force that drives the flow of the diapir towards the center of the protoplanet increases with the metallic volume fraction *f* that varies between 0 (pure silicates) and 1 (pure metal). The metal volume fraction *f* is then simply advected by the flow:

$$ \frac{Df}{Dt} = 0. \qquad (10) $$

## 4.2. Numerical model

We implement a finite volume numerical model to solve Eqs. (6), (8), (9) and (10) in axisymmetric spherical geometry. We use a stream function formulation for the equations of motion with a direct implicit inversion method (Schubert et al., 2001). Eqs. (6) and (10) are solved by an Alternating Direction Implicit (ADI) scheme (Peaceman and Rachford, 1955; Douglas, 1955). The stream function, temperature and compositional fields are described by a second-order approximation in space. To limit numerical diffusion when solving the transport equations, especially for the compositional field, we use a Total Variation Diminishing Superbee scheme (Roe, 1986; Laney, 1998) implemented in an implicit way (Srámek et al., 2010) which enables high resolution of pure advective fields. To avoid any singularity at the center of the planet, we use a staggered grid mesh (see Patankar (1980) for details) where the steam function is calculated at the corners and the temperature and compositional fields are calculated at the center of the grid cells. We use at least 200 grid points along the r-direction (i.e. the resolution *dr* = 17 km) and 400 grid points along the *θ* direction (i.e. the resolution *dθ* = 0.45°). Hence the grid-point density is slightly larger towards the center of the planet. Velocity boundary conditions are free-slip at the surface and along the symmetry axis. Thermal boundary conditions are isothermal at the surface and insulating along the symmetry axis.

## 4.3. Analytical model of the core merging process

Buoyancy-driven interactions between two deformable viscous drops has been widely studied both experimentally and numerically (Davis et al., 1989; Yiantsios and Davis, 1990; Manga and Stone, 1993). Building on this body of work, we characterize the

interaction between the preexisting core and the sinking diapir by monitoring the gap thickness, *h*, which is the thickness of the mantle layer between the base of the diapir and the CMB, measured along the axis of symmetry of the diapir (Fig. 2). We compare the results from our numerical models with analytical predictions. We consider the initial time $t_0$ to be the time at which the impact occurs.

During sinking, a diapir will deform the surrounding mantle over a radial distance that depends on the diameter of the diapir, the effective viscosities of the mantle and core, the gap thickness and the density difference across the CMB. When the diapir is small in comparison with the mantle depth and far from the core-mantle boundary, interactions with the core are negligible (Fig. 2, left). If we also neglect deformation of the density interface forming the upper boundary of the mantle the diapir sinks in the “Stokes regime” with a velocity given by Hadamard (1911) and Rybczynski (1911).

$$ v = - \frac{dh}{dt} = \frac{a_0 \Delta \rho_0 g_0 R_{Fe}^2}{\lambda_0 \eta_0}, \qquad (11) $$

where $a_0$ is a geometrical constant. In the case of a diapir sinking in an infinite medium, $a_0$ is given by the Hadamard-Rybczynski equation and varies from 4/15 = 0.27 (isoviscous) to 1/3 = 0.33 for an inviscid sphere (Hadamard, 1911; Rybczynski, 1911).

Integrating Eq. (11) with $h_0$ as the initial gap thickness,

$$ h(t) = h_0 \left( 1 - \frac{t}{t_s} \frac{R_{Fe}}{h_0} \right). \qquad (12) $$

Here $t_s$ is the characteristic time for diapir descent in this regime and is given by

$$ t_s = \frac{R_{Fe}}{v} = \frac{\lambda_0 \eta_0}{a_0 \Delta \rho_0 g_0 R_{Fe}} \qquad (13) $$

For the impactor size range studied here, the diapir is not necessarily small in comparison to the mantle and the Stokes regime described above is not valid. As the diapir descends and $h \sim R_{Fe}$, mantle flow interacts increasingly with the CMB (Fig. 2, middle). The sinking velocity in this “Early time” regime is modified from Eq. (11) and has the form (Yiantsios and Davis (1990)):

$$ v(t) = - \frac{dh}{dt} = \frac{a_1 \Delta \rho_0 g_0 R_{Fe}^{3/2} h(t)^{1/2}}{\lambda_0 \eta_0} \qquad (14) $$

<figure><img src="image_179.png" /><figcaption>Fig. 2. Dynamics of core merging for the three regimes described in Section 4.3 and computed for $R_{Fe}$ = 300 km sinking in an isoviscous differentiated Mars size planet (silicate material is represented in green and metallic material is represented in red). These results are obtained from the numerical model described in the article. Solid black lines illustrate the streamlines resulting from the diapir sinking. In the left figure, the planetary core is not present to illustrate its negligible influence on the sinking dynamics when the diapir is far from it. In this case, the viscous stress scales with $1/R_{Fe}$ (Hadamard, 1911; Rybczynski, 1911). In the middle figure, the interaction with the core results in viscous stresses that scale with $1/h$. In the right figure where $h < R_{Fe}$, interactions between the diapir and the merging core increase and the lubrication force scales with $1/h$ (Yiantsios and Davis, 1990). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
where $a_1$ is a geometrical constant that we determine from our numerical models. Integration of Eq. (14) with the same initial condition leads to

$$h(t) = h_0 \left(1 - \frac{t}{t_c}\right)^2 \qquad (15)$$

and identifies a characteristic time for diapir descent that is governed primarily by the mantle viscosity (Yiantsios and Davis, 1990):

$$t_c = \frac{2\lambda_0\eta_0}{a_1\Delta\rho_0g_0R_{Fe}} \left(\frac{h_0}{R_{Fe}}\right)^{1/2} \qquad (16)$$

In contrast to Eq. (13), this sinking time depends on the ratio $(R_{Fe}/h_0)$. When $h \to 0$, $v \to 0$ (Eq. (14)) and the merging time increases towards infinity as a result of the large lateral pressure gradients required to drive viscous mantle out of the way of the sinking diapir. Indeed, Eq. (15) is strictly valid only where $t \le t_c$. At “Late time” ($t/t_c \to \infty$) the gap becomes very thin and is expected to evolve with a form (Jones and Wilson (1978) and Yiantsios and Davis (1990)):

$$h(t) = a_2 t^{-1/3} \qquad (17)$$

where $a_2$ is a function of $\lambda_0$ and the density difference across the CMB. We obtain $a_2$ from our numerical results. Comparison of Eqs. (12), (15) and (17) shows that as the diapir approaches and deforms the CMB, the rate of descent will steadily decline.

## 4.4. Numerical results

For the range of plausible martian impactors, $R_{Fe}$ ranges from 100 km to about 650 km. In this range, the smallest values of $R_{Fe}$ lead to an initial Stokes regime, after which the diapir motion or gap evolution is successively described by the Early Time and Late Time regimes. Intermediate and large values of $R_{Fe}$ lead to an initial state described by the Early Time regime, after which the motion of the diapir is described by Late Time regime. As the diapir size increases, the relative amount of time spent in the Late Time regime increases. Because the motion of the diapir in the Stokes's regime is very fast, we consider an initial gap thickness $h_0$ that is independent of the size of the impactor, and we model only the Early Time and Late Time regimes. We set $h_0 = 580$ km, which allows us to consider $R_{Fe}$ in the range 100 km to a maximum value of $(R_{Mars} - R_c - h_0)/2 = 560$ km. Finally, we note that for small diapirs

$(R_{Fe} < 100$ km), the Stokes regime dominates during the sinking, and conversely for very large diapirs ($R_{Fe} > 650$ km) the Late Time Regime may be the only sinking regime. Neither of these two cases is considered here.

### 4.4.1. Initial setup

We use the numerical model to characterize the dynamics of core merging after a giant impact on a Mars size body. First, we consider the initial martian core temperature to be homogeneous (i.e. the thickness of the thermal boundary layer $\delta = 0$) and we use the estimated present-day temperature value $T = T_{core} = 2000$ K (Williams and Nimmo, 2004). The initial thermo-chemical conditions are shown schematically in Fig. 3 in which the hatched regions represent the metallic phase from the impactor and from Mars. The initial temperature is represented with a gray scale.

### 4.4.2. Evolution of the gap thickness

We first monitor the thermal evolution following an impact in the absence of core merging (i.e. we consider only the thermal re-equilibration since no core merging occurs). This unphysical case is used as a control case to separate the influence of the impact heating from the diapir merging on the thermal evolution of the planet. Fig. 4 shows the thermal evolution after a 600 km impact on a Mars size body with homogeneous viscosity ($\eta_0 = 10^{22}$ Pa s). After the impact heating, the hot anomaly advects and spreads beneath the planetary surface in a few million years. Then, the hot anomaly and the target planet's core cool by diffusion.

We now consider the impactor's core. Fig. 5 shows the thermal and chemical evolution after a 600 km impact (i.e. $R_{Fe} = 300$ km) on a differentiated Mars with homogeneous viscosity ($\eta_0 = 10^{22}$ Pa s). Following the impact, a metallic diapir from the impactor forms rapidly and sinks towards the pre-existing core, while undergoing significant viscous heating. Motion of the diapir displaces surrounding mantle and deforms the CMB (Figs. 2 and 5). Merging and thermal mixing of the hot diapir into the core occurs eventually, causing the average temperature of the core to rise (Fig. 5), leading in turn, to an enhanced CMB heat flux that causes the formation and rise of a plume (Fig. 5d). The plume rising from the CMB is not the result of the onset of thermal convection in the mantle to remove the heat from the core. It is a consequence of the diapir sinking and especially of the viscous dissipation that occurs along the sinking path and generates a localized hot channel.

<figure><img src="image_180.png" /><figcaption>**Fig. 3.** Sketch to show repartitioning of the post impact heating in the martian mantle. The solid black line represents the planetary surface. Melting is restricted to the material limited by the dashed line. The impactor's core is represented by hatched lines.</figcaption></figure>

<figure><img src="image_181.png" /><figcaption>**Fig. 4.** Non-dimensional temperature evolution as a function of time. In this control case, we do not consider the impactor's core (computed for a uniform viscosity with $R = 3400$ km, $R_{imp} = 600$ km, $\eta_0 = 10^{22}$ Pa s and 300 × 600 grid points). To convert this temperature field in Kelvin, $T(K) = T_0 + T\Delta T_0$.</figcaption></figure>
<figure><img src="image_182.png" /><figcaption>**Fig. 5.** Non-dimensional temperature (left) and composition (right) at times $t = 0$ (a), $t = 0.6$ Myr (b), $t = 6$ Myr (c) and $t = 38.8$ Myr (d) (computed for a uniform viscosity with $R = 3400$ km, $R_{\text{imp}} = 600$ km, $R_{\text{Fe}} = 300$ km, $\eta_0 = 10^{22}$ Pa s, $h_0 = 580$ km and $300 \times 600$ grid points). As the diapir sinks, its velocity decreases because of the increasing influence of the pre-existing core on the dynamics (see Section 4.3 for details). After the merging (c), the hot metallic fraction from the impactor spreads at the top of the martian core (d). To convert this temperature field in Kelvin, $T(K) = T_0 + T\Delta T_0$.</figcaption></figure>

The plume rising does not happen in the absence of the diapir sinking (Fig. 4).

In Fig. 6, we show the temporal evolution of the gap thickness after a 600 km impact and compare it with theoretical predictions from Eqs. (15) and (17) to obtain characteristic coalescence times $t_c$.

<figure><img src="image_183.png" /><figcaption>**Fig. 6.** Evolution of the gap thickness between the merging core and the CMB with $R_{\text{imp}} = 600$ km (i.e. $R_{\text{Fe}} = R_{\text{imp}}/2 = 300$ km). Theoretical evolutions for early times from Eq. (15) and late times from Eq. (17) are shown with red and green dashed lines respectively. From this model, we obtain a coalescence time, $t_c = 1.55$ Myr and $a_2 \sim 1.3 \times 10^9$ m s$^{1/3}$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

and values for $a_2$. We repeat this for a range of plausible impactor sizes. Consistent with Eq. (16), our numerical values for $t_c$ decrease in proportion to $R_{\text{Fe}}^{3/2}$ (red circles, Fig. 7). In Fig. 7, we fit the numerical $t_c$ values with the analytical prediction (black dashed line, Fig. 7) from Eq. (16) and obtain $a_1 = 1/9$.

### 4.4.3. Influence of the temperature-dependence of the mantle viscosity on the core merging timescale

The viscosity contrast between the sinking diapir and its environment is a key parameter in understanding the core merging dynamics. Qualitatively, as the metallic diapir sinks, shear heating occurs at the interface of the diapir (Samuel et al., 2010) and both the temperature of the diapir and the surrounding mantle can increase. Depending on the effective viscosity contrast between the metallic diapir and the mantle, the mean temperature increase within the sinking diapir can reach a maximum value $\approx 2\Delta T_0$ (Monteux et al., 2009). This temperature excess drives a heat flux from the diapir to the surrounding mantle, leading to the progressive growth of a thermal boundary layer (Fig. 5). Because the mantle viscosity declines exponentially with increasing temperature, mantle shear will be increasingly concentrated within the hottest, low viscosity part of this thermal boundary layer (Morris, 1982; Thayalan et al., 2006). Thus, heat transfer from the diapir will reduce both the effective viscosity retarding diapir motion and the radial length scale over which viscous deformation in the mantle occurs. The timescale for merging will consequently decline as a result of two processes. First, the lubricating effect of an enveloping layer of low mantle viscosity will cause metallic diapirs to descend more rapidly (cf. Eq. (14)). Second, localization of shear within only the narrow, hottest part of the thermal boundary layer will enable diapirs to more closely approach the CMB before interacting with the core.

We explore these effects quantitatively in Fig. 8, which shows analytical and numerical solutions for the evolution of the gap thickness $h$ for a range of background mantle viscosities and diapir/mantle viscosity ratios. In the isoviscous case ($\lambda_0 = 1$), when the average viscosity decreases in the mantle (from $\eta_0 = 10^{22}$ Pa s, black solid line to $\eta = 10^{20}$ Pa s, gray solid line), the diapir sinks faster (see Eq. (14)). Consequently, the sinking time values obtained by fitting our numerical results with Eq. (15) decrease by two orders of magnitude (see Table 2).

Theoretical, experimental and numerical models (Morris, 1982; Jellinek et al., 2003; Thayalan et al., 2006) show that whereas in the isoviscous case radial mantle deformation scales with the size of the diapir, $R_{\text{Fe}}$, in the temperature-dependent case deformation is confined to the hottest, low viscosity part of the thermal boundary layer with a length scale typically $\sim 0.1 R_{\text{Fe}}$. As a consequence, the value of $a_1$ obtained by fitting our numerical results with Eq. (16) increases (see Table 2) and the sinking time value for ($\eta_0 = 10^{22}$ Pa s, $\lambda_0 = 0.01$) is smaller than for ($\eta_0 = 10^{20}$, $\lambda_0 = 1$).

## 5. Effect on a preexisting dynamo

We now discuss the consequences of a giant-impact and the subsequent core merging process on the martian dynamo. We first revisit the core heat flux conditions required for dynamo generation. Next, we explore the influence of core merging on the efficiency of core cooling, and thus dynamo action, for different pre-impact mantle thermal regimes.

### 5.1. Conditions required for dynamo generation

As discussed in Monteux et al. (2011), three conditions are required to generate a dynamo on a growing planet:
<figure><img src="image_184.png" /><figcaption>**Fig. 7.** Characteristic coalescence time, $t_c$ (left panel) as a function of the impactor core radius $R_{Fe}$. Results from numerical experiments (with uniform viscosity $\eta_0 = 10^{22}$ Pa s and $h_0 = 580$ km) are represented with red circles. Theoretical fit from Eq. (16) is shown with the black dashed line ($a_1 = 1/9$). The corresponding errors between the numerical and theoretical models are plotted in the right panel. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_185.png" /><figcaption>**Fig. 8.** Time evolution of the gap thickness between the sinking core and the pre-existing core for four different rheologies. Solid black line and solid gray line represent results from uniform viscosity models with $\eta_0 = 10^{22}$ Pa s and $\eta_0 = 10^{20}$ Pa s respectively. Solid red line and solid green line represent results from temperature dependent viscosities with $\lambda_0 = 0.1$ and $\lambda_0 = 0.01$ respectively (with $\eta = \eta_0 \lambda_0^T$ and $\eta_0 = 10^{22}$ Pa s). Theoretical predictions are shown in the corresponding color with dashed lines for the coalescence theory (Eq. (15)). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 2**
Sinking time values obtained fitting numerical experiments with theoretical predictions (Eq. (15)) for different values of $\eta_0$ and $\lambda_0$ (with $R_{imp} = 600$ km, $h_0 = 580$ km).</caption><thead><tr><th></th><th>η<sub>0</sub> = 10<sup>22</sup> Pa s<br>λ<sub>0</sub> = 1</th><th>η<sub>0</sub> = 10<sup>20</sup> Pa s<br>λ<sub>0</sub> = 1</th><th>η<sub>0</sub> = 10<sup>22</sup> Pa s<br>λ<sub>0</sub> = 0.1</th><th>η<sub>0</sub> = 10<sup>22</sup> Pa s<br>λ<sub>0</sub> = 0.01</th></tr></thead><tbody><tr><td>t<sub>c</sub></td><td>1.56 × 10<sup>3</sup> kyr</td><td>15.6 kyr</td><td>79 kyr</td><td>7.5 kyr</td></tr><tr><td>a<sub>1</sub></td><td>0.11</td><td>0.11</td><td>0.21</td><td>0.23</td></tr></tbody></table>

1. The mean heat flux $q_{CMB}$ out of the core must exceed the adiabatic value $q_A$ such that convection can occur (Stevenson et al., 1983). This condition is

$$q_{CMB} > q_A = \frac{k_{Fe} \alpha_{Fe} g_c T_c}{C_{p,Fe}} \qquad (18)$$

where $k_{Fe}$, $\alpha_{Fe}$, $T_c$ and $C_{p,Fe}$ are respectively the thermal conductivity, the coefficient of thermal expansion, the temperature and the heat capacity of the metallic phase and $g_c$ is the gravity at the surface of the protocore. The largest uncertainties in the

calculation of $q_A$ are related to $T_c$ and $k_{Fe}$. Assuming reasonable values, convection in the martian core should start once $q_{CMB} > q_A = 5 - 19$ mW m$^{-2}$ (Nimmo and Stevenson, 2000) (see Table 1 for values).

2. The ratio of the rate at which gravitational potential energy is released by convection to the rate of ohmic dissipation, $\Phi$, must exceed a critical value (Buffett, 2002):

$$\frac{4\pi R_c^2 q_{CMB}}{\Phi} > \frac{1}{\epsilon_T} \qquad (19)$$

where $\epsilon_T$ is the Carnot-style efficiency for thermal convection. Here, we do not consider the effect of compositional convection or the presence of an inner core. Assuming that the characteristic scale of core flow leading to magnetic field generation is the radius of the protocore, $\Phi$ can be approximated as (Buffett, 2002):

$$\Phi = \frac{4}{3} \pi R_c \left( \frac{v \bar{B}^2}{\mu} \right), \qquad (20)$$

where $\bar{B}$ is the average strength of the magnetic field within the core, $v$ is the magnetic diffusivity and $\mu$ is the magnetic permeability (see Table 1 for values). In the absence of constraints on $\bar{B}$ for early planets, we assume a current Earth-like value of 2.5 mT (Kuang and Bloxham, 1997) that is independent of the protocore size or the planetary radius. The efficiency of thermal convection is given by Buffett (2002)

$$\epsilon_T = \frac{0.8\pi}{3} \frac{\alpha_{Fe} G \rho_{Fe} R_c^2}{C_{p,Fe}} \left( 1 - \frac{q_A}{q_{CMB}} \right), \qquad (21)$$

with $\rho_{Fe}$ the density of the metallic phase and $G$ the gravitational constant.

Hence Eqs. (19)–(21) along with values from Table 1 lead to:

$$q_{CMB} - q_A > \frac{\sqrt{B^2 C_{p,Fe}}}{0.8\pi G \mu \alpha_{Fe} \rho_{Fe} R_c^3} \approx 8 \times 10^{-2} \text{ mW m}^{-2} \qquad (22)$$

3. The structure of the convective motions carrying magnetic field lines must be sufficiently complicated to favor self-sustaining dynamo action. A measure of this complexity is that the magnetic Reynolds number (Christensen and Aubert, 2006) must exceed a critical value:

$$Re_m = \frac{UL}{v} > Re_m^{crit} = O(10 - 10^2) \qquad (23)$$
Here, L and U are the characteristic length and velocity scales for the flow within the protocore and v is the magnetic diffusivity of the metal phase. The natural length scale in the problem is the depth of the convecting iron layer, however the velocity scale depends on the leading order force balance (Christensen, 2010). As the rotation rate of growing planets is potentially time-dependent and poorly constrained, a convenient and reasonable choice is based on a balance between inertial and buoyancy forces and is (Stevenson, 2003):

$$U \sim \left( \frac{(q_{CMB} - q_A) \alpha_{Fe} g_c R_c}{\rho_{Fe} C_{p,Fe}} \right)^{1/3}, \quad (24)$$

Taking $L = R_c$ and $Re_m^{crit} = 10^2$, the combination of Eqs. (23) and (24) leads to the condition

$$q_{CMB} - q_A > \frac{\rho_{Fe} C_{p,Fe}}{\alpha_{Fe} g_c R_c} \left( Re_m^{crit} \frac{v}{R_c} \right)^3 \approx 10^{-4} \text{ mW m}^{-2} \quad (25)$$

Among the three criteria above, the first is typically considered a necessary condition for a thermally-driven dynamo (subadiabatic dynamos are possible if compositional buoyancy effects enter the problem). It is not a sufficient condition for dynamo action, as indicated by the other two conditions. However, conditions 2 and 3 require that the heat flux $q_{CMB}$ at the CMB be only about 1% greater than the adiabatic heat flux $q_A$. Thus, we take $q_A$ to be a reasonable critical value for dynamo generation (Stevenson et al., 1983; Nimmo and Stevenson, 2000; Williams and Nimmo, 2004).

## 5.2. Thermal consequences of the core merging

### 5.2.1. CMB cooling

To understand the influence of the merging process on thermally-driven dynamo action, we first specify the thermal structure and cooling regime of the core prior to merging, which depends on the previous accretion and core formation histories (Monteux et al., 2009; Samuel et al., 2010). Two end-member models for the initial interior thermal state of Mars are the isothermal case with the core and the mantle having the same temperature, and a situation where the core is much hotter than the mantle. In the case of a pre-impact core hotter than the mantle, the thickness $\delta$ of the thermal boundary layer (TBL), across which temperature declines linearly from the core to the mantle value is governed by the mode of heat loss to the mantle. Strong cooling to the mantle leads to thin thermal boundary layers and vigorous convection, weak cooling leads to relatively more sluggish flow or no convection at all. To simplify this problem we consider three pre-impact core thermal states: an isothermal core and mantle (uniform temperature), and two models where the core is hotter than the mantle with $\delta = 0.1 R_c$ (strong core cooling) and $\delta = 0.4 R_c$ (weak core cooling and no core convection).

Our results show that one major consequence of large impacts is to heat the mantle and reduce the temperature difference between the core and mantle (where there is one initially) and, thus the CMB heat flux (Roberts et al., 2009). Where a metallic diapir forms and descends such that it becomes hotter than the pre-existing core, merging can lead to the spread of a thin, hot, gravity current along the CMB. If the timescale for spreading $t_{spread}$ is much smaller than the timescale for thermal diffusion $t_{diff}$ across the current then this stratification can insulate the core against the mantle. This can in turn drive a heat flux into the core itself and inhibit dynamo action (Arkani-Hamed and Olson, 2010).

In the uniform temperature case, following merging and coalescence (Fig. 9, top-left panel, solid circles) an increase in the CMB heat flux occurs at $\sim 10t_c$ because the production of thermal stratification increases the temperature difference to the mantle. The

magnitude of this temperature difference increases with rate of spreading relative to the rate at which thermal diffusion smoothes temperature variations. Thus, the largest diapirs spreading in the lowest viscosity cores are characterized by the greatest increases in CMB heat flux (green curve). By contrast, in the weak and strong core cooling cases, initial mantle heating causes a reduction in the CMB heat flux following coalescence that is followed, in turn, by an increase in this heat flux until $t/t_{diff} \to 1$. The magnitude of this increase is in proportion to the strength of the resulting thermal stratification. As in the uniform temperature case, the greatest CMB heat fluxes then correspond with the largest diapirs and lowest viscosity cores (green curves). Notably, the smallest and most slowly spreading diapir in the strong cooling regime (Fig. 9, red curve) produces no significant stratification because $t_{spread}/t_{diff} \gg 1$.

### 5.2.2. Merging, stratification and thermal mixing: Potential effect on an early Mars dynamo

Material from our lowest viscosity and largest diapirs spreads along the CMB as a thin, axisymmetric gravity current. In the limit that this spreading is under the zero Reynolds number conditions of our numerical simulations, the radius of the current will increase at a rate that is approximately proportional to time $t^{1/8}$ (Huppert, 1982) to form a hot, stratified layer that insulates the core and inhibits core convection and dynamo action (Fig. 9, green curves) (Arkani-Hamed and Olson, 2010). This thermal stratification persists until the thermal anomaly introduced by the gravity current has diffused into the overlying mantle and underlying core. Where weak or no stratification is produced, merging and spreading weakens core convection as a result of mantle heating but will not necessarily destroy a pre-existing dynamo (Fig. 9, blue and red curves).

In reality, the gravity current driven as a result of the merging of a $\sim 400$ km radius diapir carrying $\sim 800$ K temperature variations can spread in a very low viscosity core under finite-or even high-Re conditions. In the high-Re limit, at large times where the current is thin in comparison to its radius, the radius of the current will increase in proportion to $t^{1/2}$ (Kotsovinos, 2000) and may involve extensive mixing across the density interface forming the base of the flow, depending on the magnitude of the temperature difference, the thickness of the flow and the rotation rate of Mars (Linden, 1977; Linden, 1979; Fleury et al., 1991; Kotsovinos, 2000).

A careful consideration of the efficiency of thermal mixing across this density interface requires a different numerical model and is beyond the scope of this paper. However, it is useful to explore the end-member case of perfect thermal mixing in order to determine the maximum core heat flux that is likely to occur after merging, spreading and thorough mixing. To obtain a crude estimate of the effects of complete thermal mixing of the hot gravity current into the underlying core, we stir in the gravity current at $t = t_{spread}$ and then determine the mean temperature of the core $\bar{T}_c$ and the corresponding CMB heat flux (Fig. 10, solid square). Because spreading and mixing in a high-Re limit will occur on a timescale that is much smaller than $t_{diff}$ (i.e., for $t_{spread} \ll t_{diff}$), we explore only the effect of stirring in the stratified layer corresponding to the largest, lowest viscosity diapir in Fig. 9 (Fig. 10). In Fig. 10, the discontinuity in the heat flux evolution is the consequence of imposing a new thermal state. The main effect of an imposed mixing event at $t = t_{spread}$ is to reduce the average core temperature to a value intermediate between the initial and fully stratified cases (Fig. 10). Whereas the peak in the CMB heat flux corresponds to stratification that would likely destroy the dynamo (cf., Fig. 9), the lower, post-mixing, heat flux can drive a dynamo for several diffusion times, which corresponds to a few hundred million years.
<figure><img src="image_186.png" /><figcaption>**Fig. 9.** Time evolution of the mean CMB heat flux in the stratified case as a function of the viscosity and impactor size. The characteristic diffusion time $t_{diff} = d^2/\kappa$, where the steady-state thickness $d$ is defined where the spreading gravity current extends over 25% of the CMB. In the three panels different initial thermal states before impact are considered. The two horizontal black dashed lines represent the theoretical dynamo criterion for $q_{CMB}$ ($q_{CMB} > q_A = 5 - 19$ mW m$^{-2}$). In the top-left panel, pre-impact temperature is assumed to be constant ($T = T_0$). In the bottom panels, the core is initially hotter than the mantle ($T = T_{core} > T_0$). The black dotted lines represent $q_{CMB}$ in the absence of an impact. In the bottom left panel the CMB thermal boundary layer thickness $\delta = 0.4R_{Fe}$ and in the bottom right panel $\delta = 0.1R_{Fe}$. Red, blue and green solid lines represent the evolution with ($R_{imp} = 500$ km, $\eta_0 = 10^{22}$ Pa s), ($R_{imp} = 800$ km, $\eta_0 = 10^{22}$ Pa s), and ($R_{imp} = 800$ km, $\eta_0 = 10^{20}$ Pa s) respectively. In the top-left panel, the blue dashed line is the result of the same model as the solid blue line but with $v_{imp} = 2v_{esc}$ (i.e. the post-impact temperature increase is $\approx 4$ times larger). Filled circles represent the coalescence times $t_c/t_{diff}$ from Eq. (16). In the weak and strong core cooling cases, the mean heat flux decreases for $t_c/t_{diff} < t/t_{diff} < 10t_c/t_{diff}$ because of the post-impact mantle thermal anomaly that inhibits the CMB heat loss. The combination of the spreading and cooling mantle thermal anomaly, and of the thermal core stratification leads to an increase of the CMB heat flux after $t \sim 10t_c/t_{diff}$. The maximum average heat flux increase at the CMB is reached when $R_{imp} = 800$ km and when $\eta_0 = 10^{20}$ Pa s. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_187.png" /><figcaption>**Fig. 10.** Time evolution of the mean CMB heat flux in the perfect mixing case as a function of the viscosity and impactor size. The characteristic diffusion time $t_{diff} = d^2/\kappa$, where the steady-state thickness $d$ is defined where the spreading gravity current extends over 25% of the CMB. Here $t_{diff} = 27.8$ Myr. Perfect thermal mixing is imposed within the core when $t = t_{spread}$ (vertical dashed line). The two horizontal black dashed lines represent the theoretical dynamo criterion for $q_{CMB}$ ($q_{CMB} > q_A = 5 - 19$ mW m$^{-2}$). The green solid line represent the evolution with ($R_{imp} = 800$ km, $\eta_0 = 10^{20}$ Pa s, $t_{spread} = 4$ Myr). The filled circle represents the coalescence times ($t_c/t_{diff}$) from Eq. (16), the filled square represents the spreading time ($t_{spread}/t_{diff}$) and the filled triangle represents the relevant mean heat flux after thermal mixing. The right panel shows the corresponding non-dimensional temperature field right after the perfect mixing. To convert this temperature field in Kelvin, $T(K) = T_0 + T\Delta T_0$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
## 6. Conclusions

Giant impacts may have played an important role in the thermal and magnetic evolutions of terrestrial planets. The giant impact potentially forming the martian dichotomy could have supplied a significant volume of iron that merged with the preexisting core. We model the dynamics of the sinking of a large metallic diapir within a silicate mantle and the influence of the deformation of the preexisting core on the merging process. We obtain characteristic timescales of the coalescence time between the two cores as a function of the viscosity of the mantle and the size of the sinking diapir. Our results are in agreement with theoretical and experimental models of buoyancy-driven motions sinking toward a deformable surface. We also show that the thermal consequences of the merging process can influence the generation of a martian dynamo. Assuming perfect thermal mixing, the dynamo generation might be favored by large impacts, low viscosity of the preexisting core and strong preexisting core cooling. The dynamo can last as long as several hundred million years depending on the mixing process within the core.

The novelty of our models is that we characterize and understand theoretically the sinking of a large volume of iron through a silicate mantle as well as the interactions between the two cores during the merging. Our thermal models also include the dissipation of energy via viscous heating which allows us to characterize the thermal consequences of the impactor's core sinking within the mantle and the pre-existing core. The results illustrate that the core merging is a process with a characteristic time that is essentially governed by the mechanical properties of the mantle. The dynamics and thermal consequences of the core merging process are important aspects of the giant impacts that potentially delivered large volumes of iron during the late accretionary history of terrestrial planets.

The heating induced by the impact within the mantle and the viscous dissipation during diapir sinking enhances a thermal anomaly within the mantle. These processes strongly affect the heat budget of the mantle and may have long term consequences on a planet's thermal evolution. The impact heating within the mantle leads to a dichotomy of the core cooling: the core heat flux at the CMB below the impact site is different than the CMB heat flux on the opposite side of the core. This mantle heterogeneity might affect the dynamo and the resulting paleomagnetic field recorded at the surface (Stanley et al., 2008; Sreenivasan and Jellinek, 2012). Moreover, recent dynamo simulations show that the efficiency of a mantle heterogeneity centered at the geographical pole in producing a south-north dichotomy is much higher than that of an heterogeneity centered at the equator in producing an east-west dichotomy (Amit et al., 2011). Hence, our models may have interesting implications for the understanding of the structure of the past martian dynamo.

We illustrate the importance of better constraining the thermal states of the cores, the size of the largest metallic diapir that remains after a giant impact, and the depth at which the impactor's core is buried before experiencing viscous deformation and the core merging process. We assume that the impact occurs with the martian escape velocity which is the minimum impact velocity; larger impact velocities can substantially increase the post-impact temperature in the mantle and in the merging core prior to its sinking. Hence, in this respect our results provide a lower bound on the temperature increase in the planetary core after the diapir merging and on the CMB heat flux increase (see Fig. 9, top left panel).

In our models, we make the strong assumption that the sinking impactor's core, the target silicate mantle and the target iron core are in the infinite Prandtl number regime. In reality, both the molten metallic cores have viscosities that are 10–20 orders of

magnitude smaller than the mantle viscosity. Accounting for such large realistic viscosity contrasts may have consequences on both the sinking and the merging dynamics. Indeed, considering a more realistic viscosity contrast between the sinking metallic diapir and the silicate mantle may affect the heat dissipation during the sinking. Viscous heating will preferentially occur in the surrounding mantle reducing the maximal temperature reached within the merging diapir and the energy brought to the preexisting core (Monteux et al., 2009, 2011; Samuel et al., 2010). However, considering an inviscid sinking metallic diapir should not change significantly the diapir sinking velocity which is mainly governed by the mantle viscosity (Hadamard, 1911; Rybczynski, 1911; Yientsios and Davis, 1990). Inertial and Coriolis forces can even favor dynamo generation during the sinking of the molten metallic diapir (Monteux et al., 2012). Concerning the merging dynamics, the injection of a large volume of molten iron in a turbulent and rotating preexisting core is a process that is difficult to handle numerically. To the first order, the merging dynamics is probably bounded between two extreme cases. In the infinite Prandtl number approximation, the stratification at the end of the merging is unavoidable and the core merging process will stop the preexisting dynamo. However, turbulence and rotational forces may enhance efficient thermal mixing within the preexisting core which may reactivate a dynamo previously killed by the impact heating.

Future missions will provide important additional constraints on interior structure and remanent crustal magnetic field of Mars that will help constrain interior thermal evolution and dynamo models. The MAVEN mission will provide low-altitude vector magnetic field observations at locations other than those sampled at low altitudes by Mars Global Surveyor. The InSight mission will provide constraints on the outer core radius, mantle seismic velocity structure (and hence on temperature and composition), and heat flow at a single lander location. The inference (from seismology) of any thermal dichotomy between the northern hemisphere and the southern hemisphere could be the consequence of a large core merging event. Another key to testing such models may consist of measurements of the degree of equilibration between the mantle and the core material of Mars. During the core merging following a giant impact, the metal may sink as unequilibrated blobs (Dahl and Stevenson, 2010; Kleine and Rudge, 2011). Hence, some memory of the core merging process might have been recorded in the martian mantle. An inferred dichotomy in the chemical composition of the martian mantle between the northern lowlands and the southern highlands could also be the consequence of a giant impact and of the subsequent core merging process.

## Acknowledgments

J. Monteux is funded by Agence Nationale de la Recherche (Accretis decision No. ANR-10-PDOC-001-01). M. Jellinek and C.L. Johnson acknowledge support from the Natural Science and Engineering Research Council of Canada. We thank J. Roberts and anonymous reviewer for thoughtful and constructive comments.

## References

* Acuña, M. et al., 1999. Global distribution of crustal magnetization discovered by the Mars Global Surveyor MAG/ER experiment. Science 284, 790–793.

* Agee, C.B., 1997. Melting temperatures of the Allende meteorite: Implications for a Hadean magma ocean. Phys. Earth Planet. Inter. 100, 41–47.

* Amit, H., Christensen, U.R., Langlais, B., 2011. The influence of degree-1 mantle heterogeneity on the past dynamo of Mars. Phys. Earth Planet. Inter. 189, 63–79. http://dx.doi.org/10.1016/j.pepi.2011.07.008.

* Andrews-Hanna, J.C., Zuber, M.T., Banerdt, W.B., 2008. The Borealis basin and the origin of the martian crustal dichotomy. Nature 453, 1212–1215. http://dx.doi.org/10.1038/nature07011.

* Arkani-Hamed, J., Olson, P., 2010. Giant impact stratification of the martian core. Geophys. Res. Lett. 37 (L02), 201. http://dx.doi.org/10.1029/2009GL041417.
Asphaug, E., 2010. Similar-sized collisions and the diversity of planets. Chem. Erde/ Geochim. 70, 199–219. http://dx.doi.org/10.1016/j.chemer.2010.01.004.

Benz, W., Slattery, W.L., Cameron, A.G.W., 1988. Collisional stripping of Mercury's mantle. Icarus 74, 516–528. http://dx.doi.org/10.1016/0019-1035(88)90118-2.

Buffett, B.A., 2002. Estimates of heat flow in the deep mantle based on the power requirements for the geodynamo. Geophys. Res. Lett. 29 (12), 1–4.

Canup, R.M., 2004. Simulations of a late lunar-forming impact. Icarus 168, 433–456. http://dx.doi.org/10.1016/j.icarus.2003.09.028.

Christensen, U.R., 2010. Dynamo scaling laws and applications to the planets. Space Sci. Rev. 152, 565–590.

Christensen, U.R., Aubert, J., 2006. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields. Geophys. J. Int. 166, 97–114.

Citron, R.I., Zhong, S., 2012. Constraints on the formation of the martian crustal dichotomy from remnant crustal magnetism. Phys. Earth Planet. Inter. 212, 55–63. http://dx.doi.org/10.1016/j.pepi.2012.09.008.

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

Ćuk, M., Stewart, S.T., 2012. Making the Moon from a fast-spinning Earth: A giant impact followed by resonant despinning. Science 338, 1047–1052. http://dx.doi.org/10.1126/science.1225542.

Dahl, T.W., Stevenson, D.J., 2010. Turbulent mixing of metal and silicate during planet accretion and interpretation of the Hf-W chronometer. Earth Planet. Sci. Lett. 295, 177–186. http://dx.doi.org/10.1016/j.epsl.2010.03.038.

Davies, G.F., 1982. Ultimate strength of solids and formation of planetary cores. Geophys. Res. Lett. 9, 1267–1270.

Davis, R.H., Schonberg, J.A., Rallison, J.M., 1989. The lubrication force between two viscous drops. Phys. Fluids 1, 77–81. http://dx.doi.org/10.1063/1.857525.

Douglas, J., 1955. On the numerical integration of $\frac{\partial^2 y}{\partial x^2} + \frac{\partial^2 y}{\partial z^2} = \frac{\partial u}{\partial t}$ by implicit methods. J. Soc. Ind. Appl. Math. 3, 42–65. http://dx.doi.org/10.1137/0103004.

Elkins-Tanton, L.T., Parmentier, E.M., Hess, P.C., 2003. Magma ocean fractional crystallization and cumulate overturn in terrestrial planets: Implications for Mars. Meteorit. Planet. Sci. 38, 1753–1771.

Ernst, C.M. et al., 2010. Exposure of spectrally distinct material by impact craters on Mercury: Implications for global stratigraphy. Icarus 209, 210–223. http://dx.doi.org/10.1016/j.icarus.2010.05.022.

Fleury, M., Mory, M., Hopfinger, E.J., Auchure, D., 1991. Effects of rotation on turbulent mixing across a density interface. J. Fluid Mech. 223, 165–191. http://dx.doi.org/10.1017/S0022112091001386.

Frey, H.V., 2006. Impact constraints on the age and origin of the lowlands of Mars. Geophys. Res. Lett. 33, L08S02. http://dx.doi.org/10.1029/2005GL024484.

Gerya, T.V., Yuen, D.A., 2007. Robust characteristics method for modelling multiphase visco-elastic thermo-mechanical problems. Phys. Earth Planet. Inter. 163, 83–105.

Hadamard, J., 1911. Mouvement permanent lent d'une sphère liquide et visqueuse dans un liquide visqueux. C.R. Acad. Sci. 152, 1735–1738.

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

Höink, T., Schmalzl, J., Hansen, U., 2005. Formation of compositional structures by sedimentation in vigorous convection. Phys. Earth Planet. Inter. 153, 11–20.

Hood, L.L., Richmond, N.C., Pierazzo, E., Rochette, P., 2003. Distribution of crustal magnetic fields on Mars: Shock effects of basin-forming impacts. Geophys. Res. Lett. 30 (6), 1–4.

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.

Jellinek, A.M., Gonnermann, H.M., Richards, M.A., 2003. Plume capture by divergent plate motions: Implications for the distribution of hotspots, geochemistry of mid-ocean ridge basalts, and estimates of the heat flux at the core-mantle boundary. Earth Planet. Sci. Lett. 205, 361–378. http://dx.doi.org/10.1016/S0012-821X(02)01070-1.

Johnson, C.L., Phillips, R.J., 2005. Evolution of the Tharsis region of Mars: Insights from magnetic field observations. Earth Planet. Sci. Lett. 230, 241–254.

Jones, A.F., Wilson, S.D.R., 1978. The film drainage problem in droplet coalescence. J. Fluid Mech. 87, 263–288. http://dx.doi.org/10.1017/S0022112078001585.

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

Ke, Y., Solomatov, V.S., 2009. Coupled core–mantle thermal evolution of early Mars. J. Geophys. Res. (Planets) 114 (13), 1–12.

Kleine, T., Rudge, J.F., 2011. Chronometry of meteorites and the formation of the Earth and Moon. Elements 7, 41–46. http://dx.doi.org/10.1016/j.gca.2006.06.004.

Kotsovinos, N.E., 2000. Axisymmetric submerged intrusion in stratified fluid. J. Hydraul. Eng. 126, 446–456. http://dx.doi.org/10.1061/(ASCE)0733-9429(2000)126:6(446).

Kuang, W., Bloxham, J., 1997. An Earth-like numerical dynamo model. Nature 389, 371–374.

Laney, C.B., 1998. Computational Gasdynamics. Cambridge University Press, Cambridge.

Lee, D.C., Halliday, A.N., 1997. Core formation on Mars and differentiated asteroids. Nature 388, 854–857.

Lillis, R.J., Frey, H.V., Manga, M., 2008a. Rapid decrease in martian crustal magnetization in the Noachian era: Implications for the dynamo and climate of early Mars. Geophys. Res. Lett. 35 (L14), 203.

Lillis, R.J., Frey, H.V., Manga, M., Mitchell, D.L., Lin, R.P., Acuña, M.H., Bougher, S.W., 2008b. An improved crustal magnetic field map of Mars from electron reflectometry: Highland volcano magmatic history and the end of the martian dynamo. Icarus 194, 575–596.

Linden, P.F., 1977. The flow of a stratified fluid in a rotating annulus. J. Fluid Mech. 79, 435–447. http://dx.doi.org/10.1017/S0022112077000251.

Linden, P.F., 1979. Mixing in stratified fluids. Geophys. Astrophys. Fluid Dyn. 13, 3–23. http://dx.doi.org/10.1080/03091927908243758.

Manga, M., Stone, H.A., 1993. Buoyancy-driven interactions between two deformable viscous drops. J. Fluid Mech. 256, 647–683. http://dx.doi.org/10.1017/S0022112093002915.

Marinova, M.M., Aharonson, O., Asphaug, E., 2008. Mega-impact formation of the Mars hemispheric dichotomy. Nature 453, 1216–1219. http://dx.doi.org/10.1038/nature07070.

Marinova, M.M., Aharonson, O., Asphaug, E., 2011. Geophysical consequences of planetary-scale impacts into a Mars-like planet. Icarus 211, 960–985. http://dx.doi.org/10.1016/j.icarus.2010.10.032.

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

Monteux, J., Ricard, Y., Coltice, N., Dubuffet, F., Ulvroda, M., 2009. A model of metal-silicate separation on growing planets. Earth Planet. Sci. Lett. 287, 353–362.

Monteux, J., Jellinek, A.M., Johnson, C.L., 2011. Why might planets and moons have early dynamos? Earth Planet. Sci. Lett. 310, 349–359. http://dx.doi.org/10.1016/j.epsl.2011.08.014.

Monteux, J., Schaeffer, N., Amit, H., Cardin, P., 2012. Can a sinking metallic diapir generate a dynamo? J. Geophys. Res. (Planets) 117 (16), E10005. http://dx.doi.org/10.1029/2012/JE004075.

Morris, S., 1982. The effects of a strongly temperature-dependent viscosity on slow flow past a hot sphere. J. Fluid Mech. 124, 1–26. http://dx.doi.org/10.1017/S0022112082002389.

Nimmo, F., Kleine, T., 2007. How rapidly did Mars accrete? Uncertainties in the Hf W timing of core formation. Icarus 191, 497–504. http://dx.doi.org/10.1016/j.icarus.2007.05.002.

Nimmo, F., Stevenson, D.J., 2000. Influence of early plate tectonics on the thermal evolution and magnetic field of Mars. J. Geophys. Res. 105, 11969–11980.

Nimmo, F., Hart, S.D., Korycansky, D.G., Agnor, C.B., 2008. Implications of an impact origin for the martian hemispheric dichotomy. Nature 453, 1220–1223. http://dx.doi.org/10.1038/nature07025.

O'Keeffe, J.D., Ahrens, T.J., 1977. Impact-induced energy partitioning, melting, and vaporization on terrestrial planets. In: Merril, R.B. (Ed.), Lun. Planet. Sci. Conf., vol. 8, pp. 3357–3374.

Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Series in Computational Methods in Mechanics and Thermal Sciences. Hemisphere Publishing Corp., Washington, DC, 210p.

Peaceman, D.W., Rachford, H.H., 1955. The numerical solution of parabolic and elliptic differential equations. J. Soc. Ind. Appl. Math. 3, 28–41. http://dx.doi.org/10.1137/0103003.

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

Ratcliff, J.T., Tackley, P.J., Schubert, G., Zebib, A., 1997. Transitions in thermal convection with strongly variable viscosity. Earth Planet. Inter. 102, 201–212.

Reese, C.C., Solomatov, V.S., 2010. Early martian dynamo generation due to giant impacts. Icarus 207, 82–97. http://dx.doi.org/10.1016/j.icarus.2009.10.016.

Ricard, Y., Šrá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. http://dx.doi.org/10.1016/j.icarus.2011.11.038.

Roberts, J.H., Barnouin, O.S., 2012. The effect of the Caloris impact on the mantle dynamics and volcanism of Mercury. J. Geophys. Res. (Planets) 117, E02007. http://dx.doi.org/10.1029/2011JE003876.

Roberts, J.H., Zhong, S., 2006. Degree-1 convection in the martian mantle and the origin of the hemispheric dichotomy. J. Geophys. Res. (Planets) 111 (10), E06013. http://dx.doi.org/10.1029/2005JE002668.

Roberts, J.H., Lillis, R.J., Manga, M., 2009. Giant impacts on early Mars and the cessation of the martian dynamo. J. Geophys. Res. (Planets) 114 (13), E04009. http://dx.doi.org/10.1029/2008JE003287.

Roe, P.L., 1986. Characteristic-based schemes for the Euler equations. Ann. Rev. Fluid Mech. 18, 337–365.

Rosenblatt, P., 2011. The origin of the martian moons revisited. Astron. Astrophys. Rev. 19, 44. http://dx.doi.org/10.1007/s00159-011-0044-6.

Rubie, D.C., Melosh, H.J., Reid, J.E., Liebske, C., Righter, K., 2003. Mechanisms of metal-silicate equilibration in the terrestrial magma ocean. Earth Planet. Sci. Lett. 205, 239–255.

Rybczynski, W., 1911. über die fortschreitende bewegung einer flüssigen kugel in einen medium. Bull. Acad. Sci. Cracovie 1, 40–46.

Safronov, V.S., 1978. The heating of the Earth during its formation. Icarus 33, 3–12. http://dx.doi.org/10.1016/0019-1035(78)90019-2.

Samuel, H., Tackley, P.J., 2008. Dynamics of core formation and equilibration by negative diapirism. Geochem. Geophys. Geosyst. 9, 6011–6026.

Samuel, H., Tackley, P.J., Evonuk, M., 2010. Heat partitioning in terrestrial planets during core formation by negative diapirism. Earth Planet. Sci. Lett. 290, 13–19.

Schubert, G., Turcotte, D.L., Olson, P., 2001. Mantle Convection in the Earth and Planets. Cambridge University Press.
Senshu, H., Kuramoto, K., Matsui, T., 2002. Thermal evolution of a growing Mars. *J. Geophys. Res.* 107, 1–13.

Shannon, M.C., Agee, C.B., 1996. High pressure constraints on percolative core formation. *Geophys. Res. Lett.* 23, 2717–2720.

Sleep, N.H., 1994. Martian plate tectonics. *J. Geophys. Res.* 99, 5639. http://dx.doi.org/10.1029/94JE00216.

Smith, J.V., 1979. Mineralogy of the planets: A Voyage in space and time. *Mineral. Mag.* 43, 1–89.

Srámek, O., Ricard, Y., Dubuffet, F., 2010. A multiphase model of core formation. *Geophys. J. Int.* 181, 198–220. http://dx.doi.org/10.1111/j.1365-246X.2010.04528.x.

Sreenivasan, B., Jellinek, A., 2012. Did the Tharsis plume terminate the martian dynamo? *Earth Planet. Sci. Lett.* 349, 209–217. http://dx.doi.org/10.1016/j.epsl.2012.07.013.

Stanley, S., Elkins-Tanton, L., Zuber, M.T., Parmentier, E.M., 2008. Mars' paleomagnetic field as the result of a single-hemisphere dynamo. *Science* 321, 1822–1825. http://dx.doi.org/10.1126/science.1161119.

Stevenson, D.J., 1989. Formation and early evolution of the Earth. In: Peltier, W.R. (Ed.), *Mantle Convection and Plate Tectonics*, pp. 818–868.

Stevenson, D.J., 2001. Mars' core and magnetism. *Nature* 412, 214–219.

Stevenson, D.J., 2003. Planetary magnetic fields. *Earth Planet. Sci. Lett.* 208, 1–2.

Stevenson, D.J., Spohn, T., Schubert, G., 1983. Magnetism and thermal evolution of the terrestrial planets. *Icarus* 54, 466–489.

Thayalan, V., Jellinek, A.M., Lenardic, A., 2006. Recycling the lid: Effects of subduction and stirring on boundary layer dynamics in bottom-heated planetary mantle convection. *Geophys. Res. Lett.* 33, L20318. http://dx.doi.org/10.1029/2006GL027668.

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

Watters, W.A., Zuber, M.T., Hager, B.H., 2009. Thermal perturbations caused by large impacts and consequences for mantle convection. *J. Geophys. Res. (Planets)* 114, E02001. http://dx.doi.org/10.1029/2007JE002964.

Wilhelms, D.E., Squyres, S.W., 1984. The martian hemispheric dichotomy may be due to a giant impact. *Nature* 309, 138–140. http://dx.doi.org/10.1038/309138a0.

Williams, J., Nimmo, F., 2004. Thermal evolution of the martian core: Implications for an early dynamo. *Geology* 32, 97–100.

Yiantsios, S.G., Davis, R.H., 1990. On the buoyancy-driven motion of a drop towards a rigid surface or a deformable interface. *J. Fluid Mech.* 217, 547–573. http://dx.doi.org/10.1017/S0022112090000842.

Yoder, C.F., Konopliv, A.S., Yuan, D.N., Standish, E.M., Folkner, W.M., 2003. Fluid core size of Mars from detection of the solar tide. *Science* 300, 299–303. http://dx.doi.org/10.1126/science.1079645.

Yoshino, T., Walter, M.J., Katsura, T., 2003. Core formation in planetesimals triggered by permeable flow. *Nature* 422, 154–157.

Ziethe, R., Spohn, T., 2007. Two-dimensional stokes flow around a heated cylinder: A possible application for diapirs in the mantle. *J. Geophys. Res.* 112, 1–13.
RESEARCH ARTICLE

10.1002/2013JE004587

**Key Points:**

* Following a giant impact, the impact parameters govern the merging temperature

* For plausible viscosity contrasts the merging time occurs in less than 1 kyr

* Core merging only delays the re-initiation of the dynamo for a very short time

**Correspondence to:**

J. Monteux,
julien.monteux@univ-nantes.fr

**Citation:**

Monteux, J., and J. Arkani-Hamed (2014), Consequences of giant impacts in early Mars: Core merging and Martian dynamo evolution, *J. Geophys. Res. Planets*, 119, 480–505, doi:10.1002/2013JE004587.

Received 28 NOV 2013
Accepted 29 JAN 2014
Accepted article online 6 FEB 2014
Published online 6 MAR 2014

# Consequences of giant impacts in early Mars: Core merging and Martian dynamo evolution

Julien Monteux¹ and Jafar Arkani-Hamed²

¹Laboratoire de Planétologie et de Géodynamique, Université de Nantes, Nantes, France, ²Department of Physics, University of Toronto, Toronto, Ontario, Canada

**Abstract** A giant impact is an increasingly popular explanation for the formation of the northern lowland on Mars. It is plausible that at the impact time both Mars and the impactor were differentiated with solid silicate mantles and liquid iron cores. Such a large impact likely resulted in merging of the cores of both bodies, a process which will have implications on the thermal state of the planet. We model the evolution of the Martian mantle following a giant impact and characterize the thermochemical consequences of the sinking of an impactor's core as a single diapir. The impact heating and the viscous heating induced during the core merging may affect the early thermal state of Mars during several tens of million years. Our results show that large viscosity contrasts between the impactor's core and the surrounding mantle silicates can reduce the duration of the merging down to 1 kyr but do not modify the merging temperature. When the viscosity contrast between the diapir and the surrounding silicates is larger than a factor of 1000, the descent of the diapir can lead to some entrainment of the relatively shallow silicates to deepest regions close to the core-mantle boundary. Finally, the direct impact heating of Martian core leads to thermal stratification of the core and kills the core dynamo. It takes on the order of 150–200 Myr to reinitiate a strong dynamo anew. The merging of the impactor's core with the Martian core only delays the reinitiation of the dynamo for a very short time.

## 1. Introduction

Giant impacts were common in the later stage of accretion of terrestrial planets. The Earth is likely formed by accretion of a few dozen moons to Mars-size planetary embryos (see the review by Chambers (2004)). A Mars-size impact on Earth may have resulted in the formation of the Moon [*Hartmann and Davis, 1975*; *Cameron and Benz, 1991*; *Canup, 2004*]. An oblique collision of a large body with a mass about one sixth of the Mercury's has likely stripped away a significant part of Mercury's mantle [*Smith, 1979*; *Benz et al., 1988*]. As an alternative to endogenic formation models [e.g., *Zhong and Zuber, 2001*; *Reese and Solomatov, 2006*], an 800–1300 km size impact has been proposed to explain the formation of northern lowland of Mars [*Wilhelms and Squyres, 1984*; *Nimmo et al., 2008*; *Andrews-Hanna et al., 2008*; *Marinova et al., 2008*], and the lunar farside highlands were likely formed by accretion of a companion moon [*Jutzi and Asphaug, 2011*].

Part of an impactor may be retained by the planet, depending on the impact angle, impact velocity and the impactor to planet mass ratio. The accretion efficiency is close to 1 for impacts occurring at low velocities (i.e., the impact velocity $v_{imp} <$ twice the escape velocity of the planet $v_{esc}$) with impact angles of ~45° or less relative to the local vertical direction, whereas for impacts occurring at higher velocities (i.e., $v_{imp} > 2v_{esc}$) and at higher impact angles the efficiency is negligible [*Asphaug, 2010*]. The impact-induced shock wave propagation heats the planet and results in a large thermal anomaly buried deep in the mantle of the planet with a volume comparable to the volume of the impactor [e.g., *Tonks and Melosh, 1993*; *Pierazzo et al., 1997*; *Arkani-Hamed, 2005*; *Monteux et al., 2007*; *Arkani-Hamed and Ghods, 2011*]. This thermal energy added to the preimpact temperature of a planet may result in melting and iron-silicate segregation in the mantle.

Recent studies have explored the effects of large impacts on the Martian mantle dynamics [e.g., *Watters et al., 2009*; *Ghods and Arkani-Hamed, 2011*; *Roberts and Arkani-Hamed, 2012*] and on the cessation of dynamo action [e.g., *Roberts et al., 2009*; *Arkani-Hamed and Olson, 2010a, 2010b*; *Arkani-Hamed, 2012*]. At the time of a giant impact, both Mars and the impactor were probably differentiated [*Yoshino et al., 2003*; *Neumann et al., 2012*]. It is likely that the entire impactor's iron core merged with the preexisting core of Mars. The processes involved in the merging dynamics are strongly dependent on the rheology of the impacted mantle where the iron core of the impactor sinks and on the characteristic size of the sinking core. These processes can involve sedimentation of metallic rain at the base of a magma ocean [*Rubie et al., 2003*; *Höink et al., 2005*]
followed by a large diapir, which sinks through the underlying solid mantle [Tonks and Melosh, 1992]. The impact heating can have severe consequences on the thermal state of the planet and on the cessation of the core dynamo [e.g., Arkani-Hamed and Olson, 2010a, 2010b]. Also, the hot and molten iron brought by the merging of an impactor's core can envelop the preexistent core and generate dynamo if the layer is thicker than about 30 km [Reese and Solomatov, 2010].

Monteux et al. [2013] investigated the dynamics and the thermal effects of the core merging by monitoring the evolution of a spherical diapir sinking through a solid mantle. The authors used some idealized models of Mars and the impact heating as first-order estimates. For example (1) the iron diapir is assumed to start as a spherical body. During the impact process, the molten iron of the impactor is mixed with the partially molten parts of the impactor's and planet's silicates. However, due to its higher density compared to the silicates, the molten iron segregates and pounds at the bottom of the partially molten zone in the mantle while conforming to the shape of the zone. Hence, the initial shape of the iron diapir is a thin spherical cap with a flat top and an almost spherical bottom, certainly far from being spherical. Because the dynamics of a sinking body is strongly affected by its shape the spherical shape assumption has significant effects on the sinking time of the diapir. (2) The impact creates a spherical isobaric zone with an almost constant temperature, and temperature decays exponentially outside this zone. By choosing the decay scaling factor of ~($r_o/r$)⁴.⁴, where $r_o$ and $r$ denote the radius of the isobaric sphere and the distance from the center of the sphere, the authors severely underestimated the impact-induced temperature increase outside the isobaric core. This resulted in much colder thus higher viscosity lower mantle beneath the impact site in their temperature-dependent viscosity models, hence hampering the sinking of the iron diapir. (3) The mantle viscosity was assumed either constant (isoviscous) or temperature dependent. The mantle viscosity is actually temperature and pressure dependent, decreases with increasing temperature and increases with increasing pressure. Also the viscosity reduces substantially upon partial melting. Moreover, in reality, the rheology also depends on the strain rate and plasticity can occur [Golabek et al., 2009]. The latter is not taken into account here.

The aim of the present study is to monitor the thermal evolution of the Martian mantle following a giant impact and, in particular, the dynamics of the sinking iron diapir of the impactor. We relax the assumptions made by Monteux et al. [2013] by using more realistic models. This provides a means to assess the adverse effects of the assumptions. For example, (1) we start with the molten iron of the impactor as an almost spherical cap, numerically determined by the shape of the partially molten impact-heated region in the upper mantle of the Mars models, (2) we use a more realistic scaling law for the decay of temperature outside the isobaric zone derived on the basis of hydrocode numerical simulations of an impact [Pierazzo et al., 1997; Watters et al., 2009], and (3) we use not only constant and temperature-dependent viscosity models but also incorporate the pressure- and melt-dependent viscosity models. Furthermore, we investigate the effects of other factors such as (1) the higher viscosity contrast between the iron diapir and the surrounding mantle, (2) the higher impact velocity of 10 km/s which is within the range of impact velocities of Mars [Neukum and Wise, 1976] (on Mars, $v_{esc} = 5$ km/s), (3) different physical properties of Martian mantle which control the Equation Of State (EOS) for propagation of shock waves, hence the ability of the mantle to convert the energy from the impact-induced shock wave to internal heating, and (4) the influence of melting on the heat budget and mantle dynamics. In comparison with the models from Monteux et al. [2013], the preimpact thermal state, the postimpact thermal state, the postimpact chemical repartition of the impactor's core, and the mantle viscosity are considered more realistically. We also investigate in detail the effect of an impact on the core dynamo of Mars. This includes two aspects: (1) the direct heating of the core by the impact-induced shock wave, and (2) the impactor core merging on the Martian core. The first aspect, which is the most important factor, was not considered by Monteux et al. [2013], but the second aspect was also investigated by the authors.

The second section deals with the preimpact Martian model, while the third section presents the postimpact temperature and chemical repartitions. The core merging dynamics is described in the fourth section, and the results are presented in the fifth, sixth, and seventh sections. The effects of impacts on the core dynamo of Mars are investigated in the following section. The concluding remarks are relegated to the last section.

## 2. Preimpact Martian Model

This section describes the preimpact models of Mars considered in this study. The models are spherically symmetric with radially dependent physical parameters.
## 2.1. Pressure and Gravity Profiles

The Martian topographic dichotomy is one of the oldest features on Mars, formed sometime between 4.5 and 3.7 Ga [McGill and Dimitriou, 1990; Frey et al., 2002; Solomon et al., 2005]. Hence, at the time of the giant impact that might be responsible for the Martian dichotomy, both the impactor and Mars were probably differentiated. Hf/W chronology suggests that core formation happened during the first 10–30 Myr of Mars’ history [Lee and Halliday, 1997; Kleine et al., 2002; Nimmo and Kleine, 2007; Dauphas and Pourmand, 2011]. Such a rapid process involves extensive melting, which is potentially enhanced by radiogenic heating as a result of the decay of short-lived radionuclides ²⁶Al and ⁶⁰Fe [Yoshino et al., 2003; Neumann et al., 2012], impact heating [Tonks and Melosh, 1992; Senshu et al., 2002; Monteux et al., 2009], and gravitational energy conversion during the core formation resulting in metal/silicate separation [Stevenson, 1989; Ricard et al., 2009]. The metal/silicate separation can occur via a wide range of phenomena such as percolation [Shannon and Agee, 1996; Walte et al., 2007; Bagdassarov et al., 2009; Walte et al., 2011], the sedimentation of metallic rain through a magma ocean [Rubie et al., 2003; Höink et al., 2005] or a large diapir sinking through a solid mantle [Tonks and Melosh, 1992; Monteux et al., 2009]. Whatever the mechanism, Mars’ internal structure, characterized by a ≈ 1700 km radius liquid iron core [Yoder et al., 2003], was mostly established within ≈ 10 Myr of the planet’s formation [Kleine et al., 2002].

Considering a simple two-layered density model, $ρ_{Fe}$ for the core (with radius $R_{core}$) and $ρ_{Si}$ for the mantle, the integration of

$$dP = -ρgdr \qquad (1)$$

leads to the following gravity and pressure profiles in the core,

$$g(r) = \frac{4}{3} \pi G \rho_{Fe} r \qquad (2)$$

$$P(r) = \frac{2}{3} \pi G \rho_{Fe}^2 (R_{core}^2 - r^2) + \frac{2}{3} \pi G \rho_{Si}^2 (R^2 - R_{core}^2) + \frac{4}{3} \pi G \rho_{Si} R_{core}^3 (\rho_{Fe} - \rho_{Si}) \left(\frac{1}{R_{core}} - \frac{1}{R}\right) \qquad (3)$$

and in the mantle,

$$g(r) = \frac{4}{3} \pi G \left( r \rho_{Si} + \frac{R_{core}^3}{r^2} (\rho_{Fe} - \rho_{Si}) \right) \qquad (4)$$

$$P(r) = \frac{2}{3} \pi G \rho_{Si}^2 (R^2 - r^2) + \frac{4}{3} \pi G \rho_{Si} R_{core}^3 (\rho_{Fe} - \rho_{Si}) \left(\frac{1}{r} - \frac{1}{R}\right) \qquad (5)$$

where $ρ_{Fe}$ and $ρ_{Si}$ are assumed constant, $g$ is the gravitational acceleration, $r$ is the distance to the center of the planet, and $G$ is the gravitational constant. Figure 1 presents the density, gravity, and pressure profiles used in our models.

## 2.2. Temperature Field

The preimpact thermal state of Mars strongly depends on the short-lived radiogenic heating, the accretion processes and the dissipation of potential energy during the core formation [e.g., Senshu et al., 2002; Golabek et al., 2009; Šrámek et al., 2010]. 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.

The preimpact Martian temperature field is assumed only radially dependent. From the surface, the temperature increases rapidly through the top thermal boundary layer and then increases adiabatically and gradually down to the bottom thermal boundary layer, where it again increases rapidly toward the core-mantle boundary (CMB). We use the temperature profile from Roberts and Arkani-Hamed [2012] with the CMB temperature of $T_{CMB} = 2000$ K as the preimpact temperature profile (see Figure 1). The temperature increases adiabatically in the core as

$$\frac{dT}{dr} = - \frac{\alpha_{Fe} g(r) T}{c_{p,Fe}} \qquad (6)$$
<figure><img src="image_188.png" /><figcaption>Figure 1. The preimpact density, gravity, pressure, and temperature profiles of the Mars model.</figcaption></figure>

With the linear increase of gravity in the core (equation 2) and assuming constant thermal expansion coefficient $\alpha_{Fe}$ and heat capacity $c_{p,Fe}$, the integration of equation 6 leads to [e.g., Coradini et al., 1983]

$$T_{core}(r) = T_{CMB} \exp\left(\left(\frac{2\pi\alpha_{Fe}G\rho_{Fe}}{3c_{p,Fe}}\right)(R_{core}^2 - r^2)\right) \quad (7)$$

Included in Figure 1 is the preimpact temperature profile used in our models.

### 3. Postimpact Temperature and Chemical Repartition

After a giant impact, the shock wave propagation results in a large thermal anomaly in the upper mantle of the planet [e.g., Tonks and Melosh, 1993; Senshu et al., 2002; Monteux et al., 2007; Watters et al., 2009; Arkani-Hamed and Olson, 2010a, 2010b; Roberts and Arkani-Hamed, 2012]. Depending on the impact velocity and the efficiency of conversion from the kinetic energy to heat, the subsequent temperature increase may exceed the melting temperature of the mantle material [Pierazzo et al., 1997; Ghods and Arkani-Hamed, 2011; Monteux et al., 2011]. In the following section we describe our model of large impact heating of the Martian mantle and core.

#### 3.1. The Impact-Induced Temperature Increase

Following a giant impact, a nearly uniform shock pressure $P_s$ is generated inside the so-called isobaric sphere of radius $r_c$, [Croft, 1982; Melosh, 1989]

$$P_s = \rho_{Si}(C_0 + S u_o) u_o \quad (8)$$

where $\rho_{Si}$ and $C_0$ are the preshock density and acoustic velocity of the silicate mantle, $u_o$ is the particle velocity in the isobaric sphere, which is taken to be one half the impact velocity $v_{imp}$ assuming similar target and impacting materials, and $S$ is a constant. $D_{imp}$ is the impactor diameter and $v_{imp}$ is in km/s [Arkani-Hamed and Olson, 2010a, 2010b] of $r_c = 0.0525 D_{imp} v_{imp}$. The shock pressure decays with distance from the impact site outside the isobaric region. Several different models have been proposed for the shock pressure distribution outside the isobaric region [Ahrens and O'Keefe, 1987; Melosh, 1989; Pierazzo et al., 1997; Mitani, 2003]. In this study, we adopt the average model of Pierazzo et al. [1997]. Accordingly, the shock pressure $P_s$ at a distance $d$ from the isobaric center is calculated by

$$P_s(d) = P_s(r_c) \cdot (r_c/d)^n \quad d > r_c; n = -1.84 + 2.61 \log v_{imp} \quad (9)$$

The axisymmetric shock pressure in the Martian interior is calculated at 1 × 1 km grid intervals in a spherical coordinate system with the symmetry axis passing through the impact site and the planet's center. The
<table><caption>**Table 1.** C₀ and S Values for the Dunite</caption><thead><tr><th>C₀ (m/s)</th><th>S</th><th>Reference</th></tr></thead><tbody><tr><td>6600</td><td>0.86</td><td>McQueen et al. [1967], P &lt; 44 GPa</td></tr><tr><td>5600</td><td>1.2</td><td>McQueen et al. [1967], P &gt; 45 GPa</td></tr><tr><td>5500</td><td>1.80</td><td>Ahrens and Johnson [1995]</td></tr><tr><td>5140</td><td>1.12</td><td>Trunin [2001], 29.8 &lt; P &lt; 102.1 GPa</td></tr><tr><td>5120</td><td>1.27</td><td>Trunin [2001], 30.6 &lt; P &lt; 223.7 GPa</td></tr><tr><td>7800</td><td>0.2</td><td>Melosh [1989], 44 &lt; P &lt; 73 GPa</td></tr><tr><td>4400</td><td>1.5</td><td>Melosh [1989], P &gt; 73 GPa</td></tr></tbody></table>

pressure reduction near the surface due to interference of the direct and reflected waves is taken into consideration after modifying Melosh's [1989] rectangular coordinate algorithm to a spherical coordinate algorithm [Shahnas and Arkani-Hamed, 2007; Louzada and Stewart, 2009]. At the core-mantle boundary the incident pressure wave is partly reflected to the mantle and partly transmitted to the core [Arkani-Hamed and Ghods, 2011]. The interference of the direct

and reflected waves in the mantle near the core-mantle boundary has only minor effects on the shock pressure and is neglected.

As the impact-induced shock wave propagates in the planet it creates a high pressure behind the shock front and sets the material in motion. The pressure behind the shock front is determined using the conservation of momentum across the front,

$$P_s - P_o = \rho_o u U. \qquad (10)$$

where *u* and *U* are the particle and shock wave velocities, respectively. *P* is the total pressure (shock pressure plus the preshock lithostatic pressure) behind the shock front. *P*<sub>o</sub> and *ρ*<sub>o</sub> are the preshock lithostatic pressure and density ahead of the shock front, which are radially dependent (Figure 1). We also use an empirical, linear Equation of State (EOS) to relate the shock wave velocity to the particle velocity,

$$U = C_0 + S u. \qquad (11)$$

The acoustic velocity C₀ (= (K/ρ₀)¹/², K is the bulk modulus) and the parameter S are pressure and temperature dependent in general. We assume they are constant in our models, with different values in the core and the mantle. Table 1 is a list of measurements for dunite, which is probably a reasonable candidate for the mantle material. The S values of 1.8 by Ahrens and Johnson [1995] and 0.2 by Melosh [1989] are well off compared to the others. In our models C₀ ranges between 5140 and 6965 m/s and S between 0.8 and 1.12. The C₀ = 6965 (m/s) is the average value in the upper mantle of the Earth, from the base of the crust down to 670 km transition zone which is comparable to Martian interior as far as temperature and pressure is concerned, determined using the PREM model of the Earth [Dziewonski and Anderson, 1981].

The shock-induced temperature increase ΔT is determined following the “foundering” model of Watters et al. [2009] (see Figure 2, bottom left),

$$\Delta T = \left\{ P_{\delta}(1-1/f)/(2 \rho_o S) - (C_0/S)^2 (f - pt \ln f - 1) \right\} / C_p, \qquad (12)$$

where

$$f = -P_{\delta}/\left(\beta \left\{1 - \left[ (2 P_{\delta}/\beta) + 1 \right]^{1/2} \right\}\right), \qquad (13)$$

$$P_{\delta} = P_{s} - P_{o}, \qquad (14)$$

$$\beta = (C_0^2 \rho_o)/(2 S), \qquad (15)$$

P₀ is the preshock lithostatic pressure, ρ₀ stands for the preshock density, and Cₚ is the specific heat of the mantle.

## 3.2. Postimpact Melting and Thermal Mixing

The temperature field below the impact site and before the core merging is obtained by adding the preimpact temperature profile T₀(r) and the impact-induced temperature increase ΔT(r,θ). Considering that the preimpact mantle is solid, the postimpact temperature can exceed the solidus temperature of silicates within the Martian mantle, consume latent heat, and generate a large partially to completely molten zone. The latent heat of melting is accounted for in the heat budget. For simplicity, we consider that the solidus (Tₘₐ) and the liquidus (Tₗₐ) temperatures of silicates are constant with depth and Tₘₐ = 1800 K and Tₗₐ = 2300 K [Ghods and Arkani-Hamed, 2007]. Truly, Tₘₐ and Tₗₐ are also pressure dependent [Katz et al., 2003; Roberts and Barnouin, 2012]. The iron core of the impactor is assumed to be molten prior to the impact and no latent heat is consumed to melt it.
Pre-impact

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

For a given mantle volume, the potential temperature $T(r, \theta) = T_0(r) + \Delta T(r, \theta)$ must be modified to take into account the latent heat of melting. The fraction of melt $f_{melt}$ is calculated as follows. For a potential temperature larger than $(T_{liq} + L_{Si}/C_{pSi})$ melting is complete ($f_{melt} = 1$) and the latent heat consumption leads to a temperature reduction from $T$ to $T'$ following

$$T'(r, \theta) = T(r, \theta) - L_{Si}/C_{pSi}. \quad (16)$$

For a potential temperature lower than the $T_{sol}$, $f_{melt} = 0$ and no latent heat is consumed

$$T'(r, \theta) = T(r, \theta). \quad (17)$$

Where the potential temperature $T$ ranges between $T_{sol}$ and $T_{liq} + L_{Si}/C_{pSi}$, the latent heat consumption linearly increases within the partially molten region and

$$T'(r, \theta) = \frac{[C_p T(r, \theta)(T_{liq} - T_{sol}) + T_{sol}L_{Si}]}{[C_{pSi}(T_{liq} - T_{sol}) + L_{Si}]} \quad (18)$$

After the impact heating and melting, mechanical mixing due to segregation of iron occurs within the partially molten region leading to a homogenization of the temperature where partial melting occurs [Senshu et al., 2002; Deguen et al., 2011]. The efficiency of mixing increases with the size of the planet and

Post-impact

<figure><img src="image_190.png" /><figcaption>**Figure 2.** Preimpact temperature (in K)/compositional fields, shock-induced temperature increase (in K), and postimpact temperature (in K)/ compositional fields after a 750 km diameter impact with $C_0 = 6000$ m/s and $S = 0.86$ (Reference case). In the right bottom figure, the blue colored material represents partially molten silicates.</figcaption></figure>

gravitational energy released per impactor mass [Pierazzo et al., 1997; Deguen et al., 2011]. Hence, in the partially to completely molten volume, called hereafter the segregation volume, both silicates and iron are in thermal equilibrium. As a rough estimate we determine the temperature inside the segregation volume by volume averaging the temperatures of the partially and completely molten material.

The final postimpact temperature field strongly depends on $D_{imp}$, $C_0$ and $S$ (see Table 2). Figure 2 (second column) shows the preimpact and postimpact temperature fields for $D_{imp} = 750$ km, $C_0 = 6000$ m/s, and $S = 0.86$ model,

<table><caption>**Table 2.** Free Parameters Range Studied in Our Models</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>C<sub>0</sub></td><td>Acoustic velocity of the mantle</td><td>5140, 6000, and 6965 m/s</td></tr><tr><td></td><td>Acoustic velocity of the core</td><td>3800 m/s</td></tr><tr><td>S</td><td>EOS Constant of the mantle</td><td>0.8, 0.86, and 1.12</td></tr><tr><td></td><td>EOS Constant of the core</td><td>1.6</td></tr><tr><td>v<sub>imp</sub></td><td>Impact velocity</td><td>10 km/s</td></tr><tr><td>D<sub>imp</sub> (=2 R<sub>imp</sub>)</td><td>Impactor diameter</td><td>500, 750, and 1000 km</td></tr><tr><td>η<sub>0</sub></td><td>Reference viscosity</td><td>10<sup>20</sup>–10<sup>22</sup> Pa s</td></tr><tr><td>ρ<sub>0</sub></td><td>Mean density</td><td>4060 kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>Si</sub></td><td>Mantle density</td><td>3500 kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>Fe</sub></td><td>Core density</td><td>7500 kg m<sup>-3</sup></td></tr><tr><td>c<sub>T</sub></td><td>Temperature-induced viscosity contrast</td><td>10–1000</td></tr><tr><td>c<sub>Z</sub></td><td>Pressure-induced viscosity contrast</td><td>10</td></tr><tr><td>c<sub>1</sub></td><td>Temperature-induced viscosity contrast (if f<sub>melt</sub> = 0)</td><td>10–1000</td></tr><tr><td>c<sub>2</sub></td><td>Melt-induced viscosity contrast (if f<sub>melt</sub> > 0)</td><td>30–3000</td></tr></tbody></table>