© ESO 2018

1 Introduction

Fracturing is prevalent on many scales on Comet 67P/Churyumov-Gerasimenko (hereafter 67P) when observed by Rosetta’s OSIRIS imaging instrument. Large fractures extending up to several hundred metres in length are obvious in several places on exposed consolidated terrains (see Thomas et al. 2015; El-Maarry et al. 2015a for a detailed description of terrain types), such as in the lineaments on the Hathor cliffs or in the ~500 m feature near the neck (El-Maarry et al. 2015b). On smaller scales (metres to decimetres), irregular fracture networks are seen globally, covering many of the consolidated and brittle terrain types, as well as many cliff faces. Fractured cliffs are commonly associated with mass wasting and landslides (Pajola et al. 2015), and fracturing may represent an important erosion process. Pervasive, semi-regular networks of intersecting fractures have been likened to thermal contraction crack polygons, also seen on Earth and Mars (El-Maarry et al. 2015b,a; Auger et al. 2015, 2018). Polygonal fractures and other patterns often overlay each other, implying different formation mechanisms or timeframes. Morphologically, most individual fractures resemble tensile, “opening” features (El-Maarry et al. 2015b), but a population of fractures associated with the neck region have been identified as a possible shear zone (Matonti et al. 2016). Fractures are also present on the surface of the boulders which litter the cometary surface, with examples down to a few metres in size covered in fractures that are from tens of centimetres to metres long (Pajola et al. 2015) and often surrounded by fragmented remains. An erosion sequence of chunks of consolidated cliff material into ever smaller boulders is implied; Pajola et al. (2015) have suggested a balance between this process and sublimation, leading to the observed steady-state distribution in boulder sizes. At the smallest scales, the Philae lander’s CIVA camera saw tens of cm to sub-cm fractures covering the boulders and blocks making up the Abydos landing site (Poulet et al. 2016). They are noted to be neither linear nor regularly polygonal, but with a morphology suggestive of a slow stress evolution.

Such a variety of fracture types likely have a range of formation mechanisms; rotational, collisional, and thermal mechanisms have all been proposed. Aligned fracture systems in the neck may be associated with torques from rotation- or orbit-induced stresses or, alternatively, with stresses from a putative collision and merger of the two comet lobes (El-Maarry et al. 2015b; Matonti et al. 2016). Desiccation stresses, from changes in volume caused by loss of volatiles, are important on Earth and may also play a role where outgassing is occurring on the comet. Other features, especially the polygonal networks, are more likely caused by thermal stresses, as evidenced by their uniformity over the nucleus and their tensile morphology, and supported by the large changes in temperature expected on cometary surfaces. El-Maarry et al. (2015b) compute diurnal temperature changes of ~230 K for 67P at perihelion and this is similar in magnitude to the Moon, where thermal weathering has been shown to be important (Molaro et al. 2015). Alí-Lagoa et al. (2015) also invoke thermal weathering to explain activity on the comet starting before perihelion. Thermal stresses may well be the dominant mechanism for the formation of polygonal and irregular fractures (those not associated with the neck) on both large (Thomas et al. 2015; El-Maarry et al. 2015a) and small scales (Poulet et al. 2016) on 67P, and the remainder of this paper will focus on these processes.

Thermal stresses are induced by changes in temperature causing the expansion and contraction of material (in both uniform samples and mixtures with varying thermal expansion coefficients) leading to breakdown in a process known as thermal insolation weathering. If stresses exceed a material’s tensile strength at any time, failure will occur immediately, leading to a thermal shock fracture. Alternatively, repeated stress cycles below the failure stress can lead to gradual weakening and the development of micro-fractures, leading to eventual failure in a process known as thermal fatigue (Hall & Thorn 2014; Delbo et al. 2014). Thermal fatigue follows existing material weaknesses, often leading to fracturing parallel to the surface (flaking), whereas thermal shock cuts across existing weaknesses, forming distinctive rectilinear or polygonal fracture networks. Links between the two are poorly understood, but it is unlikely that thermal fatigue enhances the probability of fractures developing by thermal shock (Hall & Thorn 2014).

Thermal fatigue may be the dominant mechanism for the formation of regolith on asteroids (Delbo et al. 2014; Dombard et al. 2010). When subjected to diurnal temperature variations of the order of 100 K, rocks larger than a few centimetres are broken down by the gradual growth of fractures. Thermal fatigue may also be important on comets (Alí-Lagoa et al. 2015), but large temperature variations mean that immediate material breakdown by thermal shock may dominate. Kuehrt (1984) modelled Halley’s comet as solid water ice and showed that seasonal thermal stresses at perihelion can reach tens of MPa, easily exceeding the tensile strength. These results were robust when extended from merely elastic stresses to a viscous description (Tauber & Kuehrt 1987), and similar when considering inhomogeneous bodies, showing that comets should be thermally fractured and that this may be important for erosion and outgassing.

On Earth and Mars, stresses exceeding a material’s tensile strength over large, uniform areas lead to networks of fractures which intersect in polygonal patterns. Each fracture acts to dissipate stress to some distance around itself, depending on the material properties and fracture depth, thus determining the distance to the next fracture. Polygon size is therefore related to the material and thermal environment in the subsurface (Lachenbruch 1962). Thermal contraction crack polygons are seen in terrestrial glaciers (e.g. Marchant et al. 2002), bedrock (e.g. El Maarry et al. 2010), and ice-bonded soils (e.g. Mellon 1997; Lachenbruch 1962). Auger et al. (2018) mapped polygons on the surface of 67P, finding them on both comet lobes, in the northern and southern hemispheres, and with no particular correlation with latitude or other features. Their morphology closely resembled terrestrial and Martian features, strongly suggesting thermal fracturing as the origin (see e.g. their Fig. 12). On Earth, polygon formation is enhanced by the presence of water ice (because of its large coefficient of thermal expansion and because it further opens and erodes fractures by freeze-thaw weathering), but the presence of ice is not a prerequisite (Hall & Thorn 2014). Water ice would appear to be required on 67P in order to bond the otherwise poorly consolidated nucleus material together to form a layer competent enough for high stresses to develop. Auger et al. (2015, 2018) suggest that a hard sintered layer within a few centimetres of the surface would be consistent with the observed polygons of a few metres in size. We discuss the evidence for such a hard layer in Sect. 2. In Sect. 3, we introduce a thermo-mechanical model of the cometary surface, with the aim of further quantifying the thermal stresses on 67P and their variations. The results of this model are discussed in Sect. 4, followed by the implications and conclusions in Sect. 5.

2 Evidence for a hard layer

The presence of thermal contraction crack polygons implies a competent layer within centimetres or metres of the cometary surface. Additional evidence also supports this hypothesis and we discuss it now.

2.1 Constraints from laboratory experiments

Since the 1950s, several laboratory simulation experiments have been performed in order to understand the physical processes occurring on cometary nuclei. These experiments consist of producing centimetre-thick porous and fluffy samples made of particles of H₂O (and CO₂) ice and dust (silicates, organics, etc.) maintained in a cold thermal-vacuum chamber while irradiated with a Sun-simulator. The irradiation of these ice-dust samples triggers the sublimation of the ice, producing a dehydrated dust mantle at the surface. Underneath this mantle, a hardening of the remaining ice was observed in many experiments (Kochan et al. 1989; Ratke & Kochan 1989; Seiferlin et al. 1995; Thiel et al. 1989; Poch et al. 2016a,b). This consolidation can be explained by the re-condensation of H₂O/CO₂ between the individual ice and dust particles and/or by sintering: both processes build bonds connecting the ice particles together and cause solidification of the sample. Measurements of the hardness of these ice-dust crusts, performed after the KOSI (Kometensimulation, DLR) experiments, yielded values of 0.15–5 MPa (Gruen et al. 1993; Kochan et al. 1989). Recently, new simulation experiments performed at the Laboratory for Outflow Studies of Sublimating Materials (LOSSy) of the University of Bern showed that the hardening depends on the way the ice and the dust are initially mixed together: sublimation of an intimate mixture of ice and dust at the particle level (inter-mixture) leads to a higher solidification than a mixture where the dust particles are included in the grains of ice (intra-mixture) (Poch et al. 2016b).

2.2 Constraints from modelling

Kossacki et al. (2015) modelled the sintering of water ice grains in the upper layers of 67P using a thermophysical model of a mixture of ice and mineral grains under a dusty mantle (thickness of 0–16 cm), subjected to solar heating. Sublimation and recondensation of the water ice reshaped the grains, enhancing their surface contact area and increasing the tensile and compressive strengths. Depending on the initial grain size, thickness of the dust layer, and insolation conditions, Kossacki et al. (2015) found a sintered layer extending up to several metres below the dust mantle after a few orbital periods. Within this layer the strength was enhanced, relative to an unconsolidated mixture, by several orders of magnitude, so that it was comparable to a few times less than the ~MPa strength of solid water ice. Kossacki et al. (2015) use 1, 2.5, and 10 MPa tensile strengths for, respectively, pure poly-crystalline ice at 240– 270 K and at 100 K, and an ice-rock mixture of roughly 20% rock by mass at 100 K. This hardening occurred at both the locations studied (the equator and the Agilkia landing site) and is expected elsewhere whenever the insolation is high enough (making it less efficient in shadowed areas) and the initial grain size sufficiently small (tens of microns in radius or smaller).

2.3 Constraints from the lander

Two Philae instruments, SESAME-CASSE and MUPUS-PEN, provide direct clues about the mechanical properties of the near-surface layer of the comet at the Abydos landing site. At Abydos, CASSE listened to the hammering of the MUPUS PEN during its insertion phase. An evaluation of Rayleigh wave velocities at this site yields a surface layer of 10–50 cm thickness which has a Young’s modulus of at least 7 MPa, with a tentative upper limit of 980 MPa. This layer has a shear modulus in the lower third of the variability range of terrestrial snow, while its density possibly exceeds that of pore-free water ice. Both can be achieved by introducing a regolith component of loose grains. Underneath this layer, shear wave velocity is significantly reduced (32% or more), likely due to a further reduced shear modulus (Knapmeyer et al. 2017). This reduction might be the result of an incomplete sintering of the ice (and thus the absence of cementing of the non-ice components). Although there is no direct relation between elastic moduli and strength, some laboratory data on snow and porous ice/dust mixtures can be used to correlate Young’s modulus with uniaxial compressive strength. Jessberger & Kotthaus (1989) showed that for porous ice/dust mixtures (KOSI ice) and temperatures down to 123 K the elastic modulus is about two orders of magnitude higher than the compressive strength. Assuming that this relation is also applicable at Abydos, the compressive strength of the upper hard layer would be between 70 kPa and 9 MPa. This range is compatible with the lower limit of 2 MPa for the compressive strength deduced by Spohn et al. (2015) from the fact that the MUPUS PEN was not able to penetrate into the ground at Abydos. This is also consistent with the results from experiments and modelling above (bearing in mind that compressive strength is generally greater than tensile for relevant materials).

Finally, differences between the results of the SESAME permittivity probe and the CONSERT Rosetta-to-Philae radiowave experiment can be related to an increase in the porosity (i.e. volume fraction of vacuum) from the near-surface of the small lobe of the nucleus sounded by SESAME-PP to the deep material sounded by CONSERT (Kofman et al. 2015; Brouet et al. 2016). These results further support the presence of a hard layer at the location of Philae.

Together, these results point to a hard sintered or recondensed layer within tens of centimetres of the cometary surface and extending down of the order of a metre, with tensile and compressive strengths around a MPa and a Young’s modulus several orders of magnitude higher. Compared to this, loosely bonded granular ice should have a significantly lower tensile strength, of the order of hundreds of kPa down to ~1 kPa (Kossacki et al. 2015).

3 Modelling thermal stresses

3.1 Thermal model

We computed the temperature inside the nucleus as a function of time to derive the seasonal temperature trends. We used a spherical shape model for the nucleus, orientated accordingly to its pole with RA = 69.57^∘ and Dec = 64.01^∘ and a rotational period of P = 12.40 h (Jorda et al. 2016). Our thermal model takes into account the solar insulation, the thermal emission, and the sublimation of water ice on the nucleus surface, and the heat conductivity inside it (Groussin & Lamy 2003). We computed the temperature on the surface and inside the nucleus over one complete revolution, taking into account the diurnal and seasonal changes in insolation with heliocentric distance. To ensure convergence, we used a time step of 12.4 s and ran the thermal model over five complete revolutions, for a total of 8.1 millions time steps. We then computed, for each nucleus rotation (every 12.40 h), the average diurnal temperature at each depth interval, i.e. the seasonal trend. Our model has 2000 depth intervals of one-fifth of a diurnal skin depth, given by $\xi =I∕c\rho \sqrt{P∕\pi }$, where I is thermal inertia in units of J m⁻² K⁻¹ s^−1∕2, and c = 1000 J K⁻¹ and ρ = 532 kg m⁻³ (Jorda et al. 2016) are taken as the heat capacity and density of the cometary material. Overall, this allows us to examine the stress profile down to a total depth of between ~4 and ~20 m, depending on the chosen value of I between 10 and 1000 J m⁻² K⁻¹ s^−1∕2 (Sects. 4.1 and 4.2).

With the above thermal model, we could not use a realistic shape model for 67P, since even a very low-resolution shape model with a few thousand facets would require more than one month of computer time for each value of the thermal inertia. However, to overcome this numerical issue, we used a second thermal model, which takes into account the real shape of 67P and associated projected shadows and self-heating, but neglects thermal inertia (I = 0 J m⁻² K⁻¹ s^−1∕2). For this purpose, we used a decimated version of shape model SHAP7 (Preusker et al. 2017) with 6000 facets. For each facet, we computed the temperature on the surface of the nucleus over one complete revolution, taking into account the diurnal and seasonal changes in insulation. We then computed, for each facet, the average diurnal temperature, i.e. the seasonal trend. Because of the null thermal inertia, we can only calculate stress in an infinitesimally thin surface layer and the absolute values derived from here are unrealistic. This is a conceptual limitation of the model; nonetheless, relative stresses between facets should still reflect real differences in seasonal conditions and will be used to study variations in stress across the nucleus surface (Sect. 4.3).

3.2 Stress model

We follow the method of Mellon (1997), who modelled the formation of polygons in Martian permafrost. This represents an advancement over the semi-analytical approach of Kuehrt (1984) and Tauber & Kuehrt (1987), who used fixed rather than temperature varying parameters and only considered pure water ice and a single effective temperature profile. As described in Mellon (1997) ice, rocks, and frozen soils respond to applied stresses in a Maxwellian, viscoelastic way, i.e. with elastic and viscous deformations acting in series (see their Fig. 2 for a schematic representation). The total time-dependent strain can then be found by summing these two terms and, with the addition of a third term describing the thermal expansion and contraction, a thermo-viscoelastic model of the material is developed. Differentiating by time then leads to an expression for the total strain rate, $\stackrel{}{ϵ}$, as $$\stackrel{}{ϵ}={\stackrel{}{ϵ}}_e+{\stackrel{}{ϵ}}_T+{\stackrel{}{ϵ}}_v,$$(1) with the subscripts e, T, and v referring to the elastic, thermal, and viscous terms, respectively. Next, the assumption of a plain geometry extending in all horizontal directions is made. In this geometry the total horizontal strain rate must sum to zero as the material is restricted from expanding or contracting in this direction, and a resultant horizontal stress is induced. Using the standard elastic and thermal expressions and a creep model for the third term, as in Mellon (1997), and ignoring the second-order terms, we come to the equation $$\stackrel{}{ϵ}=0=\frac{1-\nu }{E}\stackrel{}{\sigma }+\alpha Ṫ+\text{sign}\left(\sigma \right)A_0{\mathrm{\text{e}}}^{-Q∕\mathrm{RT}}{\left(\frac{\sigma }{2}\right)}^n,$$(2)with the three terms, elastic, thermal, and viscous, as above. Here, σ is the resultant stress in the horizontal direction; ν is the material’s Poisson’s ratio; E and α are its Young’s modulus and coefficient of thermal expansion; T is the temperature; A₀, Q, and n are viscous material properties; and R = 8.31 J kg⁻¹mol⁻¹ is the gas constant. The sign function in the viscous term is an addition from Maloof et al. (2002) and ensures than the term acts in the correct direction, independent of the chosen n exponent. The convention adopted here is the same of that of Mellon (1997): positive values of σ indicate tension, while negative values indicate compressive stresses.

Equation (2) then represents the non-linear material response to a time-varying temperature profile, and must be solved numerically for σ. When this tensile stress exceeds the material’s tensile strength, fracturing will occur. We solve the equation using scientific Python’s fsolve routine for each time- and depth-step of a one-dimensional temperature profile. Before examining the stresses on 67P, however, we first confirm the validity of our numerical solutions by reproducing the results of Mellon (1997). Figure 1 shows the temperature and stress at the Viking 2 landing site on Mars at a depth of 0.5 m, from their Fig. 3. We use this temperature data as an input to our own code and, using identical material properties, produce the dashed stress curve shown in the lower panel. The overall shape and intensity of the stress is well reproduced; compressive (negative) stresses are encountered during summer as temperatures rise, while tensile (positive) stresses dominate with the falling winter temperatures. Small differences between the curves (RMS residuals of 746 kPa) may be due to a difference in time-step, which is unspecified in Mellon (1997).

Having validated the code, we now run it for a seasonal cycle of the comet thermal model. This is repeated a number of times (up to several hundred) to ensure convergence. At large depths, where temperature variations are very slow, the code can still fail to fully converge. However, the trend here is for ever decreasing stress with time; furthermore, we check each run carefully to make sure that the depths we are interested in, where the stresses change from tensile to compressive, are not affected.

For our baseline mode, we adopt the material parameters used by Mellon (1997) for Martian permafrost. These are for a linear mixture of ice and rock with 45% ice volume, and are shown in Table 1. Comparing the values in Table 1 to those in the discussion in Sect. 2, it can be seen that our Young’s modulus value (with range ~ 9– 20 GPa) is higher than the upper limits suggested by Philae. The measured material properties of the comet’s surface are not well constrained,however, so we choose our frozen soil parameters partly for ease of comparison with previous studies. We are also justified by the fact that, as noted in Mellon (1997) and references therein, ice, snow (i.e. low-density granular ice) and ice-bonded soils have similar rheologies. Any hard layer containing significant amounts of ice should therefore exhibit the same type of viscoelastic behaviour. The stratigraphy of the comet is also poorly constrained, so we use a simple one-layer model in order to quantify the maximum induced stresses. We vary each of the main parameters in Sect. 4.4 to investigate the sensitivity of the results to these assumptions. In terms of thermal inertia, the Rosetta results imply low bulk values of I = 10–50 J m⁻² K⁻¹ s^−1∕2 (Gulkis et al. 2015), but a hard ice-enhanced layer may have higher thermal inertia, closer to that of pure water-ice (2000 J m⁻² K⁻¹ s^−1∕2; Ferrari & Lucas 2016). We therefore run our model with values of I = 10, 50, 200, and 1000 J m⁻² K⁻¹ s^−1∕2.

Fig. 1 Reproduction of the results of Fig. 3 from Mellon (1997). Top panel: temperature over a Martian year for the Viking 2 landing site at a depth of 0.5 m. The solid curve is a digitisation of the Mellon (1997) data with the grey shading representing a ± 5 K uncertainty. Bottom panel: induced thermal stresses from Mellon (1997; solid) and from our code (dashed), including the uncertainty from the temperature digitalisation. RMS residuals between the Mellon (1997) results and ours are 746 kPa.

4 Modelling results

4.1 Temperature maps

The output of the spherical nucleus thermal model is shown for I = 50 J m⁻² K⁻¹ s^−1∕2 in Fig. 2 as maps of temperature with time and depth for each of the latitudes examined. The seasonal thermal wave can be seen propagating downwards from the surface to some depth, below which temperature variations are negligible and the material remains near the initial temperature of 30 K. The phasing of the seasonal change varies between the northern and southern hemispheres because of the comet’s obliquity. The spin axis is aligned such that at aphelion the northern pole points roughly towards the sun and the southern pole away from it, and vice versa at perihelion. Therefore, as the comet approaches the Sun and insolation increases, northern latitudes become progressively warmer before experiencing a sudden temperature drop as they are orientated away from the Sun (winter), whereas southern latitudes receive very little sunshine until a sudden burst at closest approach, leading to a sharp temperature peak (summer). Because of this, surface temperatures in the southern hemisphere reach greater maximum values of just over 200 K, while peak temperatures in the north are roughly 190 K, and less near the equator. Despite these differences, the depth to which the seasonal thermal wave penetrates is similar at all latitudes: ~ 100 ξ, or 0.5–1 m for I = 50 J m⁻² K⁻¹ s^−1∕2. Similar results are also seen when varying I, but with correspondingly different magnitudes of depth, as discussed below.

Fig. 2 Temperature maps with depth and time for each of the labelled latitudes on comet 67P, produced by our 1D thermal model with a thermal inertia of I = 50 J m−2 K−1 s−1∕2. Zero/one years is aphelion and 0.5 yr is perihelion.

4.2 Stress model results

Figure 3 shows maps of the thermal stresses induced by the temperatures of Fig. 2. The general trend is for compressive stresses to follow the increasing progression of temperature in summer, followed by tensile stresses as the material cools. Variations of magnitude and timing can be seen with latitude. Tensile stresses remain high (up to tens of MPa) for much of the cycle, particularly in the cold southern hemisphere, where the depth of the transition from tension to compression is roughly the same in all cases at ~ 0.25 m. In the north, this depth is greater at the pole, around 0.45 m, and decreases towards the equator to only 0.1 m or less. Below around 0.5 m the model failed to converge, due to a combination of small temperature variation and high stiffness (E), and we exclude these depths from the plot.

The highlighted contours show the 1, 2.5, and 10 MPa values for the tensile strength of an ice-bonded hard layer from Kossacki et al. (2015). Therefore in our baseline case, thermal stresses are high enough to exceed the tensile strength of water-ice and ice-rock mixtures, and cause fracturing to depths of tens of centimetres over most of the cometary surface. In low northern latitudes the situation is less clear cut and fracturing may occur to at least between 1 and 10 cm for weak materials (with tensile strengths below 1 MPa).

Figures 4–6show the results for thermal inertia values of I = 10, 200, and 1000 J m⁻² K⁻¹ s^−1∕2, respectively.The shapes of the stress patterns are similar in each case, but shifted up or down in depth by the differing diurnal skin depths. For I = 10 J m⁻² K⁻¹ s^−1∕2, tensile stresses are only high enough to overcome the likely tensile strengths down to a few centimetres,and not at all near the equator, but for I = 200 J m⁻² K⁻¹ s^−1∕2 fracturing can be expected down to at least a metre in most cases and to at least several metres for I = 1000 J m⁻² K⁻¹ s^−1∕2.

These results are summarised in Fig. 7, which shows the maximum seasonal stress as a function of latitude for the four thermal inertia values. The highlighted contours again show values of tensile strength from Kossacki et al. (2015). The top panel shows the maximum stress encountered at the surface for the four values of thermal inertia, plus an additional curve computed using no thermal inertia, I = 0 J m⁻² K⁻¹ s^−1∕2. The band of low stress centred just north of the equator is clearly seen in all cases, with a width of between around 80^∘ and 20^∘ depending on I. This represents the region with the narrowest seasonal temperature range (as shown in Fig. 2) due to the particular orientation of the comet’s obliquity. A trend for increasing stress with thermal inertia can also be seen near the equator; this results from the higher average diurnal temperatures produced by higher thermal inertia during the summer. The results of the zero thermal inertia model show stresses that are everywhere equal to or lower than the I = 50–1000 J m⁻² K⁻¹ s^−1∕2 curves, and generally lower than the I = 10 J m⁻² K⁻¹ s^−1∕2 curve, except at the equator itself. A model without thermal inertia can therefore be considered as a lower estimate of thermal stress.

Fig. 3 Thermal stresses with depth and time for each of the labelled latitudes on comet 67P, with I = 50 J m−2 K−1 s−1∕2. Highlighted contours are likely tensile strengths, from Kossacki et al. (2015), of 1, 2.5, and 10 MPa. Zero/one years is aphelion and 0.5 yr is perihelion.

Fig. 4 Thermal stresses. Identical to Fig. 3, but with a thermal inertia of I = 10 J m−2 K−1 s−1∕2.

Fig. 5 Thermal stresses. Identical to Fig. 3, but with a thermal inertia of I = 200 J m−2 K−1 s−1∕2.

Fig. 6 Thermal stresses. Identical to Fig. 3, but with a thermal inertia of I = 1000 J m−2 K−1 s−1∕2.

Fig. 7 Top: maximum surface thermal stress over one revolution for different thermal inertia values. Bottom: maximum thermal stress maps with depth and latitude over one revolution for the four non-zero thermal inertia values.

4.3 Stress results on the shape model

We compare these results to the simplified thermal stress model using the comet shape, and to the observed locations of polygons. Figure 8 shows the maximum seasonal stress found at the surface using this approach. Thermal stresses broadly follow the pattern of Fig. 7, with the greatest stresses near the south pole, intermediate stress near the north pole, and the lowest values near the equator and in the low-latitude northern regions. Significant “patchiness” is seen because of the effects of local topography and the low resolution of the shape model used. Absolute stress values are not shown here due to the lack of thermal inertia making them unrealistic; the figure is mainly presented to compare stresses in different regions in a qualitative manner. As described in the previous section, however, this zero thermal inertia model can be considered as a minimum stress estimate compared to models with realistic thermal inertia values.

Also shown in Fig. 8 are the locations of the polygons measured by Auger et al. (2018). They are found all over the comet, at varying latitudes and on facets with low and high stress, with no apparent correlation to peak thermal stress. No significant differences are noted between the distribution of stress on all facets and those containing polygons. As noted in Auger et al. (2018) there is a bias against the detection of polygons in the southern hemisphere, and therefore the highest stress regions, due to incomplete image coverage here at the time of the survey. The presence of polygons everywhere suggests stresses exceeding the material strength at all latitudes and, since this only occurs for the higher thermal inertia models described above, we infer I ≳ 50 J m⁻² K⁻¹ s^−1∕2 for all the regions where polygons are observed.

Fig. 8 Maximum thermal stresses on the surface of the 67P shape model for the I = 0 J m−2 K−1 s−1∕2 model. Values are for qualitative comparison only. Black points are the locations of the thermal contraction crack polygons measured in Auger et al. (2018). No significant differences are noted between the distribution of stresses on all facets and those containing polygons.

4.4 Sensitivity to material parameters

We now examine the sensitivity of these results to changes in the material parameters, many of which are poorly constrained.

In Fig. 9 we examine the effects of changing the ice content of the modelled ice-rock mixture from our baseline 45% ice volume, using the equatorial I = 50 J m⁻² K⁻¹ s^−1∕2 case as an example. This is achieved by modifying the thermal model as well as the thermal expansion parameter, α, and the viscous term, as shown in Table 1 and described in Mellon (1997) and Mellon et al. (2008). A trend for decreasing stress with ice content is seen, with no fracturing predicted to occur for the 25% ice volume case. For even lower ice content, thermal expansion is small and the model begins to break down and experience convergence problems at smaller and smaller depths. Figure 10 shows the same test for a latitude of − 60^∘. The same trend is observed, but here the stresses are high enough that fracturing should still occur even with small amounts of ice.

The trends in these plots are driven by the large thermal expansion coefficient of ice. Holding the other parameters constant and varying α on its own produces a similar effect, a roughly linear relation between peak stress and thermal expansion, while varying the viscous constant, A₀ , on its own has little effect. Similarly, the viscous exponent, n, appears to have little effect when varied between the different values of Mellon (1997) and Mellon et al. (2008). This is consistent with the results of Tauber & Kuehrt (1987), who suggested viscous effects to be of lesser importance than elastic ones. Varying ice content will, additionally, change the Young’s modulus and tensile strength. As shown in the contours in Figs. 3–7 above, tensile strength is two to three times higher for ice-rock mixtures than for pure water ice (Kossacki et al. 2015). Mellon (1997) describe the baseline Young’s modulus used here as conservative, and note that coarse-grained mixtures are found to have values up to five times higher than pure ice, whereas Philae’s measurements of the hard layer suggest a value an order of magnitude lower. Holding the other parameters constant and varying E by an order of magnitude in either direction produced roughly corresponding changes in peak stress, suggesting that the two effects are approximately balanced (i.e. greater rock content will decrease the Young’s modulus, and therefore the induced stresses, but also reduce tensile strength by the same order of magnitude),leaving the thermal expansion coefficient to once again dominate.

From this we can see that the primary influences on thermal stress are thermal inertia, ice content, and the material strength and stiffness (tensile strength and Young’s modulus) of the hard layer. If thermal inertia varies across the comet then the depths to which high stresses penetrate will be different in different regions. This would affect the depths of crack propagation and distance between adjacent fractures. The other mechanical parameters are interconnected, so that a strong, stiff material may be more resistant to thermal expansion, while a weaker material (lower E) will experience lower stress, but will also require lower stresses to fracture (lower tensile strength). Unfortunately, the properties of partially bonded granular materials and the relations between these parameters are poorly understood. However, we can say that a hard layer with thermal inertia ≳ 50 J m⁻² K⁻¹ s^−1∕2, and with significant quantities of water ice (≳ 45% at the equator, less elsewhere) should fracture, whatever its locations on the comet. If no such layer exists, and the material is purely unconsolidated and granular, then it is difficult to see how fracturing can be supported with our model.

Fig. 9 Thermal stress for the equator of 67P with I = 50 J m−2 K−1 s−1∕2 and three different values of ice volume fraction.

Fig. 10 Thermal stress for a latitude of − 60∘ with I = 50 J m−2 K−1 s−1∕2 and three different values of ice volume fraction.

5 Discussion and conclusions

5.1 Summary

The results of our thermo-viscoelastic model show stresses of up to several tens of MPa induced by seasonal temperature changes on a spherical model of comet 67P down to depths of between a centimetre and a metre. These are seen at most locations on the comet, with the exception of a band of low stress in the middle and low latitudes of the northern hemisphere (Fig. 7), which is particularly pronounced for low values of thermal inertia (I ≲ 10 J m⁻² K⁻¹ s^−1∕2). Introducing local topography, in the form of 67P’s complicated shape (but neglecting thermal inertia), limits these regions of low surface stress to patchy areas, mostly still in the north and equatorial region (Fig. 8). On the spherical model, stress increases with increasing thermal inertia and, for values of I = 50–1000 J m⁻² K⁻¹ s^−1∕2 (and ice volume fractions ≳ 45%), exceeds reasonable material tensile strengths of ~1 MPa at all latitudes, leading to fracturing.

Additionally, stresses vary with ice content, Young’s modulus, and thermal expansion coefficient. If these material parameters, as well as thermal inertia, vary across the surface then fracturing will not be uniform. Stress decreases with decreasing ice content and at zero ice-fraction the unconsolidated material is probably not competent enough to support fractures. Between this and the lower end of our simulations (25% ice) the minimum ice content for fracturing is hard to quantify, and varies with thermal inertia and latitude (see Figs. 9 and 10). The strongest constraints are found in the equatorial region, where 45% ice is required. For comparison, Mellon (1997) note that soils with low ice content behave similarly to dry soil (and ice with small impurities retains ice-like rheology), while thermal contraction polygons have been observed in terrestrial surfaces ranging from nearly pure ice (Marchant et al. 2002) to permafrost with 10% or less ice content (Bockheim 1995). Various Rosetta measurements (Fulle et al. 2016) point to a global dust-to-ice ratio for 67P’s nucleus of around 6, or an average of 15% water-ice, insufficient to cause fracturing at all latitudes in our model (though probably sufficient at the poles). A higher than average ice content may well be expected, however, in the hardened layer of interest. Oklay et al. (2017), for example, saw ice-enhanced regions on the surface and boulders with ~48% and ~25% ice content, respectively, suggesting that enhancements in an ice-rich subsurface layer to the volume fractions examined here are not unreasonable.

To summarise, a hard layer, at a depth of a few centimetres to metres and with sufficient thermal inertia (I ≳ 50 Jm⁻² K⁻¹ s^−1∕2) and ice content, will undergo thermal fracturing globally. Near the poles this layer may have approximately the global average ice content, but near the equator it must be significantly more icy (≳45%) and/or have a higher thermal inertia for fracturing to occur.

5.2 Comparison with observed thermal contraction crack polygons

Thermal contraction polygons are found all over 67P’s surface (Fig. 8). There appears to be little correlation between peak stress, as predicted by our model, and polygon location.

Wherever polygons are observed, we predict thermal inertia and ice content values of I ≳ 50 J m⁻² K⁻¹ s^−1∕2 and ≳45% (for the equator), respectively, at the time of formation. For the polar regions, lower thermal inertia and ice content constraints are placed because of the higher temperature ranges and stresses. Higher inertias are expected on the consolidatedterrains than on the dust covered surfaces, and this is indeed where polygons are found. The presence of polygons on consolidated, high thermal inertia terrains is consistent with current seasonal thermal stresses, but we cannot also rule out ancient polygons resulting from differing thermal environments due to past orbital or obliquity changes.

In terms of polygon size, examples on Earth and Mars typically have a distance between adjacent cracks of three (Lachenbruch 1962) to ten (El Maarry et al. 2010) times their fracture depth. Fracture propagation is complex (Lachenbruch 1962) and not well studied in a vacuum (Molaro et al. 2015), but is enhanced in cold, brittle materials under low hydrostatic pressure, as occurs on the comet. Under such conditions, fractures can extend to depths 3– 15 times that of the high-stressed region (Maloof et al. 2002). From Fig. 7, therefore, thermal fractures of 67P penetrate to between 3 cm and 15 m deep, with expected polygons of between 0.09 and 45 m across, but likely ~ 1– 20 m for reasonable (I = 50 J m⁻² K⁻¹ s^−1∕2) thermal inertia values. This is consistent with the measured polygons of mean size 3 m and 90% between one and five metres (Auger et al. 2018). Observed polygon size is relatively uniform across the comet, suggesting little variation in the depth or mechanical properties of the hard layer. This is, however, slightly puzzling when compared to the variation in the depth of the stressed layer shown in Fig. 7. It is possible that variations in topography and thermal environment change the stress locally (as shown in Fig. 8) confusing any smooth variation of polygon size with latitude.

5.3 Implications for erosion

Thermal fracturing on the metre-scale should be an important erosion mechanism. Gradual weakening of material by thermal fatigue will break down boulders into smaller and smaller pieces (Pajola et al. 2015), whilst the presence of fractures in cliff walls will weaken the material, making collapse under its own weight more likely. Debris fields at the bottom of many cliffs suggest this is a common process (El-Maarry et al. 2015a; Groussin et al. 2015), while individual collapses have been linked with outburst activity (Pajola et al. 2017). Vincent et al. (2016) suggest that intact but fractured cliffs allow heat to penetrate further into the underlying volatile materials, increasing sublimation and driving activity. Our model supports and validates the idea of weakening material with thermal stresses.

We extend these ideas of eroding cliff activity by first considering the formation and fracturing of a hard layer. Figure 11 shows a cartoon of a cliff topped by a dust mantle. If the dust-layer is thick enough, then it will insulate the material immediately below it from large temperature changes, but the exposed consolidated material of the cliff wall will suffer the full seasonal temperature profile, as shown in Fig. 2. Volatiles in the outer layers will sublimate and begin to flow outwards through the cliff face as the comet heats up. During the comet’s retreat from the Sun, the cliff surface will begin to cool again which, in some cases, can lead to a temperature inversion where the surface layers are cooler than the interior. In Fig. 11 we show temperature profiles extracted from our thermal model for the south polar case, where clear temperature inversions can be seen shortly after perihelion. Sublimating water ice will recondense in these regions (the upper ~10 cm in this case) and may form a hard layer which will experience high stresses and thermal fracturing on the next, and subsequent, orbital cycles. Once fractured, the cliff is more vulnerable to partial collapse, leaving an overhang, full collapse into a debris field, and the sort of enhanced gas flow and activity proposed by Vincent et al. (2016). This process could be occurring all over the comet, leading to the gradual retreat of cliffs and removal of material across the surface.

It should also be noted that the same mechanisms will apply to other objects, similar to 67P. Other Jupiter-family comets undergo comparable seasonal temperature changes and, therefore, the presence of hard, ice layers and thermal fracturing isexpected on them as well. Thermal contraction polygons may well be found on other comets, if observed with sufficient resolution (Auger et al. 2018), and the linked processes of fracturing, cliff erosion and activity may also be a common feature (Vincent et al. 2016).

Fig. 11 Model for the collapse of cliffs driven by thermal stresses. Solar insolation heats exposed interior material on the face of cliffs, driving the sublimation of volatiles, whilst flat, dust covered surfaces are insulated from large temperature changes. Temperature inversions can occur as the seasonal thermal wave penetrates to depth, even as the surface begins to cool again, allowing recondensation of the volatiles in the top layers. Recondensed water-ice forms a hard layer which is vulnerable to thermal fracturing, leading eventually to cliff collapse. Temperature profiles are taken from the I = 50 J m−2 K−1 s−1∕2, latitude = −90∘ case, for times in days since perihelion. Grey boxes show the regions, within five skin depths of the surface, where diurnal temperature changes dominate.

5.4 Perspectives

In this workwe have considered two models: the simplified cases of a uniform material with depth on a spherical nucleus, and a material without thermal inertia on a shape model with complex topography. A more complete model could combine both of these to simulate temperature with depth, thermal inertia, and topography, but as noted this would be very computationally expensive.

We also did not investigate diurnal temperature changes, which are similar in magnitude to the seasonal changes, but on shorter timescales, leading to high rates of change. The limited penetration of the diurnal temperature wave (skin depth of 1– 2 cm; Gulkis et al. 2015), however, should limit the induced stresses to only the very top layers. Thus, diurnal stresses may well be important in the weathering of the outer centimetres of boulders and consolidated material, but are unlikely to play a role in the formation of polygons and other metre-scale fractures. Fracturing has indeed been observed at these small scales (Poulet et al. 2016) and may be evidence for such a scenario. The fracture networks observed by CIVA do not resemble the regular polygons found at larger scales. However, the stress pattern may be modified by pre-existing weaknesses in the micro-structure (which can affect polygon shape; Lachenbruch 1962; Auger et al. 2018) or by other stresses, such as from desiccation. Desiccation-formed fractures on Mars are similar in morphology to thermal contraction polygons (El Maarry et al. 2010), but the interactions of different stresses from seasonal and diurnal temperature changes, thermal fatigue, and thermal shock, and of desiccation, have not been well studied.

Finally, the work might be extended by following the finite element modelling of Mellon et al. (2008) to further investigate the relationship between stress, ice content, and polygon size. Nevertheless, uncertainties in material properties would involve a large parameter space.

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 686709. This work was supported by the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number 16.0008-2. The opinions expressed and arguments employed herein do not necessarily reflect the official view of the Swiss Government. We thank the reviewers for their useful comments.

References

All Tables

All Figures

Fig. 1 Reproduction of the results of Fig. 3 from Mellon (1997). Top panel: temperature over a Martian year for the Viking 2 landing site at a depth of 0.5 m. The solid curve is a digitisation of the Mellon (1997) data with the grey shading representing a ± 5 K uncertainty. Bottom panel: induced thermal stresses from Mellon (1997; solid) and from our code (dashed), including the uncertainty from the temperature digitalisation. RMS residuals between the Mellon (1997) results and ours are 746 kPa.

In the text

	Fig. 2 Temperature maps with depth and time for each of the labelled latitudes on comet 67P, produced by our 1D thermal model with a thermal inertia of I = 50 J m−2 K−1 s−1∕2. Zero/one years is aphelion and 0.5 yr is perihelion.
In the text

	Fig. 3 Thermal stresses with depth and time for each of the labelled latitudes on comet 67P, with I = 50 J m−2 K−1 s−1∕2. Highlighted contours are likely tensile strengths, from Kossacki et al. (2015), of 1, 2.5, and 10 MPa. Zero/one years is aphelion and 0.5 yr is perihelion.
In the text

	Fig. 4 Thermal stresses. Identical to Fig. 3, but with a thermal inertia of I = 10 J m−2 K−1 s−1∕2.
In the text

	Fig. 5 Thermal stresses. Identical to Fig. 3, but with a thermal inertia of I = 200 J m−2 K−1 s−1∕2.
In the text

	Fig. 6 Thermal stresses. Identical to Fig. 3, but with a thermal inertia of I = 1000 J m−2 K−1 s−1∕2.
In the text

	Fig. 7 Top: maximum surface thermal stress over one revolution for different thermal inertia values. Bottom: maximum thermal stress maps with depth and latitude over one revolution for the four non-zero thermal inertia values.
In the text

	Fig. 8 Maximum thermal stresses on the surface of the 67P shape model for the I = 0 J m−2 K−1 s−1∕2 model. Values are for qualitative comparison only. Black points are the locations of the thermal contraction crack polygons measured in Auger et al. (2018). No significant differences are noted between the distribution of stresses on all facets and those containing polygons.
In the text

	Fig. 9 Thermal stress for the equator of 67P with I = 50 J m−2 K−1 s−1∕2 and three different values of ice volume fraction.
In the text

	Fig. 10 Thermal stress for a latitude of − 60∘ with I = 50 J m−2 K−1 s−1∕2 and three different values of ice volume fraction.
In the text

Fig. 11 Model for the collapse of cliffs driven by thermal stresses. Solar insolation heats exposed interior material on the face of cliffs, driving the sublimation of volatiles, whilst flat, dust covered surfaces are insulated from large temperature changes. Temperature inversions can occur as the seasonal thermal wave penetrates to depth, even as the surface begins to cool again, allowing recondensation of the volatiles in the top layers. Recondensed water-ice forms a hard layer which is vulnerable to thermal fracturing, leading eventually to cliff collapse. Temperature profiles are taken from the I = 50 J m−2 K−1 s−1∕2, latitude = −90∘ case, for times in days since perihelion. Grey boxes show the regions, within five skin depths of the surface, where diurnal temperature changes dominate.

In the text