Serpentinization in the Thermal Evolution of Icy Kuiper Belt Objects in the Early Solar System

In this part of the paper we focus solely on the algorithm that is able to follow the serpentinization process itself. We assume that the reaction described by Equation (1) is the dominant process for serpentinite production in the planetesimals, and we take only this reaction into account when we calculate the heat released, following Góbi & Kereszturi (2017). The planetesimal used in the model is made of seven components: silicate rocks are forsterite (Mg-rich olivine), enstatite (Mg pyroxen), hydrated rock (serpentinite), nonreactive solid material, and H₂O in the three-phase state—liquid, solid, and void space filled with H₂O vapor. The melting point of H₂O (268 K) is obtained from the properties of a saturated solution of MgSO₄ (Kargel 1998), which was used in models in earlier studies too (Cohen & Coker 2000; Góbi & Kereszturi 2017). While MgSO₄ is quite common in some chondrites (Zolensky et al. 1999), in the case of TNOs other solutes could also be considered. However, the effect of any salt or salt combination on the melting temperature of water is probably small.

To be comparable with the Góbi & Kereszturi (2017) model, the serpentinization algorithm itself is tested in a selected layer of a sphere, with the homogeneous and isotropic distribution of the components within that layer. Due to the simple forward Euler scheme, the model remains sensitive to the choice of the time step. We explore this in Section 2.2 to find the optimal (largest allowable) time step for our calculations. Apart from the improvements discussed below, this is the same scheme as was used in Góbi & Kereszturi (2017) (see the model scheme in Figure 1).

Zoom In Zoom Out Reset image size

Heat capacity. Heat properties of materials and the average values of the physical properties of planetesimals depend on the temperature and composition (see Section A.1). We calculated the heat capacity of H₂O in a different, more accurate way than in the base model. After fixing it, the final heat production was higher by a few kelvin. However, the difference remained within the uncertainty range, even at longer timescales.

Lithospheric pressure. Several quantities depend on the total pressure, which is calculated as the sum of lithospheric and vapor pressures (see Section A.3).

Góbi & Kereszturi (2017) used a simplification to determine the lithospheric pressure owing to the small size (R = 15 km) of their main test objects, as the contribution of the lithospheric pressure to total pressure is insignificant for small objects. This was replaced by Equation (A33), which gives reliable results for larger objects too. Using this method, the reaction rate becomes lower, which makes a more significant difference for larger objects. We compare our results with the Góbi & Kereszturi (2017) lithospheric pressure values in Figure 2.

Zoom In Zoom Out Reset image size

Interfacial liquid water. To obtain the heat gain, we need to know the exact amount of liquid water and the latent heat of water vaporization (Section A.4). The latent heat was calculated in a way different from the base model, and it caused a few tenths of kelvin difference in the final temperature, with a very small amount of vapor formed under typical conditions. In those cases when the temperature is lower than the melting point of H₂O, ice is present instead of liquid water, which would not allow serpentinization to take place. However, Góbi & Kereszturi (2017) pointed out that the presence of microscopic-scale interfacial water is possible on the surface of olivine grains at low temperatures, making serpentinization possible, and it may gradually melt the icy surrounding of silicate particles owing to the heat produced in this exothermic reaction. The model of Góbi & Kereszturi (2017) allowed the formation of more interfacial water than would have actually been possible at the given temperature. They did not take into account the rate of serpentinization, although this effect is significant in the temperature range below the melting point of bulk ice. In our model we consider four values for the amount of water: (i) the actual amount of ice in the layer examined, (ii) the maximum value of interfacial water that can be formed, (iii) the amount of water that is able to react in a specific step, and (iv) the amount of ice that can be melted by the serpentinization heat (see Equation (A42)). The minimum of these four values determines the actual amount of interfacial water that can be formed.

In the Góbi & Kereszturi (2017) paper the improper handling of the interfacial water also caused an issue in selecting the proper simulation time step: if larger time steps for the same time span were used, the calculated initial temperatures required for the same serpentinization level were lower and the final temperatures higher, causing a greater overall heat production in the simulations. The reason behind this phenomenon was that in the case of decreased time steps the serpentinization reaction started later in time. By incorporating the new condition in the melting ice calculation, the simulations did not produce more interfacial water any more than was possible at the given temperature, and there was no difference in the amount of heat produced as a function of timescale at the same time.

To demonstrate the capabilities of our improved algorithm to follow the serpentizination process, we used the same setup as Góbi & Kereszturi (2017), but with the corrections listed above. We consider the deepest 100 m radius, without taking other heat sources and heat/material transfer into account. In each case the amount of nonreacting material of the test object was constant, 14% of the volume of planetesimals. We investigated how the different parameters of the model can influence the outcome compared with the Góbi & Kereszturi (2017) results.

Initial temperature. The output of any evolutionary model that includes chemical reactions like serpentinization is very sensitive to the choice of initial temperature. We examined the temperature increase (ΔT) during the reaction, as well as the time needed to consume 90% of the reagent material during the serpentinization reaction (t₉₀) at different initial temperature (T_ini) values (Figure 3). The temperature dependence of chemical reactions is known in general: the higher the initial temperature, the faster the reaction. The dependence of the temperature increase on the initial temperature is a result of the temperature dependence of the heat capacity of the constituent minerals (equations for calculating the heat capacities can be found in Section A.1).

Zoom In Zoom Out Reset image size

The serpentinization time t₉₀ is significantly longer than in the previous studies (Figure 3), and the serpentinization rate is higher toward higher initial temperatures. Mainly due to the accurate calculation of the lithospheric pressure and the consideration of the maximum value of subfreezing interfacial water, t₉₀ increased. In the case of a planetesimal with a radius of 15 km, the process can be up to 3–4 times longer than in the previous calculations, depending on the initial temperature.

This is due to the fact that only a small amount of interfacial liquid water is present in the system below the freezing point and its production depends on the amount of water reacted, and thus the process is slow. When the temperature reaches the melting point, the heat production is entirely used to melt the ice. Above the melting point, the reaction speeds up owing to the accessibility of large amounts of liquid water. Starting the reaction below the melting temperature, most of the heat is consumed by melting the ice, and this results in an overall smaller temperature increase.

Importance of microscopic liquid water. Interfacial liquid water is an important component, as it promotes the progress of the reaction in the early stage when the temperature is below the melting point of ice. Compared with previous results, it can be seen that in the subfreezing initial temperature range the process is slower and the heat production rate is also lower (Figure 3). The reaction rate increases after the melting point.

Size of planetesimals. The reaction rate increases with increasing pressure, and it depends on both the size of the test object through the lithospheric pressure and the vapor pressure (Section A.3). The pressure dependence of serpentinization was studied previously (Martin & Fyfe 1970; Wegner & Ernst 1983; Cohen & Coker 2000; Jones & Brearley 2006), and it was found to be nearly linear in the range of 1–200 MPa, which corresponds to the pressure expected in the size range in our investigation.

The initial temperature required for serpentinization to consume 90% of the materials in 10,000 yr (T₉₀) reaches values similar to those determined in the previous work (Góbi & Kereszturi 2017); the deviation is only a few kelvin (see Figure 4). The largest deviation occurs at the size/pressure value where the initial temperature falls below the freezing point.

Zoom In Zoom Out Reset image size

There is an increasing trend in temperature change (ΔT) with increasing size (see Figure 4) due to heat capacities that are lower owing to the lower initial temperature (T₉₀). When the required initial temperature falls below the freezing point, the heat production rate decreases significantly, as some of the heat produced is used to melt ice. From this point on, the ΔT again shows an increasing trend as a function of object size.

Olivine to H₂ O ratio and initial porosity. The component ratio of olivine to water is one of the most important initial parameters of the serpentinization reaction, which significantly influences both the course and the outcome of the reaction. In the case of higher olivine-to-water ratios, the reaction can be faster than in the cases of higher water content (see, e.g., Figure 1 in Góbi & Kereszturi 2017). For each olivine-to-water ratio value, a higher T₉₀ is needed (see, e.g., Figure 4), but the difference is only a few kelvin. The rate of heat production is closely related to heat capacity, which increases with water content, as the heat capacity of water is higher than that of rocky components. This slows down the reaction toward higher water content.

Time step. As we use a simple forward Euler scheme in our model, the results are expected to be sensitive to the time step applied. Choosing a large time step (Δt = 10–100 yr, or larger) results in considerable instability in the calculations, i.e., the final results (e.g., T₉₀ or ΔT we investigated above) do not show a consistent trend and depend strongly on the time step chosen. This annoying effect disappears when the time step is decreased, and the calculations become stable for Δt ≲ 1 yr (see Figure 5). To avoid this problem, we used a time step of Δt = 0.5 yr. Reducing the time step further did not cause a considerable change in the final results.

Zoom In Zoom Out Reset image size

The main source of heat within planetesimals is the radioactive decay of both short- and long-lived radionuclides. In Figure 6 we plot the total energy output as a function of time for some short-lived radionuclides. Because serpentinization is a rapid process, only a few × 10⁴ yr, we investigated the early stages of the heat development of planetesimals. During this period, as shown in Figure 6, the most significant portion of the heat production from radioactive decay was provided by ²⁶Al (Lugaro et al. 2018), and only this isotope was considered in our thermal evolution model as a heat source, with the following parameters: half-life is t_1/2 = 0.717 Myr, total energy is TE = 5.07 × 10⁶ J kg⁻¹, and we assume that ²⁶Al was homogeneously distributed in the early solar system (Lugaro et al. 2018).

Zoom In Zoom Out Reset image size

There are two major methods of heat transfer inside a planetary body: thermal conduction and convection. Both of them were considered, and the thermal radiation from the surface was also included in the thermal evolution model.

Thermal radiation. The power output was calculated from the Stefan-Boltzmann law using a constant emissivity factor ( = 0.9) and assuming an ambient temperature of T_amb = 50 K, which is a typical surface temperature in the outer solar system (Equation (A48)).

Thermal conduction. The heat conduction rate was estimated with Equation (A49), in all boundaries of all layers following Hussmann et al. (2006).

Thermal convection. If the temperatures within the layers are sufficiently high and water is present in a liquid state, convective motions can occur in the porous media (Hewitt et al. 2014). The requirement for starting convection is that the Rayleigh number (Ra) exceeds a specific value. Ra is calculated in a porous material as

$$\begin{eqnarray}&&\mathrm{Ra}=\displaystyle \frac{\rho \,{\beta }_{0}\,{\rm{\Delta }}T\,{k}_{0}\,g\,{dr}}{\eta \,\alpha },\end{eqnarray} \tag{ 3 }$$

where ρ is the density in the layer; β₀ is the thermal expansion coefficient, and we used β₀ = 10⁻³ as a safe upper limit for any of the possible constituents; ΔT is the temperature difference across distance dr, which is the thickness of the layer; and k₀ is the permeability, and we used two values, 10⁻¹² m² following Cohen & Coker (2000) and 10⁻⁷ m² following Bear (1972). The results were the same because the Ra did not reach the critical value required for the start of convection with any of the permeability values. g is the local gravitational acceleration, α = k/(ρ · c_p ) is the thermal diffusivity (k is thermal conductivity, and c_p is specific heat capacity), and η is the dynamic viscosity of the fluid:

$$\begin{eqnarray}&&\eta (T)={\eta }_{0}\,\exp (25\times ({T}_{\mathrm{melt}}/T-1)),\end{eqnarray} \tag{ 4 }$$

where η₀ = 10¹³ Pa s is the melting point viscosity (Hussmann et al. 2006) and T_melt is the melting point. A critical Rayleigh number of 1000 was used for the onset of convection following Hussmann et al. (2006). However, this critical value can be as low as 40 (Nield & Bejan 1999) in a porous medium. In our simulations these Ra values have not been reached at T ≲ T_melt.

The main purpose of our present work is to follow the evolution of the serpentinization process and to determine the conditions and onset timescales this process can operate at. As we show below, before and during the active serpentinization stage the temperature in all of our simulations is close to or below the melting point of ice, and the heat transfer is governed by these "icy" conditions. During serpentinization, the melted water is quickly consumed by the chemical reaction, preventing convection (which would be caused by the otherwise significantly decreased viscosity of liquid water). After the end of the serpentinization stage, radiogenic decay can still provide extra heat to further melt the ice. In these conditions the low viscosity of liquid water may lead to convection, generating a faster heat transfer within the inner layers, and may also lead to rearrangements in the structure of these layers. The consideration of this process, as well as, e.g., that of the diffusion of water into solid grains, is beyond the scope of this paper.

Serpentine is a primary candidate to explain the reflectance spectra of TNOs (e.g., Protopapa et al. 2009), and it has also been considered as primordial material in impact simulations aimed to explain the formation of trans-Neptunian binary systems (e.g., the Pluto-Charon system; Canup 2005). Therefore, it is an intriguing question how serpentine can form under the conditions in the trans-Neptunian region, where objects have typically high ice contents and possibly high porosity, and how far it can contribute to the heat budget and thermal evolution of these planetesimals.

We consider a spherical body with spherical symmetry (all variables depend on the radius only), with a size-dependent number of layers. We assumed a homogeneous and isotropic composition at the start and a common T_ini in the whole planetesimal. The effect of deviations from the latter assumption is investigated in Section 3.3. Heat production by radiogenic decay and heat transfer is considered as described in Sections 2.3 and 2.4.

Composition. The low density of TNOs is due partly to their composition, as the H₂O content of these distant objects is significantly higher than for those in the inner solar system, but it is also due to their higher porosity. In the case of smaller TNOs, when the radius does not exceed ∼150 km, the porosity can reach 60% (Bierson & Nimmo 2019). Such a high porosity is only possible if they formed after the decay of ²⁶Al to maintain their high porosity, as there is no internal transformation by heat at this stage. In the absence of accurate knowledge of the internal composition, a simplified composition was used, which has been used in previous serpentinization studies (see Góbi & Kereszturi 2017; Cohen & Coker 2000; see also the Appendix).

Porosity. Porosity is a very important factor in this model; the initial temperature (T_ini) and the serpentinization time are also strongly dependent on it. To estimate the porosity, we used the calculations from Yasui & Arakawa (2009), where they determine a size/pressure-dependent porosity for small icy bodies. The equation defined for the range of Regime 3 (P > 2 MPa) was applied owing to the pressure conditions in our objects:

$$\begin{eqnarray}&&\phi ={a}_{3}{P}^{{b}_{3}},\end{eqnarray} \tag{ 5 }$$

where P is the lithospheric pressure (in MPa) in a middle layer in the measured bodies and a₃ and b₃ are constants. We used the approximate values of a₃ = 0.5 and b₃ = −0.2 (see Yasui & Arakawa 2009). With these calculations, the porosity values are between 13% and 34% in the examined size range for the 42:58 rock/water ratio specified earlier (see Figure 7).

Zoom In Zoom Out Reset image size

Temperature. In the early solar system the initial temperature of the objects was determined by the size, composition, accretion heat, and the heat produced by the decay of short-lived radioactive nuclei. The heating from solar irradiation may be important closer to the Sun, but it is negligible in the outer regions. In our model, the initial temperature was taken to be the heat from the accretion (Hanks & Anderson 1969):

$$\begin{eqnarray}&&{T}_{\mathrm{acc}}=\frac{3}{5}\frac{{GM}}{R{C}_{p}},\end{eqnarray} \tag{ 6 }$$

where G is the gravitational constant, C_p is the average heat capacity, M is the mass, and R is the radius of the planetesimal (see Figure 7). At t = 0 the planetesimal has a homogeneous temperature distribution with T_acc. Note that since C_p is temperature dependent, Equation (6) is an implicit equation for T_acc, and it is solved iteratively; it also determines the shape of the T_acc curve presented in Figure 7.

In a next step we tested the sensitivity of the radioactive decay heat product on the time step and layer thickness. In these tests we used a rock/water ratio of 42:58, which corresponds to an olivine-to-water ratio of 0.12. Radiogenic heat production was not very sensitive to the variation of time step and layer thickness.

In all simulation configurations, we used a nonreactive rock content of 14% with a density of 3630 kg m⁻³ as in earlier studies (Góbi & Kereszturi 2017). In each case, the olivine-to-water ratio was taken as 0.12, which corresponds to the assumed rock/water ratio of 42:58 in TNOs. The calculations are performed with the size-dependent porosity (Equation (5)), the initial temperature is calculated from the accretion heat (Equation (6)), and we choose a time step of 0.5 yr and a layer thickness of 20 km. The starting time t = 0 was considered to be the end of the accretion.

Objects with a radius below 150 km are expected to have high porosity (up to 60%), and the initial (accretion) temperature would remain very low (<42 K), as indicated by our calculations. These objects may have formed after the depletion of ²⁶Al (as discussed in Bierson & Nimmo 2019), and the serpentinization reaction could likely not produce a significant amount of heat. Therefore, we did not consider these objects in our further calculations. Very large objects (R > 1000 km) were also excluded from the study size range because they are likely formed by multiple accretion events and undergo fast chemical differentiation early in their evolution.

We examined the conditions under which serpentinization can produce significant heat over a few tens of thousands of years in the size range of 150 km ≤ R ≤ 1000 km. We found that for objects smaller than 600 km the initial temperature from the accretion will be so low that the reaction cannot start or will be very slow, even in the core, in the absence of radiogenic decay. This appears in Figure 8 as a "jump" in temperature at this specific size. In larger objects, the process is already able to produce a notable amount of heat without the contribution of radioactive elements. Above an initial temperature of ∼150 K, assuming the previously determined material composition and pressure-dependent porosity (Figure 7), the reaction rate increases significantly (Figure 8). As the size increases, both the initial temperature and the internal pressure increase, resulting in an increase in the rate of serpentinization. This results in the extension of the central region, where the serpentinization process can produce significant heat in a few tens of thousands of years. The top panel of Figure 9 shows the ratio of the internal volume in which 90% olivine has been consumed in the first 100,000 yr for a planetesimal with a specific size. For radii R ≥ 650 km the serpentinized zone extends almost to the surface, while it remains in the core below R = 600 km, in agreement with the slow reaction rate indicated by the small temperature increase seen in Figure 8. The bottom panel of Figure 9 shows an example of how serpentinization proceeds inside an object of R = 620 km in the first 100,000 yr. As shown also in the top panel, a maximum of 480 km radius of the serpentinized region is reached, corresponding to a maximum volume ratio of 46%. This R ≈ 600 km critical size limit for the efficient progress of the serpentinization process was obtained without the consideration of radiogenic decay, and it already indicates that serpentinization can be an efficient process for the large objects (the dwarf planets of the trans-Neptunian region) even if they formed after 10⁷ yr or later, when the heat produced by radiogenic decay is significantly lower (see Figure 6). This also suggests that the size limit of efficient serpentinization is smaller when an additional heat source—radiogenic decay—is considered.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To investigate this, we modeled the thermal evolution of test objects with radii between 150 km < R < 500 km considering all heat sources (accretion heat, radiogenic decay) in addition to serpentinization itself. In Figure 10 we present the thermal evolution of three objects with radii 200, 340, and 500 km. The serpentinization process starts to produce a significant amount of heat when the temperature of the layer exceeds a critical value of ∼180 K, at time t₀ after the start of the simulation. This t₀ depends strongly on the size of the object (see Figure 11, middle panel) and decreases with size, in a similar way to the end time of the reaction (t₉₀), represented by vertical lines in Figure 10. In these simulations serpentinization completes fully in all the studied layers, even at the melting point, and therefore it does not contribute to the future thermal evolution. We note that all our simulations have been run beyond the time of completion of the reaction. In our simulations serpentinization reduces the time needed to reach the melting temperature of water by a factor of 1.6–2.6, depending on the size of the object (Figure 11, bottom panel). The contribution of the different processes to the temperature of the object at the end of serpentinization (t₉₀) is presented in Figure 11 (top panel). At smaller sizes radiogenic decay contributes the most, and the importance of accretion heat increases notably with size. The heat obtained from serpentinization (50–80 K) becomes more important with growing size and exceeds the contribution from radiogenic decay for R ≈ 500 km. In large enough objects (where the critical reaction-starting temperature of ∼180 K is reached) serpentinization proceeds quickly, and the whole process is finished in ∼10⁴ yr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results of Wakita & Sekiya (2011) show that the formation time is very important because 2.4 Myr after Ca-Al-rich inclusions formation the icy planetesimals may not reach the melting temperature of ice. We also examined the effect of formation time on the thermal evolution, considering the same setup as above.

In the bottom right panel of Figure 10 we demonstrate the effect of a delayed start (by 0.7 Myr, the half-life of ²⁶Al) and hence reduced radiogenic heat. In this simulation the t₀, t₁, and t₂ timescales are notably (by a factor of ∼2) longer.

Despite these longer timescales, serpentinization can still fully proceed, and the temperature reaches the melting point of ice. As is expected, the contribution from radiogenic decay to the final temperature is smaller (11 K in this specific case), and the relative contribution of the temperature increase due to serpentinization is higher.

In Figure 12 we show the effect of a late formation of the planetesimals, and consequently a late start of the serpentinization process, for different starting times. The late formation results in a reduced amount of heat from radiogenic decay, due to the depletion of ²⁶Al. These simulations were performed for R = 200 km objects, assuming that the formation of the planetesimal happened at a t = [0,1,...7] × t_1/2 after the onset of isotopic decay, when the abundance of ²⁶Al was at its maximum, as considered in our previous simulations. Serpentinization needs a much longer time to start, either in the core or in layers closer to the surface, for later formation times. While this t₀ is on the order of a few thousand years for an early formation, it is several million years for a formation at t ≈ 5 Myr. The t₉₀ timescale, i.e., the time when 90% of the potentially serpentine-forming material in the core of the planetesimal is consumed, is also notably longer, while the ratio of t₀ to t₉₀ remains roughly the same. Altogether the lower temperatures due to the smaller radiogenic heat slow down serpentinization considerably.

Zoom In Zoom Out Reset image size

We examined what changes occur when the initial temperature distribution is inhomogeneous and the surface is warmer (Figure 13). In this test we used an object with a radius of 240 km, and the initial temperature distribution was set in a way that the surface temperature was 5/3 times the central temperature (Hanks & Anderson 1969). Despite the lower starting temperature, the serpentinization process is the fastest in the core owing to the higher lithospheric pressure. The reaction takes longer in the outer layers, and as the initial temperature is higher, the final temperature will also be higher, but the temperature difference is smaller than in the homogeneous case.

Zoom In Zoom Out Reset image size

We also examined the case when the initial temperature distribution was inhomogeneous and the formation occurred in two steps (Figure 14). First, an object formed with a radius of 160 km, and τ_f = 10,000 or 100,000 yr later the formation was completed and the object reached a final radius of 320 km. In the time between the two formation/accretion events the object may have warmed up from both radiogenic decay and serpentinization. We assumed two values for the time of the first formation event: t_s = 0, the maximum ²⁶Al heat production date, and t_s = 3 t_1/2, i.e., three ²⁶Al half-lives later. For τ_s = 100,000 yr and t_s = 0 the serpentinization could go along all the way before the second accretion event in the center and left the core at a high temperature. Due to the lack of further reactants, the core was heated by radiogenic decay only from this point on (almost straight blue curve in Figure 14, bottom left). The final object is built on this warm core in the second accretion event. Also in our other three cases, the results show that serpentinization is faster in the core independently of the initial conditions. The difference between the core and the outer layers is more significant for t_s = 0, and, as expected, everything occurs significantly later for t_s = 3 t_1/2 yr.

Zoom In Zoom Out Reset image size

Heat capacity (c_p ) and density (ρ) calculations were obtained from Cohen & Coker (2000). These physical properties of all components are temperature (T) dependent:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\,{C}_{p}(\mathrm{olivine}) & = & -11.32+13.58\,x\\ & & -\,4.25\,{x}^{2}+0.44\,{x}^{3}\end{array}\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\,{C}_{p}(\mathrm{enstatite}) & = & -8.620+10.39\,x\\ & & -3.00\,{x}^{2}+0.28\,{x}^{3}\end{array}\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\,{C}_{p}(\mathrm{water}) & = & 8.250-4.180\,x\\ & & +1.12\,{x}^{2}-0.076\,{x}^{3},\end{array}\end{eqnarray} \tag{ A3 }$$

where $x=\mathrm{log}\,T$. These expressions are valid from approximately 50 to 500 K:

$$\begin{eqnarray}\begin{array}{rcl}{C}_{p}(\mathrm{serpentinite}) & = & 1145+0.048\,T\\ & & -2.65\times {10}^{7}\,{T}^{-2}\quad (T\,\gt \,273\,{\rm{K}})\\ \mathrm{log}\,{C}_{p}(\mathrm{serpentinite}) & = & -0.59-1.51\,x\\ & & +\,(2.82\,{x}^{2})-(0.66\,{x}^{3})\\ & & (T\,\lt \,273\,{\rm{K}})\end{array}\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}&&{C}_{p}(\mathrm{vapor})=1730.54+0.45\,T\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}\begin{array}{rcl}{C}_{p}(\mathrm{ice}) & = & 152.46+7.12\,T\quad (T\,\gt \,150\,{\rm{K}})\\ & = & 126.89+7.50\,T\quad (150\,{\rm{K}}\,\gt \,T\,\gt \,95\,{\rm{K}})\\ & = & -49.97+9.5\,T\quad (95\,{\rm{K}}\,\gt \,T\,\gt \,50\,{\rm{K}}).\end{array}\end{eqnarray} \tag{ A6 }$$

For a specific layer of the object, we calculated the average heat capacity from the components' heat capacities and masses (m):

$$\begin{eqnarray}\begin{array}{rcl}{C}_{p}(\mathrm{layer}){m}_{\mathrm{lay}} & = & {C}_{p}(\mathrm{oli}){m}_{\mathrm{oli}}\\ & & +{C}_{p}(\mathrm{ens}){m}_{\mathrm{ens}}+{C}_{p}(\mathrm{ice}){m}_{\mathrm{ice}}\\ & & +{C}_{p}(\mathrm{wat}){m}_{\mathrm{wat}}+{C}_{p}(\mathrm{vap}){m}_{\mathrm{vap}}\\ & & +{C}_{p}(\mathrm{ser}){m}_{\mathrm{ser}}+{C}_{p}(\mathrm{oli}){m}_{\mathrm{nre}}.\end{array}\end{eqnarray} \tag{ A7 }$$

In the cases of olivine, enstatite, serpentinite, and nonreactive material, constant values of density were used (see in Table 2). For H₂O we used

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{\mathrm{wat}} & = & -221+13.1\,T-0.050\,7\,{T}^{2}\\ & & +\,8.49\times {10}^{-5}\,{T}^{3}-5.48\times {10}^{-8}\,{T}^{4}\end{array}\end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{\mathrm{ice}} & = & -46.9\,x+1032.71\quad (T\,\gt \,137\,{\rm{K}})\\ & = & -1.32\,x+935.32\quad (T\,\lt \,137\,{\rm{K}})\end{array}\end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray}&&{\rho }_{\mathrm{vap}}={P}_{\mathrm{vap}}\,0.018\,{T}^{-1}\,{R}_{g}^{-1}\end{eqnarray} \tag{ A10 }$$

where P_vap is the vapor pressure of the vapor and R_g is the gas constant. For the examined layer, we calculated the average density from the components' masses and volumes (V):

$$\begin{eqnarray}&&{\rho }_{\mathrm{lay}}=\displaystyle \frac{{m}_{\mathrm{oli}}+{m}_{\mathrm{ens}}+{m}_{\mathrm{ice}}+{m}_{\mathrm{wat}}+{m}_{\mathrm{vap}}+{m}_{\mathrm{ser}}+{m}_{\mathrm{nre}}}{{V}_{\mathrm{oli}}+{V}_{\mathrm{ens}}+{V}_{\mathrm{ice}}+{V}_{\mathrm{wat}}+{V}_{\mathrm{ser}}+{V}_{\mathrm{voi}}+{V}_{\mathrm{nre}}},\end{eqnarray} \tag{ A11 }$$

where V_voi is the void space that is filled with H₂O vapor.

This part of the algorithm considers a specific layer in a body, at a specific depth. The volume of the examined layer is as follows:

$$\begin{eqnarray}&&{V}_{\mathrm{lay}}=\left(\displaystyle \frac{4}{3}\,\pi \,{r}_{\mathrm{lt}}^{3}\right)-\left(\displaystyle \frac{4}{3}\,\pi \,{r}_{\mathrm{lb}}^{3}\right),\end{eqnarray} \tag{ A12 }$$

where V_lay is the volume of the examined layer, r_lt is the top of the layer, and r_lb is the bottom of the layer.

The initial values of volumes and masses of components (marked with (0)) were calculated from the porosity (ϕ), olivine-to-water ratio (n_oli/n_wat), and the fraction of nonreactive material (n_nre):

$$\begin{eqnarray}&&{V}_{\mathrm{voi}}(0)={V}_{\mathrm{lay}}\,\phi \end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray}&&{V}_{\mathrm{nre}}={V}_{\mathrm{lay}}\,{n}_{\mathrm{nre}}\end{eqnarray} \tag{ A14 }$$

$$\begin{eqnarray}&&{m}_{\mathrm{nre}}={V}_{\mathrm{nre}}\,{\rho }_{\mathrm{nre}}\end{eqnarray} \tag{ A15 }$$

$$\begin{eqnarray}&&{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}(0)={n}_{\mathrm{oli}}/{n}_{\mathrm{wat}}\,140.71/18\end{eqnarray} \tag{ A16 }$$

$$\begin{eqnarray}&&{m}_{\mathrm{lay}}={V}_{\mathrm{lay}}\,{\rho }_{\mathrm{lay}}\end{eqnarray} \tag{ A17 }$$

$$\begin{eqnarray}\begin{array}{rcl}{V}_{\mathrm{ens}}(0) & = & ({V}_{\mathrm{lay}}-({V}_{\mathrm{nre}}+{V}_{\mathrm{voi}}))\\ & & \times \,\displaystyle \frac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}/1.3/{\rho }_{\mathrm{ens}}}{\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}}{{\rho }_{\mathrm{oli}}}+\tfrac{1}{{\rho }_{\mathrm{ice}}}+\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}/1.3}{{\rho }_{\mathrm{ens}}}}\quad (T\,\lt \,{T}_{\mathrm{melt}})\\ & = & ({V}_{\mathrm{lay}}-({V}_{\mathrm{nre}}+{V}_{\mathrm{voi}})\\ & & \times \,\displaystyle \frac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}/1.3/{\rho }_{\mathrm{ens}}}{\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}}{{\rho }_{\mathrm{oli}}}+\tfrac{1}{{\rho }_{\mathrm{wat}}}+\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}/1.3}{{\rho }_{\mathrm{ens}}}}\quad (T\,\gt \,{T}_{\mathrm{melt}})\end{array}\end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray}\begin{array}{rcl}{V}_{\mathrm{ice}}(0) & = & ({V}_{\mathrm{lay}}-({V}_{\mathrm{nre}}+{V}_{\mathrm{voi}}))\\ & & \times \,\displaystyle \frac{1/{\rho }_{\mathrm{ice}}}{\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}}{{\rho }_{\mathrm{oli}}}+\tfrac{1}{{\rho }_{\mathrm{ice}}}+\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}/1.3}{{\rho }_{\mathrm{ens}}}}\quad (T\,\lt \,{T}_{\mathrm{melt}})\\ & = & 0\quad (T\,\gt \,{T}_{\mathrm{melt}})\end{array}\end{eqnarray} \tag{ A19 }$$

$$\begin{eqnarray}\begin{array}{rcl}{V}_{\mathrm{wat}}(0) & = & 0\quad (T\,\lt \,{T}_{\mathrm{melt}})\\ & = & ({V}_{\mathrm{lay}}-({V}_{\mathrm{nre}}+{V}_{\mathrm{voi}}))\\ & & \times \,\displaystyle \frac{1/{\rho }_{\mathrm{wat}}}{\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}}{{\rho }_{\mathrm{oli}}}+\tfrac{1}{{\rho }_{\mathrm{wat}}}+\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}/1.3}{{\rho }_{\mathrm{ens}}}}\quad (T\,\gt \,{T}_{\mathrm{melt}})\end{array}\end{eqnarray} \tag{ A20 }$$

$$\begin{eqnarray}\begin{array}{rcl}{V}_{\mathrm{oli}}(0) & = & ({V}_{\mathrm{lay}}-({V}_{\mathrm{nre}}+{V}_{\mathrm{voi}}))\\ & & \times \,\displaystyle \frac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}/{\rho }_{\mathrm{oli}}}{\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}}{{\rho }_{\mathrm{oli}}}+\tfrac{1}{{\rho }_{\mathrm{ice}}}+\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}/1.3}{{\rho }_{\mathrm{ens}}}}\quad (T\,\lt \,{T}_{\mathrm{melt}})\\ & = & ({V}_{\mathrm{lay}}-({V}_{\mathrm{nre}}+{V}_{\mathrm{voi}}))\\ & & \times \,\displaystyle \frac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}/{\rho }_{\mathrm{oli}}}{\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}}{{\rho }_{\mathrm{oli}}}+\tfrac{1}{{\rho }_{\mathrm{wat}}}+\tfrac{{m}_{\mathrm{oli}}/{m}_{\mathrm{wat}}/1.3}{{\rho }_{\mathrm{ens}}}}\quad (T\,\gt \,{T}_{\mathrm{melt}}).\end{array}\end{eqnarray} \tag{ A21 }$$

The melting point of ice (T_melt = 268 K) is obtained considering a saturated solution of MgSO₄.

The masses of components (enstatite, ice, water, olivine) are calculated in the zeroth time step in the following way:

$$\begin{eqnarray}&&{m}_{j}(0)={V}_{j}(0){\rho }_{j},\end{eqnarray} \tag{ A22 }$$

where subscript j refers to the materials. In all other time steps these masses are obtained as

$$\begin{eqnarray}&&{m}_{\mathrm{ser}}={m}_{\mathrm{ser}}^{i-1}+\displaystyle \frac{277.1\,{\rm{\Delta }}{n}_{\mathrm{ser}}^{i-1}}{1000}\end{eqnarray} \tag{ A23 }$$

$$\begin{eqnarray}&&{m}_{\mathrm{oli}}={m}_{\mathrm{oli}}^{i-1}-({\rm{\Delta }}{n}_{\mathrm{ser}}^{i-1}\times 0.140\,71)\end{eqnarray} \tag{ A24 }$$

$$\begin{eqnarray}&&{m}_{\mathrm{ens}}={m}_{\mathrm{ens}}^{i-1}-({\rm{\Delta }}{n}_{\mathrm{ser}}^{i-1}0\times 0.100\,39)\end{eqnarray} \tag{ A25 }$$

$$\begin{eqnarray}&&{m}_{\mathrm{ice}}=\max (0,{m}_{\mathrm{ice}}^{i-1}-{\rm{\Delta }}{m}_{\mathrm{ice}}^{i-1})\end{eqnarray} \tag{ A26 }$$

$$\begin{eqnarray}&&{m}_{\mathrm{wat}}={m}_{\mathrm{wat}}^{i-1}+{\rm{\Delta }}{m}_{\mathrm{ice}}^{i-1}-2\,{\rm{\Delta }}{n}_{\mathrm{ser}}^{i-1}\times 0.018)\end{eqnarray} \tag{ A27 }$$

$$\begin{eqnarray}&&{m}_{\mathrm{vap}}=\displaystyle \frac{0.018\,{P}_{\mathrm{vap}}\,{V}_{\mathrm{voi}}}{{R}_{g}\,T},\end{eqnarray} \tag{ A28 }$$

where the i − 1 superscript refers to the value of the variable in the previous time step,

Δn_ser is the serpentine produced and Δm_ice is the changing mass of ice due to the melting. We calculated the volumes of components from the masses and densities:

$$\begin{eqnarray}&&V=m/\rho \end{eqnarray} \tag{ A29 }$$

$$\begin{eqnarray}&&{V}_{\mathrm{voi}}={V}_{\mathrm{voi}}^{i-1}+\displaystyle \frac{3}{{10}^{6}}\times {\rm{\Delta }}{n}_{\mathrm{ser}}^{i-1}.\end{eqnarray} \tag{ A30 }$$

Amount of substance of water (n_wat ) and olivine (n_oli ):

$$\begin{eqnarray}&&{n}_{\mathrm{wat}}=1000\times {m}_{\mathrm{wat}}/18\end{eqnarray} \tag{ A31 }$$

$$\begin{eqnarray}&&{n}_{\mathrm{oli}}=1000\times {m}_{\mathrm{oli}}/140.71.\end{eqnarray} \tag{ A32 }$$

To obtain the total pressure, the sum of the lithospheric (P_lit) and vapor pressures (P_vap), we first calculate the gravitational acceleration (g(l)) and then P_lit in that specific layer:

$$\begin{eqnarray}\begin{array}{rcl}g(l) & = & G\displaystyle \frac{{m}_{r}}{{r}^{2}}\\ {P}_{\mathrm{lit}} & = & \displaystyle \sum _{l=r}^{R}{\rho }_{\mathrm{lay}}(l)g(l){dr},\end{array}\end{eqnarray} \tag{ A33 }$$

where G is the gravitational constant, r is the radius of the specific layer, m_r is mass within the radius r, R is the radius of the planetesimal, and dr is the thickness of the layer.

From the Clausius−Clapeyron relation, the vapor pressure has the approximate form

$$\begin{eqnarray}&&{P}_{\mathrm{vap}}={P}_{0}\,{e}^{\tfrac{{T}_{0}}{T}},\end{eqnarray} \tag{ A34 }$$

where P₀ = 3.58 × 10¹² Pa and T₀ = − 6140 K in the presence of ice, while it is changed to P₀ = 4.7 × 10¹⁰ Pa and T₀ = − 4960 K in the presence of water (Grimm & Mcsween 1989). In the case of mixed phases of H₂O, the mass-weighted average of P_vap is used (Góbi & Kereszturi 2017):

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{vap}}(0) & = & \displaystyle \frac{{m}_{\mathrm{ice}}\,{P}_{0}\,{e}^{\tfrac{{T}_{0}}{T}}+{m}_{\mathrm{wat}}\,{P}_{0}\,{e}^{\tfrac{{T}_{0}}{T}}}{{m}_{\mathrm{ice}}+{m}_{\mathrm{wat}}}\\ {P}_{\mathrm{vap}} & = & \displaystyle \frac{{m}_{\mathrm{ice}}^{i-1}\,{P}_{0}\,{e}^{\tfrac{{T}_{0}}{T}}+{m}_{\mathrm{wat}}^{i-1}\,{P}_{0}\,{e}^{\tfrac{{T}_{0}}{T}}}{{m}_{\mathrm{ice}}^{i-1}+{m}_{\mathrm{wat}}^{i-1}}.\end{array}\end{eqnarray} \tag{ A35 }$$

The reaction rate depends linearly on pressure and exponentially on temperature. By knowing the serpentinization rate (k_r ), the amount of serpentine produced Δn_ser can then be calculated in moles:

$$\begin{eqnarray}&&{k}_{r}=4383\,\displaystyle \frac{{P}_{\mathrm{lit}}+{P}_{\mathrm{vap}}}{{10}^{8}}\,{e}^{\tfrac{-3463}{T}}\end{eqnarray} \tag{ A36 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{n}_{\mathrm{ser}}=\min \left({n}_{\mathrm{oli}}(1-{e}^{-{kr}\times {\rm{\Delta }}t}),\displaystyle \frac{{n}_{\mathrm{wat}}}{2}\right).\end{eqnarray} \tag{ A37 }$$

When olivine is in excess, then n_wat/2 is used as the limiting reagent.

To determine the extent of heat production, it is necessary to know the heat of reaction:

$$\begin{eqnarray}&&{\rm{\Delta }}{h}_{\mathrm{ser}}={\rm{\Delta }}{H}_{r}\,{\rm{\Delta }}{n}_{\mathrm{ser}},\end{eqnarray} \tag{ A38 }$$

where ΔH_r = 69 kJ mol⁻¹ is the reaction enthalpy (Robie & Waldbaum 1968).

The vaporation heat is calculated as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{H}_{v} & = & 3713997.2-7822.6569\,T\\ & & +\,17.613373\,{T}^{2}-0.019018061\,{T}^{3}\end{array}\end{eqnarray} \tag{ A39 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{h}_{\mathrm{vap}}={\rm{\Delta }}{H}_{v}\,{\rm{\Delta }}{m}_{\mathrm{vap}},\end{eqnarray} \tag{ A40 }$$

where ΔH_v is the latent heat of vaporization of water, T is capped at 400 K, and ${\rm{\Delta }}{m}_{\mathrm{vap}}={m}_{\mathrm{vap}}-{m}_{\mathrm{vap}}^{i-1}$.

In those cases when the initial temperature is lower than the melting point of the ice, a certain part of ice can transform into interfacial water (Anderson et al. 1973):

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}\,w & = & 0.2618+0.5519\,\mathrm{ln}\,S\\ & & -\,1.449\,{S}^{-0.264}\,\mathrm{ln}({T}_{\mathrm{melt}}-T),\end{array}\end{eqnarray} \tag{ A41 }$$

where w is the content of microscopic liquid water (in g/100 g soil) and S is specific surface area (100 m² g⁻¹ based on Anderson et al. 1973).

In the next step, we calculated the amount of ice that is converted to water (Δm_ice) considering the actual amount of ice in the layer, the possible amount of interfacial water that can be formed, the amount of water that is able to react in this specific step, and the amount of ice that can be melted by serpentinization:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{m}_{\mathrm{ice}}(0) & = & \min \left({m}_{\mathrm{ice}},w\times \displaystyle \frac{{m}_{\mathrm{oli}}+{m}_{\mathrm{ens}}+{m}_{\mathrm{nre}}}{100\times {\rho }_{\mathrm{wat}}}\right)\\ & & \,(T\,\lt \,{T}_{\mathrm{melt}})\\ & = & \displaystyle \frac{{\rm{\Delta }}{h}_{\mathrm{ser}}-{\rm{\Delta }}{h}_{\mathrm{vap}}}{L}\ \ \ (T\,\gt \,{T}_{\mathrm{melt}}\quad \mathrm{and}\quad {\rm{\Delta }}{h}_{\mathrm{ser}}\,\gt \,0)\\ & = & 0\ (T\,\gt \,{T}_{\mathrm{melt}}\quad \mathrm{and}\quad {\rm{\Delta }}{h}_{\mathrm{ser}}\,\lt \,0)\\ {\rm{\Delta }}{m}_{\mathrm{ice}} & = & \min \left({m}_{\mathrm{ice}},w\times \displaystyle \frac{{m}_{\mathrm{oli}}+{m}_{\mathrm{ens}}+{m}_{\mathrm{nre}}+{m}_{\mathrm{ser}}}{100\,{\rho }_{\mathrm{wat}}},\right.\\ & & \left.2\,{\rm{\Delta }}{n}_{\mathrm{ser}}\times 0.018\right),\displaystyle \frac{{\rm{\Delta }}{h}_{\mathrm{ser}}-{\rm{\Delta }}{h}_{\mathrm{vap}}}{L})\quad (T\,\lt \,{T}_{\mathrm{melt}})\\ & = & \displaystyle \frac{{\rm{\Delta }}{h}_{\mathrm{ser}}-{\rm{\Delta }}{h}_{\mathrm{vap}}}{L}\quad (T\,\gt \,{T}_{\mathrm{melt}}\quad \mathrm{and}\quad {\rm{\Delta }}{h}_{\mathrm{ser}}\,\gt \,0)\\ & = & 0\ (T\,\gt \,{T}_{\mathrm{melt}}\quad \mathrm{and}\quad {\rm{\Delta }}{h}_{\mathrm{ser}}\,\lt \,0),\end{array}\end{eqnarray} \tag{ A42 }$$

where L is the latent heat of the ice-water phase transition

$$\begin{eqnarray}&&{\rm{\Delta }}{h}_{\mathrm{ice}/\mathrm{wat}}=L{\rm{\Delta }}{m}_{\mathrm{ice}}.\end{eqnarray} \tag{ A43 }$$

The temperature increase (ΔT_ser) due to serpentinization is calculated as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{T}_{\mathrm{ser}}(0) & = & \displaystyle \frac{{\rm{\Delta }}{h}_{\mathrm{ser}}}{{C}_{p}(\mathrm{layer})\times {m}_{\mathrm{lay}}}\\ {\rm{\Delta }}{T}_{\mathrm{ser}} & = & \displaystyle \frac{{\rm{\Delta }}{h}_{\mathrm{ser}}-{\rm{\Delta }}{h}_{\mathrm{ice}/\mathrm{wat}}-{\rm{\Delta }}{h}_{\mathrm{vap}}}{{C}_{p}(\mathrm{layer})\times {m}_{\mathrm{lay}}}.\end{array}\end{eqnarray} \tag{ A44 }$$

We considered solely the ²⁶Al isotope when calculating the heat from the radiogenic decay. The radiogenic heat production rate (Q_rad) is obtained by the following equations (Desch et al. 2009):

$$\begin{eqnarray}\begin{array}{rcl}\lambda & = & \displaystyle \frac{\mathrm{ln}(2)}{{t}_{1/2}}\\ {q}_{\mathrm{rad}} & = & \mathrm{TE}\times \lambda \\ {Q}_{\mathrm{rad}} & = & {q}_{\mathrm{rad}}\,{e}^{(-t\times \lambda )}\,{m}_{\mathrm{rock}}\,{dt},\end{array}\end{eqnarray} \tag{ A45 }$$

where λ is the decay constant, TE is the total energy, and q_rad is the radiogenic heat production rate (see Table 2).

Thermal conduction.—The thermal conductivity of the rocky material components was considered to be independent of the temperature (Robertson 1988; Cohen & Coker 2000). We calculated the thermal conductivity of the different phases of H₂O in the following way:

$$\begin{eqnarray}\begin{array}{rcl}{k}_{\mathrm{ice}} & = & 9.828\,\exp (-0.0057\,T)\\ {k}_{\mathrm{vap}} & = & -0.0143+1.02\times {10}^{-4}T\\ {k}_{\mathrm{wat}} & = & -0.581+6.34\times {10}^{-3}T\\ & & -\,7.93\times {10}^{-6}{T}^{2}\quad (T\,\lt \,410\,{\rm{K}})\\ {k}_{\mathrm{wat}} & = & 0.9721(-0.142+4.12\times {10}^{-3}T\\ & & -\,5.01\times {10}^{-6}{T}^{2})\quad (T\,\gt \,410\,{\rm{K}}).\end{array}\end{eqnarray} \tag{ A46 }$$

When determining the average thermal conductivity of the investigated layer, a composite rock with a homogeneous material distribution was considered. We calculated the parallel bulk rock conductivity (k_p , its grains arranged in a parallel orientation to the direction of heat flow) and the series conductivity (k_s , its grains arranged in a layered sequence perpendicular to the heat flow direction). The mean values of k_p and k_s fit well with the observed values, especially for more porous rocks (Robertson 1988):

$$\begin{eqnarray}\begin{array}{rcl}{k}_{p} & = & {{nV}}_{1}{k}_{1}+{{nV}}_{2}{k}_{2}+{{nV}}_{3}{k}_{3}+...\\ \displaystyle \frac{1}{{k}_{s}} & = & \displaystyle \frac{{{nV}}_{1}}{{k}_{1}}+\displaystyle \frac{{{nV}}_{2}}{{k}_{2}}+\displaystyle \frac{{{nV}}_{3}}{{k}_{3}}+...\\ {k}_{\mathrm{lay}} & = & \displaystyle \frac{{k}_{p}+{k}_{s}}{2},\end{array}\end{eqnarray} \tag{ A47 }$$

where nV₁, nV₂, nV₃, ... are fractional volumes of components (nV_component = V_component/V_lay).

Heat radiation.—On the surface of the body the radiated heat is calculated as

$$\begin{eqnarray}&&{Q}_{\mathrm{th}}=\epsilon \sigma \,({T}_{\mathrm{surf}}^{4}-{T}_{\mathrm{amb}}^{4})\,A\,{dt},\end{eqnarray} \tag{ A48 }$$

where is the emissivity factor ( = 0.9); T_surf is the temperature of the surface; T_amb is the ambient temperature, which is assumed to be 50 K, a typical surface temperature of airless bodies due to solar irradiation in the outer solar system; and A is the radiating surface area.

Thermal conduction.—The heat transferred by thermal conduction is calculated as

$$\begin{eqnarray}&&{Q}_{\mathrm{cond}}=A\,{k}_{\mathrm{lay}}\,\displaystyle \frac{{\rm{\Delta }}T}{{dr}}\,{dt},\end{eqnarray} \tag{ A49 }$$

where ΔT is the temperature difference between the adjacent layers.

Thermal convection.—If the Rayleigh number exceeded a critical value, we also calculated the Nusselt number and considered the convective heat flow in our calculations:

$$\begin{eqnarray}\begin{array}{rcl}{Nu} & = & {\left(\displaystyle \frac{\mathrm{Ra}}{{\mathrm{Ra}}_{\mathrm{crit}}}\right)}^{\beta }\\ {Q}_{\mathrm{conv}} & = & \displaystyle \frac{{k}_{\mathrm{lay}}}{{dr}}{Nu}\,A(T(l+1)-T(l)){dt},\end{array}\end{eqnarray} \tag{ A50 }$$

where β is a dimensionless value that can take values between 0.25 and 1/3 depending on the geometry and boundary conditions. As a reference value, we use β = 0.3 after Hussmann et al. (2006). We used Ra_crit = 1000 as the critical Rayleigh number.

The final heat equation is

$$\begin{eqnarray}\begin{array}{rcl}{Q}_{\mathrm{sum}} & = & {\rm{\Delta }}{h}_{\mathrm{ser}}-{\rm{\Delta }}{h}_{\mathrm{ice}/\mathrm{wat}}-{\rm{\Delta }}{h}_{\mathrm{vap}}\\ & & +\,{Q}_{\mathrm{rad}}-{Q}_{\mathrm{cool}}+{Q}_{\mathrm{cond}}+{Q}_{\mathrm{conv}},\end{array}\end{eqnarray} \tag{ A51 }$$

where Q_cool = Q_th on the surface and in the inner layers Q_cool = Q_cond at the upper boundary of the examined layer.