  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 ($\gtrsim50$ 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 $\sim0.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.
