Comet 67P/Churyumov–Gerasimenko: Hardening of the sub-surface layer

From the operation of the MUPUS thermal probe Spohn et al. (2015) concluded that the material of the nucleus of 67P/Churyumov–Gerasimenko is likely to have a high strength, at least locally at the Philae landing site. In this work we consider the derived strength of the material in order to constrain its granulation. For this purpose we performed numerical simulations of the long-term sintering of ice–dust granular mixtures of different granulation, covered by a dust mantle. The dust mantle has a thickness of 0–16 cm, and a (pore size and temperature-dependent) thermal conductivity. According to our simulations a hardened layer at least a meter thick forms beneath the dust only when the grains are tens of microns in radius, or smaller. 2015 Elsevier Inc. All rights reserved.


Introduction
After the final touch-down of the Rosetta lander Philae the MUPUS PEN probe was deployed in an attempt to measure, among other parameters, the local strength of the nucleus of Comet 67P/ Churyumov-Gerasimenko. This probe, part of the MUPUS package (MUltiPUrpose Sensors for Surface and Sub-Surface Science) is equipped with 16 temperature sensors to be inserted into the soil by a hammer mechanism. After deployment, MUPUS PEN probe was nominally commanded to start the hammering sequence. The progress of the hammering was measured by a sensor. The depth sensor recorded an initial progress of about 20 mm and then its readings oscillate for 3.5 h by a total of 5 mm, suggesting that no further penetration took place. The findings of the MUPUS experiment have been described in detail by Spohn et al. (2015).
This work is an attempt to estimate the tentative thickness and strength of a layer of icy surface material covered by a dust mantle. This material is likely to have undergone a process of hardening due to sintering of ice grains in the ice-dust mixture. In the case of the MUPUS experiment, the measured strength should be considered as the crushing strength, as in the laboratory measurements described in Grün et al. (1993). Laboratory studies of the sintering process under cometary-like conditions include KOSI, a series of comet simulation experiments. It was found that after 45 h' irradiation of a fluffy dust-ice-sample with 1-1.6 solar constants a porous dust/ice crust with a crushing strength exceeding 5 MPa was formed below a thermally insulating dust mantle (Grün et al., 1993). The tensile strength is typically an order of magnitude lower than the compressive strength and is more temperature sensitive (Schulson and Duval, 2009;Kimberly et al., 2012). According to Schulson and Duval (2009), the ratio between the compressive strength and the tensile strength is 8 at À10°C and 20 at À50°C. On comets, tensile strength controls the ejection of grains, or agglomerates of grains. Ejection is possible only when the gas pressure beneath exceeds the tensile strength of the material. In our work we discuss both tensile strength and crushing strength. Our model parameters were chosen to reflect the conditions at the Philae landing site on Comet 67P/Churyumov-Gerasimenko.
The sintering of ice and other materials has been investigated by numerous authors. Unfortunately, most researchers investigated the sintering of ice mainly in the context of terrestrial snow, e.g. Salm (1982) and Colbeck (1998Colbeck ( , 2001. Thus, the corresponding equations describe granular ice with neutral gas in its pores, with a gas pressure well above the triple point. Pores of cometary ice contain vapor, but not a neutral atmospheric gas. Therefore, a more suitable model for cometary ice has to use sintering formulas derived theoretically, without the assumption that any neutral atmospheric gas is present in the pores, e.g. Swinkels and Ashby (1980) and Kossacki et al. (1994)

Physical approach
Comet 67P/Churyumov-Gerasimenko has a long and complex orbital history (Krolikowska, 2003). We focus our attention on the time span between the last major change of orbit in 1959, when the heliocentric distance at perihelion decreased from more than 3 AU to 1.29 AU, and the 2015 perihelion (8.5 orbital periods) (Krolikowska, 2003). In this time frame the orbital elements, and orientation of the rotation axis in space, are nearly constant. Thus, we consider them to be fixed values, while we varied the properties of the nucleus material.
A brief description of the method used in this work is given below.
-Calculations are performed for one selected location on the nucleus. Although the exact location of Philae's final landing site is unknown, it is likely to be closer to the equator than the initial landing site, i.e. the site of the first touch down. Two locations are considered for the purpose of simulations: one at the latitude of the initial landing site, 15°N, and one at the equator.
-The physical structure of the surface at Philae's final landing site is currently not known, neither is the exact orientation of Philae. For simplicity, we therefore assume the surface to be horizontal and smooth in this work. If the surface investigated by MUPUS is significantly inclined, and shadowed most of the day, the material should sinter slower than indicated by the presented simulations. Thus, when our simulations predict negligible sintering, the real process should happen even slower. -In our model, the uppermost layer of the nucleus is composed of dust only. This layer is hereafter called the dust mantle. Underneath the dust there is a mixture of crystalline H 2 O grains and dust. The model allows taking into account other components, such as amorphous H 2 O, and CO ice, but these are not considered in this work. -Granulation of the dust mantle can be different from the underlying material. -The material underneath the dust is composed either of ice grains with dust cores (Model A), or of agglomerates of grains (Model B). In the latter case the agglomerates are of two types; they are either made out of dust grains, or out of ice grains. -In Model B, the sizes of agglomerates and their internal porosities w int remain the same, and porosity w describes the space between agglomerates. The internal porosity of agglomerates affects their density.  -The dust grains do not sinter. In the material underneath the dust mantle, the sintering process leads to a depth dependence of the hardness. The rate of sintering of ice grains depends on their sizes. Small ice gains sinter much faster than large ones. Thus, individual ice particles may sinter to form agglomerates much faster than the agglomerates can clump together through sintering. The rate of sintering is crucially dependent on the local temperature. In the uppermost centimeters near the surface of a cometary nucleus the temperature can vary by several tens of degrees°due to the time-dependent absorption of the solar energy. The presence of a dust mantle may have a complex influence on the temperature. A low-albedo dust mantle of low thermal conductivity may warm up to a temperature higher than clean exposed ice. However, a thick, low-conductivity dust mantle will significantly reduce the heat flux transported downward through the dust mantle into the ice. The process of sintering is described in greater detail in Section 2.2. -Generally, the strength of cometary granular material is likely to be proportional to the ratio between the grain-to-grain contact areas and the cross sections of the grains (Skorov and Blum, 2012;Thomas et al., 1994).

Mathematical formulation
The model used in this study is an improved version of the program used by Kossacki and Szutowicz (2008) to simulate the emission of water from Comet 9P/Tempel 1 due to sublimation beneath the dust mantle and, most recently, to simulate observed recession of the surface of the same comet (Kossacki, 2015). The most important feature of this model is the inclusion of the sintering of ice grains, and the temperature-dependent sublimation coefficient of ice. The previous version of the model has been described by Kossacki (2015).
The model used in this work includes: -temperature dependence of the thermal conductivity of dust; -evolution of the nucleus cohesion due to sintering of ice grains (the so-called Kelvin effect which modifies grain-to-grain contact areas but does not affect the degree of compaction); -changes of porosity due to sublimation/condensation in the medium; -crystallization of amorphous water ice (if present); -sublimation of the CO ice (including explosive sublimation); -sublimation of H 2 O ice covered by the dust mantle, taking into account the temperature-dependent sublimation coefficient (Kossacki and Leliwa-Kopystynski, 2014); -diffusion of the vapor through the dust and the resulting recession of the surface; -illumination dependent on the local orientation of the surface; -the material beneath the dust mantle can be either: ice grains with non-volatile cores (Model A), or agglomerates of grains (Model B). When the surface is not covered by dust, the surface temperature T s is given by the equation: where S c is the solar constant, R h is the actual heliocentric distance in AU, A is the surface albedo and a is the zenith angle of the Sun, i.e.
the angular position of the Sun in relation to the local normal at a given point on the surface. r denotes the Stefan-Boltzmann constant, and the emissivity. The term rT is the temperature gradient at the surface. k denotes thermal conductivity, which will be described later in this section. H is the latent heat of sublimation and F s is the flux of molecules subliming and subsequently escaping into space.
Equations for the surface temperature T s of the dust mantle, and the temperature beneath the dust are: Here, the index d indicates the properties of the dust, i.e. k d refers to the thermal conductivity of dust, and F sd is the flux of molecules subliming and diffusing into space through the dust mantle. The index i indicates the dust-ice interface. The temperature gradients are positive when the temperature increases versus depth.
The thermal conductivity of the dust mantle, k d , and the thermal conductivity of the underlying material, k, are temperature dependent. The equation for the thermal conductivity of dust is Here k d0 is the thermal conductivity, which characterizes the heat transport within solid matrix of dust grains, but not in void space between grains (a value usually applicable at low temperatures, when heat transport by radiation is less relevant). The radiative term in Eq. (4) is important only at high temperature. When pore radius r p ¼ 100 lm and T ¼ 250 K, and k d0 ¼ 2 mW m À1 K À1 the radiative heat transport is <20% of conductive heat transport.
The ice-dust material underlying the dust mantle has thermal conductivity The terms on the right hand side of Eq. (5) describe: heat transport through the solid matrix of grains, and heat transport through the pore space by vapor. k g is the thermal conductivity of the grains. At low temperatures the second term on the right hand side of Eq. (5) can be neglected. For more details see Kömle and Steiner (1992). When the grains are composed of solid ice k g refers to the thermal conductivity of pure ice. When the grains have dust cores, (Haruyama et al., 1993). The symbols k c and k i denote the thermal conductivity of the mineral cores and of the bulk ice, respectively. The flux F s of vapor molecules escaping to space can be calculated in different ways. When the dust mantle is thick compared to the size of pores (thicker than ten pore radii), one can use the so-called long tube model. In this case (Mekler et al., 1990). Here r d is the characteristic radius of the pores, v dm is the volume fraction of the dust in the mantle (porosity of the dust layer w d ¼ ð1 À v dm Þ), D d is the thickness of the dust, and s is the tortuosity of the pores. l is the molar mass of water, R g denotes the universal gas constant, and p sat is the gas pressure at phase equilibrium at the given temperature T. Eq. (7) ignores factors such as the shapes of ice crystals, the structure of their surfaces, and the presence of impurities. All these factors can be included in the temperature dependent sublimation coefficient a s (Kossacki et al., 1999;Gundlach et al., 2011;Kossacki and Leliwa-Kopystynski, 2014) and references therein. After including the sublimation coefficient, Eq. (7) becomes When the dust mantle is thinner than ten pore radii Eq. (8) should be replaced with the modified Clausing formula When compared to the original version of the Clausing formula the sublimation coefficient is added. The value of the Hertz factor h evolves over time, due to the sintering of grains. The rate of sintering depends on the temperature, the radii of grains r g , and also on the current value of the Hertz factor h itself. In the case of coarse-grained material, sintering is slow and h is always close to its initial value, but the sintering of small grains can be very fast, depending on the temperature. There is a number of known sintering mechanisms. For H 2 O ice under cometary conditions the most efficient one is the so-called Kelvin effect, i.e. the transport of molecules via the vapor phase from adjacent grains onto the connecting neck. The difference between local surface curvatures causes a difference of the phase equilibrium pressure. At the neck it is lower than at the grains. Thus, molecules sublime from the grains and condense on the neck. It is important to note that this type of sintering mechanism does not lead to any  changes of porosity. The radii of the grains also remain unchanged. However, the grains change their shapes. The rate of increase of the contact area as described by h is (Swinkels and Ashby, 1980;Kossacki et al., 1994). The symbol r n denotes the radius of the grain-to-grain contact area (radius of the sinter neck). The parameters c and X are the surface energy and molar volume, respectively. For water ice, X 2 c ¼ 4 Á 10 À11 m 4 mole À2 J. The dimensionless variable S describes the fraction of the surface of a grain from which the material is removed (source area) and transported onto the sinter neck (sink area). Eq. (10) was developed for sintering under isothermal conditions. Presence of a thermal gradient results in sublimation of ice at the warmed surface and subsequent migration of molecules toward the cold, deeper layers. In such a case the molecules condense both on the grains and on the connecting necks. Thus, the ratio of the cross-sections, i.e. the Hertz factor, can be expected to undergo small changes. Indeed, Eq. (10) was tested using laboratory experiments dealing with comet-like samples under conditions far from isothermal. Radiative warming of the surface and cooling of the bottom resulted in the presence of thermal gradients exceeding 10 K cm À1 (Kossacki et al., 1997). In real cometary nuclei gradients of that magnitude should not be expected, except in a dust mantle of very low thermal conductivity. Note that there is a range of other sintering mechanisms which are not considered in this paper because they are of no importance in our case. These are: surface diffusion from a surface source, lattice diffusion from a surface source, grain boundary diffusion from a boundary source, lattice diffusion from a boundary source, and lattice diffusion from dislocation sources. A parameter requiring some attention in this context is the density. In Model A the sizes of grains composing the dust mantle are the same as the grains in the underlying icedust material. The average density is where . i is the bulk density of ice, and . d is the bulk density of dust.
In Model B the material beneath the dust mantle consists of agglomerates of the same sizes, but composed either of dust grains, or ice grains. Thus, the average density of the ice-dust material is where v ai , and v ad denote the volume fractions of different types of agglomerates in the material. The symbol w int denotes the internal porosity of agglomerates. The void space between agglomerates is given by The key feature of this model is the calculation of the depthdependent strengthening of the material due to the sintering of grains. As outlined above, strength may have different meanings. We consider the tensile strength and the crushing strength. Skorov and Blum (2012) calculated the strength of layers composed of dust agglomerates as a product of the packing density U a of aggregates, their intrinsic tensile strengths, and the ratios between the contact areas and the cross sections of aggregates (i.e. the Hertz factor). Thus, when the material is composed of aggregates of ice particles, which sinter together, two different Hertz factors need to be considered: h int which describes bonding of particles within an aggregate, and h a which describes effective contacts between aggregates. Thus, the equation for the tensile strength of the material should be where w int is the internal porosity of aggregates, and w a ¼ 1 À U a (the packing density of aggregates). In our Model B the material underlying the dust mantle is composed of two types of aggregates. The aggregates of dust particles do not change their strength, while the ice aggregates undergo strengthening. Thus, after some time dust aggregates can be much weaker than ice aggregates. We estimate the effective strength assuming that the dust aggregates within the ice-dust material behave as voids. In that case Eq. (14) becomes However, it should be noted that our approach to estimate the strength of the material composed of different aggregates needs verification by laboratory measurements. When the material is composed of solid particles (rather than aggregates) Eq. (14) becomes where r Tg is the tensile strength of a single grain. In our Model A the grains contain not only ice, but also mineral cores. Thus, r Tg is larger than for pure ice.
The crushing strength of a granular material r C is related to that of a compact material r Ci by (Thomas et al., 1994).
The simulations were performed for a layer 20 m in thickness. At the bottom of the model grid the temperature is determined by the zero flux boundary condition. The initial temperature is a free parameter.

Parameters
The orbital parameters of the model nucleus were chosen to represent 67P/Churyumov-Gerasimenko. The semimajor axis is 3.5029 AU, the eccentricity is 0.6319, and the orbital period is 6.5563 years. The dust mantle is characterized by five parameters: the albedo A, the emissivity , the porosity w d , the thermal conductivity at low temperature k d0 and the thickness of the dust layer D d .
The observational values of albedo are very low. The normal albedo derived from VIRTIS measurements is 0:060 AE 0:003 at 0.55 lm (Capaccioni et al., 2015). Similarly, geometric albedo of the entire nucleus derived from OSIRIS data is 5:9 AE 0:2 at 550 nm (Sierks et al., 2015). In this work the surface albedo A of the dust covered surface is 0.01, to obtain an upper estimate of the energy absorbed by the surface. The emissivity is 0.9.
The remaining parameters are not well known, but some estimates are available. According to Krause et al. (2011), a dust layer created in vacuum by random ballistic deposition may have a very large porosity of w d ¼ 0:85 and a thermal conductivity as low as k d ¼ 2:6 mW m À1 K À1 . Krause et al. used small silicate grains in their experiments. The grains were 1.5 lm in radius, with a bulk thermal conductivity k d ¼ 1:4 W m À1 K À1 and a specific heat c d ¼ 840 kg m À3 . In this work we take into account the temperature dependence of the thermal conductivity of the dust. The parameter k d0 which describes the thermal conductivity at low temperature, is within the range 1-10 mW m À1 K À1 . The values of k d , and c d are the same as reported by Krause et al. (2011).
The material underneath the dust layer is characterized by the volume fractions of the respective components, their bulk parameters, and the Hertz factor. The volume fraction of ice is v i ¼ 0:4, and the volume fraction of the dust component is v d ¼ 0:1. The volume v d is constant, but v i evolves due to the sublimation and condensation of ice. However, the changes of v i are small due to small temperature gradient beneath the dust mantle. The properties of the dust particles in the ice-dust matrix are the same as in the dust mantle. The initial value of the Hertz factor h 0 is within the range 0.001-0.01.

Temperature and density
In Fig. 1 we show profiles of the temperature versus depth at perihelion, at local noon. The grain radius r g is 1.5 lm; the radius of pores in the dust layer r d is 3 lm; the initial Hertz factor h 0 is 0.001; and the initial temperature is 80 K. The curves are drawn for orbital periods: 1, 2, and 4. It can be seen that an almost isothermal layer forms beneath the dust mantle, with temperature gradients of the order of only 1 K m À1 . The diurnal changes of the temperature are visible at depth smaller than 0.2 m. Fig. 2 is analogous to Fig. 1. The radii of grains and pores are 10 times larger: the grain radius r g is 15 lm and the radius of pores in the dust layer r d is 30 lm. Comparison of Figs. 1 and 2 indicates, that when the material is composed of larger grains the temperature beneath the dust mantle is significantly higher. This can be due to slower sintering of larger grains and hence slower thickening of sintered layer of enhanced thermal conductivity. We investigated also evolution of the average density due to the downward migration of vapor, as well as erosion of the surface due to the escape of subliming molecules. Both effects were negligible, except the situation when the dust mantle was not present on the surface. In that case the volume fraction of ice increased due to condensation of vapor, but less than 5%. Small increase of density was due to erosion of the surface. When the material was fine grained, r g ¼ 1:5 lm, the surface receded by 1-10 cm per orbital period, depending on the albedo.

Strengthening
In Fig. 3 we show the Hertz factor h versus depth and the rate of sintering, dh=dt, versus depth. The model parameters are the same as in the case of Fig. 1. Comparison of Figs. 1 and 3 indicates, that the Hertz factor grows in the layer, where the temperature gradient is small. Further below, where the temperature gradient is high, the temperature is so low that the sintering rate is insignificant. The sintering is slow also just beneath the dust mantle, where the Hertz factor has already exceeded the value 0.3.
In Fig. 4 we show the Hertz factor h versus depth and dh=dt, versus depth. The model parameters are the same as in the case of Fig. 2. Comparison of Figs. 2 and 4 indicates, that the efficiency of sintering is high only in a layer characterized by a small temperature gradient, as in the case of fine grained material.
In Fig. 5 we show the thickening of the hardened sub-dust layer for a dust mantle that is 4 cm thick. The hardened layer is considered as the layer where the Hertz factor exceeds a threshold value h th , which is 0.33. This value is approximately the limit of efficiency of the sintering mechanism considered in this work. The location assumed in our model is at the latitude 15°N. Plotted is the thickness of this hardened layer versus time for cases in which: the grain radius r g is 1.5 lm, and 15 lm; the radius of pores in the dust layer r d is 3 lm and 30 lm; the initial Hertz factor h 0 is 0.001 and 0.01. The initial temperature T 0 is 80 K and 40 K. T 0 depends on the heliocentric distance where the comet formed which is why we considered two values for T 0 . It can be seen that the granulation of the material is very important. When the material is fine-grained (r g ¼ 1:5 lm) and unconsolidated (h 0 ¼ 0:001), the hardened layer grows in thickness to 5 meters within just 2-3 orbital periods, depending on the initial temperature. When the radii of grains and pores are ten times larger, or the initial Hertz factor is hundred times higher, the sintered layer grows only to 2 m thickness during the considered period. When r g ¼ 15 lm and k d = 2:6 Â 10 À3 þ 4r d rT 3 , or when r g ¼ 1:5 lm and k d ¼ 1 Â 10 À3 þ 4r d rT 3 the layer of h > 0:33 does not form during 8.5 orbital periods. In the latter case the time of 10.5 orbital periods is needed for the formation of a 10 cm thick hardened layer. The initial Hertz factor h 0 , as seen in Fig. 5 has only minor significance for the formation of a hardened subsurface layer. The hardened subsurface layer grows in thickness at almost the same rate at both 40 K and 80 K initial temperature. Fig. 6 shows results obtained when the model location is at the equator. The profiles can be directly compared to these in the upper panel of Fig. 5. When the sintering is fast, the exact latitude has minor significance. At the equator the initial thickening is faster than at the latitude 15°, but later it slows down and the overall progress of the sintering over the considered period is similar. When the sintering is slow its progress at the equator can be more prominent than at the latitude 15°N.
Imagery shows that the surfaces of cometary nuclei can be smooth, or rugged. We therefore have to consider a scenario in which the surface is shadowed. Fig. 7 demonstrates the influence of local shadowing. The model parameters are the same as in the case of Fig. 5, but the solar flux is reduced by 25%. The thickening of the hardened layer is noticeably slower, which emphasizes the significance of the granulation of the material.
In Fig. 8 the thickness of the hardened sub-dust layer versus time is shown when the dust mantle is not present, or very thick, namely 8 cm and 16 cm. The latitude is b ¼ 15 N. The hardened layer of h > 0:33 forms both under a very thick dust mantle, and under uncovered surface. However, when the dust is 16 cm thick, and k d0 is as low as 2.6 mW m À1 K À1 the rate of sintering is so low, that in the considered time span the Hertz does not grow to the threshold value at any depth. Fig. 9 shows results obtained when the entire model volume of the nucleus, i.e. both the dust mantle and the underlying medium, is composed of agglomerates of small grains. We have plotted the thickness of the hardened layer versus time. In Model B the material underlying the dust mantle is characterized by two Hertz factors, h int and h a . Thus, we consider material as hardened when both coefficients exceed the threshold value. The aggregates have internal porosity w int ¼ 0:3, and the porosity w a ¼ 0:3, so the average density of the material is the same as in Model A (Fig. 5).

Model B
The values of the remaining model parameters are: r g = 15 lm, r d = 30 lm, h 0 = 0.001, and T 0 ¼ 80 K. A comparison of the profiles shown in Figs. 8 and 9 indicates that replacing the solid grains with aggregates of grains reduces the thickening rate of the hardened layer by more than an order of magnitude. When D d ¼ 16 cm thick, and k d ¼ 10 Â 10 À3 þ 4r d rT 3 , the hardened layer grows in thickness to more than 5 m in 3 orbital periods, or only to 0.4 m in 8.5 orbital period depending on the structure of the material. Table 1 gives the state of the hardened sub-dust layer of Comet 67P/C-G after 8.5 orbital revolutions . We give the following values: Hertz factor just beneath the dust mantle, dimensionless tensile strength, dimensionless crushing strength (only for Model A), and depth where the Hertz factor exceeds 0.33. The strengths are normalized to the appropriate bulk strengths of the material, i.e. the H 2 O ice with dust grains. The selected value of the Hertz factor is close to the efficiency limit of the considered sintering mechanism which is most effective under cometary conditions, i.e. the transport of molecules via the vapor phase from adjacent grains onto the connecting neck. Note that this does not preclude further progress of sintering due to the other mechanisms mentioned in the previous section. As the aim of this work is to determine what material structure is necessary for sintering-induced hardening of the material only consideration of the first stage of sintering was needed and we focused on the dominant sintering mechanism during this stage.

Summary and discussion
The mechanical properties of solid H 2 O ice have been investigated by various authors, but in most cases at high temperatures. In the temperature range 240-270 K pure poly-crystalline ice has a tensile strength of about 1 MPa (Hawkes and Mellor, 1972;Currier and Schulson, 1982;Kimberly et al., 2012). Both tensile strength and compressive strength are temperature dependent (Schulson and Duval, 2009;Kimberly et al., 2012). At a temperature of about 100 K the tensile strength of poly-crystalline ice is about 2.5 MPa (Kimberly et al., 2012). The compressive strength of solid H 2 O ice made of 1 mm grains is >30 MPa at 200 K, and >70 MPa at 100 K (Schulson and Duval, 2009). Thus, at the temperature 100 K, the granular ice of porosity w ¼ 0:5, with a Hertz factor h ¼ 0:33, the tensile strength is r T ¼ 0:4 MPa, while the compressive strength is r C ¼ 8 MPa. In the case of ice of loosely bound grains, h ¼ 0:001, the tensile strength should be only 1.25 kPa and the compressive strength 25 kPa. For comparison, the estimated tensile strength of the nucleus of Comet 9P/Tempel 1 is within the range 0-12 kPa (Holsapple and Housen, 2007), or above 10 kPa (Reach et al., 2010). It should be noted that the presence of dust, or sand, in the ice may significantly increase the strength. According to Kimberly et al. (2012) the strength of ice with a 20% mass fraction of basalt and ammonium sulfate at 251 K is two times higher than of pure ice, while the strength of ice with urea is three times higher than of pure ice. This indicates that at 100 K ice with an admixture may have a tensile strength of approximately 10 MPa.
In order to better simulate the topography at the Philae landing site, we simulated the evolution of the material when the solar flux is reduced by 25%. In both cases, also when the illumination was deliberately overestimated, our simulations indicate that the coarse-grained material should not undergo significant sintering. However, the failed attempt to penetrate the subsurface later by MUPUS-PEN probably indicates high strength.
Putting this into the context of our model, the conclusion is that the original material must have been fine-grained rather than coarse in order to produce this hardened layer.

Conclusions
We performed numerical simulations of the long-term sintering of ice-dust granular mixtures of different granulation in order to constrain the granulation of the Rosetta target Comet 67P/Chury umov-Gerasimenko. We find that a hardened layer at least a meter thick forms beneath the dust. This hardening can only be observed in the models when the ice grains are smaller than a few tens of microns in radius.