The Conductive Cooling of Planetesimals With Temperature‐Dependent Properties

Modeling the planetary heat transport of small bodies in the early Solar System allows us to understand the geological context of meteorite samples. Conductive cooling in planetesimals is controlled by thermal conductivity, heat capacity, and density, which are functions of temperature (T). We investigate if the incorporation of the T‐dependence of thermal properties and the introduction of a nonlinear term to the heat equation could result in different interpretations of the origin of different classes of meteorites. We have developed a finite difference code to perform numerical models of a conductively cooling planetesimal with T‐dependent properties and find that including T‐dependence produces considerable differences in thermal history, and in turn the estimated timing and depth of meteorite genesis. We interrogate the effects of varying the input parameters to this model and explore the nonlinear T‐dependence of conductivity with simple linear functions. Then we apply non‐monotonic functions for conductivity, heat capacity, and density fitted to published experimental data. For a representative calculation of a 250 km radius pallasite parent body, T‐dependent properties delay the onset of core crystallization and dynamo activity by ∼40 Myr, approximately equivalent to increasing the planetary radius by 10%, and extend core crystallization by ∼3 Myr. This affects the range of planetesimal radii and core sizes for the pallasite parent body that are compatible with paleomagnetic evidence. This approach can also be used to model the T‐evolution of other differentiated minor planets and primitive meteorite parent bodies and constrain the formation of associated meteorite samples.

chondritic meteorites contain primitive material including solids that condensed from hot gas in the Solar Nebula (MacPherson, 2014). Understanding the geological context of differentiated meteorites and their parent bodies' thermal evolution allows constraints to be placed on the formation, differentiation, and eventual breakup of planetesimals, and on the early evolution of the Solar System. In this context, models of conductive cooling of differentiated primary parent bodies are frequently used to aid the interpretation of meteorite samples. In this study we investigate the importance of including temperature dependent thermal properties in such models. We use a pallasite parent body as an example to illustrate the influence that including T-dependent properties can have on understanding the origin of meteorite samples.
One approach to understanding the formation of meteorites is to analyze the thermal processing experienced by meteorite samples and to compare this to estimated temperature conditions within the parent body using thermal evolution models. Heat flow in conductively cooling planetesimals is controlled by the material properties of their constituent minerals -thermal conductivity (k, W m −1 K −1 ), density (ρ, kg m −3 ), and heat capacity (C, J kg −1 K −1 ), in addition to the boundary conditions imposed and the geometry of the planetesimal. Large temperature gradients are expected in planetesimals, with typical surface temperatures of ∼250 K rising to ∼1800 K at the center (Bryson et al., 2015;Scheinberg et al., 2016). Planetesimals experience much lower internal pressures than planets: the center of a 250 km body with an olivine mantle and an iron core would be at ∼ 300 MPa, in comparison to Earth's central pressure of 364 GPa (Dziewonski & Anderson, 1981;Scheinberg et al., 2016). If k, ρ, and C are assumed constant, they can be expressed in terms of diffusivity    k c . This is a common approximation made in conductive cooling models of differentiated planetesimals with olivine mantles, despite temperature and pressure dependence (Bryson et al., 2015;Fu et al., 2014;Haack et al., 1990;Tarduno et al., 2012). While the finite difference methods frequently used in these models can be applied to systems involving T-dependent properties, the heat conduction equation becomes nonlinear and more expensive to solve when T-dependent k is included (Özısık, 1993). Bulk rock conductivity decreases by 40%-60% of its value at room temperature in mantle rocks when temperature increases from room temperature to 1273 K, while conductivity increases by approximately 4% with an increase in pressure of 1 GPa (Hofmeister, 1999;Seipold, 1998;Wen et al., 2015). Due to the weaker dependence of conductivity on pressure, and the low pressure gradients expected in planetesimals, in this paper we will focus on the temperature dependence of material properties.
Previous models of planetesimal thermal evolution take various approaches to the incorporation of k, ρ, and C. These models address different stages of planetesimal evolution, depending on the meteorite group of interest, and can be broadly grouped into two classes. Models focusing on the accretion, early heating and melting of asteroids and planetoids investigate the origin of primitive meteorites (Allan & Jacobs, 1956;Elkins-Tanton et al., 2011;Hevey & Sanders, 2006), while conductive cooling models examine the post-peak-T period following recrystallization and capture the genesis of extensively differentiated meteorites such as pallasites (Bryson et al., 2015;Ghosh & McSween, 1998;Haack et al., 1990;Nichols et al., 2016;Scheinberg et al., 2016;Tarduno et al., 2012). Models in the first class, for example those investigating the ordinary chondrite parent body, often employ temperature-dependent diffusivity from Yomogida and Matsui (1983): κ = A + B/T, where A and B are terms that describe the degree of compaction of the parent body (Akridge et al., 1998;Bennett & McSween, 1996;Harrison & Grimm, 2010). Ghosh and McSween (1999) highlight the importance of incorporating a temperature-dependent specific heat capacity in the modeling of primitive asteroids, recording a decrease in peak temperatures and corresponding change in closure temperatures when T-dependent C is used, but k and ρ are held constant.
The second class of models, which address conductive cooling in differentiated planetesimals such as the primary pallasite parent body (Bryson et al., 2015;Ghosh & McSween, 1998;Nichols et al., 2016;Scheinberg et al., 2016), generally assume mantle k, ρ, and C are independent of temperature. When experimentally investigating the effect of Fe content on olivine conductivity, Zhang et al. (2019) comment on the inclusion of T-dependent and composition-dependent k in their COMSOL TM models and note that the inclusion of variable properties have a non-negligible effect on the thermal evolution of a silicate sphere. However, the focus of the study is olivine forsterite content and the impact of olivine composition on the thermal evolution of planetary bodies, and T-dependence is not systematically explored. The implications of neglecting T-dependent k, ρ, and C on the interpretation of meteorite parent body models is not understood.
Meteorites that display remnant magnetization can inform us about the magnetic field present in the environment of their parent body, which in turn allows us to estimate when an internal dynamo may have been active (Scheinberg et al., 2016). The pallasite parent body has been chosen as an example for this study as previous research has tied paleomagnetism identified in meteorite samples to the period of core crystallization in the parent body (Bryson et al., 2015Nichols et al., 2016;Tarduno et al., 2012). In order for the metal portion of a pallasite meteorite to record a convectional core dynamo, it must cool through the tetrataenite chemical ordering temperature of the metal portion while the core is crystallizing (Bryson et al., 2015;Scheinberg et al., 2016). Modifying the material properties of the body affects whether this condition is met. The geochemical and petrological heterogeneity exhibited across pallasite meteorites has been used to argue for multiple parent bodies or alternatively different environments and depths of formation within a single parent body (Boesenberg et al., 2012;McKibbin et al., 2019). Paleomagnetism places an easily testable constraint on models to investigate the importance of including T-dependent properties when deciding parent body geometry, the formation depth of pallasite meteorites, and the number of parent bodies involved in formation.
Before we address the specific example of the pallasite parent body, we outline the approach used to incorporate T-dependent properties in models of conductive cooling of planetesimals and show how, even in simple cases, this can have an important influence on their thermal history. We first address the model and numerical scheme in section 2, before exploring the sensitivity of the model to different parameters with mantle k, C, and ρ as independent of T and investigating the range of parameters used in the literature. We then address the incorporation of a non-linear term when T-dependent k is included by using a series of simple linear functions for k(T) in section 3.2. We implement T-dependent functions for k, C, and ρ in section 3.3, and attempt to recreate these results by averaging the values for k, C, and ρ radially and through time and then using these mean values in the constant model. Finally, we discuss the relevance to modeling the pallasite parent body.

Methods
To investigate the effect of including temperature-dependent properties in the thermal evolution of planetesimals, we used the 1D in radius r heat conduction equation with a non-linear term to allow for temperature dependence of k, ρ, and C (Carslaw & Jaeger, 1986;Özısık, 1993). As in Bryson et al. (2015), the layered body is composed of three primary materials: a metallic FeS core which is initially molten, a solid olivine mantle, and an insulating megaregolith layer (see Figure 1). Assuming a purely conductive mantle following magma-ocean solidification, in which convective heat transport is neglected, the temperature T in the mantle satisfies the differential heat conduction equation in spherical geometry: where t is time. The non-linear term arises due to the T-dependence of k. The insulating megaregolith layer is given a constant diffusivity lower than that of the mantle as in Bryson et al. (2015). Pressure and self-gravitation are not incorporated into the current model. The boundary and initial conditions are chosen as follows: Figure 1. Not to scale. General model setup, both before and during core solidification, displaying the functions relevant to different regions. Core radius is defined as a fraction of the total planetary radius, which includes the megaregolith layer. The megaregolith has a constant κ.
While core is molten While core freezes where r p is the planetesimal radius, r c is the core radius, T surf is the constant surface temperature, T c is the core temperature, and T init is the initial temperature, implying a homogeneous initial interior temperature distribution at t 0 ; the code can accommodate a heterogeneous initial temperature array but this is not used in this study. A Dirichlet boundary condition has been applied to the surface as in Bryson et al. (2015) instead of a radiative condition as used by Ghosh and McSween (1998), assuming the temperature at the surface of the planetesimal is constant and that of the ambient Solar Nebula. While a radiative boundary condition is a closer approximation to the physical system, a simpler fixed-temperature boundary condition has been found to produce negligible difference in inner-Solar System asteroidal models (Hevey & Sanders, 2006;Moskovitz & Gaidos, 2011).
The boundary condition for the base of the mantle depends on the core temperature. Because of our focus on the effect of T-dependent properties of the mantle, we follow the previous simplified core models of Bryson et al. (2015) and Tarduno et al. (2012) and assume the core is initially entirely liquid and vigorously convecting, and that on cooling it behaves as if it were pure iron or as an FeS mixture with eutectic composition. We discuss the implications of this simplified core model in section 4. The core temperature is updated by considering the total energy extracted across the core-mantle boundary (CMB). The energy transferred during a small time increment δt is where A c is the surface area of the core, r c is the radius of the core, and k cmb is the thermal conductivity at the base of the mantle at the CMB, that is, k CMB = k m (T (r c , t)). As E CMB = ρ c V c C c ΔT where V c is the total volume of the core and ΔT is change in temperature, the change in the core temperature in one-time increment (ΔT C ) is: The temperature at the base of the mantle is then updated by adding ΔT to the temperature at the previous timestep: The core cools as the mantle conducts heat to the surface, and is assumed to solidify when T c reaches the melting temperature of the FeS core (T l , in this case T l = 1200 K; Bryson et al., 2015). Once the core begins to freeze, the temperature is constant at T l as latent heat is extracted across the CMB. The liquid and solid fraction act identically during this process and partitioning of elements is not addressed during freezing. The core solidifies entirely once the total latent heat associated with crystallization has been extractedwhen E CMB during the solidification period exceeds E l , where the total latent heat of the core is: where m c is the mass of the core and L c the specific latent heat of fusion of the core (Bryson et al., 2015;Tarduno et al., 2012).

Numerical Implementation
We solve the conduction equation numerically for the mantle using an explicit finite difference scheme with first order differences in time and second order in space. Equation 1 can be rewritten with the temperature at radius r and time t denoted by where δt and δr are the constant timestep and radius step, and k is evaluated at   t t r T . A consequence of this discretization is that temperature dependent properties lag if evaluated at t − δt. A more accurate method is to evaluate k as: and similarly for C and ρ (Özısık, 1993). To reduce the error associated with variable k not being centered in time, we chose a sufficiently small δt such that We compared this with a selection of runs using the more accurate but computationally expensive method above for k t and Cρ t , and the differences in results were negligible.
The maximum timestep allowable for stability in the Forward-Time Central-Space (FTCS) scheme must satisfy Von Neumann stability criteria in 1D: , with the largest diffusivity of the scheme being chosen for the most restrictive conditions (Charney et al., 1950;Crank & Nicolson, 1947). For a constant spatial grid with δr = 1,000 m, δt = 1 × 10 11 s was sufficient to meet this criterion for the most restrictive cases with large κ. An adaptive grid was not used due to the first-order nature of the problem being addressed.
In order to assess accuracy, this numerical solution, with constant k, C, and ρ, was compared to the analytical solution for a sphere given by equation 6.18 in Crank (1979) with an initial uniform temperature T i and a constant surface temperature T s : where r = 0 at the center of the sphere and κ is a constant diffusivity, given by    k C (see supporting information). We also verified that we can reproduce the results of Bryson et al. (2015) when using the same input parameters.

Meteorite Formation Depth
The FeNi portion of pallasite meteorites records the cooling rate of the sample at 800 K ). This measurement is intrinsic to the meteorite sample and independent of parent body modeling. For a given cooling model, the intersection between the contour that matches the measured cooling rate of the meteorite sample and the 800 K isotherm gives a formation depth for pallasite material within the planetesimal. Then, the time when this depth passes through the tetrataenite chemical ordering temperature (593 K) and is magnetically recording can be compared to the timing of core crystallization to see if it occurs while the core is freezing, thus potentially recording core dynamo activity (Bryson et al., 2015).
To illustrate the implications of this study on the pallasite parent body, we calculate the formation depths of two pallasite meteorite samples, Imilac and Esquel, which have published cooling rates and remnant magnetization (Bryson et al., 2015;. We use the cooling rates applied by Bryson et al. (2015), calculated from cloudy-zone particle size ).

Parameter Choices for the Pallasite Parent Body
We selected parameters from previous models of planets, planetesimals, and asteroids in the literature and experimental results from geochemistry and mineral physics studies as detailed in Table 1. For many of these parameters, we have chosen both a reference value relevant to our example case of the pallasite parent body, and a range of values used in other models of differentiated planetesimals with different assumptions regarding geometry and composition. We have chosen a reference initial temperature that ensures a solid mantle that conductively cools, and a reference surface temperature that reflects the average midplane temperature of the circum-Solar disk at 2.5 AU, 1 Myr after Solar System formation (Hevey & Sanders, 2006). Reference values related to the megaregolith, the core and the boundary conditions are from Bryson et al. (2015), while mantle olivine properties have been chosen from experimental results and other planetesimal models (Su et al., 2018;Xu et al., 2004, see Table 1 for further citations). We have chosen r p = 250 km as our reference value so that paleomagnetic recording occurs while the core is crystallizing for both samples (sections 2.2 and 3).
Initially, we allowed models to run for 400 million years. We increased the run time if it did not capture the period of core solidification, for example in cases with larger radii. The core reverts to an isothermal state following the solidification period. This simplified approximation of a highly conductive metallic core is sufficient for the example application in this study, for which the post core-solidification period is not of interest.

Incorporation of Temperature Dependent Properties
In solids at low temperatures (T < θ D , the Debye temperature), heat capacity increases from zero at 0 K as C v ∼ AT 3 , where C v is specific heat capacity at a constant volume and A is a constant (Debye, 1912). At high temperatures (T > θ D ), heat capacity is weakly dependent on temperature and can be approximated with a constant value (Petit & Dulong, 1819). This results in approximately 30% increase in C in olivine over the temperature range commonly modeled for planetesimals ( Figure 2).
In electrically insulating solids such as mantle silicates, heat is primarily transferred through lattice or phonon conduction. As temperature increases, the mean energy per phonon also increases due to the change in phonon specific heat. At lower temperatures (T < θ D ), the inelastic phonon relaxation time is constant as scattering is primarily due to crystal defects or boundaries. This results in k ∝ T 3 due to the T-dependence of C (Hofmeister, 1999;Poirier, 2000). When phonon momentum exceeds a threshold at high temperatures, phonon-phonon Umklapp scattering acts to reduce k, producing a  1 k T dependency (Poirier, 2000). This non-monotonic behavior is illustrated for olivine in Figure 2.
A change in density with temperature can be linked to thermal expansion by the coefficient of expansivity, α: ρ = ρ 0 − αρ 0 (T − T 0 ), where ρ 0 is a reference density at T 0 , commonly room temperature (∼ 295 K). Density is less temperature dependent than C or k, and is combined with heat capacity in Figure 2 as volumetric heat capacity, both as a constant and as a T-dependent function to illustrate the scale of its effect.
In order to fully understand the effect of including temperature dependence in our model, we constructed a simple linear function for conductivity before investigating the more complex equation based on experimental results (Equation 14): where k 0 is a reference conductivity at 0 K and β controls the temperature dependence, and can be set as positive or negative. β and k 0 must be chosen such that k does not become negative over the temperatures explored in the body. In order to contrast a T-dependent conductivity with simply setting the average conductivity higher or lower, functions with both positive and negative β were chosen to approximate the same mean conductivity over radius and time. Additionally, the cases were run with and without the non-linear term. Both ρ and C were held constant to isolate the effect of the conductivity. The megaregolith layer maintains a constant κ for all model runs including those with fully variable k, ρ, and C, as after initial rapid equilibration with the surface temperature, this layer has a constant temperature. The core properties have also been kept constant.
For this study, we have chosen the function used for heat capacity in olivine from Su et al. (2018), based on lattice vibration theory from Berman and Brown (1985) and fit to experimental data from Isaak (1992): 1343 2.887 10 6.166 10 995.1 .
note that this is valid for the range of temperatures T surf -T init . We do not explore temperatures close to 0 K. The expression for thermal expansivity is also taken from Su et al. (2018) based on the functional fit by Fei (2013) and using experimental data from Suzuki (1975): As the lower temperatures modeled (∼250 K) are rarely of interest in terrestrial mineral physics and are less accessible to experimental studies, we constructed a simple conductivity function for olivine spanning 250-1800 K. As discussed above, conductivity is controlled by different processes at high-and low-temperatures, resulting in different temperature dependencies. For the high-T region, we used the experimentally derived curve from Xu et al. (2004): where a = 0.032 GPa −1 (experimentally derived) and P = 4 GPa. As T-dependence of k at temperatures ≪θ D is similar to that of C, a function identical in shape to Equation 11 but normalized such that C = 1 at T > θ D was used for the low-T region. As this low-T curve is constant and equal to one above θ D , it can be multiplied by Equation 13 to fill in the low-T region without altering the higher-T experimental results. Our resultant function is differentiable and non-monotonic: .
While the pressures inside the planetesimal are ≪ 4 GPa, changing pressure to < 1 GPa in Equation 13 increases conductivity in our composite function by < 0.3 W m −1 K −1 at all temperatures. As this is outside of the calibration range of the experiments by Xu et al. (2004) we have chosen not to include this adjustment as it may not be physically realistic and pressure effects are not the focus of this study, and instead use a and P as quoted by Xu et al. (2004). These functions are illustrated in Figure 2.

Results
The model produces arrays of temperature and cooling rate through time and radius. For any radius r, the linear, geometric, and non-linear (if applicable) terms of the heat conduction equation can be plotted against time. Model outputs that are important to the interpretation of meteorites include the initiation and duration of core crystallization, the depth within the parent body from which the meteorite was derived and when this occurred, and the peak cooling rates reached. In the specific case of the pallasite parent body, the calculated depth of formation can then be tracked to see if this region of the parent body passes through the temperature where magnetism is recorded while the core is solidifying, thus potentially recording core dynamo activity.

Constant k, ρ, and C
The model was run with constant k, ρ, and C for both the reference parameters in Table 1 and the end-member values quoted, if applicable. In addition, parameters were varied by ±10% of the reference value to gauge the sensitivity of the model to different inputs. The full results of these parameter explorations are tabulated in the supplementary information.  the temperature anomaly propagates through the mantle to deeper regions with a time delay determined by the diffusion timescale.
The slope of T(r) from the base of the mantle to the surface is negative for the duration of the model run.
Initially, T(r) is convex upwards but flattens over time and becomes convex downwards as the linear term changes sign: initially within the body   2 2 T r is negative for all radii and increases with time, becoming positive at the boundaries first, with this change in sign propagating towards the middle of the mantle. When the core is removed to approximate a solid sphere, this effect is only seen to propagate downwards from the surface boundary as the breaking effect of the core on the cooling of the mantle is not present. The geometric term then drives further cooling after this point (Figure 3).
When the core reaches 1200 K and begins to freeze, the temperature at the CMB is held constant. The fixed core temperature reduces the cooling rate in the mantle sharply; in the deeper regions of the mantle    T t drops towards zero as the mantle reaches the same temperature as the core. The effect is less pronounced in the shallow regions as the cooling rate has already slowed significantly and is approaching zero.
The body cools rapidly at the surface, with shallow depths quickly equilibrating with the constant surface temperature (Figure 4). High temperatures are maintained for longer deeper within the body due to the overlying insulating mantle. Using the cooling rates applied by Bryson et al. (2015), calculated from cloudyzone particle size  we calculated source depths of 64 km for Esquel and 57 km for Imilac, approximately midway through the mantle (Figure 4 and Table 2).
The geometry of the body is a strong controlling factor on the cooling rate and timing of core crystallization ( Table 2). The planetary radius has the largest effect: increasing the total radius by 10% slows the cooling of the planetesimal at depth and delays the onset of core crystallization by 38 Myr. When the core fraction is increased by 10%, the core begins to freeze 5 Myr earlier as there is less insulating mantle, but takes 4 Myr longer to freeze fully due to its increased size. While the average cooling rate of the body drops sharply for all cases on initiation of core solidification, the effect is more pronounced when the core fraction is increased as the cooling rate of the core dominates the overall cooling rate. Increasing the insulating megaregolith thickness by 1 km while maintaining a 250 km total radius does not delay the onset of core crystallization, but does increase the duration of the solidification period by 1 Myr. Increasing the megaregolith thickness further does delay core solidification, with a 20 km thick megaregolith causing a 73 Myr delay when compared to the reference case (see supplementary information). The resulting changes in the calculated source region depth for pallasite meteorites is also shown in Table 2.
Increasing k by 10% accelerates the cooling in the body, causing the core to begin solidifying 15 Myr earlier.
Increasing ρ or C by 10% has the opposite effect, and delays the onset of core crystallization by 8 Myr. Table 2 also shows the results of setting k = 4 W m −1 K −1 and 1 W m −1 K −1 , which reflect the end-member expected values if k varied with T (see Figure 2). Between these two cases, there is a 198 Myr difference in the timing of the start of core solidification. The core begins to freeze at 132 Myr and the freezing period lasts 53 Myr when k = 4 W m −1 K −1 , while the core begins to freeze at 330 Myr when k = 1 W m −1 K −1 . An increase in conductivity results in deeper source regions for the pallasite meteorites, with the Esquel and Imilac source regions moving 13 and 10 km deeper respectively when k = 4 W m −1 K −1 , while both move ∼22 km shallower when k = 1 W m −1 K −1 .

Simple Linear T-Dependent Conductivity
In this section we explore k(T) in the form k = k 0 + (βT) with ρ and C held constant. For the examples shown in Figure 5 and summarized in Table 3, we chose β = ±0.0025 W m −1 K −2 and k 0 such that k = 3.0 W m −1 K −1 at the mean temperature of the reference case with constant k, ρ, and C (with megaregolith thickness set to 0 km - Table 3) to isolate the effect of T-dependence. The model was run both with and without the non-linear term in Figures 5a and 5b. When compared to the constant case with k = 3 W m −1 K −1 , allowing k to vary with T changes the timing and duration of the core crystallization period (see Table 3). For β = 0.0025 W m −1 K −2 and k 0 = 1.1125 W m −1 K −1 (panel (a), Figure 5), the onset of core crystallization is 19 Myr earlier than for the constant case (Table 3); in the early stages of the model run the average cooling rate throughout the body is higher than the constant case due to higher initial conductivity in the mantle (panel (c) of Figure 5). After ∼80 Myr (before the core begins to freeze), the average cooling rate throughout the body drops below the constant case, resulting in a 3 Myr longer core-crystallization period. The duration of core crystallization is close to that of the constant case as, during this time period, the variable conductivity is similar to the fixed conductivity of the constant case (panel (c), Figure 3).
When the nonlinear term is neglected (panel (b), Figure 5), core crystallization initiates 46 Myr earlier than in the constant reference case, due to increased cooling rates despite a lower average conductivity. The nonlinear term is always positive and slows cooling if β > 0, reducing the peak cooling rates experienced at this depth and the average cooling rates in the mantle.  The calculated source depth of the Imilac and Esquel meteorites for this model setup are shown in both plots, using the cooling rates applied by Bryson et al. (2015), calculated from cloudy-zone particle size . Temperature contours highlight the tetrataenite formation temperature when paleomagnetism can be recorded (593 K) and the temperature for which the sample's cooling rates were measured (800 K), while cooling rate contours show the measured cooling rates for both samples. the mantle with time. The nonlinear term in this case is negative, owing to the negative sign of dK dt , and it amplifies the initial peak cooling rates at the depth examined (panel (d), Figure 5); however, the overall average cooling rate of the body is initially lower due to the low conductivity (Figure 5f). When the nonlinear MURPHY QUINLAN ET AL.   term is neglected, the core begins to solidfy 146 Myr later than in the constant case, and solidification takes 24 Myr longer. As the core does not freeze at the midpoint between the initial and surface temperatures, the nonlinear terms for positive and negative β are not symmetric.
In summary, positive β leads to earlier onset of core freezing and a longer duration of core freezing, while negative β results in later onset of freezing and a shorter freezing period. For both ± β the change in onset time when compared to the constant case is much larger than the change in the duration of core freezing, as there is a much greater difference between constant and variable k earlier in the model than during core solidification (Figures 5c and 5d). Even for linear conductivity functions with shallow slopes, the conductivity structure of the mantle is very different to that of the constant case and the temporal dependence of this structure has implications for the timing of events within the body that cannot be approximated by changing the value of k in the constant case. Inclusion of the nonlinear term is essential as neglecting it can result in large over-or under-estimations of core crystallization onset time (for negative β, neglecting the nonlinear term results in 119 Myr delay in the onset of core crystallization). The implications of these results on the pallasite parent body are investigated using the experimentally derived functions in the next section.

Temperature-Dependent Properties: Using Experimental Functions
The fully variable case, using the default parameters in Table 1 and the k(T), C(T), and ρ(T) functions (Equations 11, 12 and 14), resulted in a 40 Myr delay in the onset of core crystallization but only 3 Myr longer period of core crystallization when compared to the reference case with constant properties (Figure 6). The temperature distribution in the shallow mantle is similar to that of the constant reference case, but the interior stays hotter for longer when T-dependent properties are used ( Figure 6). The fully variable case requires deeper source regions for the pallasite meteorite samples than the reference case, with a depth of 61 km calculated for Imilac and 68 km for Esquel (Table 4).
When discussing simple linear functions for k(T), we have demonstrated that cases with constant and variable properties should be correctly calibrated in order to make meaningful comparisons. In order to do so, we measured the average temperature in the mantle of the fully variable case and used this to calculate new constant values of k, C, and ρ using Equations 11, 12 and 14. The mean temperature of the mantle over the 400 Myr of the model lifetime was 780 K, giving k = 2.8 W m −1 K −1 , ρ = 2,945 kg m −3 , and C = 996 J kg −1 K −1 . The model with constant properties was then rerun with these updated values for k, ρ, and C, to more closely approximate the results from the fully variable model. In this section, this new model using updated constant k, ρ, and C is referred to as the constant mean values case, and the results are shown in Table 4.
In the fully variable case (Figure 7), the nonlinear term is negative and enhances the overall cooling rate at the depths displayed for all times shown (up to 400 Myr), as the slope of the function for k is negative MURPHY QUINLAN ET AL.  for all T > 300 K (Figure 2). A thin insulating layer in the shallow mantle forms where T < 300 K and the nonlinear term is positive. The core begins to freeze 211 Myr after model initiation, and takes 61 Myr to fully solidify. The constant mean values case does not replicate this result: with constant k, ρ, and C, the core begins to solidify at 189 Myr and takes 53 Myr to fully freeze (Table 4). In addition, the constant mean values case requires shallower source regions for the pallasite meteorites Imilac and Esquel: 53 and 60 km respectively (Table 4). Qualitatively, the fully variable case is similar to the case with linear k and negative β in section 3.2: the core begins to freeze later but takes a shorter time to fully crystallize than the constant mean values case (Tables 3 and 4). However, the insulating layer in the shallow mantle with a positive nonlinear term cannot be replicated by the simple linear case and so the fully variable case must be used for quantitative results. When the non-linear term is set to zero, again the fully variable model behaves similarly to the β < 0 linear case (Table 4).
When the different properties are allowed to vary in turn, T-dependent C produces the smallest deviation in core crystallization timing from the constant mean values case, as at high T (temperatures such as those experienced by the planetesimals prior to and during core crystallization), C is approximately constant MURPHY QUINLAN ET AL.

10.1029/2020JE006726
13 of 18 The calculated source depth of the Imilac and Esquel meteorites for this model setup are shown in both plots, using the cooling rates applied by Bryson et al. (2015), calculated from cloudy-zone particle size . Temperature contours highlight the tetrataenite formation temperature when paleomagnetism can be recorded (593 K) and the temperature that corresponds to the sample's measured cooling rates (800 K), while cooling rate contours show the measured cooling rates for both samples.
( Figure 2). Including variable ρ results in a 9 Myr delay in the onset and 2 Myr longer duration of core crystallization in comparison to the constant mean values case, while including only variable k results in an 11 Myr delay in the onset and a 4 Myr shorter duration of core crystallization. Variable ρ produces the shallowest meteorite source regions of the three properties while variable k produces the deepest (Table 4). Including just one T-dependent property cannot replicate the fully variable model.

Discussion and Conclusion
Including T-dependent thermal properties changes the temperature structure in the modeled planetesimal: predictions of mantle temperature can differ by 50 K over tens of millions of years even when the best estimates for constant k, ρ, and C are used ( Figure 8). This results in significant changes in the timing and duration of core crystallization: the onset of core solidification is 22 Myr later, a delay of 12%, while the core solidified 3% faster. The delay in onset of core crystallization is equivalent to increasing the radius of the planetesimal by 10% with constant parameters, but increasing r p extends the period of solidification by 13% (Table 2). We use the example of a pallasite parent body to illustrate these results: including T-dependent properties delays the onset of core crystallization and results in deeper source regions for pallasite meteorites than when constant k, ρ, and C are used (Figure 8). In this example, T-dependent k, ρ, and C result in a hotter deep mantle but cooler shallow mantle, which cannot be replicated by constant values (Figure 8).
Including T-dependent properties also affects whether or not samples are predicted to preserve remnant magnetization from a core dynamo: while in the constant reference case, both the Imilac and Esquel meteorite source depths cool through 593 K during core solidification, the Imilac region cools down below 593 K before core solidification when variable k, ρ, and C or mean constant values based on the variable case are used ( Note. Summary of key results. Timing of core crystallization period given in millions of years after model start (Myr) and formation depth of meteorites given in km. a Reference case with constant k = 3 W m −1 K −1 , ρ = 3,341 kg m −3 , and C = 819 J kg −1 K −1 . b Constant case here differs from the reference case: values for k, ρ, and C are calculated at the mean T in the fully variable case: k = 2.8 W m −1 K −1 , ρ = 2,945 kg m −3 , and C = 996 J ,kg −1 K −1 . c Case with T-dependent k, ρ, and C.
d T-dependent properties, but with nonlinear term neglected. e One property allowed to vary with T with other properties held at mean values as in ( b ). Table 4 Variable k, ρ, and C Figure 7. Results for the reference case with T-dependent k, ρ, and C. The components of the heat conduction equation are shown at a depth of (a) 42 km (one third of the thickness of the mantle) and (b) 84 km (two thirds). The cooling rate is multiplied by −1 to illustrate how it balances the other components to add to zero. The green area defines the period of core crystallization when T-dependent properties are used, while the pink area highlights the period of core crystallization from the mean constant case for comparison. constant mean case for this example, the input values for k, ρ, and C require the fully variable case to be run initially in order to be calculated.
In our example of a 250 km radius parent body, Imilac forms only ∼5 Myr before the core begins to crystallize and so can be accounted for by error in the measurement of the cooling rate from this sample (Bryson et al., 2015;. However, larger discrepancies in timing can be found for different cooling rates, parent body radii, megaregolith thickness or core fraction (Figure 9). Including T-dependent properties narrows the range of input parameters that allow meteorite samples to potentially record paleomagnetic signatures. This provides a simple criterion for testing different parameter combinations: whether the meteorite source region cools through the tetrataenite chemical ordering temperature during core solidification. As shown in Figure 9, when constant k, ρ, and C are used, megaregolith thicknesses anywhere between 0-12 km satisfy the above criteria for a planetesimal of 250 km radius and a core that is 50% of r p , while a megaregolith layer of 4-8 km is required when T-dependent properties are used. If the core fraction is reduced to 30% of r p , a 250 km body with megaregolith between 0-8 km can accommodate both meteorite samples, whereas no suitable combination of parameters can be chosen when T-dependent k, ρ, and C are used. Similarly, no suitable parent body with a 250 km radius and a core fraction of 70% r p can be found if T-dependent properties are used, whereas if these values are taken as constant, then a planetesimal with a radius of 300 km including an 8 km thick megaregolith can produce the cooling rates and required timings in both meteorites. Nichols et al. (2016) find that two additional pallasites, Marjalahti, and Brenham, record a weak magnetic field and argue that these samples cooled through the tetrataenite formation temperature before the onset of core crystallization. This timing stipulation could provide an additional constraint on the allowable physical parameters in the model. However, for the range of parameters explored in Figure 9, Marjalahti and Brenham record a magnetic remanence before core crystallization for all cases except a selection already ruled out by Esquel and Imilac recording a remanence after core crystallization. Therefore, in this case they do not provide an additional constraint on the timing of core crystallization, but may be useful for different parameter searches.
One limitation of this work comes from the simplified approach to modeling the core. We assumed a dynamo is generated during bottom-up eutectic solidification and have neglected T-dependent properties in this region. In reality, core solidification in planetesimals is likely to be complex and strongly dependent on the bulk sulfur content of the core. Bulk S content is difficult to estimate from iron meteorite samples due to the incompatibility of S in solid iron, and the resulting low-S composition of these samples . Bryson et al. (2019) predict a core dynamo driven by compositional convection at earlier times than suggested by the model we present, while the initially non-eutectic core composition evolves towards a eutectic composition. They argue that eutectic solidification would not be expected to generate a core dynamo, as sulfur would not be expelled during solidification as the inner core adopts an FeS composition. Within this framework, Esquel and Imilac are instead predicted to experience a magnetic field in the period "before core solidification" in Figure 9, but this period of non-eutectic solidification cannot be easily quantified without a fuller treatment of the core, which is not warranted by our focus on the importance of temperature dependent conductivity on the cooling of the mantle. Furthermore, core solidification fronts may initiate at the core-mantle boundary, resulting in top-down solidification through dendritic growth (Scheinberg et al., 2016). Top-down crystallization would allow a dynamo to be generated during eutectic solidification, and both modes of solidification have been inferred for differentiated planetesimals based on MURPHY QUINLAN ET AL.  k, ρ, and C through time and radius are equal. Period of core crystallization is shown in dashed white for the constant cases, and in green for the variable case. Symbols mark the source regions for the Imilac and Esquel meteorites as they pass through the 593 K isotherm; white circles show the results from the constant cases, while green shows the result when variable properties are used. We use the cooling rates applied by Bryson et al. (2015), calculated from cloudy-zone particle size . iron meteorite cooling rates Yang et al., 2008;Yang, Goldstein, Michael, et al., 2010).
Following crystallization, the core is assumed to return to an isothermal state due to the high conductivity of the material. For the pallasite example case, this is an acceptable simplification as it is the times preceding and during the core solidification period that are of interest. For other applications it may be required to restart the model with the core included in the iterative solution with a Neumann boundary condition at the center, as used for approximating the analytical solution. The effects of pressure and gravity have also been neglected due to the low pressure gradient expected within the body as discussed in section 1.
In conclusion, T-dependent properties can significantly impact the output of planetesimal cooling models, even if the model results are being used qualitatively or to judge the relative timing of processes within the body, such as whether meteorite formation regions cool through specific temperatures before, during or after the period of core crystallization. The inclusion of T-dependent k, ρ, and C results in later crystallization of the core (∼40 Myr later than the constant reference case and ∼20 Myr later than the updated constant case) and deeper meteorite formation depths due to suppressed cooling rates in the mantle. This result cannot be replicated with constant values for k, ρ, and C, even when these values are chosen to match the mean values of each through time and radius in the variable model. If T-dependent κ is included without a nonlinear term, the reduction in cooling rates through the body is overestimated, resulting in core solidification 33 Myr after the variable case and 73 Myr after the constant case. These results are shown with relevance to the pallasite parent body. The parameter space which satisfies the cooling rate criteria for the material which formed the Imilac and Esquel meteorites shrinks when T-dependent mantle properties are included; it follows that if more samples are investigated the parameter space will shrink further. Future work could use this more restrictive parameter space to address the ongoing debate over the number of required pallasite parent bodies and potentially place a minimum constraint on the number of bodies required. T-dependent properties should also be addressed for other planetesimals and meteorite parent bodies where conduction is involved, for example the ordinary chondrite parent body, where peak temperatures and the inferred parent body radius may be incorrectly calculated. Figure 9. Planetary radius, core size, and megaregolith thickness investigation for the constant k, ρ, and C case, and the fully variable case. The color and symbol denote whether or not the Imilac and Esquel meteorite source region cooled through 593 K during core crystallization ±10 Myr: green triangles mark models where this criteria was met. Red crosses denote models where the meteorite cooled through 593 K after core crystallization, whereas blue squares show where this happened before the core began to crystallize. Gray markers note that no matches for the meteorite cooling rates at 800 K were found, implying the meteorite could not have formed in that body. Where both samples have different results, Imilac is shown on the left and Esquel on the right. We use the cooling rates applied by Bryson et al. (2015), calculated from cloudy-zone particle size .