Analysis of layering-related linear features on comet 67P/Churyumov-Gerasimenko
Birko-Katarina Ruzicka, Luca Penasa, Hermann Boehnhardt, Andreas Pack, Benoit Dolives, Fabrice Souvannavong, Emile Remetean
aa r X i v : . [ a s t r o - ph . E P ] J a n MNRAS , 1–6 (2018) Preprint 3 January 2019 Compiled using MNRAS L A TEX style file v3.0
Analysis of layering-related linear features on comet67P/Churyumov-Gerasimenko
Birko-Katarina Ruzicka, , ⋆ Luca Penasa, Hermann Boehnhardt, Andreas Pack, Benoit Dolives, Fabrice Souvannavong, and Emile Remetean Max-Planck-Institut f¨ur Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 G¨ottingen, Germany Center of Studies and Activities for Space, CISAS, ‘G. Colombo’, University of Padova, Via Venezia 15, 35131 Padova, Italy Georg-August-Universit¨at G¨ottingen, Geowissenschaftliches Zentrum, Abteilung Isotopengeologie, Goldschmidtstraße 1, 37073 G¨ottingen, Germany Magellium, Imaging and Applications Department, 24 Rue Herm`es, BP 12113, 31521 Ramonville Saint-Agne Cedex, France Engineering for Future Missions Department, Centre National d’ ´Etudes Spatiales, Bpi 1712, 18 av. E. Belin, 31401 Toulouse Cedex 9, France
This article has been accepted for publication in MNRAS, Published by Oxford University Press, with the DOI 10.1093/mnras/sty3079
ABSTRACT
We analysed layering-related linear features on the surface of comet 67P/Churyumov-Gerasimenko (67P) to determine the internal configuration of the layerings within thenucleus. We used high-resolution images from the OSIRIS Narrow Angle Camera on-board the Rosetta spacecraft, projected onto the SHAP7 shape model of the nucleus,to map 171 layering-related linear features which we believe to represent terrace mar-gins and strata heads. From these curved lineaments, extending laterally to up to 1925m, we extrapolated the subsurface layering planes and their normals. We furthermorefitted the lineaments with concentric ellipsoidal shells, which we compared to the es-tablished shell model based on planar terrace features. Our analysis confirms that thelayerings on the comet’s two lobes are independent from each other. Our data is notcompatible with 67P’s lobes representing fragments of a much larger layered body.The geometry we determined for the layerings on both lobes supports a concentricallylayered, ‘onion-shell’ inner structure of the nucleus. For the big lobe, our results arein close agreement with the established model of a largely undisturbed, regular, con-centric inner structure following a generally ellipsoidal configuration. For the smalllobe, the parameters of our ellipsoidal shells differ significantly from the establishedmodel, suggesting that the internal structure of the small lobe cannot be unambigu-ously modelled by regular, concentric ellipsoids and could have suffered deformationalor evolutional influences. A more complex model is required to represent the actualgeometry of the layerings in the small lobe.
Key words: comets: general – comets: individual: 67P/Churyumov-Gerasimenko –methods: data analysis – techniques: miscellaneous
In-situ images of the nucleus surfaces of Jupiter-familycomets (e.g., 9P, 81P, 103P) have previously been used tospeculate about a layered structure of these nuclei (e.g.,Thomas et al. 2007; Bruck Syal et al. 2013; Cheng et al.2013). Belton et al. (2007) proposed a model in which theirinterior consists of a core overlain by layerings that werelocally and randomly piled onto the core through colli-sions (the ‘talps’ or ‘layered pile’ model). More recently,Belton et al. (2018) suggested that the layerings could havebeen formed by fronts of self-sustaining amorphous to crys- ⋆ E-mail: [email protected] talline ice phase-change propagating from the nucleus sur-face to the interior.The ongoing debate about the origin of layerings incometary nuclei might benefit from a more comprehen-sive understanding of their geometry and orientation.This understanding was greatly improved through thedata collected by ESA’s Rosetta mission to comet67P/Churyumov-Gerasimenko (67P). Images taken bythe OSIRIS camera system, surpassing the spatial res-olution of previous missions by more than an order ofmagnitude, resolved features exposed across most of thenucleus surface which we interpret as layerings. Repetitivestaircase patterns are formed by laterally persistent cliffsseparating planar ‘terrace’ surfaces. The cliff faces display © Ruzicka et al. parallel linear grooves, which are reminiscent of sedimentaryoutcrops on Earth where differential erosion carves suchgrooves into layerings of alternating hardness. The dust-freewalls of deep pits reveal quasi-parallel sets of lineaments,oriented roughly perpendicular to the local gravity vectorand extending to depths of at least a hundred meters belowthe present-day nucleus surface (Vincent et al. 2015).A first systematic analysis of the orientation of 67P’slayerings was conducted by Massironi et al. (2015). Theycreated a series of geologic cross-sections of the comet nu-cleus based on the orientation of planes fitted to morpho-logically flat areas (‘terraces’) on a shape model of the nu-cleus surface. Massironi et al. (2015) concluded that the twolobes of the nucleus are independently-formed bodies withan ‘onion-shell’ layered inner structure that formed beforethe two bodies merged in a gentle collision to form the nu-cleus of 67P.Using a similar approach, albeit on a much smaller num-ber of measurements, Rickman et al. (2015) suggested thatmorphological ridges and other features on opposing sides ofthe comet’s lobes can be connected by planar features. Theyinterpreted this correlation as evidence for a semi-planar,pervasive internal layering. This would suggest that the twolobes of 67P are fragments of a much larger body.Penasa et al. (2017) modeled the inner layered structureof comet 67P by fitting concentric ellipsoidal shells to a totalof 483 terraces on both lobes, providing a first simplified3D geological model of their inner structure. By comparingthe orientation of the surface planes with the two modelellipsoids, they suggest that the inner structure of the twolobes can be explained by a set of concentric ellipsoids.The studies by Massironi et al. (2015), Rickman et al.(2015), and Penasa et al. (2017) made use of terraces as aproxy for the underlying orientation of the layers. The ad-vantage of such an approach is that terraces are ubiquitouson the cometary body, that they can easily be mapped, andthat their orientation can be easily estimated by means ofa best fitting procedure applied to the vertices of a shapemodel (Massironi et al. 2015). A downside of the terrace ap-proach is that results may be biased due to airfall and masswasting processes from nearby cliffs (Figure 1). Penasa et al.(2017) acknowledge this limitation and estimate that thismight introduce an error of up to ∼ ° to their normals.Here we use layering-related linear features instead ofterraces to avoid possible bias due to depositional processes.We studied two types of linear features (Figure 1):i) morphological edges along terrace margins, adjacentto cliffs (referred to as ‘terrace margins’) andii) lineaments on hill slopes and cliff faces (referred toas ‘strata heads’, in accordance with Lee et al. 2016) pro-duced by the intersection of the layers with the topographicsurface. Both types of features appear to be erosional conse-quences of discontinuities between the individual layerings,which are possibly related to planes of different physicalproperties within layered materials. The linear traces markthe locations where the inner bedding planes intersect withthe topography at a high angle (Massironi et al. 2015). Figure 1.
Schematic illustration of the two types of layering-related linear features we analyzed in this paper: i) morpholog-ical terrace margins (red curves) and ii) strata heads visible aslineaments on cliff faces and hill slopes (blue curves). Beddingorientations derived from these linear features are not affected bydeposits on the terraces.
For the three-dimensional representation of the nucleusof 67P we used the SHAP7 shape model, which is ob-tained from a stereo-photogrammetric (SPG) analysis of im-ages taken by the OSIRIS Narrow Angle Camera (NAC;Keller et al. 2007) onboard Rosetta. The model covers thewhole nucleus and consists of about 44 million facets, hasa mean accuracy of 0.3 m at a horizontal sampling ofabout 1-1.5 m, and vertical accuracy at the decimeter scale(Preusker et al. 2017).Most of our mapping was conducted on high-resolutionimages of the nucleus surface. We selected suitable OSIRISNAC images, publicly available and retrieved from the ESAPlanetary Science Archive (PSA), according to these crite-ria: Images taken with the OSIRIS NAC orange filter F22(which have a high signal-to-noise ratio; images taken at aspacecraft-comet distance of <
30 km (resulting in a pixelresolution of 0.5 to 1.5 m at the image centre); images wherethe illumination is both sufficiently bright and also in a suit-able direction to show layering-related features in good con-trast. Initially, we selected only images that were calibratedfor geometrical distortion due to the internal camera geom-etry (CODMAC level 4; Tubiana et al. 2015). We later de-cided to supplement these with geometrically uncalibratedimages (level 3) in order to improve coverage of the nucleus.We minimized the impact of this distortion on our data byrestricting our mapping efforts to the central area of uncal-ibrated images, where the distortion decreases to near-zero.
Mapping of linear features
Large-scale morphological terrace margins were mapped di-rectly on the shape model by manually tracing the featuresas polylines in
CloudCompare (Girardeau-Montaut 2014).Mapping the finer morphological details of strata headsrequired a higher resolved mapping medium, for which weprojected two-dimensional OSIRIS images onto the three-
MNRAS , 1–6 (2018) ayering-related linear features on comet 67P dimensional shape model. We used a customized versionof the software philae localization workshop (PLW,Remetean et al. 2016) for the measurement process of thelinear features. For each image, we manually selected a setof reference points, consisting of prominent landmarks visi-ble on both the image and the shape model (e.g., large boul-ders, sharp corners), and used the software to spatially alignthe image with the shape model by minimizing the distancebetween corresponding reference points. Using a set of 20 ref-erence points we achieved a root-mean square error (RMSE)between 3.1 and 9.9 m (7.0 m on average) for the alignment,depending on how many clearly defined landmarks are visi-ble on the image that is to be aligned. Increasing the numberof reference points did not further improve the RMSE. Thepoint of view onto the projected image can then be changedin 3D to allow mapping from an optimum viewing angle.Subsequently, we manually selected nodes along eachfeature of interest. The nodes each have three-dimensionalXYZ point-coordinates in the comet-fixed Cheops referenceframe. For ease of handling we exported the nodes as onepolyline per feature. To aid visualization of the layering ori-entations, we determined best matching local plane solutionsfor each linear feature. The plane solutions consist of thenormal vector ( n ) and the reference point of the plane ( nb ),which is the centroid of the nodes used to fit the plane. Italso corresponds to the location of the base of n . We foundthese plane solutions by applying a weighted least squaresplane-fitting routine to the nodes of the feature (Planefit,Schmidt 2012) in matlab (The MathWorks 2017, release2017a)).We assessed the uncertainty of the plane normal vec-tors n by means of a Monte Carlo simulation: The valuesof coordinates X , Y , and Z of each node were varied, by anormally-distributed random value within the error of thenodes, around the measured coordinates X m , Y m , and Z m .Those planes where n was poorly defined (variance of MonteCarlo results in at least one direction of ≥ ° ) were dis-carded from the pool of mapped features. Ellipsoidal model fitting
To evaluate whether linear features can be used as astandalone product to produce three-dimensional geologi-cal models, we made use of the same model defined byPenasa et al. (2017). The model represents the layering ofeach lobe as a scalar field: f ( c x , c y , c z ) = k (1)Function f () is defined such that for any constant value of k a single contour surface with ellipsoidal shape is determined,while the value of the scalar field represents a metric in thestratigraphic space ( in the centre of the ellipsoidal modeland increasing toward the outermost layers).Function f () is completely defined by eight parameters: c x , c y , and c z for the centre of the ellipsoidal model, b and c for the axial ratios and finally α, β , and γ for orientingthe concentric ellipsoids in space. The parameters can bedetermined by maximizing the accordance of the model withthe provided constraints. The orientation of the terraces,were used by Penasa et al. (2017) to provide observations ofthe gradient of the function f () in a specific point p , thus providing a modeling constraint of the type: ∇ f ( p ) = n (2)where n is the normal to the surface element located in thepoint p .In this work we instead tested the use of polylines formodelling purposes. Each polyline is formed by segmentswhich are expected to lie on a contour surface of the modeland are thus tangential to the ellipsoidal shell passing forthat location. Each segment can be defined by a pair ofpoints p and p describing a direction in space, which canbe represented as a unit vector: ˆ t = p − p k p − p k (3)where ˆ t represents a constraint of tangent type (e.g.Hillier et al. 2014, on this subject): h∇ f ( p ) , t i = (4)where p is the location of the observation (i.e. a referencepoint for the location of the segment). From these observa-tions an angular misfit of the observation in respect to anymodel can be obtained by computing: ϑ = arccos (cid:18)(cid:13)(cid:13)(cid:13)(cid:13) ∇ f ( p )k∇ f ( p )k × ˆ t (cid:13)(cid:13)(cid:13)(cid:13)(cid:19) (5)By minimizing the squares of the angles provided byEquation 5 for each segment composing the mapped poly-lines, we were able to determine the most-likely parametersfor the ellipsoidal model approximating the observations inthis work. We employed a weighting strategy to ensure thateach polyline contributes equally to the obtained solution.For this we divided each squared residual by the total num-ber of segments of the specific polyline. Finally, we used abootstrap strategy, based on the resampling of the polylines,to estimate the standard errors associated with each param-eter. We used the PLW software to align a total of 34 OSIRISimages (
Table A2 in the supplementary material), cover-ing most of the nucleus surface. On these images we mapped171 linear features, of which 31 are terrace margins and 140are strata heads. The mapped features are distributed ap-proximately evenly between 67P’s big and small lobe. Thefeatures extend laterally for 863 m on average, ranging from185 to 1925 m. The feature length is calculated from thecumulative length of all segments connecting the polylinenodes. Most features contain between 9 and 20 segments(14 on average), and the average segment length is 38 m.The uncertainty of each node results cumulatively fromthe error of the shape model ( ± <
10 m.
MNRAS , 1–6 (2018)
Ruzicka et al.
Table 1.
Exemplary excerpt of the normal vectors n and theirbases nb to the plane solutions of the mapped strata heads onthe big lobe. r is the radial distance of each vector base from thegravitational centre of the lobe. Values are in [km] in the comet-fixed Cheops reference frame. The complete tables are availablein the supplementary material. nb x nb y nb z n x n y n z r0.316 0.695 0.616 0.632 -0.425 0.648 1.3-0.311 0.740 0.696 0.840 -0.117 0.530 1.40.473 0.981 -1.320 0.013 0.131 -0.991 1.90.489 0.536 -1.257 0.919 0.025 -0.393 1.70.290 -0.074 -1.237 0.776 -0.084 -0.625 1.5... ... ... ... ... ... ... The Monte Carlo error analysis showed that affectingthe node coordinates X , Y , and Z by random amounts be-tween zero and 10 m results in a variance of the normals n by less than 10 ° for most features ( Figure A1 in thesupplementary material). For a small number of features,whose polylines have a low curvature, n is more stronglyaffected. As expected in such cases, the uncertainties showsignificant directional asymmetry and have large amplitudesonly perpendicular to the main extension of the lineament.As Figure A1(B) shows, these cases are the exception inour data, and the node uncertainty does not have a sub-stantial effect on the bulk of the layering orientations wereconstructed from the linear features.An exemplary excerpt of the parameters of the planesolutions to the mapped features are listed in Table 1 (thecomplete
Table A1 is available in the supplementary ma-terial). The 3D orientation of vectors n is illustrated inFigure 2. Qualitatively, by visual impression the normals areoriented perpendicular to the nucleus surface. For the bigand small lobe separately, the normals are pointing outwardfrom the respective lobe’s gravitational centre.For the purpose of comparison, we fitted our own set ofellipsoids to the polylines of our linear features. The param-eters we obtained for the best-fitting ellipsoidal model aresummarized in Table 2, next to the parameters achieved byPenasa et al. (2017). The parameters are consistent for thebig lobe within the achieved uncertainties, but we found anotable a misfit of the results for the small lobe. Our smalllobe model is offset by 0.32 km from their model (whichamounts to ca. 40% of the lobe’s semi-minor axis, accordingto Jorda et al. (2016)), there is a minor difference in the el-lipsoidal axis ratio, and a major mismatch in the rotationalangles.Finally, we compared the orientation of the plane nor-mals n to the orientation of the ellipsoid surface for bothPenasa et al. (2017) and our model. Figure 3 shows the an-gular misfit between n and the corresponding normal to theellipsoid surface n ell , both at location nb . Again, we findan overall low angular misfit for the big lobe (median 16.7 ° ,panel A), and a larger misfit for the small lobe (median 19.6 ° with larger percentiles, panel B), indicating that the linearfeatures are less congruent with the ellipsoid on the smalllobe. Penasa et al. (2017) did not observe this dichotomy. Figure 2.
The SHAP7 shape model of comet 67P, showing thefeature normals n at the bases nb for all mapped terrace margins(blue arrows) and strata heads (pink arrows). The dotted linesrepresent the mapped nodes along each feature. Coordinates arein the Cheops reference frame. A: ‘Front’ view towards positivey-values; B: ‘Back’ view rotated around the z-axis by 180 ° . Full-resolution images are available in the supplementary materials. While the layering-related linear features are not affected bysedimentation, we cannot rule out that the subset of terracemargins might be influenced by erosion and cliff collapse(cf. e.g. Pajola et al. 2017). However, we took care to mini-mize this effect by accurately following the external border ofthe mapped edges, including any niches resulting from localbreakoffs, thus preserving the orientation of the underlyinglayering.The orientation of our feature normals n particularlyin the ‘neck’ region of the nucleus (Figure 2), supports thewidely accepted findings of previous works that the orien-tation of the layerings on the big and small lobe of 67Pare independent from each other (e.g. Massironi et al. 2015;Davidsson et al. 2016); it does not support a common enve-lope structure surrounding both lobes, nor the interpretationthat both lobes represent fragments of a much larger, layeredbody (as proposed by Rickman et al. 2015).We find that the internal structure of the nucleus, asdeduced from the orientation of layerings mapped at thesurface, is in agreement with the ‘onion-shell’ model pro-posed by Massironi et al. (2015) and the concentric ellip- MNRAS , 1–6 (2018) ayering-related linear features on comet 67P Table 2.
Maximum likelihood estimates of the ellipsoidal model parameters and relative 2 σ errors, modelled to our data and comparedto Penasa et al. (2017). Distances in km in the Cheops reference frame, and Tait-Bryan angles in degrees. BL means 67P’s big lobe, SLis the small lobe; b and c are the axial ratios with respect to the a-axis ( a = 1).Values for which this study differs from Penasa et al.(2017) in excess of the uncertainties are highlighted in bold in the table.Parameter This study (171 linear features) Penasa et al. (2017) (483 terraces)BL 2 σ SL 2 σ BL 2 σ SL 2 σ c x -0.55 0.12 c y c z -0.15 0.11 α β γ Figure 3.
Histogram of the angular misfit between normal vec-tors n to the surface features and corresponding normal vectors n ell of the ellipsoid model. The misfit for our linear featuresis plotted in orange, the misfit for the terraces in Penasa et al.(2017) is plotted in blue. A: Features on 67P’s big lobe, B: Fea-tures on the small lobe. soidal shell model by Penasa et al. (2017). Particularly forthe big lobe, our results match those of Penasa et al. (2017)closely. We understand this as confirmation that the biglobe has a regular, concentric inner structure that followsa generally ellipsoidal configuration. Nevertheless, our datacannot confirm that the layerings are indeed connected intoglobally coherent shells. Approximating the minimal lobecircumferences as 6 km and 8.6 km, respectively, our mea-surements (including polylines of almost 2 km length) mightbe mistaken to span a substantial portion of the nucleussurface. However, most of our polylines intentionally have a strong curvature and represent features with a continu-ous lateral extent of no more than a few hundred metres.This leaves room for the possibility of a discontinuously lay-ered structure, as would be a consequence of e.g. the ‘talps’model (Belton et al. 2007) or layering by thermal processes(Belton et al. 2018). Our findings would be also compatiblewith either concept.For the small lobe, our results differ significantly fromthose of the other authors. Neither are the orientations of ourproposed layerings clearly compatible with the ellipsoidalshell model based on its terraces, nor do the parameters ofour own ellipsoidal model agree with those of the terrace-based model. This disagreement could either be explainedby the circumstance that Penasa et al. (2017) mapped theirterraces exclusively on the shape model, and thereby mighthave included some planar areas that are unrelated to thelayerings. Another possible explanation is that the smalllobe’s inner structure has been affected by processes of evo-lution or deformation, and thus cannot be unambiguouslymodelled by regular, concentric ellipsoids. In this case, amore complex model is required to represent the real geom-etry of the layerings.
ACKNOWLEDGEMENTS
BKR was financially supported by the International MaxPlanck Research School (IMPRS) hosted by the Max PlanckInstitute for Solar System Research (G¨ottingen). This workbenefited from insightful discussions with Matteo Massironiand his research group at the University of Padua. Geo-metrically calibrated images were made available by theOSIRIS team. The PLW software was provided by the Cen-tre National d’´Etudes Spatiales (CNES) in collaborationwith Magellium. Borys Dabrowski is thanked for his helpwith the matlab coding.
REFERENCES
Belton M. J. S., et al., 2007, Icarus, 191, 573Belton M. J., Zou X.-D., Li J.-Y., Asphaug E., 2018, Icarus, 314,364Bruck Syal M., Schultz P. H., Sunshine J. M., A’Hearn M. F.,Farnham T. L., Dearborn D. S. P., 2013, Icarus, 222, 610MNRAS , 1–6 (2018)
Ruzicka et al.
Cheng A. F., Lisse C. M., A’Hearn M., 2013, Icarus, 222, 808Davidsson B. J. R., Sierks H., G¨uttler C., Marzari F., Pajola M.,Rickman H., 2016, Astronomy & Astrophysics, 592Girardeau-Montaut D., 2014, Cloudcompare, a 3D Point Cloudand Mesh Processing Free Software. EDF R&D, TelecomParisTech.Hillier M. J., Schetselaar E. M., de Kemp E. A., Perron G., 2014,Mathematical Geosciences, 46, 931Jorda L., et al., 2016, Icarus, 277, 257Keller H. U., et al., 2007, Space Science Reviews, 128, 433Lee J.-c., et al., 2016, MNRAS, pp 1–17Massironi M., et al., 2015, Nature, 526, 402Pajola M., et al., 2017, Nature Astronomy, 1, 0092Penasa L., et al., 2017, Monthly Notices of the Royal Astronomical Society,14, 1Preusker F., et al., 2017, Astronomy and Astrophysics, 1, 1Remetean E., Dolives B., Souvannavong F., Germa T., GinestetJ., Torres A., Mousset T., 2016, Acta Astronautica, 125, 161Rickman H., et al., 2015, Astronomy & Astrophysics, 583, A44Schmidt V., 2012, Planefit,
The MathWorks I., 2017, MATLABThomas P. C., et al., 2007, Icarus, 191, 51Tubiana C., et al., 2015, Astronomy and Astrophysics, 583, A46Vincent J.-B., et al., 2015, Nature, 523, 63
APPENDIX A: SUPPLEMENTARY MATERIAL
Supplementary material is available online. It includes
Fig-ures 1 , and in full resolution, and Figure 2 as a mat-lab figure, and the complete table of feature normals (
TableA1 ). Table A2 lists all OSIRIS NAC images that were usedfor mapping, with the corresponding CODMAC calibrationlevels.In addition, we included
Figure A1 , which is a visual-ization of the results of the Monte Carlo analysis of our planenormals, as described in the text. A: Plot of all Monte Carlosolutions. Each set of green arrows represents the spread ofthe normal vectors of a single linear feature. B: The MonteCarlo solutions of most features show low variability, as seenin this histogram of the standard deviations (2 σ ). This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000