<figure><img src="image_378.png" /><figcaption>**Fig. 1.** Schematic representation at scale of the impact pressure field, as considered in this study, for two impactor sizes: $d_i = 5$ km and 20 km and two impactor velocities: $v_i = 11$ km s$^{-1}$ and 15 km s$^{-1}$, hitting a NH$_3$-enriched icy layer of thickness $h_c$. The impact-induced pressure field remains nearly constant in a region beneath the impact site, called isobaric core, characterized by its radius $R_{ic}$ and the depth of its center below the surface, $z_{ic}$. The comparison between case (1) and case (2) shows the influence of the impact velocity on the impact-induced pressure field: owing to higher impact velocity in case (1), the pressure inside the isobaric core is higher and the decay pressure zone is larger. As shown by comparing case (2) and case (3), a greater impactor increases the size of the isobaric core but does not affect the peak pressure.</figcaption></figure>

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

Physical constants for the numerical modelisation.</caption><thead><tr><th>Parameter</th><th>Name</th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td>R<sub>T</sub></td><td>Titan's radius</td><td>2575 km</td><td></td></tr><tr><td>g<sub>T</sub></td><td>Titan's gravity</td><td>1.35 m s<sup>-2</sup></td><td></td></tr><tr><td>v<sub>esc</sub></td><td>Escape velocity</td><td>2639 m s<sup>-2</sup></td><td></td></tr><tr><td>M<sub>TA</sub></td><td>Titan's present day atmospheric mass</td><td>9.2 &middot; 10<sup>18</sup> kg</td><td>Niemann et al. (2010)</td></tr><tr><td>&rho;<sub>t</sub></td><td>Density of the icy crust (target)</td><td>917 kg m<sup>-3</sup></td><td>Stewart (2005)</td></tr><tr><td>&rho;<sub>i</sub></td><td>Density of the impactor</td><td>917 kg m<sup>-3</sup></td><td>Stewart (2005)</td></tr><tr><td>C<sub>t</sub></td><td>Bulk sound wave velocity in water ice</td><td>3610 m s<sup>-1</sup></td><td>Stewart (2005)</td></tr><tr><td>S<sub>t</sub></td><td>Non dimensional constant</td><td>0.92</td><td>Stewart (2005)</td></tr></tbody></table>

channel with a density that is 1–2 orders of magnitude lower than the ambient air density. The heated gas flow expands outward through this channel and is possibly ejected into space. (ii) When hitting the surface, an expanding plume develops above the impact site, where the gas molecules are also heated and can be accelerated upward to escape the gravity of Titan (see Shuvalov (2009) for more details). We describe these atmospheric erosion processes using the parameterization of Shuvalov (2009), initially built up for the Earth. This parameterization is valid for a broad range of atmospheric densities and impactor sizes compatible with LHB population and it has been adapted here to Titan.

According to hydrodynamic simulations of Shuvalov (2009), the normalized eroded mass $\chi_a$ is related to the impact erosion efficiency $\xi$ defined respectively as:

$$ \chi_a = \frac{m_a}{m_i} \frac{v_{esc}^2}{v_i^2 - v_{esc}^2}, \qquad (8) $$

$$ \xi = \frac{d_i^3 \rho_i}{H^3 \rho_{atm}} \frac{(v_i^2 - v_{esc}^2)}{v_{esc}^2} \frac{\rho_t}{(\rho_t + \rho_i)}, \qquad (9) $$

with $m_a$ the atmospheric mass escaped during the impact, $m_i$ the impactor mass, $d_i$ the impactor diameter and $v_{esc}$ the escape velocity. By fitting simulation results obtained on a broad range of $\chi_a$ and $\xi$, Shuvalov (2010) proposed the following polynomial relationship:

$$ \log_{10}(\chi_a) = -6.375 + 5.239(\log_{10}(\xi)) - 2.121(\log_{10}(\xi))^2 \\ + 0.397(\log_{10}(\xi))^3 - 0.037(\log_{10}(\xi))^4 \\ + 0.0013(\log_{10}(\xi))^5. \qquad (10) $$

To adapt this parameterization to Titan, we compute $\xi$, then using Eq. (10) we find $\chi_a$ and finally $m_a$ by replacing the parameters appropriate for Titan in each expression. The range of $\xi$ values we encounter during the simulations is comparable to the range explored in Shuvalov (2009). In Fig. 2 we compare the atmospheric erosion obtained with the parametrization of Shuvalov (2009) to the hydrocode simulations of impact-induced erosion of Korycansky and Zahnle (2011), a study devoted to Titan. Fig. 2 presents the total amount of atmosphere that is eroded on Titan following the parameterization of Shuvalov (2009), as a function of the impactor diameter and for three different atmospheric pressures (15 bar, 1.5 bar, 0.01 bar) and an impact velocity of 10 km s$^{-1}$. For the case with 1.5 bar, we compare our calculations to the values obtained by Korycansky and Zahnle (2011) from hydrocode simulations. The parameterization of Shuvalov (2009) consider an eroded mass averaged over the impact angle. When compared to the simulations of Korycansky and Zahnle (2011) for fixed impact angle of 45° and 75°, and impactor diameters of 10 km, the atmospheric mass loss is comparable for the two studies. For larger impactors, there are some discrepancies between the two models, but the magnitudes remain comparable. Fig. 2 also illustrates that the erosion decreases rapidly with decreasing impactor sizes. This results in a more efficient deceleration of the smallest impactors by atmospheric drag, leading to a smaller
<figure><img src="image_379.png" /><figcaption>**Fig. 2.** Eroded atmospheric mass as a function of the impactor diameter for the model of Shuvalov (2009, 2010), for a $v_0 = 10$ km s$^{-1}$. Korycansky and Zahnle (2011) values are given for a present-day Titan's atmosphere, for a 2D plane-polar geometry calculations on a 800 × 1280 grid and $v_0 = 10$ km s$^{-1}$.</figcaption></figure>

impact-induced plume. For impactors larger than ~50 km, the eroded atmospheric mass is limited by an asymptotic value corresponding to the atmospheric mass above the plane tangent to the satellite's surface.

## 2.4. Evolution by multi-impact

We implement the parametrizations described above into a multi-impact approach, using an impactor flux characteristic of the LHB event. Then, we integrate the evolution of the volatile reservoir throughout the LHB-era.

### 2.4.1. Multi-impact approach

To generate our LHB impactor population, we carried out a Monte-Carlo calculation. Sampling procedures for the impactor diameter, velocity and coordinates of the impact sites are adapted from Zahnle (2001) and Zahnle et al. (2003). Impactor diameters are chosen from a size distribution characteristic of the LHB and scaled to match the crater record on lapetus (Charnoz et al., 2009). The velocity distribution follows a Gaussian distribution (Korycansky and Zahnle, 2011). Assuming an isotropic impact flux, the latitude and longitude of the impact site are randomly drawn. Parameters used to define the distributions are summarized in Table 3.

To limit the computation time, we assume a lower limit of 1 m for impactor diameters. We also prescribe an upper limit of 100 km. Above this value, the validity of the scaling laws used here becomes questionable because they do not include physical parameters that may influence the impact processes at large scales, such as the curvature of the satellite and its internal structure (Louzada and Stewart, 2009; Bruesch and Asphaug, 2004).

<caption>**Table 3**</caption>

Impactor sampling parameters. The Gaussian distribution of the impactor velocity $v_0$ is taken from Korycansky and Zahnle (2005).

<table><thead><tr><th>Parameter</th><th>Distribution</th><th>Interval of sampling</th></tr></thead><tbody><tr><td>d<sub>i</sub> – impactor diameter</td><td>From Charnoz et al. (2009)</td><td>1 m to 100 km</td></tr><tr><td>v<sub>0</sub> – initial velocity of the impactor</td><td>Gaussian</td><td>11.3 km s<sup>-1</sup> ± 4 km s<sup>-1</sup> (1σ)</td></tr><tr><td>λ – longitude</td><td>Linear</td><td>0–360°</td></tr><tr><td>φ – latitude</td><td>cos(φ)</td><td>-90° to 90°</td></tr></tbody></table>

For each impact, we compute the velocity decrease due to atmospheric drag using Eq. (2), we determine the impact pressure field and corresponding degassing using Eqs. (1), (5) and (7) and the atmospheric mass loss using Eqs. (8)–(10). The new atmospheric equilibrium as well as the crustal NH$_3$ distribution are re-evaluated after each impact. The simulations are stopped when the total mass of impactors delivered to Titan reaches a value of ~3 × 10$^{20}$ kg (Barr et al., 2010).

### 2.4.2. Subsurface distribution and numerical evolution of NH$_3$ reservoirs

The distribution of NH$_3$ at the time of the LHB is mostly governed by Titan's degree of differentiation, which is strongly related to its accretional history (e.g. Monteux et al., 2014). NH$_3$ could be either homogeneously distributed within an undifferentiated interior or mostly located in a chemically differentiated crust. To explore different possibilities, we consider different initial NH$_3$ distributions and simplified evolutions of NH$_3$ reservoirs summarized in the three following cases:

* Case A: we consider a uniform and constant distribution of NH$_3$ all along the simulation. This simplified case considers an infinite reservoir of NH$_3$.

* Case B: the reservoir is assumed initially uniform throughout a given thickness (up to 100 km), and the remaining fraction of ammonia is reevaluated after each impact as a function of depth. This one-dimensional approach allows us to take into account the averaged progressive depletion in ammonia as a function of depth during the simulation.

* Case C: a reservoir of limited depth (less than 5 km) is modeled as a spherical shell with a 0.1° spatial resolution in latitude and longitude. The spatial distribution of NH$_3$ is re-evaluated in each cell after an impact, assuming that the lost fraction of NH$_3$ is uniform over the cell volume. This two-dimensional approach allows us to record the local depletion in NH$_3$, but do not take into account the variations with depth. Fig. 3 shows an example of such distribution in the end of the simulation (no. 14, see Table 5). A full 3D evolution of the reservoir, including lateral and depth variation would be too much CPU time and therefore has not been considered here.

Case A considering an infinite NH$_3$ reservoir is used only to compare our results to those of Sekine et al. (2011). Case B can be used either for undifferentiated or differentiated interior, depending on the assumed thickness of NH$_3$ reservoir. Case C is valid only for relatively thin NH$_3$ reservoir (≤5 km), and therefore represents a differentiated interior case with a crust enriched in NH$_3$. For these three cases, the volume affected by impact is computed from Eqs. (1), (3), (4) and (5). The only difference between these three cases concerns the evolution of NH$_3$ distribution for each reservoir after an impact: no NH$_3$ depletion in Titan's interior for the case A, NH$_3$ depletion with depth for the case B and NH$_3$ depletion at Titan's surface, as displayed for case C in Fig. 3(a).

As the surface and atmospheric conditions on early Titan are unconstrained, we performed numerous simulations (49 in total) to cover a wide range of possible parameter values, such as the surface albedo, the size and concentration of crustal NH$_3$ reservoir, and the mass and composition of the pre-LHB atmosphere (see Table 4). For the surface albedo, we considered values ranging between 0.1 and 0.7. Except Enceladus and Triton, most of the icy moons has albedos comprised between 0.1 and 0.6. The present-day surface albedo on Titan is estimated between 0.1 and 0.2, its dark surface being mostly the consequences of deposition of photochemically-produced organics. In absence of methane in the past, no photochemical byproducts would accumulate, and
<figure><img src="image_380.png" /><figcaption>Fig. 3. (a) Distribution of the NH₃ at Titan's surface after the simulation no. 14 in Table 5: no initial atmosphere was considered here; (b) distribution of cumulated impact sizes during the LHB for the same simulation; and (c) distribution of affected terrains after the LHB for the map (a) and the simulation no. 26 in Table 5.</figcaption></figure>

<table><caption>**Table 4**
List of parameters explored in this study and the corresponding range of values.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Range of values explored in this study</th></tr></thead><tbody><tr><td>Surface albedo</td><td>A<sub>s</sub></td><td>0.1–0.7</td></tr><tr><td>Initial atmospheric mass</td><td>M<sub>ini</sub><sup>A</sup></td><td>0–10 M<sub>TA</sub></td></tr><tr><td>CH<sub>4</sub> fraction in the impactor</td><td>X<sub>CH<sub>4</sub></sub><sup>A</sup></td><td>0% or 1%</td></tr><tr><td>Initial CH<sub>4</sub> fraction in the atmosphere</td><td>X<sub>CH<sub>4</sub></sub><sup>A</sup></td><td>0–10%</td></tr><tr><td>Mass concentration of NH<sub>3</sub> in the icy crust</td><td>X<sub>NH<sub>3</sub></sub></td><td>1–5%</td></tr><tr><td>Mass concentration of NH<sub>3</sub> in the impactor</td><td>X<sub>NH<sub>3</sub></sub><sup>A</sup></td><td>1%</td></tr><tr><td>Thickness of the NH<sub>3</sub>-enriched icy layer</td><td>h<sub>c</sub></td><td>∞ for case A<br>5–50 km for case B<br>1–5 km for case C</td></tr></tbody></table>

therefore its surface may have been much brighter, that's why we considered values as high as 0.7. We also investigated different scenarios for NH₃ enrichment in the outer layers corresponding to different differentiation state at the time of LHB, with ammonia mass fraction varying between 1% and 5%. As initial conditions, we considered either no atmosphere or the presence of an pre-LHB atmosphere with a mass up to 10 times present-day Titan's atmospheric mass and with methane mass fraction up to 10%. The corresponding results are presented in Table 5 for simulations with no preexisting atmosphere and in Table 6 for the simulations with pre-LHB atmospheres.

### 3. Results

In the following, we first present the atmospheric mass balances between atmospheric loss and supply for individual impacts. Then we investigate the possible emergence of an atmosphere

from an initially air-less body and after the evolution of initially massive atmosphere for multiple impacts characteristic of the LHB.

#### 3.1. Balance between atmospheric loss and supply

Fig. 4 presents the balance between N₂ crustal degassing and atmospheric erosion for individual impacts as a function of the impactor diameter for three different atmospheric pressures (Fig. 4a), impactor velocities (Fig. 4b) and crustal reservoirs (Fig. 4c). Our calculations show that the maximal erosion rate is obtained for impactor diameters of 2.5, 12 and 30 km for atmospheric pressures of 0.015, 1.5 and 15 bars, respectively, explaining the minima in the N₂ balance (Fig. 4a).

For impact diameters smaller than these optimal values, the efficiency of erosion rapidly decays, as the balance is mostly controlled by the degassing efficiency. Owing to both decrease in atmospheric erosion and increase in N₂ degassing relatively to the impactor mass, the balance becomes positive for diameters lower than 1.3, 6, and 13 km for atmospheric pressures of 0.015, 1.5 and 15 bars, respectively (see Fig. 4a). However, as impactors smaller than 13 km represent only 11% of the total mass of the LHB impact flux, erosion should prevail.

For even smaller impactor diameters, the atmospheric drag significantly decelerates the impactors leading to reduced impactor velocities and hence to lower shock pressures and degassing rates. Below diameters of 1.2 and 6 km for surface pressures of 1.5 and 15 bars respectively (Fig. 4a), the atmospheric balance is zero because both degassing and impact-induced erosion become negligible. As decreasing impactor velocity leads to a lower degassing efficiency and lower atmospheric erosion, the atmospheric balance curves in Fig. 4b are more and more flattened with decreasing velocity and the minima and maxima shift toward larger impactors. Naturally, as shown in Fig. 4c the balance is positive for a wider range of diameters if larger reservoirs of NH₃ are considered.
<caption>**Table 5**</caption>

List of parameters and results for an initially air-less Titan simulations. $M_A$ is the final atmospheric mass, $M_c$ is the final mass of condensed volatiles at Titan's surface, $M_{TA}$ is the present-day Titan's atmospheric mass.

<table><caption>**Table 6**</caption><thead><tr><th>No.</th><th>A<sub>s</sub></th><th>CH<sub>4</sub></th><th>NH<sub>3</sub><br>reservoir</th><th>X<sub>NH<sub>3</sub></sub><br>(%)</th><th>h<sub>c</sub><br>(km)</th><th>M<sub>A</sub>/M<sub>TA</sub><br>(% of CH<sub>4</sub>)</th><th>M<sub>c</sub>/M<sub>TA</sub><br>(% of CH<sub>4</sub>)</th></tr></thead><tbody><tr><td>1</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case A</td><td>1</td><td>-</td><td>0.06</td><td>0</td></tr><tr><td>2</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case A</td><td>1.5</td><td>-</td><td>0.14</td><td>0</td></tr><tr><td>3</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case A</td><td>2</td><td>-</td><td>0.24</td><td>0</td></tr><tr><td>4</td><td>0.3</td><td>no CH<sub>4</sub></td><td>Case A</td><td>2</td><td>-</td><td>0.31</td><td>0</td></tr><tr><td>5</td><td>0.4</td><td>no CH<sub>4</sub></td><td>Case A</td><td>2</td><td>-</td><td>0.30</td><td>0</td></tr><tr><td>6</td><td>0.5</td><td>no CH<sub>4</sub></td><td>Case A</td><td>2</td><td>-</td><td>0.27</td><td>0.1<br>(liquid)</td></tr><tr><td>7</td><td>0.6</td><td>no CH<sub>4</sub></td><td>Case A</td><td>2</td><td>-</td><td>0.14</td><td>0.95<br>(liquid)</td></tr><tr><td>8</td><td>0.7</td><td>no CH<sub>4</sub></td><td>Case A</td><td>2</td><td>-</td><td>0.06</td><td>1.65<br>(solid)</td></tr><tr><td>9</td><td>0.2</td><td>1% in the<br>impactor</td><td>Case B</td><td>5</td><td>5</td><td>0.03 (26%)</td><td>0</td></tr><tr><td>10</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case B</td><td>1</td><td>50</td><td>0.03</td><td>0</td></tr><tr><td>11</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case B</td><td>1.5</td><td>50</td><td>0.06</td><td>0</td></tr><tr><td>12</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case B</td><td>2</td><td>50</td><td>0.10</td><td>0</td></tr><tr><td>13</td><td>0.5</td><td>no CH<sub>4</sub></td><td>Case B</td><td>2</td><td>50</td><td>0.14</td><td>0</td></tr><tr><td>14</td><td>0.7</td><td>no CH<sub>4</sub></td><td>Case B</td><td>2</td><td>50</td><td>0.06</td><td>0.91<br>(solid)</td></tr><tr><td>15</td><td>0.2</td><td>1% in the<br>impactor</td><td>Case B</td><td>5</td><td>50</td><td>0.93 (6%)</td><td>0</td></tr><tr><td>16</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case B</td><td>1</td><td>100</td><td>0.04</td><td>0</td></tr><tr><td>17</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case B</td><td>1.5</td><td>100</td><td>0.08</td><td>0</td></tr><tr><td>18</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case B</td><td>2</td><td>100</td><td>0.17</td><td>0</td></tr><tr><td>19</td><td>0.2</td><td>1% in the<br>impactor</td><td>Case B</td><td>5</td><td>100</td><td>1.43 (5%)</td><td>0</td></tr><tr><td>20</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case C</td><td>5</td><td>5</td><td>0.02</td><td>0</td></tr><tr><td>21</td><td>0.2</td><td>no CH<sub>4</sub></td><td>Case C</td><td>2</td><td>5</td><td>0.006</td><td>0</td></tr><tr><td>22</td><td>0.2</td><td>1% in the<br>impactor</td><td>Case C</td><td>5</td><td>5</td><td>0.03 (24%)</td><td>0</td></tr><tr><td>23</td><td>0.3</td><td>1% in the<br>impactor</td><td>Case C</td><td>5</td><td>5</td><td>0.04 (27%)</td><td>0</td></tr><tr><td>24</td><td>0.5</td><td>1% in the<br>impactor</td><td>Case C</td><td>5</td><td>5</td><td>0.03 (16%)</td><td>0.15<br>(94%, liquid)</td></tr><tr><td>25</td><td>0.5</td><td>no CH<sub>4</sub></td><td>Case C</td><td>5</td><td>5</td><td>0.02</td><td>0</td></tr><tr><td>26</td><td>0.7</td><td>no CH<sub>4</sub></td><td>Case C</td><td>5</td><td>5</td><td>0.03</td><td>0</td></tr><tr><td>27</td><td>0.7</td><td>no CH<sub>4</sub></td><td>Case C</td><td>2</td><td>5</td><td>0.01</td><td>0</td></tr></tbody></table>

List of parameters and results for simulations with an initial atmosphere. $X_{CH_4}^t$ and $X_{CH_4}^A$ represent the bulk fraction of CH$_4$ in the impactor and in the initial atmosphere, respectively. In all simulations the NH$_3$ reservoir is the one described by the case C, with $X_{NH_3} = 5\%$ and $h_c = 5$ km.

<table><thead><tr><th>No.</th><th>A<sub>s</sub></th><th>M<sub>ini</sub>/<br>M<sub>TA</sub></th><th>X<sub>CH<sub>4</sub></sub><br>(%)</th><th>X<sup>A</sup><sub>CH<sub>4</sub></sub><br>(%)</th><th>M<sub>A</sub>/M<sub>TA</sub> (% of<br>CH<sub>4</sub>)</th><th>M<sub>c</sub>/M<sub>TA</sub> (% of<br>CH<sub>4</sub>)</th></tr></thead><tbody><tr><td>28</td><td>0.2</td><td>1</td><td>1</td><td>5</td><td>0.03 (24%)</td><td>0</td></tr><tr><td>29</td><td>0.2</td><td>1</td><td>0</td><td>0</td><td>0.02</td><td>0</td></tr><tr><td>30</td><td>0.2</td><td>2</td><td>1</td><td>5</td><td>0.04 (39%)</td><td>0.05 (97%)</td></tr><tr><td>31</td><td>0.2</td><td>2</td><td>0</td><td>0</td><td>0.02</td><td>0</td></tr><tr><td>32</td><td>0.2</td><td>5</td><td>1</td><td>5</td><td>0.18 (9%)</td><td>0.17 (81%)</td></tr><tr><td>33</td><td>0.2</td><td>5</td><td>0</td><td>0</td><td>0.77</td><td>1.27</td></tr><tr><td>34</td><td>0.2</td><td>7</td><td>1</td><td>5</td><td>0.76 (5%)</td><td>0.15 (59%)</td></tr><tr><td>35</td><td>0.2</td><td>8</td><td>1</td><td>5</td><td>1.42 (11%)</td><td>0</td></tr><tr><td>36</td><td>0.2</td><td>10</td><td>1</td><td>5</td><td>2.91 (9%)</td><td>0</td></tr><tr><td>37</td><td>0.2</td><td>10</td><td>0</td><td>0</td><td>0.8</td><td>0.25</td></tr><tr><td>38</td><td>0.1</td><td>1</td><td>1</td><td>5</td><td>0.03 (25%)</td><td>0</td></tr><tr><td>39</td><td>0.3</td><td>1</td><td>1</td><td>5</td><td>0.04 (32%)</td><td>0.8 (97%)</td></tr><tr><td>40</td><td>0.3</td><td>5</td><td>1</td><td>5</td><td>0.19 (5%)</td><td>0.22 (74%)</td></tr><tr><td>41</td><td>0.3</td><td>1</td><td>0</td><td>0</td><td>0.02</td><td>0</td></tr><tr><td>42</td><td>0.5</td><td>1</td><td>1</td><td>5</td><td>0.03 (17%)</td><td>0.25 (94%)</td></tr><tr><td>43</td><td>0.5</td><td>5</td><td>1</td><td>5</td><td>0.17 (0.05%)</td><td>0.31 (67%)</td></tr><tr><td>44</td><td>0.5</td><td>1</td><td>0</td><td>0</td><td>0.03</td><td>0</td></tr><tr><td>45</td><td>0.7</td><td>1</td><td>0</td><td>0</td><td>0.06</td><td>0.8</td></tr><tr><td>46</td><td>0.2</td><td>1</td><td>1</td><td>10</td><td>0.04 (62%)</td><td>0</td></tr><tr><td>47</td><td>0.2</td><td>2</td><td>1</td><td>10</td><td>0.04 (35%)</td><td>0.09 (97%)</td></tr><tr><td>48</td><td>0.2</td><td>5</td><td>1</td><td>10</td><td>0.15 (10%)</td><td>0.23 (85%)</td></tr><tr><td>49</td><td>0.2</td><td>10</td><td>1</td><td>10</td><td>2.75 (13%)</td><td>0</td></tr></tbody></table>

Nevertheless even in the most optimistic case (Case A, $X_{NH_3} = 5\%$), a negative balance is obtained for diameters from 6.5 to 27 km.

## 3.2. Simulations with no initial atmosphere

In the following set of simulations, we assume an initially air-less Titan with various NH$_3$ reservoirs and impactors containing 1% of NH$_3$ and either 0% or 1% of CH$_4$. All of the simulations with no initial atmosphere are presented in Table 5. Fig. 5 shows the influence of the albedo and NH$_3$ reservoir on the accumulation of volatiles on Titan's surface. For Case A with $X_{NH_3} = 2\%$ (Fig. 5a and b), the total mass of N$_2$ at the end of the simulations never exceed 0.4 M$_{TA}$ as long as the albedos remains below 0.5. For $A_s = 0.5$, even if the saturation pressure is reached and N$_2$ starts condensing at half of the simulation (Fig. 5a), the accumulation of volatiles at Titan's surface remains limited. The effect of N$_2$ accumulation at the surface becomes important only for $A_s \ge 0.6$, as N$_2$ condensation starts directly at the beginning of the simulations. For this

<figure><img src="image_381.png" /><figcaption>**Fig. 4.** Balance between the mass of N$_2$ degassed from NH$_3$ crustal reservoir by impact ($M_d$) and the atmospheric mass lost by impact-induced erosion ($M_c$), normalized by the impactor mass ($m_{imp}$), as a function of impactor diameters for three different values of (a) surface pressures, (b) impact velocities before entering the atmosphere, (c) initial reservoirs of NH$_3$, $X_{NH_3}$ is the mass fraction of the NH$_3$ contained in the NH$_3$-enriched surface layer and $h_c$ is the thickness of this layer. For comparison, in all three figures, the green curve has the same parameters. (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_382.png" /><figcaption>**Fig. 5.** Evolution of atmospheric pressure (a, c) and $M_{\text{tot}}$, the total mass of N$_2$ in the atmosphere and condensed at the surface (b, d) during the LHB for an initially air-less Titan for cases A (infinite reservoir; a–b) and B ($h_c = 50$ km, c–d) with ammonia fraction $X_{\text{NH}_3}$ of 2%, and three different values of surface albedos ($A_s = 0.2, 0.5$ and $0.7$). Simulations displayed on panels a and b correspond to simulations nos. 3, 6 and 8 in Table 5, on panels c and d to simulations nos. 12, 13 and 14, respectively.</figcaption></figure>

most favorable case where no depletion of the crustal NH$_3$ reservoir is considered, the mass of N$_2$ accumulated at the end of the simulation, mostly in the form of N$_2$ ices, can exceed 1.5 $M_{\text{TA}}$. This result is comparable to the result obtained by Sekine et al. (2011) for a cold Titan, an infinite reservoir of crustal NH$_3$ (2%), and impactor radius of 20 km.

For simulations including a progressive depletion of NH$_3$ in the outer layer owing to successive impacts (Case B), an atmospheric mass comparable to the present one is obtained only if the enriched layer is at least 50 km thick. For $A_s = 0.7$, ~2% of ammonia in this layer is sufficient to generate nearly 1 $M_{\text{TA}}$ (simulation no. 14), while for $A_s = 0.2$, 5% of ammonia is needed (simulations nos. 15 and 19).

For more realistic thicknesses of ammonia-enriched crust ($h_c = 5$ km and $X_{\text{NH}_3} = 5\%$), the final atmospheric pressures are considerably lower (see Figs. 5 and 6), and the accumulated volatile mass barely exceeds 0.04 $M_{\text{TA}}$ (see Table 5, simulation nos. 20, 21, 25, 26 and 27). Adding the contribution of methane brought by the impactors has only a moderate effect on the atmosphere evolution (simulations nos. 22, 23 and 24 in Table 5).

### 3.3. Simulations with an initial atmosphere

In the following set of simulations, we consider preexisting N$_2$-CH$_4$ atmospheres with composition between 0% and 10% of CH$_4$ and masses ranging from 1 to 10 present-day Titan's atmospheric mass ($M_{\text{TA}}$), at the beginning of the LHB. Fig. 7 shows that an atmosphere with an initial mass of 1 $M_{\text{TA}}$ is severely eroded by impacts: after half of the LHB mass has accumulated, about 85% of the initial atmosphere is lost. For a surface albedo of 0.2, and initial temperature of 90 K, all volatiles are in the gaseous phase. The surface liquid condensates when the surface pressure and temperature reach

<figure><img src="image_383.png" /><figcaption>**Fig. 6.** Evolution of the atmospheric pressure at Titan's surface during the LHB, for an initially air-less Titan, the case C and 0.7, 0.5 and 0.2 albedos of surface (simulations nos. 20, 25, 26 from Table 5).</figcaption></figure>

values of 1.08 bar and 87 K, respectively. At its maximum, 24% of the remaining volatiles are condensed at the surface. Most of the condensed liquids is composed of CH$_4$ (67%), while the atmosphere is dominated by N$_2$ (95%). After this maximum, the liquid reservoir progressively evaporates as the atmosphere is continuously eroded by impacts.

For albedos of 0.1, 0.2, 0.3 and 0.5, and similar atmospheric parameters, comparable atmospheric evolutions are obtained (see Fig. 8). The difference between these simulations (nos. 38, 39 and 42 in Table 6) is the mass of volatiles condensed at the surface. For increasing albedos, the condensed volatile reservoir tends to remain longer at the surface of the satellite. For a 0.1 surface albedo, $1.8 \times 10^{20}$ kg of impactors are sufficient to erode the atmosphere and evaporate the surface liquids. For surface albedo of 0.3,
<figure><img src="image_384.png" /><figcaption>**Fig. 7.** Evolution of the atmospheric and condensed volatile mass during the LHB for an initial atmospheric mass equivalent to the present-day Titan's atmospheric mass (simulation no. 28 in Table 6). The vertical line marks the beginning of the condensation of volatiles at the surface.</figcaption></figure>

the surface liquids are maintained during the LHB, while for the simulation with an albedo of 0.5 the surface condensed reservoir keeps increasing all along the simulation.

A strong erosion of the atmosphere during the LHB is observed for all initial atmospheric masses tested in this study. Whatever the initial atmospheric mass, when surface pressure reaches ~1.1 bar, volatiles condense at the surface in the form of a liquid layer (Fig. 9). Liquid reservoirs are mostly composed of CH₄, while the atmospheric compositions are dominated by N₂. For the set of simulations shown in Fig. 9, an atmosphere mass comparable to the present-day mass is preserved only if the initial atmospheric mass is about 7 times the present-day one.

In Fig. 10 we compare the evolution of volatile reservoirs for different initial atmospheric CH₄ fractions. Evolutions of the volatile reservoirs for initial 5% and 10% of CH₄ atmosphere are similar. They only differ in the initiation times of the volatiles condensation. Pure N₂ reservoir experiences less severe erosion. Unlike the previous cases, the liquids at the surface in the N₂ pure case are conserved and the atmospheric surface pressure remains constant all along the simulation. As the atmospheric surface pressure is lower than the atmospheric surface pressure of N₂–CH₄ cases, the erosion is less efficient (see Fig. 4(a)). As a result, a less massive initial volatile reservoir is required to lead to a volatile mass comparable to a present-day one.

Fig. 11 presents an extreme case with a very high albedo ($A_s = 0.7$) and without CH₄, where N₂ remains frozen during the entire LHB. In these conditions, Titan's surface remains in the frozen state during the entire simulation. Due to the low surface pressure, the impact-induced atmospheric erosion is much lower and the atmosphere is globally maintained trough the LHB.

<figure><img src="image_385.png" /><figcaption>**Fig. 8.** Evolution of volatiles masses for 0.1, 0.2, 0.3 and 0.5 surface albedos.</figcaption></figure>

<figure><img src="image_386.png" /><figcaption>**Fig. 9.** Evolution of reservoirs of volatile (atmospheric + condensed masses) during the LHB for several initial atmospheric masses. Solid lines represent the evolution of atmospheric reservoirs, while dotted lines represent the evolution of the liquid surface reservoir.</figcaption></figure>

<figure><img src="image_387.png" /><figcaption>**Fig. 10.** Evolution of the volatile reservoirs for an initial atmosphere ×5 $M_{\text{TA}}$, for three different CH₄ abundances, for a surface albedo of 0.2, and a case C reservoir with $X_{\text{NH}_3}$ = 5% and $h_c$ = 5 km.</figcaption></figure>

<figure><img src="image_388.png" /><figcaption>**Fig. 11.** Evolution of the atmospheric and condensed N₂ during the LHB for a surface albedo of 0.7, for an initial N₂ mass equivalent to a present-day Titan's atmospheric mass.</figcaption></figure>

## 4. Discussion

Following the scenario proposed by Sekine et al. (2011), we tested the possible formation of a N₂-rich atmosphere during the LHB by impact-induced conversion and outgassing of a NH₃ crustal reservoir. To expand the study of Sekine et al. (2011) that considered constant impactor sizes equal to either 20 or 30 km, we considered in our study a distribution of impactor sizes and velocities representative of LHB. To estimate the degassed mass in the impacted crust, we used the scaling laws of Kraus et al. (2011)
derived from hydrocode simulations using a new 5-phase model equation of state for H₂O (Stewart, 2005). Using these scaling laws, we found that an impactor of 20 or 30 km in radius hitting the surface at a velocity of 11 km s⁻¹ degases a crustal volume equal to about five times the impactor mass, while Sekine et al. (2011) predicted a factor of eight for similar conditions. Despite these differences, we obtained results relatively similar to those of Sekine et al. (2011) for initially airless Titan with a high surface albedo (A_s ≥ 0.7) and a constant ammonia fraction (~2%) in the interior (Case A). By taking into account the progressive ammonia depletion in the outer layer owing to successive impacts (Case B), which was not considered in Sekine et al. (2011), we confirmed that, for high surface albedo (A_s ≥ 0.7), a N₂ mass equivalent to the present-day mass may be generated if the outer layer contained about 2% of ammonia over a thickness of at least 50 km. For lower albedos, equivalent generated masses require higher ammonia content: for example, if A_s = 0.2, an ammonia content of 5% is needed. Having a few percent of ammonia over a thickness of 50–100 km might be possible only if Titan was undifferentiated at the time of LHB, which requires a cold accretion. If the outer layer melts during accretion, which is suggested by thermo-accretional modeling (Monteux et al., 2014), most of the ammonia initially incorporated in Titan would migrate in the internal ocean, leaving the outer ice shell depleted in ammonia. Ammonia might be reinjected in a secondary crust formed by cryovolcanic activity (Choukroun and Grasset, 2010), however it seems unlikely that this enriched crust exceeds 5–10 km. In case of a crust enriched in ammonia (Case C), we showed that no atmosphere can be generated and that preexisting atmospheres are severely eroded by impacts during the LHB, independently of the initial atmospheric masses or methane fractions. Our simulations showed that the impact-induced erosion appears to be a controlling process in the evolution of Titan's atmosphere. Our results rely on the parameterization of Shuvalov (2009), dedicated to Earth conditions and where the impactor diameters range between 1 and 30 km. As suggested by the comparison with the simulations of Korycansky and Zahnle (2011) dedicated to Titan's atmosphere (see Fig. 2), it is possible that the erosion efficiency predicted by Shuvalov's scaling laws is overestimated for larger impactors. To test the sensitivity of our model to the erosion parameterization, we arbitrarily divided the erosion efficiency by 2, 3 and 5, independently of impactor diameters. These tests showed that the erosion can be counterbalanced if the erosion efficiency is divided by three and if the most favorable degassing scenario (case C, X_{NH_3} = 5%, h_c = 5 km) is considered. These tests indicate that our results may change only if there is a significant overestimation of the erosion efficiency by Shuvalov (2009). A lower impactor flux may also limit the atmospheric erosion. It was suggested by Nimmo and Korycansky (2012) that the total impactor mass delivered during the LHB could be ~10 times smaller than the value used here, estimated from the Nice model (Charnoz et al., 2009). Obviously in these conditions, the erosion will be strongly reduced for preexisting massive atmosphere, reducing the requirement on the mass for the primitive atmosphere. However, it would still be impossible to generate a massive atmosphere.

Impact-induced erosion of volatiles condensed at the surface may be another source of uncertainties. Indeed, in our approach we do not consider this process. For the simulation where the initial atmosphere is seven times more massive than the present day atmosphere, the fraction of volatiles in the liquid phase never exceeds 16%. Therefore, its influence on the surface erosion remains negligible. Surface liquid fractions may play a non negligible role in simulations with moderate initial atmospheric masses (1–5 times the present-day Titan's atmospheric mass, see Figs. 7 and 9). In these cases, we do not preserve a massive atmosphere

at the end of the simulations. Including the surface erosion will only enhance the volatile loss, therefore it should not change our conclusions. In cases of a high surface albedo (>0.7), most of the atmospheric mass condenses at the surface and is preserved through the LHB, because low atmospheric pressure leads to a lower impact-induced atmospheric erosion. These are the only cases where the surface erosion could have a non-negligible influence on the volatile mass at Titan's surface and this result has to be examined with caution. Indeed, the study of Nimmo and Korycansky (2012) showed that the surface erosion by impact may eject a non negligible mass of surface material during the LHB. Therefore, even in the favorable case with a high surface albedo, a significant fraction of condensed N₂ may be still lost by successive impacts.

Another critical aspect in the evolution of the atmospheric balance is the volatile mass supplied to the atmosphere during the impact. In our study we consider only the methane supplied by the impactor, no crustal CH₄ degassing is estimated. Nevertheless, methane could be present in a significant amount in the form of clathrates hydrates (Osegovic and Max, 2005; Tobie et al., 2006, 2012). The pressures and temperatures reached during the high-velocity impacts are high enough to destabilize the methane clathrate and, hence, to release CH₄ to Titan's atmosphere. Our simulations show that Titan's surface layers are significantly affected by the impacts at the end of the LHB. An important fraction of methane that might be stored in the crust should be released during the LHB. The non-degassed terrains represent 30% of Titan's surface in the case of an initial atmosphere ten times more massive than the present-day Titan's atmosphere and only 0.02% for an initially air-less body (see Fig. 3c). Even though the destabilization of CH₄ clathrate reservoir is different from crustal degassing implemented here, we expect that a large part of CH₄ should be degassed by impact during the LHB. If a clathrate-rich crust was already formed at the time of the LHB, most of this crustal reservoir would be destabilized leading to accumulation of several hundred of meters of liquid methane at the surface. Other gas compounds present in comets or in the crust in similar or higher proportions than NH₃ and CH₄, in particular CO or CO₂, could also affect the atmospheric balance. The low CO abundance in today's atmosphere (Baines et al., 2006; de Kok et al., 2007) suggests that either a small amount of CO was brought to Titan, which seems in contradiction with cometary composition, or that a very efficient process (still to be determined) destroyed most of the primordial CO. CO₂ brought by the LHB impactors, degassed from the crust or possibly converted from CO during impacts (Ishimaru et al., 2011), will remain in the condensed phase and should therefore not affect significantly the radiative balance. However, further studies are needed to evaluate the influence of these chemical species and their fate during the LHB era.

As shown in Section 3.1, the atmospheric mass balance is also sensitive to the impactor diameter, and therefore our results depend on the assumed size distribution. Here, we have used the size distribution of Charnoz et al. (2009) relying on the Nice model and scaled to lapetus' crater distribution. Considering the 1.5 bar atmosphere case from Fig. 4 a. and using this size distribution, only 0.03% of the total mass is brought by impactors smaller than 6 km, that is the critical impactor diameter above which the atmospheric mass balance become negative. However, this fraction should be at least 50% in order to limit erosion and generate an atmosphere. Increasing this tendency would imply a significant change in the assumed size distribution. Secondary debris resulting from a giant impact in the system of Saturn (e.g. Levison et al., 2011; Asphaug and Reufer, 2013) could possibly provide a source of small impactors (Nakamura and Fujiwara, 1991). As suggested by Farinella et al. (1997) and Dobrovolskis and Lissauer (2004), Hyperion may have resulted from such giant impacts. The total mass of debris
resulting from proto-Hyperion disruption and possibly re-accreted by Titan is evaluated ranging between ~10¹⁹ and several times 10²⁰ kg, which is comparable to the estimated LHB cumulative mass. These debris have a significantly lower velocity (<4 km s⁻¹) (Farinella et al., 1997) resulting in less atmospheric erosion but also limited degassing as the impact pressure barely exceeds 8 GPa (obtained for a velocity of 3.4 km s⁻¹). Large impactors might also be more numerous than assumed here, and have more complex effects than what we considered in the present model. As shown by Zahnle et al. (2014), impactors larger than $d_i > 13$ km with an impact velocity of 11 km s⁻¹ may completely break the ice crust, resulting in a major water flooding event and strong degassing from volatiles contained in the ocean (Tobie et al., 2012). Such catastrophic events would strongly affect the atmospheric balance, possibly leading to the formation of the atmosphere. This suggests that a few large impact events may be more favorable to the generation of an atmosphere than a continuous flux of kilometric-size impactors. A full approach, solving consistently the impact dynamics and the coupling between the subsurface ocean, ice shell and atmosphere, is required to investigate the consequences of such large impacts on Titan's atmospheric balance.

## 5. Conclusion

We investigate the supply and loss of $N_2$ and $CH_4$ by impacts on Titan. We show that an atmosphere with a mass equivalent to the present-day is unlikely to form during the LHB era. We investigate the sensitivity of the emergence of Titan's atmosphere on the $NH_3$ crustal reservoir and on the surface conditions. We show that a $N_2$-rich atmosphere can be generated during the LHB only if Titan contained a few percents of ammonia over at least 50 km beneath its surface, implying an undifferentiated interior. We also show that preexisting $N_2$ atmospheric mass equivalent to the present-day atmosphere mass could survive to the LHB only for very high surface albedos resulting in very low atmospheric pressures and if the impact-induced loss of $N_2$ condensed at the surface is limited. Further modeling efforts are needed to understand the coupled evolution of crustal, surface and atmospheric reservoir of N-bearing species and the consequence for the emergence of atmosphere on early Titan.

Except for high surface albedos ($A_s \ge 0.7$) where most of the released volatiles are frozen at the surface, the high-velocity impacts during the LHB led to a strong atmospheric erosion. The atmospheric loss by impact is controlled by the equilibrium atmospheric pressure, which depends mostly on the surface albedo and the fraction of methane. For albedo of 0.2 and initial methane fraction between 0% and 5% in the atmosphere the pre-LHB atmosphere should be between ~5 and 7 times more massive than at present to sustain an atmosphere equivalent to the present-day one.

If Titan was differentiated at the time of the LHB and had a surface albedo lower than 0.6, our results imply that either a massive atmosphere was formed on Titan during its accretion or that the nitrogen-rich atmosphere was generated after the LHB. The first hypothesis involves a warm accretion with an efficient conversion of $NH_3$ into $N_2$, due to either photochemistry (Atreya et al., 1978) or impact-induced chemistry (McKay et al., 1988; Ishimaru et al., 2011). This is likely that a massive atmosphere interacting with a water ocean formed at the end of accretion due to cumulative impact heating (Monteux et al., 2014). However it is still unclear if the conditions to produce several bars of $N_2$ from $NH_3$ or other N-bearing molecules was met. The second hypothesis implies an efficient outgassing process and an internal conversion of N-bearing compounds into $N_2$, as suggested by (Glein et al.,

2009; Tobie et al., 2012). Here again, though outgassing is possible, the internal production of $N_2$ remains unconstrained. Further studies will be needed to constrain these two possible scenarios.

## Acknowledgments

We thank two anonymous reviewers, whose comments improved the manuscript. The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 259285). J. Monteux was funded by Agence Nationale de Recherche (Accretis decision ANR-10-PDOC-001-01). B. Charnay acknowledges support from an appointment to the NASA Postdoctoral Program at NAI Virtual Planetary Laboratory, administered by Oak Ridge Affiliated Universities.

## References

* Alibert, Y., Mousis, O., 2007. Formation of Titan in Saturn's subnebula: Constraints from Huygens probe measurements. Astron. Astrophys. 1060, 1051–1060.

* Asphaug, E., Reifer, A., 2013. Late origin of the Saturn system. Icarus 223, 544–565.

* Atreya, S.K., Donahue, T., Kuhn, W., 1978. Evolution of a nitrogen atmosphere on Titan. Science 201, 611–613.

* Atreya, S.K., Lorenz, R.D., Waite, J.H., 2009. Volatile origin and cycles: Nitrogen and methane. In: Titan from Cassini-Huygens, pp. 177–199.

* Baines, K. et al., 2006. On the discovery of CO nighttime emissions on Titan by Cassini/VIMS: Derived stratospheric abundances and geological implications. Planet. Space Sci. 109, 1552–1562.

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

* Barr, A.C., Citron, R.I., Canup, R.M., 2010. Origin of a partially differentiated Titan. Icarus 209, 858–862.

* Bockelée-Morvan, D., Crovisier, J., 2004. The composition of cometary volatiles. In: Comets II, pp. 391–424.

* Bruesch, L.S., Asphaug, E.I., 2004. Modeling global impact effects on middle-sized icy bodies: Applications to Saturn's moons. Icarus 168, 457–466.

* Charnay, B. et al., 2014. Titans past and future: 3D modeling of a pure nitrogen atmosphere and geological implications. Icarus 241, 269–279.

* Charnoz, S. et al., 2009. Did Saturn's rings form during the Late Heavy Bombardment? Icarus 199, 413–428.

* Choukroun, M., Grasset, O., 2010. Thermodynamic data and modeling of the water and ammonia-water phase diagrams up to 2.2 GPa for planetary geophysics. J. Chem. Phys. 133, 144502.

* Croft, S., 1982. A first-order estimate of shock heating and vaporization in oceanic impacts. Geol. Soc. Am. 190, 143–152.

* Davison, T., Collins, G., Ciesla, F., 2010. Numerical modelling of heating in porous planetesimal collisions. Icarus 208, 468–481.

* de Kok, R. et al., 2007. Oxygen compounds in Titan's stratosphere as observed by Cassini CIRS. Icarus 186, 354–363.

* de Niem, D. et al., 2012. Atmospheric erosion and replenishment induced by impacts upon the Earth and Mars during a heavy bombardment. Icarus 221, 495–507.

* Dobrovolkskis, A.R., Lissauer, J.J., 2004. The fate of ejecta from Hyperion. Icarus 169, 462–473.

* Farinella, P., Marzari, F., Matteoli, S., 1997. The disruption of Hyperion and the origin of Titan's atmosphere. Astron. J. 113, 2312–2316.

* Fortes, A. et al., 2007. Ammonium sulfate on Titan: Possible origin and role in cryovolcanism. Icarus 188, 139–153.

* Genda, H., Abe, Y., 2005. Enhanced atmospheric loss on protoplanets at the giant impact phase in the presence of oceans. Nature 433, 842–844.

* Glein, C.R., Desch, S.J., Shock, E.L., 2009. The absence of endogenic methane on Titan and its implications for the origin of atmospheric nitrogen. Icarus 204, 637–644.

* Gomes, R. et al., 2005. Origin of the cataclysmic Late Heavy Bombardment period of the terrestrial planets. Nature 435, 466–469.

* Grasset, O., Sotin, C., Deschamps, F., 2000. On the internal structure and dynamics of Titan. Planet. Space Sci. 48, 617–636.

* Griffith, A., Zahnle, K.J., 1995. Influx of cometary volatiles to planetary moons: The atmospheres of 1000 possible Titans. J. Geophys. Res. 100, 16907–16922.

* Griffith, C. et al., 2013. Titan's evolving climate. In: Comparative Climatology of Terrestrial Planets, pp. 1–27.

* Hersant, F. et al., 2008. Interpretation of the carbon abundance in Saturn measured by Cassini. Planet. Space Sci. 56, 1103–1111.

* Ishimaru, R. et al., 2011. Oxidizing proto-atmosphere on Titan: Constraint from $N_2$ formation by impact shock. Astrophys. J. 741, L10.

* Kargel, J.S., 1991. Brine volcanism and the interior structures of asteroids and icy satellites. Icarus 94, 368–390.

* Korycansky, D.G., Zahnle, K.J., 2005. Modeling crater populations on Venus and Titan. Planet. Space Sci. 53, 695–710.

* Korycansky, D.G., Zahnle, K.J., 2011. Titan impacts and escape. Icarus 211, 707–721.
Kraus, R.G., Senft, L.E., Stewart, S.T., 2011. Impacts onto H2O ice: Scaling laws for melting, vaporization, excavation, and final crater size. *Icarus* 214, 724–738.

Levison, H. et al., 2011. Ridge formation and de-spinning of Iapetus via an impact-generated satellite. *Icarus* 214, 773–778.

Linstrom, P.J., Mallard, W.G. (Eds.), 2015. NIST Chemistry WebBook, NIST Standard Reference Database Number 69. National Institute of Standards and Technology, Gaithersburg MD, p. 20899, http://webbook.nist.gov.

Lorenz, R.D., Lunine, J.I., McKay, C.P., 1997. Titan under a red giant Sun: A new kind of habitable Moon. *Geophys. Res. Lett.* 24, 2905–2908.

Lorenz, R.D., McKay, C.P., Lunine, J.I., 1999. Analytic investigation of climate stability on Titan: Sensitivity to volatile inventory. *Planet. Space Sci.* 47, 1503–1515.

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

Mandt, K.E. et al., 2012. The 12C/ 13C Ratio on Titan from Cassini INMS measurements and implications for the evolution of methane. *Astrophys. J.*, 749

Mandt, K.E. et al., 2014. Protosolar ammonia as the unique source of Titan's nitrogen. *Astrophys. J.* 788, L24.

McKay, C.P., Scattergood, T., Pollack, J.B., 1988. High-temperature shock formation of N₂ and organics on primordial Titan. *Nature* 332, 520–522.

McKay, C.P., Lorenz, R., Lunine, J., 1999. Analytic solutions for the antigreenhouse effect: Titan and the early Earth. *Icarus* 61, 56–61.

Melosh, H., Vickery, A., 1989. Impact erosion of the primordial atmosphere of Mars. *Nature* 338, 487–489.

Monteux, J. et al., 2014. Can large icy moons accrete undifferentiated? *Icarus* 237, 377–387.

Morbidelli, A. et al., 2005. Chaotic capture of Jupiter's Trojan asteroids in the early Solar System. *Nature* 435, 462–465.

Mousis, O., Gautier, D., Bockelée-Morvan, D., 2002. An evolutionary turbulent model of Saturn's subnebula: Implications for the origin of the atmosphere of Titan. *Icarus* 156, 162–175.

Nakamura, A., Fujiwara, A., 1991. Velocity distribution of fragments formed in a simulated collisional disruption. *Icarus* 92, 132–146.

Niemann, H.B. et al., 2010. Composition of Titan's lower atmosphere and simple surface volatiles as measured by the Cassini–Huygens probe gas chromatograph mass spectrometer experiment. *J. Geophys. Res.* 115, E12006.

Nimmo, F., Korycansky, D.G., 2012. Impact-driven ice loss in outer Solar System satellites: Consequences for the Late Heavy Bombardment. *Icarus* 219, 508–510.

Osegovic, J.P., Max, M.D., 2005. Compound clathrate hydrate on Titan's surface. *J. Geophys. Res.* 110, E08004.

Owen, T.C., 1982. The composition and origin of Titan's atmosphere. *Planet. Space Sci.* 30, 833–838.

Pham, L., Karatekin, O., Dehant, V., 2011. Effects of impacts on the atmospheric evolution: Comparison between Mars, Earth, and Venus. *Planet. Space Sci.* 59, 1087–1092.

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

Pierazzo, E., Melosh, H., 2000. Understanding oblique impacts from experiments, observations, and modeling. *Annu. Rev. Earth Planet. Sci.* 28, 141–167.

Rousselot, P. et al., 2014. Toward a unique nitrogen isotopic ratio in cometary ices. *Astrophys. J.* 780, L17.

Sekine, Y. et al., 2011. Replacement and late formation of atmospheric N₂ on undifferentiated Titan by impacts. *Nat. Geosci.* 4, 359–362.

Senshu, H., 2002. Thermal evolution of a growing Mars. *J. Geophys. Res.* 107, 5118.

Shinnaka, Y. et al., 2014. ¹⁴NH₂/¹⁵NH₂ ratio in comet C/2012 S1 (ISON) observed during its Outburst in 2013 November. *Astrophys. J.* 782, L16.

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

Shuvalov, V., 2010. Atmospheric erosion induced by oblique impacts. In: *Lunar and Planetary Science Conference*.

Stewart, S.T., 2005. Shock properties of H2O ice. *J. Geophys. Res.* 110, E03005.

Svetsov, V.V., 2000. On the efficiency of the impact mechanism of atmospheric erosion. *Sol. Syst. Res.* 34, 441–453.

Svetsov, V.V., 2007. Atmospheric erosion and replenishment induced by impacts of cosmic bodies upon the Earth and Mars. *Sol. Syst. Res.* 41, 28–41.

Thompson, W.R., 1985. In: Rolfe, E., Battrick, B. (Eds.), Phase Equilibria in N₂-Hydrocarbon Systems: Applications to Titan. *ESA Special Publication*, pp. 109–119.

Tobie, G., Lunine, J.I., Sotin, C., 2006. Episodic outgassing as the origin of atmospheric methane on Titan. *Nature* 440, 61–64.

Tobie, G., Gautier, D., Hersant, F., 2012. Titan's bulk composition constrained by Cassini–Huygens: Implication for internal outgassing. *Astrophys. J.* 752, 125.

Tsiganis, K. et al., 2005. Origin of the orbital architecture of the giant planets of the Solar System. *Nature* 435, 459–461.

Waite, J.H. et al., 2009. Liquid water on Enceladus from observations of ammonia and 40Ar in the plume. *Nature* 460, 487–490.

Zahnle, K.J., 2001. Differential cratering of synchronously rotating satellites by ecliptic comets. *Icarus* 153, 111–129.

Zahnle, K.J. et al., 2003. Cratering rates in the outer Solar System. *Icarus* 163, 263–289.

Zahnle, K.J., Korycansky, D.G., Nixon, C.A., 2014. Transient climate effects of large impacts on Titan. *Icarus* 229, 378–391.
<figure><img src="image_389.png" /></figure>

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

# Consequences of large impacts on Enceladus' core shape

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

CrossMark

J. Monteux<sup>a,b,*</sup>, G.S. Collins<sup>c</sup>, G. Tobie<sup>b</sup>, G. Choblet<sup>b</sup>

## ARTICLE INFO

**Article history:**

Received 13 May 2015

Revised 31 August 2015

Accepted 22 September 2015

Available online 8 October 2015

**Keywords:**

Enceladus

Impact processes

Cratering

Interiors

Accretion

## ABSTRACT

The intense activity on Enceladus suggests a differentiated interior consisting of a rocky core, an internal ocean and an icy mantle. However, topography and gravity data suggests large heterogeneity in the interior, possibly including significant core topography. In the present study, we investigated the consequences of collisions with large impactors on the core shape. We performed impact simulations using the code iSALE2D considering large differentiated impactors with radius ranging between 25 and 100 km and impact velocities ranging between 0.24 and 2.4 km/s. Our simulations showed that the main controlling parameters for the post-impact shape of Enceladus' rock core are the impactor radius and velocity and to a lesser extent the presence of an internal water ocean and the porosity and strength of the rock core. For low energy impacts, the impactors do not pass completely through the icy mantle. Subsequent sinking and spreading of the impactor rock core lead to a positive core topographic anomaly. For moderately energetic impacts, the impactors completely penetrate through the icy mantle, inducing a negative core topography surrounded by a positive anomaly of smaller amplitude. The depth and lateral extent of the excavated area is mostly determined by the impactor radius and velocity. For highly energetic impacts, the rocky core is strongly deformed, and the full body is likely to be disrupted. Explaining the long-wavelength irregular shape of Enceladus' core by impacts would imply multiple low velocity (<2.4 km/s) collisions with deca-kilometric differentiated impactors, which is possible only after the LHB period.

© 2015 Elsevier Inc. All rights reserved.

## 1. Introduction

Despite its small size (R = 252 km), Saturn's moon Enceladus is one of the most geologically active body of the Solar System. Its surprising endogenic activity is characterized by a very active province at the South Pole, from which eruptions of water vapor and ice grains emanating from warm tectonic ridges have been observed by the Cassini spacecraft (Porco et al., 2006; Hansen et al., 2006; Waite et al., 2006; Spencer et al., 2006). This activity is associated with a huge heat power estimated between 5 and 15 GW from thermal emission (Spencer and Nimmo, 2013), which implies a warm interior, consistent with a liquid water layer underneath the ice shell and a differentiated interior (Nimmo et al., 2007; Schubert et al., 2007). Models of tidal dissipation may explain why the activity is concentrated at the poles, where dissipation is predicted to be maximal (Tobie et al., 2008; Běhounková et al., 2010). However, there is still no satisfactory

explanation for why this activity is located only in the south, and not in the north.

Based on the global shape data which show a depression at the south pole (Thomas et al., 2007), it has been proposed that the ocean may be located only in the southern hemisphere (Collins and Goodman, 2007), thus explaining why the activity would be concentrated at the south (Tobie et al., 2008). Gravity and shape data indicate that such an ocean would be at depths of about 30–40 km and extend up to south latitudes of about 50° (less et al., 2014). It has been proposed that the dichotomy between the north and south hemispheres may be the result of asymmetry in core shape (McKinnon, 2013). Due to the low pressure and moderate temperature expected in Enceladus' core, large topography anomalies may indeed be retained on very long periods of time (McKinnon, 2013) and may explain why convection-driven activities in the ice shell is confined only to the south polar terrain (Showman et al., 2013). Besides the south polar depression, core topography anomalies could explain, at least partly, the existence of other big depressions observed at moderate latitudes (between 15°S and 50°N) and uncorrelated with any geological boundaries (Schenk and McKinnon, 2009).

<footnote><sup>a</sup> Laboratoire Magmas et Volcans, Université Blaise Pascal, CNRS, IRD, Clermont-Ferrand, France</footnote>

<footnote><sup>b</sup> Laboratoire de Planétologie et de Géodynamique, LPG Nantes, CNRS UMR6112, Université de Nantes, France</footnote>

<footnote><sup>c</sup> Impacts and Astromaterials Research Centre, Department of Earth Science and Engineering, Imperial College London, United Kingdom</footnote>

<footnote>* Corresponding author at: Laboratoire Magmas et Volcans, Université Blaise Pascal, CNRS, IRD, Clermont-Ferrand, France.</footnote>
McKinnon (2013) proposed three hypotheses to explain the possible irregularity of Enceladus' rocky core: accretional melting of the outer region of the icy moon associated with a degree-one instability; accretion of icy protomoons around irregular rock chunks; and collisional merger of two previously differentiated protomoons. Here we test the latter hypothesis by investigating the consequences of the collision of a large differentiated impactor on the shape of Enceladus' core. Collisions with large differentiated bodies were likely at the end of satellite accretion, during the final assemblage phase (e.g. Asphaug and Reufer, 2013). Large impact basins on other saturnian moons (e.g. Iapetus (Giese et al., 2008), Mimas (Schenk, 2011), Titan (Neish and Lorenz, 2012)) and other Solar System bodies (e.g. Vesta (Schenk et al., 2012)) could represent remnant evidences of such collisions. Large impacts occurring at the end of the accretion and after, during the rest of the satellite's evolution, likely influenced the internal structure and especially the shape of its rocky core. It is also important to determine the conditions under which Enceladus would have survived disruption by collisions with deca-kilometric objects, which would place constraints on its accretion and the subsequent impact history.

To constrain the consequences of large-scale impacts on Enceladus, we simulated head-on collisions of differentiated impactors with diameter ranging between 50 and 200 km using the iSALE2D shock physics code (Wünnemann et al., 2006; Collins et al., 2004; Davison et al., 2010). From these simulations, we tracked the evolution of rock fragments coming from the impactor and the impact-induced modification of Enceladus's core shape. In particular, we quantified the sensitivity in these outcomes to key model parameters, such as impactor velocity and radius, as well as structure and mechanical properties of Enceladus' interior (porosity, strength, temperature profile, core size, presence of an internal ocean). In Section 2, we describe our numerical modelling approach; in Section 3 we present our results. We discuss our results in the context of the presence of a water ocean in Section 3.3. Conclusions are highlighted in Section 4.

## 2. Impact modeling

To model the thermo-mechanical evolution of material during an impact between two differentiated icy bodies, we use iSALE2D (Wünnemann et al., 2006; Collins et al., 2004). This numerical model is a multi-rheology, multi-material shock physics code based on the SALE hydrocode (Amsden et al., 1980) that has been extended and modified specifically to model planetary-scale impact crater formation (e.g., Amsden et al., 1980; Melosh et al., 1992; Ivanov et al., 1997; Collins et al., 2004; Wünnemann et al., 2006; Davison et al., 2010). In our simulations, the target structure and the impactor were simplified to two- or three-layer spherical bodies consisting of a rocky core, an icy mantle and for the three-layer case an internal ocean. Interpretation of gravity data collected by the Cassini spacecraft indicates that the core density could be as low as 2400 kg m⁻³, corresponding to a core radius of about 200 km (less et al., 2014). However, as Enceladus appears to be relatively far from hydrostatic equilibrium (less et al., 2014), there are still significant uncertainties on the core radius and density. The low core density inferred from gravity data suggests that the rocky core might be significantly porous, with pores filled by water ice and/or liquid water, and that a significant fraction of the core may consist of hydrated silicate minerals. Currently, iSALE2D does not have provision to describe the behavior of an ice-rock or water-rock mixture. In our simulations, for simplicity, we assume complete segregation of the rock and ice-water phase into discrete layers and we consider dunite as representative of the rock phase (with density ρs = 3330 kg m⁻³). We reduce the

density of the core by including some initial porosity φ (defined as the ratio of pore volume to total volume) within it, varying from 0% to 50%, corresponding to radius varying between typically 160 km and 200 km. Assuming a core made of pure dunite, a radius as large as 200 km is consistent with a core porosity of about 50%, which is at the upper end of the estimated porosity in large asteroids (Lindsay et al., 2015). A significant fraction of the core may also consist of hydrated minerals such as serpentine. In this case a 200 km core radius would imply a lower porosity. For simplicity, we consider only dunite as core materials and vary the porosity up to values of 50%. We also test the possible effect of porosity in the ice shell by considering values up to 20% as suggested by Besserer et al. (2013).

In our models, we consider the extreme case where the pores of both ice and rocks consist of voids, and are not filled with secondary materials (i.e. water or ice in rock pores). The difference between saturated porosity (with ice or liquid water) and voids may lead to differences in terms of mechanical and thermal properties. This aspect will be discussed in the last section. The effect of both rock and ice porosity is treated using the ε-α porosity compaction model (Wünnemann et al., 2006; Collins et al., 2013), which accounts for the collapse of pore space by assuming that the compaction function depends upon volumetric strain. For sake of simplicity, we assume that the impactor material has an identical composition and porosity to those of the target.

The impact velocity $v_{imp}$ can be decomposed into two contributions:

$$v_{imp} = \sqrt{v_{esc}^2 + v_{\infty}^2} \quad (1)$$

where $v_{esc}$ is the escape velocity of the impacted planet and $v_{\infty}$ is the velocity of the impactor at a distance much greater than that over which the gravitational attraction of the impacted planet is important. The escape velocity of Enceladus is $v_{esc} = 240$ m/s. As we consider collisions with relatively large objects ($R_{imp} = 25–100$ km), we limit our analysis to moderate relative velocities, varying between $v_{esc}$ and $10 \times v_{esc}$, in order to limit the impact-induced deformation of the satellite and avoid full disruption (Benz and Asphaug, 1999; Asphaug, 2010). Moreover, this low-velocity impact regime is representative of the collisional environment at the end of the accretion. Indeed, N-body simulations from Dwyer et al. (2013) show that random impact velocity of proto-satellites mostly ranges between $v_{esc}$ and $5v_{esc}$.

We approximated the thermodynamic response of the icy material using the Tillotson EoS for Ice as in Bray et al. (2008) and of the rocky material using the ANEOS EoS for dunite material as in Benz et al. (1989) and Davison et al. (2010) (see Table 1 for parameter values). Standard strength parameters for dunite were used to form the static strength model for the rocky core (Collins et al., 2004; Davison et al., 2010). The static strength model for ice used in iSALE was derived from low temperature, high pressure laboratory data and accounts for the material strength dependence on pressure, damage and thermal softening (Bray et al., 2008). We also explored the effect on our results of the cohesion of the damaged material (referred to here as $Y_i$ for ice and $Y_s$ for rocks), which represents the minimum zero-pressure shear strength of cold material (strength is reduced to zero at the melt temperature). The minimum strength values considered in our models range between 10–500 kPa for ice and 100–10⁴ kPa for silicate material. The Tillotson EoS for ice is severely limited in its applicability for hypervelocity impact; it includes no solid state or liquid phase changes. However, as we limit here our analysis to low velocity encounters (240 < $v_{imp}$ < 2400 m s⁻¹), thought to be dominant at the end of the accretion, as shown in our simulations, no significant ice melting occurs and the use of Tillotson EoS is a reasonable approximation. We also used the Tillotson EoS for the liquid water.
<table><caption>**Table 1**

Typical parameter values for numerical models.</caption><tbody><tr><td>Enceladus radius</td><td><i>R</i></td><td>250 km</td></tr><tr><td>Rocky core radius</td><td><i>R</i><sub>core</sub></td><td>160–200 km</td></tr><tr><td>Icy mantle thickness</td><td><i>δ</i><sub>m</sub></td><td>50–90 km</td></tr><tr><td>Surface gravity field</td><td><i>g</i><sub>0</sub></td><td>0.113 m s<sup>-2</sup></td></tr><tr><td>Escape velocity</td><td><i>v</i><sub>esc</sub></td><td>240 m/s</td></tr><tr><td>Impactor radius</td><td><i>R</i><sub>imp</sub></td><td>25–100 km</td></tr><tr><td>Impact velocity</td><td><i>v</i><sub>imp</sub></td><td>240–2400 m/s</td></tr><tr><td colspan="3"><em>Mantle properties (Ice)</em></td></tr><tr><td>Initial density</td><td><i>ρ</i><sub>i</sub></td><td>820 kg m<sup>-3</sup></td></tr><tr><td>Equation of state type</td><td></td><td>Tillotson</td></tr><tr><td>Poisson</td><td></td><td>0.33</td></tr><tr><td>Porosity</td><td></td><td>0–20%</td></tr><tr><td>Minimum strength</td><td><i>Y</i><sub>i</sub></td><td>10–500 kPa</td></tr><tr><td colspan="3"><em>Core properties (Dunite)</em></td></tr><tr><td>Rocky core density</td><td><i>ρ</i><sub>s</sub></td><td>3330 kg m<sup>-3</sup></td></tr><tr><td>Equation of state type</td><td></td><td>ANEOS</td></tr><tr><td>Poisson</td><td></td><td>0.25</td></tr><tr><td>Porosity</td><td></td><td>0–50%</td></tr><tr><td>Minimum strength</td><td><i>Y</i><sub>s</sub></td><td>100 kPa–10 MPa</td></tr></tbody></table>

Material weakening during impact may also be achieved by acoustic fluidization and/or thermal softening (Melosh and Ivanov, 1999), the latter of which is especially efficient for large-scale events (Potter et al., 2012). Our simulations including acoustic fluidization that assumed typical block-model parameters favored in other works showed no significant effect on simulation results (see also discussion section). Hence, for simplicity and to reduce the number of free parameters, we chose to neglect acoustic fluidization. We do, however, include the effect of temperature on shear strength using the temperature-strength relationship proposed by Ohnaka (1995) and described by Collins et al. (2004) and we set the thermal softening coefficient in this expression to 1.2 as suggested by Bray et al. (2008). Since we consider the thermal softening during the impact, the thermal structure of Enceladus before the impact is probably a key parameter governing the post-impact state. However, the early temperature profile for such a small body is poorly constrained. Accretionary models seem to favor a cold accretion with inner temperatures close to the equilibrium temperature (Schubert et al., 1981; Monteux et al., 2014). To test the influence of the initial thermal conditions, we consider three different pre-impact temperature profiles for the impacted moon: constant temperature, conductive profile, two-layered advective profile. The impactor is assumed to have a constant temperature with $T = 100$ K.

Owing to the axisymmetric geometry of iSALE2D, we consider only head-on collisions (impact angle of 90° to the target tangent plane). The role of impact angle is left to future studies. To limit computation time, a 1-to-2 km spatial resolution is used, which is sufficient to describe the deflection of the rock core surface. The gravity is calculated from the density structure. For the largest and fastest impacts, we use iSALE2D's self-gravity gravity model (Collins et al., 2011) to correctly assess the gravity field as the body is strongly deformed and the center of mass of the target moves upon the collision. As this self-gravity model is expensive in terms of computational time, we limit our post impact monitoring to the time needed to deform the rocky core (i.e. we consider that the fall-back of icy material and the icy-mantle slumping has only a very minor effect on the morphology of the rocky core). For all the impacts characterized here, this corresponds to the first hour after the impact.

# 3. Numerical results

## 3.1. Non-porous models

Fig. 1 shows three characteristic simulations: ($v_{imp} = 10v_{esc}$, $R_{imp} = 25$ km), ($v_{imp} = 10v_{esc}$, $R_{imp} = 75$ km) and ($v_{imp} = v_{esc}$,

$R_{imp} = 75$ km). After such events, a large volume of Enceladus' mantle is displaced or escapes the orbit of the icy moon. To get a quantitative measure of deformation induced by the impact event, we monitor the plastic strain experienced by the impacted material. In particular, we calculate the total plastic strain which is the accumulated sum of plastic shear deformation, regardless of the sense of shear (Collins et al., 2004). As represented in Fig. 1, the icy material is highly disturbed by the impact and most of the plastic deformation occurs in this layer. For the largest impact velocities (Fig. 1, left and middle), deformation also occurs at the top of the rocky core and leads to the formation of a depression. The material removed from the depression is displaced in a very small uplift of the core, surrounding the depression.

For small impact velocities (Fig. 1, right), the icy mantle is also highly deformed but the impactor's rocky core is trapped within the ice layer. In this low-velocity case, the deformation of the target's core and the impact melt production are minor but the surrounding ice is warmed up. Hence, over a longer time scale governed by a Stokes' flow, the impactor's core gently spreads over the surface of the pre-existing rocky core favoring the formation of successive fragmented silicate layers (Roberts, 2015). Depending on the impactor size and impact velocities, our simulations show that core merging occurs into three distinct regimes (Fig. 2):

(1) For small impactors and impact velocities close to ~$v_{esc}$, the impactor's core is simply buried within Enceladus' icy mantle at a depth that scales with the penetration depth $p$ (Orphal et al., 1980; Murr et al., 1998):

$$p/R_{imp} = A v_{imp}^{2/3} \qquad (2)$$

where $A$ is a function of the bulk sound velocity, the geometry and density difference between the impactor and the target.

(2) For higher impact velocities or larger impactors, the kinetic energy increases and hence penetration of the impactor's core through the target ice mantle is facilitated. When the impactor penetration depth, $p$ (2), exceeds the ice-mantle thickness, $\delta_m$, the impactor induces a deflection of the core boundary (Fig. 2), the amplitude of which depends on the impactor energy remaining after crossing the ice mantle. For $p \sim \delta_m$ or slightly larger, the impactor core spreads above the target's core (leading to a positive core-topography anomaly defined as the difference of post- and pre-impact core radii below the impact site). (3) However, if more energy is available, $p > \delta_m$ and the core is strongly deformed, possibly leading to severe deformation of the satellite, as illustrated in Fig. 2 for impactors larger than 75 km and/or impact velocities $\ge 10v_{esc}$. It has to be noted that, as we limit our post impact monitoring to one hour, for the most energetic impact cases with large impact velocities ($\ge 6$ km/s) and large impactor radii ($\ge 75$ km) the rocky material excavated from Enceladus' core and orbiting around the moon is still moving with significant velocity at the end of the simulation.

The thermal softening is an efficient process for large-scale events (Potter et al., 2012). This process is strongly dependent on the pre-impact temperature field that is unfortunately poorly constrained. To test the influence of the pre-impact thermal state, we consider three different pre-impact temperature profiles for the impacted moon (Fig. 3): constant temperature (with $T \sim 100$ K), conductive profile (with a temperature gradient value of 1 K/km), two-layered convective profile (with a core temperature of 450 K and a mantle temperature of 250 K). As illustrated in Fig. 3, a hotter temperature profile in the icy shell strongly enhances the ice flow back and the refill of the core depression. One hour after the impact, a large cavity remains open in the icy mantle for the constant and cold temperature case. For the two-layered convective case where the mantle temperature is close to the melting temperature of ice, the icy mantle rapidly flows back leading to a huge jet
<figure><img src="image_392.png" /><figcaption>**Fig. 1.** Material repartition (left column) and total plastic deformation (right column) as a function of time (from top to bottom) on Enceladus for 3 impact cases: ($v_{\text{imp}} = 10v_{\text{esc}}$, $R_{\text{imp}} = 25$ km) (left), ($v_{\text{imp}} = 10v_{\text{esc}}$, $R_{\text{imp}} = 75$ km) (center) and ($v_{\text{imp}} = v_{\text{esc}}$, $R_{\text{imp}} = 75$ km) (right). In these models, the grid resolution is 1 km in all directions. Here both the rocky core and the icy material are considered as nonporous materials.</figcaption></figure>

of ice at the impact site. However, even if considering three pre-impact thermal states significantly modifies the post-impact dynamics of the icy mantle, this only weakly affects the depth of the depression within the rocky core that ranges between 12 and 15 km (Fig. 3). Hence, in the following, we consider models with a constant pre-impact temperature field.

## 3.2. Influence of ice and rock porosity

The porosity of the material involved during the impact is known to be a key factor in both the fragmentation and disruption of the impactor and the target (Jutzi et al., 2008, 2009), and therefore it may play a role in our results. Enceladus is believed to contain a high degree of porosity, as are many other small bodies in the different populations of asteroids and comets (e.g. Lindsay et al., 2015). To explain the long-wavelength topography of Enceladus, recent models also invoke porosity values ranging between 20% and 30% within the icy mantle of Enceladus (Besserer et al., 2013). We monitored the rocky core deformation as a function of the icy mantle porosity with porosities ranging between 0% and 20%. Similar to the simulations with different initial thermal

conditions (Fig. 3), the dynamics of post-impact ice flow in the deep cavity depends significantly on the porosity, as it affects the ice mechanical properties (Fig. 4). When the ice porosity equals 20% and because the compacted ice is thermally softened, the icy material (which is heated by impact to temperatures up to 250 K) re-fills the impact induced cavity in less than one hour.

Nevertheless, as illustrated in Fig. 5, the effect of the icy mantle porosity on the post-impact core morphology is rather small, at least for initial porosities ranging from 0% to 20% and for impact parameters leading to moderate core deformation ($v_{\text{imp}} = 10v_{\text{esc}}$ and $R_{\text{imp}} = 25$ km). Fig. 6 shows the depth of the impact-induced core depression as a function of the mantle porosity. According to this figure, the depth of the depression ranges between 8 and 13 km. As mentioned earlier (see Fig. 4), the major influence of the mantle porosity is its ability to flow back and refill the core depression. As the impacted ice is severely deformed and compacted during the shockwave propagation, the impact will increase locally the porosity and the temperature of the icy mantle below the impact site.

Figs. 6 and 7 show that the influence of core porosity on core deformation is larger than the corresponding influence of mantle
<figure><img src="image_393.png" /><figcaption>**Fig. 2.** Rocky core morphology as a function of the impactor size and the impact velocity ($v_{esc} = 240$ m/s). In these models the porosity of the icy material is zero. For each morphology, the red circle represents the pre-impact spherical shape of the impacted core. The dashed black line represents Eq. (2) with $A = 2$. Above this critical theoretical line, the impact induced topography is positive. The dotted black line represents Eq. (2) with $A = 1$. Above this critical theoretical line, very highly deformed cores are formed and acoustic fluidization may contribute to their final shape. However the deformation is too large and probably not compatible with the Enceladus morphology. We limit our post impact monitoring to one hour which means that for large impact velocities ($\ge 6$ km/s) and large impactor radii ($\ge 75$ km) the rocky material excavated from Enceladus' core and orbiting around the moon is still moving with significant velocity at the end of the simulation. (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_394.png" /><figcaption>**Fig. 3.** Material repartition one hour after the impact (bottom) for three different pre-impact temperature profiles (top) (with $v_{imp} = 10v_{esc}$, $R_{core} = 160$ km and $R_{imp} = 25$ km). The color of the temperature profile corresponds to the color of the rectangle surrounding the material repartition snapshot. (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_395.png" /><figcaption>**Fig. 4.** Material repartition as a function of the icy mantle porosity one hour after the impact ($v_{imp} = 10 v_{esc}$, $R_{imp} = 25$ km). The rocky core is represented in gray while the icy material is represented in white. In these models, the grid resolution is 1 km in all directions.</figcaption></figure>

<figure><img src="image_396.png" /><figcaption>**Fig. 5.** Rocky core morphology as a function of the icy mantle porosity (with $R_{core}$ = 160 km). For each morphology, the red circle represents the pre-impact spherical shape of the impacted core. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

porosity. Indeed, increasing the porosity of the core from 0% to 50% (and thus increasing its radius from 160 to 200) increases the maximum depth of the depression caused by the impact from ~13 km to ~31.5 km. To explain this feature, two effects shall be invoked. The first one is that increasing the rocky core porosity increases its size to maintain its mass. Hence, the top of the rocky core is closer to the surface and the impactor penetration depth needed to deform the rocky core is reduced accordingly. The second one is that porosity can enhance the rocky core deformation because the core material is less dense and easier to compact. To decipher between these two effects we ran a non-consistent model with a non-porous 200 km rocky core radius surrounded by a 50 km thick icy mantle (Fig. 8, first column). At the end of this model, the depression depth is 18.5 km (compared to 31.5 km when the rocky core porosity is 50% and to 13 km when the rocky core has a radius of 160 km) meaning that both increasing the core size and the porosity favor deeper impact-induced depressions. This also suggests that density/compaction has a greater influence than core radius on the depth of the impact-induced core depression. We also ran a model with a 50% porosity 160 km rocky core radius (Fig. 6) where the obtained depression depth is 15 km (close to value obtained in the non-porous case). In this non-consistent case, a 8 km-thick ice block is trapped between the impactor and the target's core that prevents the formation of a deeper cavity. We should, however, keep in mind that in our simulations, we consider

void porosity, while in reality pores should be filled by liquid water or water ice, which would affect compaction. The results presented here should be considered as an estimation of the maximal effect associated to impact-induced porosity compaction.

### 3.3. Influence of minimum strength values and water ocean

In all the models described above, the minimum strength values were set to $Y_i = 500$ kPa for ice and $Y_s = 10$ MPa for silicate material. These values represent the upper range of the plausible values since recent estimates of the strength of the surface of Comet Tempel-1 obtained minima strength values in the order of 1–10 kPa (Richardson and Melosh, 2013). For the minimum strength of the rocky mantle, this value is also likely to range between the strength of the lunar soil (1 kPa) to the strength of the terrestrial soil (<100 kPa) (Mitchell et al., 1972; Lambe and Whitman, 1979). We have tested the influence of these two parameters using lower values, $Y_i = 10$ kPa and $Y_s = 100$ kPa. As illustrated in Fig. 9 (second column) (called “highly deformable”), decreasing the minimum strength of both the ice and the rocky materials tends to increase the deformability of the rock core leading to both a deeper and wider depression. Ultimately, for a 200 km radius rocky core with 50% porosity (Fig. 9, second columns), the depth of the depression can reach 54.5 km. Here again, the conditions in term of porosity and strength are rather extreme, and the objectives of
<figure><img src="image_397.png" /><figcaption>**Fig. 6.** Depth of the impact induced depression as a function of the rocky core porosity (black circles) and as a function of the icy mantle porosity (red squares) ($v_{imp} = 10v_{esc}$ and $R_{imp} = 25$ km). The vertical line for 0% porosity represents the range of depression depths obtained when considering a 100% icy (lower value) and a 100% rocky (upper value) impactor. The black filled circle at 50% porosity represents the unrealistic case with a core radius of 160 km (while in the other cases the core radius increases with porosity). (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_398.png" /><figcaption>**Fig. 7.** Rocky core morphology as a function of the rocky core porosity. For each morphology, the red circle represents the pre-impact spherical shape of the impacted core. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

this simulation are to illustrate the maximal depression depth that could be generated by a large impact on Enceladus.

Figs. 8 (third and fourth columns) and 9 (third column) show that the presence of a deep water ocean (considered as an inviscid fluid with a density of 910 kg/m³) above the rocky core tends to reduce the impact-induced deflection of the core surface. Liquid water and water ice have comparable compressibility, water being slightly more compressible. The main difference concerns their resistance to shear. Liquid water has no strength (and is considered a completely inviscid material in the simulation), while ice has some strength. In the presence of liquid water, there is complete mechanical decoupling of shear deformation between the water and the core, whereas in the latter case shear stresses exist at the ice-core boundary. In the presence of the water ocean, the lateral extent of the morphology anomaly as well as its depth are smaller than without an ocean. Indeed, for $R_{core} = 160$ km, the depth of the impact induced cavity decreases from 13 km without an ocean to 3.5 km with an ocean. For $R_{core} = 200$ km and $\phi = 50\%$, the depth of the impact induced cavity decreases from 31.5 km

without an ocean to 22.5 km with an ocean. This tends to illustrate that it is easier to enhance post-impact negative topography anomalies in the absence of a water ocean. Including a thick sub-surface water ocean has the opposite effect of increasing the impact velocity or the impactor size, because it concentrates deformation in the ice mantle above, decoupling it from the rocky core below. On the other hand, the presence of the ocean seems to enhance the plastic strain in the deepest part of the core (Fig. 8, third and fourth columns). In parallel to compaction, impact-induced fracturing is likely to generate a porosity increase (via the dilatancy process) (Collins, 2014) that could in return favor fluid circulation within the deformed rocky core.

## 4. Discussion and conclusion

In order to investigate the morphological consequences of collisions between differentiated impactors and Enceladus, we performed numerical impact simulations for impactor radii and
<figure><img src="image_399.png" /><figcaption>**Fig. 8.** Material repartition (left column) and total plastic deformation (right column) as a function of time (from top to bottom) on Enceladus for $R_{core} = 200$ km, ($v_{imp} = 10v_{esc}$ and $R_{imp} = 25$ km). We consider 4 models: a non-consistent non-porous rocky core (first column), a porous rocky core with a porosity of 50% (second column), a non-consistent non-porous rocky core overlaid by a 20 km thick water ocean (third column) and a porous rocky core with a porosity of 50% overlaid by a 20 km thick water ocean (fourth column). In these models, the grid resolution is 1 km in all directions.</figcaption></figure>

velocities ranging between 10% and 40% Enceladus' radius and 1–10 times Enceladus' escape velocity (0.24–2.4 km/s), and for various assumptions for the structure and mechanical properties of Enceladus' interior. Our results showed that the dynamical response of the icy mantle to the impact is strongly dependent on the assumed thermo-mechanical properties for the ice. However, the icy mantle response has minor effects on the impact-induced deformation of the rock core. Only the presence of an internal ocean between the icy mantle and the rock core can significantly limit the rock core deformation.

Our simulations showed that the main controlling parameters for the post-impact shape of Enceladus' rock core are the impactor radius and velocity. We have identified three regimes: (1) For low energy impacts (≤1.5–2 × 10²³ J), the impactors do not pass completely through the icy mantle and the core surface remains unmodified. The rock core of the impactors are deformed by the impact events, but remains trapped within the icy mantle. The impactor core embedded in the icy mantle would then progressively sink and spread, leading to a positive core topographic anomaly. (2) For more energetic impacts, the impactors completely penetrate though the icy mantle and hit the core surface. The impact leads to a negative core topography surrounded by a positive anomaly of smaller amplitude. The depth and lateral extent of the excavated area is mostly determined by the impactor radius and velocity. The shape of the excavated area can be significantly enhanced for high core porosity and very low material strengths, but its amplitude and extent remain primarily determined by the impactor parameters. In this regime, accounting for the acoustic

fluidization does not change the final core morphology (not shown here). (3) For even more energetic impacts, the core is very strongly deformed, which does not appear to be compatible with Enceladus' core morphology (see Fig. 2). Our simulations of these events do not follow the full evolution of the impact scenario so we cannot predict the final core structure; however, it is likely that some of these events lead to full body disruption and that, in non-disruptive impacts, acoustic fluidization may contribute to the final shape of the rocky core and would therefore need to be included to analyze possible outcomes.

For impact velocities higher than 2.4 km s⁻¹(10 × vesc), moderate deformation of the core is possible only for impactors smaller than 25 km. During the Late Heavy Bombardment, high-velocity collisions with impactors exceeding 20 km is likely and therefore, as recently highlighted by Movshovitz et al. (2015), full disruption and re-accretion of the satellite may have occurred possibly several times during this period. This implies that any large impact leaving a long-wavelength signature on the core shape should have taken place after the Late Heavy Bombardement. This also requires relatively low velocity impacts, and therefore encounter with planetocentric bodies rather than with heliocentric bodies. Alternatively, as proposed by Charnoz et al. (2011), Enceladus may have formed late during the history of the Saturn system, thus limiting the risk of full disruption. Following the model of Charnoz et al. (2011), Enceladus may have accreted from a swarm of differentiated embryos emerging from the outer edge of a massive ring system. In such a model, multiple low velocity collisions between decametric differentiated impactors and a growing Enceladus are expected.
<figure><img src="image_400.png" /><figcaption>**Fig. 9.** Rocky core morphology for different pre-impact core radii ($R_{core}$ = 160 km (top) and 200 km (bottom)). First and third columns: $Y_i$ = 500 kPa and $Y_s$ = 10 MPa, second column (“highly deformable”) $Y_i$ = 10 kPa and $Y_s$ = 100 kPa. In the third column we consider a water ocean (with a thickness of 20 km) above the rocky core. For each morphology, the red circle represents the pre-impact spherical shape of the impacted core. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

The irregular core shape of Enceladus, as constrained from Cassini gravity and topography data (McKinnon, 2013; Lefèvre et al., 2015), may constitute a record of this accretional process.

Various processes will probably alter the core topography after an impact event, so that the amplitude of core deflection predicted in our simulations should be considered as an upper limit. Rock fragments would be likely transported by the ice flow back to the impact cavity, filling partly the impact-induced depression. Even if the core is relatively cold, topography relaxation may occur to some extent, especially for low-strength rock material. Prolonged water interactions may also partly erode the topography and again reduce the topography anomaly. Detailed modelling of the subsequent topography evolution is beyond the scope of the present study, and will require future modeling effort. The 2D nature of our simulations also optimizes the amplitude of impact-induced core deflection as only head-on collisions can be considered. It is known that impact angle affects the strength and distribution of the shock wave generated in the impact and therefore the perturbed region (e.g. Pierazzo and Melosh, 2000). For more oblique impacts, the impactor kinetic energy will be more efficiently transferred to the icy mantle, leading to a more efficient deformation of the icy mantle and a larger amount of escaping materials (e.g. Korycansky and Zahnle, 2011). The volume of icy mantle affected by the impact, which is already large for head-on collisions as shown with our 2D simulations, will be further increased. Another limitation of our modelling approach is the assumption regarding the mechanical properties of the rock core. We considered dunite with various degree of void porosity as representative of the rock core composition, since a relatively well-defined equation of state exists for this material (Davison et al., 2010). Based on the interpretation of the Cassini gravity data, which suggest a low density core (2400 kg m⁻³, less et al. (2014)), the rock core may contain a significant fraction of highly

hydrated minerals, as well as free water or/and ice in rock pores. Currently, we are not able to consider a mixture of ice and rocks for both the impactor's core and the target's core. However, to estimate an upper limit of the deformation, we have performed a run corresponding to our classical impact model ($v_{imp} = 10v_{esc}$ and $R_{imp} = 25$ km) with 100% ice-filled pores (i.e. a core made of pure ice). In this unrealistic case (not shown here), the impactor's core is eventually buried at a depth of ~170 km (i.e. 80 km below the core-mantle boundary) which is far larger than the depth of the depression (~30 km) obtained with a 50% porous rocky core. This limitation also stands for the structure of the impactor's core that is likely to have remained undifferentiated in the context of an early formation. To estimate the influence of the impactor's degree of differentiation, we have also considered the $v_{imp} = 10v_{esc}$ case with a 25 km radius impactor made of pure ice and an impactor made of pure dunite. In the first case, the impact induces a flattening of ~0.4 km at the core's surface below the impact site (see Fig. 6). In the second case, the impact induces a flattening of ~23.2 km. This result, even if performed for an unrealistic water ice content, suggests the ice/rock ratio in the core may play a strong influence on the response of the core to large impacts. This suggests that the results presented here should be considered valid only for differentiated interior models with rock-dominated core and a relatively small porosity content (<10–20%). Future works are required to constrain more precisely the effect of hydrated minerals and mixture with high ice-water/rock ratio in the interior.

Large impacts are likely to modify the ice/rock ratio by eroding significantly the shallower part of the impacted moon. Our results show that vertical impacts with $v_{imp} > 6v_{esc}$ and $R_{imp} > 75$ km, can erode up to half the ice volume from the impacted body (Fig. 1, second column). Several factors such as a hot, porous pre-impact mantle and the presence of a deep water ocean increase the ability of the icy mantle to deform. Hence, these parameters are also likely
to influence the post-impact ice/rock ratio by decreasing the fraction of ice in the post-impact moon. The impact angle is another key parameter that governs the fraction of escaped material (e.g. Korycansky and Zahnle, 2011). However, to limit the computational time and as we have restricted our study to vertical impacts, monitoring the long-term evolution of the ice/rock ratio is beyond the scope of our study.

Despite the limitations, the simulations we performed highlight the crucial role played by impacts on the evolution of Enceladus. Besides explaining the irregular shape of the core, impacts also provide efficient mechanisms to enhance thermo-chemical exchanges between the deep interior and the surface. For models with an internal water ocean, we can see that a large volume of the ocean is temporarily exposed to the surface, thus potentially releasing a large fraction of volatile initially stored dissolved in the ocean. Large impacts cause a strong damage of the ice on a very large portion of the icy mantle, which will likely have consequences on the subsequent convective mantle dynamics and interaction with the fractured surface. These also lead to a large plastic strain in the rock core underneath the impact site, which may enhance macroporosity. This would promote fluid circulation throughout a large fraction of the core, favoring serpentinization processes (Malamud and Prialnik, 2013) and hydrothermal activities (e.g. Hsu et al., 2015). Further modeling efforts will be needed to understand the consequences of such impact events on the long-term evolution of Enceladus. Lastly, the effects of large impacts are not confined to Enceladus. Similar effects are very likely on the other moons of Saturn as well as on other planetary objects, such as Ceres (e.g. Davison et al., 2015; Ivanov, 2015) and Pluto (e.g. Bray and Schenk, 2015) for which impact bombardment has probably played a key role in their evolution.

## Acknowledgments

J. Monteux was funded by Agence Nationale de la Recherche (Accretis decision No. ANR-10-PDOC-001-01). The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 259285). We thank James Roberts and the anonymous reviewer for constructive comments. We gratefully acknowledge the developers of iSALE, including Kai Wünnemann, Boris Ivanov, Jay Melosh, and Dirk Elbeshausen. We also thank Katarina Miljkovich and Tom Davison for useful discussions.

## References

* Amsden, A.A., Ruppel, H.M., Hirt, C.W., 1980. SALE: Simplified ALE Computer Program for Fluid Flow at All Speeds, Tech. Rep. LA-8095. Los Alamos Natl. Lab., Los Alamos, NM.

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

* Asphaug, E., Reufer, A., 2013. Late origin of the Saturn system. Icarus 223, 544–565. http://dx.doi.org/10.1016/j.icarus.2012.12.009.

* Benz, W., Asphaug, E., 1999. Catastrophic disruptions revisited. Icarus 142, 5–20. http://dx.doi.org/10.1006/icar.1999.6204.

* Benz, W., Cameron, A.G.W., Melosh, H.J., 1989. The origin of the moon and the single impact hypothesis. III. Icarus 81, 113–131. http://dx.doi.org/10.1016/0019-1035(89)90129-2.

* Besserer, J. et al., 2013. Convection-driven compaction as a possible origin of Enceladus's long wavelength topography. J. Geophys. Res. (Planets) 118, 908–915. http://dx.doi.org/10.1002/jgre.20079.

* Bray, V.J., Schenk, P.M., 2015. Pristine impact crater morphology on Pluto – Expectations for new horizons. Icarus 246, 156–164. http://dx.doi.org/10.1016/j.icarus.2014.05.005.

* Bray, V.J. et al., 2008. The effect of target properties on crater morphology: Comparison of central peak craters on the moon and Ganymede. Meteoritics Planet. Sci. 43, 1979–1992. http://dx.doi.org/10.1111/j.1945-5100.2008.tb00656.x.

* Běhounková, M. et al., 2010. Coupling mantle convection and tidal dissipation: Applications to Enceladus and Earth-like planets. J. Geophys. Res. (Planets) 115, E09011. http://dx.doi.org/10.1029/2009JE003564.

* Charnoz, S. et al., 2011. Accretion of Saturn's mid-sized moons during the viscous spreading of young massive rings: Solving the paradox of silicate-poor rings versus silicate-rich moons. Icarus 216, 535–550. http://dx.doi.org/10.1016/j.icarus.2011.09.017.

* Collins, G.S., 2014. Numerical simulations of impact crater formation with dilatancy. J. Geophys. Res. Planets 119, 1–20. http://dx.doi.org/10.1002/2014JE004708.

* Collins, G.C., Goodman, J.C., 2007. Enceladus' south polar sea. Icarus 189, 72–82. http://dx.doi.org/10.1016/j.icarus.2007.01.010.

* Collins, G.S., Davison, T.M., Ciesla, F.J., 2011. Numerical simulations of sub-catastrophic porous planetesimal collisions. Proc. Lunar Sci. Conf. 42, 1933.

* Collins, G.S., Melosh, H.J., Wünnemann, K., 2013. Improvements to the $\epsilon-\alpha$ porous compaction model for simulating impacts into high-porosity Solar System objects. Int. J. Impact Eng. 38, 434–439.

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

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

* Davison, T.M. et al., 2015. Impact bombardment of ceres. Proc. Lunar Sci. Conf. 46, 2116.

* Dwyer, C.A. et al., 2013. The influence of imperfect accretion and radial mixing on ice:rock ratios in the Galilean satellites. Icarus 225, 390–402. http://dx.doi.org/10.1016/j.icarus.2013.03.025.

* Giese, B. et al., 2008. Enceladus: An estimate of heat flux and lithospheric thickness from flexurally supported topography. Geophys. Res. Lett. 35, L24204. http://dx.doi.org/10.1029/2008GL036149.

* Hansen, C.J. et al., 2006. Enceladus' water vapor plume. Science 311, 1422–1425. http://dx.doi.org/10.1126/science.1121254.

* Hsu, H.-W. et al., 2015. Ongoing hydrothermal activities within Enceladus. Nature 519, 207–210. http://dx.doi.org/10.1038/nature14262.

* Iess, L. et al., 2014. The gravity field and interior structure of Enceladus. Science 344, 78–80. http://dx.doi.org/10.1126/science.1250551.

* Ivanov, B.A., 2015. Ceres: possible records of giant impacts. Proc. Lunar Sci. Conf. 46, 1077.

* Ivanov, B.A., Deniem, D., Neukum, G., 1997. Implementation of dynamic strength models into 2D hydrocodes: Applications for atmospheric breakup and impact cratering. Int. J. Impact Eng. 20, 411–430.

* Jutzi, M., Benz, W., Michel, P., 2008. Numerical simulations of impacts involving porous bodies. I. Implementing sub-resolution porosity in a 3D SPH hydrocode. Icarus 198, 242–255. http://dx.doi.org/10.1016/j.icarus.2008.06.013.

* Jutzi, M. et al., 2009. Numerical simulations of impacts involving porous bodies. II. Comparison with laboratory experiments. Icarus 201, 802–813. http://dx.doi.org/10.1016/j.icarus.2009.01.018.

* Korycansky, D.G., Zahnle, K.J., 2011. Titan impacts and escape. Icarus 211, 707–721. http://dx.doi.org/10.1016/j.icarus.2010.09.013.

* Lambe, T., Whitman, R., 1979. Soil mechanics, SI version. In: Series in Soil Engineering. Wiley.

* Lefèvre, A. et al., 2015. Enceladus' internal ocean constrained from Cassini gravity and topography data. In: European Planetary Science Conference.

* Lindsay, S.S. et al., 2015. Composition, mineralogy, and porosity of multiple asteroid systems from visible and near-infrared spectral data. Icarus 247, 53–70. http://dx.doi.org/10.1016/j.icarus.2014.08.040.

* Malamud, U., Prialnik, D., 2013. Modeling serpentinization: Applied to the early evolution of Enceladus and Mimas. Icarus 225, 763–774. http://dx.doi.org/10.1016/j.icarus.2013.04.024.

* McKinnon, W.B., 2013. The shape of Enceladus as explained by an irregular core: Implications for gravity, libration, and survival of its subsurface ocean. J. Geophys. Res. (Planets) 118, 1775–1788. http://dx.doi.org/10.1002/jgre.20122.

* Melosh, H.J., Ivanov, B.A., 1999. Impact crater collapse. Annu. Rev. Earth Planet. Sci. 27, 385–415. http://dx.doi.org/10.1146/annurev.earth.27.1.385.

* Melosh, H.J., Ryan, E.V., Asphaug, E., 1992. Dynamic fragmentation in impacts – Hydrocode simulation of laboratory impacts. J. Geophys. Res. 97, 14,73. http://dx.doi.org/10.1029/92JE01632.

* Mitchell, J.K. et al., 1972. Mechanical properties of lunar soil: Density, porosity, cohesion and angle of internal friction. In: Metzger, A.E. et al. (Eds.), Proc. Lunar Sci. Conf., vol. 3, p. 3235.

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

* Movshovitz, N. et al., 2015. Disruption and reaccretion of midsized moons during an outer Solar System Late Heavy Bombardment. Geophys. Res. Lett. 42, 256–263. http://dx.doi.org/10.1002/2014GL062133.

* Murr, L. et al., 1998. The low-velocity-to-hypervelocity penetration transition for impact craters in metal targets. Mater. Sci. Eng. 256, 166–182.

* Neish, C.D., Lorenz, R.D., 2012. Titan's global crater population: A new assessment. Planet. Space Sci. 60, 26–33. http://dx.doi.org/10.1016/j.pss.2011.02.016.

* Nimmo, F. et al., 2007. Shear heating as the origin of the plumes and heat flux on Enceladus. Nature 447, 289–291. http://dx.doi.org/10.1038/nature05783.

* Ohnaka, M., 1995. A shear failure strength law of rock in the brittle-plastic transition regime. Geophys. Res. Lett. 22, 25–28. http://dx.doi.org/10.1029/94GL02791.

* Orphal, D.L. et al., 1980. Impact melt generation and transport. In: Bedini, S.A. (Ed.), Proc. Lunar Sci. Conf. vol. 11, pp. 2309–2323.
Pierazzo, E., Melosh, H.J., 2000. Melt production in oblique impacts. *Icarus* 145, 252–261. http://dx.doi.org/10.1006/icar.1999.6332.

Porco, C.C. et al., 2006. Cassini observes the active south pole of Enceladus. *Science* 311, 1393–1401. http://dx.doi.org/10.1126/science.1123013.

Potter, R.W.K. et al., 2012. Constraining the size of the South Pole-Aitken basin impact. *Icarus* 220, 730–743. http://dx.doi.org/10.1016/j.icarus.2012.05.032.

Richardson, J.E., Melosh, J.H., 2013. An examination of the Deep Impact collision site on Comet Tempel 1 via Stardust-NExT: Placing further constraints on cometary surface properties. *Icarus* 222, 492–501. http://dx.doi.org/10.1016/j.icarus.2012.04.019.

Roberts, J.H., 2015. The fluffy core of Enceladus. *Icarus* 258, 54–66. http://dx.doi.org/10.1016/j.icarus.2015.05.033.

Schenk, P., 2011. Geology of Mimas? In: *Lunar and Planetary Science Conference, Lunar and Planetary Inst. Technical Report*, vol. 42, p. 2729.

Schenk, P.M., McKinnon, W.B., 2009. One-hundred-km-scale basins on Enceladus: Evidence for an active ice shell. *Geophys. Res. Lett.* 36, L16202. http://dx.doi.org/10.1029/2009GL039916.

Schenk, P. et al., 2012. The geologically recent giant impact basins at Vesta south pole. *Science* 336, 694. http://dx.doi.org/10.1126/science.1223272.

Schubert, G., Stevenson, D.J., Ellsworth, K., 1981. Internal structures of the Galilean satellites. *Icarus* 47, 46–59. http://dx.doi.org/10.1016/0019-1035(81)90090-7.

Schubert, G. et al., 2007. Enceladus: Present internal structure and differentiation by early and long-term radiogenic heating. *Icarus* 188, 345–355. http://dx.doi.org/10.1016/j.icarus.2006.12.012.

Showman, A.P., Han, L., Hubbard, W.B., 2013. The effect of an asymmetric core on convection in Enceladus' ice shell: Implications for south polar tectonics and heat flux. *Geophys. Res. Lett.* 40, 5610–5614. http://dx.doi.org/10.1002/2013GL057149.

Spencer, J.R., Nimmo, F., 2013. Enceladus: An active ice world in the Saturn system. *Annu. Rev. Earth Planet. Sci.* 41, 693–717. http://dx.doi.org/10.1146/annurev-earth-050212-124025.

Spencer, J.R. et al., 2006. Cassini encounters Enceladus: Background and the discovery of a south polar hot spot. *Science* 311, 1401–1405. http://dx.doi.org/10.1126/science.1121661.

Thomas, P.C. et al., 2007. Shapes of the saturnian icy satellites and their significance. *Icarus* 190, 573–584. http://dx.doi.org/10.1016/j.icarus.2007.03.012.

Tobie, G., Čadek, O., Sotin, C., 2008. Solid tidal friction above a liquid water reservoir as the origin of the south pole hotspot on Enceladus. *Icarus* 196, 642–652. http://dx.doi.org/10.1016/j.icarus.2008.03.008.

Waite, J.H. et al., 2006. Cassini ion and neutral mass spectrometer: Enceladus plume composition and structure. *Science* 311, 1419–1422. http://dx.doi.org/10.1126/science.1121290.

Wünnemann, K., Collins, G.S., Melosh, H.J., 2006. A strain-based porosity model for use in hydrocode simulations of impacts and implications for transient crater growth in porous targets. *Icarus* 180, 514–527. http://dx.doi.org/10.1016/j.icarus.2005.10.013.
<figure><img src="image_401.png" /></figure>

# Water and the Interior Structure of Terrestrial Planets and Icy Bodies

J. Monteux<fnref n="1" /> · G.J. Golabek<fnref n="2" /> · D.C. Rubie<fnref n="2" /> · G. Tobie<fnref n="3" /> · E.D. Young<fnref n="4" />

Received: 30 July 2016 / Accepted: 11 January 2018
© Springer Science+Business Media B.V., part of Springer Nature 2018

**Abstract** Water content and the internal evolution of terrestrial planets and icy bodies are closely linked. The distribution of water in planetary systems is controlled by the temperature structure in the protoplanetary disk and dynamics and migration of planetesimals and planetary embryos. This results in the formation of planetesimals and planetary embryos with a great variety of compositions, water contents and degrees of oxidation. The internal evolution and especially the formation time of planetesimals relative to the timescale of radiogenic heating by short-lived $^{26}$Al decay may govern the amount of hydrous silicates and leftover rock-ice mixtures available in the late stages of their evolution. In turn, water content may affect the early internal evolution of the planetesimals and in particular metal-silicate separation processes. Moreover, water content may contribute to an increase of oxygen fugacity and thus affect the concentrations of siderophile elements within the silicate reservoirs of Solar System objects. Finally, the water content strongly influences the differentiation rate of the icy moons, controls their internal evolution and governs the alteration processes occurring in their deep interiors.

**Keywords** Water · Differentiation · Internal evolution · Terrestrial planets · Icy moons

The Delivery of Water to Protoplanets, Planets and Satellites

Edited by Michel Blanc, Allessandro Morbidelli, Yann Alibert, Lindy Elkins-Tanton, Paul Estrada, Keiko Hamano, Helmut Lammer, Sean Raymond and Maria Schönbächler

J. Monteux
julien.monteux@uca.fr

<footnote n="1">1 CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, Université Clermont Auvergne, 63000 Clermont-Ferrand, France</footnote>

<footnote n="2">2 Bayerisches Geoinstitut, Universität Bayreuth, Universitätsstrasse 30, 95440 Bayreuth, Germany</footnote>

<footnote n="3">3 Laboratoire de Planétologie et Géodynamique, UMR-CNRS 6112, Université de Nantes, 44322 Nantes Cedex, France</footnote>

<footnote n="4">4 Department of Earth, Planetary, and Space Sciences, University of California Los Angeles, Los Angeles, CA 90095, USA</footnote>
# 1 Introduction

The Solar System displays a strong dichotomy between the inner region that is characterized by relatively dry planetary objects having a very small water fraction (<0.1%), and the outer solar system where water ice constitutes a large fraction of solid material (>20%), which is inherited from accretion processes. When compared with other planetary systems, the Solar System seems however rather unusual. Exoplanet surveys have revealed that planets intermediate in mass between Earth and Neptune are surprisingly common, but are notably absent in the Solar System (Mayor et al. 2014; Howard et al. 2012; Marcy 2014). Model mass-radius relationships indicate a great diversity of interior compositions and atmospheric extents for the Super-Earth/Mini-Neptune-planet class (e.g. Howard et al. 2012), suggesting a wide range of volatile contents (including water) and compositions.

The distribution of water in the planetary system is controlled by the temperature structure in the protoplanetary disk and dynamics and migration of planetesimals and planetary embryos (Raymond et al. 2004; Cowan and Abbot 2014; O'Brien et al. 2014; Rubie et al. 2015a). This results in the formation of planetesimals and planetary embryos with a great variety of compositions and water contents. Subsequent accretionary processes can lead to the formation of planets that range from completely dry planets to "water worlds" with more than 100 Earth-oceans of water (1 ocean = $1.5 \times 10^{21}$ kg H$_2$O) (Raymond et al. 2004; Cowan and Abbot 2014; O'Brien et al. 2014; Rubie et al. 2015a). As the evolution of water content inside a growing proto-planet is a strong function of pressure-temperature conditions, the early stages of planetary formation that govern their accretion rates and their early heat budget may have played a major role in the water distribution in the solar system. Indeed, depending on the accretion rate and on the formation time of the building blocks of the terrestrial planets relative to the timescale of radiogenic heating by the decay of short-lived $^{26}$Al, water can either be incorporated as an ice-rock mixture or in the form of hydrated silicates or be absent. Aqueous alteration appears to be common in carbonaceous chondrites (Kerridge and Bunch 1979; Zolensky et al. 1989; Krot et al. 1998). The proportion of planetesimal water incorporated into rock depends critically on whether water was mobile or stagnant (e.g., Young 2001; Young et al. 1999, 2003; Bland et al. 2009; Fu et al. 2017; Castillo-Rogez and Young 2017). In turn, the water content of Solar System objects likely played a key role in the internal evolution of the solar system bodies. For example, water dissolved in silicate minerals is known to significantly reduce mantle viscosity and melting temperature (e.g. Hirth and Kohlstedt 1996; Katz et al. 2003; Hirschmann 2006). It has therefore a potential influence on the early differentiation processes and a strong feedback on the thermal state of the mantle and the vigor of convection (e.g. Korenaga 2010).

Water is also the major constituent of the moons orbiting Jupiter (except for Io), Saturn, Uranus and Neptune, and all planetary objects beyond Neptune (e.g. Hussmann et al. 2015). In these icy objects, water is likely present either as thick ice shells or deep water oceans (e.g. Hussmann et al. 2006) and in some cases is being ejected from venting plumes composed mainly of water vapor and ice particles (e.g. Porco et al. 2006 for Enceladus, Roth et al. 2014 for Europa). The combined effect of accretional heating, radiogenic decay by short-lived isotopes, tidal heating associated with despinning, and viscous heating associated with sinking negative rock diapirs may increase the internal temperature above the melting point of ice (e.g. Tobie et al. 2013). This may result in a partial to full ice-rock separation. During this differentiation as well as at the present-day, interactions between solid rock and liquid water might have promoted chemical reactions facilitated by the presence of ammonia (Schubert et al. 2010; Sohl et al. 2010). It is also possible that aqueous alteration
might have led to serpentinization of olivine-rich rocks and to the formation of a highly hy-
drated rocky core that may result from a rapid differentiation process that implies large-scale
melting events, as proposed for Saturn’s moon Titan (e.g. Castillo-Rogez and Lunine 2010).

Finally, comets and comet-like objects are another class of objects in which ice has been
detected as a major component (80%) (Encrenaz 2008). These bodies, with radii usually
smaller than 10 km, were formed in the outer Solar System. Comets are undifferentiated
and, hence, represent weakly altered remnants of the early stages of solar-system formation.
Because of radiogenic heating from $^{26}$Al and $^{60}$Fe, the interior of comets may reach temper-
atures above the melting point of water and ultravolatile species may be lost (Prialnik et al.
2008; Mousis et al. 2017). In the following study, we do not consider the comets and we
will only focus on objects with radii larger than 25 km (see review of Prialnik et al. 2008 for
more details).

The structure of this article is as follows: We first detail the distribution of water in the protoplanetary disc and its effect on oxidation. Second, we present a discussion of the effect of internal evolution of planetesimals on water content based on numerical models. In particular we investigate the influence of planetesimal radius and formation time. We then discuss the role of water on the internal evolution of proto-planets and terrestrial bodies focusing particularly on metal/silicate separation and oxygen fugacity evolution. Finally we briefly describe the internal evolution of icy satellites from their accretion to their present-day activity and emphasize the role of water in differentiation and alteration processes.

## 2 Distribution of Water in the Protoplanetary Disc and Its Effect on Oxidation

Protoplanetary disks of gas and dust are by necessity characterized by radial increases in
temperature and pressure towards the nascent star at their center. The inner portions of the
disk are heated by two mechanisms. So-called “passive” disks are heated at their surfaces by
illumination from the central star. Viscous heating in “active” accretion disks can warm up
the inner portions of the disk to temperatures well above those imposed by passive illumina-
tion alone. The degree of viscous heating depends on the mechanism for angular transport
(accretion onto the star) (e.g., Bai and Stone 2013). The greater the effective viscosity of
the disk, the greater the heating related to accretion. In all cases, there will be specific radii
corresponding to the condensation fronts for various volatile components (referred to as
“snowlines”). From the perspective of volatiles and planet formation, the most important
of these may be that for water. The water snowline separating the inner regions of the disk
where water vapor is stable from the more distal regions of the disk where water ice is stable
is thought to exert a first-order control on the chemistry of planetesimals. The snowline in
an accretion disk is closest to the star at the midplane where the intervening mass of gas and
dust shields the disk from light from the central star. Higher in the disk, the snowline curves
to greater radial distances due to heating of the surfaces of the disk by the central star. We are
most concerned with the position of the water snowline at the midplane when contemplating
planetesimal formation. However it must be kept in mind that the time-dependent position
of the snowline in the midplane and its effects on planetesimal formation over time are still
debated (e.g., Morbidelli et al. 2016).

The bulk compositions of planetesimals and planetary embryos in the protoplanetary
disc, including water and some volatile element abundances, have been constrained by core
formation modelling (Rubie et al. 2015a, 2016). Mantle abundances of siderophile elements
<figure><img src="image_402.png" /><figcaption>**Fig. 1** The composition of initial bodies (embryos and planetesimals) as a function of heliocentric distance determined in the accretion/core formation model of Rubie et al. (2015a). At heliocentric distances less than ~1 au, 99% of all Fe is present as metal and ~18% of total Si is dissolved in the metal—which makes the bulk composition very reduced. Beyond ~1 au, compositions become increasingly oxidized with increasing heliocentric distance because the fraction of Fe present as metal decreases and the fraction of Fe present as FeO in the silicate increases. Beyond 6–7 au, compositions are fully oxidized (no metal) and contain 20 wt% water—these are bodies that originally accreted beyond the snowline. The *red arrows* indicate values of parameters that were refined by least squares</figcaption></figure>

(Ni, Co, W etc.) can be reproduced in multistage core formation models when the initial 60–70% of Earth accretes from highly reduced material and the final 30–40% from relatively-oxidized material, in combination with the effects of Si partitioning into the core (Rubie et al. 2011). The late accretion of oxidized material is also consistent with the relatively late accretion of volatile elements to the Earth (e.g. Schönbächler et al. 2010; Dauphas 2017). This scenario was developed further by Rubie et al. (2015a, 2016) by formulating a combined accretion/core formation model. This was based on “Grand Tack” N-body accretion models (Walsh et al. 2011; O’Brien et al. 2014) that started with up to ~200 embryos (distributed from 0.7 to 3.0 au) embedded in a protoplanetary disc consisting of up to 4400 planetesimals (distributed from 0.7 to 9.5 au). Each accretional collision between embryos and other embryos and planetesimals was considered to result in an episode of core formation, provided the impactor was metal-bearing. The bulk composition of each of the starting bodies was defined in terms of CI ratios of non-volatile elements, with oxygen and water contents varying as a function of heliocentric distance. Based on these bulk compositions, the compositions of equilibrated metal and silicate in each core-formation event were determined by a combined mass balance/element partitioning approach, thus enabling the evolving compositions of the mantles and cores of all growing planets to be modeled. The composition of the model Earth in each simulation was fit to Earth’s mantle composition by a least squares refinement of the metal-silicate equilibration pressure and several parameters that describe the compositions of initial bodies (Fig. 1). Best fits to Earth’s mantle composition were obtained when compositions of bodies close to the Sun are highly reduced and those further out are relatively oxidized (Fig. 1).

The gradient in Fig. 1 requires explanation in the context of the young solar protoplanetary disk. Dynamical models have shown that radial mixing within the protoplanetary disk was likely to occur especially during the late stages of planetary formation (e.g. Chambers
<figure><img src="image_403.png" /><figcaption>**Fig. 2** Schematic illustrating the enhancement of the fayalite component in olivine due to an increase in water. This scheme could have operated at moderately high temperatures in the protoplanetary disk where dust was processed</figcaption></figure>

2013; Raymond et al. 2004). Variation in oxygen fugacity across the protoplanetary gaseous disk is a longstanding question. At temperatures greater than about 500 to 600 K or so, Fe₂SiO₄, representing oxidized iron, is unstable at solar oxygen fugacities. This temperature is lower than the ca. 1400 K required for diffusion and annealing of silicate grains (Figs. 2 and 3). The solution is to pump up the oxidation state of the vapor phase by infusing it with water inside of the snowline. As the protoplanetary disk evolved, small icy bodies will have drifted inwards of the snowline, resulting in addition of H₂O vapor by subsequent sublimation of water ice. The longevity of this vapor depended critically on the competing rates of inward migration across the snowline and the outward diffusive flux of water driven by freezing of water at the snowline (e.g., Cuzzi and Zahnle 2004). The result of the competing fluxes may have been a pulse of elevated H₂O/H₂ just inside the snowline (Fig. 4). Such a localized pulse of water vapor has been observed in the TW Hydra transition disk at a distance of about 4 au from the central star (Zhang et al. 2013). The significant amount of oxidized iron in grains comprising chondrite matrix, even in rocks with limited evidence for water-rock reactions, is testament to the fact that such a process must have occurred where planetesimals were forming (Grossman et al. 2012). The amount of water required can be estimated from the equilibrium relationship

$$H_2 + 1/2O_2 \rightleftharpoons H_2O$$

for which we obtain the equation for fO₂ (e.g., Krot et al. 2000)

$$\log(f_{O_2}) \approx 2\log\left(\frac{H_2O}{H_2}\right) + 5.67 - \frac{256664}{T(\text{K})}. \quad (1)$$

Inspection of this equation shows that in order to change oxygen fugacity by ~5 log units as required to shift from solar to “planetary” oxidation states, H₂O/H₂ must have increased by a factor of ~375 at inner disk temperatures of ~1400 K (Fig. 5). Geometric considerations suggest that focusing solar abundances of water into the inner 5 au of the disk by inward transport can only raise the H₂O/H₂ ratio by about a factor of 10 (where, for example, surface density varies as 1/R where R is the radial distance from the Sun).

The picture that emerges is one in which ice accreted within planetesimals outside of the snowline may have reacted with the host rock as the planetesimals evolved in response to heating by ²⁶Al (see also Sect. 3). Liquid water reactions with rock will have oxidized the rock. Inside of the snowline there was likely a region where inward migration of ice enhanced the H₂O/H₂ ratio in the vapor phase just inside of the snowline, but further into the disk in the region of terrestrial planet formation, the ambient conditions favored reduced iron over oxidized iron as a result of higher temperatures and a decrease in gaseous H₂O toward the star. The snowline must have migrated inward with time in response to a decrease
<figure><img src="image_404.png" /><figcaption>**Fig. 3** Schematic diagram showing the relationship between oxygen fugacity and temperature for a solar gas and for the oxidation state of iron. At higher temperatures, reduced iron is favored over the fayalite component for a solar-like gas. At lower temperatures, Fe₂SiO₄ can be stabilized, but whether this reaction can actually occur at lower temperatures is in question. Typical activities for fayalite in olivine at two different temperatures are shown for illustration</figcaption></figure>

<figure><img src="image_405.png" /><figcaption>**Fig. 4** Schematic diagram illustrating the competition between advection of icy material through the snowline of the protoplanetary disk and backwards net transfer of water via diffusion driven by a buildup of water vapor inside of the snowline. The process is described by Cuzzi and Zahnle (2004) and is consistent with Spitzer and Herschel observations of TW Hydra (Zhang et al. 2013)</figcaption></figure>

in the stellar accretion rate, and the best fit shown in Fig. 1 suggests that its time-integrated “fossilized” location (see below) was somewhere between 6 and 2 au.

The position of the midplane water snowline as a function of time can be obtained from the temperature and pressure structure of a generic protoplanetary accretion disk. The *T* and *P* structure of the disk can be modeled using basic principles of dissipative radiative losses resulting from orbital torques. The torques are sustained by the viscosity of the disk that results in the inward spiral of gas towards the growing star. For dimensionless viscosity parameter *α* (e.g., Shakura and Sunyaev 1973), sound speed *c*<sub>s</sub>, and scale height *H*(*R*), the effective viscosity of the disk is *v* = *αc*<sub>s</sub>*H*(*R*). The temperature at the midplane of the disk can be written in terms of the accretion rate of the system, *ṁ*, as follows

$$T_{\text{midplane}} = \left[ \frac{\frac{3}{2} \Sigma(R) \kappa_R}{4} + 1 \right]^{1/4} \left[ \frac{3}{8\pi\sigma} \frac{GM_*\dot{m}}{R^3} \right]^{1/4} \quad (2)$$

where *κ*<sub>*R*</sub> is the Rosseland mean opacity, *Σ*(*R*) is the surface density of the disk that is a function of *R* from the central star, *G* is the gravitational constant, *M*<sub>*</sub> is the mass of the
<figure><img src="image_406.png" /><figcaption>**Fig. 5** The logarithm of oxygen fugacity versus temperature for various oxygen buffers relevant for the early solar system. The heavy dashed line shows the relationship for a solar gas. For comparison, planetary materials (Earth's mantle, basalts from various bodies) exhibit oxygen fugacities closer to the iron-wüstite (IW) and quartz-fayalite-magnetite (QFM) oxygen buffers. Note that the difference between a solar gas and those planetary materials implies an increase of more than 5 orders of magnitude in oxygen fugacity. The region shaded in grey represents oxidation states that favor Ti³⁺ over Ti⁴⁺ in calcium-aluminum-rich inclusions (CAIs) in carbonaceous chondrites. Contours for the activity ratio of the two forms of Ti in CAIs are also shown</figcaption></figure>

central star, $\dot{m}$ is the accretion rate and $\sigma$ is the Stefan-Boltzmann constant. The first term on the right-hand-side of (2) accounts for the optical depth perpendicular to the midplane and is equivalent to $T_{\text{midplane}}/T_{\text{surface}}$. The surface temperature is given by the second term. Equation (2) shows that the temperature of the disk at any radial distance $R$ from the star depends on the rate of mass accretion through the disk. The latter has been measured in young stellar objects today and varies from about $10^{-6} M_{\odot} \text{yr}^{-1}$ early in the history of stellar accretion to $10^{-8} M_{\odot} \text{yr}^{-1}$ (Hartmann 2000) later on.

Evaluation of (2) requires evaluating the surface density with distance from the star, $\Sigma(R)$. The surface density can be written in terms of the viscosity as (Armitage 2010)

$$ \Sigma(R) = \frac{\dot{m}}{3\pi\alpha c_s H(R)}, \qquad (3) $$

where scale height $H$ of the disk as a function of radial distance from the star is

$$ H(R) = \frac{c_s}{\Omega} = \sqrt{\frac{k_b T R^3}{G M_* \mu}}, \qquad (4) $$

and where $\Omega$ is orbital angular velocity, $k_b$ is Boltzmann's constant, $T$ is the temperature and $\mu$ is the average molecular mass for the gas comprising the disk. Equation (2) can be evaluated analytically with reasonable precision by assuming a constant value of 1 m²/kg for $\kappa_R$ and replacing (3) with a parameterization for $\Sigma(R)$ suggested by the mass distribution of planets in the solar system: $\Sigma(R) \sim 1200 \text{ kg/m}^2 (R/5 \text{ au})^{-3/4}$ (compare with Chiang and Youdin 2010). The result for $10^{-7} M_{\odot} \text{yr}^{-1}$ is shown in Fig. 6.

The competing rates of temperature and pressure-dependent adsorption and desorption of water onto grain surfaces in the disk define the position of the water snow line. Pressure
<figure><img src="image_407.png" /><figcaption>Fig. 6 Calculated temperature and pressure vs. distance from the central star for a solar-like protoplanetary accretion disk with a mass accretion rate of $10^{-7}$ solar masses per year</figcaption></figure>

defines the total number density $n_{total}$ of molecules since $P = n_{total}k_bT$. Midplane pressures as a function of radial distance from the star are

$$P(R) = \frac{\Sigma k_b T}{\sqrt{\pi H(R)} \mu} \qquad (5)$$

The kinetic theory of gases yields the rate constant for adsorption of molecules onto grain surfaces:

$$k_{\text{ads}} = \langle \pi r_{\text{grains}}^2 \rangle \hat{V}_{\text{gas}} n_{\text{grains}} \qquad (6)$$

where $k_{\text{ads}}$ is the rate constant (s⁻¹), $r_{\text{grains}}$ is the average radius of the dust grains, $n_{\text{grains}}$ is the number density of the dust grains (cm⁻³), and $\hat{V}_{\text{gas}} = \sqrt{8k_bT/(\pi\mu)}$ is the average gas velocity. Hasegawa et al. (1992) suggest a parameterization for the number density of grains based on total hydrogen such that $n_{\text{grains}} \sim 1.33 \times 10^{-12} n_{\text{H}}$ and $n_{\text{H}} = 2 \times n_{\text{H2}} \sim n_{\text{total}}$. Substitution of this expression into (6), together with the mass of water molecules, yields in SI units

$$k_{\text{ads}} = \langle \pi r_{\text{grains}}^2 \rangle \sqrt{1.18 \times 10^4 T (\text{K})} 1.33 \times 10^{-6} (P_{\text{H}_2}/(k_b T)). \qquad (7)$$

Equation (7) can be evaluated using midplane temperatures and pressures (2), (4), and (5)) and the fact that typical radii for dust grains relevant to the protoplanetary disk environment are 0.1 µm, or $1 \times 10^{-7}$ m. The competing rate constant for desorption of water from grains in SI units is (Willacy et al. 1998):

$$k_{\text{des}} = 1 \times 10^{12} \exp(-4815/T(\text{K})). \qquad (8)$$
<figure><img src="image_408.png" /><figcaption>**Fig. 7** Position of the water snow line as a function of mass accretion rate as defined by the ratio of the adsorption and desorption rate constants</figcaption></figure>

At equilibrium rates of adsorption and desorption of water onto dust grains are equal so we have $k_{ads}n_{H_2O,gas} = k_{des}n_{H_2O,grains}$. Therefore, $k_{ads}/k_{des} \gg 1$ represents water frozen to grains and $k_{ads}/k_{des} \ll 1$ represents where most water molecules exist in the gas phase. A unit value for $k_{ads}/k_{des}$ defines conditions corresponding to the effective position of the snow line in the midplane of the disk. Substitution of (2) and (5) into (7) and (8) yields the plot in Fig. 7 showing the calculated position of the snowline in the disk as a function of accretion rate. The effective freeze out temperature for $\dot{m} = 10^{-8} M_{\odot}/yr$ is ~190 K. The effective temperature for freeze out of water is 170 K for $\dot{m} = 10^{-7} M_{\odot}/yr$. Most importantly, we note that basic physics of accretion disks indicates that the calculated position of the snow line varies from ~5 au to 2 au when accretion rate varies from $10^{-7} M_{\odot}/yr$ to $10^{-8} M_{\odot}/yr$, respectively. This is the basis for the assertion that in a solar-like protoplanetary disk, we have a good idea of where the snowline should have been as a function of time.

Of course, radial drift complicates the simple picture described above. Morbidelli et al. (2016) point out that decoupling of radial drift velocities between dust and gas can result in a decoupling of “wet” and “dry” gas from the position of the snowline. The rates of icy particle drift inwards towards the star compared with the rates of planetesimal formation then become critical for understanding the compliments of water in protoplanetary materials, as described above in connection with Fig. 1. Morbidelli et al. (2016) suggested that the apparent position of the water snowline recorded by planetesimal formation may be a “fossilized” snowline imprinted on the disk by the gap opened up by a growing Jupiter. Therefore, the rough correspondence between the position of the shift in oxygen fugacity in Fig. 1 and the
natural position of the water snowline for disk mass accretion rates of between $10^{-7} M_{\odot}/yr$ and $10^{-8} M_{\odot}/yr$ could be a manifestation of a locking in of the position of the snowline at these accretion rates by opening of a gap in the vicinity of the snowline by proto-Jupiter.

### 3 Effect of Internal Evolution on Water

During the early evolution of the terrestrial planets and icy bodies, four main sources of energy may have contributed to their internal thermal budget (e.g. Tobie et al. 2013): (i) Tidal heating associated with despinning, (ii) radiogenic heating by the decay of short-lived radioactive isotopes, (iii) impact heating and (iv) conversion of potential energy into heat via viscous dissipation during core-mantle separation. These four sources of energy likely concern different growth stages or reservoirs and hence will be separated in time and space. For planetesimals, only short-lived radiogenic (mainly $^{26}$Al) heating is relevant (Yoshino et al. 2003; Rubie et al. 2015b). This source of energy was for several million times larger than the radiogenic heating rate due to the decay of long-lived isotopes currently expected on planets and moons. However, this huge heat production decayed rapidly. On the other hand on planetary embryos and fully accreted terrestrial planets, the release of potential energy and impact heating have to be taken into account in addition to radiogenic heating. The impact heating strongly depends on the impact velocity, on the masses of the target and the impactor as well as the way kinetic energy is converted into heat (Monteux et al. 2007). Hence the size of the growing body that governs the minimum impact velocity (i.e. the escape velocity) is a key factor. For large growing bodies, large impact events may generate significant melting within the target and trigger metal-rock or ice-rock separation potentially leading to runaway differentiation processes (Golabek et al. 2009; Ricard et al. 2009). Finally on icy moons long-lived radiogenic heating and tidal dissipation are the most important energy sources, while only the largest impact heating events can play a role (Monteux et al. 2014).

The composition of the building blocks from which planetary objects form as well as the early thermal budget should control the amount of water that can be stored in planetary embryos. To get a first order estimation of the amount of hydrous silicates and leftover rock-ice mixture available at the end of the thermo-mechanical evolution of planetesimals, we perform numerical models in 2D infinite cylinder geometry using the finite-difference marker-in-cell code I2ELVIS (Gerya and Yuen 2007; Golabek et al. 2014; Lichtenberg et al. 2016). The physical parameters and the setup employed are identical to those used in our previous studies (see Golabek et al. 2014; Lichtenberg et al. 2016).

Based on current planetesimal formation models (e.g. Johansen et al. 2015; Simon et al. 2016), we model planetesimals with radii $R_p$ ranging between 25 and 230 km and consider instantaneous formation times $t_{\text{form}}$ ranging from 0 to 3.5 Myr after formation of Ca-Al-rich inclusions (CAIs). We study here the evolution of planetesimals that formed across the snowline and were scattered later into the inner solar system (Raymond and Izidoro 2017). Thus we assume for each planetesimal a starting temperature equal to the temperature of surrounding space $T_{\text{space}} = 150$ K. Since for small objects the release of potential energy related to accretion is small (Schubert et al. 1986), the resulting temperature increase is negligible and is not considered here. Thus only short-lived radiogenic heating is relevant. For simplicity we assume that the planetesimals do not contain metal because both formation of hydrated silicates and their dehydration occur at lower temperatures than the melting temperature of the eutectic Fe-FeS mixture (1243 K) at low pressures (Brett and Bell 1969),
so core formation due to percolation is expected to start only after the breakdown of hy-
drated silicates. Therefore we assume that at the start of the model each planetesimal is
composed only of a rock-ice mixture. Here we assume that at T = 273 K the ice melts
and hydrated silicates can form. Under the low-pressure conditions inside a planetesimal the
most temperature-resistant hydrated silicates (amphibolites) break down at T ≈ 1223 K and
dehydrate (Fu and Elkins-Tanton 2014). However it should be taken into account that by
considering the breakdown of amphibolites the current models give only an upper limit for
the amount of hydrated silicates present, since for example more abundant hydrated silicates
like serpentine-phyllosilicates break down at significantly lower temperatures (573-673 K)
(Nakamura 2006; Nakato et al. 2008).

The following additional assumptions and simplifications are used in the numerical
model: We ignore here reaction kinetics and assume that both the formation of hydrated
silicates and their breakdown occur instantaneously at the given temperatures. Also we ig-
nore possible water delivery or loss processes related to impact events. Additionally we do
not consider pore water convection (Young et al. 2003) suggested for smaller planetesimals
(R<sub>P</sub> < 40–60 km) and assume that after silicate dehydration the water previously incorpo-
rated into hydrated silicates is completely lost to space. Also it should be kept in mind that
even dry silicates still contain small amounts of water (see e.g. Fu and Elkins-Tanton 2014),
thus even planetesimals losing all hydrated silicates still contain some water, which is not
considered in our numerical model.

The time-dependent amount of hydrated silicate material in each modeled planetesimal is tracked on the markers of the numerical model. The final 2D area covered by (a) rock-ice mixture and (b) hydrated silicates is converted for representation purposes analytically into a 3D volume (see Fig. 8a, b).

The general results show that the amount of both remaining rock–ice mixture and hy-
drated silicates is strongly dependent on the formation time tstart since this parameter con-
trols the abundance of short-lived ²⁶Al (t₁/₂ = 0.716 Myr), while the radius has only a mod-
est effect (see Fig. 8a, b). As expected late formed objects contain orders of magnitude more
rock–ice mixture and hydrated silicates than early-formed planetesimals. Since formation of
hydrated silicates requires some heating, more of the interior contains hydrated silicates in
late-formed, but large planetesimals, while a primordial rock–ice mixture is best preserved
in the smallest objects experiencing fast cooling. However the results also show that less
hydrated silicates form in planetesimals that accreted very late because the reduced amount
of ²⁶Al radiogenic heating is insufficient to allow for the formation of hydrated silicates
throughout the bulk of the planetesimal interior. This leads to a “sweet spot” for the forma-
tion of hydrated silicates, which is narrower for small planetesimals with R<sub>P</sub> < 80 km since
these objects cool more efficiently than large planetesimals (see Fig. 8a).

All models start with a rock–ice mixture. Because the temperature at the center of the
planetesimal reaches the maximum value, formation of hydrated silicates starts there and
over time the formation front propagates outwards. For large objects (R<sub>P</sub> ≥ 100 km) that
formed early enough (t<sub>form</sub> < 0.75 Myr) the center is also the location where dehydration
starts during later evolution and a dehydration front propagates from there towards the sur-
face. For planetesimals with sufficient ²⁶Al heating both the formation and later the dehy-
dration front are able to propagate throughout the entire planetesimal interior and within
~10⁶ years the entire planetesimal will lose all hydrated silicates (see Fig. 9). On smaller
objects full dehydration occurs only on objects that formed earlier than t<sub>form</sub> < 0.75 Myr,
while in the smallest planetesimals considered (R<sub>P</sub> ≤ 50 km) complete loss of hydrated
silicates is not achieved even inside planetesimals formed at the same time as CAIs (see
Fig. 8a).
<figure><img src="image_409.png" /></figure>

<figure><img src="image_410.png" /><figcaption>**Fig. 8** (a) Final volume of hydrated silicates $V_{\text{hydr}}$ and (b) final volume of rock-ice mixture $V_{\text{rock-ice}}$. Both volumes are scaled by the total volume of the specific planetesimal $V_p$ for models with different planetesimal radii $R_p$ and formation times $t_{\text{form}}$ assuming an initial temperature $T_{\text{space}} = 150$ K. Squares represent the numerical model results. Colours of the squares stand for various final internal states: No hydrous silicates left (red), outer hydrated shell (orange), inner hydrated shell (yellow) and hydrated interior (green). For clarity three of the possible final outcomes are sketched on top of subfigure (a) with dehydrated silicates (white), hydrated silicates (black) and rock-ice mixture (grey).</figcaption></figure>
<figure><img src="image_411.png" /><figcaption>**Fig. 9** Time-dependent fraction of hydrated silicates inside planetesimals with $t_{\text{form}} = 1$ Myr and $R_P = 25–230$ km.</figcaption></figure>

Large, but later-formed planetesimals ($0.75 \text{ Ma} \le t_{\text{form}} < 1.25 \text{ Myr}$) or early-formed, but small objects exhibit a different final outcome, namely an outer hydrated shell, where the formation front reached the object's surface, while the propagation of the dehydration front ceased at depth due to insufficient $^{26}$Al radiogenic heating. Inside large planetesimals that formed even later ($1.25 \text{ Ma} \le t_{\text{form}} < 1.75 \text{ Myr}$) neither the formation nor the dehydration front reach the surface. Under these circumstances the surface is still composed of a rock-ice mixture while the deep interior is dehydrated and an annulus at intermediate depth contains hydrated silicates. For late-formed objects (large objects with $t_{\text{form}} \ge 1.75 \text{ Myr}$ or for smaller planetesimals with $t_{\text{form}} \sim 1.5 \text{ Myr}$ after CAI formation) temperatures are never sufficiently high as to start dehydration and these objects display a hydrated interior and a rock-ice mixture closer to the surface (see Fig. 8a, b). This is in agreement with observations indicating that some asteroids contain even at present-day significant amounts of volatiles (Hsieh and Jewitt 2006; Jewitt et al. 2009; Küppers et al. 2014).

Based on these simple models, we can speculate that large ($R_P > 80$ km), late-formed planetesimals ($t_{\text{form}} > 1.7$ Myr), initially formed across the snowline and implanted into the main asteroid belt or scattered into the inner solar system during the accretion of the gas giants (Raymond and Izidoro 2017), can contribute significant amounts of water in the form of hydrated silicates to growing planetary embryos. This is in general agreement with both the classical model (e.g. Chambers 2013; Raymond et al. 2009) and the Grand Tack scenario (Walsh et al. 2011; Rubie et al. 2015a) suggesting that during the later stages of accretion planetary embryos will accrete material that formed at larger distances from the Sun.

In conclusion, we have shown that the internal evolution and especially the formation time of planetesimals relative to the timescale of radiogenic heating by short-lived $^{26}$Al decay may govern the amount of hydrous silicates and leftover rock-ice mixture available during the late stages of their evolution. In turn, water content may affect the early internal evolution of the planetesimals and in particular the metal-silicate separation processes to be discussed in the next section.
**Fig. 10** Effects of water on the melting of peridotite.

(a) Depression of the peridotite solidus temperature as a function of pressure for bulk water contents of 0, 0.05, 0.1 and 0.3 wt%. (b) Effect of water on melt fraction *F* during partial melting of peridotite shown as isobaric melting curves at 1 GPa for bulk water contents ranging from 0 to 0.3 wt%. Note that the kink at a melt fraction of ~0.3 is caused by the disappearance of clinopyroxene from the mineral assemblage (from Katz et al. 2003). Note that equivalent temperatures for a bulk chondritic composition are lower by ~200°C (Fig. 7 in Asahara et al. 2004). Also the addition of any incompatible elements (e.g. the *C* and *S* that are present in chondritic meteorites) will reduce melting temperatures

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

# 4 Effects of Water on Internal Evolution

The most important process occurring during the early internal evolution of a rocky planetary body is differentiation that results in a silicate mantle and a metallic iron-rich core. The terrestrial planets, Mercury, Venus, Earth and Mars, as well some asteroids, such as Vesta, all underwent core-mantle differentiation. In addition, based on the existence of iron meteorites, many early-formed planetesimals also experienced core-mantle differentiation. Undifferentiated material, as represented by chondritic meteorites, consists of intimate mixtures of metal, sulfide, silicate and oxide grains. The formation of metallic cores in differentiated bodies therefore requires metal and silicate to separate over significant length scales (up to 3000 km in the case of Earth but on the order of 100–500 km in the case of planetesimals). Such segregation was only physically possible at high temperatures at which at least the metal and probably also the silicates were in a molten state (Stevenson 1990; Rubie et al. 2015b). High temperatures were the result of the decay of the short-lived ²⁶Al isotope in planetesimals during the first ~3 Myr of Solar System history and were later the result of high-energy collisions that took place during planetary accretion (e.g. Rubie et al. 2015b; de Vries et al. 2016). Thus planetesimals could only differentiate early, while ²⁶Al was ex- tant, and were largely unaffected by impact-induced heating because of their low masses (see Rubie et al. 2015b, Fig. 10). In general, core-mantle differentiation that occurred as a consequence of impact-induced melting could only happen later after bodies had grown sufficiently in mass. Physical and chemical conditions differed significantly during these respec-
<figure><img src="image_413.png" /><figcaption>**Fig. 11** Schematic section through a terrestrial planet during core formation. A global magma ocean is shown overlying solid crystalline mantle below which lies the metallic proto-core. Metal is delivered by impactors, many of which are likely to have been differentiated into core and mantle. The molten impactor cores partially or completely emulsify in the magma ocean to form small droplets that sink to the base of the magma ocean and collect there to form liquid "iron ponds". These sink through the crystalline mantle as diapirs or possibly are transported through fractures (see Rubie et al. 2015b for a review). Although a possible solid crust is shown, this may not form when an atmosphere is present because its insulating effect keeps the temperature at the surface of a magma ocean above the solidus (e.g., Matsui and Abe 1986)</figcaption></figure>

tive episodes. For example, gravitational forces were much smaller during early planetesimal
differentiation and such bodies likely differentiated in the presence of nebular gases, which
dispersed before the formation of the terrestrial planets.

The presence of water in planetary bodies could have had a significant effect on both
the physical processes of core formation and their chemical consequences. For example,
the presence of even small amounts of water strongly reduces the solidus temperature of
silicates (Fig. 10a). It also increases the melt fraction at a given temperature, although this
effect is largest close to the solidus and becomes much smaller as the degree of melting
increases (Fig. 10b).

**4.1 Effects of Water on Mechanisms of Core Formation**

There are several different mechanisms by which metal and silicate could have segregated (Fig. 11) and the presence of water would have affected each of these differently. Potentially the mechanisms shown in Fig. 11 could all operate during the differentiation of an Earth-mass planet. However, as discussed towards the end of this section, core-mantle differentiation on small bodies (planetesimals) likely required significant melting (with a melt fraction of at least 30–50%) as the result of heat produced by the decay of $^{26}$Al (see also Scheinberg et al. 2015).

We first review the possible metal-silicate segregation mechanisms and how these would
be affected by the presence of water and then discuss chemical consequences of water during
core-mantle differentiation.

**Percolation of Liquid Metal Through a Polycrystalline or Partially-molten Silicate Matrix**

Liquid iron-rich metal or sulfide, (Fe, Ni)S, can percolate through a solid polycrystalline silicate matrix by porous flow provided the liquid phase forms an interconnected network along grain boundaries or grain edges rather than isolated melt pockets. Interconnection depends upon the dihedral angle, $\theta$, which is the angle between two solid-liquid boundaries that are intersected by a solid-solid boundary at a triple junction (Fig. 12; von
<figure><img src="image_414.png" /><figcaption>**Fig. 12** The dihedral angle $\theta$ is the angle between two melt-crystal interfaces where they are intersected by a crystal-crystal grain boundary at a triple junction. When $\theta < 60^{\circ}$, the melt forms an interconnected network along grain edges (*bottom left*), whereas when $\theta > 60^{\circ}$ there is no interconnectivity (unless the melt fraction is high) (from Stevenson 1990)</figcaption></figure>

Bargen and Waff 1986; Stevenson 1990; Rubie et al. 2015b). When $\theta$ is less than 60°, the liquid metal phase is interconnected at any melt fraction and can percolate efficiently through the crystalline matrix. When $\theta$ exceeds 60°, interconnection only occurs when a critical melt fraction (e.g. a few volume % or more) is exceeded, which means that metal-silicate segregation cannot occur efficiently. Experimental studies of both natural and synthetic systems relevant to core formation have shown that, in general, $\theta$ lies in the range 80–120°, at least up to pressures of ~25 GPa (Rubie et al. 2015b and references therein). Therefore percolation is unlikely to be a major mechanism of metal-silicate segregation in the terrestrial planets. However, at pressures <2–3 GPa, the dihedral angle is less than 60° if the liquid Fe contains significant concentrations of oxygen and sulfur (Terasaki et al. 2005, 2008). It is therefore possible that under oxidizing conditions percolation was an important core formation mechanism in planetesimals during heating by $^{26}$Al decay, especially when the metal contains sulfur (which reduces the melting temperature of the metal relative to that of the silicate).

As discussed in the previous section, dehydration of silicates is likely to leave only a few hundred ppm of water in planetesimals that have been heated by the decay of $^{26}$Al. However, such H$_2$O concentrations would be sufficient to lower the silicate solidus temperature by 100–150 K (Fig. 10a), which means that the percolation of liquid metal or sulfide through partially molten silicates, rather than crystalline silicates, needs to be considered. Intuitively, it might be expected that a few percent silicate liquid in the system would enhance the ability of liquid metal to percolate. However, the opposite is the case and the presence of silicate melt actually reduces the interconnectivity of metallic liquids (Holzheid et al. 2000; Rushmer and Petford 2011; Holzheid 2013; Cerantola et al. 2015; Todd et al. 2016). This is the result of the high liquid metal-liquid silicate interfacial energy that causes metallic blobs to adopt a spherical geometry so that they become isolated in silicate melt. Consequently the metallic blobs become trapped in interstices in the partial molten aggregate. Thus, once silicate liquid forms by partial melting, as temperature increases, any interconnectivity of liquid metal or sulfide is destroyed. As the temperature rises further, the silicate melt fraction will eventually become large enough to enable the metallic droplets to become mobile and to segregate gravitationally. The value of the silicate melt fraction at which this happens is poorly known but is likely to be in the range 30–50% (Stevenson 1990; Minarik et al. 1996; Holzheid et al. 2000; Costa et al. 2009; Solomatov 2015). At this degree of partial
melting a magma ocean forms and the segregation mechanism changes to that described below. Thus the presence of water is likely to significantly delay metal-silicate segregation in planetesimals by requiring a higher temperature to start than under anhydrous conditions.

**Segregation of Liquid Metal in a Silicate Magma Ocean** As described above, when the fraction of silicate melt exceeds a certain value (probably 30–50%), a magma ocean state is achieved. When the metallic cores of impacting bodies sink through a magma ocean by gravitational settling, they tend to break up into smaller masses, eventually forming small droplets through emulsification. The extent to which emulsification occurs is important because it determines the fraction of accreted metal that equilibrates chemically with silicate liquid of the magma ocean (e.g., Nimmo et al. 2010; Rudge et al. 2010; Rubie et al. 2015a). When emulsification and equilibration fail to occur, the impactor's core merges directly with the core of the target body by “core merging”. The efficiency of emulsification is uncertain, especially for the cores of giant impactors (Rubie et al. 2003; Dahl and Stevenson 2010; Samuel 2012; Deguen et al. 2014; Wacheul et al. 2014; Kendall and Melosh 2016; Landeau et al. 2016). The process is extremely difficult to study numerically because the length scales involved range from 100's km (or more) to a few mm. In the case of the Earth, it has been estimated that emulsified droplets of liquid Fe have a stable diameter of ~1 cm and sink with a velocity of ~0.5 m/s (Stevenson 1990; Rubie et al. 2003). Such parameters are a result of molten peridotite having a very low viscosity (Liebske et al. 2005)—which also results in magma oceans convecting vigorously in a turbulent regime. Because peridotite liquid is already highly depolymerized, the addition of water is unlikely to have a significant effect on its viscosity. In addition, as shown in Fig. 10b, the presence of water will have only a small effect on the temperature of magma ocean formation when the latter is defined to have a silicate melt fraction ≥0.3.

**Descent of Diapirs of Liquid Metal Through Crystalline Mantle** It has been proposed that segregated liquid metal accumulates as a pool at the bottom of a magma ocean during the differentiation of terrestrial planets. If the base of the magma ocean is separated from the proto-core by a crystalline mantle, it is likely that the iron melt pool sinks through the latter as diapirs, perhaps 1–10 km in diameter or larger (Fig. 11; e.g. Stevenson 1990; Karato and Murthy 1997; Samuel et al. 2010). Water dissolved in crystalline upper mantle rocks has a strong weakening effect on rheology and reduces the viscosity of solid-state flow by up to 2–3 orders of magnitude. For example, Hirth and Kohlstedt (1996) show that 810 ± 490 H/10⁶ Si dissolved in olivine reduces the viscosity by a factor of 500 ± 300 compared with that of dry olivine. The descent velocity of metal diapirs can be described by Stokes' Law. Since the Stokes' sinking velocity is inversely proportional to mantle viscosity, the effect of water is large and could enhance the descent velocity by up to 2–3 orders of magnitude, at least in upper mantle peridotite.

Numerical models of diapirs sinking through the solid mantle of Earth-sized bodies indicate that a wet rheology favors diapir sinking, while a dry rheology leads to preferential formation of iron-filled dikes (Golabek et al. 2009). It was shown both analytically and numerically that this changes the heat partitioning between iron and silicates (Ke and Solomatov 2009; Golabek et al. 2009). Whereas the sinking of diapirs results in heat preferentially partitioning into the silicate mantle, the propagation of iron dikes results in the preferential heating of the iron. Thus the presence of water in the crystalline mantle could affect the thermal state of a terrestrial planet after core formation.

In summary, the effect of water on the temperature and timing of magma ocean formation in accreting planetary bodies, as required for metal-silicate segregation, is very small.
In relatively oxidized planetesimals, the presence of water may eliminate the possibility of metal-silicate segregation by percolation, as proposed by Terasaki et al. (2008), so that segregation only starts at a magma ocean stage after mantle material is at least 30–50% molten. If the descent of diapirs of liquid metal through crystalline mantle occurs, water dissolved in the mantle material would greatly enhance the rate of this process. On the other hand, an absence of water could change the descent mode from sinking diapirs to diking, thus affecting the heat partitioning between iron and silicates and modifying the thermal state of the early terrestrial planet.

## 4.2 Effects of Water on the Chemistry of Core Formation

The main geochemical effect of metal-silicate segregation to form the core and mantle of a planetary body is that siderophile (metal-loving) elements, such as Ni, Co, V, Cr, W, Mo, etc., partition into liquid metal and are thus transported from the mantle into the core. Consequently, all siderophile elements are variably depleted in Earth's mantle, relative to solar relative abundances, by factors of up to almost 1000. The degree of depletion depends on the extent to which an element partitions into metal. This is described by the metal-silicate partition coefficient, $D_M^{\text{met-sil}}$, which for element M is defined simply as:

$$D_M^{\text{met-sil}} = \frac{C_M^{\text{met}}}{C_M^{\text{sil}}}$$

where $C_M^{\text{met}}$ and $C_M^{\text{sil}}$ are the molar (or wt%) concentrations of M in the metal and silicate, respectively. Most metal-silicate partition coefficients for siderophile elements depend on pressure (P) and/or temperature (T). In addition, and of great importance for the present discussion, they depend strongly on oxygen fugacity $f_{O_2}$, as shown by the redox reaction:

$$ \underset{\text{metal}}{M + (n/4)O_2 = MO_{n/2}} $$
silicate liquid

where n is the valence of element M when incorporated in silicate liquid as an oxide component. When $f_{O_2}$ is low, this reaction shifts to the left thus increasing $D_M^{\text{met-sil}}$, whereas when $f_{O_2}$ is high the reaction shifts to the right thus decreasing $D_M^{\text{met-sil}}$.

When modeling core formation, the extent to which Earth's mantle has been depleted in siderophile elements enables estimates to be made of the conditions (P, T, $f_{O_2}$) of metal-silicate equilibration during core formation (e.g. Wade and Wood 2005; Walter and Cottrell 2013; Rubie et al. 2011, 2015a, 2016). Traditionally, it has been assumed that core formation was a single-stage event, thus enabling unique values for P, T and $f_{O_2}$ to be determined (Walter and Cottrell 2013; see Table 3 in Rubie et al. 2015b). However, “single-stage” core formation is clearly unrealistic because metal is delivered to a growing planet during the entire accretion process, which in the case of Earth likely lasted for up to 100 Myr (e.g. Jacobson et al. 2014). Thus, in recent years more realistic core formation models have been formulated. Firstly, in models of “continuous” core formation (Wade and Wood 2005; Wood et al. 2006), material is accreted to Earth in small (1%) increments with the metal component being equilibrated with the silicate mantle at the base of a magma ocean. In a second class of models, core formation is multi-stage and the accretion of mostly-differentiated embryos and planetesimals to growing planets is linked to episodes of metal-silicate equilibration and core formation (Rubie et al. 2011, 2015a, 2016). In both types of models, magma ocean depth increases as planets grow
so that the pressure and temperature of metal-silicate equilibration both increase as accretion proceeds. In addition, the oxygen fugacity at which metal and silicate equilibrate has to increase by 2–3 orders of magnitude during accretion in order to reproduce the concentrations of siderophile elements in Earth’s mantle (Wade and Wood 2005; Wood et al. 2006; Rubie et al. 2011). There are at least three mechanisms that can result in, or contribute to, this increase of $f_{O_2}$, two of which are intimately related to water. Note that oxygen fugacity is determined from the compositions of co-existing metal and silicate relative to that of the iron-wüstite (IW) buffer, as follows:

$$ \Delta \text{IW} = 2 \log \left[ \frac{a_{\text{FeO}}}{a_{\text{Fe}}} \right]. $$

Here $\Delta \text{IW}$ is the deviation of $f_{O_2}$ from the IW buffer in log units, and $a_{\text{FeO}}$ and $a_{\text{Fe}}$ are the activities of FeO and Fe in the silicate and metal, respectively. Under low oxygen fugacities (e.g. $\Delta \text{IW} = -4$) the FeO content of the silicate is low (e.g. <1 wt%), whereas at higher oxygen fugacities (e.g. $\Delta \text{IW} = -2$) the FeO content becomes similar to that of Earth’s mantle (~8 wt%).

**Partitioning of Silicon into the Core** Under reducing conditions and especially at high temperature, Si becomes siderophile (e.g. Mann et al. 2009) and this element partitions into liquid core-forming metal by the reaction

$$ \begin{array}{l} \text{SiO}_2 + 2\text{Fe} \rightarrow \text{Si} + 2\text{FeO}. \\ \text{mantle core} \quad \text{core mantle} \end{array} $$

According to this reaction, for every mole of Si that partitions into the core, two moles of FeO are added to the mantle. This reaction is therefore very efficient at increasing the oxygen fugacity, as required by the core formation models of Wade and Wood (2005) and Rubie et al. (2011, 2015b).

**Oxidation of Metal by Accreted Water** At the end of core formation, the oxygen fugacity of Earth’s mantle was low (~IW – 2) whereas it is currently several orders of magnitude higher (close to the FMQ buffer). In the upper mantle, the currently high oxygen fugacity is the result of a high $Fe^{3+}$ concentration in some of the constituent minerals (O’Neill et al. 1993). The increase in $f_{O_2}$ after core formation may have been caused by the disproportionation reaction

$$ \text{Fe}^{2+} \rightarrow \text{Fe}^{3+} + \text{Fe metal} $$

which occurs as a result of the crystallization of the lower-mantle mineral bridgmanite (silicate perovskite) (Frost et al. 2004). Water can also potentially cause oxidization in planetary interiors, for example by the reaction:

$$ \begin{array}{l} \text{H}_2\text{O} + \text{Fe} \rightarrow \text{H}_2 + \text{FeO}, \\ \text{metal} \quad \text{silicate} \end{array} $$

for which the loss of hydrogen is required, or, in the absence of metal, by the dissociation of H₂O into oxygen and hydrogen, with the latter being lost. Sharp et al. (2013) proposed that such reactions were the cause of the increase in $f_{O_2}$ of Earth’s mantle following core formation.
All oxidation reactions that are caused by water require that hydrogen is lost from the system. In the absence of metal, this may be problematic because hydrogen is a strong reducing agent. Therefore, if oxidation is taking place deep in the Earth's interior, the released hydrogen has to migrate to the Earth's surface without causing any reduction of the material that it comes into contact with, which seems very unlikely. If metal is present, the hydrogen may be removed by dissolving in the metal, so that, during core formation, it would be transported to the core of the planet. It is uncertain how efficient this process would be because the partitioning of hydrogen between metal and silicate under core formation conditions is poorly known (Okuchi 1997). The oxidation of metal by water during core formation has been included in the accretion/core formation model of Rubie et al. (2015b) assuming that all released hydrogen dissolves in the metal. It was thus shown that the role of accreted water in oxidizing the mantle during core formation is insignificant and that all oxidation (up to the end of core formation) can be explained by the dissolution of Si into the core (as discussed above) combined with the effects of the accretion of oxidized bodies due to an oxidation gradient in the solar nebula (Fig. 1).

# 5 Differentiation and Evolution of Water-Rich Objects: From Icy Moons to Ocean-Planets

Most of the physical processes that are discussed to occur during the early stages of planetary formation also occur in the outer part of the Solar System during the formation of icy moons. In these bodies the mantle is composed mainly of water ice while their core is essentially made of rocky material. While the variety of internal structures among the icy moons is still debated, the essential role of the water is widely accepted. These objects are discussed in this section.

## 5.1 Variety of Internal Structures of Icy Moons

Differences in composition and internal structure exist among the major icy satellites of Jupiter and Saturn, suggesting distinct accretion and differentiation histories (e.g., Kirk and Stevenson 1987; Mueller and McKinnon 1988; Mosqueira and Estrada 2003; Barr and Canup 2008; Monteux et al. 2014). In this section we will focus particularly on the largest icy moons having a planetary scale, namely Jupiter's moons Ganymede and Callisto, and Saturn's moon Titan. The internal structure of these bodies has been constrained from the measurements of the gravity coefficients from spacecraft radio tracking during close flybys (e.g. Galileo, Cassini) (Anderson et al. 1996, 2001; Iess et al. 2010, 2012). The degree-2 gravity coefficients can then be used to infer the moment of inertia of the satellite, thus providing constraints on the density structure. The normalized moment of inertia about the spin axis, $C/(M_T R_T^2)$, is commonly calculated from the fluid Love number using the Radau–Darwin approximation. This formulation implies that the body's equipotential corresponds to the equilibrium shape due to fluid dynamic flattening. In this approximation, the degree-two gravity coefficients, $C_{20}$ and $C_{22}$, should be equal to $-10/3$ and are proportional to the fluid Love number.

In the case of Titan, the ratio between the observed $C_{20}^{\text{obs}}$ and $C_{22}^{\text{obs}}$ is close to $-10/3$ (Iess et al. 2010), suggesting that Titan is relatively close to hydrostatic equilibrium. We can therefore use the degree-two gravity coefficients to first order to estimate the fluid Love number and the corresponding moment of inertia factor. Depending on whether we use the $C_{20}$ or $C_{22}$
coefficients as a reference to determine the fluid Love number, the normalized moment of inertia varies between 0.335 and 0.342. For Ganymede and Callisto, the two coefficients have not been determined independently, and the -10/3 relationship has been prescribed when inverting the degree-two gravity coefficients. Using the degree-two coefficients, and applying the Radau-Darwin relationship, the moment of inertia factors of Ganymede and Callisto are estimated to be 0.31 and 0.355 (assuming hydrostatic equilibrium), respectively (Anderson et al. 1996, 2001). Titan's value is therefore intermediate between those of Ganymede and Callisto, suggesting a modest increase of density toward the center. The high moment of inertia factor obtained by Galileo gravity measurements ($C/(MR^2) = 0.355$) (Anderson et al. 2001) suggests an incomplete ice-rock separation within Callisto. In contrast, Ganymede has a much smaller moment of inertia ($C/(MR^2) = 0.31$) (Anderson et al. 2001) and shows signs of past endogenic activity (Pappalardo et al. 2004). Hence, Ganymede is potentially fully differentiated with an icy upper mantle, a rocky lower mantle and a metallic core, which is the origin of a relatively intense intrinsic magnetic field (Kivelson et al. 1998). With similar size and mass, Titan may be an intermediate case between Callisto and Ganymede. Based on the estimation of its moment of inertia factor, Titan's interior is more differentiated than Callisto's but probably much less than Ganymede's interior. Like Callisto, Titan might still possess a layer of ice-rock mixture between a rocky core and an outer ice-rich mantle, unless the rocky core is composed mostly of hydrated minerals (Sohl et al. 2010; Castillo-Rogez and Lunine 2010). A dense atmosphere is currently present at its surface but no intrinsic magnetic field has yet been detected on Titan. However the presence of a metallic core is plausible within Titan if enough iron was delivered during its accretion (Grasset et al. 2000). The fact that the interior of Callisto and possibly Titan may still contain a layer of ice-rock mixture suggests that the satellites have avoided significant melting during accretion and subsequent evolution.

We now discuss the apparently partially-differentiated state of Callisto and Titan, as suggested by their elevated moment of inertia factors (e.g. Anderson et al. 2001; Iess et al. 2010; Gao and Stevenson 2013). In the absence of lithospheric strength and shell thickness variations, these satellites would adopt a perfectly hydrostatic shape, corresponding to a degree-two gravity field. On Titan, the existence of a non-negligible degree-three in the gravity field as well as significant topography suggest that non-hydrostatic effects may significantly affect the estimation of the moment of inertia factor (Iess et al. 2010; Gao and Stevenson 2013; Baland et al. 2014) and that the moment of inertia factor may be significantly smaller than the value estimated from the Radau-Darwin approximation. On Callisto, similar non-hydrostatic contributions originating in the lithosphere, which could not be tested with the Galileo spacecraft data, may also affect the estimation of its moment of inertia factor (McKinnon 1997; Gao and Stevenson 2013). On these two moons, the hydrostatic dynamical flattening is relatively small as they orbit relatively far from their respective planets, and therefore the non-hydrostatic contributions need to be correctly estimated in order to accurately infer the moment of inertia factor and the density profiles of their interior. Future measurements by the JUICE mission (Grasset et al. 2013) will allow us to test how far Ganymede and Callisto are from hydrostatic equilibrium.

## 5.2 Differentiation Process in Water-Dominated Interiors

Several scenarios have been proposed to explain the wide range of extents of differentiation among the icy moons that resulted from radioactive and accretional heating. O'Rourke and Stevenson (2014) showed that although rock-ice separation may be delayed by double-diffusive convection in the ice-rock interior, ice melting due to progressive radiogenic heating and subsequent differentiation cannot be prevented. Previous studies showed that it was
possible to avoid melting if the accumulation of accretion energy was inefficient, i.e. if the energy was radiated away at a rate comparable to the accretion rate (e.g., Schubert et al. 1981; Squyres et al. 1988; Kossacki and Leliwa-Kopystynski 1993; Coradini et al. 1995; Grasset and Sotin 1996). Based on these models, the accretion timescales should be longer than 1 Myr to avoid significant melting and hence differentiation of Callisto, while an accretion timescale as short as $10^3-10^4$ yr may be possible for Ganymede. However, these timescales are dependent on the way heat deposition and cooling are treated. Barr and Canup (2010) proposed that the Ganymede-Callisto dichotomy can be explained through differences in the energy received during the Late Heavy Bombardment. Impacts would have been sufficiently energetic on Ganymede to lead to a complete ice-rock separation, but not on Callisto. More recently, Monteux et al. (2014) have developed a 3D numerical model that characterizes the thermal evolution of a satellite growing by multiple impacts, simulating the satellite growth and thermal evolution for a body radius ranging from 100 to 2000 km. Their results indicate that a satellite exceeding 2000 km in radius may accrete without experiencing significant melting only in case (i) its accretion history is dominated by small impactors (radius less than a few kilometers) and (ii) the conversion of impact energy into heat is unrealistically inefficient. In this context, they underlined that global melting for large bodies like Titan or Callisto is difficult to avoid.

Once the melting of the icy phase is initiated, ice-rock separation can proceed. This separation will lead to the accumulation of dense blocks of rock overlying a undifferentiated core consisting of a mixture of rock and ice (e.g., Friedson and Stevenson 1983; Kirk and Stevenson 1987). The dense layer of accumulated rock is gravitationally unstable, which leads to large-scale ice-rock separation. Friedson and Stevenson (1983) suggested that an increase in the silicate volume fraction from 44% for Callisto to 55% for Ganymede could have induced an increase in the internal viscosity of between 1 and 2 orders of magnitude. This increase in the internal viscosity of Ganymede might have led to a failure of its solid-state convection self-regulation. This failure would then have limited heat loss, favored partial melting of the ice and enhanced a catastrophic runaway differentiation (Kirk and Stevenson 1987). Similarly, Nagel et al. (2004) suggested that if Callisto accreted from a mixture of rock and ice and the average size of the rocks was of the order of meters to tens of meters, then the viscosity would range between $10^{13}$ and $10^{17}$ Pa·s. For such viscosities, Callisto would have experienced a gradual, but incomplete unmixing of ice and rocks. This mantle/core separation potentially involves a large amount of heating due to the viscous dissipation of the potential energy (Flasar and Birch 1973). Contrary to terrestrial planets (Golabek et al. 2009; Ricard et al. 2009), this runaway separation has not yet been modeled.

## 5.3 Metal-Rock Separation

If a sufficient amount of native iron is present within the icy moon, an iron core may eventually form during or after the ice-rock separation. In this case, the deep interior would be comprised of a probably liquid iron core, and a convecting silicate mantle. Gravity measurements performed by the Radio Science System onboard the Cassini spacecraft suggest that the deep interior of Titan is not fully differentiated (Sotin et al. 2009a; Iess et al. 2010). However, the presence of an iron core cannot be definitely ruled out, and additional measurements are required to determine whether Titan possesses an iron core like Ganymede. On Ganymede, a second differentiation event likely occurred during or after the ice-rock separation since enough heat had been stored in the metal-silicate mixture leading to the formation of an iron core. The major observable feature of this metal-silicate separation is the presence of an intrinsic magnetic field generated in its deep interior (Kivelson
<figure><img src="image_415.png" /><figcaption>**Fig. 13** Schematic representation of the structure of an icy moon (a) just after accretion and after differentiation leading to (b) a Callisto-like structure, (c) a Titan-like structure or (d) a Ganymede-like structure</figcaption></figure>

et al. 1996). While the ice-rock separation has been extensively studied, the modalities of metal-rock separation on icy moons and, as a consequence, the modalities of the generation of intrinsic magnetic fields remain poorly constrained. Another major difficulty is to constrain the composition and size of the metallic core. Depending on whether the core is pure iron or a sulfur-iron alloy, it could represent 1 to 30% of the total mass of the icy moon (Anderson et al. 1996).

Although the interior of Ganymede is currently fully differentiated, the silicate rock and the iron likely coexisted and reacted with H₂O (as ice or water) before metal/silicate separation occurred. In this context, it is difficult to envision the formation of a pure iron core because oxidation of the iron would be unavoidable. This led Crary and Bagenal (1998) to suggest that Fe₃O₄-magnetite could maintain a strong remanent magnetic field within the rocky material of Ganymede. Scott et al. (2002) have shown that the presence of water could have enhanced the formation of low-density hydrated silicates and the alloying of iron with sulfur would have resulted in FeS-dominated cores. They performed high-pressure experiments that combined chondritic chemistry with water at a pressure of 1.5 GPa, temperatures between 300°C and 800°C and a range of oxygen fugacities. Their results indicate that sulfur may be a significant constituent of Ganymede's core but they did not observe the presence of either iron oxide species or metallic iron. Hence they concluded that magnetite is not a stable phase in the chondrite-water system in contrast with the suggestion of Crary and Bagenal (1998).

Hence, Ganymede's core probably contains not just pure iron but also lighter elements. As mentioned above a candidate for the dominant core-forming assemblage is an alloy of iron (Fe) and sulfur (S) (Hauck et al. 2006). Indeed sulfur is an element that was abundant in the proto-solar system (Hillgren et al. 2000). Within Ganymede, the core sulfur concentrations could range between 12 and 20 wt%. Under the pressure conditions of Ganymede's core (≈6–10 GPa) a mixture of iron and sulfur forms a binary, eutectic sys-
tem. A major consequence of adding sulfur to liquid iron is a depression of the liquidus (e.g., Usselman 1975; Fei et al. 1997). Upon freezing, nearly pure iron (Fe) (Li et al. 2001; Kamada et al. 2010) crystallizes on the Fe-rich side of the eutectic and iron sulfide (FeS) on the Fe-poor side. Various experiments have revealed that the eutectic composition of Fe-FeS becomes increasingly Fe-rich with increasing pressure (e.g., Chudinovskikh and Boehler 2007; Fei et al. 1997, 2000; Morard et al. 2008). Buono and Walker (2011) have suggested low and even negative gradients for the melting temperatures of Fe-rich Fe-FeS alloys in the pressure range of Ganymede's core. If the melting temperature has a lower gradient than the core temperature, iron will start to crystallize at the core-mantle boundary (CMB) and a growing snow zone will form during further cooling (McKinnon 1996; Hauck et al. 2006). However, the unknown oxidation state of the interior during differentiation limits our ability to constrain whether the composition of a sulfur-bearing core is on the Fe- or FeS-rich side of the eutectic composition (e.g., Scott et al. 2002; Breuer et al. 2015). Another important element could be hydrogen. The interaction of native iron with hydrated material could have led to the addition of hydrogen to the iron core resulting in the presence of iron hydride (Ohtani et al. 2005; Shibazaki et al. 2009).

## 5.4 Evidence for Active and Past Aqueous Processes

A range of geophysical evidence indicates that several moons of Jupiter (Europa, Ganymede, Callisto) and Saturn (Titan, Enceladus, Mimas) harbor extensive liquid water reservoirs beneath their cold icy surfaces (Khurana et al. 1998; Kivelson et al. 2000; Iess et al. 2014; Tajeddine et al. 2014; Thomas et al. 2016). Detection of ammoniated phyllosilicates on Ceres (de Sanctis et al. 2015) as well as the young age of some extensional features observed on Pluto (Hammond et al. 2016; Moore et al. 2016) suggest the existence of a subsurface ammonia-rich ocean, at least in the past, inside these two dwarf planets. Based on evolution models, many other icy bodies are also predicted to harbor internal oceans underneath their cold icy shells (Hussmann et al. 2006). These massive oceans constitute the greatest volumes of potentially habitable space in the Solar System. If chemical interactions with the rocky core and the irradiated surface provide nutriments and energy sources in sufficient amounts, the development of primitive life may be envisaged in these water-rich worlds. Several observational clues suggest that complex interactions between the subsurface oceans and their warm rocky interiors are presently occurring.

The most convincing evidence for active aqueous processes is provided by the surprising observations of eruptions of water vapor and salted icy grains at Enceladus' active south pole, indicating connections with a subsurface ocean (Postberg et al. 2009, 2011). The discovery of silicon-rich nano-particles further indicates that hydrothermal interactions are currently occurring in porous warm (100–200 °C) rock and that hydrothermal products are quickly transported to the plume source (Hsu et al. 2015; Sekine et al. 2015). Moreover, analysis of gravity, shape and libration data shows that the ice shell thickness at the south pole is thinner than 5 km, with a globally average value of not more than 21 km, requiring about 25–30 GW of internal power for such a configuration to be thermally stable (Cadek et al. 2016). This provides additional constraints in favor of a warm rocky core, powered by strong tidal dissipation. However, it remains unclear how Enceladus reached such a hyperactive state.

On Europa, although there is no direct evidence yet, the occurrence of hydrothermal processes at the seafloor, at present or in recent geological history (<100 Ma), is also likely. The observation of non-ice colored and hydrated materials in many disrupted areas on the surface of Europa points to exposure of oceanic species at the surface and to the possible
existence of suboceanic hydrothermal systems (Zolotov and Kargel 2009). The likelihood of presently active hydrothermal processes on Europa mostly depends on how much heat can be generated in the interior at present. Similarly to Enceladus and Io, tidal friction may contribute to Europa's internal budget, but in absence of observational constraints, it is still unclear if tidal friction in the rocky mantle is large enough at present to maintain active hydrothermal processes (Sotin et al. 2009b).

Despite the separation of the subsurface ocean from the rocky core by a thick high pressure ice mantle, water-rock interactions at the base of the thick ice mantle and the transport of brine through the ice mantle are still possible on the largest moons like Saturn's moon Titan and Jupiter's moon Ganymede (Vance and Brown 2013; Choblet et al. 2017). It is very likely that water-rock interactions occurred on these large moons in the past, especially during the differentiation stage when a large volume of water percolated throughout rock-dominated layers toward the surface (Tobie et al. 2013) (Fig. 13). The detection of $^{40}$Ar, the decay product of $^{40}$K initially contained in silicate minerals, in Titan's atmosphere (Niemann et al. 2010) as well as the likely presence of salts in the internal ocean of Titan (Béghin et al. 2012; Baland et al. 2014; Mitri et al. 2014), Ganymede (Saur et al. 2015) and Callisto (Khurana et al. 1998) indicate that water-rock interaction occurred during the evolution of these large icy moons. However, it is still unclear for how long these water-rock interactions lasted and what effects they had on the thermo-chemical evolution of their interiors.

## 5.5 Role of Volatile Compounds in the Aqueous Alteration Within Icy Bodies

During the early evolution of Pluto, impact heating, radiogenic heat sources (short- and long-lived) and tidal heating due to the interaction with Charon have likely been significant to induce a complete differentiation of Pluto's interior (Schubert et al. 2010, Sekine et al. 2017). At present-day, because of the rapid decay of short-lived radioisotopes and since the Pluto/Charon system has reached the “dual synchronous” state in which Charon's orbital period, spin period, and Pluto's rotation period are equal, these two major heat sources are probably too weak to maintain a current water ocean within Pluto's interior (Robuchon and Nimmo 2011; Barr and Collins 2015). The surprising discovery of two quasi-circular mounds on Pluto by New Horizons, suspected to be of cryovolcanic origin (Moore et al. 2016), suggests that the interior of the dwarf planet is more active than initially anticipated and that water-volatile-rock interactions may have operated during a long period of time in its interior. The presence of significant amounts of anti-freeze compounds, such as ammonia and methanol, together with gas compounds in the form of clathrate hydrates may explain the occurrence of such activity on Pluto. This hypothesis is consistent with the detection of ammonia-rich and methane-rich areas on Charon (Grundy et al. 2016). Although ammonia has not been identified on Pluto's surface, that is dominated by $N_2$ ices (Grundy et al. 2016), its presence on Charon indicates that it may have played a key role in the chemical evolution of the Pluto-Charon system, as well as in those of other Kuiper belt objects (Hussmann et al. 2006).

For decades, both ammonia and methanol, have been suspected to strongly affect the melting and crystallization process in icy bodies (e.g. Lewis 1971; Kargel 1992). The role of ammonia is expected to be particularly strong during the early stage in the evolution of icy bodies, when the first aqueous melts are generated during the post-accretional differentiation process and at the end of their evolution when the internal ocean becomes strongly concentrated in solutes by progressive crystallization of ice (Hussmann et al. 2006; Malamud and Prialnik 2015). Even in the absence of significant heat sources, the presence of ammonia, as well as methanol (both are present in comparable amount in comets, Cochran
et al. 2015) can sustain active aqueous processes for long periods of time. The recent detection of widespread ammoniated phyllosilicates across Ceres' surface (de Sanctis et al. 2015) may provide evidence for alteration of chondritic materials by N-bearing fluids. However, in absence of reliable experimental data on aqueous alteration by concentrated N-bearing fluids, it is difficult to investigate the conditions under which such ammoniated phyllosilicates may have formed.

# 6 Conclusions and Outlook

Water content and the internal evolution of terrestrial planets and icy bodies are closely linked. Indeed, the internal evolution and especially the formation time of the planetesimals and therefore the extent of radiogenic heating governs the amount of hydrous silicates and leftover rock-ice mixtures available in the late stages of their evolution. In turn, water content may affect the early internal evolution of the planetesimals and metal-silicate separation by favoring the sinking of metal diapirs over percolation processes. On the other hand, the absence of water could change the descent mode on larger terrestrial bodies from sinking diapirs to diking, thus affecting the heat partitioning between iron and silicates and modifying the thermal state of the accreting planet. Moreover, water content can contribute to an increase of $f_{O_2}$ and thus affect the concentration of siderophile elements within the silicate reservoirs of the Solar System objects. Finally, the water content strongly influences the differentiation rate among the icy moons, controls their internal evolution and governs the alteration processes occurring in their deep interiors.

Partitioning of water between the mantle and the surface/atmosphere is also a key factor in regulating both planetary climate and internal dynamics. This partitioning is controlled by release of water by erupting melts, by hydrothermal alterations, mainly serpentinization, trapping some water in the form of hydrated minerals in the crust, and by recycling part of it deep into the mantle where it can be dissolved in nominally-anhydrous minerals. On Earth, where the water cycle is driven by plate tectonics, exactly how water partitions between surface and interior remains uncertain, but it is commonly thought to be of the order of 1. Parameterized models have recently been proposed to address the question of water cycle on other terrestrial planets, with different sizes and water contents (Crowley et al. 2011; Cowan and Abbot 2014; Schaefer and Sasselov 2015). These parameterized approaches enable some important characteristics of the water cycle to be understood, especially how surface gravity and the water volume at the surface influence the interior/surface exchange. However, they cannot describe the complexity of water transport in the mantle and its feedback on internal melting and crustal recycling. Describing correctly the coupled evolution of surface and internal reservoirs of water requires the use of thermo-chemical convection codes in 3D or at least 2D, including exchange of water between the crust and mantle, the convective transport in the mantle, its retroaction on viscosity, melting and associated outgassing.

Finally, with the discoveries of thousands of exoplanets, we are now entering a new era in comparative planetology. Even if the dataset is limited to the mass, radius and orbital characteristics of the planetary systems, it already offers a unique opportunity to test the models of planet formation and evolution. Model mass-radius relationships indicate a great diversity of interior compositions and atmospheric extents for the Super-Earth/Mini-Neptune-planet class (e.g. Howard et al. 2012), suggesting a wide range of volatile contents and compositions. Hence, planetary worlds are probably much more diverse than originally thought, with a wide range of water and other volatile contents.

<figure><img src="image_416.png" /></figure>
**Acknowledgements** We thank L. Elkins-Tanton, A. Morbidelli and two anonymous reviewers for detailed and thoughtful comments that helped to improve the manuscript. D.C.R. was supported by the European Research Council Advanced Grant “ACCRETE” (contract number 290568) and additional support was provided by the German Science Foundation (DFG) Priority Programme SPP1833 “Building a Habitable Earth” (Ru 1323/10-1). J.M. was funded by the Auvergne Fellowship program. This is Laboratory of Excellence Cler-Volc contribution no. 283. G.J.G. thanks Taras Gerya for providing the code I2MART. G.T. acknowledges support from the project ANR OASIS.

## References

J.D. Anderson, R.A. Jacobson, T.P. McElrath, W.B. Moore, G. Schubert, P.C. Thomas, Shape, mean radius, gravity field, and interior structure of Callisto. *Icarus* **153**, 157–161 (2001)

J.D. Anderson, E.L. Lau, W.L. Sjogren, G. Schubert, W.B. Moore, Gravitational constraints on the internal structure of Ganymede. *Nature* **384**, 541–543 (1996)

P.J. Armitage, *Astrophysics of Planet Formation* (Cambridge University Press, Cambridge, 2010), 294 p.

Y. Asahara, T. Kubo, T. Kondo, Phase relations of a carbonaceous chondrite at lower mantle conditions. *Phys. Earth Planet. Inter.* **143–144**, 421–432 (2004)

X-N. Bai, J.M. Stone, Wind-driven accretion in protoplanetary disks. I. Suppression of the magnetorotational instability and launching of the magnetocentrifugal wind. *Astrophys. J.* **769**, 21 (2013)

R.M. Baland et al., Titan's internal structure inferred from its gravity field, shape, and rotation state. *Icarus* **237**, 29–41 (2014)

A.C. Barr, R.M. Canup, Constraints on gas giant satellite formation from the interior states of partially differentiated satellites. *Icarus* **198**, 163–177 (2008)

A.C. Barr, R.M. Canup, Origin of the Ganymede–Callisto dichotomy by impacts during the Late Heavy Bombardment. *Nat. Geosci.* **3**, 164–167 (2010)

A.C. Barr, G.C. Collins, Tectonic activity on Pluto after the Charon-forming impact. *Icarus* **246**, 146–155 (2015)

C. Béghin, O. Randriamboarison, M. Hamelin, E. Karkoschka, C. Sotin, R.C. Whitten, J.-J. Berthelier, R. Grard, F. Simões, Analytic theory of Titan's Schumann resonance: constraints on ionospheric conductivity and buried water ocean. *Icarus* **218**(2), 1028–1042 (2012)

P.A. Bland et al., Why aqueous alteration in asteroids was isochemical: high porosity does not equal permeability. *Earth Planet. Sci. Lett.* **287**, 559–568 (2009)

R. Brett, P.M. Bell, Melting relations in the Fe-rich portion of the system Fe–FeS at 30 kb pressure. *Earth Planet. Sci. Lett.* **6**, 479–482 (1969)

D. Breuer et al., Iron snow, crystal floats, and inner-core growth: modes of core solidification and implications for dynamos in terrestrial planets and moons. *Prog. Earth Planet. Sci.* **2**, 39 (2015)

A.S. Buono, D. Walker, The Fe-rich liquidus in the Fe–FeS system from 1 bar to 10 GPa. *Geochim. Cosmochim. Acta* **75**, 2072–2087 (2011)

O. Cadek, G. Tobie, T. Van Hoolst, M. Massé, G. Choblet, A. Lefèvre, G. Mitri, R-M. Baland, M. Behounkova, O. Bourgeois, A. Trinh, Enceladus's internal ocean and ice shell constrained from Cassini gravity, shape, and libration data. *Geophys. Res. Lett.* **43**(11), 5653–5660 (2016)

J.C. Castillo-Rogez, J.I. Lunine, Evolution of Titan's rocky core constrained by Cassini observations. *Geo-phys. Res. Lett.* **37**, L20205 (2010)

J.C. Castillo-Rogez, E. Young, Origin and evolution of volatile-rich asteroids, in *Early Differentiation and Consequences for Planets*, ed. by L.T. Elkins-Tanton, B. Weiss. Cambridge Planetary Science (2017)

F. Cerantola, N. Walte, D.C. Rubie, Deformation of a crystalline olivine aggregate containing two immiscible liquids: implications for early core-mantle differentiation. *Earth Planet. Sci. Lett.* **417**, 67–77 (2015)

J.E. Chambers, Late-stage planetary accretion including hit-and-run collisions and fragmentation. *Icarus* **224**, 43–56 (2013)

E. Chiang, A.N. Youdin, Forming planetesimals in solar and extrasolar nebula. *Annu. Rev. Earth Planet. Sci.* **38**, 493–522 (2010)

G. Choblet, G. Tobie, C. Sotin, K. Kalousova, O. Grasset, Heat and melt transport in the high-pressure ice mantle of large icy moons. *Icarus* **285**, 252–262 (2017)

L. Chudinovskikh, R. Boehler, Eutectic melting in the system Fe–S to 44 GPa. *Earth Planet. Sci. Lett.* **257**(1), 97–103 (2007)

A.L. Cochran et al., The composition of comets. *Space Sci. Rev.* **197**(1–4), 9–46 (2015)

A. Coradini, C. Federico, O. Forni, G. Magni, Origin and thermal evolution of icy satellites. *Surv. Geophys.* **16**, 533–591 (1995)