© ESO, 2014

1. Introduction

Investigations of the structure of chondritic meteorites reveal significant variations in the average porosity between samples (Wilkinson et al. 2003), suggesting that they formed at different depths in the respective parent body. Furthermore, the average porosities of less than 27% (0–27% for H chondrites, 0–14% for L chondrites, and 0–18% for LL chondrites, see Wilkinson et al. 2003) are rather small compared to the average porosity of a model structure with a random loose or a random close packing of grains. Random packings should actually be expected for objects that accreted from dust in the protoplanetary nebula and that did not experience any kind of compaction, i.e. values of the porosity range between 36% and 44% (Dullien 1992).

The measured average porosities of the meteorites, however, indicate that processes such as compaction and metamorphism have modified the parent body. Compaction, i.e. the decrease of the porosity below the equilibrium value (see below), can be caused by mechanical processes (such as stress, plastic or brittle deformation, and fracture evolution) and/or by geochemical processes (such as dissolution and reprecipitation). Thus, the porosity of the random packings can be considered as initial porosity of the superficial material that accumulates at the surface during the accretion of a planetesimal.

Assuming random packing of particles that accreted at the surface of a planetary object, the porosity of that material depends on its structural parameters. The most important parameters are the grain size, the packing, the particle shape and the grain size distribution. The average or total porosity of an object is defined as the ratio of pore volume to the total bulk volume and can be expressed as a fraction or as a percentage. We note that the fraction value does not contain any information about distribution, sizes, or connectivity of the pores. Thus, objects with equal total porosity can have different physical properties. The porosity of randomly packed material is in equilibrium if it cannot be changed by tapping and shaking and without applying any pressure or increasing temperature. The value at which such an equilibrium state is achieved depends on the frictional and cohesive forces between the grains that are proportional to the exposed surface area. Per unit volume, the exposed surface area is inversely proportional to the grain size. Thus, coarse grain structures stabilise at a lower porosity than finer grain structures assuming the same mass. If the grain size increases above 50 μm, the frictional and cohesive forces decrease and a limit of the porosity is reached at a value of φ ≈ 0.399; for larger grain sizes this limit remains constant. This value represents a lower boundary for the porosity of randomly packed spheres (due to shaking and tapping only). Once non-spherical grains are considered, the porosity of the non-frictional limit changes, i.e. random packings of angular grains (e.g. cubes, cylinders, discs) attain equilibrium porosities higher than 0.399. However, a mixture of two grain sizes, for example, leads to smaller porosities than those of the samples with the involved constant grain sizes (see e.g. Coogan & Manus 1975). In contrast, the equilibrium porosity of ordered packings is independent of the grain size and depends simply on the choice of the packing. In the present study, we will use ordered packing based on the work by Kakar & Chaklader (1967) as this allows us to model self-consistently the initial porosity and the compaction process for selected packings.

In general, compaction can be considered as a two-stage process. Owing to the radial growth during the accretion, accreted particles are subsequently covered by fresh material and experience increasing lithostatic pressure. This means that the randomly or ordered packed grains can become closer and compact to the limit of φ ≈ 0.4 (see e.g. Guettler et al. 2009). This is a purely self-gravitational effect that does not require any temperature increase above the nebula temperature, and is referred to as cold pressing. Further compaction below the porosity of φ ≈ 0.4 is due to the so-called hot pressing or sintering, which can be described by plastic deformation of the grains by thermally activated creep processes. It depends on the temperature-dependent viscosity of the material and the effective stress acting on a dust grain. The effective stress depends on the geometric packing of the particles and the lithostatic pressure. Thus, compaction is in particular governed by temperature and stress. Because the pressure and thus the effective stress in a planetesimal increases with depth, shallow layers require a higher temperature than deep layers to compact by the same degree.

In the literature, the process of porosity loss by sintering in planetesimals is commonly modelled using a parametrised approach, for which compaction occurs at the fixed temperature interval of 670–700 K (e.g. Hevey & Sanders 2006; Sahijpal et al. 2007; Gupta & Sahijpal 2010; Solano et al. 2012). The approach dates back to a paper by Yomogida & Matsui (1984). In that work, the experimental results from Schwenn & Goetze (1978) were combined with a geometric approach of the calculation of the effective stress on dust grains from Kakar & Chaklader (1967) and Rao & Chaklader (1972) to model self-consistently compaction of planetesimals with radii of ≤100 km. The hot pressing experiments by Schwenn & Goetze (1978) with olivine have been reproduced to model the evolution of log (1 − φ) (where φ is the porosity; see also Eq. (9)) using an activation energy of E = 85 ± 29 kcal mol^-1 (E = 356 ± 121 kJ mol^-1). In their numerical simulations, however, Yomogida & Matsui (1984) used E = 45 kcal mol^-1 (E = 188 kJ mol^-1) for temperatures below 700 K. This value was chosen in order to obtain compaction at low temperatures to reconcile the meteoritic evidence of compaction even with inefficient heating by long-lived isotopes ⁴⁰K, ²³²Th, ²³⁵U, and ²³⁸U. Heating by ²⁶Al and ⁶⁰Fe was not considered in their study and most bodies considered by Yomogida & Matsui (1984) had maximum temperatures of ≤700 K. Yomogida & Matsui (1984) concluded that the sintering temperature, i.e. the temperature at which sintering starts to become important, is around 600−650 K. In this manner they artificially favoured compaction at low temperatures. The value of E = 45 kcal mol^-1 (E = 188 kJ mol^-1) lies clearly outside of the interval given by Schwenn & Goetze (1978) for olivine dominated material and is too low for an H-chondritic composition. Typical values of the activation energy for creep processes are 295−530 kJ mol^-1 in olivine dominated mantle rocks (Karato 2013), 270−720 kJ mol^-1 for enstatite powders (Bystricky et al. 2011; Karato 2013), and 49−130 kJ mol^-1 for ice (Goldsby & Kohlstedt 2001; Durham et al. 2001).

Considering the assumptions in Yomogida & Matsui (1984), their conclusions about the sintering temperature are valid. However, their findings have been simplified and generalised in other studies. Hevey & Sanders (2006) referred to Yomogida & Matsui (1984) and incorporated compaction at 700 K in their simulations. Later, Sahijpal et al. (2007) and Gupta & Sahijpal (2010) included a gradual reduction of porosity in the temperature range of 670−700 K in their models and Solano et al. (2012) used the same approach to model compaction of Vesta. Apart from the models by Henke et al. (2012) and Neumann et al. (2012), who consider compaction of chondritic planetesimals as a thermally activated creep process (and possibly the calculations of Akridge et al. 1997; Senshu 2006), most studies use the simplified parametrised approach based on Yomogida & Matsui (1984) and allow compaction in a rather small temperature window at ≈700 K (e.g. Akridge et al. 1998; Hevey & Sanders 2006; Sahijpal et al. 2007; Gupta & Sahijpal 2010) to study the insulating effect of regolith on the thermal evolution. Other studies which do not consider compaction explicitly, but discuss its effects on the results, still refer to some of the publications mentioned above (e.g. Elkins-Tanton et al. 2011; Moskovitz & Gaidos 2011).

In a previous study (Neumann et al. 2012), we investigated compaction of planetesimals with radii of ≤120 km by hot pressing using a creep-related approach combined with accretion. We compared the thermal evolution of initially consolidated bodies with that of initially porous ones and concluded that the maximum temperature in initially highly porous planetesimals is significantly higher than in consolidated bodies. This effect is particularly pronounced (temperature increase by up to 300%) for small planetesimals with radii of less than a few km. We also concluded that compaction starts at temperatures around 600 K (at high pressures) and could proceed to the temperatures above 1200 K (at low pressures). These findings are consistent with the results of Henke et al. (2012). Nevertheless, a systematic study is lacking that examines whether compaction achieved by the reduction of the porosity in a fixed temperature interval is an acceptable approximation and that discusses the implications of this approach.

In the present study, we compare two approaches to model compaction of planetesimals via hot pressing. The first one is based on the plastic deformation of dust grains and is referred to as the “creep related approach” in the following. The second one is based on the linear reduction of the porosity from its initial value to zero in the temperature interval 650−700 K, and is referred to as the parametrised approach. We show that the parametrised approach overestimates compaction by hot pressing due to the usage of a specific temperature interval. Furthermore, we investigate in which way ordered packing of grains, grain size, and activation energy influence the compaction process.

2. Model

We use the model of thermal evolution and differentiation of planetesimals described in Neumann et al. (2012) and Neumann et al. (2013). The energy balance is solved in a spherically symmetric planetesimal, using a 1D heat conduction equation where density, thermal conductivity, heat source density, and radius are functions of the porosity. The body is heated by ²⁶Al and ⁶⁰Fe and is initially porous (we note that in the current study we adopt a different abundance of ⁶⁰Fe than in our previous papers; see Table 2). The porosity changes with time because of compaction via hot pressing. The hot pressing is implemented based on the experimental results from Schwenn & Goetze (1978) and on the geometric deformation theory from Kakar & Chaklader (1967) and Rao & Chaklader (1972) (also utilised by Yomogida & Matsui 1984). In Neumann et al. (2012, 2013) we considered one possible (orthorhombic) packing of equally sized spheres and a constant initial porosity of φ₀ = 0.4 for the computation of the effective stress between dust grains in a compacting system.

A similar model, but with an initial porosity of φ₀ = 0.5, was used by Henke et al. (2012). In the present study, we consider four basic systems of packing in order to compare the resulting thermal evolution and porosity distribution with the ones from the parametrised approach, i.e. compaction occurs in the temperature interval 650 − 700 K. While the thermal evolution including the dependence of density, thermal conductivity, heat source density, and radius on porosity and the setup of the numerical model was described thoroughly in Neumann et al. (2012, 2013), we describe here in more detail the theoretical basis on the modelling of the porosity loss.

2.1. Theoretical approach to the modelling of compaction

One of the most common forms of plastic flow is power law creep, which describes the relation between the steady-state strain rate dε/dt and the effective stress during hot pressing, $\frac{\mathrm{d}\mathit{\epsilon }}{\mathrm{d}\mathit{t}}\mathrm{=}\mathit{C}{\mathit{\sigma }}_{\mathrm{1}}^{\mathit{n}}\mathit{,}$(1)where σ₁ is the effective stress acting on the material during creep, the stress exponent n is a material constant, the factor C = C(T) is a function of temperature, and t is time. The strain rate during densification of powder can be given in terms of the rate of change of the relative density D = 1 − φ (or of the porosity φ) as $\frac{\mathrm{d}\mathit{\epsilon }}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathrm{d}\mathit{l}}{\mathit{l}\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{1}}{\mathit{D}}\frac{\mathrm{d}\mathit{D}}{\mathrm{d}\mathit{t}}\mathrm{·}$(2)This is an approximation of the strain rate during densification of a powder compact of the height l in a die and can be applied in the first approximation to small planetary bodies. By setting $\frac{\mathrm{1}}{\mathit{D}}\frac{\mathrm{d}\mathit{D}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathit{C}{\mathit{\sigma }}_{\mathrm{1}}^{\mathit{n}}$(3)and $\frac{\mathrm{1}}{\mathit{D}}\frac{\mathrm{d}\mathit{D}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{\partial }\log \mathrm{(}\mathrm{1}\mathrm{-}\mathit{\phi }\mathrm{)}}{\mathit{\partial t}}$(4)we obtain $\frac{\mathit{\partial }\log \mathrm{(}\mathrm{1}\mathrm{-}\mathit{\phi }\mathrm{)}}{\mathit{\partial t}}\mathrm{=}\mathit{C}{\mathit{\sigma }}_{\mathrm{1}}^{\mathit{n}}\mathit{,}$(5)where C and n need to be determined experimentally and σ₁ needs to be computed with respect to the geometry of the material on a small scale. Thereby, σ₁ changes with the relative density during compaction. It is widely accepted that the effective stress during hot pressing is not equal to the applied stress σ₀ and depends on the porosity φ. Often, an inversely linear dependence σ₁ = σ₀D^-1 on the relative density is used, suggested by Farnsworth & Coble (1966) and so any specific intrinsic geometry of the material is neglected.

Fig. 1 Relative stress σ1/σ0 for the adopted four basic systems of packing plotted against the relative density D = 1 − φ (after Rao & Chaklader 1972, their Fig. 6). The lines are truncated at the respective critical density Dc (see Table 1), above which the effective stress is computed using a linear relationship between the value at Dc and the lithostatic pressure.

The intrinsic geometry is accounted for in the approach of Rao & Chaklader (1972), ${\mathit{\sigma }}_{\mathrm{1}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{0}}{\mathit{\alpha }}_{\mathrm{1}}^{-1}{\mathrm{(}{\mathit{D}}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{\beta }}_{\mathrm{1}}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{-1}\mathit{,}$(6)which we adopt in the present study. Here, σ₀ is the lithostatic pressure inside a planetesimal, the parameters α₁ and β₁ are constant and reflect the assumptions on the packing of dust grains, and r is the numerical grain radius (varying between r₀ and r_c; see Table 1 and Appendix A for more details). Figure 1 shows the relative stress σ₁/σ₀ as a function of the relative density D calculated according to Eq. (6). Prior to compaction (having a high porosity and a low relative density), the value of the effective stress σ₁ is several orders of magnitude higher than that of the lithostatic pressure. During the ongoing consolidation of the material, the relative stress approaches the value 1 meaning that the effective stress converges against the lithostatic pressure if the relative density approaches 1 (or the porosity approaches 0).

In the study by Kakar & Chaklader (1967) four basic systems of packing are considered, defined by the geometric arrangement of spherical grains and by the coordination number (i.e. the number Z of the contact points, an average value for equally sized spheres). In the simple cubic packing (SCP) every grain has four neighbours in a horizontal layer and two additional neighbours in the layers above and below. Orthorhombic (ORP) packing corresponds to six neighbours in a layer (plus two above and below). The rhombohedral packing (RHP) has coordination number 12 and the body-centred cubic packing (BCCP) 8. A geometric theory was introduced by Kakar & Chaklader (1967) in order to treat the deformation of the spheres involved in any of these four packings and thus the loss of porosity in a dust sample. We note that according to this theory, the deformation that corresponds to BCCP is a two-stage process. We refer to the two stages as BCCP1 and BCCP2 (see Appendix A for the detailed information). Based on the above model, Rao & Chaklader (1972) introduced a method for the computation of the effective stress acting on a dust particle (see Eq. (6)). Combination of the effective stress from Eq. (6) with Eq. (5) yields $\frac{\mathit{\partial }\log \mathrm{(}\mathrm{1}\mathrm{-}\mathit{\phi }\mathrm{)}}{\mathit{\partial t}}\mathrm{=}\mathit{C}{\left(\frac{{\mathit{\sigma }}_{\mathrm{0}}}{{\mathit{\alpha }}_{\mathrm{1}}\left({\mathit{D}}^{\frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\beta }}_{\mathrm{1}}^{\frac{\mathrm{2}}{\mathrm{3}}}{\mathit{r}}^{\mathrm{2}}\mathrm{-}\mathrm{1}\right)}\right)}^{\mathit{n}}\mathrm{·}$(7)There are two parameters C = C(T) and n which depend on the temperature and on the material properties (such as activation energy). Experimental data for olivine is available from the laboratory experiments of creep of olivine during hot-pressing performed by Schwenn & Goetze (1978). They provide a fit obtained for spherical olivine grains with radii ≤50 μm: $\mathit{C}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\mathit{B}{\mathit{b}}^{\mathrm{-}\mathit{m}}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{E}}{\mathrm{ℛ}\mathit{T}}}\mathrm{,}\mathit{n}\mathrm{\approx }\frac{\mathrm{3}}{\mathrm{2}}\mathrm{·}$(8)Here, B is a constant factor between 5.06 × 10^-19 and 1.71 × 10^-18 m³Pa^−3/2s^-1, m is the grain size exponent, E is the activation energy in kJ mol^-1, R is the gas constant in kJ mol^-1 K^-1, and b is the grain radius in m (see Table 2 for the adopted values). Consequently, we obtain $\frac{\mathit{\partial }\log \mathrm{(}\mathrm{1}\mathrm{-}\mathit{\phi }\mathrm{(}\mathit{r}\mathrm{\left)}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}\mathit{B}{\mathit{b}}^{-3}{\mathit{\sigma }}_{\mathrm{1}}^{\frac{\mathrm{3}}{\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{E}}{\mathrm{ℛ}\mathit{T}}}\mathit{.}$(9)Apart from different numbers of contact points, all four packings have different initial porosities φ₀ (or initial densities D₀) and critical porosities φ_c (or densities D_c) at which the geometric model is no longer valid, and different parameters to be used in Eq. (6). These parameters and the parameters that will be described in Appendix A are summarised in Table 1.

2.2. Structural evolution of a planetesimal

We approximate the continuous accretion of a body by considering instantaneous formation at a certain time t₀ relative to the formation of the calcium-aluminium-rich inclusions (CAIs), at which accretion is completed. The planetesimals accrete as porous aggregates and undergo compaction via hot pressing due to increasing temperature in the presence of stress. Assuming a terminal radius $\overline{\mathit{R}}$, which corresponds to a body with the porosity equal to zero, the initial radius R_p(t₀) at the accretion onset time t₀ is calculated using the prescribed uniform initial porosity φ₀. In the case of instantaneous formation the initial radius is larger than the terminal radius and is given by ${\mathit{R}}_{\mathrm{p}}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\frac{\overline{\mathit{R}}}{{\left(\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{0}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}}\mathrm{·}$(10)The radius changes with time as a result of sintering converging against the terminal radius $\overline{\mathit{R}}$ for falling average porosity $\overline{\mathit{\phi }}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ at the time t: ${\mathit{R}}_{\mathrm{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{0}}}{\mathrm{1}\mathrm{-}\overline{\mathit{\phi }}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}{\mathit{R}}_{\mathrm{p}}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathit{.}$(11)However, the terminal radius can only be reached if the body becomes completely consolidated, which is usually not the case. Thus, a planetesimal experiences a macroscopic change (shrinkage) during the microscopic change (compaction) and the depth at which a particle is located changes with the radius and the porosity distribution.

Fig. 2 Porosity after the cooling of the body below 300 K as a function of depth. Left panel: $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$, models Bcreep (solid lines) and Bpar (dashed lines). Right panel: $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{5}\mathrm{km}\mathit{,}\mathrm{2.3}\mathrm{Ma}\mathrm{\right)}$, models Screep (solid lines) and Spar (dashed lines). For each model four different packings are considered. From top to bottom (at the depth of zero km) the lines (both solid and dashed) correspond to SCP, ORP, BCCP, and RHP. We note that for the single packings the different porosities in the primordial layer result in a variation of the radius (see Table 3). In addition, for a chosen packing with different final average porosities (i.e. different thickness of the porous layer), the associated solid and dashed lines have different lengths. This is not obvious here because the dashed lines are overlaid by the solid lines at depth.

3. Results

Compaction of a planetesimal via hot pressing is a direct result of heating by the radioactive decay of the short-lived isotopes ²⁶Al and ⁶⁰Fe. The initial concentrations of these isotopes in a planetesimal depend on the formation time relative to the formation of the CAIs. Because of the rapid decline in the concentration of these heat sources, a late formation time reduces the amount of heating significantly. The conductive heat loss through the surface counteracts the internal heating by radiogenic heat sources. The efficiency of this cooling depends on the surface-to-volume ratio of a body, i.e. the smaller the body, the more efficient the cooling. Another factor that influences the thermal evolution of a planetesimal is the presence of porous material which has a lower thermal conductivity compared to consolidated material and reduces conductive cooling. Assuming equal composition and material properties, a rough measure of the compaction potential of a planetesimal is the combination of its size (terminal radius $\overline{\mathit{R}}$; see Eq. (11)) and formation time (t₀, relative to the formation of the CAIs). In the following we will consider two bodies for which we choose $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{5}\mathrm{km}\mathit{,}\mathrm{2.3}\mathrm{Ma}\mathrm{\right)}$ (case S refers to small body) and $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$ (case B refers to big body). For both cases, we compare the extent of compaction using either the creep-related approach described in Sect. 2 (models S_creep and B_creep; see Table 3), or the parametrised approach, assuming linear decrease of the porosity in the temperature range 650–700 K (models S_par and B_par; see Table 3). The radii were chosen in order to simulate a relatively small and a larger planetesimal, whereas the corresponding formation times were chosen such that the bodies would neither compact completely, nor would they not compact at all. The assumed composition and the physical parameters correspond to an H chondrite (see Table 2).

For each case, i.e. S_creep, S_par, B_creep, and B_par, we investigate the thermal and structural evolution assuming four different packings of dust grains with the corresponding initial porosity and radius, as well as the corresponding computation of the effective stress (for the creep related approach, see also Appendix A). All other parameters are fixed and if not stated otherwise are identical to the parameters used in Neumann et al. (2012).

We further vary the activation energy E and the grain size b for case B_creep; both parameters strongly influence the compaction process. The activation energy is varied between 235 kJ mol^-1 and 477 kJ mol^-1 (see Sect. 3.2) corresponding to the measured activation energy for olivine (Schwenn & Goetze 1978). Olivine is assumed to dominate the rheology in ordinary chondrites. As a reference value we use E = 356 kJ mol^-1 (Schwenn & Goetze 1978). The grain size of the primordial material before compaction is an unknown parameter. It has been constrained by determining the size of the matrix grains in chondritic meteorites. These grains tend to be smaller than 5μm (Hutchison 2004). However, these grains may result from coarsening of smaller particles. Another constraint is provided by the interplanetary dust particles that may serve as a substitute of protoplanetary nebular particles. According to Rietmeijer (1993), the radii of interplanetary dust particles scatter from to , having a typical size of 0.1−0.2μm. We use as a reference value a grain radius of b = 1μm in the present study, but also investigate the influence of the grain size on compaction by varying b between 10^-5 m and 10^-7 m (see Sect. 3.3). The possibility of a non-trivial particle radius distribution that would act to reduce porosity by filling the gaps between larger grains with smaller ones is neglected here, as well as the possibility of irregular shapes, which would act to reduce the contact areas thereby increasing the porosity.

If not otherwise stated, the parameters b = 10^-6 m and E = 356 kJ mol^-1 are used in the models.

3.1. Parametrised vs. creep related compaction

Figure 2 shows final porosity profiles obtained for model B_creep and B_par (left panel), and S_creep and S_par (right panel). Four initial structures (packings) are considered, for which the initial (uniform) porosities vary between 26% and 47% (see Table 1).

The planetesimal starts as a porous body that heats up as a result of the decay of ²⁶Al and ⁶⁰Fe. The temperature profile in such a body with a homogeneous structure and heat source distribution is monotonous with a maximum at the centre. Thus, porosity reduction starts in the central region and the compaction front proceeds from the inside outwards. Thereby, the radius and the bulk volume decrease, while the bulk density increases. Eventually, a stable porosity profile is established (see Fig. 2). The porosity as a function of radius varies nearly exponentially between 0 in the completely compacted region and the initial porosity φ₀ in an outer layer consisting of primordial material. This insulating blanket of a variable thickness is retained at the surface, where, because of low temperature and pressure, compaction does not occur. While for decreasing porosity the heat source density increases and promotes a stronger heating, the decreasing thickness of the insulating porous layer and the shrinking radius enhance cooling of the body. The interplay of these three factors determines the temperature evolution.

For the parametrised approach of compaction (cases S_par and B_par), both bodies consolidate completely in the centre and have a primordial porous layer of a variable thickness at the surface. The change from the consolidated region to the porous layer is very abrupt. This behaviour can be explained by the temperature profile, which is very flat apart from the porous part where the temperature decreases rapidly toward the surface. The specific shape of the profile is even more pronounced due to the insulating effect of the outer low-conductivity region. As a consequence, when the central region exceeds the temperature of 700 K, the entire body except the outer ≲10 km (for B_par) or ≲4 km (for S_par) experiences similar heating and compacts to φ = 0, the porous blanket D_p being rather thin (5.5−9.5 km for B_par and 1.1−3.4 km for S_par; see Table 3). Using a different interval for sintering, e.g. 670−700 K, the same scenario would apply in general. The porous layer will be slightly thicker and the porosity profile will have a steeper slope between the compacted and porous regions. The lack of strong insulation and the decline of the heat production inhibit further temperature increase. Thus, for S_par and B_par, T_max varies between 706−734 K and 734−736 K, respectively, and remains below 800 K in all cases.

In contrast, for the creep related approach, the porosity changes nearly exponentially with time and temperature. Compaction of the planetesimal is less efficient because a temperature of 700 K is not sufficient to achieve φ = 0 even for the larger body with a reference radius of $\overline{\mathit{R}}\mathrm{=}\mathrm{50}$ km. Thus, the low conductivity of the porous material provokes further heating of the interior, until temperatures of up to 925 K (case S_creep) or only up to 739 K (case B_creep) are reached. Even for those temperatures, the porosity reduction is either moderate (case B_creep), or nearly non-existent (case S_creep, the profile corresponding to the packing RHP in Fig. 2). The outer porous layer remains significantly thicker than for the parametrised approach with values between 20−25 km for B_creep and more than 2 km (or even comprising the entire body) for S_creep. As a consequence, using the creep related approach for hot pressing results in much higher average porosities after compaction in contrast to the parametrised approach. For case S_creep, the average porosity is a factor of 1.2−1.9 larger than for S_par, and for case B_creep it is a factor of 1.8−2.8 larger than for B_par (we note that different factors correspond to different packings).

It should also be noted that the choice of the packing has a larger influence on small bodies with low pressure than on larger bodies with higher pressure. This can be seen in the final structure (Fig. 2 and Table 3) and therefore also has a consequence on the temperature evolution. The maximum temperatures vary by less than 40 K for different packings in the case of S_par, B_creep, and B_par, but for S_creep variation of the maximum temperature reaches 120 K. This finding is the result of the low pressure in small bodies, resulting in a strong variation of the porosity profiles for different packings.

Table 3 lists initial and final properties for all cases. The parameter I_T defines the temperature interval of compaction in the centre of the body under consideration. For creep-related models, the lower value of I_T corresponds to the temperature at which the porosity changes by 1% of its initial value φ₀, while the upper value corresponds to the temperature at which the porosity decreases below 10^-5. For parametrised models, the lower and upper values of I_T indicate the temperatures at which compaction starts and is completed, respectively (by definition, 650–700 K). If only one value is given, then the centre did not compact completely. In general, for the large body case (i.e. B_creep with an initial central pressure of 2.1−3.3 MPa depending on the value of φ₀), I_T varies around 650−740 K. For the small body case (i.e. S_creep with an initial central pressure of 0.021−0.033 MPa), the interval I_T varies around 740−925 K and for the packings RHP and BCCP, compaction in the centre stops even before reaching φ = 0. In the latter cases, the maximum temperature of 808 K or 869 K were not sufficiently high to reach φ = 0. It should be noted that I_T has been constrained for the centre of the respective body and the differences in I_T between the parametrised and the creep related cases are even larger in the more shallow regions.

To some extent, the timescale of the compaction can be estimated by using the relation ${\mathit{t}}_{\mathrm{comp}}\mathrm{\approx }\mathit{\eta }{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}$, where η = B^-1b³e ^E/RT (i.e. by inverting Eq. (9)). In the following we demonstrate this in two examples: B_creep with ORP (see Fig. 6, left panel) and S_creep with RHP (see Fig. 3).

Fig. 3 Central temperature (solid line) and central porosity (dashed line) as functions of time for the case Screep with RHP.

B_creep with ORP: the central porosity changes by 1% of its initial value of φ₀ = 0.3954 at the temperature of 655 K (see Table 3). The central pressure σ₀ at this time is 2.5 MPa. At the relative density of ≈0.61, the effective stress is σ₁ ≈ 200σ₀ ≈ 500 MPa. For this effective stress and for the temperature of T = 655 K, we obtain t_comp ≈ 52Ma. Subsequently, the temperature increases, while the effective stress and the porosity decrease. At t ≈ 6.6 Ma and T ≈ 722 K the values of φ ≈ 0.3, σ₁ ≈ 10σ₀ and σ₀ ≈ 2.9 MPa follow. Thus, t_comp is reduced to 8.6 Ma. In fact, the pores around the centre close at t = 12.8 Ma ≤6.6Ma + 8.6Ma = 15.2 Ma predicted by the last estimate (see Fig. 6; we note that the temperature increases further up to 739 K).

S_creep with RHP: the central porosity changes by 1% of its initial value of φ₀ = 0.2595 at the temperature of 740 K (see Table 3). The central pressure σ₀ at this time is 0.033 MPa. At the relative density of ≈0.7431, the effective stress is σ₁ ≈ 100σ₀ ≈ 3.3 MPa. For this effective stress and for the temperature of T = 740 K, we obtain t_comp ≈ 53 Ma. Subsequently, the temperature increases, while the effective stress and the porosity decrease. At t ≈ 3.34 Ma and T ≈ 808 K the values of φ ≈ 0.22, σ₁ ≈ 13σ₀, and σ₀ ≈ 0.0334 MPa follow. Thus, t_comp is reduced to 8.6 Ma. However, the temperature falls right after reaching 808 K rapidly (see Fig. 3 below) and compaction stalls.

3.2. Variation of the activation energy

In the following, we vary the creep activation energy for model B_creep within the interval E = 356 ± 121 kJ mol^-1 (Schwenn & Goetze 1978). The general compaction scenario is identical to the one described above. Figure 4 shows the final porosity profiles and Table 4 lists the final properties of the planetesimal obtained for different values of the activation energy. The final average porosity is ≈3% for E = 235 kJ mol^-1 and increases monotonously with increasing value of E. For the maximum value of E = 477 kJ mol^-1 compaction is almost negligible and the final porosity remains ≈φ₀. Although the interior structure differs significantly due to a variation in E, the maximum temperature in the centre increases only moderately with increasing E by ≈2% from 727 K to 740 K. This is caused by the decline in the heat production with time. Even if the interior remains more porous with increasing E and could heat up efficiently because of its low thermal conductivity, not enough energy is available to strongly heat the interior. Smaller activation energy implies compaction onset at a lower temperature and a shift of I_T to lower temperatures. Consequently, a larger part of the interior consolidates completely, and the thickness of the porous layer is smaller than for higher values of E, which also implies a weaker isolation and stronger cooling leading to a slightly lower T_max. The interval I_T experiences a considerable shift. While for E = 235 kJ mol^-1 compaction in the centre takes place at T ≥ 440 K, for E = 356 kJ mol^-1 it takes place at T ≥ 650.

Fig. 4 Porosity as a function of depth for $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$ for model Bcreep assuming ORP. The lines correspond to the adopted values of the activation energy E within the interval inferred by Schwenn & Goetze (1978).

For the extreme case of E = 477 kJ mol^-1 the threshold temperature for compaction has not been reached and thus the porosity change is infinitesimal.

Fig. 5 Porosity as a function of depth for $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$ for model Bcreep assuming ORP. The lines correspond to the adopted values of the grain size b. The grain size varies by two orders of magnitude between b = 10-7 and b = 10-5 m.

Fig. 6 The central temperature (solid lines) and the central porosity (dashed lines) as functions of time. Left panel: case $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$ with ORP and the associated effective stress σ1. Right panel: same as left panel but lithostatic pressure σ0 used instead of the effective stress σ1.

3.3. Variation of the grain size

In this section we address the influence of the grain size while keeping all other parameters fixed. Figure 5 shows the final porosity profiles. The other output parameters such as average final porosity ${\overline{\mathit{\phi }}}_{\mathrm{f}}$, final planetesimal radius R_f, thickness of the porous layer D_p, central maximum temperature T_max, and the temperature interval I_T in which compaction took place in the centre, are listed in Table 5. All parameters are identical to those used for the case B_creep with ORP, except the grain size b. There is a straightforward relationship between the adopted value of the grain radius and the extent of compaction. For decreasing b the final average porosity and the maximum temperature decrease as well. Assuming b = 10^-5 m (i.e. b = 10μm), the porosity starts decreasing at a temperature of 716 K, but as the maximum temperature is only 740 K, compaction stops at φ = 0.33 in the centre. For the reference case with b = 10^-6 m, a substantial part of the interior becomes completely consolidated. With decreasing grain size, i.e. b = 10^-7 m for the case B_creep, compaction is even more effective, taking place between 597 K and 694 K in the centre. Although the temperature interval I_T for b = 10^-6 m is closer to 650−700 K, at smaller depths the interval changes towards higher temperatures. Eventually, a porous outer layer with a thickness of only ≈5 km is established for b = 10^-7 m. Similarly to the effect of the activation energy, a smaller value of b favours creep at lower temperatures. This leads to a thinner insulating layer and better cooling, resulting in lower temperatures and thus lower T_max.

In conclusion, it is easier to achieve consolidation of a fine-grained material via hot pressing than for coarse-grained material. For the latter, higher temperatures than those that are possible for the pair $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$ are needed to change the porosity further (e.g. through an earlier formation time). We have explored this aspect with some test runs by increasing the grain radius even further up to the size of 1 cm appropriate for chondrules (Weisberg et al. 2006). For that extreme case, some residual porosity was present even at super-solidus temperatures of the metal phase and even of the silicate phase (T ≥ 1425 K) of ordinary chondritic material. This means that at moderate pressures and for certain values of activation energy, compaction takes place simultaneously with the production of partial melt (see Discussion).

3.4. Effective stress

The compaction model used considers that for porous material the effective stress can be several orders of magnitude higher than the stress caused simply by lithostatic pressure (see Fig. 2). In Fig. 6 we show the influence of the efffective stress on the evolution of the central temperature (solid lines) and of the central porosity (dashed lines). The left panel corresponds to the case $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$ with ORP and the associated effective stress σ₁, whereas the right panel corresponds to a simulation where the lithostatic pressure was used instead of the effective stress. While in the first case the central porosity changes by 1% of its initial value φ₀ at the temperature of 655 K and reaches zero at 739 K (see Table 3), in the second case compaction requires a higher threshold temperature and starts at 732 K. Compaction is, however, completed for 739 K as in the case using the effective stresses. Additional output parameters of this simulation are ${\overline{\mathit{\phi }}}_{\mathrm{f}}\mathrm{=}\mathrm{34.40}$%, R_f = 57.57 km, D_p = 41.20 km, and T_max = 739 K (for the same case with a proper σ₁ we had ${\overline{\mathit{\phi }}}_{\mathrm{f}}\mathrm{=}\mathrm{21.15}$%, R_f = 54.12 km, D_p = 23.70 km, and T_max = 739 K; see Table 3). Although the maximum temperature is the same, the smaller values of σ₀ compared to those of σ₁ (see Fig. 1) lead to a less effective compaction. Thus, the insulation of the centre is more effective (D_p = 41.20 km) and the temperature of ≥700 K in the centre can be maintained for more than 60 Ma. This allows the centre to compact to φ = 0 after ≈40 Ma (in the standard case with σ₁ only ≈8 Ma were needed, compare the compaction time t_comp in Table 3).

The compaction timescale can be roughly estimated as in 3.1. If not the effective stress, but the central lithostatic pressure is used, for T = 739 K and σ₀ = 2.5 MPa we obtain t_comp ≈ 87 Ma. However, the pressure increases upon compaction up to the value of 4.6 MPa. For this value, t_comp ≈ 35 Ma follows. This agrees very well with the results obtained from the simulations.

The porosity and temperature in Fig. 6 show that compaction by hot pressing is influenced strongly by the interplay of the temperature and the applied stresses and their temporal evolution. A lack of pressure can be compensated by heating far above some critical temperature (onset of compaction at 732 K). Furthermore, although for t₀ = 3.2 Ma a higher maximum temperature than 739 K cannot be reached, a drastically prolonged heating of the centre due to an effective insulation by the porous shell is possible, leading eventually to the closure of the pores.

4. Summary and discussion

Here we compare hot pressing via the creep related approach (Kakar & Chaklader 1967; Rao & Chaklader 1972; Schwenn & Goetze 1978) with the simplified parametrised approach used by several studies (e.g. Hevey & Sanders 2006; Sahijpal et al. 2007; Gupta & Sahijpal 2010; Solano et al. 2012). Using the creep related approach, hot pressing is triggered when a powder assemblage reaches a certain temperature. This critical (or threshold) temperature is a function of the effective stress for otherwise fixed parameters. The higher the effective stress, the lower the critical temperature. Thus, the threshold temperature varies with depth (pressure) within a planetesimal. Since compaction is a creep process, to fully compact the material, this threshold temperature must either be sustained for a sufficiently long time, or higher temperatures are necessary to enhance the compaction rate. Thus, compaction in a thermally evolving planetesimal takes place in a certain temperature interval for a given stress. In our study, we found that for a central pressure of 200−300 MPa (equivalent to a planetesimal with a radius of 500 km) the temperature interval is 597−690 K, for a pressure of 2−3 MPa (equivalent to a planetesimal with a radius of 50 km) it is 645−739 K, and for a pressure of 0.02 − 0.03 MPa (equivalent to a planetesimal with a radius of 5 km) it is 740−925 K, assuming an activation energy of 356 kJ mol^-1 (Table 2).

In conclusion, only for a rather narrow parameter space or pressure range in a planetesimal, compaction occurs approximately between 650 K and 700 K as used in the parametrised approach. Depending on the pressure, initial grain size, activation energy and initial porosity, the temperature interval is instead 500 − 1000 K. For some extreme values of E and b, the upper boundary exceeds even the solidus temperature of the silicates. Significant differences in the interior structure and temperature evolution are the consequence; with the simplified approach compaction is overestimated predicting final porosities that are too low and thus a cooling that is unrealistically fast.

For the small body of 5 km radius and a formation time t₀ of 2.3 Ma the porous layer is at least 1 km thinner and for a body of 50 km radius and t₀ = 3.2 Ma at least 12 km, if the parametrised approach is used. Differences in the maximum temperatures between the two approaches can be observed, especially for small bodies (radius ≈). The maximum temperature varies by up to 32% as the consequence of different extents of compaction.

We note that although porosity plays a substantial role in the thermal evolution and interior structure of undifferentiated planetesimals, in the case of high levels of partial melt and efficient differentiation, effects associated with porosity and compaction would be painted over in the subsequent evolution.

The parameter that has the strogest influence is the activation energy E of the material. A small value of E shifts the onset of compaction to lower temperatures, and a high E to higher temperatures. For our cases with E = 477 kJ mol^-1, some left-over porosity is observed in the centre directly before the melting of the metal-rich phase begins. This in fact requires a simultaneous treatment of porosity loss and melting. Experiments by e.g. Kohlstedt (2002) and Jin et al. (1994) suggest that a few percent of melt can reduce the creep strength by more than one order of magnitude because of the fast diffusion of the melt (see also Karato 2013). Therefore, it is likely that in these specific cases partial melting provokes quasi-instantaneous compaction on a geological timescale.

The activation energy is in general related to the mineralogical composition. For planetesimals with a “wet” rheology, i.e. containing hydrated silicates or ice, as suggested for carbonaceous chondritic compositions that are dominated by serpentine, the activation energy is lower than for dry material and thus lower compared to the ordinary chondritic composition assumed in the current study (Hilairet et al. 2007). For the ordinary chondritic composition it is assumed that olivine dominates the rheology. On the other hand, the creep activation energy of an enstatite chondritic composition seems to be higher than that of the ordinary chondritic composition (Karato 2013). This would also suggest that a carbonaceous chondritic body will compact at lower temperatures than the ordinary chondritic and the enstatitic ones at similar pressure. It should be noted, however, that the pre-factor B, the stress exponent n, and the grain size exponent m need to be adjusted to the rheological properties of the material.

The shift of the compaction temperature interval to lower temperatures in a body consisting of a carbonaceous chondritic assemblage also implies that compaction and dehydration of hydrated minerals (≈298–1073 K, Castillo-Rogez 2011) occur simultaneously. The shrinkage due to hot pressing overlaps with the shrinkage resulting from the density changes of the silicate minerals due to the dehydration on the one hand and by the filling of the pores with the water freed by the dehydration process on the other hand. Furthermore, dehydration of certain minerals proceeds in some cases simultaneously with the hydration of other silicate minerals, which could be accompanied by expansion. Thus, compaction of planetesimals that contained water in some form is an even more complex issue than compaction of rocky bodies.

Four different packings have been compared on the compaction process. The assumption on the packing of the dust grains influences the initial porosity and is particularly important for small bodies with $\overline{\mathit{R}}\mathrm{=}\mathrm{𝒪}\mathrm{\left(}\mathrm{1}\mathrm{km}\mathrm{\right)}$; a variation in T_max by up to 16% has been observed. With increasing body radii the offset becomes gradually smaller and for a planetesimal radius of 50 km and larger the choice of the suitable packing of grains is less important.

We have further examined the influence of the effective stress on the compaction history. The results indicate that neglecting the effective stress, which acts on the contact points of the grains in case of porous material, and assuming only the lithostatic pressure leads to a longer timescale of compaction and to a higher threshold temperature at which compaction starts. Thus, compaction is underestimated when neglecting the effective stress. This finding might have implications for a number of studies that calculate or estimate the porosity reduction of a porous layer and only consider the lithostatic pressure for compaction (e.g. Wieczorek et al. 2013; Nimmo et al. 2003; Fu & Elkins-Tanton 2014). The timescale of the porosity loss t_comp varies in our calculations between and (see Table 3), which argues against assumptions of fixed compaction timescales utilised in some planetesimal studies.

In the current approach, we neglect compaction prior to hot pressing and assume that the initial porosity φ₀ of the material is limited by the porosities of the random close and random loose packing. The initial porosity, however, may be reduced by a purely mechanical mechanism called cold pressing, which is driven by the applied pressure. Cold pressing does not impose any further requirements and can act at any temperature. Porosity evolution due to cold pressing is an important topic in the study of dust coagulation and planetesimal formation in protoplanetary discs and was studied in theory in e.g. Ormel et al. (2007), Suyama et al. (2008), Zsom et al. (2010, 2011), Kataoka et al. (2013a). This kind of compaction was also discussed and investigated through the laboratory experiments of Guettler et al. (2009) and the numerical experiments of Kataoka et al. (2013b). The latter authors derived a relation for the change of the relative density D with the applied stress σ₀. In the limit σ₀ → ∞, the maximum relative density of D_max = 0.58 (or minimal porosity ${\mathit{\phi }}_{\min }^{\mathrm{cold}}\mathrm{=}\mathrm{0.42}$ is obtained. In conclusion, it is only possible to compact porous material by cold pressing to the limit of φ = 0.42. Since SCP with φ₀ = 0.4764 is the only model with an initial porosity of φ₀ ≥ 0.4, we neglect this kind of compaction in the present study.

As a further complication, however, recent dynamic cold pressing experiments of Beitz et al. (2013) indicate that this process might have been important during accretion. Beitz et al. (2013) used aluminium cylinders as projectiles to compact the chondrite-analog samples with the initial porosities of 39 − 70% in a velocity range of 1200−165 m s^-1, down to the porosities of 1−30%. Taking the projectile velocity as velocity of the accreting dust particles suggests that a smaller initial porosity φ₀ should be used than the values corresponding to the packings. Moreover, this cold pressing experiment does most certainly result in crushing of the dust grains, meaning that the geometric model for the computation of the effective stress would be not applicable without further adjustments. Here future studies are needed in order to assess how this might change compaction and the structural and thermal evolution of a planetesimal.

Online material

Appendix A: Geometric model for the computation of the effective stress

The geometric approach to the deformation theory of hot pressing and the conclusions for the plastic flow have been investigated by Kakar & Chaklader (1967) and Rao & Chaklader (1972). However, the detailed derivation of the formulae was not provided. Thus, we present the theory in detail and improve some aspects that were treated in a simplified manner in the studies mentioned above.

In the study by Kakar & Chaklader (1967) four basic systems of packing are considered, defined by the geometric arrangement of the grains and by the coordination number (this is the number Z of the contact points, an average value for equally sized spheres). In the simple cubic packing (SCP) every grain has four neighbours in a horizontal layer and two additional neighbours in the layers above and below. Orthorhombic (ORP) packing corresponds to six neighbours in a layer (plus two above and below). The rhombohedral packing (RHP) has coordination number 12 and the body-centred cubic packing (BCCP) 8. Although in a unidirectional force field the body-centred cubic packing is an unstable arrangement, we consider this case for the sake of completeness. Furthermore, although stability varies with the packing type, we assume that the symmetry is sustained throughout the whole compaction process. Assuming homogeneous and isotropic composition of the grains and symmetrical spread upon application of stress, for different packings the spheres will deform to polyhedra upon complete compaction, e.g. cubes for SCP, hexagonal prisms for ORP, rhombic dodecahedra for RHP, and tetrakaidecahedra for BCCP. During the deformation, necks develop at the contact points. The form of the grains will change, whereas the volume will remain constant.

For each kind of packing a specific unit cell of a certain form which is occupied by a dust particle is considered. The form of a unit cell is defined by the polyhedra mentioned above. From Figs. A.2–A.4, the volume V of such a unit cell can be computed. For SCP, this is the volume of a cube, and, since $\mathit{y}\mathrm{=}\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}}$, we obtain $\mathit{V}\mathrm{=}{\left(\mathrm{2}\mathit{y}\right)}^{\mathrm{3}}\mathrm{=}\mathrm{8}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}\mathit{.}$(A.1)For ORP, the volume of a hexagonal prism is $\mathit{V}\mathrm{=}\frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{2}}{\left(\frac{\mathrm{4}}{\sqrt{\mathrm{3}}}\mathit{y}\right)}^{\mathrm{2}}\left(\mathrm{2}\mathit{y}\right)\mathrm{=}\mathrm{4}\sqrt{\mathrm{3}}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}\mathit{.}$(A.2)For RHP we compute $\mathit{V}\mathrm{=}\frac{{\left(\mathrm{2}\mathit{y}\right)}^{\mathrm{3}}}{\sqrt{\mathrm{2}}}\mathrm{=}\mathrm{4}\sqrt{\mathrm{2}}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}\mathit{,}$(A.3)and for BCCP $\mathit{V}\mathrm{=}\frac{\mathrm{8}\sqrt{\mathrm{2}}}{\mathrm{27}}{\left(\sqrt{\mathrm{6}}{\mathit{y}}_{\mathrm{1}}\right)}^{\mathrm{3}}\mathrm{=}\frac{\mathrm{32}\sqrt{\mathrm{3}}}{\mathrm{9}}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{{\mathit{a}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}}^{\frac{\mathrm{3}}{\mathrm{2}}}\mathit{.}$(A.4)Thus, the volume of the unit cell can be expressed as the product of a constant β₁ (the pre-factor of the bracket) and the power ${{}^{\mathrm{(}}\mathit{r}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\mathrm{)}}^{}}^{\mathrm{3}\mathit{/}\mathrm{2}}$.

Fig. A.1 General deformation geometry at any contact point between two grains.

The deformation of the BCC packing proceeds in two stages. For this reason we describe the theory for the other cases and the first stage of the BCCP deformation first, and then consider the second deformation stage of the BCC packing separately.

Prior to the hot pressing, the grain touches the surface of the unit cell at Z points where the neighbouring grains touch each other (Z₁ for BCCP). During compaction, the grains deform plastically at the contact points and develop necks. The cross-section of such a neck between two grains has a circular form and is called the face (the contact points have evolved to contact cross-sections). Still, the volume of the particle needs to be kept constant. Thus, the current radius of the grain is assumed to increase and the volume of those parts of the grain now lying outside of the unit cell (the caps) need to be subtracted from the current sphere volume (we note that the initial sphere volume is normalised to 1): $\mathrm{1}\mathrm{=}\mathrm{const}\mathrm{.}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}^{\mathrm{3}}\mathrm{-}\mathit{Z}{\mathit{V}}_{\mathrm{cap}}\mathit{,}$(A.5)where r is the current grain radius and ${\mathit{V}}_{\mathrm{cap}}\mathrm{=}\mathit{\pi }\frac{{\mathit{h}}^{\mathrm{2}}}{\mathrm{3}}\mathrm{(}\mathrm{3}\mathit{r}\mathrm{-}{\mathit{h}}^{\mathrm{)}}$(A.6)is the volume of a cap. Using the relation h = r − y defined by the geometry of the packings (Fig. A.1), we obtain the normalised volume of the particle for the first three packings as $\begin{array}{ccc}\mathrm{1}& \mathrm{=}& \frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}^{\mathrm{3}}\mathrm{-}\mathit{Z}{\mathit{V}}_{\mathrm{cap}}\\ & \mathrm{=}& \frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}^{\mathrm{3}}\mathrm{-}\frac{\mathit{Z\pi }}{\mathrm{3}}{\left(\mathit{r}\mathrm{-}\mathit{y}\right)}^{\mathrm{2}}\left(\mathrm{2}\mathit{r}\mathrm{+}\mathit{y}\right)\\ & \mathrm{=}& \frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}^{\mathrm{3}}\mathrm{-}\frac{\mathit{Z\pi }}{\mathrm{3}}\left(\mathrm{2}{\mathit{r}}^{\mathrm{3}}\mathrm{-}\mathrm{2}{\mathit{r}}^{\mathrm{2}}\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}}\right)\\ & \mathrm{=}& \frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}^{\mathrm{3}}\mathrm{-}\frac{\mathit{Z\pi }}{\mathrm{3}}\mathrm{2}{\mathit{r}}^{\mathrm{3}}\mathrm{+}\frac{\mathit{Z\pi }}{\mathrm{3}}\mathrm{(}\mathrm{2}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{)}\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}}\\ & \mathrm{=}& \frac{\mathit{Z\pi }}{\mathrm{3}}\mathrm{(}\mathrm{2}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{)}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{-}\frac{\mathrm{(}\mathrm{2}\mathit{Z}\mathrm{-}{\mathrm{4}}^{\mathrm{)}}\mathit{\pi }}{\mathrm{3}}{\mathit{r}}^{\mathrm{3}}\\ & \mathrm{=}& \mathit{\alpha }\mathrm{(}\mathrm{2}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{)}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{-}\mathit{\beta }{\mathit{r}}^{\mathrm{3}}\\ & \mathrm{=}& \end{array}$(A.7)This defines the geometric constants α and β from Table 1 as α = Zπ/ 3 and β = (2Z − 4)π/ 3. From Eq. (A.7) the current radius of the particle is computed:

Fig. A.2 Unit cells and their intersections with an inscribed grain for the SCP and the ORP (top row from left to right) along with the top views of the unit cells and of the grain at its initial and critical (final) stages of deformation (bottom row).

$\mathit{r}\mathrm{=}{\left[\mathit{\alpha }\left(\mathrm{2}\mathrm{+}{\left(\frac{\mathit{a}}{\mathit{r}}\right)}^{\mathrm{2}}\right){\left(\mathrm{1}\mathrm{-}{\left(\frac{\mathit{a}}{\mathit{r}}\right)}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{-}\mathit{\beta }\right]}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}}\mathrm{·}$(A.8)Given the ratio a/r, the current radius r of the sphere during deformation can be found (in arbitrary units corresponding to unit volume of the sphere). The ratio a/r is computed from the expression for the bulk density D, defined as the ratio of the volume of the particle to the volume of the unit cell. This last value is β₁(r² − a²)^3/2 (see Table 1). From the normalised grain volume 1 and the volume of the unit cell we obtain $\mathit{D}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\beta }}_{\mathrm{1}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}\mathrm{·}$(A.9)Prior to compaction, i.e. for a = 0, the radius equals the initial radius (r = r₀), hence $\mathit{D}\mathrm{=}{\mathit{D}}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\beta }}_{\mathrm{1}}{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}$(A.10)holds, thus, $\begin{array}{ccc}\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{0}}}{\mathrm{1}\mathrm{-}\mathit{\phi }}& \mathrm{=}& \frac{{\mathit{D}}_{\mathrm{0}}}{\mathit{D}}\mathrm{=}{\left[\frac{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\right]}^{\frac{\mathrm{3}}{\mathrm{2}}}\mathrm{=}{\left[{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{2}}\mathrm{-}{\left(\frac{\mathit{a}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{2}}\right]}^{\frac{\mathrm{3}}{\mathrm{2}}}\\ & \mathrm{=}& \end{array}$(A.11)This is solved for (a/r)²: ${\left(\frac{\mathit{a}}{\mathit{r}}\right)}^{\mathrm{2}}\mathrm{=}\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{r}}_{\mathrm{0}}}{\mathit{r}}\right)}^{\mathrm{2}}{\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{0}}}{\mathrm{1}\mathrm{-}\mathit{\phi }}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}\mathrm{·}$(A.12)For the computation of the effective stress, we also need the area of the cross-section of the unit cell. In terms of the geometries shown in Figs. A.2–A.4, the cross-section of a unit cell is the area of its projection onto the horizontal plane. For SCP this is a square with the edge length 2y, for ORP a regular hexagon with the edge length $\mathrm{(}\mathrm{2}\mathit{/}\sqrt{\mathrm{3}}\mathrm{)}\mathit{y}$, for RHP a square with the edge length 2y, and for BCCP an irregular octagon with four edges of the length $\sqrt{\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{y}}_{\mathrm{1}}$, four edges of the length $\mathrm{(}\mathrm{2}\mathit{/}\sqrt{\mathrm{3}}\mathrm{)}{\mathit{y}}_{\mathrm{1}}$ and the internal angle of 135 degrees. The areas of those polygons can be written as γy² (or as $\mathit{\gamma }{\mathit{y}}_{\mathrm{1}}^{\mathrm{2}}$ for BCCP; see Table 1 for the values of γ).

Fig. A.3 Unit cell of the RHP and its intersections with an inscribed grain (left panel), the top right view of the unit cell and of the grain at its initial and critical (final) stage of deformation (central panel), and a single face of the unit cell at the critical stage of the deformation (right panel).

Fig. A.4 Unit cells and their intersections with an inscribed grain for the two deformation stages of the BCCP (top row, first and second panels); the faces of the unit cell at the second critical deformation stage (top row, third and fourth panels); topview of the unit cell before the first critical stage (bottom row, left panel) and after the first but before the second critical stage (bottom row, central panel); and the arrangement of the spheres at the first critical deformation stage as seen in the plane defined by the centres of the three spheres (bottom row, right panel).

If a load m acts on a face then the stress on the face is equal m/ (γy²). The total stress on the sphere is Zm/ (γy²). In a planetesimal, the applied stress is the lithostatic pressure σ₀ and is equal to the stress on the face: ${\mathit{\sigma }}_{\mathrm{0}}\mathrm{=}\frac{\mathit{m}}{\mathit{\gamma }{\mathit{y}}^{\mathrm{2}}}\mathrm{·}$(A.13)From this the load acting on the cross-section of the unit cell follows as m = σ₀γy² and can be divided into two components. It is equal to the sum of the load acting on the sphere through N non-vertical flat faces πa²cos(θ)σ₁ lying above the equator of the unit cell and of the load acting on the void space cσ₂, $\mathit{\gamma }{\mathit{y}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{=}\mathit{N\pi }{\mathit{a}}^{\mathrm{2}}\cos \mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{\sigma }}_{\mathrm{1}}\mathrm{+}\mathit{c}{\mathit{\sigma }}_{\mathrm{2}}\mathit{,}$(A.14)where πa² is the area of the neck, θ is the inclination angle of the faces, c is the cross-section of the porous part of the unit cell, γ is a geometric constant which depends on the packing (see Table 1), σ₁ is the stress acting along the z-axis through the neck, and σ₂ = 0 is the stress on the void space. The last equation is solved for σ₁ to obtain ${\mathit{\sigma }}_{\mathrm{1}}\mathrm{=}\frac{\mathit{\gamma }}{\mathit{N\pi }\cos \mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{\left(\frac{\mathit{y}}{\mathit{a}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{=}\mathrm{:}\frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{1}}}\frac{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}}{\mathit{\sigma }}_{\mathrm{0}}\mathit{.}$(A.15)Above we have shown $\mathit{D}\mathrm{=}{\mathit{\beta }}_{\mathrm{1}}^{-1}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}$. From this we compute $\begin{array}{ccc}\mathit{D}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{\beta }}_{\mathrm{1}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}\\ & \mathrm{\iff }& {\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{D}}^{\frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\beta }}_{\mathrm{1}}^{\frac{\mathrm{2}}{\mathrm{3}}}}\\ & \mathrm{\iff }& \frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}}\mathrm{=}{\mathit{\beta }}_{\mathrm{1}}^{\frac{\mathrm{2}}{\mathrm{3}}}{\mathit{D}}^{\frac{\mathrm{2}}{\mathrm{3}}}\\ & \mathrm{\iff }& \frac{{\mathit{r}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}}\mathrm{=}{\mathit{r}}^{\mathrm{2}}{\mathit{\beta }}_{\mathrm{1}}^{\frac{\mathrm{2}}{\mathrm{3}}}{\mathit{D}}^{\frac{\mathrm{2}}{\mathrm{3}}}\\ & \mathrm{\iff }& \frac{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}}\mathrm{=}{\mathit{r}}^{\mathrm{2}}{\mathit{\beta }}_{\mathrm{1}}^{\frac{\mathrm{2}}{\mathrm{3}}}{\mathit{D}}^{\frac{\mathrm{2}}{\mathrm{3}}}\\ & \mathrm{\iff }& \mathrm{1}\mathrm{+}\frac{{\mathit{a}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}^{\mathrm{2}}}\mathrm{=}{\mathit{r}}^{\mathrm{2}}{\mathit{\beta }}_{\mathrm{1}}^{\frac{\mathrm{2}}{\mathrm{3}}}{\mathit{D}}^{\frac{\mathrm{2}}{\mathrm{3}}}\\ & \mathrm{\iff }& \end{array}$(A.16)Thus, ${\mathit{\sigma }}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{1}}}\frac{\mathrm{1}}{{\mathit{D}}^{\frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\beta }}_{\mathrm{1}}^{\frac{\mathrm{2}}{\mathrm{3}}}{\mathit{r}}^{\mathrm{2}}\mathrm{-}\mathrm{1}}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{·}$(A.17)To obtain the particle volume in the second deformation stage of the BCCP, one needs to subtract Z₁ = 8 caps of the volume ${\mathit{V}}_{\mathrm{cap,}\mathrm{1}}\mathrm{=}\mathit{\pi }\frac{{\mathit{h}}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{3}}\mathrm{(}\mathrm{3}\mathit{r}\mathrm{-}{\mathit{h}}_{\mathrm{1}}\mathrm{)}$ and additionally Z₂ = 6 caps of the volume ${\mathit{V}}_{\mathrm{cap,}\mathrm{2}}\mathrm{=}\mathit{\pi }\frac{{\mathit{h}}_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{3}}\mathrm{(}\mathrm{3}\mathit{r}\mathrm{-}{\mathit{h}}_{\mathrm{2}}\mathrm{)}$: $\begin{array}{ccc}\mathrm{1}& \mathrm{=}& \frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}^{\mathrm{3}}\mathrm{-}{\mathit{Z}}_{\mathrm{1}}{\mathit{V}}_{\mathrm{cap,}\mathrm{1}}\mathrm{-}{\mathit{Z}}_{\mathrm{2}}{\mathit{V}}_{\mathrm{cap,}\mathrm{2}}\\ & \mathrm{=}& \frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}^{\mathrm{3}}\mathrm{-}{\mathit{Z}}_{\mathrm{1}}\mathit{\pi }\frac{{\mathit{h}}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{3}}\mathrm{(}\mathrm{3}\mathit{r}\mathrm{-}{\mathit{h}}_{\mathrm{1}}\mathrm{)}\mathrm{-}{\mathit{Z}}_{\mathrm{2}}\mathit{\pi }\frac{{\mathit{h}}_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{3}}\mathrm{(}\mathrm{3}\mathit{r}\mathrm{-}{\mathit{h}}_{\mathrm{2}}\mathrm{)}\\ & \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}}\frac{{\mathit{Z}}_{\mathit{i}}\mathit{\pi }}{\mathrm{3}}\mathrm{(}\mathrm{2}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{{\mathit{a}}_{\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{{\mathit{a}}_{\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{-}\frac{\mathrm{(}\mathrm{2}{\mathit{Z}}_{\mathrm{1}}\mathrm{+}\mathrm{2}{\mathit{Z}}_{\mathrm{2}}\mathrm{-}{\mathrm{4}}^{\mathrm{)}}\mathit{\pi }}{\mathrm{3}}{\mathit{r}}^{\mathrm{3}}\\ & \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}}{\mathit{\alpha }}_{\mathrm{0}\mathit{,i}}\mathrm{(}\mathrm{2}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{{\mathit{a}}_{\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{{\mathit{a}}_{\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{-}\mathit{\beta }{\mathit{r}}^{\mathrm{3}}\\ & \mathrm{=}& \end{array}$(A.18)This defines the geometric constants ${\mathit{\alpha }}_{\mathrm{0}\mathit{,}\mathrm{1}}\mathrm{=}\frac{{\mathit{Z}}_{\mathrm{1}}\mathit{\pi }}{\mathrm{3}}$, ${\mathit{\alpha }}_{\mathrm{0}\mathit{,}\mathrm{2}}\mathrm{=}\frac{{\mathit{Z}}_{\mathrm{2}}\mathit{\pi }}{\mathrm{3}}$ and $\mathit{\beta }\mathrm{=}\frac{\left(\mathrm{2}{\mathit{Z}}_{\mathrm{1}}\mathrm{+}\mathrm{2}{\mathit{Z}}_{\mathrm{2}}\mathrm{-}\mathrm{4}\right)\mathit{\pi }}{\mathrm{3}}$ for BCCP. Similarly to the other cases, Eq. (A.18) implies $\mathit{r}\mathrm{=}{\left[\sum _{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}}{\mathit{\alpha }}_{\mathrm{0}\mathit{,i}}\left(\mathrm{2}\mathrm{+}{\left(\frac{{\mathit{a}}_{\mathit{i}}}{\mathit{r}}\right)}^{\mathrm{2}}\right){\left(\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{a}}_{\mathit{i}}}{\mathit{r}}\right)}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{-}\mathit{\beta }\right]}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}}\mathit{.}$(A.19)The power (a₁/r)² is computed analogously to Eq. (A.12): ${\left(\frac{{\mathit{a}}_{\mathrm{1}}}{\mathit{r}}\right)}^{\mathrm{2}}\mathrm{=}\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{r}}_{\mathrm{0}}}{\mathit{r}}\right)}^{\mathrm{2}}{\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{0}}}{\mathrm{1}\mathrm{-}\mathit{\phi }}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}\mathrm{·}$(A.20)For (a₂/r)² we consider the first critical stage c₁ when the cap of the second kind with the radius a₂ starts to form and express the relative density in the terms of y₂,

$\mathit{D}\mathrm{=}\frac{\mathrm{1}}{\frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{8}}{\mathit{\beta }}_{\mathrm{1}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}\mathit{,}$(A.21)because ${\mathit{y}}_{\mathrm{1}}\mathrm{=}\mathrm{(}\sqrt{\mathrm{3}}\mathit{/}\mathrm{2}\mathrm{)}{\mathit{y}}_{\mathrm{2}}$. Then $\begin{array}{ccc}\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{{\mathrm{c}}_{\mathrm{1}}}}{\mathrm{1}\mathrm{-}\mathit{\phi }}& \mathrm{=}& \frac{{\mathit{D}}_{{\mathrm{c}}_{\mathrm{1}}}}{\mathit{D}}\mathrm{=}\frac{\frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{8}}{\mathit{\beta }}_{\mathrm{1}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}{\frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{8}}{\mathit{\beta }}_{\mathrm{1}}{\mathit{r}}_{{\mathrm{c}}_{\mathrm{1}}}^{\mathrm{3}}}\mathrm{=}\frac{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{3}}{\mathrm{2}}}}{{\mathit{r}}_{{\mathrm{c}}_{\mathrm{1}}}^{\mathrm{3}}}\\ & \mathrm{=}& {\left(\frac{\mathit{r}}{{\mathit{r}}_{{\mathrm{c}}_{\mathrm{1}}}}\right)}^{\mathrm{3}}{\left[\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{a}}_{\mathrm{2}}}{\mathit{r}}\right)}^{\mathrm{2}}\right]}^{\frac{\mathrm{3}}{\mathrm{2}}}\end{array}$(A.22)and ${\left(\frac{{\mathit{a}}_{\mathrm{2}}}{\mathit{r}}\right)}^{\mathrm{2}}\mathrm{=}\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{r}}_{{\mathrm{c}}_{\mathrm{1}}}}{\mathit{r}}\right)}^{\mathrm{2}}{\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{{\mathrm{c}}_{\mathrm{1}}}}{\mathrm{1}\mathrm{-}\mathit{\phi }}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}$(A.23)follow. The relation between a₁ and a₂ is obtained from Eqs. (A.20) and (A.23) as ${\mathit{a}}_{\mathrm{2}}\mathrm{=}{\left[\frac{\mathrm{1}\mathrm{-}{\mathrm{(}{\frac{{\mathit{r}}_{{\mathrm{c}}_{\mathrm{1}}}}{\mathit{r}}}^{\mathrm{)}}}^{\mathrm{2}}{\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{{\mathrm{c}}_{\mathrm{1}}}}{\mathrm{1}\mathrm{-}\mathit{\phi }}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}}{\mathrm{1}\mathrm{-}{\mathrm{(}{\frac{{\mathit{r}}_{\mathrm{0}}}{\mathit{r}}}^{\mathrm{)}}}^{\mathrm{2}}{\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{0}}}{\mathrm{1}\mathrm{-}\mathit{\phi }}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{a}}_{\mathrm{1}}\mathrm{=}\mathit{c}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}{\mathit{a}}_{\mathrm{1}}\mathit{.}$(A.24)The load on the octagonal cross-section of the unit cell of the BCCP is $\mathit{\gamma }{\mathit{y}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathrm{1}}^{\mathrm{2}}\cos \mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{\sigma }}_{\mathrm{1}}\mathrm{+}\mathit{\pi }{\mathit{a}}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}\mathit{,}$(A.25)where σ₁ is the effective stress on the faces of the first kind, σ₂ on the second kind, and $\cos \mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$. Assuming equal load per face, we obtain $\mathit{\pi }{\mathit{a}}_{\mathrm{1}}^{\mathrm{2}}\cos \mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{\sigma }}_{\mathrm{1}}\mathrm{=}\mathit{\pi }{\mathit{a}}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}$ and from this $\mathit{\gamma }{\mathit{y}}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{=}\frac{\mathit{N}}{\cos \mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\mathit{\pi }{\mathit{a}}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{1}}\mathit{,}$(A.26)where N = 5. Thus, ${\mathit{\sigma }}_{\mathrm{1}}\mathrm{=}\frac{\mathit{\gamma }}{\mathit{N\pi }\cos \mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{\left(\frac{{\mathit{y}}_{\mathrm{1}}}{{\mathit{a}}_{\mathrm{1}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{1}}}\frac{{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{a}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{a}}_{\mathrm{1}}^{\mathrm{2}}}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{·}$(A.27)The geometric relations described above are valid from the beginning of the deformation until the faces start to touch each other. In the beginning, for the normalised grain volume 1 the initial grain radius r₀ is equal to (3/4π)^{1 / 3}. The stage when the equations lose their validity corresponds to a certain particle radius r_c, the critical radius. This can be computed as

follows. For the first four cases, the critical face radius a can be expressed in terms of the critical particle radius as a_c = δr_c with the constants δ (equal to $\mathrm{1}\mathit{/}\sqrt{\mathrm{2}}$, 1/2, 1/2, and 1/2, respectively). Substituting this into Eq. (A.8) we obtain ${\mathit{r}}_{\mathrm{c}}\mathrm{=}{\left[\mathit{\alpha }\mathrm{(}\mathrm{2}\mathrm{+}{\mathit{\delta }}^{\mathrm{2}}\mathrm{)}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\delta }}^{\mathrm{2}}\mathrm{)}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{-}\mathit{\beta }\right]}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}}\mathit{.}$(A.28)From Eq. (A.10) we obtain the relative initial densities D₀ (see Table 1). The critical density D_c is obtained by substituting r_c and a_c = δr_c into Eq. (A.9): ${\mathit{D}}_{\mathrm{c}}\mathrm{=}{\mathit{\beta }}_{\mathrm{1}}^{-1}{\mathrm{(}{\mathit{r}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{-}{\mathit{\delta }}^{\mathrm{2}}{{\mathit{r}}_{\mathrm{c}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}}\mathrm{=}{\mathit{r}}_{\mathrm{c}}^{-3}{\mathit{\beta }}_{\mathrm{1}}^{-1}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\delta }}^{\mathrm{2}}\mathrm{)}}^{\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}}\mathit{.}$(A.29)At the transition from the first to the second deformation stage of the BCC packing, contact points between the particle and the six rectangular faces of the unit cell form. Using δ₁ = 1/2, the particle radius at this moment of deformation r_c1 and the first critical density D_c1 follow from Eqs. (A.8) and (A.9), or using δ₁ = 0.5 and a₂ = 0, from Eqs. (A.19) and (A.21), respectively. At the end of the second stage, all fourteen faces come in contact. Here, with ${\mathit{\delta }}_{\mathrm{1}}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$ and δ₂ = 1/3, the second critical radius r_c2 and the second critical density D_c2 follow from Eqs. (A.19) and (A.21).

After the respective critical density is reached, the geometric approach to the computation of the effective stress σ₁ can no longer be applied (i.e. for D ≥ D_c or $\mathit{D}\mathrm{\ge }{{}^{\mathrm{(}}\mathit{D}_{\mathrm{c}}{}^{\mathrm{)}}}_{\mathrm{2}}$). In this case we use a continuous linear relationship between ${}^{\mathrm{(}}\mathit{D}_{\mathrm{c}}\mathit{,}{\mathit{\sigma }}_{\mathrm{1}}{{}^{\mathrm{(}}\mathit{D}_{\mathrm{c}}{\mathrm{)}}^{}}^{\mathrm{)}}$ and ${\mathrm{(}}^{}\mathrm{1.0}\mathit{,}{\mathit{\sigma }}_{\mathrm{0}}{\mathrm{)}}^{}$.

Acknowledgments

We thank the reviewer Satoshi Okuzumi for his constructive and helpful suggestions and comments on the paper. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) Priority Programme 1385 “The First 10 Million Years of the Solar System – a Planetary Materials Approach” and by the Helmholtz Association through the research alliance “Planetary Evolution and Life”.

References

All Tables

All Figures

	Fig. 1 Relative stress σ1/σ0 for the adopted four basic systems of packing plotted against the relative density D = 1 − φ (after Rao & Chaklader 1972, their Fig. 6). The lines are truncated at the respective critical density Dc (see Table 1), above which the effective stress is computed using a linear relationship between the value at Dc and the lithostatic pressure.
In the text

Fig. 2 Porosity after the cooling of the body below 300 K as a function of depth. Left panel: $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$, models Bcreep (solid lines) and Bpar (dashed lines). Right panel: $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{5}\mathrm{km}\mathit{,}\mathrm{2.3}\mathrm{Ma}\mathrm{\right)}$, models Screep (solid lines) and Spar (dashed lines). For each model four different packings are considered. From top to bottom (at the depth of zero km) the lines (both solid and dashed) correspond to SCP, ORP, BCCP, and RHP. We note that for the single packings the different porosities in the primordial layer result in a variation of the radius (see Table 3). In addition, for a chosen packing with different final average porosities (i.e. different thickness of the porous layer), the associated solid and dashed lines have different lengths. This is not obvious here because the dashed lines are overlaid by the solid lines at depth.

In the text

	Fig. 3 Central temperature (solid line) and central porosity (dashed line) as functions of time for the case Screep with RHP.
In the text

	Fig. 4 Porosity as a function of depth for $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$ for model Bcreep assuming ORP. The lines correspond to the adopted values of the activation energy E within the interval inferred by Schwenn & Goetze (1978).
In the text

Fig. 5 Porosity as a function of depth for $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$ for model Bcreep assuming ORP. The lines correspond to the adopted values of the grain size b. The grain size varies by two orders of magnitude between b = 10-7 and b = 10-5 m.

In the text

Fig. 6 The central temperature (solid lines) and the central porosity (dashed lines) as functions of time. Left panel: case $\mathrm{\left(}\overline{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{50}\mathrm{km}\mathit{,}\mathrm{3.2}\mathrm{Ma}\mathrm{\right)}$ with ORP and the associated effective stress σ1. Right panel: same as left panel but lithostatic pressure σ0 used instead of the effective stress σ1.

In the text

	Fig. A.1 General deformation geometry at any contact point between two grains.
In the text

	Fig. A.2 Unit cells and their intersections with an inscribed grain for the SCP and the ORP (top row from left to right) along with the top views of the unit cells and of the grain at its initial and critical (final) stages of deformation (bottom row).
In the text

	Fig. A.3 Unit cell of the RHP and its intersections with an inscribed grain (left panel), the top right view of the unit cell and of the grain at its initial and critical (final) stage of deformation (central panel), and a single face of the unit cell at the critical stage of the deformation (right panel).
In the text

Fig. A.4 Unit cells and their intersections with an inscribed grain for the two deformation stages of the BCCP (top row, first and second panels); the faces of the unit cell at the second critical deformation stage (top row, third and fourth panels); topview of the unit cell before the first critical stage (bottom row, left panel) and after the first but before the second critical stage (bottom row, central panel); and the arrangement of the spheres at the first critical deformation stage as seen in the plane defined by the centres of the three spheres (bottom row, right panel).

In the text