which is designated as our reference model for the rest of the study. The impact angle should also influence the shape of the postimpact thermal anomaly and, hence, the shape of the segregation volume where melting occurs [Pierazzo and Melosh, 2000]. Indeed, increasing the impact angle from 0° (head-on impact) to larger values (oblique impacts) reduces the volume of the impact melt. Here we consider the upper limit case of a vertical impact where the efficiency of melt production is optimal.

### 3.3. Postimpact Chemical Repartition

Assuming that the bulk chemistry of the impactor is the same as Martian (12.5% volume of the metallic phase), a metallic volume with an equivalent radius $R_{Fe} = R_{imp}/2$ is retained when the planet is hit by an impactor with a radius $R_{imp}$. The retained impactor's core is denser than the silicates and rapidly sinks and fills the bottom of the segregation volume [Monteux et al., 2009; Deguen et al., 2011]. This leads to the formation of a complex shape diapir that starts to sink toward the Martian core while deforming the relatively colder and more viscous surrounding mantle (see Figure 2, third column).

## 4. Dynamic Model of Core Merging

Once buried below the surface, the dense metallic core of the impactor sinks toward the center of the target planet. Dissipation of the work done by buoyancy forces driving motion of this diapir occurs and causes heating in the mantle and/or the diapir.

### 4.1. Physical Model

To investigate the dynamics of the metallic diapir, we follow the procedure adopted by Monteux et al. [2009] and assume that the mantle deforms in a diffusion creep limit. Briefly, the conservation of energy leads to

$$ \frac{DT}{Dt} = \frac{\nabla^2 T}{Ra_\chi} + Di \left( \frac{\eta}{\eta_0} B\Omega - v_r \left( T + \frac{T_0(r)}{\Delta T_0} \right) r \right) \quad \text{with} \quad \frac{D}{Dt} = \frac{\partial}{\partial t} + (v \cdot \nabla) \quad (19) $$

$T$, $t$, $r$, and $v_r$ are dimensionless temperature, time, radius, and radial velocity, respectively. $Ra_\chi = \frac{\Delta\rho_0 g_0 R^3}{\kappa\eta_0}$ is the compositional Rayleigh number and $Di = \frac{\eta_0 g_0 R}{\rho_0 C_p}$ is the dissipation number (with $\kappa$ the heat diffusivity and $\overline{\rho_0 C_p}$ the average volume specific heat of the impacted body). $B = \frac{\Delta\rho_0}{\rho_0 \alpha \Delta T_0}$ is the chemical/thermal buoyancy ratio, where $\rho_0$ and $\alpha$ are the average density and thermal expansion of the planet. The density difference, which drives the flow, is a function of temperature and chemical composition. Because $B > 1$, the influence of the chemical difference between iron and silicates is preponderant in the density difference. The temperature-, pressure-, and melt-dependent viscosity is $\eta$, and $\eta_0$ is a constant reference viscosity. The dimensionless dissipation function $\Omega$ expresses the conversion of potential energy into heat through viscous dissipation,

$$ \Omega = 2\epsilon : \epsilon \quad (20) $$

where $\epsilon$ is the dimensionless shear strain rate tensor. The preimpact temperature $T_0(r)$ is radially dependent. Assuming an incompressible medium with an infinite Prandtl number, the conservation of mass and momentum are expressed as

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

and

$$ -\nabla P + \nabla \cdot \left( \frac{\eta}{\eta_0} [\nabla v + [\nabla v]] \right) + \left( \frac{T}{B} - f \right) re_r = 0 \quad (22) $$

where $v$ and $P$ are the nondimensional velocity and pressure and $e_r$ is the radial unit vector. The buoyancy force that drives the diapir toward the center of the planet increases with the volume fraction of metal $f$ that varies between 1 (pure metal) and 0 (pure silicates). The $f$ is simply advected by the flow

$$ \frac{Df}{Dt} = 0. \quad (23) $$

### 4.2. Viscosity Model

After ponding at the bottom of the segregation zone, the dense metallic phase sinks toward the core of the impacted planet. How this occurs is debated. Proposed mechanisms include percolation through a solid mantle [Shannon and Agee, 1996; Stevenson, 1990], settling of centimeter-sized droplets in a magma ocean
[Rubie et al., 2003] and large mantle diapirs [Ricard et al., 2009; Senshu et al., 2002; Tonks and Melosh, 1992]. In addition, it has been proposed that metallic diapirs can sink toward the centers of the planets via fracturing and diking [Tonks and Melosh, 1992; Solomatov, 2000; Stevenson, 2003] or via viscous deformation [Honda et al., 1993; Monteux et al., 2009] depending on the thermal state of the planets' mantle. In the frame of the impactor's core sinking as a single large metallic diapir, the rheology of the Martian mantle plays a key role since it governs its ability to deform during the merging. Under the conditions of a massive iron diapir sinking through solid mantle the large deviatoric stress generated by the diapir could lead to a nonlinear rheology [Samuel and Tackley, 2008; Golabek et al., 2009], elastoplastic deformations [Gerya and Yuen, 2007] or even to fracturing if they exceed the ultimate strength of rocks which is ~1–2 GPa [Davies, 1982].

However, several lines of evidence suggest the viscous deformation considered in the present study to be a more likely mechanism: (1) Detailed study of the distance from Mars to Phobos, the motion of the shadow of Phobos on the surface of Mars measured by Mars Orbiter Laser Altimeter, and the Viking I Lander observation of the motion of the shadow of Phobos passing over it led Bills et al. [2005] to suggest that the mantle of Mars at present is more deformable than that of the Earth, (2) in its early history, because of the combination of the heating from accretion and short-lived radioisotopes, and from conversion of potential energy during core formation, the Martian mantle was probably much hotter and softer than today [Senshu et al., 2002; Breuer et al., 2010], and more importantly (3) the iron diapir sinks through the already impact-heated region of the mantle below the segregation zone. Included in Figure 2 (second row, middle) is the temperature field immediately after a 750 km diameter impact (other temperature fields obtained from the different impact parameters investigated in this study will be represented later in Figures 5, 7, and 8). The temperature directly beneath the segregation zone is much closer to the solidus of the mantle, implying that the metallic core likely sinks via diffusional flow deformation (Newtonian flow) [Honda et al., 1993; Monteux et al., 2013]. Hence, excluding power law rheologies and plasticity seems to be reasonable as a first-order approximation. However, as the mantle viscosity in early Mars may be sufficiently large to support large stresses, more sophisticated models considering these two effects might be envisioned. Due to the numerical limitations of our model, this is currently beyond the scope of the present study.

To account for the viscosity contrast between the hot molten material and the cold solid material, we adopt a temperature-, melt-, and pressure-dependent viscosity model:

$$ \eta(T, P)/\eta_0(T_0, P_0) = \exp(-c_T T + c_Z P). \quad (24) $$

where $\eta_0 = \eta(T_0, P_0)$ is the viscosity of the cold and shallow undifferentiated material at the start of the experiment, $c_Z$ is a constant, and $c_T$ depends on the melt fraction of the material. For completely solid, $c_T$ equals $c_1$, and $c_2$ for completely molten regions (see Table 2). In the partially molten region $c_T$ is linearly interpolated between $c_2$ and $c_1$,

$$ c_T = (c_2 - c_1)/(T_{\text{liq}} - T_{\text{sol}})T(r, \theta) + c_1 - (c_2 - c_1)/(T_{\text{liq}} - T_{\text{sol}})T_{\text{sol}}. \quad (25) $$

This type of viscosity model 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 dependency of the viscosity for two reasons. First, the compositional contrasts are sharp and may lead to steep viscosity fronts, difficult to handle numerically. Second, among the parameters that might affect the viscosity (composition, temperature, melting degree, and pressure) experimental results suggest that the viscosity contrast is more sensitive to the melting degree and temperature. Ultimately, the viscosity contrast between hot molten iron (or fully molten silicates) and solid cold silicates can reach 20 orders of magnitude [Vocadlo et al., 2000].

To model the influence of the viscosity contrast between hot molten iron and cold solid silicates on the core merging dynamics, we use much smaller contrast. This approach is encouraged by the fact that the viscosity of the surrounding silicate material mainly governs the merging velocity [Hadamard, 1911; Batchelor, 1967]. Diapirs more than 10 times less viscous than the wall rock will behave like inviscid diapirs [Weinberg and Podladchikov, 1994]. However, shear-heating partitioning is strongly dependent on the viscosity contrast between the two media [Ke and Solomatov, 2009; Monteux et al., 2009]. In our study we consider two end-member viscosity models: models where the viscosity remains uniform and constant and models where the hot diapir is initially $10^3$ times less viscous than the coldest mantle. Further viscous heating during the core merging can eventually increase the viscosity contrasts up to $10^4$. Incorporating larger viscosity contrasts leads to numerical instabilities.
## 4.3. Numerical Model

The above dynamic equations are solved in an axisymmetric spherical geometry using 200 × 400 to 400 × 800 grid points in the *r* and *θ* directions. An implicit inversion method [Schubert et al., 2001] is adopted for the momentum equation and an Alternating Direction Implicit scheme [Peaceman and Rachford, 1955; Douglas, 1955] for the energy and mass conservation. The numerical diffusion in solving the transport equations, especially for the compositional field, is constrained using total variation diminishing Superbee scheme [Roe, 1986; Laney, 1998] in an implicit procedure [Šrámek et al., 2010], which enables a high resolution of pure advective fields. Free-slip velocity is imposed at the surface and along the axis of symmetry. The core-mantle boundary is a mobile interface that can deform in response to the merging of the impactor's core. Thermal boundary conditions are isothermal at the surface and insulating along the symmetry axis. The spherically symmetric boundary conditions are applied at the center of the planet.

## 4.4. Simplifications

A number of approximations are made in the numerical models to simplify the analysis of the core merging process:

1. We assume that the heat capacities per unit volume of iron and silicates are equal and representative of a mixture: $ρC_p = ρ_{Fe}C_{p,Fe} = ρ_{Si}C_{p,Si}$. Actually, $ρ_{Fe}C_{p,Fe} = 4500 \text{ kJ K}^{-1} \text{ m}^{-3}$ and $ρ_{Si}C_{p,Si} = 4200 \text{ kJ K}^{-1} \text{ m}^{-3}$. We use $ρC_p = 4312 \text{ kJ K}^{-1} \text{ m}^{-3}$.

2. The characteristic diffusion time ($\propto R_{Fe}^2/\kappa$) during the core merging is much larger than the characteristic sinking time ($\propto \eta_0 / \Delta\rho_0 g R_{Fe}$) even with the metallic thermal conductivity $k_{Fe}$ (= 40 W m$^{-1}$ K$^{-1}$) [Stacey and Anderson, 2001] larger than the silicate thermal conductivity $k_{Si}$ (= 3–4 W m$^{-1}$ K$^{-1}$) [Hofmeister, 1999]. Hence, the heat diffusion during the sinking is negligible and the thermal conductivities for metallic and silicate material are assumed to be equal.

3. The heat partitioning during the metal/silicate separation is strongly dependent on the mechanisms of iron segregation. If liquid iron is transported to the center of the planet through channels, the potential energy related to core formation would largely go into the core [Ke and Solomatov, 2009]. However, when the liquid metal is transported via a large diapir and depending on the viscosity contrast between the diapir and the surrounding silicates, the potential energy of the diapir that is released via viscous heating heats the surrounding silicates causing further melting [Ke and Solomatov, 2009; Monteux et al., 2009; Samuel et al., 2010]. Also, crystallization of the molten silicate occurs especially near the cold Martian surface. These two processes can influence the heat budget and affect the temperature evolution. Moreover, the viscous heating can have consequences on the rheology of the material surrounding the diapir and modify the sinking dynamics. In our dynamic model of core merging, we do not consider any consumption or release of latent heat during the sinking of the diapir. However, in the melt-dependent viscosity models, the viscosity increases abruptly by a factor of $c_2$ when temperature within a grid cell decreases below the solidus. On the opposite, when the temperature within a grid cell exceeds $T_{liq}$ its viscosity is reduced abruptly by a factor $c_2$ (see Table 2 for values).

4. During the sinking of the large volume of the impactor's iron core, the gravity field locally and temporally changes slightly. We neglect the changes of the gravity during the sinking.

5. The impactor's iron core, the silicate mantle, and the iron core of the planet are considered in the infinite Prandtl number regime. In reality, the molten metallic cores have viscosities that are about 20 orders of magnitude smaller than the mantle viscosity. Such large but realistic viscosity contrasts may have consequences on both the sinking and the merging dynamics and may affect the heat dissipation during the sinking [Ke and Solomatov, 2009; Monteux et al., 2009, 2011; Samuel et al., 2010]. However, in the case of Newtonian rheology as considered here, increasing the viscosity contrast between the hot liquid diapir and the surrounding silicates may not significantly change the sinking velocity of the diapir which is mainly controlled by the mantle viscosity [Hadamard, 1911; Rybczynski, 1911; Yiantsios and Davis, 1990]. Diapirs more than 10 times less viscous than the wall rock will behave like inviscid bubbles and further decrease in viscosity will not alter the behavior [Hadamard, 1911; Rybczynski, 1911; Weinberg and Podladchikov, 1994]. When considering a power law rheology, a large viscosity contrast between the diapir and the mantle might affect the terminal velocity of the diapir and induce differences with the theoretical velocities predicted with a Newtonian rheology [Weinberg and Podladchikov, 1994].
<figure><img src="image_191.png" /><figcaption>**Figure 3.** (left to right) Thermochemical readjustment for the model with $D_{\text{Imp}} = 750$ km, $C_0 = 6000$ m s⁻¹, $S = 0.86$, and $\eta_0 = 10^{21}$ Pa s (isoviscous) at certain times after the impact. In each panel the left part shows the temperature distribution, and the right part shows the iron fraction. The blue colored material in the right part represents partially molten silicates.</figcaption></figure>

## 5. Reference Model (Influence of the Initial Shape)

Figure 3 shows the thermochemical evolution for the reference model where $D_{\text{Imp}} = 750$ km, $C_0 = 6000$ m s⁻¹, $S = 0.86$, and $\eta = 10^{21}$ Pa s (isoviscous). The dense iron diapir rapidly sinks toward the planet's core. During the sinking, the potential energy of the diapir is converted into heat via viscous dissipation and heats both the mantle and the diapir. The diapir becomes slightly hotter than the Martian core and spreads on the core in a short time. However, because of the small density difference between the molten and solid silicates, it takes much longer time for the hot and molten silicates to ascend and spread beneath the rigid lithosphere of the planet, potentially forming a magma ocean. The spreading is followed by the initiation of convection in the Martian mantle where a cold slab sinks toward the CMB from the edge of the hot magma ocean.

During sinking, the metallic diapir deforms the surrounding mantle over a radial distance that depends on the diameter of the diapir, the effective viscosities of the mantle and the diapir, the distance between the diapir and the planet's core (i.e., the gap thickness), and the density difference across the CMB. The core of the planet has a minor effect on the dynamics of the diapir at a large gap thickness and the diapir sinks with a Stokes-like velocity [Batchelor, 1967; Monteux et al., 2013]. When the gap thickness is small, the mantle flow that enables the diapir sinking is confined between two mechanical boundary layers, one at the bottom of the diapir and the other at the top of the core. The horizontal shear stress in these boundary layers hampers the sinking of the diapir, and the sinking velocity decreases.

In Figure 4 (top, black lines), we monitor the gap thickness for our reference case. The gap thickness is calculated using the composition profile along the symmetry axis and particularly the difference in positions where the second derivative is zero. Because of some numerical diffusion of the compositional field, the final gap thickness of our models is slightly larger than the vertical spatial resolution of our models. As illustrated also in Figure 4 (top, black lines), the final gap thickness decreases for higher numerical resolutions. We consider that the time when the final gap thickness is reached corresponds to the “merging time” $t = t_m$ (i.e., the two cores are about to enter physically in contact). The characteristic coalescence time $t_c$ in case of a spherical diapir is [Monteux et al., 2013]

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

where $\lambda_0$ is the mantle viscosity/diapir viscosity ratio, $h_0$ is the initial gap thickness, and $a_1$ is a geometrical constant that is typically on the order of 0.1 for the isoviscous case and 0.2 for the T-dependent viscosity [Monteux et al., 2013]. Calculating the characteristic coalescence time from equation 26, we obtain a characteristic sinking time $t_c = 0.36$ Myr. From Figure 4 (top, black lines), we observe that our numerical results (where the merging time $\approx 0.36-0.47$ Myr) are in agreement with the theoretical coalescence time of the spherical diapir and that the initial shape of the diapir does not significantly affect the sinking time [Weinberg and Podladchikov, 1994]. An elongated diapir could increase the surface contact between the iron and the silicated material, potentially increasing the chemical degree of equilibration. However, the sinking time is too short for such a chemical equilibration.

## 6. Influence of $C_0$, $S$, and $D_{\text{Imp}}$

This section characterizes the influence of $C_0$, $S$, $D_{\text{Imp}}$ for different viscosity models. The characteristics of these models and the main results are listed in Table 3. In particular we monitor the merging temperature $T_m$ which is
<figure><img src="image_192.png" /><figcaption>**Figure 4.** Gap thickness evolution for models with different (C₀, S) sets of (top) parameters, (middle) impactor diameters *D*<sub>Imp</sub>, and (bottom) reference viscosities η₀.</figcaption></figure>

the average temperature within the sinking diapir when $t = t_m$. The initial gap thickness is a function of the postimpact segregation volume. Hence, it depends on the (C₀, S) set of parameters and on the impactor size. We monitor the thermochemical evolution after a large impact on Mars using different sets of parameters C₀, S, and D<sub>Imp</sub>.

Figure 5 (first column) illustrates the influence of C₀ and S on the postimpact segregation volume and hence on the thermochemical state before the merging process. For the three sets of (C₀, S) examined, the initial temperature within the segregation volume slightly varies and is ≈ 2100 ±10 K. When C₀ increases and S decreases a smaller segregation volume is generated after the impact. On the other hand, when C₀ decreases and S increases, the size of the segregation volume increases. Hence, the impactor's core is initially closer to the CMB (Figure 5, top). Consequence, the core merging occurs in a short time (Figure 4, top). Within the range of C₀, S studied here the sinking time (i.e., the time needed for the two cores to enter in contact) can decrease by ≈ 35% (from 0.38 Myr for (C₀ = 6965 m s⁻¹, S = 0.8) to 0.25 Myr for (C₀ = 5140 m s⁻¹, S = 1.12)) (Figure 6, right). As C₀ and S control the shape of the segregation volume and especially the initial gap thickness between the two cores, these parameters play also a significant role on the temperature of the diapir immediately before merging. Within the range of C₀, S studied here the merging temperature can increase by ≈ 10% (from 2130 K for (C₀ = 5140 m s⁻¹, S = 1.12) to 2350 K for (C₀ = 6965 m s⁻¹, S = 0.8)) which can have important consequences on the preexisting core cooling after merging (Figure 6, left).

The Martian topographic dichotomy is related to an impact ranging from 800 to 1300 km in diameter [Marinova et al., 2008; Andrews-Hanna et al., 2008; Nimmo et al., 2008]. In our models, we consider that the impactor's core is collected at the bottom of the segregation volume. When the impactor is too large (*D*<sub>Imp</sub> > 1000 km), the segregation volume is in contact with the Martian core, and the merging dynamics is governed by the viscosities of the molten silicates and the liquid iron. Hence, the merging dynamics involves turbulent flows and our models are no more valid. Therefore, we limit our study to impactor diameters ranging from 500 km to 1000 km. Moreover, as detailed above, decreasing C₀ and increasing S would increase the size of the segregation volume, and for *D*<sub>Imp</sub> = 1000 km, C₀ < 6965 m s⁻¹ and S > 0.8, the
<table><caption>**Table 3.** Model Specifics</caption><thead><tr><th>Model</th><th>Impactor Diameter (Shell Thickness) (km)</th><th>C<sub>0</sub> (m/s)</th><th>S</th><th>Reference Viscosity (Pa s)</th><th>Viscosity Contrast</th><th>Merging Time (Myr)</th><th>Merging Temperature (K)</th></tr></thead><tbody><tr><td>M0 (ref.)</td><td>750 (0.8)</td><td>6000</td><td>0.86</td><td>10<sup>21</sup></td><td>1</td><td>0.36</td><td>2286.7</td></tr><tr><td>M1</td><td>750 (0.8)</td><td>6965</td><td>0.8</td><td>10<sup>21</sup></td><td>1</td><td>0.37</td><td>2352.5</td></tr><tr><td>M2</td><td>750 (0.8)</td><td>5140</td><td>1.12</td><td>10<sup>21</sup></td><td>1</td><td>0.24</td><td>2137.6</td></tr><tr><td>M3</td><td>1000 (1.8)</td><td>6965</td><td>0.8</td><td>10<sup>21</sup></td><td>1</td><td>8.9 10<sup>-2</sup></td><td>2127.7</td></tr><tr><td>M4</td><td>500 (0.2)</td><td>6965</td><td>0.8</td><td>10<sup>21</sup></td><td>1</td><td>1.3</td><td>2392.9</td></tr><tr><td>M5</td><td>750 (0.8)</td><td>6000</td><td>0.86</td><td>10<sup>20</sup></td><td>1</td><td>3.6 10<sup>-2</sup></td><td>2287.3</td></tr><tr><td>M6</td><td>750 (0.8)</td><td>6000</td><td>0.86</td><td>10<sup>22</sup></td><td>1</td><td>3.6</td><td>2280.5</td></tr><tr><td>M7</td><td>750 (0.8)</td><td>6000</td><td>0.86</td><td>10<sup>21</sup></td><td>10 (T-P)</td><td>0.1</td><td>2281.6</td></tr><tr><td>M8</td><td>750 (0.8)</td><td>6000</td><td>0.86</td><td>10<sup>21</sup></td><td>10 (T-P melt)</td><td>6.3 10<sup>-2</sup></td><td>2290.1</td></tr><tr><td>M9</td><td>750 (0.8)</td><td>6000</td><td>0.86</td><td>10<sup>21</sup></td><td>100 (T-P)</td><td>1.3 10<sup>-2</sup></td><td>2277.4</td></tr><tr><td>M10</td><td>750 (0.8)</td><td>6000</td><td>0.86</td><td>10<sup>21</sup></td><td>100 (T-P melt)</td><td>7.5 10<sup>-3</sup></td><td>2323.2</td></tr><tr><td>M11</td><td>750 (0.8)</td><td>6000</td><td>0.86</td><td>10<sup>21</sup></td><td>1000 (T-P)</td><td>1.5 10<sup>-3</sup></td><td>2292.5</td></tr><tr><td>M12</td><td>750 (0.8)</td><td>6000</td><td>0.86</td><td>10<sup>21</sup></td><td>1000 (T-P melt)</td><td>1.0 10<sup>-3</sup></td><td>2301.3</td></tr></tbody></table>

segregation volume also overlaps the CMB. Hence, to measure the influence of the impactor size, we only consider $C_0 = 6965 \text{ m s}^{-1}$ and $S = 0.8$. Figure 7 shows the merging dynamics for three impactor diameters. When $D_{\text{imp}}$ increases, the segregation volume and the metallic volume of the impactor increase. However, the mean temperature of the segregation volume is not significantly affected (differences of ≈ 10 K between $D_{\text{imp}} = 500 \text{ km}$ and $D_{\text{imp}} = 1000 \text{ km}$) (Figure 6, left). The major difference resides in the volume of the postimpact segregation zone of the planet which is much more significant for $D_{\text{imp}} = 1000 \text{ km}$ (see Figure 7, third row). When the impactor diameter increases, the volume of the segregation zone and the size of the merging diapir increase, resulting in the reduction of $h_0$ and in the increase of $R_{\text{Fe}}$, which reduce the sinking time from ≈ 1 Myr for $D_{\text{imp}} = 500 \text{ km}$ to ≈ 0.1 Myr for $D_{\text{imp}} = 1000 \text{ km}$ (Figure 6, right). As $D_{\text{imp}}$ controls the initial distance between the two cores, it also influences the temperature of the diapir before merging [Monteux et al., 2009]. Hence, when $D_{\text{imp}}$ increases from 500 to 1000 km, the impactor's core becomes initially closer to the CMB and $h_0$ decreases from 935 km to 120 km (Figure 4, middle). Hence, its available potential energy decreases, and the merging temperature decreases by ≈ 12% (from 2400 to 2130 K) (see Figure 6, left).

<figure><img src="image_193.png" /><figcaption>**Figure 5.** Thermochemical readjustment as a function of time for models with three ($C_0$, $S$) sets of parameters. We consider here that the viscosity $\eta_0$ is uniform and constant.</figcaption></figure>
<figure><img src="image_194.png" /></figure>

**Figure 6.** (left) Merging temperatures and (right) merging times as a function of the viscosity contrast. Green filled solid triangles represent results for uniform viscosities, with $C_0 = 6965$ m/s and $S = 0.8$ and for different impactor diameters (indicated close to the corresponding symbol). The open symbols represent the values obtained with $C_0 = 6000$ m/s and $S = 0.86$ for T-P dependent viscosities (black squares) and for T-P melt-dependent viscosities (red circles).

Finally, by controlling the segregation volume $D_{\text{imp}}$, $C_0$ and $S$ indirectly control the final radial extension of the magma ocean spreading, and the position of the cold slab that sinks from the surface toward the CMB (see Figures 5 and 7).

# 7. Influence of the Viscosity Model

## 7.1. Influence of the Average Viscosity

The mantle viscosity is a key parameter in the core merging process that governs the ability of the mantle to deform and accommodate the metallic diapir sinking, while deforming at farther distances to compensate for the accommodation. The viscosity also governs the spreading of the hot silicate material in the upper mantle. When consider dynamical models with uniform viscosities (i.e., $\eta_0$ is constant and $c_T = c_Z = 0$), modifying the mean viscosity does change the heat repartitioning at the end of the merging and affects the characteristic

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

**Figure 7.** Thermochemical readjustment as a function of time for models with three impactor sizes. We consider here that the viscosity $\eta_0$ is uniform and constant. $C_0 = 6965 \text{ m s}^{-1}$ and $S = 0.8$.
merging and spreading times. In the laminar flow regime, the terminal sinking velocity U of a metallic diapir sinking in an infinite media obeys the Hadamard velocity law [Hadamard, 1911; Batchelor, 1967; Monteux et al., 2009]. As a consequence

$$U = \frac{2}{9} \left( \frac{\eta_d + \eta_s}{\eta_d + \frac{2}{3} \eta_s} \right) \frac{\Delta \rho_0 g R_{\text{Fe}}^2}{\eta_s} \quad (27)$$

where $\eta_d$ is the viscosity inside the diapir and $\eta_s$ is the viscosity of the surrounding material. In the isoviscous case (Figure 4, bottom) when the average viscosity of the mantle decreases the diapir sinks faster (see equation 27), and the sinking time is significantly reduced. For example, in an isoviscous model, $\eta_d = \eta_s$, increasing the viscosity from $10^{20}$ Pa s (Figure 4, bottom, red line) to $10^{22}$ Pa s (Figure 4, bottom, green line) leads to an increase of the sinking times from ≈ 0.1 Myr to ≈ 10 Myr (Table 3). Similarly, an increase of the average viscosity by 2 orders of magnitude increases the spreading time of the hot mantle thermal anomaly from ≈ 1 Myr to ≈ 100 Myr, even if the thermal anomaly cools by diffusion. However, the changes in the average uniform viscosity have less effects on the merging temperatures that range between 2280 and 2290 K at $t = t_m$ (Table 3). For large viscosities, as the characteristic merging and spreading times are longer, heat diffusion plays a significant role and smoothes the temperature contrasts between the hot merging channel and the surrounding silicates, at the CMB, and between the sinking slab and the colder surrounding mantle.

## 7.2. Influence of Varying Viscosity

The viscosity contrast between the sinking diapir and its surroundings is a key parameter in the core merging dynamics and the repartitioning of the viscous heating [Ke and Solomatov, 2009; Samuel et al., 2010; Monteux et al., 2011, 2013]. Qualitatively, as the metallic diapir sinks, shear heating occurs at the interface of the diapir and the surroundings [Samuel et al., 2010; Monteux et al., 2011], leading to the temperature increase of both the diapir and the surrounding mantle. Depending on the effective viscosity contrast between the metallic diapir and the mantle, the mean temperature of the diapir sinking in a Mars-size planet can increase by a few hundred Kelvins [Monteux et al., 2009]. This temperature increase is maximal when the diapir and the surrounding mantle are isoviscous [Monteux et al., 2009]. 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 is similar to viscosity variations of order $10^2$ [Jellinek et al., 2003; Thayalan et al., 2006], the effects of which we can explore parametrically.

To characterize the influence of a variable viscosity, we monitor the thermochemical readjustment following a 750 km diameter impact and the subsequent core merging for three viscosity models: isoviscous, T-P dependent and T-P melt dependent. During the impactor's core sinking the viscous heating enhances the temperature and decreases the viscosity of the material immediately surrounding the diapir. The viscosity decrease is mainly localized in the sinking channel where the highest temperatures are reached. In this region, the temperature of the already impact-heated mantle beneath the segregation zone can exceed the solidus temperature and lead to partial melting and a significant viscosity decrease. For the T-P melt-dependent viscosity model (Figure 8, third row) and using a reference viscosity $\eta_0 = 10^{21}$ Pa s, the core merging is much faster than for the isoviscous model and occurs in less than 1 kyr when the viscosity contrast is 1000 between the hottest and the coldest material (see Table 3 for values). The spreading of the hot silicate mantle beneath the lithosphere is also facilitated by the increasing viscosity contrast. However, it is achieved within a much longer timescale that is governed by the surrounding mantle viscosity.

The radial deformation of the mantle in the isoviscous case scales with the size of the diapir [e.g., Morris, 1982; Jellinek et al., 2003; Thayalan et al., 2006], whereas in the temperature-dependent case the deformation is confined to the hottest, lowest viscosity material surrounding the diapir with a length scale typically ≈ 0.1 R$_{Fe}$. Moreover, the horizontal shear stress in the boundary layer between the diapir and the CMB that retards the sinking of the diapir is very small for the T-P and T-P melt-dependent viscosity. Both these effects increase the sinking velocity when increasing the viscosity contrast between the diapir and the surrounding mantle. Figure 9 shows the evolution of the gap thickness for three viscosity models, isoviscous (with $\eta_0 = 10^{21}$ Pa s), T-P dependent (with $c_Z = 10$ and $c_T = 10–1000$), and T-P melt dependent (with $c_Z = 10$, $c_1 = 10–1000$, and $c_2 = 30–3000$). The highest viscosity contrasts between the hot diapir and the surrounding mantle are obtained in the T-P melt-dependent viscosity model (Figure 8); hence, the coalescence time is shortest (see Figure 6, right).
<figure><img src="image_196.png" /><figcaption>**Figure 8.** Thermochemical readjustment as a function of time for an isoviscous, a *T-P* dependent (with a viscosity contrast of 100) and a *T-P* melt-dependent (with a viscosity contrast of 100 to 300) viscosity models. (rows 2 and 3) The right part of each panel shows the evolution of the viscosity during the core merging process.</figcaption></figure>

Figure 6 summarizes the influence of the viscosity contrast on the merging times (right) and the temperature of the merging diapir (left) (see values in Table 3). This figure emphasizes the weak influence of the viscosity contrast on the temperature of the impactor's core before merging with the Martian core. The temperature ranges between 2280 and 2330 K while the viscosity contrast increases by 3 orders of magnitude. As detailed above and illustrated in Figure 6 (left) the merging temperature is much more sensitive to the impactor size, while the viscosity contrast is a key parameter when considering the merging time. Figure 6 (right) shows that an increase of the viscosity contrast by 3 orders of magnitude decreases the merging time from 0.36 Myr to ≈ 1 kyr. Here the influence of the impactor size is weaker according to the results of our models.

<figure><img src="image_197.png" /><figcaption>**Figure 9.** Gap thickness evolution for three different viscosities: uniform (black solid line), *T-P* dependent (colored solid lines), and *T-P* melt dependent (colored dashed lines). For *T-P* and *T-P* melt-dependent viscosities, three viscosity contrasts are represented: 10 (red), 100 (green), and 1000 (blue).</figcaption></figure>

**8. Thermochemical Consequences of the Core Merging**

This section addresses the consequences of the core merging in terms of thermochemical equilibration.

**8.1. The Core-Mantle Equilibrium**

After the impact and the formation of the segregation volume, mechanical mixing due to segregation of iron occurs within the partially molten region and can enhance thermochemical equilibrium between the iron from the impactor and molten silicates. The degree of
<figure><img src="image_198.png" /><figcaption>**Figure 10.** (left to right) Close-up view of the time evolution of the compositional field within the core merging area using the T-P melt viscosity model with a viscosity contrast of (top) 10, (middle) 100, and (bottom) 1000. Solid lines represent the streamlines in the reference frame of Mars.</figcaption></figure>

equilibration decreases with increasing characteristic size of the metallic diapirs [Ulrová et al., 2011; Deguen et al., 2011; Samuel, 2012]. For centimeter-sized droplets sinking through a turbulent magma ocean, the thermochemical equilibration is efficient and rapidly achieved [Samuel, 2012]. However, during the core merging, the metallic phase may sink as blobs and lead to metal/silicate disequilibrium [Dahl and Stevenson, 2010; Kleine and Rudge, 2011]. In our models, the impactor's metallic phase is sinking as a large single diapir and chemical equilibrium should depend on the complexity of the flow structure within and outside the diapir [Ulrová et al., 2011]. However, from the sinking times obtained in our models (< 1 kyr for the most realistic viscosity contrast) (see Table 3), it seems difficult to envision any chemical reequilibration during the diapir sinking. As detailed above, increasing the viscosity contrast between the hot sinking diapir and the surrounding material drastically reduce the merging time but also leads to an elongation of merging core's tail because of the large-scale mantle flow and of the increase in the ability of the diapir to deform (Figure 10). If any chemical equilibration might occur, this thin tail is a good candidate.

The flow induced in the mantle during the sinking can also lead to some entrainment of the silicates from the upper to the lower mantle. This entrainment increases with the diapir size (i.e., with the impactor size) and with the viscosity contrast between the metallic diapir and the mantle (Figure 10). Once the two cores have fully merged, the flow induced by the thermal readjustment within the mantle tends to bring the deep silicate material upward. The accuracy of this phenomenon increases as the viscosity contrast between the diapir and the silicates increases (Figure 10). However, according to the streamlines shown in Figure 10 (right column) the molten silicate material entrained downward during the diapir sinking does not seem to be
<figure><img src="image_199.png" /><figcaption>**Figure 11.** Surface heat flux as a function of time and latitude for different impact properties and average viscosities.</figcaption></figure>

retrained upward after the merging. Hence, core merging subsequent to a giant impact results in the mixing of silicate material in the mantle.

## 8.2. Consequences for the Martian Cooling

### 8.2.1. Mantle Cooling

The impact heating and the viscous dissipation associated with the impactor's core sinking modify the thermal state of the Martian mantle. During and after the merging, a thermal readjustment occurs in the mantle modifying the surface heat flux. The readjustment of the postimpact thermal anomaly within the mantle and its spreading beneath the Martian lithosphere leads to an increase of the surface heat flux. During the early advective stage, corresponding to the flattening of the hot mantle anomaly, the radial extent of the high heat flux patch increases with time. This early advective stage is followed by a later stage of motionless diffusive cooling where the surface heat flux progressively decreases and the mantle thermal anomaly diminishes [Monteux et al., 2007]. The radial extent of the thermal anomaly, its thickness and the duration of the advective stage are strongly dependent on the size of the thermal anomaly and on the physical properties of the planet. Figure 11 shows the time evolution of the surface heat flux ($f_s$), as a function of the colatitude for different $C_0$, $S$, $D_{\text{imp}}$, and $\eta_0$ models. Low viscosities and large thermal anomalies favor the spreading of the anomaly over large distances. Hence, the largest sizes for these patches ($\theta \approx 60^{\circ}$) are reached for small
viscosities ($\eta_0$ = $10^{20}$ Pa s) and large impactor sizes ($D_{imp}$ = 1000 km). The time of the transition between the advective and the diffusive stages strongly depends on the mantle viscosity. For $\eta_0$ = $10^{21}$ Pa s this time is typically on the order of 10 Myr which is comparable with the models from Monteux et al. [2007].

T-P and T-P melt-dependent viscosities also influence the evolution of the surface heat flux. Large viscosity contrasts between the hot anomaly and the relatively colder mantle enable easier spreading of the anomaly below the surface [Koch and Koch, 1995; Monteux et al., 2007]. As the resistance to internal shearing decreases with increasing viscosity contrast, the horizontal velocity becomes more significant for low viscosities. As a result, for the T-P and T-P melt-dependent viscosity models, the radial extent increases by ≈ 10%, while the advection time decreases by a factor of ≈ 2 compared to the isoviscous models [Monteux et al., 2007]. Long-term monitoring of the surface heat flux can be very time consuming especially when considering large viscosity contrasts. Hence, we did not monitor the long-term surface heat flux in the nonuniform viscosity cases. However, the temperature dependence of the viscosity should affect our results by a minor amount because the readjustment is mostly controlled by the viscosity far from the segregation volume.

### 8.2.2. Core Cooling

During the sinking, viscous heating leads to a temperature increase in both the diapir and the surrounding mantle. For a uniform viscosity and for diapir larger than ≈ 70 km in diameter, viscous coupling between the metallic diapir and the surrounding material occurs and a fraction of the gravitational energy of the diapir is converted to heat up the diapir [Monteux et al., 2009]. However, when the diapir is less viscous than the mantle, viscous heating is restricted to the surrounding mantle and is concentrated at the diapir's poles where strain rates are greatest [Samuel et al., 2010; Monteux et al., 2011]. In this case, the temperature of the impactor's core does not increase significantly during the sinking and remains close to its initial presinking temperature at the base of the segregation zone [Monteux et al., 2009].

The temperature of the diapir when merging with the Martian core is also dependent on the distance traveled during the sinking and increases with the initial distance between the diapir (i.e., the initial depth of the segregation zone) and the CMB [Monteux et al., 2009]. Hence, the merging temperature of the diapir may be larger or smaller than the core temperature, depending on the size of the diapir, the distance traveled by the diapir, and the viscosity contrast with the surrounding mantle. For instance, when the segregation volume is large (Figure 5, first row and Figure 7, third row), the initial distance between the diapir and the CMB is not large enough for the diapir to heat up significantly, and the merging diapir is colder than the Martian core. For other sets of parameters the merging core is hotter than the Martian core (Table 3).

Following the merging, the impactor's core is trapped between the hot impacted core material and the hot mantle material from the channel where the diapir has sunk. The positive buoyancy of these two materials between which the impactor's core is trapped tends to spread the impactor's core beneath the CMB (Figure 12). However, the mixing dynamics that occurs within the Martian core after the merging has to be considered here with caution. The injection of a large volume of molten iron of the impactor in a rotating core of the planet is a process that involves Coriolis and inertial forces that are neglected here, because of our infinite Prandtl number approximation. We address this issue in section 10.

## 9. Effects on the Core Dynamo

A giant impact not only excavates the near surface and heats the upper mantle of a planet as discussed above but also heats the core of the planet in two distinct stages. In the first stage, the shock wave created by the impact propagates in the core, heating it differentially within about 1 h. In the second stage, the merging of the impactor iron diapir modifies the thermal state of the planet's core. Here we investigate the effects of these two stages on the core dynamo of Mars.

The first stage has been studied by many investigators [e.g., Arkani-Hamed and Olson, 2010a, 2010b; Ghods and Arkani-Hamed, 2011]. Here we follow the procedure adopted by Arkani-Hamed and Olson [2010b]. Briefly, the shock wave leads to a temperature increase within the core of the planet, much stronger in the region directly beneath the impact site (Figure 2, bottom left). The low-viscosity rotating liquid core cannot sustain lateral variations of temperature and the core overturns, resulting in a stably stratified temperature which increases with radius. The thermal stratification diminishes the possible preexisting core convection, hence the core dynamo, within a few kyr [Arkani-Hamed and Olson, 2010a]. Here we assume it occurs immediately after the impact. Shortly after, the juxtaposition of the superheated stratified core to the relatively colder
<figure><img src="image_200.png" /><figcaption>**Figure 12.** Close-up view of the temperature field within the core merging area for two impactor sizes, two sets of $C_0$ and $S$, and viscosity of $10^{21}$ Pa s.</figcaption></figure>

mantle results in convection in the outer part of the core which generates a dynamo as it thickens over time. Here we first present the results for our reference model ($D_{\text{imp}} = 750$ km, $C_0 = 6000$ m s⁻¹, $S = 0.86$ and $\eta = 10^{21}$ Pa s) in some detail and then provide the major results for the other models.

Included in Figure 2 (bottom left) is the impact-induced temperature increase inside our reference model, showing the differentially heated mantle and core. We note that the temperature in the uppermost part of the core directly beneath the impact site is higher than that in the adjacent mantle. This is partly due to shock pressure jump as the shock wave inters the core and partly because of low-specific heat of the iron core, 600 J/kg/K, compared to that of the silicate mantle, 1200 J/kg/K. Also, the impact heating of the core is appreciable only within ~40° colatitude relative to the symmetric axis passing the impact site and the center of the planet, the area which covers only 12% of the core's surface. About 88% of the core is juxtaposed to the base of the mantle with preimpact temperature. Once the core stratifies and the hottest part of the core is placed directly beneath the core-mantle boundary, the core essentially cools almost globally. The thermal evolution of the core is calculated by numerically solving the 1-D enthalpy equation and following the procedure adopted by Arkani-Hamed and Olson [2010b]. The radial grid interval is taken to be 100 m to allow accurate determination of the thickness of the thermal boundary layer in the convecting outer core. The domain of calculation includes the entire core plus a 100 km thick thermal boundary layer at the base of the mantle. The preimpact temperature is assumed adiabatic inside the core with 2000 K at the core-mantle boundary. It linearly decreases in the overlying mantle layer to 1500 K at the top of the layer. The spherically symmetric temperature distribution inside the core immediately after the impact is obtained by first adding the preimpact temperature to the impact-induced temperature increase and then allowing the thermal stratification. The liquid core has a kinetic viscosity of $10^4$ m²/s and it rotates with a 24 h period (see Table 2 for physical properties of the core). The core viscosity must be much smaller. However, reducing the viscosity by 8 orders of magnitude has only minor effects on the core cooling [Arkani-Hamed and Olson, 2010b], because the core cooling is mainly controlled by the overlying solid mantle.

To illustrate the effects of the impactor core merging on the Martian core dynamo two scenarios are calculated for each models listed in Table 3: one without considering the merging and the other with merging. Figure 13 shows the thermal evolution of the core and the overlying mantle for 300 Myr after the impact where no merging is considered. The positive temperature gradient deep in the core retains stable condition, prohibiting convection, while the upper most part of the core convects and maintains adiabatic temperature as it thickens and penetrates to deeper parts of the core. Although the very high temperature at the top of the stratified core diminishes rapidly, the impact heating remains appreciable for much longer time. Figure 14 shows the time variations of the thickness of the convecting outer core, the heat flux at the core-mantle boundary, the magnetic
<figure><img src="image_201.png" /></figure>

**Figure 13.** Thermal evolution of Martian core after an impact, where no core merging is considered. The temperature is shown at 30 Myr intervals. The numbers on the curves denote time after the impact in Myr. The curve "Initial Adiabat" is the preimpact temperature, and "0+ Stratified" denotes the temperature after the first stage of stratification. $D_{\text{imp}}$ is the impactor diameter in kilometer, $C_0$ is the acoustic velocity of the mantle in m/s, $S$ is the constant in the EOS, and $\eta_0$ is the mantle viscosity assumed isoviscous mantle.

Reynolds number and the mean magnetic intensity inside the convecting outer core for 300 Myr after the impact. Similar to the results obtained by other investigator [e.g., Arkani-Hamed and Olson, 2010b; Arkani-Hamed and Ghods, 2011], the impact certainly cripples the possible preimpact core dynamo and it takes around 150–200 Myr for the convecting outer core to generate a strong core dynamo, taking the threshold value of 20 for the magnetic Reynolds number. We note that the outer core is not thick enough to generate an appreciable dynamo within the first ~20 Myr after the impact, despite very high heat flux at the core-mantle boundary.

Due to the positive buoyancy relative to the core of Mars, the impactor diapir spreads on the core (Figure 12) and creates a super-heated iron layer of ~0.8 km

thickness. Figure 15 compares the thermal evolution of the core for the first 10 Myr with and without the diapir merging. At the time of merging the iron diapir is about 100 K hotter than the uppermost part of Mars' core. However, the iron layer cools very rapidly such that the difference between the two models almost disappears within less than 2 Myr. This is better illustrated in Figure 16 where the thickness of the convecting outer core, the heat flux at the core-mantle boundary, and the magnetic Reynolds number and the mean magnetic field intensity of the convecting outer core are displayed within the first 1 Myr after the impact for both models. The merging iron layer reinitiates convection at about 0.37 Myr, and its effects diminish rapidly within the first 1 Myr.

The effects of physical parameters of the mantle and the size of the impactor are also estimated for both scenarios, with and without merging. Figures 17 and 18 show the effects of the diapir merging on the core dynamo for different viscosity models of the mantle, different $C_0$ and $S$ values as well as different size of the impactor. The

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

**Figure 14.** The evolution of the thickness of the convecting outer core, the heat flux at the core-mantle boundary, the magnetic Reynolds number, and the mean magnetic intensity in the convecting zone corresponding to the thermal evolution of the core shown in Figure 13.
<figure><img src="image_203.png" /><figcaption>**Figure 15.** Thermal evolution of Martian core within 10 Myr after an impact. (top) With no core merging and (bottom) with core merging that occurs at 0.36 Myr after the impact. The blue and magenta curves in Figure 15 (bottom) show the temperature distribution immediately before and immediately after the core merging. See Figure 13 for details.</figcaption></figure>

isoviscous mantle with a viscosity of 10²² Pa s, and to a lesser extent 10²¹ Pa s, delays the merging time and has appreciable effects on the dynamo. It takes longer time for the iron diapir of the smaller impactor ($D_{\text{imp}} = 500$ km), which initiates at the base of a smaller segregation zone, to sink through the mantle before merging to the core, hence has appreciable effect on the core dynamo because the Martian core cools prior to the merging. The $C_0 = 6965$ (m/s) and $S = 0.8$ model results in an appreciable merging effects. However, the effects diminish within a few Myr in all models. It is plausible to conclude that the main effect of an impact on the core dynamo is due to direct heating of the core (see Figure 14), which was not considered by Monteux *et al.* [2013], and diapir merging has almost negligible effects. This is largely because the iron layers produced on the Martian core by the diapirs are very thin, and no dynamo

can be generated inside the layers. It requires a thickness of about 30 km for an iron layer to generate dynamo [e.g., Reese and Solomatov, 2010]. The iron layers can only delay the initiation of the dynamo for a very short time.

## 10. Limitations of Our Models

The stratification of the differentially heated Martian core occurs while the impactor's core diapir is still descending in the mantle. The core merging incorporates some new material within the Martian core with a

<figure><img src="image_204.png" /><figcaption>**Figure 16.** The evolution of the thickness of the convecting outer core, the heat flux at the core-mantle boundary, the magnetic Reynolds number, and the mean magnetic intensity in the convecting zone within the first 1 Myr. The red curves are for the no merging scenario, and the black ones are for the merging scenario.</figcaption></figure>
<figure><img src="image_205.png" /></figure>

**Figure 17.** Effects of the mantle viscosity on the Reynolds number and the mean magnetic intensity in the convecting outer core for a model with 750 km impactor diameter, $C_0 = 6000$ m/s, and $S = 0.86$. Eta at the top denotes the mantle viscosity: either isoviscous with $10^{20}$, $10^{21}$, or $10^{22}$ Pa s, or temperature- and pressure-dependent ($f(P,T)$), or temperature, pressure, and melt dependent ($f(P,T,M)$) with a viscosity contrast of 1000. The red curves are for the no merging scenario, and the black ones are for the merging scenario.

temperature that depends on the size of the impactor, the segregation volume and the viscosity contrast between the diapir and the mantle. In any case, the core merging leads to lateral temperature variations and may initiate restratification of the core.

Our present numerical model of spreading the iron diapir on the Martian core does not include the rotation effect (i.e., no Coriolis forces) and assumes an infinite Prandtl number (i.e., no inertial forces) in the momentum equation. Hence, we do not obtain the stratification that likely occurs immediately after the impact as discussed in the previous section, and the second restratification that would occur after the complete core merging. Consequently, materials from our lowest viscosity and largest diapir spread beneath the CMB as a thin, axisymmetric gravity current (see Figure 12). The spreading of the hot core anomaly stops when diffusive cooling overcomes advective transport leading to a partial stratification of the core [Monteux et al., 2013]. This thermal stratification persists until the thermal anomaly introduced by the gravity current diffuses into the overlying mantle and the underlying part of the core. The dynamics of thermal mixing under realistic core viscosities will be the subject of a separated study.

## 11. Conclusions

Giant impacts may have significantly influenced the mantle dynamics of Mars. We investigated the postimpact thermochemical readjustment after one single giant impact occurring during the early stages of

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

**Figure 18.** Effects of the acoustic velocity, $C_0$, and EOS constant, $S$, of the mantle on the magnetic Reynolds number and mean magnetic intensity of the convecting outer core. The mantle has a constant viscosity of $10^{21}$ Pa s. $D$ denotes the impactor diameter in kilometer. The red curves are for the no merging scenario, and the black ones are for the merging scenario.
the Martian history. In particular, we studied the thermal consequences of the sinking of an impactor's core in the Martian mantle as a function of the impactor size and for different mantle viscosity models. Our results show that the core merging is a fast process that occurs in a time mainly governed by the mantle viscosity. Although this study is the continuation of the study by Monteux et al. [2013], several improvements have been achieved in the characterization of the impact heating, the initial postimpact geometry of the impactor's core and the viscosity models. Indeed, from our results we can now separate the influences of the impactor size ($D_{\text{imp}}$), the ability of the impacted mantle to be shock heated ($C_0$ and S) and the mantle viscosity (type and viscosity contrast) on the time required for sinking and on the temperature of the impactor's core when merging with the Martian one. We can now conclude that the viscosity does not play a key role in the thermal evolution during the merging (only modifying the merging temperature by ≈ 1 % for the range of viscosities adopted here). However, it plays an important role on the sinking time which is very sensitive to the viscosity contrast between the sinking diapir and the surrounding material that can decrease from 0.3 Myr to 1 kyr. From the range of parameters considered here (with viscosity contrast up to 3 orders of magnitude), the merging time is probably too short to envision any chemical equilibration during the sinking. Concerning the thermal evolution, the impactor diameter $D_{\text{imp}}$ and the ($C_0$, S) values seem to be the key parameters that govern both the merging time (by controlling the diapir size and initial depth) and the merging temperature that can vary by 10% in the range of parameters studied here (from 2130 to 2350 K). In our models, we assume that the impact occurs with a moderate impact velocity (about twice the Martian escape velocity); larger impact velocities can substantially increase the postimpact temperature in the merging core prior to its sinking which strengthen the importance of the impact parameters on the merging temperature.

Our study also highlights that, even if the core merging processes could reduce the efficiency of chemical equilibration during the core formation in terrestrial planets [Kleine and Rudge, 2011], the flow dynamics generated within the mantle during the merging process for large viscosity contrasts could enhance some mixing between upper and deep silicate material even in the absence of mantle convection. Because the characteristic timescales of thermal readjustment within the mantle are long, the hot postimpact mantle anomaly creates a surface heat flux anomaly as it spreads below the Martian surface. The spreading stage lasts ≈ 10 Myr and is followed by a diffusive cooling that occurs within a longer timescale. More importantly, we can now make a better estimate of the influence of the core merging process on the Martian dynamo. The main effect of a large impact on the core dynamo is due to direct heating of the core, and diapir merging has almost negligible effects. The spreading of the merging diapir as a hot and thin layer surrounding the preexisting Martian core can only delay the initiation of the dynamo for a very short time.

In our models, the impact angle is vertical which maximizes the postimpact melt volume, hence the segregation volume in which the metal/silicates separation occurs [Pierazzo and Melosh, 2000]. Oblique impacts could influence the core merging dynamics by reducing the initial position of the metallic diapir and modifying the shape of the postimpact mantle thermal anomaly. The asymmetry introduced in the thermal state of the mantle by an oblique impact cannot be included in our current axisymmetric model and has to be modeled in a 3-D model in the future. However, an oblique impact has less effect on the cessation and reactivation of the core dynamo, except for the amount of impact heating of the core. This is because the core stratifies almost immediately after the impact and retains a spherically symmetric temperature distribution.

## Acknowledgments

J.M. is funded by Agence Nationale de la Recherche (Accretis decision ANR-10-PDOC-001-01), and J.A.H. is supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada. We also thank the two anonymous reviewers and the Associate Editor for valuable comments, which helped to improve the manuscript significantly.

## References

* Ahrens, T. J., and M. L. Johnson (1995), Shock wave data for rocks, in *Mineral Physics and Crystallography, A Handbook of Physical Constants*, vol. 3, edited by T. J. Ahrens, pp. 35–44, Amer. Geophys. Union, Washington, D. C.

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

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

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

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

* Arkani-Hamed, J. (2012), Life of the Martian dynamo, *Phys. Earth Planet. Inter.*, 196-197, 83–96.

* Arkani-Hamed, J., and P. Olson (2010a), Giant impact stratification of the Martian core, *Geophys. Res. Lett.*, 37, L02201, doi:10.1029/2009GL041417.

* Arkani-Hamed, J., and P. Olson (2010b), Giant impacts, core stratification, and failure of the Martian dynamo, *J. Geophys. Res.*, 115, E07012, doi:10.1029/2010JE003579.

* Asphaug, E. (2010), Similar-sized collisions and the diversity of planets, *Chemie der Erde / Geochem.*, 70, 199–219, doi:10.1016/j.chemer.2010.01.004.

* Bagdassarov, N., G. Solferino, G. J. Golabek, and M. W. Schmidt (2009), Centrifuge assisted percolation of Fe–S melts in partially molten peridotite: Time constraints for planetary core formation, *Earth Planet. Sci. Lett.*, 288(1–2), 84–95, doi:10.1016/j.epsl.2009.09.010.
Batchelor, G. (1967), *An Introduction to Fluid Dynamics*, Cambridge Univ. Press, Cambridge, U. K.

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

Bills, B. G., G. A. Neumann, D. E. Smith, and M. T. Zuber (2005), Improved estimate of tidal dissipation within Mars from MOLA observations of the shadow of Phobos, *J. Geophys. Res.*, **110**, E07004, doi:10.1029/2004JE002376.

Breuer, D., S. Labrosse, and T. Spohn (2010), Thermal evolution and magnetic field generation in terrestrial planets and satellites, *Space Sci. Rev.*, **152**, 449–500.

Cameron, A. G. W., and W. Benz (1991), The origin of the Moon and the single impact hypothesis IV, *Icarus*, **92**, 204–216.

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

Chambers, J. E. (2004), Planetary accretion in the inner solar system, *Earth Planet. Sci. Lett.*, **223**, 241–252.

Coradini, A., C. Federico, and P. Lanciano (1983), Earth and Mars—Early thermal profiles, *Phys. Earth Planet. Inter.*, **31**, 145–160.

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

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 Planet. Sci. Lett.*, **295**, 177–186.

Dauphas, N., and A. Pourmand (2011), Hf-W-Th evidence for rapid growth of Mars and its status as a planetary embryo, *Nature*, **473**, 489–493.

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.

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

Dziewonski, A. M., and D. L. Anderson (1981), Preliminary reference Earth model, *Phys. Earth Planet. Inter.*, **25**, 297–356.

Frey, H. V., J. H. Roark, K. M. Shockey, E. L. Frey, and S. E. H. Sakimoto (2002), Ancient lowlands on Mars, *Geophys. Res. Lett.*, **29**(10), 1384, doi:10.1029/2001GL013832.

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.

Ghods, A., and J. Arkani-Hamed (2007), Impact-induced convection as the main mechanism for formation of mare basalts, *J. Geophys. Res.*, **112**, E03005, doi:10.1029/2006JE002709.

Ghods, A., and J. Arkani-Hamed (2011), Effects of the Borealis impact on the mantle dynamics of Mars, *Phys. Earth Planet. Inter.*, **188**, 37–46.

Golabek, G. J., T. V. Gerya, B. J. P. Kaus, R. Ziethe, and P. J. Tackley (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.

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

Hofmeister, A. M. (1999), Mantle values of thermal conductivity and the geotherm from phonon lifetimes, *Science*, **283**, 1699–1706.

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.

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

Jellinek, A. M., H. M. Gonnermann, and M. A. Richards (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, doi:10.1016/S0012-821X(02)01070-1.

Jutzi, M., and E. Asphaug (2011), Forming the lunar farside highlands by accretion of a companion moon, *Nature*, **476**, 69–72.

Katz, R. F., M. Spiegelman, and C. H. Langmuir (2003), A new parameterization of hydrous mantle melting, *Geochem. Geophys. Geosyst.*, **4**(9), 1073, doi:10.1029/2002GC000433.

Ke, Y., and V. S. Solomatov (2009), Coupled core-mantle thermal evolution of early Mars, *J. Geophys. Res.*, **114**, E07004, doi:10.1029/2008JE003291.

Kleine, T., and J. F. Rudge (2011), Chronometry of meteorites and the formation of the Earth and Moon, *Elements*, **7**, 41–46.

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.

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

Laney, C. B. (1998), *Computational Gas Dynamics*, Cambridge University Press, Cambridge.

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

Louzada, K. L., and S. T. Stewart (2009), Effects of planet curvature and crust on the shock pressure field around impact basins, *Geophys. Res. Lett.*, **36**, L15203, doi:10.1029/2009GL037869.

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

McGill, G. E., and A. M. Dimitriou (1990), Origin of the Martian global dichotomy by crustal thinning in the Late Noachian or early Hesperian, *J. Geophys. Res.*, **95**, 12,595–12,605.

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

Melosh, H. J. (1989), *Impact Cratering: A Geologic Process*, Oxford Univ. Press, New York.

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

Monteux, J., N. Coltice, F. Dubuffet, and Y. Ricard (2007), Thermo-mechanical adjustment after impacts during planetary growth, *Geophys. Res. Lett.*, **34**, L24201, doi:10.1029/2007GL031635.

Monteux, J., Y. Ricard, N. Coltice, F. Dubuffet, and M. Ulvrová (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 (2011), 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., A. M. Jellinek, and C. L. Johnson (2013), Dynamics of core merging after a mega-impact with applications to Mars' early dynamo, *Icarus*, **226**, 20–32.

Morris, S. (1982), The effects of a strongly temperature-dependent viscosity on slow flow past a hot sphere, *J. Fluid Mech.*, **124**, 1–26, doi:10.1017/S0022112082002389.
Neukum, G., and D. V. Wise (1976), Mars: A standard crater curve and possible new time scale, *Science*, **194**, 1381–1387.

Neumann, W., D. Breuer, and T. Spohn (2012), Differentiation and core formation in accreting planetesimals, *Astron. Astrophys.*, **543**, A141, 23.

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

Nimmo, F., S. D. Hart, D. G. Korycansky, and C. B. Agnor (2008), Implications of an impact origin for the Martian hemispheric dichotomy, *Nature*, **453**, 1220–1223, doi:10.1038/nature07025.

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

Pierazzo, E., and H. J. Melosh (2000), Melt production in oblique impacts, *Icarus*, **145**, 252–261.

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

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

Reese, C. C., and V. S. Solomatov (2006), Fluid dynamics of local Martian magma oceans, *Icarus*, **184**, 102–120.

Reese, C. C., and V. S. Solomatov (2010), Early Martian dynamo generation due to giant impacts, *Icarus*, **207**, 82–97, doi:10.1016/j. icarus.2009.10.016.

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.

Roberts, J. H., and J. Arkani-Hamed (2012), Impact-induced mantle dynamics on Mars, *Icarus*, **218**, 278–289, doi:10.1016/j.icarus.2011.11.038.

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

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

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

Rubie, D. C., H. J. Melosh, J. E. Reid, C. Liebske, and K. Righter (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 einem Medium, *Bull. Acad. Sci. Cracovie*, **1**, 40–46.

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. J. Tackley (2008), Dynamics of core formation and equilibration by negative diapirism, *Geochem. Geophys. Geosyst.*, **9**, Q06011, doi:10.1029/2007GC001896.

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

Schubert, G., D. L. Turcotte, and P. Olson (2001), *Mantle convection in the Earth and planets*, Cambridge University Press, New York.

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

Shahnas, H., and J. Arkani-Hamed (2007), Viscous and impact demagnetization of Martian crust, *J. Geophys. Res.*, **112**, E02009, doi:10.1029/ 2005JE002424.

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

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

Solomatov, V. S. (2000), Fluid dynamics of a terrestrial magma ocean, in *Origin of the earth and moon*, edited by R. M. Canup et al., pp. 323–338, University of Arizona Press, Tucson.

Solomon, S. C., et al. (2005), New perspectives on ancient Mars, *Science*, **307**, 1214–1220.

Š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.

Stacey, F. D., and O. L. Anderson (2001), Electrical and thermal conductivities of Fe-Ni-Si alloy under core conditions, *Phys. Earth Planet. Inter.*, **124**, 153–162.

Stevenson, D. J. (1989), Formation and early evolution of the Earth, in *Mantle Convection and Plate Tectonics*, edited by W. R. Peltier, pp. 818–868 Routledge, Abingdon, U. K.

Stevenson, D. J. (1990), Fluid dynamics of core formation, in *Origin of the Earth*, edited by H. E. Newsom and J. H. Jones, pp. 231–249, Oxford Univ, New York.

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

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

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.

Trunin, R. F. (2001), Experimental data on shock compression and adiabatic expansion of condensed matter, *Academic Publication*, RFNC-VNIIEF, SAROV.

Ulvrová, M., N.oltice, Y. Ricard, S. Labrosse, F. Dubuffet, J. Vélimský, 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.

Vocadlo, L., D. Alfè, G. D. Price, and M. J. Gillan (2000), First principles calculations on the diffusivity and viscosity of liquid Fe–S at experi- mentally accessible conditions, *Phys. Earth Planet. Inter.*, **120**, 145–152.

Walte, N. P., J. K. Becker, P. D. Bons, D. C. Rubie, and D. J. Frost (2007), Liquid-distribution and attainment of textural equilibrium in a partially- molten crystalline system with a high-dihedral-angle liquid phase, *Earth Planet. Sci. Lett.*, **262**(3-4), 517–532.

Walte, N. P., D. C. Rubie, P. D. Bons, and D. J. Frost (2011), Deformation of a crystalline aggregate with a small percentage of high-dihedral-angle liquid; implications for core-mantle differentiation during planetary formation, *Earth Planet. Sci. Lett.*, **305**, 1–2.

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

Weinberg, R. F., and Y. Podladchikov (1994), Diapiric ascent of magmas through power law crust and mantle, *J. Geophys. Res.*, **99**, 9543–9559, doi:10.1029/93JB03461.

Wilhelms, D. E., and S. W. Squyres (1984), The Martian hemispheric dichotomy may be due to a giant impact, *Nature*, **309**, 138–140, doi:10.1038/309138a0.
<figure><img src="image_207.png" /></figure>

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

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

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

Zhong, S. J., and M. T. Zuber (2001), Degree-1 mantle convection and the crustal dichotomy on Mars, *Earth Planet. Sci. Lett.*, **189**, 75–84.

Ziethe, R., and T. Spohn (2007), Two-dimensional stokes flow around a heated cylinder: A possible application for diapirs in the mantle, *J. Geophys. Res.*, **112**, B09403, doi:10.1029/2006JB004789.
<figure><img src="image_208.png" /></figure>

ELSEVIER

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

# Giant impacts, heterogeneous mantle heating and a past hemispheric dynamo on Mars

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

CrossMark

Julien Monteux<sup>a,b,c,*</sup>, Hagay Amit<sup>a</sup>, Gaël Choblet<sup>a</sup>, Benoit Langlais<sup>a</sup>, Gabriel Tobie<sup>a</sup>

## ARTICLE INFO

**Article history:**

Received 22 July 2014

Received in revised form 12 November 2014

Accepted 16 December 2014

Available online 12 January 2015

**Keywords:**

Mars

Dynamo

Magnetic field

Core-mantle boundary

Impact

Heat flux

## ABSTRACT

The martian surface exhibits a strong dichotomy in elevation, crustal thickness and magnetization between the southern and northern hemispheres. A giant impact has been proposed as an explanation for the formation of the Northern Lowlands on Mars. Such an impact probably led to strong and deep mantle heating which may have had implications on the magnetic evolution of the planet. We model the effects of such an impact on the martian magnetic field by imposing an impact induced thermal heterogeneity, and the subsequent heat flux heterogeneity, on the martian core-mantle boundary (CMB). The CMB heat flux lateral variations as well as the reduction in the mean CMB heat flux are determined by the size and geographic location of the impactor. A polar impactor leads to a north-south hemispheric magnetic dichotomy that is stronger than an east-west dichotomy created by an equatorial impactor. The amplitude of the hemispheric magnetic dichotomy is mostly controlled by the horizontal Rayleigh number $Ra_h$ which represents the vigor of the convection driven by the lateral variations of the CMB heat flux. We show that, for a given $Ra_h$, an impact induced CMB heat flux heterogeneity is more efficient than a synthetic degree-1 CMB heat flux heterogeneity in generating strong hemispheric magnetic dichotomies. Large $Ra_h$ values are needed to get a dichotomy as strong as the observed one, favoring a reversing paleo-dynamo for Mars. Our results imply that an impactor radius of ~1000 km could have recorded the magnetic dichotomy observed in the martian crustal field only if very rapid post-impact magma cooling took place.

© 2014 Elsevier B.V. All rights reserved.

## 1. Introduction

Giant impacts have strongly influenced the internal structure and dynamics of the terrestrial planets during the primordial stages of their evolutions (Hartmann and Davis, 1975; Benz et al., 1988; Asphaug et al., 2006; Andrews-Hanna et al., 2008; Marinova et al., 2008; Nimmo et al., 2008; Jutzi and Asphaug, 2011). These events are plausible explanations for remarkable features of the solar system such as the small volume of Mercury's mantle relative to its core (Benz et al., 1988; Gladman and Coffey, 2009), the Earth-Moon system (Canup, 2004) and the topographic martian and lunar hemispheric dichotomies (Marinova et al., 2008; Nimmo et al., 2008; Jutzi and Asphaug, 2011). Giant

impacts have also been invoked to explain the initiation or cessation of the dynamos of the terrestrial planets and moons (Roberts et al., 2009; Arkani-Hamed and Olson, 2010; Reese and Solomatov, 2010; Monteux et al., 2013; Monteux and Arkani-Hamed, 2014). In these models, the impactors' radii typically range between 100 and 1000 km. These impacts deliver a large amount of heat to the deep mantle, which is likely to strongly affect the efficiency of core cooling and in turn the dynamo activity. Although there is a higher probability that a giant impact will fall on low-latitudes of the planetary surface (Le Feuvre and Wieczorek, 2011), true polar wander events can ultimately place the resulting thermal anomaly at high-latitudes of the Core Mantle Boundary (CMB). Moreover, large impacts could be responsible for significant resurfacing and reset the magnetization of the pre-impact material (Langlais and Thébault, 2011; Lillis et al., 2013).

On Earth, the influence of lower mantle thermal heterogeneity on core magnetohydrodynamics has been extensively studied

<footnote><sup>a</sup>Laboratoire de Planétologie et Géodynamique, Université de Nantes, CNRS, UMR-6112, 2 rue de la Houssinière, 44322 Nantes Cedex, France</footnote>

<footnote><sup>b</sup>ISTerre, Université de Grenoble 1, CNRS, Grenoble, France</footnote>

<footnote><sup>c</sup>Laboratoire Magmas et Volcans, Université Blaise Pascal – CNRS – IRD, OPGC, 5 rue Kessler, 63038 Clermont Ferrand, France</footnote>

<footnote>* Corresponding author at: Laboratoire Magmas et Volcans, Université Blaise Pascal – CNRS – IRD, OPGC, 5 rue Kessler, 63038 Clermont Ferrand, France.
E-mail address: j.monteux@opgc.univ-bpclermont.fr (J. Monteux).</footnote>
using numerical dynamos with imposed non-uniform outer boundary conditions. It has been shown that heterogeneous CMB heat flux causes a deviation from axisymmetry in the core flow (Aubert et al., 2007), in the time-average paleomagnetic field (Olson and Christensen, 2002) and in locations of intense magnetic flux patches on millennial time-scales (Bloxham, 2002; Amit et al., 2010). It may also explain the emergence of intense magnetic flux patches in the equatorial region (Amit and Choblet, 2012) and may even yield field locking (Gubbins et al., 2007; Willis et al., 2007). Heterogeneous CMB heat flux may also recover the lateral variations in the inner-core boundary seismic properties (Aubert et al., 2008; Amit and Choblet, 2009). Finally, reversal frequency and the trajectory of the paleomagnetic dipole axis during reversals may also be governed by the heterogeneous lower mantle (Glatzmaier et al., 1999; Kutzner and Christensen, 2004; Olson et al., 2010, 2013; Olson and Amit, 2014).

Heterogeneous mantle control has also been proposed to explain some features of planetary magnetic fields. Cao et al. (2014) found that high equatorial CMB heat flux breaks the core flow symmetry and produces north-south asymmetric magnetic fields which may explain the observed field of Mercury (Anderson et al., 2012). Stanley (2010) argued that temperature differences in the surrounding envelope of the convective zone of Saturn axisymmetrize its magnetic field. It has also been proposed that CMB heterogeneity may have controlled the shape of the current Martian magnetic field (Stanley et al., 2008). Mars is characterized by a striking magnetic field dichotomy, which is correlated with the topographic dichotomy. The Northern Lowlands are mostly devoid of significant magnetic fields. In contrast the southern highlands exhibit large and in some places intense magnetic field anomalies, up to 1500 nT at 90 km altitude as measured by Mars Global Surveyor (Acuña et al., 1998). This is two orders of magnitude larger than the crustal magnetic field on Earth. In terms of magnetized material, this suggests a thick (40 km) and intensely magnetized (up to 12 A/m) lithosphere to produce the observed magnetic field (Langlais et al., 2004), or any combination of a thinner lithosphere and a more intense magnetization (e.g., Parker, 2003).

The martian magnetic dichotomy can be explained using two end-members scenarios. In the external scenario, the dynamo was equally strong in both hemispheres, and the resulting magnetization was equally strong in both hemispheres. Then the magnetization of the northern hemisphere was removed or erased after the dynamo cessation, e.g., by a giant impact (Nimmo et al., 2008) or volcanic activity (Lillis et al., 2008). Alternatively a significant magnetization was never recorded in the northern hemisphere because surface conditions, lithological or alteration processes were different from those in the southern hemisphere (Rochette, 2006; Quesnel et al., 2009; Chassefière et al., 2013). In the internal scenario, the magnetization is strong only in the southern hemisphere because the dynamo was hemispheric to begin with (Langlais and Amit, 2008; Stanley et al., 2008; Amit et al., 2011).

Such an hemispheric dynamo could have been driven by CMB heat flux heterogeneity possibly caused by a very large-scale mantle convection pattern (Harder and Christensen, 1996; Zhong and Zuber, 2001; Elkins-Tanton et al., 2003, 2005; Ke and Solomatov, 2006; Roberts and Zhong, 2006) or by a giant impact (Roberts et al., 2009). In this study we propose a model for the magnetic field dichotomy in which the dynamo hemisphericity (internal origin) is related to a large impact (external origin). For that purpose, we model heterogeneous CMB heat flux resulting from giant impact heating and investigate its influence on the core dynamo by imposing it as a static, laterally-varying outer boundary condition on numerical dynamo models. In this approach, the CMB heat flux pattern and amplitude, as well as the reduction in the mean

heat flux with respect to a reference pre-impact value, are determined by the impactor size, using a synthetic description of the impact heating zone. In Section 2 we describe our method. The results are presented in Section 3. Discussion, post-impact time evolution and applicability of our results to Mars are given in Section 4. Conclusions and possible planetary applications are highlighted in Section 5.

## 2. Method

### 2.1. Impact heating at the Core Mantle Boundary

Large impacts brought to Mars a formidable amount of energy that is a function of the impactor mass and velocity, the latter strongly depending on the impacted planet radius $R$. After a large collision on a Mars-size body, a significant fraction of this energy is deeply buried as heat within the mantle and leads to a local temperature increase $\Delta T_0$ below the impact site. The size and the shape of the post-impact thermal anomaly depend on several parameters such as the size of the impactor, the impact velocity and angle, and the structure of the martian mantle. Increasing the size of the impactor leads to an increase of the heated volume while increasing the impact angle from 0 (head-on impact) to larger values (oblique impacts) reduces the maximal depth reached by the post-impact thermal anomaly (Pierazzo et al., 1997; Pierazzo and Melosh, 2000). Here for simplicity, we consider that the volume of the thermal anomaly only scales with the size of the impactor and we consider the case of a head-on impact. Hence, the post-impact thermal anomaly in our models is approximately uniform within a spherical volume (termed isobaric core) with radius $R_{ic}$ that is 1 to 1.44 times larger than the radius of the impactor $R_{imp}$ (Pierazzo et al., 1997; Senshu et al., 2002; Monteux et al., 2013).

On Mars, the impactor size invoked to explain the topographic dichotomy ranges between 320 and 1350 km (Marinova et al., 2008; Nimmo et al., 2008). This has to be compared to the size of the martian mantle. Based on solar tidal deformations, the martian core radius has been estimated between 1520 and 1840 km (Yoder et al., 2003). For simplicity, we assume a core radius of 1700 km, which implies a mantle thickness of about 1700 km. Hence, considering that $R_{ic} = 1.44R_{imp}$, the post-impact spherical thermal anomaly is likely to overlap the CMB for $R_{imp} > 500$ km. For an impactor radius of $R_{imp} = 1200$ km, the disruption of the impacted planet will only occur when the impact velocity reaches values of ~100 km/s which is much larger than the impact velocity $v_{imp}$ considered here ($v_{imp} = 5$ km/s) (Tonks and Melosh, 1992; Reese et al., 2010). In our models, we consider that the impactor radius ranges between 600 and 1000 km bearing in mind that larger impactors with larger impact angles could have similar thermal consequences at the CMB (Pierazzo et al., 1997; Pierazzo and Melosh, 2000).

As the volume of the isobaric core is governed by the size of the impactor, the magnitude of the temperature increase can be directly related to the impact velocity. Making the conservative hypothesis that the impact velocity is close to the martian escape velocity and that the volume of the isothermal sphere is 3 times larger than the impactor (Senshu et al., 2002; Monteux et al., 2013), the energy balance accounting for heating and melting of both the impactor and impacted material may lead to a uniform spherical temperature increase of ~400 K in the martian mantle (Monteux et al., 2013). Away from the isothermal sphere, the temperature decreases rapidly with distance $r$ as $(R_{ic}/r)^m$ with $m$ typically ranging between 4 and 5 (Senshu et al., 2002; Monteux et al., 2007).

Geochemical evidence and crater densities indicate that the martian topographic dichotomy could have formed within the first 50 Myr of Solar System formation and that the martian northern
hemisphere has been low and stable for nearly all of Mars' history (Zuber, 2001; Frey et al., 2002; Solomon et al., 2005; Marinova et al., 2008). Hence, the impact-driven temperature increase is superimposed to the pre-impact thermal state of the martian mantle that strongly depends on the short-lived radiogenic heating, the accretion processes and the dissipation of gravitational energy during the core formation (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. For simplicity, we assume here a 1D radially dependent pre-impact mantle temperature field. The choice of this specific temperature field is not crucial as long as the impact heating of the mantle is predominant. As shown later, the main parameter affecting the dynamo is the amplitude of the heat-flux heterogeneity at the CMB. Before the impact heating, we consider a simplified temperature profile as in Monteux et al. (2013) with a CMB temperature of $T = 2000$ K and a convective mantle temperature of $\sim 1600$ K (Roberts and Arkani-Hamed, 2012). It should be noted that in reality, the pre-impact thermal state of the martian mantle was probably much more complicated than the one used in our models, with lateral heterogeneities as well as radial variations (including thermal boundary layers and pressure dependence). For simplicity, we do not consider the changes of mantle properties with depth such as the pressure increase and the corresponding adiabatic heating. Because we consider here that the mantle temperature above the CMB is uniform and equal to 1600 K while the core temperature is 2000 K, adding more complexity should slightly decrease the amplitude of the heat flux heterogeneity. However, since the early martian thermal state is poorly constrained, our simple model may be considered as a reasonable first step to understand the influence of giant impacts on planetary dynamos.

The superposition of the large impact-driven temperature increase leads to a significant perturbation of the pre-impact homogeneous CMB mean heat flux $q_0^h$. To obtain the post-impact heat flux at the martian CMB for a given impactor radius, we used an impact heating model similar to the one described in Monteux et al. (2013). A uniform spherical temperature anomaly rapidly decreasing with distance is superimposed on the martian pre-impact temperature field. The impact heating is followed by a thermal re-adjustment within a characteristic time that is governed by the rheology of the mantle surrounding the impact-induced thermal anomaly. This characteristic time is $\tau_{imp} \sim 10$ Myr (Monteux et al., 2007). The characteristic magnetic diffusion time is $\tau_\lambda = r_0^2/\lambda$ where $r_0$ is the core radius and $\lambda$ is the magnetic diffusivity (e.g., Bloxham and Jackson, 1991). For Mars $r_0 = 1700$ km and combined with updated estimates of the electrical conductivity of molten iron in Earth's core conditions of $\lambda \sim 0.5 \text{ m}^2/\text{s}$ (Pozzo et al., 2012), these estimates give $\tau_\lambda \sim 180$ kyr. Since $\tau_{imp} \gg \tau_\lambda$, the post impact mantle temperature field may be considered constant in our numerical dynamo models.

The shock wave also leads to a temperature increase within the core of the impacted planet, much stronger in the region directly beneath the impact site. The low-viscosity rotating liquid core cannot sustain lateral variations of temperature and the core overturns, resulting in a stably stratified temperature which increases with radius. In the case of a homogeneous CMB heat flux, the thermal stratification occurs within a few kyr (Arkani-Hamed and Olson, 2010). This stratification kills the possible pre-existing core convection, and hence the core dynamo. Then it can take up to several tens of Myr to remove this stratification by conductive core heat loss (Arkani-Hamed and Olson, 2010). However, a significant fraction of the mantle right above the CMB may have experienced melting which facilitated the core cooling. The time needed to remove the core excess heat by convection is $10^3–10^4$ years when

considering a molten layer above the CMB (Monteux et al., 2011). It is therefore likely that after the impact the martian dynamo died during the time needed to generate and remove the post-impact core stratification (i.e., $\sim 10^3–10^4$ years). Then, the dynamo probably re-started, while the impact-driven lower mantle anomaly was still in place (during $\tau_{imp} \sim 10$ Myr). For simplicity, we do not consider here the core impact heating and the subsequent rapid thermal readjustment. Our dynamo models are convectively unstable throughout the shell, corresponding to the state of the system after core stratification has been removed and while the CMB heterogeneity was still in place. We discuss this aspect in the conclusion section.

As the post-impact mantle temperature reaches the temperature of the core, the heat flux $q$ is nearly zero where the isobaric core overlaps with the CMB. Away from the isobaric core and along the CMB, the heat flux increases rapidly to its pre-impact mean value $q_0^h$ (Fig. 1). The reduction in the mean CMB heat flux due to the impact corresponds to the relative CMB surface that is heated by the impactor and is defined by $q_0^h$ as

$$q_0^r = \frac{q_0^h - q_0}{q_0^h} \qquad (1)$$

where $q_0$ is the post-impact mean heat flux. The amplitude of the heat flux heterogeneity is commonly given by $q^*$ (Olson and Christensen, 2002) with

$$q^* = \frac{q_{max} - q_{min}}{2q_0} \qquad (2)$$

where $q_{max}$ and $q_{min}$ are the maximal and minimal values of the post-impact heat flux, respectively. In Fig. 1, we show $q_0^r$ and $q^*$ obtained from our post-impact heating model as a function of the impactor size. As the impactor size increases, the extension of the isobaric core on the CMB increases, i.e., $q_0$ decreases and therefore $q_0^r$ and $q^*$ increase. An impactor with a radius smaller than $R_{imp} \sim 500$ km has a negligible effect on the heated surface of the CMB. For small impactors a small portion of the CMB is heated so $q_0 \approx q_0^h$, and since $q_{max} = q_0^h$ and $q_{min} = 0$ always hold, for small impactors $q^* \approx 0.5$. An impactor radius of $R_{imp} \sim 1000$ km decreases the mean CMB heat flux by 26.5% and produces an heterogeneity with amplitude $q^* \sim 0.67$. In our models, we limit the impactor radius to 1000 km and consider a head-on impact which may represent cases with larger impactors and smaller impact angles (Pierazzo et al., 1997; Pierazzo and Melosh, 2000).

Next we expand the impact-driven CMB heat flux pattern in terms of spherical harmonic coefficients (Fig. 2). In order to use it as an outer boundary condition for heat flux in our dynamo models, we performed a spherical harmonic expansion truncated at $l_{max} = 20$ (see Fig. 2 for $R_{imp} = 800$ km). First we fit the CMB heat flux with the following analytical expression in terms of the angular distance $\omega$ from the impactor's center:

$$q_{\omega} = \exp \left[ -\frac{1}{n} \left( \frac{\omega}{\omega_0} \right)^{n} \right] \qquad (3)$$

The Gauss-like function (3) avoids undesirable Gibbs effects associated with discontinuous gradients on the edge of the imprint of the isobaric core on the CMB. The best fit parameters found for the three impactor radii are $\omega_0 = 15^\circ$ and $n = 3$ for $R_{imp} = 600$ km, $\omega_0 = 30^\circ$ and $n = 7$ for $R_{imp} = 800$ km and $\omega_0 = 48^\circ$ and $n = 7$ for $R_{imp} = 1000$ km.

## 2.2. Numerical dynamo models

We solve the set of self-consistent non-dimensional Boussinesq magnetohydrodynamics equations for dynamo action due to thermal convection of an electrically conducting fluid in a rotating
<figure><img src="image_211.png" /><figcaption>Fig. 1. (a) martian mantle post-impact temperature increase for R<sub>imp</sub> = 800 km. (b) Conventional amplitude of heat flux heterogeneity on the CMB q' (black line, Eq. (2)) and the mean heat flux reduction q'<sub>0</sub> (red line, Eq. (1)) vs. impactor size R<sub>imp</sub>. In the grey domain (R<sub>imp</sub> < 500 km), the impactor is too small for isobaric core to reach the CMB. (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_212.png" /><figcaption>Fig. 2. (a) Post-impact CMB heat flux (black line), analytical fit (red line) and spherical harmonic expansion (green line) vs. angular distance from the center of the impact driven heated area. (b) Imposed CMB heat flux anomalies resulting from an impactor of radius 800 km falling on the north pole (left) or the equator (right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

spherical shell (for governing equations and more details see Amit et al., 2011). We use the code MAGIC by Wicht (2002). We analyze numerical dynamos with rigid insulating boundary conditions. The models differ in the imposed outer boundary heat flux pattern and the amplitude of its variation. A summary of model parameters, outer boundary heat flux patterns and main results is given in Table 1.

It is likely that Mars has no solid inner core up to the present day (e.g., Schubert and Spohn, 1990; Breuer et al., 2010) and, as a consequence, convection in the early martian dynamo was purely thermal (Amit et al., 2011; Dietrich and Wicht, 2013), driven by secular cooling and perhaps by radioactive heating. This convection mode is highly sensitive to CMB heat flux heterogeneity and may thus break internal dynamo symmetries with relatively moderate heterogeneity amplitudes (Hori et al., 2014). Due to numerical singularity at the center of the planet, we retain in our dynamo models a small inner core with a radius r<sub>i</sub>/r<sub>o</sub> = 0.2 of the outer core radius. We impose zero heat flux on the inner boundary so that the inner core is convectively passive. Overall, results by Aubert et al. (2009) and Hori et al. (2010) suggest that such a relatively small and passive inner core has little effect on the dynamo models.

Four internal non-dimensional parameters control the dynamo action. The heat flux Rayleigh number (Olson and Christensen, 2002) represents the strength of buoyancy force driving the convection relative to retarding forces

$$Ra = \frac{\alpha g_0 q_0 D^4}{k \kappa \nu} \qquad (4)$$

where α is thermal expansivity, g<sub>0</sub> is gravitational acceleration on the outer boundary at radius r<sub>o</sub>, q<sub>0</sub> is the mean heat flux across

the outer boundary, D is shell thickness, k is thermal conductivity, κ is thermal diffusivity and ν is kinematic viscosity. The Ekman number represents the ratio of viscous and Coriolis forces

$$E = \frac{v}{\Omega D^2} \qquad (5)$$

The Prandtl number is the ratio of kinematic viscosity to thermal diffusivity

$$Pr = \frac{v}{\kappa} \qquad (6)$$

and the magnetic Prandtl number is the ratio of kinematic viscosity to magnetic diffusivity λ

$$Pm = \frac{v}{\lambda} \qquad (7)$$

A fifth non-dimensional number is the amplitude of the outer boundary heat flux heterogeneity which is expressed by q* (see Eq. (2)).

In all cases a volumetric homogeneous heat source ϵ compensates for the loss of heat through the outer boundary according to

$$-4\pi r_0^2 Pr \left[ \frac{\partial T}{\partial r}(r_0) \right] = \frac{4}{3} \pi (r_0^3 - r_i^3) \epsilon \qquad (8)$$

where [...] denotes averaging over the outer boundary surface S. In terms of the non-dimensional variables $[\frac{\partial T}{\partial r}(r_0)] = 1$, so for r<sub>i</sub>/r<sub>o</sub> = 0.2 (the geometry used in the study) the source term is ϵ ≃ 2.42. We use moderate amplitudes of CMB heat flux heterogeneity to avoid violation of the Boussinesq approximation on which the dynamo code relies.
<caption>**Table 1**</caption>

Summary of numerical dynamo models parameters and resulted magnetic hemispheric dichotomies. In all cases $r_i/r_0 = 0.2$, $Pr = 1$ and $Pm = 3$. The upper part includes reversing dynamos with $E = 3 \cdot 10^{-4}$; The lower part includes non-reversing dynamos with $E = 1 \cdot 10^{-4}$. The imposed CMB heat flux patterns are either degree-1 single harmonics ($Y_1^0$ or $Y_1^1$) denoted as cases Y, or impact driven centered at the pole or at the equator denoted as cases I. All Y cases are from Amit et al. (2011) except cases Y2 and Y3 which are from this study. The radius of the impactor is $R_{imp}$. The conventional amplitude of the heat flux anomaly $q^*$ is defined as the ratio of the peak-to-peak difference to twice the mean (Olson and Christensen, 2002). The reduction in the mean heat flux due to the impact is $q_0^*$. $Rm$ is the magnetic Reynolds number based on the rms velocity in the volume of the shell. SN and EW denote south-north and east-west, 'cnt' and 'rnd' subscripts denote continuous and random crust formation. East-west ratios are calculated with respect to the hemisphere centered at longitude 0° (center of large heat flux anomaly for $Y_1^1$ and equatorial I cases) to the hemisphere centered at longitude 180°. For reversing cases only ratios based on random crust formation are given, for $Y_1^0$ or polar I cases only south-north ratios are given, for $Y_1^1$ or equatorial I cases only east-west ratios are given. As reference cases, the homogeneous CMB heat flux models from Amit et al. (2011) give SN and EW values of unity. The average chron duration $\tau_{ch}$ is given in kyrs for the reversing cases.

<table><thead><tr><th>Case</th><th>Pattern</th><th>R<sub>imp</sub> [km]</th><th>Ra</th><th>q<sup>*</sup></th><th>q<sup>*</sup><sub>0</sub></th><th>Rm</th><th>SN<sub>cnt</sub></th><th>SN<sub>rnd</sub></th><th>EW<sub>cnt</sub></th><th>EW<sub>rnd</sub></th><th>τ<sub>ch</sub></th></tr></thead><tbody><tr><td>I1</td><td>Polar</td><td>600</td><td>2.90 · 10<sup>7</sup></td><td>0.51</td><td>0.04</td><td>386</td><td>–</td><td>1.18</td><td>–</td><td>–</td><td>40.9</td></tr><tr><td>I2</td><td>Polar</td><td>800</td><td>2.69 · 10<sup>7</sup></td><td>0.56</td><td>0.13</td><td>423</td><td>–</td><td>1.54</td><td>–</td><td>–</td><td>22.6</td></tr><tr><td>I3</td><td>Polar</td><td>1000</td><td>2.21 · 10<sup>7</sup></td><td>0.66</td><td>0.26</td><td>439</td><td>–</td><td>2.25</td><td>–</td><td>–</td><td>18.3</td></tr><tr><td>Y1</td><td>Y<sub>1</sub><sup>0</sup></td><td>–</td><td>3 · 10<sup>7</sup></td><td>0.8</td><td>–</td><td>488</td><td>–</td><td>2.54</td><td>–</td><td>–</td><td>12.2</td></tr><tr><td>I4</td><td>Polar</td><td>1000</td><td>4 · 10<sup>7</sup></td><td>0.66</td><td>0.26</td><td>572</td><td>–</td><td>2.58</td><td>–</td><td>–</td><td>15.4</td></tr><tr><td>Y2</td><td>Y<sub>1</sub><sup>0</sup></td><td>–</td><td>2 · 10<sup>7</sup></td><td>0.8</td><td>–</td><td>412</td><td>–</td><td>2.35</td><td>–</td><td>–</td><td>14.9</td></tr><tr><td>I5</td><td>Equatorial</td><td>800</td><td>2.69 · 10<sup>7</sup></td><td>0.56</td><td>0.13</td><td>350</td><td>–</td><td>–</td><td>–</td><td>1.22</td><td>28.2</td></tr><tr><td>Y3</td><td>Y<sub>1</sub><sup>0</sup></td><td>–</td><td>1 · 10<sup>7</sup></td><td>0.8</td><td>–</td><td>303</td><td>–</td><td>1.68</td><td>–</td><td>–</td><td>26.3</td></tr><tr><td>Y4</td><td>Y<sub>1</sub><sup>1</sup></td><td>–</td><td>3 · 10<sup>7</sup></td><td>0.8</td><td>–</td><td>405</td><td>–</td><td>–</td><td>–</td><td>1.14</td><td>17.5</td></tr><tr><td>I6</td><td>Polar</td><td>800</td><td>8.97 · 10<sup>6</sup></td><td>0.56</td><td>0.13</td><td>99</td><td>1.12</td><td>1.13</td><td>–</td><td>–</td><td>∞</td></tr><tr><td>I7</td><td>Polar</td><td>1000</td><td>7.35 · 10<sup>6</sup></td><td>0.66</td><td>0.26</td><td>182</td><td>1.59</td><td>1.57</td><td>–</td><td>–</td><td>∞</td></tr><tr><td>Y5</td><td>Y<sub>1</sub><sup>0</sup></td><td>–</td><td>1 · 10<sup>7</sup></td><td>0.3</td><td>–</td><td>114</td><td>1.16</td><td>1.16</td><td>–</td><td>–</td><td>∞</td></tr><tr><td>Y6</td><td>Y<sub>1</sub><sup>0</sup></td><td>–</td><td>1 · 10<sup>7</sup></td><td>0.5</td><td>–</td><td>157</td><td>1.49</td><td>1.50</td><td>–</td><td>–</td><td>∞</td></tr><tr><td>I8</td><td>Equatorial</td><td>800</td><td>8.97 · 10<sup>6</sup></td><td>0.56</td><td>0.13</td><td>104</td><td>–</td><td>–</td><td>1.01</td><td>1.01</td><td>∞</td></tr><tr><td>Y7</td><td>Y<sub>1</sub><sup>1</sup></td><td>–</td><td>1 · 10<sup>7</sup></td><td>0.5</td><td>–</td><td>133</td><td>–</td><td>–</td><td>1.01</td><td>1.01</td><td>∞</td></tr></tbody></table>

Most dynamo models fall into two categories. In the first, the radial field on the CMB is dominated by an axial dipole component, but the field does not reverse. In the second, the field is multipolar and dipole reversals occur (Kutzner and Christensen, 2002). Earth-like models that are both dipole-dominated and reversing are only found in a narrow transitional regime of parameters space (Olson, 2007; Wicht et al., 2009). Following Amit et al. (2011), we consider cases from both regimes.

## 2.3. Hemispheric magnetic dichotomy monitoring

The kinetic energy of the impactor is dissipated as a result of the irreversible work done by shock waves in damaging crustal rocks as well as heating and melting the target material (Pierazzo et al., 1997; Senshu et al., 2002; Reese and Solomatov, 2006; Monteux et al., 2011). After the excavation process, a significant fraction of the material molten by impact is redistributed heterogeneously at the surface of the impacted planet (Marinova et al., 2008, 2011). The distribution of the molten material is governed by the impact parameters such as the impact velocity and angle. At low impact velocities (6–10 km/s) and oblique impact angles (30–60°), 50–70% of the impact-induced melt distribution might be contained within the area of impact and 25–30% might be deposited at the antipode of the impact site (Marinova et al., 2008). The impact-induced molten material might be redistributed over a thickness that ranges between 30 and 50 km (Marinova et al., 2008, 2011). The cooling and the crystallization of the molten material leads to the formation of the impact induced crust potentially recording the anomalous dynamo.

We consider two end member crust formation scenarios (Amit et al., 2011; Langlais and Thébault, 2011) as illustrated in Fig. 3. In the case of continuous and homogeneous crust formation, each part of the martian crust is formed by a large number of incremental and superposed additions (e.g., lava flows, sills, dykes) over an extended period of time. Each new layer records the magnetic field at its time of cooling below the Curie Temperature, and the present crustal field at a specific location results from the vertical superposition of the magnetization vectors of the various layers (Fig. 3, bot-

tom line). If this crust formation scenario occurs while the dynamo reverses (Fig. 3, bottom right), the present local crustal field would thus represent the intensity of a long-term time-average martian paleomagnetic field. In the extreme case of periodic inversions and layers with equal thicknesses, the superposed opposite sign magnetization vectors could eventually cancel each other. The other end-member model assumes a random crust formation where crustal units are formed in relatively rapid events (Fig. 3, top line). Individual crustal blocks are created randomly both in space and time and each block acquires a magnetization which only depends on the dynamo field at the time of cooling. The present local field in this scenario thus results from adjacent magnetization vectors, or in a probabilistic way, from the time-average of the paleofield intensity. It is likely that neither one of these end-member scenarios represent what actually occurred on Mars, and rather that the actual way in which the magnetized part of the martian crust formed may be intermediate between the two scenarios.

We follow the statistical measures proposed by Amit et al. (2011) corresponding to these two end-member crust formation scenarios. In the context of a continuous homogeneous crust formation scenario, we calculate the ratio of intensities of the time-average field at the planet surface

$$SN_{\text{cnt}} = [\langle |\vec{B}| \rangle]_{\text{sh}} / [\langle |\vec{B}| \rangle]_{\text{nh}} \quad (9)$$

where $SN$ denotes the ratio between the rms surface average in the southern and northern hemispheres and the subscript 'cnt' denotes continuous crust formation. In the context of a random crust formation scenario, we calculate the magnetic dichotomies based on the time-average of the magnetic field intensity

$$SN_{\text{rnd}} = [\langle |\vec{B}| \rangle]_{\text{sh}} / [\langle |\vec{B}| \rangle]_{\text{nh}} \quad (10)$$

where the subscript 'rnd' denotes random. Eqs. (9) and (10) are applied for the east-west dichotomies $EW_{cnt}$ and $EW_{rnd}$ by replacing the summations with the appropriate hemispheres.

The relevance of the continuous or random crust formation scenario depends on the relation between the cooling time of the Impact Induced Molten Material (IIMM) and the typical magnetic
<figure><img src="image_213.png" /><figcaption>**Fig. 3.** Schematic illustration of the random (top) and continuous (bottom) crust formation scenarios and the corresponding recorded magnetization without (left) or with (right) magnetic reversals. Blue/red denotes negative/positive radial magnetic paleofields, respectively, so alternating colors correspond to paleomagnetic reversals. Cubes represent a vertical cut through the martian crust that has recorded strong (deep blue or deep red) or weak (light red or light blue) magnetic paleofields. The weak magnetic paleofield is restricted to the impacted pole while the strong magnetic paleofield is restricted to the opposite one. The current magnetic field observed at the surface results from the vertically integrated magnetization over the cube. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

timescales. The cooling time of the IIMM strongly depends on the presence of an atmosphere that may prevent the IIMM from rapid cooling (e.g., Lebrun et al., 2013). In the case of molten material induced by a giant impact, a significant fraction of the atmosphere (if any) can be eroded from the impacted body (Shuvalov, 2009) which strongly enhances the excess heat removal and decreases the cooling and solidification times. Hence, the duration of the partially molten stage decreases from more than 1 million years with an atmosphere to ~1000 years when no atmosphere is present (Lebrun et al., 2013). The two crust formation scenarios give identical dichotomies for non-reversing dynamos (Amit et al., 2011). In the reversing case, if the crust was formed over a period much longer than a typical magnetic timescale, the continuous scenario is relevant (Dietrich and Wicht, 2013). Conversely, if the crust was formed faster than a typical magnetic timescale (i.e., without any atmosphere), the random scenario is relevant (see Fig. 3).

### 3. Influence of impact heating on core dynamo

Various CMB heat flux models corresponding to varying impactor sizes were examined. We considered two extreme geographic locations of the impactor, with a center falling on the geographical pole or on the equator. The impactor introduces warm material to the lower mantle, thus locally decreases the CMB heat flux. Apart from pattern dependence, the reduction in mean CMB heat flux (expressed by the Ra number compared to its reference homogeneous case) as well as the amplitude of the heat flux heterogeneity $q^*$ are all dependent in a self-consistent manner on the radius of the impactor. A larger impactor warms a larger part of the CMB, thus reduces the mean CMB heat flux $q_0$ more and produces larger $q^*$ (see Table 1).

Fig. 4 shows the time-average rms radial field at the CMB and the time-average field intensity at the surface of Mars for case I2 with an impactor of radius 800 km falling on the north geographic

pole. Here and elsewhere all magnetic field values are given in units of $\sqrt{\rho\mu_0\lambda^2\Omega}$ where $\rho$ is the fluid density and $\mu_0$ is permeability of free space. The impactor yields an hemispheric field similar to those obtained with synthetic $Y_1^0$ CMB heat flux patterns (Stanley et al., 2008; Amit et al., 2011; Dietrich and Wicht, 2013).

Fig. 5 shows the time-average zonal temperature and flow in case I2. The reduced CMB heat flux in the north pole region caused by the impactor results in a relatively warmer fluid there. The colder fluid in the more vigorously convecting southern hemisphere is associated with fluid downwelling at high-latitudes that concentrates the magnetic field. This produces a south-north magnetic hemispheric dichotomy. In addition, the boundary driven thermal wind flow exhibits a large one cell meridional circulation with surface flow going southward, carrying weak magnetic flux from north to south. These dynamical features are also in agreement with those obtained with a synthetic $Y_1^0$ CMB heat flux pattern (Stanley et al., 2008; Amit et al., 2011).

When the impactor is falling on the equatorial plane, some magnetic dichotomy may be expected between eastern and western hemispheres. Indeed Fig. 6 shows that in case I5 such a dichotomy is obtained. Note that the east-west dichotomy in this case is significantly weaker than the south-north dichotomy in the corresponding polar impactor case I2 with the same impactor radius of 800 km. However, compared to the east-west dichotomy obtained with a synthetic $Y_1^1$ CMB heat flux pattern (Amit et al., 2011), the impactor driven east-west dichotomy is significantly stronger (compare cases I5 and Y4 in Table 1).

Fig. 7 shows two snapshots of the radial field at the CMB and the intensity at the surface of Mars for case I3 with the largest impactor studied here (radius of 1000 km) falling on the north geographic pole as well as the corresponding long term time-average surface intensity. The differences between the two snapshots attest to the chaotic time dependence of this model. Nevertheless, in both snapshots the relatively small-scale radial field at the CMB is strongest at high latitudes of the southern hemisphere. The surface
<figure><img src="image_214.png" /><figcaption>**Fig. 4.** Time-average magnetic field properties on the CMB (left) and on the surface of Mars (right) in case I2. The CMB field is upward continued to the surface of Mars as a potential field.</figcaption></figure>

<figure><img src="image_215.png" /><figcaption>**Fig. 5.** Time-average zonal temperature (left) and flow (right) in case I2. In the right subplot colors denote azimuthal flow and streamlines denote meridional circulation (solid/dashed are anti-clockwise/clockwise, respectively). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

intensity is large scale, and also peaks at the polar region of the southern hemisphere. In this case, the time-average south-north dichotomy is 2.25 (see Table 1), very close to the lower bound estimation of 2.4 obtained by Amit et al. (2011) based on observations of the martian crustal magnetic field (Langlais et al., 2004).

Finally we examine the influence of the strength of the internal core convection. Case I4 is identical to case I3 except for its Ra value which is larger. The south-north magnetic dichotomy in this case is 2.58 (Table 1), demonstrating that more vigorous core convection produces more hemispherical fields.

Our models show that a polar impactor leads to a stronger north-south hemispheric magnetic dichotomy than an equatorial

impactor to an east-west dichotomy. We also find that in the non-reversing regime of parameters, the magnetic field dichotomy that can be recorded in the cooling IIMM is independent of the crust formation scenario (random or continuous). Reversing dynamos that convect stronger produce stronger magnetic dichotomies than non-reversing dynamos. These results are in agreement with those obtained by Amit et al. (2011) for degree-1 heat flux patterns.

The amplitude of the CMB heat flux heterogeneity (Fig. 2) is commonly measured by the peak to peak lateral variation (Olson and Christensen, 2002). However, when the pattern is not smooth and very localized, with imbalance between areas of positive and negative anomalies, $q^*$ is inadequate. This is exactly the case in the impact driven heat flux patterns considered here. For example, for very small impacts the amplitude is negligible but $q^*$ approaches 0.5 (Fig. 1b). We quantify the impact driven heterogeneity amplitude by $q_0^r$ (Eq. (1)) which measures the extent of the CMB surface affected by the impact heating. Therefore, in order to compare the efficiency of magnetic hemispheric dichotomy generation by impact driven patterns with the efficiency by synthetic degree-1 patterns, $q_0^r$ in the first must be compared with $q^*$ in the latter. In addition we propose that the governing parameter controlling the amplitude of the hemispheric magnetic dichotomy is the horizontal Rayleigh number $Ra_h$ (Willis et al., 2007). We define $Ra_h$ as:

$$Ra_h = \begin{cases} q^*Ra & \text{in Y cases} \\ q_0^rRa & \text{in I cases} \end{cases} \qquad (11)$$

To adequately compare the synthetic $Y_1^0$ CMB heat flux cases from Amit et al. (2011) and the polar impact driven CMB heat flux cases from this study, we plot in Fig. 8 the increase of the south-north dichotomy $SN-1$ as a function of $Ra_h$ (11). This figure shows that for a given lateral forcing, an impact induced CMB heat flux heterogeneity is significantly more efficient than a synthetic $Y_1^0$ CMB heat

<figure><img src="image_216.png" /><figcaption>**Fig. 6.** As in Fig. 4 for case I5 (note different color scales). (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_217.png" /><figcaption>**Fig. 7.** Two arbitrary snapshots of the radial magnetic field $B_r$ on the CMB (a and b) and surface intensity of the magnetic field (c and d) in case I3. The time-average surface intensity is shown in (e). Here the spherical harmonic expansion is truncated at $\ell_{max} = 10$. Note differences in color scales among the maps. The CMB field is upward continued to the surface of Mars as a potential field. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

flux heterogeneity in generating hemispheric magnetic dichotomies. Note that the curve of the impact cases is elevated with respect to that of the $Y_1^0$ cases, although the slopes of both curves seem comparable. Fig. 8 also shows that the impactor size needed to generate the observed martian dichotomy (Amit et al., 2011) is about 1000 km, which is within the estimated range of impactor sizes (Marinova et al., 2008). The better efficiency of impact driven hemispheric dynamos over the synthetic cases is even more pronounced in the equatorial cases (see Table 1).

4. Implication for crustal magnetic dichotomy

Our results show that large $Ra_h$ (and therefore large $Ra$) are needed in order to get a dichotomy as strong as the observed. As large $Ra$ are likely to be associated with reversing dynamos (e.g., Kutzner and Christensen, 2002), the possibility to record a dichotomy in the crust depends on the duration of a chron relative to the duration of the crystallization of the IIMM induced by the impact. If a slow crystallization rate was coupled with short magnetic chrons (Dietrich and Wicht, 2013), then our models cannot explain the observed martian dichotomy. However, if a relatively rapid magma cooling was coupled with a relatively low reversal frequency then the impact induced magnetic dichotomy may have been recorded within the martian crust and explain the observed signal. Hence, the relevant magnetic timescale is a duration of a typical chron.

In our non-reversing dynamo models the chron duration is effectively infinite and the hemispheric magnetic dichotomies may be recorded by the crustal magnetization (Fig. 3, left). However, the amplitudes of the dichotomies in these cases are too low to explain the observed hemispheric crustal magnetic dichotomy on Mars (Table 1). Dynamo models with large amplitude

hemispheric dichotomies tend to reverse frequently, in agreement with the findings of Dietrich and Wicht (2013). For example, in case I3 a typical chron persists for ~18 kyr. In order for such a model to record an hemispheric magnetic dichotomy at the crust, very fast crust formation is required.

Dietrich et al. (2013) studied hemispherical dynamos in the framework of classical mean field theory. For $q^* > 0.6$ Dietrich and Wicht (2013) found hemispheric $\alpha\Omega$ dynamos with fast oscillations over periods of ~10 kyrs. In our models the amplitude of CMB heat flux heterogeneity is moderate with $q^* \sim 0.6$, and more importantly $q_0^r \ll 0.6$. Our dynamos are therefore of the $\alpha^2$-type with chaotic (reversing or non-reversing) behavior. In addition, the moderate $q_0^r$ does not change the dynamo regime from stable to reversing. In these dynamos the duration of a chron (or the reversal frequency) depends on the level of turbulence in the core (Olson and Amit, 2014), which is in general unknown (even for the Earth). Chron duration varied immensely over Earth's history between 40 kyrs and 40 Myrs (Merrill et al., 1998).

Fig. 9 illustrates the interplay between the relevant time-scales: the crystallization time $\tau_c$, the spreading time $\tau_{\text{imp}}$ and the duration of a chron. In the context of a martian giant impact, the associated molten thickness $\delta$ ranges between 30 and 50 km and the melt fraction is ~20% (Marinova et al., 2008). Hence, the crystallization and cooling times of the shallow IIMM induced by one impact should be rapid: even if the complete solidification timescale ultimately depends on the poorly constrained evolution of the post-impact transient atmosphere (Abe, 1997), crystallization of this local molten material should occur in less than 1000 years (Reese and Solomatov, 2006). The subsequent evolution of the partially molten material could involve isostatic readjustment of a deep, initially hemispheric, retained melt region and lateral spreading as a gravity current (Reese et al., 2011). Overall, the surface
<figure><img src="image_218.png" /></figure>

**Fig. 8.** South-north dichotomy as a function of the horizontal Rayleigh number $Ra_h$ (Eq. (11)). The black symbols are the values obtained (with $r_i/r_o = 0.2$) for the synthetic $Y_1^0$ CMB heat flux patterns (with $Ra_h = q^*_0 Ra$). The red symbols are the values obtained in cases of polar impact driven CMB heat flux patterns (with $Ra_h = q_0^*_0 Ra$). The size of the red symbols increases with the size of the impactor. The corresponding power law fits are plotted with dashed lines. The green horizontal line represents the martian value of $SN - 1$ from Amit et al. (2011) based on observations of the martian crustal magnetic field (Langlais et al., 2004). All degree-1 cases are from Amit et al. (2011) except cases Y2 and Y3 from this study (see also Table 1). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

temperature of the IIMM falls below the Curie temperature in a timescale that strongly depends on factors such as the characteristics of the above atmosphere or the viscosity of the solid mantle if isostatic rebound is involved. This instant characterizes the beginning of the time interval when the solidified crust starts to record the magnetic field. The cold front then propagates downwards from the surface and a secondary front most probably develops at the base of the IIMM whose nature depends on the thermal state of the pre-impact crust on top of which the IIMM is superimposed (Fig. 9). The crustal material originating in the impact event stops recording the magnetic field precisely when the innermost region is cooler than the Curie temperature. In summary, given these somewhat overlapping ranges of timescales for the crust formation and magnetic chrons, we argue that both end-member crust formation scenarios are in principle possible, bearing in mind that reality may be somewhere in between.

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

**Fig. 9.** Schematic illustration of temporal evolution of core dynamo activity (top panel) and of the IIMM temperature (bottom panel) after a giant impact. The C.A.I. (Calcium-Aluminium rich Inclusions) are the oldest objects in the solar system. In the bottom panel, the blue line represents the temperature evolution at the interior of the IIMM while the green line represents the top of the IIMM (i.e., the post impact martian surface). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Our study also emphasizes the importance of the duration of the post-impact thermal anomaly at the CMB ($\tau_{imp}$). If this timescale is shorter than the time needed by the IIMM to cool down to temperatures smaller than the Curie temperature $\tau_c$, the hemispheric dynamo has ended before the crust was able to record it. However, if $\tau_{imp}$ is larger than $\tau_c$, the post-impact magnetic field can be recorded within the crust. The beginning of this record starts at $\tau_{c,top} \sim 1$ kyr and ends at $\tau_{c,int} \sim 10–100$ Myr, this time interval being mostly controlled by heat diffusion. The deep post-impact thermal re-adjustment occurs within a characteristic time that is governed by the rheology of the mantle surrounding the post-impact thermal anomaly. For a mantle viscosity comparable to that of the present day Earth, say around $10^{21}$ Pa.s, this characteristic spreading time is $\tau_{imp} \sim 10–100$ Myr (Monteux et al., 2007; Watters et al., 2009) which is much larger than the characteristic magnetic diffusion time $\tau_\lambda \sim 30$ kyrs (e.g., Bloxham and Jackson, 1991) or the time needed by the top of the IIMM to reach the Curie temperature $\tau_{c,top}$ (Fig. 9).

## 5. Conclusion

Can the observed martian hemispheric magnetic dichotomy be the consequence of a giant impact that has led to an internal hemispheric magnetic field? Our results show that a $\sim 1000$ km radius impactor can generate a sufficiently large scale thermal anomaly at the CMB so that a hemispheric dichotomy comparable to the observed is generated. An impact induced CMB heat flux heterogeneity is more efficient than a synthetic degree-1 CMB heat flux heterogeneity in generating strong hemispheric magnetic dichotomies. This magnetic dichotomy is stronger for a polar impact than for an equatorial impact. This result reconciles the giant impact induced scenario evoked to explain the martian topographic dichotomy (Andrews-Hanna et al., 2008; Marinova et al., 2008; Nimmo et al., 2008) and the heterogeneous CMB heat flux scenario proposed to explain why the magnetization is strong only in the southern hemisphere (Stanley et al., 2008) via a more realistic CMB heat flux pattern than the previously used synthetic degree-1. Hence, we propose here that an external event at the martian surface may have produced an internal hemispheric dynamo.

From our results, three relevant timescales arise: the cooling time of the impact induced molten material, the duration of a magnetic chron and the duration of the post-impact thermal anomaly at the CMB. The first is mainly governed by the characteristics of the martian atmosphere, the second is governed by the core dynamics while the third is governed by the mantle dynamics. If a relatively rapid magma cooling was coupled with a relatively low reversal frequency and with a stable CMB heat flux heterogeneity, then the impact induced magnetic dichotomy may have been recorded within the martian crust and explain the observed signal.

According to some studies (e.g., Arkani-Hamed and Olson, 2010), a giant impact might have led to a thermal stratification at the top of the core that terminated the dynamo. However, core stratification can be removed much faster by convection of molten material, over a timescale orders of magnitude shorter than the longevity of mantle thermal anomalies. Hence our scenario of giant impact leading to heterogenous mantle heating and to a hemispheric dynamo on Mars may follow an episode of dynamo shutdown.

Our impact driven CMB heat flux heterogeneity model may also be applied to model dynamos of other planets. Indeed, giant impacts were common in the later stage of accretion of terrestrial planets. The Earth is likely formed by accretion of a few dozen moon to Mars-size planetary embryos (see 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 Ward, 1976; Canup, 2004). The Moon was also probably hit by a large planetesimal at the end of its formation (Jutzi and Asphaug, 2011). Finally, an oblique collision of a large body with a mass about one sixth of Mercury's has likely stripped away a significant part of its mantle (Smith, 1979; Benz et al., 1988). Furthermore, these three terrestrial objects have or have had an internally generated magnetic field (Stevenson, 2003). Hence, giant impacts have potentially strongly influenced their internal dynamics and dynamo activities.

## Acknowledgments

J. Monteux was funded by Agence Nationale de la Recherche (Accretis Decision No. ANR-10-PDOC-001-01). We thank Johannes Wicht for providing the MAGIC code. We thank two anonymous reviewers for their careful and insightful reviews that greatly improve the manuscript. The authors also thank N. Mangold and J.P. Lorand for useful discussions.

## References

* Abe, Y., 1997. Thermal and chemical evolution of the terrestrial magma ocean. *Phys. Earth Planet. Inter.* 100, 27–39.

* Acuña, M., Connerney, J., Wasilewski, P., Lin, R., Anderson, K., Carlson, C., McFadden, J., Curtis, D., Mitchell, D., Reme, H., Mazelle, C., Sauvaud, J., d'Uston, C., Cros, A., Medale, J., Bauer, S., Cloutier, P., Mayhew, M., Winterhalter, D., Ness, N., 1998. Magnetic field and plasma observations at Mars: initial results of the Mars global surveyor mission. *Science* 279 (5357), 1676–1680.

* Amit, H., Choblet, G., 2009. Mantle-driven geodynamo features – effects of post-Perovskite phase transition. *Earth Planets Space* 61, 1255–1268.

* Amit, H., Choblet, G., 2012. Mantle-driven geodynamo features – effects of compositional and narrow D" anomalies. *Phys. Earth Planet. Inter.* 190, 34–43.

* Amit, H., Aubert, J., Hulot, G., 2010. Stationary, oscillating or drifting mantle-driven geomagnetic flux patches? *J. Geophys. Res.* 115, B07108.

* Amit, H., Christensen, U., Langlais, B., 2011. The influence of degree-1 mantle heterogeneity on the past dynamo of Mars. *Phys. Earth Planet. Inter.* 189, 63–79.

* Anderson, B.J., Johnson, C.L., Korth, H., Winslow, R.M., Borovsky, J.E., Purucker, M.E., Slavin, J.A., Solomon, S.C., Zuber, M.T., McNutt Jr., R.L., 2012. Low-degree structure in Mercury's planetary magnetic field. *J. Geophys. Res. (Planets)* 117.

* 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.

* Arkani-Hamed, J., Olson, P., 2010. Giant impact stratification of the Martian core. *Geophys. Res. Lett.* 37, L02201.

* Asphaug, E., Agnor, C.B., Williams, Q., 2006. Hit-and-run planetary collisions. *Nature* 439, 155–160.

* Aubert, J., Amit, H., Hulot, G., 2007. Detecting thermal boundary control in surface flows from numerical dynamos. *Phys. Earth Planet. Inter.* 160, 143–156.

* Aubert, J., Amit, H., Hulot, G., Olson, P., 2008. Thermo-chemical wind flows couple Earth's inner core growth to mantle heterogeneity. *Nature* 454, 758–761.

* Aubert, J., Labrosse, S., Poitou, C., 2009. Modelling the paleo-evolution of the geodynamo. *Geophys. J. Int.* 179, 1414–1428.

* Benz, W., Slattery, W.L., Cameron, A.G.W., 1988. Collisional stripping of Mercury's mantle. *Icarus* 74, 516–528.

* Bloxham, J., 2002. Time-independent and time-dependent behaviour of high-latitude flux bundles at the core-mantle boundary. *Geophys. Res. Lett.* 29, 1854.

* Bloxham, J., Jackson, A., 1991. Fluid flow near the surface of the Earth's outer core. *Rev. Geophys.* 29, 97–120.

* Breuer, D., Labrosse, S., Spohn, T., 2010. Thermal evolution and magnetic field generation in terrestrial planets and satellites. *Space Sci. Rev.* 152, 449–500.

* Cameron, A.G.W., Ward, W.R., 1976. The origin of the Moon. *Lun. Planet. Inst. Conf. Abs.* 7, 120–122.

* Canup, R.M., 2004. Simulations of a late lunar-forming impact. *Icarus* 168, 433–456.

* Cao, H., Aurnou, J.M., Wicht, J., Dietrich, W., Soderlund, K.M., Russell, C.T., 2014. A dynamo explanation for Mercury's anomalous magnetic field. *Geophys. Res. Lett.* 41, 4127–4134.

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

* Chassefière, E., Langlais, B., Quesnel, Y., Leblanc, F., 2013. The fate of early Mars' lost water: the role of serpentinization. *J. Geophys. Res. (Planets)* 118, 1123–1134.

* Dietrich, W., Wicht, J., 2013. A hemispherical dynamo model: implications for the Martian crustal magnetization. *Phys. Earth Planet. Int.* 217, 10–21.

* Dietrich, W., Schmitt, D., Wicht, J., 2013. Hemispherical Parker waves driven by thermal shear in planetary dynamos. *Europhys. Lett.*, 104

* 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.

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

* Frey, H.V., Roark, J.H., Shockey, K.M., Frey, E.L., Sakimoto, S.E.H., 2002. Ancient lowlands on Mars. *Geophys. Res. Lett.* 29, 1384.

* Gladman, B., Coffey, J., 2009. Mercurian impact ejecta: Meteorites and mantle. *Meteorit. Planet. Sci.* 44, 285–291.

* Glatzmaier, G., Coe, R., Hongre, L., Roberts, P., 1999. The role of the Earth's mantle in controlling the frequency of geomagnetic reversals. *Nature* 401, 885–890.

* 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.

* Gubbins, D., Willis, P.W., Sreenivasan, B., 2007. Correlation of Earth's magnetic field with lower mantle thermal and seismic structure. *Phys. Earth Planet. Inter.* 162, 256–260.

* Harder, H., Christensen, U., 1996. A one-plume model of Martian mantle convection. *Nature* 380, 507–509.

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

* Hori, K., Wicht, J., Christensen, U., 2010. The effect of thermal boundary conditions on dynamos driven by internal heating. *Phys. Earth Planet. Inter.* 182, 85–97.

* Hori, K., Wicht, J., Dietrich, W., 2014. Ancient dynamos of terrestrial planets more sensitive to core-mantle boundary heat flows. *Planet. Space. Sci.* 98, 30–40.

* Jutzi, M., Asphaug, E., 2011. Forming the lunar farside highlands by accretion of a companion moon. *Nature* 476, 69–72.

* Ke, Y., Solomatov, V.S., 2006. Early transient superplumes and the origin of the Martian crustal dichotomy. *J. Geophys. Res. (Planets)* 111, 10001.

* Kutzner, C., Christensen, U., 2002. From stable dipolar towards reversing numerical dynamos. *Phys. Earth Planet. Inter.* 131, 29–45.

* Kutzner, C., Christensen, U.R., 2004. Simulated geomagnetic reversals and preferred virtual geomagnetic pole paths. *Geophys. J. Int.* 157, 1105–1118.

* Langlais, B., Amit, H., 2008. The past martian dynamo. *Science* 321, 1784–1785.

* Langlais, B., Thébault, E., 2011. Predicted and observed magnetic signatures of martian (de)magnetized impact craters. *Icarus* 212, 568–578.

* Langlais, B., Purucker, M., Mandea, M., 2004. Crustal magnetic field of Mars. *J. Geophys. Res.* 109, E02008.

* Lebrun, T., Massol, H., Chassefière, E., Davaille, A., Marcq, E., Sarda, P., Leblanc, F., Brandeis, G., 013. Thermal evolution of an early magma ocean in interaction with the atmosphere. *J. Geophys. Res. (Planets)* 118, 1155–1176.

* Le Feuvre, M., Wieczorek, M.A., 2011. Nonuniform cratering of the Moon and a revised crater chronology of the inner Solar System. *Icarus* 214, 1–20.

* 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. *Geophys. Res. Lett.* 35, 14203.

* Lillis, R.J., Stewart, S.T., Manga, M., 2013. Demagnetization by basin-forming impacts on early Mars: contributions from shock, heat, and excavation. *J. Geophys. Res. (Planets)* 118, 1045–1062.

* Marinova, M.M., Aharonson, O., Asphaug, E., 2008. Mega-impact formation of the Mars hemispheric dichotomy. *Nature* 453, 1216–1219.

* Marinova, M.M., Aharonson, O., Asphaug, E., 2011. Geophysical consequences of planetary-scale impacts into a Mars-like planet. *Icarus* 211, 960–985.

* Merrill, R., McElhinny, M., McFadden, P., 1998. The Magnetic Field of the Earth: Paleomagnetism, the Core, and the Deep Mantle. *Academic Press, San Diego, California, USA*.

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

* 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., Jellinek, A.M., Johnson, C.L., 2011. Why might planets and moons have early dynamos? *Earth Planet. Sci. Lett.* 310, 349–359.

* Monteux, J., Jellinek, A.M., Johnson, C.L., 2013. Dynamics of core merging after a mega-impact with applications to Mars' early dynamo. *Icarus* 226, 20–32.

* 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.

* Olson, P., 2007. Gravitational dynamos and the low frequency geomagnetic secular variation. *Proc. Natl. Acad. Sci.* 104, 20159–20166.

* Olson, P., Amit, H., 2014. Magnetic reversal frequency scaling in dynamos with thermochemical convection. *Phys. Earth Planet. Inter.* 229, 122–133.

* Olson, P., Christensen, U., 2002. The time averaged magnetic field in numerical dynamos with nonuniform boundary heat flow. *Geophys. J. Int.* 151, 809–823.

* Olson, P., Coe, R.S., Driscoll, P.E., Glatzmaier, G.A., Roberts, P.H., 2010. Geodynamo reversal frequency and heterogeneous core-mantle boundary heat flow. *Phys. Earth Planet. Inter.* 180, 66–79.

* Olson, P., Deguen, R., Hinnov, L.A., Zhong, S., 2013. Controls on geomagnetic reversals and core evolution by mantle convection in the Phanerozoic. *Phys. Earth Planet. Inter.* 214, 87–103.

* Parker, R., 2003. Ideal bodies for Mars magnetics. *J. Geophys. Res. Planets* 108 (E1).

* Pierazzo, E., Melosh, H.J., 2000. Melt production in oblique impacts. *Icarus* 145, 252–261.

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

* Pozzo, M., Davies, C., Gubbins, D., Alfe, D., 2012. Thermal and electrical conductivity of iron at Earth's core conditions. *Nature* 485, 355–358.
Quesnel, Y., Sotin, C., Langlais, B., Costin, S., Mandea, M., Gottschalk, M., Dyment, J., 2009. Serpentinization of the martian crust during Noachian. Earth Planet. Sci. Lett. 277, 184–193.

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.

Reese, C.C., Orth, C.P., Solomatov, V.S., 2010. Impact origin for the Martian crustal dichotomy: Half emptied or half filled? J. Geophys. Res. (Planets) 115 (14), 5004–+.

Reese, C.C., Orth, C.P., Solomatov, V.S., 2011. Impact megadomes and the origin of the martian crustal dichotomy. Icarus 213, 433–442.

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

Roberts, J., Zhong, S., 2006. Degree-1 convection in the Martian mantle and the origin of the hemispheric dichotomy. J. Geophys. Res. 111, E06013.

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.

Rochette, P., 2006. Crustal magnetization of Mars controlled by lithology or cooling rate in a reversing dynamo? Geophys. Res. Lett. 33, 2202.

Schubert, G., Spohn, T., 1990. Thermal history of Mars and the sulfur content of its core. J. Geophys. Res. 95, 14095–14104.

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

Shuvalov, V., 2009. Atmospheric erosion induced by oblique impacts. Meteorit. Planet. Sci. 44, 1095–1105.

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

Solomon, S.C., Aharonson, O., Aurnou, J.M., Banerdt, W.B., Carr, M.H., Dombard, A.J., Frey, H.V., Golombek, M.P., Hauck, S.A., Head, J.W., Jakosky, B.M., Johnson, C.L.,

McGovern, Neumann, G.A., Phillips, R.J., Smith, D.E., Zuber, M.T., 2005. New Perspectives on Ancient Mars. Science 307, 1214–1220.

Šrámek, O., Ricard, Y., Dubuffet, F., 2010. A multiphase model of core formation. Geophys. J. Int. 181, 198–220.

Stanley, S., 2010. A dynamo model for axisymmetrizing Saturn's magnetic field. Geophys. Res. Lett. 37, 5201.

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.

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

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.

Wicht, J., 2002. Inner-core conductivity in numerical dynamo simulations. Phys. Earth Planet. Inter. 132, 281–302.

Wicht, J., Stellmach, S., Harder, H., 2009. Numerical models of the geodynamo: from fundamental Cartesian models to 3D simulations of field reversals. In: Glassmeier, H., Soffel, H., Negendank, J. (Eds.), Geomagnetic Field Variations – Space-time Structure, Processes, and Effects on System Earth. Springer, Berlin, pp. 107–158.

Willis, A.P., Sreenivasan, B., Gubbins, D., 2007. Thermal core mantle interaction: exploring regimes for locked dynamo action. Phys. Earth Planet. Inter. 165, 83–92.

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.

Zhong, S., Zuber, M.T., 2001. Degree-1 mantle convection and the crustal dichotomy on Mars. Earth Planet. Sci. Lett. 189, 75–84.

Zuber, M.T., 2001. The crust and mantle of Mars. Nature 412, 220–227.
<figure><img src="image_220.png" /></figure>

REVIEW

Open Access

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

CrossMark

# Towards more realistic core-mantle boundary heat flux patterns: a source of diversity in planetary dynamos

Hagay Amit<fnref n="1" />*, Gaël Choblet<fnref n="1" />, Peter Olson², Julien Monteux³, Frédéric Deschamps⁴, Benoit Langlais<fnref n="1" /> and Gabriel Tobie<fnref n="1" />

**Abstract**

Mantle control on planetary dynamos is often studied by imposing heterogeneous core-mantle boundary (CMB) heat flux patterns on the outer boundary of numerical dynamo simulations. These patterns typically enter two main categories: Either they are proportional to seismic tomography models of Earth's lowermost mantle to simulate realistic conditions, or they are represented by single spherical harmonics for fundamental physical understanding. However, in reality the dynamics in the lower mantle is much more complicated and these CMB heat flux models are most likely oversimplified. Here we term *alternative* any CMB heat flux pattern imposed on numerical dynamos that does not fall into these two categories, and instead attempts to account for additional complexity in the lower mantle. We review papers that attempted to explain various dynamo-related observations by imposing alternative CMB heat flux patterns on their dynamo models. For present-day Earth, the alternative patterns reflect non-thermal contributions to seismic anomalies or sharp features not resolved by global tomography models. Time-dependent mantle convection is invoked for capturing past conditions on Earth's CMB. For Mars, alternative patterns account for localized heating by a giant impact or a mantle plume. Recovered geodynamo-related observations include persistent morphological features of present-day core convection and the geomagnetic field as well as the variability in the geomagnetic reversal frequency over the past several hundred Myr. On Mars the models aim at explaining the demise of the paleodynamo or the hemispheric crustal magnetic dichotomy. We report the main results of these studies, discuss their geophysical implications, and speculate on some future prospects.

**Keywords:** Magnetic field; Dynamo; Core-mantle boundary; Heat flux

**Introduction**

**General**

Various geophysical observations have been recovered by imposing heterogeneous heat flux patterns on the outer boundary of numerical dynamo simulations. For the Earth, the most popular prescribed heat flux patterns are either proportional to seismic velocity anomalies obtained from tomography models of the lowermost mantle to mimic realistic conditions, or spherical harmonic degree and order 2 which is the dominant term in these tomography models (Aubert et al. 2007, 2008; Bloxham 2002; Christensen and Olson 2003; Davies et al. 2008;

Glatzmaier et al. 1999; Gubbins et al. 2007, 2011; Heimpel and Evans 2013; Kutzner and Christensen 2004; Olson and Amit 2014; Olson and Christensen 2002; Olson et al. 2010; Sreenivasan 2009; Sreenivasan and Gubbins 2011; Takahashi et al. 2008). For other planets, various single harmonics are often used, for example degree-1 for Mars (Amit et al. 2011a; Dietrich and Wicht 2013; Stanley et al. 2008) and degree-2 for Mercury and Saturn (Cao et al. 2014; Stanley 2010, respectively), reflecting the low-resolution information on planetary deep interiors. Various single harmonics were imposed for the geodynamo as well in order to obtain fundamental physical understanding of the complex dynamo simulations (e.g., Glatzmaier et al. 1999).

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

<footnote>*Correspondence: Hagay.Amit@univ-nantes.fr</footnote>

<footnote n="1">¹CNRS UMR 6112, Université de Nantes, Laboratoire de Planétologie et de Géodynamique, 2 rue de la Houssinière, Nantes, F-44000, France
Full list of author information is available at the end of the article</footnote>
Inferring core-mantle boundary (CMB) heat flux from tomography models relies on the assumption that seismic and thermal anomalies in the lowermost mantle are highly correlated. However, this idealized linear relation may not hold in reality. Observational and numerical studies of the lower mantle suggest that the structures and dynamics there are far more complex. First, non-thermal contributions (e.g., compositional, phase changes) perturb the idealized linear seismic-thermal mapping (Bull et al. 2009; Deschamps et al. 2012; Lay et al. 2008; Nakagawa and Tackley 2008; Tackley 2011; Trampert et al. 2004). Second, while the spatial resolution inherent to tomographic models in the lowermost mantle (e.g., Lekic et al. 2012) exceeds the resolution of the observed core field, strong lateral gradients of shear wave velocities are not captured. As revealed by detailed waveform analyses at several locations, boundaries of large low shear-wave velocity provinces (LLSVPs) might however be very sharp (e.g., Ni et al. 2005; To et al. 2005). Similarly, synthetic single harmonic patterns for planetary mantle heterogeneity, for example for the CMB conditions that prevailed during the paleo dynamo of Mars, are oversimplified. Even if the Martian anomaly was indeed large-scale (Elkins-Tanton et al. 2005; Harder and Christensen 1996; Ke and Solomatov 2006; Roberts and Zhong 2006; Roberts et al. 2009), it was likely more complex than a single harmonic pattern. Note finally that simulating the geodynamo at earlier times requires information about the time-dependent mantle convection pattern, which is not witnessed by present-day tomography (Yoshida and Santosh 2011; Zhang and Zhong 2011; Zhang et al. 2010; Zhong et al. 2007).

Here we term any CMB heat flux pattern that is neither tomographic nor single harmonic as *alternative*. For the Earth, the alternative patterns include those that account for non-thermal effects or non-resolved small scales in various ways (Amit and Choblet 2009, 2012; Amit et al. in press) and those that model the mantle heterogeneity in the past hundreds of Myr in various ways (Amit and Olson 2015; Olson et al. 2013). For Mars, the alternative patterns are obtained by modeling localized mantle heating induced by a mantle plume (Sreenivasan and Jellinek 2012) or a giant impact (Kuang et al. 2014; Monteux et al. 2015).

Figure 1a illustrates some proposed dynamical scenarios of Earth's lower mantle that result in alternative CMB heat flux patterns (see also, e.g., Lay et al. 2008). Subducting slabs bring cold material to the CMB and locally increase the heat flux. Time-dependent plate motions therefore change the location of these large CMB heat flux regions. Further complications arise from the presence of post-perovskite phase transition (or even its double crossing) in the slab reservoirs of the lowermost mantle. LLSVPs may represent thermochemical piles of a not only hot

but also compositionally different material (e.g., Trampert et al. 2004). Growth and collapse of these thermochemical piles, as observed in analog experiments (Davaille 1999), may yield temporal changes in the total heat loss through the CMB and the amplitude of the lateral heterogeneity. Ultra low velocity zones (ULVZs) at the edges of LLSVPs (McNamara et al. 2010) form hot regions that are too narrow to be captured by global large-scale tomography models.

Figure 1b illustrates the effects of localized mantle heating on Mars' CMB heat flux. Giant impacts heat the mantle and reduce the CMB heat flux right below the impact site. Likewise, mantle plumes that originate at the lowermost mantle and may have produced extensive volcanism such as Tharsis are also potential zones of reduced CMB heat flux. It may also be naively expected that core heating by impacts would yield a CMB heat flux increase. However, core impact heating leads to the formation of a hot layer at the top of the core (within 10 kyr) which prevents the core from cooling (Arkani-Hamed 2012; Arkani-Hamed and Olson 2010; Roberts and Arkani-Hamed 2014). Thus, core impact heating could inhibit the dynamo generation during a timescale that is governed by the efficiency of the surrounding mantle to extract the impact induced core thermal anomaly (from several kyr up to 100 Myr).

Numerical simulations of mantle convection may shed light on the relation between the observed seismic anomalies and the CMB heat flux heterogeneity. Nakagawa and Tackley (2008) assumed a linear equation of state to relate their models' thermal, mineral, and chemical anomalies with seismic anomalies. Their purely thermal case gives nearly perfect linear relation between seismic and thermal anomalies, whereas adding post-perovskite effects breaks the linearity but maintains a unique relation (Fig. 2a). The situation is more complicated when chemical anomalies are accounted for; when the density contrast is large enough, a parallel branch appears and the thermal-seismic relation becomes non-unique (Fig. 2b). Implementation of the results of Nakagawa and Tackley (2008) might be problematic because due to limited vertical resolution global tomographic models do not sample precisely the CMB but rather a vertically averaged region above it.

The outline of this review paper is as follows. In the “Alternative models of core-mantle boundary heat flux” subsection, we introduce motivations for alternative CMB heat flux models. Related dynamo properties inferred from observations are listed in subsection “Relevant observed planetary dynamo properties”. In the section “Methods”, we recall the governing equations and control parameters of numerical dynamo simulations. We describe the incorporation of mantle control effects by imposing an outer thermal boundary condition. In the section “Review”, we describe the results of papers that
<figure><img src="image_223.png" /></figure>

imposed alternative outer boundary conditions on numerical dynamos. For each of these papers, we specify the procedure used to construct the alternative CMB heat flux model and we report how the perturbation from the conventional reference pattern affected the resulting planetary dynamo features. We discuss the results in the subsection “Discussion”, and we propose some future prospects in subsection “Future prospects”.

**Alternative models of core-mantle boundary heat flux**

The existence of the post-perovskite phase transition (Murakami et al. 2004; Oganov and Ono 2004) was invoked to explain the seismic heterogeneity of the

lowermost mantle (Lay et al. 2006). Post-perovskite is needed to explain core-reflected (Hernlund et al. 2005) and diffracted (Cobden et al. 2012) waves. It may also be a good candidate for the D" discontinuity, at least in regions where P- and S-waves have opposite signs (Cobden and Thomas 2013). In contrast, in probabilistic tomography models the thermal distribution is similar with and without post-perovskite (Mosca et al. 2012). Mantle convection simulations that incorporate post-perovskite effects indicate that its presence affects mantle dynamics (Cizkova et al. 2010; Nakagawa and Tackley 2011) and distorts the thermal-seismic relation from its linearity (Nakagawa and Tackley 2008).
<figure><img src="image_224.png" /></figure>

<figure><img src="image_225.png" /><figcaption>**Fig. 2** Scatter plots of CMB heat flux anomalies vs. seismic shear velocity anomalies at the lowermost mantle based on mantle convection simulations. **a** Isochemical cases without (red) and with (green) post-perovskite. **b** Thermochemical cases with various density contrasts $\delta\rho_m$ at the CMB. From Nakagawa and Tackley (2008)</figcaption></figure>

Sharp margins of LLSVPs were suggested by local studies of the lowermost mantle (Ni et al. 2002; To et al. 2005; Wang and Wen 2007). Dynamical models indicate that ULVZs likely correspond to the hottest regions of LLSVPs (McNamara et al. 2010). Thin hot narrow regions on the edges of the LLSVPs were identified in two end-member models of lower mantle dynamics (Bull et al. 2009). The low spatial resolution of global tomography models excludes such short wavelength features (e.g., Ritsema et al. 2007).

LLSVPs are better explained if a compositional component is also present (Deschamps et al. 2012; Mosca et al. 2012; Trampert et al. 2004). These structures could be primitive or may consist of recycled crustal material, although the seismic sensitivity of mantle material indicates that the latter is less likely (Deschamps et al. 2012). By taking into account seismic normal modes, Ishii and Tromp (1999) provided an independent constraint on mantle density distributions, which can be used to resolve thermal and compositional contributions to seismic anomalies, as done by probabilistic tomography studies (Mosca et al. 2012; Trampert et al. 2004).

One obvious motivation to consider alternative CMB heat flux patterns is the time-dependence of mantle dynamics. Global plate motion models indicate that the distribution of continents has evolved from a spherical harmonic degree-1 heterogeneity before the breakup of supercontinent Pangea to its present-day degree-2 structure (Seton et al. 2012). Some mantle convection models constrained by reconstructed time-dependent plate motions suggest that the supercontinent formation and breakup was accompanied by substantial change in mantle convection pattern (Zhang and Zhong 2011; Zhang et al. 2010; Zhong et al. 2007). According to these dynamical models, the low seismic velocity structures below

the Pacific and Africa predate and postdate, respectively, the Pangea breakup (McNamara and Zhong 2005; Zhang et al. 2010). These time-dependent mantle convection patterns would inevitably correspond to time-dependent CMB heat flux with patterns distinctively different from that inferred from present-day tomography models.

In contrast, the pattern of LLSVPs has been interpreted as a root mantle structure that prevailed for several hundreds of Myr (Burke 2011; Burke et al. 2008; Dziewonski et al. 2010). Evidence for long-lived mantle structure comes from hotspot reconstructions (Torsvik et al. 2006) and paleomagnetic pole locations that indicate substantial true polar wander in the Mesozoic (Courtillot and Besse 1987). The dynamical origin of this (stationary) degree-2 pattern may be two mantle superplumes (Romanowicz and Gung 2002; Torsvik et al. 2010), plume clusters (Schubert et al. 2004) or dense chemical piles above the CMB (McNamara and Zhong 2005; Tackley 2002; Tan and Gurnis 2007). Permanent or not, the thermochemical piles are clearly dynamic and respond to changes in the global mantle circulation, through variations in their footprints above the CMB and their height in the lower mantle. Variations in either or both of these would affect the mean CMB heat flux and the amplitude of its heterogeneity.

Alternative CMB heat flux models were also constructed for planets other than the Earth. It has been argued that the emergence of the Tharsis volcanic complex coincided with the termination of the Martian dynamo (Jellinek et al. 2008; Johnson and Phillips 2005; Lillis et al. 2008). Mars' mantle may contain a single huge plume (Harder and Christensen 1996; Srámek and Zhong 2010). If Tharsis emerged from such a deep mantle plume, then it may indeed have affected the thermal heterogeneity on Mars' CMB.
Another possible scenario is localized mantle heating by giant impacts. This process enhances the iron/silicate separation and modifies the heat repartition between the mantle and the core. Thus, large impacts potentially affected the efficiency of core cooling and dynamo activity during (Monteux et al. 2011, 2012) or after terrestrial core formation (Glenn Sterenborg and Crowley 2013; Ke and Solomatov 2009; Monteux and Arkani-Hamed 2014; Monteux et al. 2013; Reese and Solomatov 2010). This large-scale heating also strongly influences mantle dynamics (Reese et al. 2002; Roberts and Arkani-Hamed 2012; Roberts and Barnouin 2012; Watters et al. 2009) and hence the heat repartition in the deepest part of terrestrial bodies. The onset or demise of terrestrial dynamos were modeled by giant impacts (Arkani-Hamed and Olson 2010; Kuang et al. 2014; Monteux and Arkani-Hamed 2014; Monteux et al. 2013; Roberts et al. 2009). Large impacts heated Mars' mantle below the impact site, with the heated volume depending mostly on the impactor's size and velocity (Pierazzo et al. 1997; Tonks and Melosh 1993). The impact-driven heating may be modeled by a uniform spherical temperature increase with rapid temperature decrease away from the isothermal sphere (Monteux et al. 2007, 2013; Senshu et al. 2002). Hence, for large impacts, the CMB heat flux may be affected.

Giant impacts may also have affected the morphology of the Martian paleomagnetic field. These events were invoked to explain the Martian and lunar topographic hemispheric dichotomies (Jutzi and Asphaug 2011; Marinova et al. 2008; Nimmo et al. 2008). Numerical dynamos with a synthetic $Y_1^0$ heat flux pattern lead to magnetic fields concentrated at the large heat flux hemisphere (Amit et al. 2011a; Dietrich and Wicht 2013; Stanley et al. 2008). Giant impacts may produce a similar effect with a more realistic CMB heat flux pattern (Monteux et al. 2015).

## Relevant observed planetary dynamo properties

Alternative (or any) heterogeneous CMB heat flux is expected to influence the dynamo on very long timescales (Olson and Christensen 2002), but some statistical preference can already be detected on shorter timescales (Amit et al. 2010a; Bloxham 2002). Stronger heat flux heterogeneity amplitudes increase the probability of a high correlation between a dynamo property and its mantle-driven long-term average (Aubert et al. 2007). If the heterogeneity amplitude is very large compared to the mean superadiabatic CMB heat flux, locking may occur and a mantle signature may appear in a snapshot (Willis et al. 2007).

Because neither the dynamo equations nor the geometry of planetary outer cores contain any special longitude, any persistent deviation of a dynamo property from

axisymmetry can only be explained by boundary heterogeneity. In contrast there are special latitudes, obviously the geographic equator and poles, but also the latitude of the tangent cylinder (Aurnou et al. 2003). However, there is no preferred hemisphere (neither in the dynamo equations nor in the geometry of planetary outer cores), so any persistent deviation of a dynamo property from equatorial symmetry may also require boundary heterogeneity. Below we list the main dynamo-related observations that motivated imposing heterogeneous CMB heat flux, in particular alternative patterns.

The primary dynamo observation is the observed magnetic field at the surface, from which a truncated (or large-scale) part can be downward continued to model the radial magnetic field on the CMB. Many studies attempted to reproduce the locations of high-latitude intense geomagnetic flux patches. The existence of two of these prominent geomagnetic field features at each hemisphere (Jackson et al. 2000) is consistent with the $Y_2^2$ dominance of the lowermost mantle seismic heterogeneity (e.g., Masters et al. 2000). The longitudes of the flux patches and the large seismic shear velocity structures seem to be correlated (Gubbins 2003), although some patches exhibit significant mobility (Amit et al. 2011b; Korte and Holme 2010). The time-average radial field of the dynamo model of Aubert et al. (2008) decently reproduces the longitudes of the flux patches in the paleomagnetic field (Kelly and Gubbins 1997). However, a shift between the locations of the large heat flux structures and the patches was observed, with the shift value depending on the dynamo internal parameters (Aubert et al. 2007; Olson and Christensen 2002; Takahashi et al. 2008). Alternatively, some geomagnetic (or paleomagnetic) spectral information, e.g., the number of patches at each hemisphere (corresponding to the dominant order), may be used as an observational constraint (Amit and Choblet 2009).

While high-latitude intense flux patches near the edge of the tangent cylinder (Aurnou et al. 2003) are a robust feature of dipole-dominated numerical dynamos (Christensen et al. 1998), the observed intense geomagnetic flux patches at low-latitude (Jackson 2003) are less characteristic of these models. The morphological criteria for Earth-like dynamos proposed by Christensen et al. (2010) identify well magnetic field models that resemble the geomagnetic field, but these criteria contain no information about the latitudinal distribution of magnetic flux. Alternative CMB heat flux heterogeneity offers a ready explanation for localized low-latitude intense flux patches (Amit and Choblet 2012; Amit et al. in press). However, such flux patches are not necessarily connected to boundary heterogeneity. Mobile low-latitude flux patches are seen in numerical dynamos without boundary heterogeneity when a strong enough influence
of rotation on outer core convection is present (Wicht et al. 2011).

Other geodynamo properties apart from the magnetic field may also be used as observational constraints to dynamo models with heterogeneous CMB heat flux. For example, core flow models inverted from the geomagnetic secular variation (Holme 2007) show active Atlantic vs. quiet Pacific hemispheres (Amit and Olson 2006; Aubert et al. 2013; Gillet et al. 2009; Holme and Olsen 2006; Hulot et al. 2002). In addition, the seismic velocity, anisotropy, and attenuation of the inner core exhibit an east-west hemispheric dichotomy (Tanaka and Hamaguchi 1997), although these three seismic properties may sample different parts of the inner core (Deuss 2014; Irving and Deuss 2011). Aubert et al. (2008) recovered the observed patterns of the paleomagnetic field and the inner core hemispheric dichotomy in the corresponding time-average quantities of their dynamo model. The inner core hemispheric dichotomy may indeed be mantle-driven (Aubert et al. 2008; Gubbins et al. 2011), or alternatively independent inner core dynamics (Alboussiére et al. 2010; Monnereau et al. 2010) may produce anomalous bouyancy flux that would affect the outer core dynamics and the properties of the CMB (Aubert et al. 2013; Olson and Deguen 2012).

On much longer timescales of hundreds of Myr, the enigmatic variability of the geomagnetic reversal frequency may constrain changes in mantle conditions. In particular, the existence of long superchrons of tens of Myr in which no reversal was found in the paleomagnetic record in contrast to periods of hyper-reversing behavior (Gradstein et al. 2012; Merrill et al. 1998) is difficult to reconcile without time-dependent CMB conditions (Amit and Olson 2015; Biggin et al. 2012; Olson et al. 2013). In mantle convection simulations constrained by reconstructed plate motions the timing of the Kiaman Reversed Superchron (KRS) at around 270–315 Ma coincides with the low mean and equatorial CMB heat flux (Zhang and Zhong 2011). However, fast plate motions during the most recent Cretaceous Normal Superchron (CNS) may result in large CMB heat flux and frequent reversals (Olson et al. 2013).

Strong crustal magnetic field measured by Mars Global Surveyor is suggestive of a core dynamo that ceased at some time during the Noachian (Acuña et al. 2001; Langlais and Purucker 2007; Lillis et al. 2008; Milbury et al. 2012). Localized mantle heating may have affected the efficiency of the dynamo generation on Mars. A strong mantle plume rooted at the CMB could have terminated the Martian dynamo (Sreenivasan and Jellinek 2012). Alternatively, if the Martian dynamo was subcritical (Kuang et al. 2008), a moderate CMB heat flux anomaly induced by a giant impact could have stopped the dynamo (Kuang et al. 2014).

It has been proposed that the hemispheric dichotomy in Mars' crustal magnetic field was caused by heterogeneous CMB heat flux when the Martian paleo dynamo operated (Stanley et al. 2008). The magnitude of this dichotomy (Amit et al. 2011a) is recovered only with strongly connecting reversing dynamo models (Dietrich and Wicht 2013). Thus, in order to record hemispheric magnetization, the Martian crust should have cooled over a shorter period than the typical chron time (Lebrun et al. 2013), which poses a constraint on the feasibility of the hemispherical dynamo scenario for explaining the hemispheric dichotomy in Mars' crustal magnetic field (Dietrich and Wicht 2013; Monteux et al. 2015).

## Methods

The following set of self-consistent non-dimensional Boussinesq (i.e., incompressible) magnetohydrodynamics equations for dynamo action due to thermochemical convection of an electrically conducting fluid in a rotating spherical shell were solved (e.g., Olson et al. 1999):

$$ 
\begin{aligned}
E \left( \frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \nabla \vec{u} - \nabla^2 \vec{u} \right) + 2\hat{z} \times \vec{u} + \nabla P \\
= Ra \frac{\vec{r}}{r_o} T + \frac{1}{Pm} (\nabla \times \vec{B}) \times \vec{B}
\end{aligned} 
\quad (1) $$

$$ \frac{\partial \vec{B}}{\partial t} = \nabla \times (\vec{u} \times \vec{B}) + \frac{1}{Pm} \nabla^2 \vec{B} \quad (2) $$

$$ \frac{\partial T}{\partial t} + \vec{u} \cdot \nabla T = \frac{1}{Pr} \nabla^2 T + \epsilon \quad (3) $$

$$ \nabla \cdot \vec{u} = 0 \quad (4) $$

$$ \nabla \cdot \vec{B} = 0 \quad (5) $$

where $\vec{u}$ is the velocity, $\vec{B}$ is the magnetic field, $T$ is co-density (Braginsky and Roberts 1995), $t$ is time, $\hat{z}$ is a unit vector in the direction of the rotation axis, $P$ is modified pressure which contains centrifugal effects, $\vec{r}$ is the position vector, $r_o$ is the core radius, and $\epsilon$ is co-density source or sink. Four internal non-dimensional parameters in (1)–(5) control the dynamo action. The Ekman number represents the ratio of viscous and Coriolis forces

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

where $\nu$ is kinematic viscosity. The Prandtl number is the ratio of kinematic viscosity to thermal diffusivity $\kappa$

$$ Pr = \frac{\nu}{\kappa} \quad (7) $$

and the magnetic Prandtl number is the ratio of kinematic viscosity to magnetic diffusivity $\lambda$

$$ Pm = \frac{\nu}{\lambda} \quad (8) $$

The strength of buoyancy force driving the convection relative to retarding forces is represented by the
Rayleigh number. When the CMB heat flux is strongly superadiabatic, Ra can be defined based on the mean CMB heat flux $q_0$ as (Olson and Christensen 2002)

$$Ra = \frac{\alpha g_0 q_0 D^4}{k \kappa \nu} \qquad (9)$$

where $\alpha$ is thermal expansivity, $g_0$ is gravitational acceleration on the outer boundary at radius $r_o$, D is shell thickness, and k is thermal conductivity. The heat flux based Ra (9) was used with either generic thermochemical convection (Amit and Choblet 2009; Amit et al. in press; Amit and Choblet 2012) or purely thermal convection (Monteux et al. 2015). Alternatively, when convection is dominantly chemical, Ra can be defined based on the rate of light elements release into the outer core due to inner core freezing as (Olson 2007a)

$$Ra = \frac{\beta g_0 D^5 \dot{\chi}}{\kappa \nu^2} \qquad (10)$$

where $\dot{\chi}$ is the time rate of change of the light element concentration (mixing ratio) in the outer core due to inner core growth. The mixing ratio based Ra (10) was invoked with dominantly (or purely) chemical convection (Amit and Olson 2015; Olson et al. 2013). The final internal control parameter is $\epsilon$, the sink (or source) term that appears in the co-density equation (Christensen and Wicht 2007), which parameterizes the effects of mixing of light elements in the outer core due to inner core freezing, secular cooling of the outer core, and radioactivel heat sources.

Various alternative CMB heat flux models were imposed on the outer boundary of the dynamo models. In the “Review” section below, we describe in details the procedures for their constructions. The non-dimensional heat flux anomaly amplitude is classically defined in numerous studies by the peak-to-peak difference normalized by twice the mean value (e.g., Olson and Christensen 2002) as

$$q^* = \frac{q_{\max} - q_{\min}}{2q_0} \qquad (11)$$

When the CMB heat flux is localized (11) appears to be inappropriate, hence various alternative measures were proposed. Sreenivasan and Jellinek (2012) proposed their control parameter f as the peak-to-background difference normalized by the background heat flux $H_B$

$$f = \frac{q_B - q_{\min}}{q_B} \qquad (12)$$

They also defined the ratio of the average heat flux, either over the plume footprint or over the entire CMB, to the background heat flux, i.e. $\bar{q}_{\text{plume}}/q_B$ or $\bar{q}_{\text{CMB}}/q_B$, respectively (Sreenivasan and Jellinek 2012). Kuang et al. (2014) defined their control parameter $\epsilon_I$ as the relative post-impact heterogeneity amplitude, as in (12). In

contrast, according to the model of mantle heating by impact of Monteux et al. (2015), in the vicinity of the impact, the lower mantle is heated to the core temperature and the CMB heat flux becomes zero, whereas away from it the lower mantle temperature is unchanged. Therefore, the control parameter $q_0^r$ of Monteux et al. (2015) is based on the portion of the CMB heated by the impact

$$q_0^r = \frac{q_0^h - q_0}{q_0^h} \qquad (13)$$

where $q_0^h$ and $q_0$ are the pre- and post-impact mean CMB heat flux values, respectively.

To dimensionalize the dynamo models output, D, $D^2/\nu$, and $\sqrt{\rho\mu_0\lambda\Omega}$ were used to scale length, time, and magnetic field, respectively, where $\rho$ is the fluid density and $\mu_0$ the permeability of free space. Co-density is scaled either by $Dq_0/k$ when (9) is used or by $\rho\beta D^2 \dot{\chi}/\nu$ when Ra is defined by (10). For studying the present-day geodynamo an Earth-like inner-outer core radii ratio of $r_i/r_o = 0.35$ was applied (Amit and Choblet 2009, 2012; Amit et al. in press). A smaller ratio of $r_i/r_o = 0.2$ was applied for Mars (Monteux et al. 2015) to approach no inner core conditions while avoiding numerical instability. The past geodynamo was modeled with a variable aspect ratio based on thermal history modeling (Olson et al. 2013). Numerical dynamos with rigid insulating boundary conditions were analyzed. Overall, the models differ in the imposed outer boundary heat flux pattern, the amplitude of its variation and the convection vigor and style—all of these choices depend on the particular geophysical application.

Parametric studies of numerical dynamos found two main types of solutions: non-reversing dipole-dominated and reversing multipolar (Christensen and Aubert 2006; Kutzner and Christensen 2002; Olson and Christensen 2006). Careful choice of control parameters may yield some “Earth-like” overlap of dipole-dominated reversing dynamos (Olson 2007a; Wicht et al. 2009). In this review, studies that focus on morphological features of the present-day field had dipole-dominated non-reversing models as a reference state, whereas in studies of reversal frequency either reversing dynamos or models of both types were analyzed. Both types of solutions were examined in the context of the paleo dynamo of Mars. Table 1 summarizes the control parameters, convection styles, inner to outer core ratios, and dynamo types in the papers reviewed here.

## Review

To our best knowledge, only eight papers have examined the effects of alternative (i.e., neither tomographic nor single harmonic) CMB heat flux on planetary dynamos (Amit and Choblet 2009, 2012; Amit and Olson 2015; Amit et al. in press; Kuang et al. 2014; Monteux et al. 2015; Olson et al. 2013; Sreenivasan and Jellinek
<table><caption>**Table 1** Summary of dynamo models control parameters, geometry, convection styles and reversibility in the following papers: AC09 (Amit and Choblet 2009), AC12 (Amit and Choblet 2012), ADC15 (Amit et al. in press), ODHZ13 (Olson et al. 2013), AO15 (Amit and Olson 2015), SJ12 (Sreenivasan and Jellinek 2012), KJRF14 (Kuang et al. 2014) and MACLT15 (Monteux et al. 2015). In all papers Pr = 1. Control parameters correspond to the definitions in this paper (see "Methods" section), so the values may differ from the original papers if different definitions were used (in particular for the Ekman number). Alternative definitions of CMB heat flux heterogeneity amplitudes $q_{eff}^*$ (termed here *effective*) include f of SJ12 (12), $\epsilon_i$ of KJRF14 (also 12) and $q_i^*$ of MACLT15 (13). The Rayleigh number for SJ12 and KJRF14 is given in terms of the critical value for the onset of convection $Ra_c$. The time-dependent inner core size in ODHZ13 is determined based on a core evolution model</caption><thead><tr><th>Paper</th><th>E</th><th>Ra</th><th>Pm</th><th>q*</th><th>q<sub>eff</sub>*</th><th>r<sub>i</sub>/r<sub>o</sub></th><th>Conv</th><th>Rev?</th></tr></thead><tbody><tr><td>AC09</td><td>3 · 10<sup>-4</sup></td><td>10<sup>6</sup></td><td>3</td><td>0.5</td><td>–</td><td>0.35</td><td>TC</td><td>No</td></tr><tr><td>AC12</td><td>3 · 10<sup>-4</sup></td><td>(0.5 – 1) · 10<sup>6</sup></td><td>3 – 7</td><td>0.5 – 0.8</td><td>–</td><td>0.35</td><td>TC</td><td>No</td></tr><tr><td>ADC15</td><td>3 · 10<sup>-4</sup></td><td>(1 – 3) · 10<sup>6</sup></td><td>3 – 10</td><td>0.7 – 1.34</td><td>–</td><td>0.35</td><td>TC</td><td>No</td></tr><tr><td>ODHZ13</td><td>5.75 · 10<sup>-3</sup></td><td>(1.5 – 4) · 10<sup>4</sup></td><td>20</td><td>0.67 – 1.1</td><td>–</td><td>Variable</td><td>C</td><td>Both</td></tr><tr><td>AO15</td><td>(1 – 6) · 10<sup>-3</sup></td><td>(0.6 – 18) · 10<sup>5</sup></td><td>20</td><td>0.65 – 2.22</td><td>–</td><td>0.35</td><td>C</td><td>Both</td></tr><tr><td>SJ12</td><td>5 · 10<sup>-5</sup></td><td>12Ra<sub>c</sub></td><td>1</td><td>–</td><td>4 – 29</td><td>0.35</td><td>TC</td><td>No</td></tr><tr><td>KJRF14</td><td>2.6 · 10<sup>-5</sup></td><td>~ Ra<sub>c</sub></td><td>1</td><td>–</td><td>0.05 – 0.3</td><td>0.31</td><td>T</td><td>Yes</td></tr><tr><td>MACLT15</td><td>(1 – 3) · 10<sup>-4</sup></td><td>(0.735 – 4) · 10<sup>7</sup></td><td>3</td><td>0.3 – 0.8</td><td>0.04 – 0.26</td><td>0.2</td><td>T</td><td>Both</td></tr></tbody></table>

Conv convection type, T thermal, TC thermochemical, C chemical, Rev? reversibility
2012). Here we describe for each of these papers the CMB heat flux models that were imposed as outer boundary conditions on the dynamo simulations. Then, for each paper, we compare the results with those obtained with a reference heat flux model (Table 2) and we report the consequences for the resulting dynamo properties.

## Lower mantle thermal-chemical-phase heterogeneity and geodynamo morphological features

Here we review the papers by Amit and Choblet (2009), Amit and Choblet (2012), and Amit et al. (in press). As in previous studies, a reference CMB heat flux anomaly of purely thermal origin $\delta q_T$ was assumed to be linearly proportional to the lowermost mantle seismic shear velocity anomaly $\delta v_s$:

$$ \frac{\delta q_T}{q_T} \propto \frac{\delta v_s}{v_s} \qquad (14) $$

Figure 3a shows the CMB heat flux based on the tomography model of Masters et al. (2000) and (14). Nakagawa and Tackley (2008) found that when incorporating post-perovskite phase transition in their mantle convection models, the $\delta q$ to $\delta v_s$ relation becomes an inclined step function (Fig. 2a), with the spreading effect well approximated by a tanh function. Following equation (B8) of Nakagawa and Tackley (2008), Amit and Choblet (2009) attempted to parameterize the thermal-phase interpretation $\delta q_{Tp}$ as

$$ \frac{\delta q_{Tp}}{q_{Tp}} = \frac{\delta q_T}{q_T} + A \left[ 0.5 \left( \tanh \left( -\frac{1}{B} \frac{\delta q_T}{q_T} \right) - 1 \right) \right] + D \quad (15) $$

<table><caption>**Table 2** Summary of dynamo models imposed outer boundary heat flux patterns in AC09 (Amit and Choblet 2009), AC12 (Amit and Choblet 2012), ADC15 (Amit et al. in press), ODHZ13 (Olson et al. 2013), AO15 (Amit and Olson 2015), SJ12 (Sreenivasan and Jellinek 2012), and MACLT15 (Monteux et al. 2015). For each paper, the reference and alternative patterns are given. Tomographic always refers to present-day, HF1 is the time-dependent model of Zhang and Zhong (2011). For more details, see text</caption><thead><tr><th>Paper</th><th>Reference</th><th>Alternative</th></tr></thead><tbody><tr><td>AC09</td><td>Tomographic</td><td>Thermal + post-perovskite</td></tr><tr><td>AC12</td><td>γ<sub>2</sub><sup>2</sup></td><td>γ<sub>2</sub><sup>2</sup> + ridges</td></tr><tr><td>ADC15</td><td>Tomographic</td><td>Probabilistic tomography</td></tr><tr><td>ODHZ13</td><td>Tomographic</td><td>HF1</td></tr><tr><td>AO15</td><td>Homogeneous</td><td>Rotated γ<sub>2</sub><sup>0</sup></td></tr><tr><td>SJ12</td><td>Homogeneous</td><td>Plume-induced</td></tr><tr><td>KJRF14</td><td>Homogeneous</td><td>Impact-induced</td></tr><tr><td>MACLT15</td><td>Degree-1</td><td>Impact-induced</td></tr></tbody></table>

The non-dimensional parameters in (15) are the amplitude A, the tanh spreading factor B, and the offset D. For a given amplitude A, the spreading factor B is constrained by reaching the tanh asymptote while maintaining an inclined step function (rather than Z-like curve).

Figure 3b shows a heat flux anomaly map that accounts for post-perovskite according to (15). The patterns in Fig. 3a, b are similar, but the relative magnitudes of the positive/negative heat flux structures vary. When the post-perovskite effect is included, the positive structures weaken and spread over larger areas, whereas the negative structures remain strong and more localized (Fig. 3b). This is due to the particular pattern of the tomographic model that contains larger amplitude narrower regions of negative seismic velocity anomalies ($\delta v_s$, $\delta q_T < 0$) as opposed to broader regions of weaker positive anomalies (Fig. 3a). The tanh function strongly reduces intermediate $\delta v_s$ while it reduces extreme $\delta v_s$ more mildly, resulting in a spreading effect acting stronger on the positive $\delta q_T$.

Following Aubert et al. (2008), Amit and Choblet (2009) compared three long-term time-average properties of their dynamo models with relevant observations: the radial magnetic field on the outer boundary with the time-average paleomagnetic field over the past 5 Myr (Kelly and Gubbins 1997), the flow at the top of the shell with the core flow model of Amit and Christensen (2008) obtained from inversions of the historical geomagnetic secular variation (Jackson et al. 2000) and averaged over the period 1840–1990, and the inner boundary buoyancy flux with seismic anomalies of the upper inner core (Tanaka and Hamaguchi 1997). As an example, we highlight here the latter comparison (Fig. 4). The inner boundary buoyancy flux is characterized by two isolated positive low-latitude structures both in the thermal case (Fig. 4a, denoted as T) as well as when accounting for post-perovskite (Fig. 4b, denoted as Tp1). In addition, however, the ICB buoyancy fluxes contain other orders. More vigorous convection leads to stronger filtering of the CMB order 2 content and the ICB pattern remains dominantly order 1 (Aubert et al. 2008), whereas slightly supercritical convection leaves the ICB pattern with nearly non-distorted order 2 (Gubbins et al. 2011). Amit and Choblet (2009) found that the ratio of order 1 to order 2 in Fig. 4b is significantly larger than that in 4a. This suggests that accounting for post-perovskite improves the agreement of the dynamo models with the order 1 dominance in seismic properties of the upper inner core (Tanaka and Hamaguchi 1997), although the phase may vary substantially from one dynamo model to another. Amit and Choblet (2009) argued that accounting for post-perovskite also improves the recovery of the Atlantic/Pacific hemispheric dichotomy in core flow activity (Amit and Christensen 2008; Gillet et al. 2009) and the single intense
<figure><img src="image_226.png" /><figcaption>**Fig. 3** Non-dimensional CMB heat flux patterns based on the mantle tomography model of Masters et al. (2000): **a** Thermal linear interpretation T and **b** interpretation including a post-perovskite effect Tp1. From Amit and Choblet (2009)</figcaption></figure>

paleomagnetic flux patch in the southern hemisphere (Kelly and Gubbins 1997).

Next, narrow features that are not captured by global tomography models were considered. Lower mantle tomography models are dominated by a $Y_2^2$ component (e.g., Masters et al. 2000), which is classically interpreted as the largest scale temperature anomaly. Such a pattern was imposed by many authors as heat flux on the outer boundary of their numerical dynamo models (e.g., Aubert et al. 2007; Bloxham 2002; Olson and Christensen 2002; Sreenivasan 2009). Amit and Choblet (2012) used the $Y_2^2$ pattern as a reference case and superimposed on it low heat flux ridges at the four margins of hot and cold large-scale structures (Table 2). The ridges were modeled using a Gaussian geometry with an amplitude equal to that of the hot large-scale features (Fig. 5).

Figure 6 presents two close-ups of the time-average dynamo model velocity and magnetic field, one on the center of Fig. 5 (i.e., center of a ridge with a positive heat flux anomaly to its east), the other centered on a ridge with a positive heat flux anomaly to its west. In places where a positive heat flux anomaly appears east of the ridge (Fig. 6 left), the poloidal flow dominates and the fluid diverges

from the hot ridge. Stronger field appears east of the ridge where cold fluid downwelling concentrates magnetic flux. In contrast, in places where a positive heat flux anomaly appears west of the ridge (Fig. 6 right), the upwelling associated with this ridge is broken at the equator, westward toroidal flow dominates, and the magnetic field exhibits patches of comparable intensity on both sides of the ridge. Note that the fluid upwellings appear very close to the longitudes of the hot CMB heat flux ridges, in contrast to previous studies of numerical dynamos with a $Y_2^2$ CMB heat flux pattern in which a phase shift was found between the heat flux and the flow structures (Aubert et al. 2007; Olson and Christensen 2002; Takahashi et al. 2008).

Amit and Choblet (2012) interpreted the low-latitude dynamics around the ridge locations by a combined effect of westward drift in the time-average homogeneous dynamo (Aubert 2005) and the steady boundary-driven flow that roughly goes from hot to cold regions (Gubbins 2003). When cold material is east of the ridge the two flows counteract and the ridge acts as a flow barrier, whereas when cold material is west of the ridge both flows are westward and the barrier effect is diminished. Applying these inferences to the tomography model of

<figure><img src="image_227.png" /><figcaption>**Fig. 4** Time-average inner-core boundary buoyancy flux anomalies normalized by the mean buoyancy flux: **a** Model T and **b** Model Tp1. From Amit and Choblet (2009)</figcaption></figure>
<figure><img src="image_228.png" /></figure>

Masters et al. (2000), Amit and Choblet (2012) predicted (in a statistical sense) fluid upwelling at the top of the core and a barrier for azimuthal flow below central Asia and the Indian Ocean and below the American Pacific coast. Below east Asia and Oceania and below the Americas

they predicted persistent intense low-latitude geomagnetic flux patches. These predictions are in agreement with preferred longitudes of virtual geomagnetic poles during reversals both in paleomagnetic data (Hoffman 1996; Love 1998; Merrill and McFadden 1999) and in tomographic reversing dynamo models (Coe et al. 2000; Kutzner and Christensen 2004). It also resonates with the persistent production of low-latitude geomagnetic flux patches below Indonesia during the historical era (Jackson et al. 2000).

Finally, probabilistic tomography may be used to isolate the thermal component of the observed seismic anomalies. A limitation of the probabilistic tomography models published so far is that only even spherical harmonic degrees are included (Mosca et al. 2012; Trampert et al. 2004). To overcome this problem, Amit et al. (in press) mapped the ratio between shear-wave velocity $\delta lnV_s^e$ (even degrees) and thermal anomalies $\delta T^e$ (even degrees) in probabilistic tomography and applied this mapping to a conventional lower mantle seismic tomography model $\delta lnV_s^C$ (even and odd degrees) to obtain a rescaled thermal anomaly distribution $\delta T^P$ with both even and odd degrees. Mathematically they wrote

$$ \delta T^P = \frac{\delta T^e}{\delta lnV_s^e} \delta lnV_s^C \quad (16) $$

<figure><img src="image_229.png" /></figure>
To avoid numerical instabilities at regions where $\delta\ln V_s^e$ is low, at each node of the probabilistic tomography model $\delta T^e$ and $\delta\ln V_s^e$ were sampled according to their full probability distribution functions (PDFs) rather than their mean values, resulting in a full PDF of the ratio $\delta T^e/\delta\ln V_s^e$. Amit et al. (in press) applied the above procedure on the conventional seismic tomography model of Masters et al. (2000) and the probabilistic tomography model of Mosca et al. (2012) for the seismic and thermal anomalies. The resulting CMB heat flux pattern inferred from probabilistic tomography (Fig. 7) exhibits positive structures mostly concentrated at low latitude.

The dynamo models with an imposed probabilistic tomography heat flux pattern are characterized by stronger magnetic and convective activity at low-latitude than the dynamo models with conventional tomography heat flux. The timeseries of the ratio of RMS low-latitude magnetic flux to RMS high-latitude magnetic flux as well as the ratio of maximum low-latitude magnetic flux to maximum high-latitude magnetic flux (Fig. 8) show larger mean values and much larger temporal fluctuations in the probabilistic dynamo models (Fig. 8b) than those obtained with conventional tomography (Fig. 8a). For comparison, Amit et al. (in press) found that the historical geomagnetic field model *gufm1* (Jackson et al. 2000) exhibits for the period 1840–1990 average RMS and maxima ratios of 0.675 and 0.725, respectively. The higher resolution model *CHAOS-4* (Olsen et al. 2014) for the period 2000–2010 is characterized by somewhat larger average ratios of 0.797 and 0.882, and the average ratios in the lower resolution archeomagnetic field model *Cals3k.3* (Korte et al. 2009) for the past three millennia are 0.470 and 0.648. Accounting for the temporal fluctuations in the dynamo models with the CMB heat flux pattern inferred from probabilistic tomography, any of these models may occasionally recover the observed low- to high-latitude geomagnetic flux ratios.

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

### Lower mantle history and geomagnetic reversal frequency

Here we review the papers by Olson et al. (2013) and Amit and Olson (2015). Zhang and Zhong (2011) imposed a reconstructed plate tectonics model for the period 0–120 Ma (Lithgow-Bertelloni and Richards 1998) and Atlantic hemisphere continent reconstructions for older periods (Scotese 2001) as a time-dependent mechanical upper boundary condition on their mantle thermochemical convection model to infer the time-dependent CMB heat flux. According to one of the CMB heat flux models of Zhang and Zhong (2011), termed HF1, relatively large/low mean CMB heat flux corresponds to fast/slow plate reconstructions (Lithgow-Bertelloni and Richards 1998; Scotese 2001), respectively. In addition, HF1 includes temporal changes in the heat flux pattern (Fig. 9) and amplitude (Table 1), most notably a transition from degree-1 dominance during Pangea to degree-2 dominance at present.

Olson et al. (2013) imposed snapshots of HF1 (Fig. 9) on chemical convection dominated numerical dynamos. In addition, they accounted for the time-dependent core evolution that is controlled by the mean CMB heat flux. Using thermal history modeling (Labrosse 2003), Olson et al. (2013) determined the time-dependent CMB heat flux and consequently the time-dependent inner core size and inner core boundary buoyancy flux. According to these models, a few tens of Myrs after inner core nucleation outer core convection indeed becomes dominantly chemical. Very high core thermal conductivity (e.g., Pozzo et al. 2012) implies that higher CMB heat flux, radiogenic heating, or some presently unknown form of chemical differentiation in the core would be needed to power the geodynamo prior to inner core nucleation. However, the debate on core thermal conductivity is still ongoing (Zhang et al. 2015). Finally, Olson et al. (2013) also accounted for the observed nearly linear increase in length of day (Williams 2000) to derive a time-dependent rotation rate. Overall, they established the relative temporal evolution of the Rayleigh and Ekman numbers for the period corresponding to the HF1 model.

Olson et al. (2013) compared the time-dependent reversal frequency of two types of dynamo models, one with time-dependent and another with fixed present-day CMB heat flux (both cases include time-dependent core evolution). This comparison allows distinguishing between the direct effects of the CMB heterogeneity and the effects of the core evolution. Figure 10 shows comparison of the reversal frequency of the dynamo models of Olson et al. (2013) with the observed paleomagnetic reversal frequency record (Gradstein et al. 2012). The dynamo model with fixed CMB heat flux is characterized by weaker temporal variability in reversal frequency. In contrast, the dynamo model with time-dependent CMB heat
<figure><img src="image_231.png" /></figure>

<figure><img src="image_232.png" /><figcaption>**Fig. 8** Timeseries of low/high-latitude magnetic flux ratios for dynamo models with the same control parameters but different outer boundary heat flux patterns, one with a pattern from conventional tomography (a, case C2 of Amit et al. in press) and the other with a pattern inferred from probabilistic tomography (b, case P2 of Amit et al. in press). Black/red curves are RMS/maxima ratios, respectively. From Amit et al. (in press)</figcaption></figure>

flux shows larger reversal frequency fluctuations including superchrons and rather hyper-reversing periods. In particular, at 275 Ma when Pangea was at its height and prior to it, the dynamo model with time-dependent CMB heat flux recovers the low reversal frequency, whereas that is not the case for the dynamo model with fixed CMB heat flux. Olson et al. (2013) argued that the transition from low to high reversal frequency is driven by the superconductive cycle. Both dynamo models fail to reproduce the robust Cretaceous Normal Superchron at 120 Ma. However, Olson et al. (2013) demonstrated that the reversal frequency is strongly sensitive to mild changes in the CMB heat flux (mean, pattern, and its amplitude).

In contrast, the lower mantle may exhibit a quasi-stationary pattern. Amit and Olson (2015) proposed that the extreme variability in the recorded geomagnetic reversal frequency, from hyper-reversing periods to superchrons, may be explained by lower mantle piles with

a fixed lateral pattern but time-dependent height. Amit and Olson (2015) proposed a simple model of the D" layer with fixed temperature difference, conserved volume with time, and a conductive equilibrium. Based on these assumptions, they derived analytically the dependence of the mean CMB heat flux and the amplitude of its lateral heterogeneity on the pile height. They showed that superplume growth increases both the mean CMB heat flux and the amplitude of its lateral heterogeneity, whereas pile collapse decreases these two quantities.

Figure 11 shows the rotated $Y_2^0$ CMB heat flux pattern based on the model of lower mantle seismic heterogeneity by Dziewonski et al. (2010). The mean heat flux was assumed to be close to the core adiabat, i.e., the dynamo models are dominated by chemical convection. The heterogeneity pattern was generated using spherical harmonics $Y_2^0, Y_2^2$, and $Y_1^1$, with relative amplitudes of 10 : 10 : 1, respectively. Amit and Olson (2015) determined the