We simulate the stresses induced by temperature changes in a putative hard layer near the surface of comet 67P/Churyumov–Gerasimenko with a thermo-viscoelastic model . Such a layer could be formed by the recondensation or sintering of water ice ( and dust grains ) , as suggested by laboratory experiments and computer simulations , and would explain the high compressive strength encountered by experiments on board the Philae lander . Changes in temperature from seasonal insolation variation penetrate into the comet ’ s surface to depths controlled by the thermal inertia , causing the material to expand and contract . Modelling this with a Maxwellian viscoelastic response on a spherical nucleus , we show that a hard , icy layer with similar properties to Martian permafrost will experience high stresses : up to tens of MPa , which exceed its material strength ( a few MPa ) , down to depths of centimetres to a metre . The stress distribution with latitude is confirmed qualitatively when taking into account the comet ’ s complex shape but neglecting thermal inertia . Stress is found to be comparable to the material strength everywhere for sufficient thermal inertia ( \gtrsim 50 J m ^ { -2 } K ^ { -1 } s ^ { -1 / 2 } ) and ice content ( \gtrsim 45 \% at the equator ) . In this case , stresses penetrate to a typical depth of \sim 0.25 m , consistent with the detection of metre-scale thermal contraction crack polygons all over the comet . Thermal fracturing may be an important erosion process on cometary surfaces which breaks down material and weakens cliffs .