Tensile Strength of 67P/Churyumov-Gerasimenko Nucleus Material from Overhangs

We directly measure twenty overhanging cliffs on the surface of comet 67P/Churyumov-Gerasimenko extracted from the latest shape model and estimate the minimum tensile strengths needed to support them against collapse under the comet's gravity. We find extremely low strengths of around one Pa or less (one to five Pa, when scaled to a metre length). The presence of eroded material at the base of most overhangs, as well as the observed collapse of two features and implied previous collapse of another, suggests that they are prone to failure and that true material strengths are close to these lower limits (although we only consider static stresses and not dynamic stress from, for example, cometary activity). Thus, a tensile strength of a few pascals is a good approximation for the tensile strength of 67P's nucleus material, which is in agreement with previous work. We find no particular trends in overhang properties with size, over the $\sim10-100$ m range studied here, or location on the nucleus. There are no obvious differences, in terms of strength, height or evidence of collapse, between the populations of overhangs on the two cometary lobes, suggesting that 67P is relatively homogenous in terms of tensile strength. Low material strengths are supportive of cometary formation as a primordial rubble pile or by collisional fragmentation of a small (tens of km) body.


Introduction
Material strength is an important parameter in constraining the formation and evolution of comets and in explaining their morphological diversity.Low strengths support the conclusion that comets are primordial rubble piles, accreted gently (collision velocities of ∼ 1 m s −1 to tens of m s −1 ) in the early solar system, as opposed to remnants from higher velocity collisions (at hundreds of m s −1 to km s −1 ) which would undergo impact compaction, leaving them with higher strengths (Davidsson et al. 2016).The strength of cliffs and overhangs against collapse also directly relates to cometary activity as cliff collapses are an important source of jets and outbursts (Vincent et al. 2016;Pajola et al. 2017).
The tensile strength of cometary nucleus material has been estimated for a number of comets and using a number of different methods including: observations of comet breakups, from rotation (see for example Davidsson 2001) or close Solar (Klinger et al. 1989;Steckloff et al. 2014) or Jupiter (Asphaug & Benz 1996) encounters; laboratory experiments (Blum et al. 2006(Blum et al. , 2014;;Bar-Nun et al. 2007); computer modelling (Greenberg et al. 1995;Biele et al. 2009); and, in the case of 67P/Churyumov-Gerasimenko (hereafter 67P), observations of cliffs, overhangs and boul-ders.As summarised by Groussin et al. (2015), strength estimates vary over several orders of magnitude for cometary material but are generally low: in the range of pascals to tens of pascals for tensile strength (σ T ) at the metre scale, and larger for shear and compressive strengths.The discussion of scale is important because consolidated material generally shows decreasing strength at larger scales, following a power law proportional to d −q , with length scale d and an exponent of q ∼ 0.6 for water ice (Petrovic 2003).
For 67P, the high quality of OSIRIS (Optical, Spectroscopic, and Infrared Remote Imaging System) imaging allowed the examination of a variety of features at the metre and decametre scale and, from this, Groussin et al. (2015) measured the tensile strength of an overhang and a collapsed feature, estimating σ T = 1 − 3 Pa and < 150 Pa at the 5−30 m scale.When scaled to 1 m, this gives σ T = 8−39 Pa and < 1150 Pa for these features, which are located in the Imhotep and Maftet regions, on the head and body of the comet, respectively (see Thomas et al. 2015 andEl-Maarry et al. 2015 for a detailed description of 67P's regions).Vincent et al. (2017) also measured the heights of cliffs and derived 1 − 2 Pa strengths at the decametre scale across 67P's surface in order for the cliffs not to collapse.
The goal of this work is to perform a more comprehensive survey of overhanging cliffs across comet 67P's surface,

Method
We measure overhang properties using a digital terrain model (DTM or shape model) of comet 67P, constructed using the Stereo Photogrammetry (SPG) technique from Rosetta/OSIRIS images.Groussin et al. (2015) showed that the alternative SPC technique can underestimate the number of high-slope facets compared to SPG, making SPG a more appropriate technique for studying sharp topography.We use the latest DLR shape model (SPG-SHAP7, Preusker et al. 2017) with 22 million vertices, with a typical spacing of 1-2 m and orientation uncertainty of 2-5 • .For computational reasons, we use a decimated version (2 million facets, ∼ 7 m lateral vertex spacing) of the full shape model.Gravitational vectors are calculated for each facet following the method of Jorda et al. (2012), which includes centrifugal forces, and are then compared to the facet normal vectors to produce a map of gravitational slope.Locally flat regions have a gravitational slope of 0 • , while vertical cliffs are at 90 • and anything larger is an overhang.There are 15, 358 facets (0.77% of the total number) with slopes greater than 90 • and approximately 5000 greater than 100 • .A number of these are artefacts of the reconstruction but many are real overhanging features, visible in OSIRIS images.
Due to the complicated shape of the nucleus, and the presence of artefacts in the reconstructed shape model, it was deemed impossible to automatically detect and characterise every overhanging feature.Instead we identify, by eye, a number of features and investigate them in detail.The features are selected by picking the most obvious large groupings of overhanging facets on the shape model, from all over the nucleus.Selecting the largest (actually 'deepest') overhangs places the strongest constraints on strength.We examine the spatial distribution of selected features to ensure good coverage in Section 3.
For each selected feature, a local DTM is extracted from the full model, and a vertical profile of the overhang is computed in the following way.First, the negatives of the gravitational vectors are plotted on each high-slope facet.These point 'up', ie.away from the centre of gravity, and, in an overhang region, will penetrate inside the shape model before exiting it again at the top of the overhang (see Figure 1).The two facets where the gravitational vector enters and exits the nucleus define the coordinates of the 'base' and 'top' of the overhang.Selecting a third point as a facet on the 'face' of the overhang, perpendicular to the gravity vector, defines a plane, the intersection of which with the shape model can be calculated.The final step is to project the coordinates of the intersection points onto the plane itself and then rotate the data to align them with the local gravity vector.This vector will vary across the overhang but the differences over such short distances (∼ 100 m) are negligible, and the 'base' gravity vector is used for the whole profile here.A profile of the shape model along the chosen vector is then recovered.
From such a profile, an estimate of the strength of the overhang can be calculated using the equations of a can- tilevered beam.Following Tokashiki & Aydan (2010) and using the coordinate frame centred on the overhang base (as shown in Fig. 1), the maximum stress from bending is found at x = 0 and decreases along the length.The material tensile strength must be at least as large as this maximum stress to prevent immediate failure and collapse.For a beam of unit thickness and length L and height h, this can be expressed as where M = L 0 xρgh dx is the total bending moment acting on the cantilever from its own unit weight of: unit density, ρ = 537.8kg m −3 (Preusker et al. 2017), times the magnitude of the local gravity, g.For simple beam shapes, this can be integrated directly.In a rectangular beam, for example, h is constant and the bending moment evaluates to M = hρgL 2 /2 while stress becomes For a trapezium shaped beam, as shown in Fig. 1, h varies along the length but the integral can still be evaluated, so that this time the stress becomes where α = h s /h b is the ratio of the height at the far-end to the base height.These two equations are useful approximations for many overhangs and were used in the previous works (Groussin et al. 2015;Tokashiki & Aydan 2010).In our case, however, we have the full overhang profile from the shape model and can therefore integrate this directly, without having to use one of these approximations.We do this numerically by, first, interpolating the profile shape to a regularly spaced series of x positions (separated by dx), measuring h at each of these and then making the sum As can be seen from Eqs. 2 and 3, the overhang length is the most important parameter for constraining the material strength, and this remains true for the numerical integration.Therefore, for each of our overhangs we select the intersection plane to compute the profile along the deepest part of the overhang, by choosing the 'base' facet with a gravity vector penetrating most deeply into the shape model (maximum L).In most cases, this is obvious from visual inspection and in cases with several similar facets, similar strength estimates will be derived.
There is an uncertainty in the position of all coordinates measured on the shape mode, which can be conservatively estimated as the average radius of a facet.For our two million facet model this is ∼ 1.6 m.In addition, the uncertainty in the orientation of each facet has been estimated as ∆θ ≈ ±5 • (Jorda et al. 2012) and in practice this can dominate over the positional uncertainty.We therefore estimate the uncertainties in overhang proportions by rotating each profile by ±∆θ, computing the integral and measuring h, and then combining this uncertainty with the position error.Uncertainties in tensile strength are then derived using the standard error propagation formulae.

Results
Table 3 shows the results for the 20 overhanging features analysed.Overhang heights are between ∼ 10 and 100 m and depths are generally only ∼ 10 m (Fig. ??), with a single feature (number 2) having a larger, almost 40 m, depth.This feature, part of a cliff in the Babi region facing Hapi, is visually the largest overhang on the shape model but is otherwise unremarkable.
The calculated overhang strengths are very low, with a mean of 0.3 Pa and an average uncertainty of ±0.22 Pa. Figure 2 shows the distribution in measured strengths and those scaled from the feature length-scale, h, to 1 m, using the power law for ice (the scaling law for ice is used, despite the large dust content, as a first estimate and for ease of comparison with previous studies).The scaled strengths all lie between zero and 5 Pa, apart from the single large outlier of feature 2. Nonetheless, the uncertainty in sigma for this feature could easily bring its value in line with the others.We detect no particular relation between strength, scaled or raw, and overhang height, suggesting a uniform strength over the range of ∼ 10 − 100 m.
Figures 3, 4 and 5 show the locations of each measured feature on the cometary surface, along with all the highslope facets.High-slopes are typically clustered into curvilinear chains of cliffs.The only large area lacking highslopes (other than artefacts) is the smooth, dusty terrain of Hapi.Our selected overhangs are well distributed across the surface, with features seen in both hemispheres, both cometary lobes and in a variety of cometary regions.No particular trends are noted between overhang properties and location, and there are no significant differences between the average strengths on each lobe, suggesting a uniformity in the capacity of 67P's material to support overhangs.
The features have a variety of morphologies, as can be seen in their profiles and accompanying OSIRIS images in Figs. 6 -9.They range from large, shallow cliffs to protruding blocks and the lips of pit walls.Most have a shallow trapezium or triangular shape, with a number having curving top surfaces, and are not obvious 'caves', but rather are steep cliffs which lean slightly 'outwards' towards their tops, producing an overhanging region.Several features are actually double or triple overhangs (features 1, 4, 16 and 19 in Figs. 6 -9), and for these we split the integration region before calculating the moment of each area and manually summing them together.Nearly all the overhangs examined here show some evidence of collapse or erosion, such as boulders and debris fields at their base.The exceptions are features 1, 4, and 17, where the possible debris are some distance away due to their positions at the top of tall cliffs, and 3 and 20, where their presence is slightly ambiguous.
Of particular interest are features 5, 12 and 15.Feature 12 is the same as the failed section already measured in Groussin et al. (2015).The overhang in Ash near feature 5 was seen to suffer a collapse between May and December 2015 (El-Maarry et al. 2017b supplementary Fig. 2).We examine a local high-resolution DTM of the area (Preusker, F., personal communication) but were unfortunately unable to produce a profile of the collapsed region due to technical issues.Measuring directly from the DTM in a 3D model viewer, however, we estimate the overhanging segment to have L ≈ 12 m and h ≈ 31 m, giving σ T ≈ 0.4 Pa with the rectangular approximation of Eqn. 2, in-line with the other overhangs measured here.Similarly, feature 15, located in Seth and named Aswan, was seen to collapse in July 2015 (triggering an outburst, Pajola et al. 2017).They measure the overhanging section as 65 × 12 m, again consistent with the 63 × 11 m, σ = 0.33 Pa found here.

Are overhangs representative of bulk nucleus strength?
The strengths estimated above are lower limits; the material must have at least this tensile strength, in order not to immediately collapse under its own weight, but could be much stronger.The presence of collapses, however, implies that overhangs may be close to failure and that these estimates may indeed be close to the actual material strength.
Conversely, Vincent et al. (2017) and Pajola et al. (2017) argue that cliff heights are not controlled by intrinsic material strength but by external erosion, in which case observed features would not be limited by gravity and tensile strengths could be larger.Many of the features examined here have clearly fractured cliff walls (e.g.2,6,9,10,13,15,16,18), some reminiscent of thermal contraction crack No. Latitude ( polygons (Auger et al. 2015(Auger et al. , 2017)), which may imply additional processes, such as thermal stresses and sublimation, weakening them before failure occurs.Additionally, material strength can vary on a local scale.For example, results from the MUPUS and SESAME experiments on the Philae lander (Spohn et al. 2015;Knapmeyer et al. 2017) as well as theoretical (Kossacki et al. 2015) and laboratory work (Gruen et al. 1993;Kochan et al. 1989) suggest a hard, ice bonded layer with much greater strength within ∼ metres of the surface.The SESAME results in particular suggest a layer at depths of 10 − 50 cm with a tensile strength on the order of MPa (Knapmeyer et al. 2017).
In order to consider these possibilities, we compare the dimensions of the measured overhangs with the depths of a hard layer, and to those of temperatures relevant for thermal processing, using a comet thermal model described in Attree et al. (submitted).Our model takes into account a spherical nucleus (orientated according to its pole with RA = 69.57deg, DEC = 64.01deg (J2000) and with a rotational period of P = 12.40 hr; Jorda et al. 2016), solar insulation, and heat conductivity (Groussin & Lamy 2003).We compute 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 use a time step of 12.4 s and ran the thermal model over 5 complete revolutions.We then compute the maximum temperature over an orbital period, experienced at each depth interval (for 2000 depth intervals of a thickness of a fifth of a diurnal skin depth each) for three different values of thermal inertia: I = 10, 50 and 250 J m −2 K −1 s −1/2 .
Figure 10 shows the resulting thermal profiles at the equator.The horizontal lines show conservative estimates for the temperatures at which sublimation is negligible com-pared to the erosion rate at perihelion (one thousandth its value) for water ice (160 K) and CO 2 (85 K).Depths of 2.8 m, for water ice, and 8.4 m, for CO 2 , therefore define the limits to which we would expect sublimation to affect material properties.Our overhangs are typically ∼ 10 m deep and, with one exception (feature 9), all have L > 5 m.Assuming H 2 O is the dominant volatile component, we therefore expect our strength measurements to be probing material which has not been thermally processed.Further, since pressure inside the comet is on the order of tens to hundreds of Pa (Groussin et al. 2015), thermal or compressional processing should not affect material below these depths and we expect the measured tensile strength to be a good approximation for that of the bulk nucleus.
A hard, ice bonded layer will increase the average strength of near surface material.The fact that we measure such low strengths, however, implies that the hard layer is localised and does not contribute significantly to the average strength of cometary material at depth.This is consistent with an estimated hard layer thickness of 0.1 − 0.5 m, compared to the ∼ 5 − 10 m deep overhangs.Therefore, we conclude this section by reinforcing the idea that the tensile strengths measured above (∼ 0−5 Pa when scaled to metre lengths) are indeed representative of bulk cometary tensile strengths at the decametre scale.

Dynamic stresses
Thus far our analysis has only focused on the strength of overhangs to resist the static stress of their own weight, however dynamic stresses should also be present on the comet.These might include stresses induced by rotational changes (Hviid et al. 2016;Hirabayashi et al. 2016), as well as seismic events, such as cometary activity/outbursts and impacts (e.g.Thomas & Robinson 2005).
Impact features are rare on 67P's surface, and with only a few, small (∼ 10s of metres) craters detected during the duration of the Rosetta mission (El-Maarry et al. 2017a), seem unlikely to provide a comet-wide mechanism for overhang collapse.Activity related stresses, on the other hand, are very likely to have occurred during the comet's approach of the Sun, but quantifying them remains difficult due to the still poorly understood activity mechanisms.Seismic energy, released from a particular source of cometary activity, should fall off with the square of the distance times some attenuation factor (Thomas & Robinson 2005), which may be large due to the comet's fractured and porous nature.As a first order estimate, the stress in a weak elastic wave is σ = Uvρ (Melosh 1989) which, with a wave velocity of U > 80 ms −1 (as measured by experiments on Philae; Knapmeyer et al. 2017) and assuming a particle velocity of v = 1 cm s −1 , will be at least 430 Pa, greater than the strength estimates here.Continuous activity is not localised to particular regions, however, and even transient outburst features are found all-over the nucleus, suggesting it will hard to correlate the locations of activity with collapsed overhangs (the exception being where collapses trigger outbursts, such as in the Aswan case above).
Due to the reaction force on the nucleus, activity can alter cometary orbit and rotation, and several studies have investigated the stresses induced by such changes in rotation pattern on 67P's complex shape.Both Hviid et al. (2016) and Hirabayashi et al. (2016) found large stresses of up to several hundred Pa centred on the neck region.These dynamic stresses indeed exceed static stress and would lead to failure if the material has the strengths measured here.Our measured overhangs (both those with and without collapse features) have no apparent correlation with the neck, however.Additionally, the presence of collapses in features 5 and 12, located on the big lobe, far form the neck, argues for static-stress induced failure.A full analysis of the correlation between activity-and rotation-driven stress patterns and the location of overhangs could put further constraints on the material strength but is beyond the scope of this paper.Here we must limit our conclusion to the following: that static stress analysis provides lower limits to cometary material strength, but which may indeed be close to the real values due to the presence of overhangs which appear to have collapsed under their own weight.

Conclusion
We examine 20 overhanging cliffs, measuring their vertical profiles using the up-to-date SPG SHAP7 shape model (Preusker et al. 2017) of comet 67P.From this we derive lower limits for the material's tensile strength, in order to support such overhangs against gravity.Overhangs are generally shallow (most have depths ∼ 10 m) and so the resulting tensile strengths are very small; σ T ∼ 1 Pa or less at the decametre scale and ∼ 0 − 5 Pa when scaled to metre lengths (except for one outlier at ∼ 28 Pa, but with relatively large uncertainties).Nevertheless, the presence of eroded material at the base of most overhangs, and the observed collapse of two features and implied previous collapse of another, suggests that they are near to failure.Thus, σ T of a few pascals is a good estimate for the tensile strength of 67P's nucleus material, although further analysis of dynamic stresses, such as those caused by cometary activity and rotation changes, is warranted.Thermal modelling shows little material alteration at relevant depths in the subsurface, suggesting this value to be a reasonable approximation for bulk strengths at depth.This is in good agreement with previous estimates (as can be seen in the summary table in Groussin et al. 2015) from modelling (Greenberg et al. 1995;Biele et al. 2009), laboratory experiments (Bar-Nun et al. 2007;Blum et al. 2014) and some observations (Asphaug & Benz 1996), including cliff heights (Vincent et al. 2017).Other observations, such as from the breakup of sungrazing and rotating comets (Klinger et al. 1989;Davidsson 2001;Steckloff et al. 2014), suggest somewhat higher values of tens of Pa to ∼ 100 Pa.The Groussin et al. (2015) overhang results are slightly higher than those presented here because of their approximation of overhang shapes as rectangular, compared to the shape model profiles used here.
We find no particular trends in overhang properties with size, over the ∼ 10 − 100 m range studied here, or location on the nucleus.There are no obvious differences, in terms of strength, height or evidence of collapse, between the populations of overhangs on the two cometary lobes, suggesting that 67P is relatively homogenous in terms of tensile strength.
Such a low and homogeneous strength has implications for the formation and evolution of 67P.In terms of evolution, low strengths mean that material is easily eroded by sublimation, gas pressure and thermal fracturing, and is vulnerable to collapse under its own gravity.Collapses naturally explain the retreating cliffs, with debris fields and fallen boulders at their feet, seen across the comet (see for example Pajola et al. 2015;El-Maarry et al. 2017b;Pajola et al. 2017), as well as the presence of the overhangs themselves.These may form from the partial collapse of sections of cliff following preexisting weaknesses or which are further weakened by thermal fracturing (see model of erosion of Attree et al. submitted).Cometary outgassing activity has been linked to such collapses (Pajola et al. 2017) and to active cliff faces in general (Vincent et al. 2016), demonstrating that cliffs and overhangs are important areas of erosion on cometary surfaces.
Low bulk strengths support the conclusions of Davidsson et al. (2016) that 67P represents a primordial rubble pile, directly accreted from the proto-solar nebular by hierarchical aggregation or streaming instabilities, or is a colli- sional fragment from a small (tens of km) body.A fragment, or rubble pile of fragments, from the disruption of a larger (∼ 1000 km) body would inherit some of that body's properties, such as higher density and strength material from impact compaction and/or thermal processing and differentiation.Low strength is more consistent with early formation at low collision velocities (Skorov & Blum 2012) in a dynamically cold disk, whilst a homogeneity between 67P's two lobes could imply a similar formation mechanism for both.

Fig. 1 .
Fig.1.Overhang measurement method using the coordinates of three facets, at the base, top and face, of the overhang to define a plane and finding its intersection with the shape model.Aligning the profile with the local gravity vector, g, then allows the various overhang parameters to be measured, as in Tokashiki & Aydan (2010).

Fig. 2 .
Fig.2.Lower limit of tensile strengths for the measured overhangs.On the right the data have been scaled from the feature length scale (h) to the equivalent strength at the metre scale.

Fig. 4 .
Fig. 4. Location of the measured overhangs.The colour scheme shows the (unscaled) tensile strength on a log scale.

Fig. 5 .
Fig. 5. Location of the measured overhangs on the shape model.

Table 1 .
Locations (in the Cheops frame) and properties of each measured overhang.h and L are the heights and depths, as directly measured from the profiles, while σT is derived by numerically integrating the profile shape with Eqn. 1.