Airfall on Comet 67P/Churyumov–Gerasimenko

Icarus 354 (2021) 114004 Contents lists available at ScienceDirect Icarus journal homepage: www.elsevier.com/locate/icarus Airfall on Comet 67P/Churyumov–Gerasimenko Björn J.R. Davidsson a, *, Samuel Birch b, Geoffrey A. Blake c, Dennis Bodewits d, Jason P. Dworkin e, Daniel P. Glavin f, Yoshihiro Furukawa g, Jonathan I. Lunine h, Julie L. Mitchell i, Ann N. Nguyen j, Steve Squyres k, Aki Takigawa l, Jean-Baptiste Vincent m, Kris Zacny n a Jet Propulsion Laboratory, California Institute of Technology, M/S 183–401, 4800 Oak Grove Drive, Pasadena, CA 91109, USA Cornell University, 426 Space Sciences Building, Ithaca, NY 14850, USA Division of Geological & Planetary Sciences, California Institute of Technology, MC150–21, Pasadena, CA 91125, USA d Department of Physics, Auburn University, 201 Allison Laboratory, Auburn, AL 36849, USA e NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA f NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA g Department of Earth Science, Tohoku University, 6–3, Aza–aoba, Aramaki, Aoba–ku, Sendai 980–8578, Japan h Department of Astronomy and Carl Sagan Institute, Cornell University, Ithaca, NY 14850, USA i NASA Johnson Space Center, 2101 NASA Pkwy, Houston, TX 77058, USA j Jacobs, Astromaterials Research and Exploration Science, NASA Johnson Space Center, 2101 NASA Parkway, Mail Code XI3, Houston, TX 77058, USA k Blue Origin, LLC, 21218 76th Ave S, Kent, WA 98032, USA l The Hakubi Center for Advanced Research, Division of Earth and Planetary Science, Kyoto University, Kitashirakawa–Oiwakecho, Sakyo, Kyoto 606–8502, Japan m DLR Insitute of Planetary Research, Rutherfordstrasse 2, 12489 Berlin, Germany n Honeybee Robotics, 398 W Washington Blvd., Pasadena, CA 91103, USA b c A R T I C L E I N F O A B S T R A C T Keywords: 67P/Churyumov–Gerasimenko Comets, composition Comets, dust Comets, nucleus We here study the transfer process of material from one hemisphere to the other (deposition of airfall material) on an active comet nucleus, specifically 67P/Churyumov–Gerasimenko. Our goals are to: 1) quantify the thickness of the airfall debris layers and how it depends on the location of the target area, 2) determine the amount of H2O and CO2 ice that are lost from icy dust assemblages of different sizes during transfer through the coma, and 3) estimate the relative amount of vapor loss in airfall material after deposition in order to understand what locations are expected to be more active than others on the following perihelion approach. We use various numerical simulations, that include orbit dynamics, thermophysics of the nucleus and of individual coma aggregates, coma gas kinetics and hydrodynamics, as well as dust dynamics due to gas drag, to address these questions. We find that the thickness of accumulated airfall material varies substantially with location, and typically is of the order 0.1–1 m. The airfall material preserves substantial amounts of water ice even in relatively small (cm–sized) coma aggregates after a rather long (12 h) residence in the coma. However, CO2 is lost within a couple of hours even in relatively large (dm–sized) aggregates, and is not expected to be an important component in airfall deposits. We introduce reachability and survivability indices to measure the relative capacity of different regions to simultaneously collect airfall and to preserve its water ice until the next perihelion passage, thereby grading their potential of contributing to comet activity during the next perihelion passage. * Corresponding author. E-mail addresses: [email protected] (B.J.R. Davidsson), [email protected] (S. Birch), [email protected] (G.A. Blake), [email protected] (D. Bodewits), [email protected] (J.P. Dworkin), [email protected] (D.P. Glavin), [email protected] (Y. Furukawa), [email protected]. edu (J.I. Lunine), [email protected] (J.L. Mitchell), [email protected] (A.N. Nguyen), [email protected] (A. Takigawa), [email protected] (J.-B. Vincent), [email protected] (K. Zacny). https://doi.org/10.1016/j.icarus.2020.114004 Received 13 January 2020; Received in revised form 15 July 2020; Accepted 21 July 2020 Available online 1 August 2020 0019-1035/© 2020 Elsevier Inc. All rights reserved. B.J.R. Davidsson et al. Icarus 354 (2021) 114004 1. Introduction material is rapidly cleaned off. Observations by the mass spectrometer ROSINA (Hässig et al., 2015; Fougere et al., 2016a) and the near–infrared spectrometer VIRTIS (Fink et al., 2016) on Rosetta show that the two hemispheres display a strong chemical dichotomy as well. The northern hemisphere is predominantly outgassing H2O and comparably small amounts of CO2, while the southern hemisphere is a source of both water and carbon dioxide. Keller et al. (2017) interpreted the chemical dichotomy as a result of airfall. In their view, the CO2 is either missing in the solid material being ejected into the coma from the south near perihelion (i.e., the CO2 sublimation front is located at some depth) or is lost on the way during transport toward the north. Because of the substantially higher sublimation temperature of H2O compared to CO2 the water loss is significantly smaller. The airfall that is responsible for northern activity in other parts of the orbit is therefore rich in water ice but poor in supervolatiles like CO2, according to Keller et al. (2017). High–resolution imaging of smooth terrain in Ash by OSIRIS during low (~10 km) Rosetta orbits (Thomas et al., 2015a), by the Philae/ ROLIS camera during descent toward Agilkia in Ma’at (Mottola et al., 2015; Pajola et al., 2016), and by OSIRIS during the landing of Rosetta at Sais in Ma’at (Pajola et al., 2017b) revealed that the material in those locations consisted primarily of pebbles, cobbles, and boulders in the cm–m size range. In the following, such units are collectively referred to as “chunks”. The observations of similarly sized chunks in the coma (e. g., Rotundi et al., 2015; Davidsson et al., 2015), of which a substantial fraction display acceleration toward the nucleus (Agarwal et al., 2016), show that the airfall concept is viable. Smooth terrains are not featureless. Structures resembling aeolean ripples were observed in Hapi, dunes that in some cases contained pits (potentially formed by sublimation) were seen in Serqet and Maftet, and some boulders in Hapi and Maftet appeared to have wind tails (Thomas et al., 2015a), that are also seen in Ma’at (Mottola et al., 2015). Whether such features primarily form during airfall deposition or are the result of lateral transport mechanisms is unclear. They do indicate that deposition is not a trivial phenomenon and that significant local transport may take place after deposition. The features seen in smooth terrain are also not static. The first reported observations of large–scale morphological changes concerned roundish scarps that formed and expanded at ~6 m day−1 in Imhotep (Groussin et al., 2015). Later, similar scarp retreats were also seen in Hapi, Anubis, and in Seth (El-Maarry et al., 2017; Hu et al., 2017). The scarps frequently displayed strong brightening and spectral slope changes suggestive of the presence of abundant sub–surface water ice that mixed up to the surface during the retreat process. A dramatic and puzzling example of these morphological changes was the removal of the aeolean ripples in Hapi in April–July 2015 by retreating scarps, that was followed by ripple reformation in December 2015 (ElMaarry et al., 2017). OSIRIS observations allowed for a documentation of inbound removal and outbound redeposition of airfall material, albeit restricted to specific locations and time instances because of resolution, viewing geometry, and imaging cadence restrictions. A comprehensive summary is found in Hu et al. (2017). The presence of “honeycomb” features in Ash, Babi, Serqet, Seth, and particularly in Ma’at in late March 2015 (also see Shi et al., 2016) caused by the removal of overlying granular material is indicative of self–cleaning. The time–scale of their formation is uncertain because of observational biases but the earliest known honeycomb sighting occurred two months earlier. Other signs of airfall deposit removal include exposure of previously hidden outcrops and boulders, formation of shallow depressions, and smoothing of previously pitted deposits (Hu et al., 2017). The OSIRIS images were used to create three–dimensional maplets of selected terrains before and after the removal of deposits. Although it is not possible to estimate the thickness of the expelled layers by simply subtracting two maplets, insight about the erosion can still be obtained by measuring changes in the degree of roughness inferred from such maplets. Hu et al. (2017) find, for regions where changes are more easily detected, that at least 0.5 m has been removed (similar to the pixel resolution) and an estimated average The initial characterization of Comet 67P/Churyumov–Gerasimenko (hereafter 67P) by the OSIRIS cameras (Keller et al., 2007) on the Rosetta orbiting spacecraft (Glassmeier et al., 2007) revealed widespread smooth terrains on the northern hemisphere of the nucleus, primarily in Ma’at on the small lobe, in Ash, Anubis, and Imhotep on the large lobe, and in Hapi that constitutes the northern neck region between the two lobes (Sierks et al., 2015; Thomas et al., 2015b). The smooth material in the Ma’at and Ash regions formed a relatively thin coverage over partially–revealed consolidated landforms. This led Thomas et al. (2015b) to propose that the smooth material, at least in those regions, constituted a rather recent veneer of dust that had rained down from the coma. They named these deposits and the process forming them “airfall”. The presence of smooth deposits in isolated topographic lows, e.g., in Khepry (El-Maarry et al., 2015b), or in large gravitational lows found primarily in the Imhotep and Hapi regions, suggest that airfall may not be the only mechanism responsible for the formation of smooth terrain. Auger et al. (2015) propose that the vast smooth plain in Imhotep formed through mass wasting the surrounding steep cliffs and by transport downhill toward the lowland. Mass wasting from cliffs revealed by taluses, gravitational accumulation deposits, and diamictons that extends into Hapi (Pajola et al., 2019) indicate that such processes are partially responsible for the smooth material in that region as well. The presence of deposit–free regions that sharply contrast with surrounding smooth terrain, best illustrated by the Aten depression on the large lobe or the Anuket region in the neck area (El-Maarry et al., 2015b), suggests that the removal rate of airfall material through self–cleaning is locally high and that the net accumulation rate may be slow (e.g., in case the Aten depression formed recently in a massive outburst event). At the time of the Rosetta orbit insertion around 67P in August 2014 the northern hemisphere was illuminated while the southern hemisphere experienced polar night because of the orientation of the spin axis (Sierks et al., 2015). Most activity, as revealed by prominent dust jets, came from Hapi (Sierks et al., 2015; Lara et al., 2015), that also was the brightest and bluest unit in terms of spectral slope (Fornasier et al., 2015), showing that the neck region had accumulated particularly ice–rich material that was strongly active. Thomas et al. (2015b) noted that the airfall deposits on north–facing regions are more extensive than on south–facing regions and proposed that the major transport route went from Hapi to other suitably oriented regions on the northern hemisphere. During approach to the inbound equinox in May 2015, increasingly large portions of the southern hemisphere were illuminated and eventually received peak solar illumination at the perihelion passage in August 2015 while the northern hemisphere experienced polar night. Keller et al. (2015a) demonstrated that ~80% of the solar radiation is absorbed by the northern hemisphere at a low rate over an extended amount of time, while the southern hemisphere receives the remaining ~20% during a much shorter time interval. Their thermophysical modeling showed that the strongly non–linear response of sublimation to solar flux levels led to a ~3 times higher erosion in the south compared to the north. They therefore suggested that the main transport route of airfall material went from the southern hemisphere to the northern hemisphere. The views of Thomas et al. (2015b) and Keller et al. (2015a) are not in contradiction, but suggest that activity in the south around perihelion leads to airfall in the cold and inactive north. Some of this material is then redistributed from Hapi to other northern regions from the time of the outbound equinox to the shutdown of activity at larger heliocentric distance, with a pause during the aphelion passage followed by resumed redistribution when activity reappears during inbound motion. OSIRIS observations of the southern hemisphere showed a strong hemispherical dichotomy regarding smooth terrains – they were essentially missing in the south (El-Maarry et al., 2016; Lee et al., 2016; Birch et al., 2017), showing that any deposited or otherwise produced smooth 2 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 removal of 1.3 m for one specific honeycomb. Post–perihelion observations showed that the honeycombs had disappeared and were replaced by smooth terrain. Hu et al. (2017) do not quantify the thickness of the redeposited layer but a reasonable assumption is of order 0.1–1 m, so that the net orbital change is small. It is not self–evident that the northern airfall deposits currently experience net growth, or if there is currently net erosion of deposits built up at an earlier time when conditions were different than today. From the description above a basic qualitative understanding has emerged regarding the origin and behavior of smooth terrain on 67P; this new insight emphasizes the global redistribution of material on the nucleus enabled by its spin axis obliquity. This process creates a non–primordial chemical surface heterogeneity. If these processes are active on other comet nuclei as well it makes the interpretation of groundbased production rate patterns more difficult than previously envisioned. However, it is also evident that a detailed understanding of these processes requires substantial modeling efforts that can fill the gaps that are not immediately provided by available observations. Therefore, the first goal of this paper is to present a quantitative study of the south–to–north transfer process proposed by Keller et al. (2015a) and elaborated on by Keller et al. (2017), i.e., to estimate the thickness of the deposited layer built during the perihelion passage and how it varies across the northern hemisphere. Specifically, we focus on 31 target areas distributed within the Ma’at region on the small lobe and within the Ash and Imhotep regions on the large lobe. These sites were selected and investigated in the context of a Phase A study for the New Frontiers 4 candidate mission CAESAR (Comet Astrobiology Exploration SAmple Return). The sites were selected to cover a range of regions of smooth terrain in the northern hemisphere of 67P that were deemed to be accessible to the CAESAR touch and go sampling. The location of the 31 target areas are shown in Fig. 1. Our second goal is to quantify the loss of volatiles expected to take place during the transfer process, thereby enabling a comparison between the ice abundance of fresh airfall deposits with that of the sources of this material. Our third goal is to quantify the relative capacity of different airfall deposition sites to retain their water ice content after deposition, in order to better understand their capability to uphold comet activity as the comet approaches the Sun after the aphelion passage. In the context of the CAESAR mission this work contributed to locating high science value sampling sites in the smooth terrain on 67P for future Earth return. Our work relies on a pipeline of different models and methods that are described in Sec. 2. In brief, we first find Keplerian orbits that connect specific source regions and target areas that are unobstructed by the rotating nucleus. Those orbits are characterized by a specific ejection velocity vd that needs to be provided by gas drag, and are realized for a specific dust radius acrit. We calculate acrit by first determining the local surface temperature and outgassing rate, use that information to calculate the coma number density, translational temperature, and gas expansion velocity versus height, which in turn are used to evaluate the gas drag force needed for dust dynamics calculations. We also present our model for calculating the loss of H2O and CO2 from chunks in the coma, as well as the longterm nucleus thermophysics modeling. Our results are presented in Sec. 3, we discuss these results and compare with Fig. 1. The location of target areas, considered for CAESAR sample collection, on Rosetta/OSIRIS context images with file names provided in the following. Upper left: Ma’at (N20140902T084254618ID30F22.img). Upper right: eastern Ash (N20140822T054254583ID30F22.img). Lower left: western Ash (N20140902T121022552ID30F22.img). Lower right: Imhotep (N20140905T064555557ID30F22.img). 3 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Table 1 Functions and parameters with descriptions and units. Upper case Latin. Model parameters #1 Symbol Description Unit A A Bond albedo Chunk cross section – m2 Eq. (25) constant K (1 − γ)−1 A crit C CD D Dmin, Dmax, Dc D Et Er Fdrag Fg Fıȷ Cross section of maximum liftable chunk Drag coefficient Dust chunk diameter Minimum, maximum, and intermediate dust chunk diameter Molecular number flux Molecular kinetic energy Molecular rotational energy Drag force Gravitational force vector Fraction of latent heat consumed during segregation m2 – m m m−2 s−1 J J N N – F jk F F, G GN Hı View factor Fraction of dust mass at D > Dc Eq. (18) parameters Newtonian gravitational constant Energy release during crystallization – – – N m2 kg−2 J kg Qd Qj, Qs Rn R* ℛ S1 S⊙ S′ S Dust production rate Water sublimation rate Comet nucleus curvature radius Radiogenic energy production rate Reachability index Speed ratio Solar constant Normalized relative ice retainment capacity of target area Survivability index kg m−2 s−1 kg m−2 s−1 m J m−3 s−1 – – J m−2 s−1 – – Hertz–Knudsen formula kg m−2 s−1 Kt L L′ M, M1, M2 Mn Mrot P P T, Tj(x, t), Tk(x, t), Ts, T0, T1, T∞ Tfr Tr Tsurf U W, W1, W2, W∞ W Z Knudsen layer thickness Latent heat of sublimation Normalized relative ice loss of target area Mach number Comet nucleus mass Rotation matrix Nucleus or chunk rotation period Comet orbital period Temperature Freeze–out temperature Rotational temperature Surface temperature Gravitational potential Gas drift speed Speed ratio previously published airfall investigations in Sec. 4, and summarize our conclusions in Sec. 5. m J kg−1 – – kg – h yr K K K K J kg−1 m s−1 – of this information in hand, the amount of airfall material can be calculated, as described in Sec. 2.6. Section 2.7 describes the model used for modeling the thermophysical evolution of chunks in the coma, and the longterm thermophysical evolution of target sites are described in Sec. 2.8. The mathematical symbols used throughout this work are summarized in Tables 1–3. 2. Methods Section 2.1 describes the method for determining whether an unobstructed Keplerian orbit exists that connects a certain southern source region with a specific northern target area, taking the irregular shape of the rotating nucleus into account. The purpose of that calculation is to identify the dust ejection velocity vd that is needed in order to realize a particular orbit. The reliability of this analytical method based on the point–mass assumption is tested using numerical integration of orbits around irregular bodies with complex gravity fields described in Sec. 2.2. The thermophysical model described in Sec. 2.3 is needed in order to calculate the nucleus outgassing rate as function of surface location and time. The outgassing rate is used to calculate how the gas density, drift (expansion) speed, and temperature vary with height, using procedures described in Sec. 2.4. Those quantities are, in turn, necessary to determine the specific chunk size that will be accelerated to the particular vd needed to reach a given target area, as summarized in Sec. 2.5. With all 2.1. Keplerian orbits Let rs be the position vector of a source region on the southern nucleus hemisphere, in an inertial frame with its ̂ z axis aligned with the nucleus positive spin pole, at the time of dust ejection t = 0. Let vs be the dust ejection velocity vector in the same frame with one component due to gas drag acceleration along the local surface normal ̂ n s and another due to nucleus rotation, n s + ̂z × rs ω vs = vd ̂ (1) where the angular velocity is ω = 2π /P and P is the nucleus rotational period. Then the normal to the dust chunk orbital plane is h = rs × vs = {hx, hy, hz}. Let rt′ be the position vector of a target area on the northern nucleus 4 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Table 2 Functions and parameters with descriptions and units. Lower case Latin. Model parameters #2 Symbol Description Unit a, b, c a1, a2, a3 acrit ao cs d f fi gı Eq. (5) coefficients Eq. (11) coefficients Radius of maximum liftable chunk Orbital semi–major axis Specific heat capacity Distance between chunk and target location Eq. (5) coefficient Ice volume fraction Heat capacity of gaseous species ı m3 s−1 – m m J kg−1 K−1 m – – J kg−1 h h = {hx, hy, hz} ı ȷ Height above the comet nucleus surface Chunk orbit normal vector Chemical species identifier (host) Chemical species identifier (guest) m m2 s−1 – – j, k kB l m mı Facet number Boltzmann constant Spherical chunk latitudinal coordinate Water molecule mass Molecule mass of species ı – J K−1 rad kg kg p, q pı Eq. (5) coefficients Partial pressure of species ı – Pa md n, ns, n0, n1, n2 nsp ̂ ns qı qı ′ psat pv qs r rd rg rh rjk ro rs rt = {rtx, rty, rtz} sk t tpl ud vd vesc vo vj⊙, vk⊙ vs x, xj {̂ x, ̂ y , ̂z } Mass of a dust chunk Gas number density Number of chemical species Surface normal unit vector Volume mass production rate of species ı Volume mass segregation or crystallization rate Water saturation pressure Vapor pressure Dust chunk size distribution power–law index Spherical chunk radial coordinate Chunk position vector Constituent grain radius Heliocentric distance Distance between facets Orbital distance between chunk and nucleus center Source region position vector Target area inertial position vector with rt′ = rt(t = 0) Facet area Time Target area crossing time with chunk orbital plane Chunk speed Post–acceleration speed along surface normal Escape velocity Orbital speed Sunlit fraction of facet Dust chunk ejection velocity vector Depth below the comet surface Inertial system coordinate system unit vectors kg m−3 – – kg m−3 s−1 kg m−3 s−1 Pa Pa – m m m AU m m m m m2 s s m s−1 m s−1 m s−1 m s−1 – m s−1 m – (̂z parallel to the nucleus rotation axis) hemisphere in the inertial frame at the dust ejection time t = 0 (equivalently, the target area position vector in the body–fixed frame). At a later time the target area has position rt(t) = Mrot(t)rt′ = {rtx, rty, rtz} in the inertial frame because of nucleus rotation, where ⎛ ⎞ cosωt −sinωt 0 Mrot (t) = ⎝ sinωt cosωt 0⎠ (2) 0 0 1 ( ) h⋅rt tpl = 0, (3) remembering that the orbital plane necessarily must pass through the center of mass of the nucleus that hosts the origin of the coordinate system. The crossing of the target area with the dust orbital plane will repeat twice during a nucleus rotation at instances of time tpl obtained from Eq. (3) and associated equations, (√̅̅̅̅̅̅̅̅̅̅̅̅̅) 1 (4) tpl = sin−1 1 − f2 The criterion for the target area being located within the dust orbital plane is then ω 5 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Table 3 Functions and parameters with descriptions and units. Greek. Model parameters #3 Symbol α αbf αs βj, βk β− γ δ ε ζ κ λ, λ0 μj, μk μo, μ* ρ σ σc τı Φı Ψı ψ ω where ⎧ √̅̅̅̅̅̅̅̅̅̅̅̅̅ ⎪ ⎪ p p2 ⎪ ⎪ f =− ± −q ⎪ ⎪ ⎪ 2 4 ⎪ ⎪ ⎪ ⎪ ⎪ 2ac ⎪ ⎪ ⎪ p= 2 ⎪ ⎪ a + b2 ⎪ ⎨ c2 − b2 ⎪ q= 2 ⎪ ⎪ ⎪ a + b2 ⎪ ⎪ ⎪ ⎪ ⎪ a = hx rtx + hy rty ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ b = −hx rty + hy rtx ⎪ ⎪ ⎪ ⎩ c = rtz hz Description Unit Right ascension Molecular surface backflux fraction Sublimation coefficient Angle between facet normal and direction to other facet Molecular backflux ratio Gas heat capacity ratio Declination Emissivity Number of rotational degrees of freedom Heat conductivity Molecular mean free path Cosine of incidence angle Reduced mass Density Stefan–Boltzmann constant Molecular collisional cross section Phase change rate – – W m−1 K−1 m – kg kg m−3 W m−2 K−4 m2 kg m−3 s−1 Latitudinal gas diffusion rate of species ı kg m−2 s−1 Radial gas diffusion rate of species ı Porosity Angular velocity of nucleus rotation ∘ – – rad – – ∘ kg m−2 s−1 – rad s−1 location and orientation to enable successful trajectories to the 31 northern hemisphere sites. 2.2. Realistic gravity and numerical orbit integration The method in Sec. 2.1 treats the nucleus as a point mass. That allows us to consider all ~400 million combinations of source regions (~2.5 ⋅ 104), vd values (~500), and target areas (~30) at a relatively low computational cost because the approach is essentially analytical. However, the nucleus is irregular and the gravity field deviates from that of a point source. We evaluate the error introduced by our point–mass assumption by comparing the Keplerian orbits with trajectories obtained by numerical integration of the equation of motion in a realistic gravity field. The nucleus gravitational potential U at the center of each nucleus facet, the gravitational force vector Fg(x, y, z) = ∇ U at any location exterior to the body surface, and the Laplacian ∇2U (taking the value 0 outside the surface of the nucleus and −4π within it) were calculated using the method of Werner and Scheeres (1997). This method applies to bodies with an arbitrarily complex shape described by a surface represented by polyhedrons under the assumption of constant density in the interior. Validation of the implementation was made by verifying that ∇U conformed with the Newtonian solution for polyhedrons forming a sphere, and that U calculated for the highly irregular shape of comet 67P reproduces the potential map of Keller et al. (2017) in their Fig. 4. In order to represent the shape of 67P we rely on the SHAP5 version 1.5 shape model (Jorda et al., 2016) degraded to 5 ⋅ 103 facets. In order to solve the equation of motion numerically an 8th order Gauss–Jackson integrator was implemented (Fox, 1984; Berry and Healy, 2004) using a 4th order Runge–Kutta method (Burden and Faires, 1993) as a startup formula. In order to validate the implementation the orbit of Comet C/1995 O1 (Hale–Bopp) was integrated from January 1993 through 2010 including the perturbations of the eight planets. This included a 0.77 AU encounter with Jupiter on April 5, 1996 that drastically modified the orbit, e.g., changed the orbital period from P = 4221 yr to 2372 yr. The solution was compared to that obtained with the RMVS3 version of the integrator SWIFT (Levison and Duncan, 1994) and the final positions of the comet according to the two codes differed by merely 70 m. (5) Solutions exist if f is real and f ∈ [−1, 1] (i.e., the target area does cross the orbital plane). The {rs, vs} pair can be used to calculate orbital elements for the dust trajectory in ecliptic space, following standard methods (e.g., Boulet, 1991). These orbital elements can be used to calculate the location rd(tpl) of the dust chunk at tpl. The distance d = ∣ rd(tpl) − rt(tpl)∣ is a function of vd only, for any given combination of source and target location. We calculate d as function of vd, starting at vd = 0 and incrementing in steps of 0.001 m s−1. Our condition for the existence of a Keplerian trajectory that connects the source and target regions is min(d) ≤ 25 m. It means that the dust trajectory crosses the circle swept up by the target area during nucleus rotation (co–location in space) and that the dust chunk and target area arrives to this interception point simultaneously (co–location in time). Co–location in space and time is a necessary but not sufficient criterion for airfall onto the target area. It is also required that the dust chunk does not intercept another part of the nucleus on its way to the target area. Therefore, all facets on the shape model that cross the dust trajectory orbital plane are located, their crossing times are determined and the dust chunk position at that time is calculated. If a co–location in space and time takes place prior to the dust chunk arrival time to the target area, the trajectory is discarded. The trajectory search was made by considering the SHAP5 version 1.5 shape model (Jorda et al., 2016) of comet 67P degraded to 5 ⋅ 104 facets. Out of the ~2.5 ⋅ 104 facets on the southern hemisphere treated as potential source regions, there were 151 facets that had the correct 6 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 vanishing temperature gradient; ⃒ ∂Tj ⃒⃒ = 0. ∂xj ⃒x=xmax 2.3. Nucleus thermophysics The search for uninterrupted dust trajectories described in Sec. 2.1 identified 151 unique source facets on the southern hemisphere responsible for airfall at the 31 target regions. The purpose of the thermophysical nucleus modeling is to calculate the local water production rates for those source regions as a function of nucleus rotational phase during a 315 day period from the inbound equinox on May 11, 2015, throughout the perihelion passage up to the outbound equinox on March 21, 2016. At this time the southern hemisphere was illuminated and the airfall in the north would have peaked. We assume that all regions have equal potential for activity when subjected to identical illumination sequences. This is consistent with the surface H2O activity distribution of Fougere et al. (2016b) that shows little variability within the southern hemisphere. All facets on the shape model that have an unobstructed view of at least one source facet were identified (here referred to as “surrounding terrain”). The ~3 ⋅ 104 facets in surrounding terrain have the capacity of shadowing the source facets by temporarily switching off the direct solar flux, and radiation that is scattered in the visual or emitted in the infrared from these facets will illuminate the sources (such nucleus self heating elevates illumination levels at all times and are particularly important when the sources are in shadow or experiencing night–time). We used the model presented by Davidsson and Rickman (2014) in order to calculate the illumination conditions, including shadowing and self heating. For a given source facet j we evaluate whether the Sun is fully visible (vj⊙ = 1) or if the facet is partially shadowed (rounded to vj⊙ = 1/ 3 or vj⊙ = 2/3) or fully shadowed by topography (vj⊙ = 0) in 10∘ increments of nucleus rotational phase throughout the time period, applying the spin axis in equatorial right ascension and declination {α, δ} = {69.54∘,64.11∘} of 67P (Preusker et al., 2015) and taking the time–dependent nucleus rotation period into account (e.g., Keller et al., 2015b). For simplicity we assume that the combined scattering and emission of a surrounding terrain facet k equals the local incident direct illumination (i.e., there is no loss or gain because of heat conduction from or to the surface). The fraction of the radiation emanating from k that reaches a source facet j is given by the view factor (e.g., Özişik, 1985) F jk = sk cosβj cosβk . πrjk2 Solving Eqs. (7) with its boundary conditions yields Tj = Tj(x, t) and the vapor production rate then is calculated as ) ( Qj = fi αs Z Tj (0, t) (10) where fi is the volume fraction of water ice, αs is the sublimation coefficient (Kossacki et al., 1999), ( ( )) 1 1 T − a2 π αs = 1 − + tanh − a3 tan π − , (11) a1 a1 273 − a2 2 where a1 = 2.342, a2 = 150.5, and a3 = 4.353, and the Hertz–Knudsen formula is √̅̅̅̅̅̅̅̅̅̅̅̅̅ m Z = psat (T) (12) 2πkB T where the saturation pressure of water vapor (in Pa) is given by Fanale and Salvail (1984) ( ) 6141.667 psat (T) = 3.56 × 1012 exp − . (13) T 2.4. Coma structure The purpose of this Section is to present the methods used for calculating the number density n(h), translational temperature T(h), and drift (expansion) speed W(h) as functions of height h above the nucleus surface (along the local surface normal). This calculation is non–trivial because outgassing from a solid surface into vacuum yields a vapor that initially is deviating drastically from thermodynamic equilibrium (e.g., Cercignani, 2000), i.e., the velocity distribution function does not resemble the Maxwell–Boltzmann function, and hydrodynamic relations as such do not apply because they are based on a simplified Boltzmann equation where the collision integral is zero. Molecular collisions will gradually evolve the gas toward thermodynamic equilibrium, but only if molecular collisions are common. Therefore, there is a kinetic region known as the Knudsen layer above the nucleus surface that may transit into a hydrodynamic coma at some height, unless the Knudsen layer extends to infinity. We need to deal with both cases, treating the first possibility (strong sublimation) in Sec. 2.4.1 and the second (weak sublimation) in Sec. 2.4.2. For various discussions about cometary Knudsen layers see, e.g., Skorov and Rickman, 1998, 1999; Huebner and Markiewicz, 2000; Crifo et al., 2002; Davidsson and Skorov, 2004; Davidsson, 2008; Davidsson et al., 2010. In the discussion that follows, we use the following subscripts to distinguish various regions; a sub–surface region responsible for feeding the coma “s”; the bottom of the Knudsen layer “0”; the top of the Knudsen layer “1”; a location within the hydrodynamic part of the coma “2”; and infinity “∞”. We start by defining the difference between strong and weak sublimation. At a given moment the nucleus outgassing rate and surface temperature of a given facet are given by {Qs, Ts}. The sub–surface reservoir number density capable of sustaining the production rate Qs is √̅̅̅̅̅̅̅̅̅̅̅̅̅ 2π (14) ns = Qs mkB Ts (6) Here, sk is the surface area of facet k, the angles between the surface normals of j and k and their connection line are βj and βk, respectively, and rjk is the distance between the two facets. For each source facet j we solve the energy conservation equation ( ) ∂T (x, t) ∂ ∂Tj (x, t) ρc s j = κ (7) ∂xj ∂t ∂xj for the temperature Tj by using the Finite Element Method. The surface boundary condition of Eq. (7) is given by ) ∑ (S⊙ vk⊙ Aμ (t) S⊙ vj⊙ (1 − A)μj (t) k + εσ Tk (0,t)4 = εσ Tj (0, t)4 − (1 − A) F jk 2 2 rh rh k∕ =j ⃒ ( ) ∂T ⃒ + fi αs Z Tj L − κ ⃒⃒ . ∂xj xj =0 (9) (8) From left to right, the terms in Eq. (8) denote the absorbed flux from direct solar illumination; thermally emitted infrared radiation; self heating in the form of scattered optical and thermally emitted radiation from other facets; energy consumption by sublimation of near–surface ice; and heat conducted from the surface to the interior or vice versa depending on the sign of the temperature gradient. The boundary condition of Eq. (7) at the bottom of the calculational domain (typically placed at ten times the diurnal thermal skin depth below the surface) is a and the Knudsen layer bottom number density is n0 = ) 1( 1 − αbf ns 2 (15) where αbf is the surface backflux fraction (i.e., the fraction of molecules that return to the surface). The molecular mean free path is 7 B.J.R. Davidsson et al. λ= Icarus 354 (2021) 114004 1 Using the ratios in Eq. (22) we apply ⎧ n1 − n0 ⎪ h + n0 n(h) = ⎪ ⎪ Kt ⎪ ⎪ ⎨ Qs W(h) = ⎪ ⎪ mn(h) ⎪ ⎪ ⎪ ⎩ T(h) = T1 (16) σc n0 where σ c is the molecular collisional cross section. The Knudsen layer is taken to have a thickness of Kt = 30λ (90–99% of thermodynamic equilibrium is achieved at 20–200λ according to Cercignani, 2000). For a nucleus curvature radius Rn we define strong sublimation as Rn/(Rn + 30λ) > 0.9, i.e., the Knudsen layer is relatively thin compared to the nucleus radius and the degree of gas expansion within the Knudsen layer is small. Such conditions are typical on the dayside at small heliocentric distances. We define weak sublimation as Rn/(Rn + 30λ) ≤ 0.9, i.e., the gas expansion is substantial within the Knudsen layer and downstream hydrodynamic conditions may not be reached. Such conditions are typical at large heliocentric distances, and on the nightside at any point in the orbit. for heights 0 ≤ h ≤ Kt above the nucleus surface, i.e., we assume a linear reduction of the number density from n0 near the surface (Eq. 15) to n1 at Kt, evaluate the drift speed according to mass conservation, and assume that the translational temperature remains quasi–constant throughout the Knudsen layer. The last property is motivated by numerical solutions to the Boltzmann equation obtained with the Direct Simulation Monte Carlo (DSMC) technique, see Davidsson (2008). For h > Kt we apply an analytical hydrodynamic solution valid for isentropic expansion summarized by Davidsson et al. (2010), see their Eq (1)–(6), that merges smoothly with the Knudsen layer solution. Specifically, we first find the Mach number M2 at height h by numerically solving 2.4.1. Strong sublimation We first describe our applied solution for the Knudsen layer, and then for the downstream hydrodynamic coma. Classical gas kinetic theory (Anisimov, 1968; Ytrehus, 1977) shows that the temperature and number density at the top of the Knudsen layer are related to the sub–surface parameters (the so–called “jump conditions”) through √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ )2 ⎧ ( T1 1 √̅̅̅ π ⎪ ⎪ S π + 1 + S21 = − ⎪ 1 ⎪ Ts ⎪ 8 64 ⎪ ⎨ √̅̅̅̅̅ (17) T1 ⎪ ⎪ G F + ⎪ ⎪ n T ⎪ 1 ⎪ ( s ) ⎩ = ns 2exp − S12 ⎧ ⎪ ⎪ ⎨ ( ( 2 ( ) ) 2 ⎪ 2 ⎪ ⎩ G = 2S1 + 1 erfc(S1 ) − √̅̅̅ exp − S1 C = Tn1−γ (18) π (20) − 1) ⎟ 1 ⎜ 1+ ⎝ ⎠ M2 1 + 12 (γ − 1)M22 . (24) (25) 2.4.2. Weak sublimation When collisions are rare, the gas is unable to achieve thermodynamic equilibrium. An accurate description of the coma in this case requires numerical solutions to the Boltzmann equation. Such sophisticated treatment of the problem is beyond the scope of this paper and we apply a simple parameterized description of the coma that reproduce general characteristics of DSMC solutions reasonable well (see Sec. 3). The average kinetic energy of a water molecule upon release from the nucleus is Furthermore, the analysis showed that the Mach number approached unity at the top of the Knudsen layer (M1 ≈ 1), which means that T1/Ts ≈ 0.67, n1/ns ≈ 0.21 and that β− implies αbf ≈ 0.18 (see Eq. 15). However, these results apply for monatomic vapor and are not directly applicable for water molecules that have ζ = 3 and γ = 4/3. Assuming that T1/Ts ≈ 0.67 still holds, requiring that M1 = 1 for γ = 4/3, recognizing that mass conservation yields W1 = Qs/mn1 on top of the Knudsen layer while applying Eqs. (14) and (19) yields the following jump condition for water vapor, T1 = 0.67 Ts √̅̅̅̅̅̅̅̅̅̅̅̅ ⎪ Ts ⎪ ⎪ n1 = ⎩ ≈ 0.42. ns 2π γT1 γ+1 ⎞−2(γ−1) With n2 = n(h) known, W2 = W(h) follows from Eq. (26) and T2 = T (h) from Eq. (25), thus the gas parameters for the strong sublimation case can be evaluated. The ratio of the number of molecules traveling back towards the surface at the bottom relative to the top of the Knudsen layer is )√̅̅̅̅ ( √̅̅̅ 2 2S12 + 1 TT1s − 2 π S1 − √̅̅̅̅ . (21) β = F + TT1s G ⎧ ⎪ ⎪ ⎪ ⎨ 1 (γ 2 Equations (25)–(26) are inserted into Eq. (19) which is solved for n, thereby allowing that parameter to be determined, ]−1 [ ( √̅̅̅̅̅̅̅̅̅̅̅̅ ) 12 (1−γ)−1 M2 γkB C n(h) = . (27) m D Here, the heat capacity ratio γ and the number of rotational degrees of freedom ζ are related by ζ+5 . ζ+3 = ⎛ is determined by using the known {n1, T1} values and the parameter D (h) is evaluated which yields the product of number density and drift speed at the corresponding height from mass conservation, ( )2 Qs Rn D = = nW. (26) m Rn + h where erfc() is the complimentary error function and the speed ratio √̅̅̅̅̅̅̅̅ is S1 = 5/6M1 where the Mach number is √̅̅̅̅̅̅̅̅̅̅ m M=W . (19) γkB T γ= h + Rn K t + Rn )2 Next, the constant C given by where ) ( √̅̅̅ F = − π S1 erfc(S1 ) + exp − S12 (23) 3 Et = kB Ts 2 (28) and its average internal energy because of molecular rotation is ζ Er = kB Ts 2 (29) (see, e.g., Bird, 1994). If collisions are very common (λ ≪ Rn) the far downstream molecular velocity vectors are radially aligned (i.e., essentially parallel to each other) and all speed variability has been removed through energy equipartitioning. Because there are no random velocities about the local mean, T∞, λ≪Rn = 0 and all kinetic energy is in bulk drift. Therefore, an ideal hydrodynamic gas that expands towards (22) 8 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 infinity (n → 0) has T → 0 (so that Eq. 25 remains valid at arbitrarily low n). Because expansion leads to translational temperature cooling, and because collisions are common, energy is transferred from the internal to the kinetic mode and the rotational temperature cools to Tr = 0 at infinity as well. Thus, all available energy is in bulk drift kinetic energy mW∞, λ≪Rn2/2 = Et + Er and therefore √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ (3 + ζ)kB Ts W∞,λ≪Rn = . (30) m liberates dust particles from their attachment to the nucleus does not transfer momentum to the dust. For the current application we do not need to specify the mechanism that is responsible for the chunk release. However, we note that in laboratory experiments Ratke and Kochan (1989) observed oscillation of dust aggregates at a frequency of 1–100 Hz that lasted for up to 10 minutes before their final release. This suggests that gas drag ultimately is responsible not only for the acceleration of dust but also for its liberation, through a fatigue process that gradually reduces the initially strong cohesion of the granular nucleus material to the point the chunk breaks free. We evaluate the local gravitational force Fg(h) (directed towards the nucleus) at height h by using the Werner and Scheeres (1997) formalism discussed in Sec. 2.2. Therefore, the point–mass assumption is used to determine the orbital speed vd that is necessary to reach the target site, but realistic gravity is used when calculating the initial acceleration from rest to the point that vd is reached. This approach is reasonable, given that realistic gravity rapidly approaches the point–mass solution with height. We evaluate the gas drag force Fdrag(h, ud) by utilizing the local coma solution {n(h), T (h), W(h)} obtained in Sec. 2.4 for chunks of cross section A and applying (e.g., Gombosi, 1994), ⎧ 1 ⎪ ⎪ Fdrag (h, ud ) = A CD (h, ud )n(h)(W(h) − ud (h) )2 ⎪ ⎨ 2 ( ) ( ) ( ) 2 ⎪ 2W 4W 4 + 4W 2 − 1 erf (W ) + 1 exp −W 2 ⎪ ⎪ . + √̅̅̅ ⎩ CD (h, ud ) = 4 3 2W πW If collisions are very few (λ ≫ Rn) the downstream velocity vectors are still radially aligned purely because of geometric reasons. However, a molecular speed dispersion remains that collisions have not been able to remove. Therefore, the downstream flow is characterized by a non–zero translational freeze–out temperature T∞,λ≫Rn = Tfr = ∕ 0. Because most translational energy still is in bulk drift kinetic energy (Tfr is low) and because very little rotational energy has been transferred (Tr re2 mains quasi–constant) the bulk drift kinetic energy is mW∞,λ≫R /2 = Et n and therefore √̅̅̅̅̅̅̅̅̅̅̅̅ 3kB Ts W∞,λ≫Rn = . (31) m For the transition regime 0 < λ < Rn we assume that there is a linear transition from {T∞, W∞}λ≪Rn to {T∞, W∞}λ≫Rn. The temperature and drift speed change with height, starting at their near–surface values {T0, W0} and gradually approaching their downstream values {T∞, W∞}. The rate of change depends on the number of molecular collisions that take place per time unit. That number depends on the local gas density, here assumed to follow n ∝ h−2. We realize these assumptions through the expressions ⎧ [ ) ( )2 ](√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ⎪ ⎪ Rn (3 + ζmax{1 − λ0 /Rn ,0})kB Ts ⎪ ⎪ W(h) = 1 − + W0 − W 0 ⎪ ⎪ Rn + h m ⎪ ⎪ ⎪ ⎪ ⎪ [ ⎨ ( )2 ] ( ) Rn Tfr − Tfr max{1 − λ0 /Rn ,0} − T0 + T0 T(h) = 1 − ⎪ ⎪ Rn + h ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ D (h) ⎪ ⎪ ⎩ n(h) = W(h) (34) Here, erf() is the error function and W is given by √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ m W = (W(h) − ud (h) ) . 2kB T(h) Equation (33) is solved numerically for every source region, for every 10∘ advancement of the nucleus rotational phase during the considered ~11 month period around the perihelion passage, for 19 chunk size classes (with dimensions that are 0.05A crit –0.95A crit in increments of 0.05A crit , where A crit is the critical maximum liftable chunk cross section for which Fg(0) + Fdrag(0, 0) = 0 is valid at the considered location and time). We here consider chunks shaped as oblate ellipsoids with axis ratio 4, exposing their largest cross section to the gas flow at h = 0 (~2.5 times that of an equal–mass sphere), that turn with constant radial velocity to expose their smallest cross section (~0.63 times that of an equal–mass sphere) during the first 500 m of flight. The chunk shape is close to the axis ratio of 5 needed in order to match the particle speeds measured by Rosetta/GIADA according to Ivanovski et al. (2017a). The time needed for the chunks to reach the height of 500 m is similar to that needed by the aerodynamic torque to set non–spherical chunks in rotation according to simulations by Ivanovski et al. (2017b). The assumed decrease of the chunk cross section with time shortens the time needed for the particle to reach its terminal velocity, and makes Fdrag diminish in strength faster than Fg, thereby illustrating how the assumed dominance of nucleus gravity over drag (after the initial acceleration) can be motivated. We select the chunk size acrit(t) that reaches a quasi–terminal velocity of ud, max that is as close to vd as possible. This is considered the best representative of the chunk type that feed a particular target area at a given time. Note that this size changes continuously during nucleus rotation and orbital motion. (32) where λ0 = (n0σ )−1 and D (h) is evaluated as in Eq. (26). Inspired by the numerical DSMC solutions to the Boltzmann equation by Davidsson et al. (2010) we set Tfr = 20 K. The expression for W(h) in Eq. (32) can be understood as follows. The first term within the curved bracket is the downstream terminal drift speed. It is close to Eq. (30) when λ0 ≈ 0 and falls linearly with λ0 towards Eq. (31) as λ0 → Rn, and remains at that value for even more diluted gas (λ0 ≥ Rn) because of the maximum function. The expression within the curved bracket therefore measures the largest considered increase above W0, the near–surface value. The expression within the square bracket (confined to the [0, 1] interval) regulates how quickly W(h) grows from W0 at h = 0 to W∞ as h → ∞. It is an inverse–square law of height, because of the previously mentioned relations between collision frequency, rate of change of the kinetic properties, and number density. The expression for T(h) follows a similar logic. 2.5. Dust ejection 2.6. Airfall In order to calculate the velocity reached by a dust chunk of a particular size we need to solve the equation of motion along the local surface normal, md d2 h = Fg (h) + Fdrag (h, ud ) d t2 (35) In order to estimate the amount of mass that is channeled toward a certain target region it is necessary to know the total amount of dust that is being ejected from a source region at a given time, and the fraction of that mass carried by chunks of size acrit(t). The total dust mass flux is taken as Qd = 4Qs (Rotundi et al., 2015). The estimate of the fraction of this mass carried by acrit(t)–chunks requires that we specify a dust size distribution function. (33) where the chunk mass is md and the chunk speed is ud = dh/dt. We assume that the initial velocity is zero, i.e., that the mechanism that 9 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Table 4 Target area locations and thickness of accumulated layers formed by fallback of material released from the southern hemisphere. Target Lat Long #1 #2 #3 #4 #5 #6 #7 #8 #9 #10 #11 #12 #13 #14 #15 #16 40.9 34.6∘ 30.3∘ 27.4∘ 19.8∘ 9.1∘ 19.0∘ 16.4∘ 21.2∘ 13.9∘ 6.4∘ 15.3∘ 48.87∘ 56.5∘ 45.7∘ 58.9∘ ∘ 338.2 320.2∘ 332.9∘ 340.2∘ 346.0∘ 349.2∘ 339.9∘ 336.0∘ 331.7∘ 330.6∘ 333.5∘ 315.9∘ 108.0∘ 120.9∘ 135.6∘ 142.4∘ ∘ Thickness [m] Target Lat 0.86 0.77 0.24 0.40 1.82 1.48 0.19 0.00 0.62 0.00 0.27 0.41 3.64 1.30 1.80 0.41 #17 #18 #19 #20 #21 #22 #23 #24 #25 #26 #27 #28 #29 #30 #31 49.6 58.1∘ 49.2∘ 42.5∘ 33.9∘ 30.3∘ 25.4∘ 24.4∘ 19.1∘ 17.6∘ 19.2∘ 16.8∘ 11.3∘ −10.7∘ −7.2∘ For 3.7 mm ≤ acrit ≤ 1.62 cm we use the size distribution for chunks observed in Philae/CIVA images at Abydos, consisting of three segments with different slopes (Poulet et al., 2017). We extrapolate the ≤5 mm part of this distribution to smaller sizes. For 4 cm ≤ acrit ≤ 1 m we use the size distribution for chunks observed in Philae/ROLIS images at Agilkia (Mottola et al., 2015). We use this distribution to bridge the 1.62–4 cm gap in the measured data, and to extrapolate to larger sizes. With a differential size distribution (dN/dD ∝ D−qs) slope of qs = 3.8 (Mottola et al., 2015) truncated at Dmin = 0.0162 m and Dmax = 1 m the amount of mass at diameters D ≥ Dc = 0.1 m is given by ∫ Dmax 3−q 4−qs s D s dD Dmax − D4−q c F = ∫DDcmax = 4−qs (36) 4−qs 3−q D s dD Dmax − Dmin D Long Thickness [m] 158.1 171.1∘ 177.1∘ 180.5∘ 203.1∘ 200.9∘ 206.4∘ 198.8∘ 195.9∘ 202.9∘ 206.1∘ 212.1∘ 210.1∘ 118.6∘ 146.1∘ ∘ ∘ 0.77 2.95 1.44 1.22 0.87 0.42 0.79 0.04 0.85 0.51 0.00 0.51 2.07 0.11 0.09 model but here state and briefly discuss the governing equations. The energy conservation equation is given by ( ) ( ) ∂T 1 ∂ ∂T 1 ∂ sinl ∂T ρc(T) = 2 κ(ψ , T)r2 + κ(ψ , T) rsinl ∂l r ∂l ∂t r ∂r ∂r − nsp ∑ ı=4 − qı (pı , T)L ı (T) + nsp nsp ∑ ∑ ( ı=2 ȷ=5 ′ qı (T){Hı − Fıȷ L ȷ (T) } ( ) nsp ∑ ∂T Ψı ∂T gı Φı − + R* r ∂l ∂r ı=4 ) (37) with the following terms going left to right; 1) change of the internal energy; 2) radial heat conduction; 3) latitudinal heat conduction; 4) sublimation of ice and recondensation of vapor; 5) energy release during crystallization of amorphous water ice and energy consumption during its release of occluded CO and CO2, as well as energy consumption during CO segregation from its partial entrapment in condensed CO2; 6) advection during radial and latitudinal gas diffusion; 7) heating by radioactive decay. The indices in Eq. (37) refer to refractories (ı = 1), amorphous water ice (ı = 2), cubic water ice (ı = 3), hexagonal (crystalline) water ice (ı = 4), carbon monoxide (ı = 5), and carbon dioxide (ı = 6). Species denoted by ȷ are hosted within a more abundant and less volatile species denoted by ı. The upper boundary condition to Eq. (37) is given by ⃒ S⊙ (1 − A)μ(l, t) ∂T ⃒ 4 , (38) = σεTsurf − κ ⃒⃒ 2 ∂r r=Rn rh min and amounts to F = 66%. If this slope is extrapolated to Dmin = 1 mm then F = 80% for Dc = 1 cm and F = 49% for Dc = 0.1 m. It is clear that a majority of the mass is carried by chunks that are cm–sized or larger. We truncate the size distribution function at the maximum ejectable size at a given source location and time, normalize the size distribution so that it carries the entire mass being ejected, and determine the fraction of mass carried by chunks with sizes in the range acrit(t)/2–2acrit(t). The total airfall on a given target area is calculated by integrating the mass over all contributing source regions and time. We consider the activity of most target areas to be so low during the mass accumulation period that we assume the airfall has free access to the surface. The exception is target areas #30 and #31 (see Table 4) that themselves are located on the southern hemisphere and are strongly illuminated at times. We therefore calculate the maximum liftable size in those locations versus time and only accept airfall chunks that are larger (i.e., we assume that the smaller chunks will be re–ejected within a short period of time if they manage to reach the surface in the first place). Therefore, sites #30 and #31 in Imhotep are partially self–cleaning. The thickness of the accumulated layer at each site is calculated from the incident mass flux assuming that chunks build a deposit with the same bulk density (ρ = 535 kg m−3) as the nucleus (Preusker et al., 2015). and balances absorbed solar radiation with net thermal radiation losses and heat conduction. The boundary condition at the body center is a vanishing temperature gradient. The mass conservation equation for vapor (ı ≥ 4) is given by ψ mı nsp ∑ ∂nı 1 ∂( 2 ) 1 ∂ ′ Fȷı qȷ (pȷ , T), r Φı − = − 2 (Ψı sinl) + qı (pi , T) + r ∂r rsinl ∂l ∂t ȷ=2 (39) with terms left to right; 1) changes to the gas density; 2) radial gas diffusion; 3) latitudinal gas diffusion; 4) release during sublimation and consumption during condensation; 5) release of CO and CO2 during crystallization of amorphous water ice, and of CO from sublimating CO2 ice. The boundary conditions to Eq. (39) consist of a quantification of the gas venting to space at the surface and a vanishing gas flux at the center. Finally, the mass conservation equation for ice (ı ≥ 2) is given by 2.7. Thermophysics of coma chunks During their flight through the coma the airfall chunks will lose some of the ice they carry. In order to calculate the typical volatile loss from chunks and to estimate the ice abundance in fresh airfall deposits we apply the code NIMBUS (Numerical Icy Minor Body evolUtion Simulator) developed by Davidsson (2020). NIMBUS considers a rotating spherical body consisting of a porous mixture of dust and ice, and tracks the internal ice sublimation, vapor condensation, and diffusion of gas and heat in the radial and latitudinal directions over time during body rotation. We refer to Davidsson (2020) for a detailed description of the ∂ρı = −qı (pı , T) + τı (T), ∂t (40) with terms left to right; 1) changes to the amount of ice; 2) changes due to sublimation of ice and recondensation of vapor; 3) changes due to 10 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Fig. 2. Surface temperature (upper panel) and water production rate (lower panel) versus time for one of the source regions during a single nucleus rotation near inbound equinox. The circles denote minimum and maximum levels of activity. Fig. 3. Number density, translational temperature and drift speed at T = 211 K. other phase transitions (crystallization of amorphous water ice, transition from cubic to hexagonal water ice). In the current simulations we nominally only consider one volatile – hexagonal (crystalline) water ice that initially constitutes 20% of the body mass. However, one simulation also has condensed CO2 at a 5% level relative to water by number. The specific heat capacity c(T) and the heat conductivity κ(ψ , T) are mass–weighted averages of temperature–dependent values measured in the laboratory for water ice (Klinger, 1980, 1981), carbon dioxide (Giauque and Egan, 1937; Kührt, 1984) and refractories (forsterite from Robie et al., 1982 for c and ordinary chondrites from Yomogida and Matsui, 1983 for κ). The heat conductivity is corrected for porosity following Shoshany et al. (2002) and includes both solid state and radiative conduction. The volume mass production rates qı are standard expressions (e.g., Mekler et al., 1990; Prialnik, 1992; Tancredi et al., 1994) with H2O and CO2 saturation pressures taken from Huebner et al. (2006) and the gas diffusion fluxes (Φı , Ψı ) are evaluated as in Davidsson and Skorov (2002). The radiogenic heating from a number of long–lived isotopes at chondritic abundances in the refractory component is included by default but has a completely negligible influence on the results in the current application. 11 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Fig. 4. Number density, translational temperature and drift speed at T = 129 K. Fig. 5. Velocity versus height for chunks with about the right size to reach a targeted terminal velocity vd (red) needed to take them to a given target region. The upper panel shows conditions near midnight while the lower panel shows conditions near noon. 2.8. Longterm loss of volatiles 3. Results In order to better understand the fate of water ice after it has been deposited as airfall we perform thermophysical simulations for a selection of target areas. These simulations rely on the same model as described in Sec. 2.3. The major difference is that the illumination conditions at the target areas (calculated for every 10∘ of nucleus rotation) are not restricted to the orbital arc between the inbound and outbound equinoxes, but are made for the entire orbit. 3.1. The thickness of airfall debris layers In order to exemplify the application of the models presented in Sec. 2 we here describe sample results. The thermophysical modeling of the nucleus (Sec. 2.3) provides surface temperature Tsurf(t) and water production rate Qs during nucleus rotation. Fig. 2 shows an example of diurnal temperature and production curves for a southern hemisphere 12 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Fig. 6. Ejected chunk size versus time during one nucleus rotation near perihelion for a number of source regions (shown with different symbols) that all feed the same target area. Fig. 7. Ecliptic Cartesian position coordinates versus time for a Keplerian trajectory connecting a source region with target area #24 for an initial speed vd = 0.855 m s−1 along the local surface normal in addition a 0.081 m s−1 component due to nucleus rotation (dashed-dotted). The solid curves correspond to a numerical orbit integration using a realistic gravity field of a rotating nucleus for vd = 0.855 m s−1 (left) and vd = 0.800 m s−1 (right). source region1 that contribute with material to target area #3, at a time near the inbound equinox. In this particular case the surface temperature peaks near Tsurf = 211 K at local noon and falls to T = 129 K just prior to sunrise. The waviness of the temperature curve is a consequence of changes in shadowing and self heating conditions during nucleus rotation, in addition to the continuous changes of the solar height above the local horizon. The gas production rate is a strongly non–linear function of surface temperature (see Eqs. 12–13). Therefore, it varies by seven orders of magnitude at this particular location and time. This has strong implications for the properties of the near–nucleus coma. Fig. 3 shows the number density, translational temperature, and expansion velocity of the gas emanating from the same region as in Fig. 2 near noon (see the right circle in that figure). Due to the high gas production rate, this case illustrates the strong sublimation case described in Sec. 2.4.1, i.e., a Knudsen layer separates the nucleus and the hydrodynamic coma. Here, the Knudsen layer is ~2.15 m thick. The initial gas drift speed is ~303 m s−1 and the gas has an initial translational temperature of ~141 K, much lower than the surface temperature because of the kinetic jump condition (Eq. 22). As the gas recedes from the nucleus surface, the expansion results in a number density reduction which is nearly an order of magnitude at 1 km from the surface. The gas 1 This randomly selected source is located in the Seth region at longitude 105.5∘ W, latitude 21.1∘ S, at a distance 1.33 km from the nucleus center (origin of the Cheops coordinate system, see Preusker et al. 2015). 13 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 acceleration to ~605 m s−1 is a consequence of mass conservation (a thinner gas needs to flow faster to carry the same amount of mass over a boundary). The adiabatic cooling of the gas is evident and the translational temperature drops to ~66 K within the lower kilometer of the coma. The conditions in the nighttime coma are drastically different. Fig. 4 shows the number density, translational temperature, and expansion velocity of the gas emanating from the same region as in Fig. 2 just before dawn (see the left circle in that figure). Because of the low production rate, Fig. 4 illustrates the weak sublimation case described in Sec. 2.4.2. The coma never achieves thermodynamic equilibrium and has to be modeled as an infinite Knudsen layer. Over the inner 1 km of the coma, the gas translational temperature drops from ~87 K to ~34 K while the expansion velocity increases from ~238 m s−1 to ~385 m s−1. The number density is ~7 orders of magnitude lower than on the dayside, as expected. The source region in question has such a position, local rotational velocity in the inertial frame, and surface normal that it would be capable of launching chunks to a particular target area on the opposite nucleus hemisphere if these chunks are accelerated to vd ≈ 0.63 m s−1 in the vicinity of the nucleus. For comparison, the escape velocity of 67P √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ is vesc = 2GN Mn /Rn = 0.83 m s−1 , where the nucleus mass is Mn = 9.982 ⋅ 1012 kg (Pätzold et al., 2016) and the area–equivalent nucleus radius is Rn = 1932 m (Jorda et al., 2016). The size of the chunks that will reach vd ≈ 0.63 m s−1 (see Sec. 2.5) varies drastically during the day. Fig. 5 shows results from the numerical integration of Eq. (33) for both the noon and near–dawn conditions discussed previously. At noon, the gas flow is sufficiently strong to lift chunks with diameters just under 1.9 m. However, those chunks are moving too slowly to reach the particular target area in question. Chunks with ~0.3 m diameter do reach the desired velocity (Fig. 5, lower panel) and are expected to reach the target area, while chunks smaller than this will travel too fast. When the gas production rate wanes after noon, the critical chunk radius acrit capable of entering the correct transport route to the target area diminishes rapidly. Fig. 5 (upper panel) shows that only nano–sized grains would reach the desired velocity, if such particles exist. Throughout the orbit the capability of this and other source regions to launch large chunks towards the northern hemisphere changes significantly. Seen from many source regions, the Sun is circumpolar at perihelion and will continuously send large chunks to the north. For the region discussed above the solar radiation intensity actually diminishes because of its location on the nucleus. The surface temperature fluctuates between 147–206 K and gas release can eject chunks with up to 1.3 m diameter range. It would launch chunks with vd ≈ 0.63 m s−1 for diameters in the 1.4 μm–0.2 m range. These numbers serve the purpose of illustrating the drastic variations in airfall chunk sizes a given target region may experience in just a few months during perihelion approach, because of changing illumination conditions at the source during nucleus rotation and orbital motion. The substantial variability of the critical chunk size during the comet day is shown more clearly in Fig. 6. Here, acrit(t) is shown for a number of different source regions that all contribute to target area #3 (see Table 4), for a single nucleus rotation near perihelion. One source is continuously feeding cm–to dm–sized chunks to the target area because the Sun is circumpolar at that location. A couple of other sources experience day/night variations and their contributions fluctuate between micron and ~0.1 mm contributions at night to cm–to dm–sized chunks in the day. Other sources are poorly illuminated and are only capable of providing sub–micron particles. The material feeding a given target area therefore originates from a variety of locations and arrives in chunks of drastically different size. Before proceeding to the main results of this paper we first discuss the accuracy of the orbital solutions for the chunks. The orbit calculations in this paper are made assuming a nucleus point mass and corresponding Keplerian orbits (Sec. 2.1). We now compare such solutions to numerically integrated orbits in the realistic gravity field of a rotating irregular nucleus (Sec. 2.2). The dashed-dotted curves in Fig. 7 shows the {x(t), y(t), z(t)} coordinates in the ecliptic system of a Keplerian orbit that connected a particular source region with target area #24. In this particular case, airfall at that location could be realized for vd = 0.855 m s−1 according to the idealized point mass solution and a time of flight of 9.7 h. The left panel of Fig. 7 shows that the numerically integrated trajectory in a realistic gravity field rapidly diverges from the idealized one and after 9.7 h the difference amounts to 7.7 km. The main reason for this drastic difference is that the nucleus potential at the launch site in question is merely ~90% of that obtained if Fig. 8. The amount of accumulated airfall material expressed as layer thickness in meters for target areas in Serqet (blue), Ash (red and green), and Imhotep (black). 14 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 all nucleus mass was concentrated at the origin. The chunk injects into an orbit that perhaps would have been similar to that around a point mass carrying ~10% less mass than 67P. If we consider the vis–viva equation (e.g., Danby, 1989) that relates the orbital speed vo and distance ro from the origin for any elliptical orbit with semi–major axis ao, ) ( 2 1 v2o = μo (41) − ro ao among real chunks because they ought to have some density dispersion, we do not attempt such a correction but assume that our acrit–values are representative (particularly considering the factor 4 margin used in Sec. 2.6 when evaluating the fraction of the ejected mass that contribute to airfall). The final results of our simulations regarding airfall accumulation are tabulated in Table 4 and are shown graphically in Fig. 8. The thicknesses of accumulated layers during a perihelion passage vary rather substantially across the nucleus surface. Some regions collect only small amounts of material, amounting to a few centimeters or less. Most regions accumulate a few decimeters. However, some regions (particularly #5–6, 13–15, 18–20, and 29) collect substantial layers that are 1.5–3.6 m thick. The average thickness for all 31 target areas is 0.87 m. Target areas #7, 8, and 10 are the ones located closest to Agilkia, where Philae bounced. The average airfall layer thickness for those target sites is 0.27 m according to our investigation. We therefore predict that the thickness of the airfall layer is comparably thin at this part of the nucleus (about a third of the average thickness) and we note that Biele et al. (2015) estimate a thickness of ≳0.20 m for the “granular soft surface” that possibly is located “on top of a more rigid layer”. We are encouraged by the similarity between the layer thickness inferred from the Philae measurements, and the estimate based on our calculations. Target areas #1–12 are located in the Ma’at and Serqet regions on we realize that any reduction of the reduced mass μo to μo* < μo at a fixed ro could be compensated for by lowering the velocity by a factor √̅̅̅̅̅̅̅̅̅̅̅̅ μ*o /μo , thereby achieving exactly the same semi–major axis. Taking the velocity component due to nucleus rotation into account, this suggests that if the chunk velocity due to gas drag was reduced by ~5% compared to its original value (i.e., from 0.855 m s−1 to 0.800 m s−1) we would expect the numerically integrated chunk trajectory to follow the Keplerian one. The right panel of Fig. 7 shows a numerical orbit integration when vd = 0.8 m s−1 and the Keplerian orbit is indeed reproduced very well. It therefore seems like very modest modifications to the ejection speed could compensate for differences between the real and the idealized gravity field. In principle, an adjustment of the ejection velocity should be compensated for by a corresponding change of acrit. However, considering that the cross–section to mass ratio will vary Fig. 9. Properties of a D = 1 cm coma chunk after 11.4 h in the coma. Upper left: porosity. Upper right: temperature [K]. Lower left: H2O vapor pressure [Pa]. Lower right: ice abundance versus distance from the center at chunk latitude 25∘ S, normalized to the initial local value. The location of the displayed region is shown as a solid black line in the upper right panel. 15 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 the small lobe, extending from the border of Hatmehit towards the Hathor cliff. Collectively, these areas accumulate an average of 0.6 m, with the largest deposits (#5–6) found at the Hatmehit border. Target areas #13–20 are located in the eastern Ash region on the large lobe and extend from the border with Aten westward up to the antimeridian. As a group they have the largest average layer thickness of 1.7 m with the thinnest layers (#16 and #17) located at the very center. Target areas #21–29 are also located in Ash but in its western part stretching from the antimeridian up to the borders to Seth, Anubis, and Aten. The average layer thickness is 0.7 m for this group with the largest accumulation (#29) near those borders. Finally, target areas #30-31 are located in Imhotep, just south of the equatorial plane within the vast smooth region that dominates this morphological unit. Because of efficient self–cleaning the net accumulation here is the lowest, 0.1 m. the depth does not vary significantly with latitude because of the tumbling of the chunk that exposes all parts of the surface to strong solar illumination over time. At the particular snapshot shown in Fig. 9 the subsolar point has been at high northern latitudes sufficiently long for the south pole to cool off, and H2O vapor has recondensed within the dust mantle. This is hinted at by a slight decrease of porosity in the middle of the southern hemisphere dust mantle, and is seen more clearly in the lower right panel showing the ice abundance (normalized to the initial amounts) versus radial distance at latitude 25∘ S. At depth the normalized ice abundance is slightly above unity because of downward diffusion and recondensation of vapor. Near the surface the ice abundance has previously gone to zero in the dust mantle, and recently been slightly elevated again, because of the previously mentioned near–surface cooling and associated vapor recondensation. The subsolar point is moving south which means that this frost is starting to sublimate at latitudes near the equator. That gives rise to the local vapor pressure maximum seen in the lower left panel. In order to explore the robustness of our results, we made three tests where the nominal parameters in the D = 0.1 m model were changed. First, the radius of constituent grains was changed from rg = 1 μm to 100 μm, which decreases the volume sublimation rate by two orders of magnitude (q ∝ r−1 g ). The ice loss as well as internal temperatures and vapor pressures changed insignificantly. This can be understood as follows. During one second the part of the difference between absorbed flux and thermally emitted flux that is not used to heat sub–surface material must necessarily be consumed by net sublimation of ice. If an amount of energy ΔE is available, a mass Δm = ΔE/L (T) must be converted to vapor (where the latent heat is evaluated at the temperature of the sublimation front). The temperature and pressure gradients on both sides of the sublimation front will evolve in such a way that gas diffusion removes Δm from the sublimation front region. This flow gives rise to continuous reductions of the vapor pressure pv below the local saturation pressure psat, which triggers sublimation because q ∝ psat(T) − pv. The time–scale for re–establishing saturation conditions after a temporary deviation is extremely short (micro–seconds or less) because q ≫ 0 when pv < psat (and q rapidly goes to zero when pv → psat). Therefore, even if the response time to pressure changes is increased two orders of magnitude because of rg, this time is still much shorter than the characteristic time scales of gas diffusion and heat conduction, thereby making the model insensitive to the size of the icy particles. In another test the dimensions of the cylindrical tubes used to evaluate gas diffusivity were increased by two orders of magnitude from radius 1 μm to 0.1 mm and from length 10 μm to 1 mm. Note that the pore radius also affects the radiative component of heat conduction, but that parameter was still kept at 1 μm in order to isolate the dependence of the ice loss on the gas diffusion modeling. In this case, the water ice loss increased from 6.4% to 8.9%. When the gas diffusivity increases drastically because of the longer and wider pores, the main effect is a reduction of temperature and vapor pressure, in such a way that the resulting gas fluxes inward and outward are essentially the same as before. The same net sublimation Δm is required in order to consume the discrepancy between absorbed and emitted radiation that does not change much, and that goal can be met for a cooler and less pressurized gas. For example, a cell that had T = 204.4 K with the stricter tubes cooled by 15.4 K to T = 189.0 K when the tubes widened, and the vapor pressure fell from pv = 0.2234 Pa to pv = 0.0055 Pa. The sub–surface cooling slightly reduces the surface temperatures as well, meaning that less energy is lost radiatively to space and more energy is available for sublimation. That is why the ice loss becomes somewhat larger, but this increase is still modest. In the final test, the pore radii used to regulate the radiative contribution to heat conduction was increased by two orders of magnitude as well. This is expected to have a stronger effect on the ice loss because the dust mantle becomes rather warm (325 K) compared to the icy interior, and wider pores allow for the thermal radiation to penetrate more efficiently to the ice. In fact, the total ice loss increased from 6.4% to 3.2. Volatile loss during the transfer through the coma We now proceed to the question of the ice abundance in airfall deposits. Observations of dust jet switch–off after local sunset shows that the water ice sublimation front is located just a few millimeters below the nucleus surface (Shi et al., 2016). The process that liberates and ejects chunks in the cm–m class, including but not limited to various forms of cracking (El-Maarry et al., 2015a; Auger et al., 2018) and cliff collapse (Pajola et al., 2017a), combined with gas drag, should produce coma chunks with an initial ice abundance that may be similar to that of the bulk nucleus. However, some of this ice will be lost during the flight through the coma, thereby lessening the initial ice abundance in airfall deposits. In order to estimate the ice abundance in fresh airfall deposits, we used the code NIMBUS described in Sec. 2.7 to study a D = 0.1 m chunk in the coma, rotating about a fixed axis with coordinates {α, δ} = {0∘, 45∘} in the equatorial system with a period P = 20 min (individual chunks with similar sizes and rotation periods have been observed in the 67P coma, see Davidsson et al., 2015). These simulations assumed a porosity of 70%, micro–meter sized grains (influencing the volume sublimation rate q) and pore lengths and radii of 10 μm and 1 μm, respectively (influencing the gas diffusivities Φ and Ψ as well as the magnitude of radiative heat transfer). The chunk was resolved by 18 latitudinal segments and 69 cells in the radial direction using geometric progression such that surface cells were 0.5 mm thick and core cells were 1 mm thick. The diurnal skin depth was resolved by ~10 cells. After a 12 h flight this body lost a total of 4.8% of its original ice content. Because parts of the chunk received little radiation due to its spin axis orientation, and because small bodies of this type may not maintain a fixed spin axis, we then introduced a constant precession rate of α from 0∘ to 360∘ in one hour, combined with a declination nutation with a full oscillation from δ = 90∘ to −90∘ and back again during the same amount of time. Because of the more even distribution of illumination, the volatile loss increased marginally to 6.4% of the initial ice content. Next, a D = 0.01 m chunk was considered. The rotation period was reduced to 2 min but the same precession and nutation periods were applied as before. Because of the smaller size of this chunk its ice loss was more substantial. Yet after 3 h, only 37% of the ice had been lost, increasing to 45% loss after 6 h and to 56% loss after 12 h. This illustrates that the ice–loss rate slows down rapidly once the sublimation front has moved a few millimeters below the surface, and that chunks as small as ~1 cm are capable of retaining a substantial fraction of their water ice on time–scales similar to the coma flight time. Coupled with the observation that the majority of mass in smooth region deposits is bound in D≳1 cm chunks (with at least half of that mass in D≳1 cm chunks) this means that fresh northern airfall deposits are nearly as icy as the southern source regions that produced them. Fig. 9 shows the internal distributions of porosity, temperature, and H2O gas pressure for a D = 1 cm chunk about 11.4 h into the simulations. The sublimation front has withdrawn to a depth of ~1 mm and 16 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Fig. 10. Diurnal minimum and maximum temperatures as functions of time and heliocentric distance throughout the orbit of 67P at the surface and at 0.2 m below the surface for four target areas. The gray fields indicate the time period from the inbound equinox to the outbound equinox. Note that the red solid and dashed curves largely overlap. The panels to the right were first published, in a slightly different format, by Takigawa et al. (2019). 13.7%. Yet, this increase is not large enough to alter our overall conclusions – coma chunks preserve water ice rather efficiently. We do not expect these numbers to be strongly dependent on the dust–to–ice mass ratio. The water ice sublimation front is expected to withdraw to similar depths regardless of water abundance. At that critical depth, the cooling by net sublimation is reduced by the dust mantle quenching of gas diffusion to the point that the hot dust mantle dissipates the majority of absorbed energy through thermal reradiation. The fractional volatile losses in terms of volume and mass are thus expected to be similar. We now consider the presence and survival of CO2 ice in coma chunks. Davidsson et al. (2020) reproduced the observed pre–perihelion CO2 production rate curve of 67P from 3.5 AU to 1.24 AU with a NIMBUS model that had the perihelion CO2 sublimation front located ~1.9 m below the surface near the equator, and where that depth diminished with latitude to ~0.2 m at the south pole. If the suggested presence of CO2 ice at shallow depths is correct, it is conceivable that regular sublimation–driven activity occasionally can expel extremely cold chunks containing supervolatiles at perihelion. The presence of CO2 ice patches on the surface of 67P observed spectroscopically by Filacchione et al. (2016) could be the result of such temporary exposure of near–surface CO2 ice. There are other forms of activity that might produce CO2–rich chunks as well. For example, the cliff collapse in Aswan observed by OSIRIS suddenly exposed material located ~12 m below the previous surface (Pajola et al., 2017a), and though such events are rare, they could occasionally inject unusually cold chunks into the coma that still contain CO2. It is interesting and important to understand to what extent such material could be transported and mixed in with other airfall material. Therefore, in our last numerical experiment, 5% condensed CO2 relative to H2O (by number) was added to the chunk, and the initial temperature was lowered from T(t = 0) = 150 K to T(t = 0) = 50 K to avoid an explosive sublimation of the carbon dioxide. Our NIMBUS simulations of a D = 0.1 m chunk shows that it takes 2 h for all CO2 to be lost. A substantially larger chunk could potentially preserve a fraction of its CO2 ice during the transfer to the northern hemisphere. After a total of 12 h exposure to the Sun at the perihelion distance of 67P, the amount of water ice lost in the D = 0.1 m chunk is 5.4%, i.e., somewhat lower than in the simulation without CO2. This reduction happens for two reasons; 1) more energy is required to heat the body which start out 100 K colder than in the other simulations; 2) energy is needed in order to sublimate the carbon dioxide. As a consequence, the water sublimation is somewhat delayed. 17 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 temperature profile in the upper 0.2 m differed at most a few Kelvin between models with and without airfall. Because this thermal relaxation time is rather short compared to the orbital period we consider Fig. 10 a sufficiently good representation of temperatures at the target sites. A detailed comparison of outbound temperatures indicates that some regions will lose less of the airfall ice they have accumulated than others. In order to compare the sites with each other, we integrated the gas production rates from the outbound equinox to perihelion and normalized with respect to the site with the highest loss, thereby obtaining an index L′ quantifying the relative capacity to get rid of the ice. We then defined S′ = 1 − L′ as the relative capacity to retain ice, and renormalized the set to the highest retainer. Thereby we obtained an index S that we refer to as “survivability”. We also form an index ℛ by normalizing the thickness of the accumulated airfall layer to the highest in the set that we call “reachability”. If a target location simultaneously is characterized by a high reachability (it is reached by much airfall material) and high survivability (it remains comparably cool on the outbound orbital branch), then it is likely to be relatively active when the comet approaches perihelion on the inbound orbital branch. In this system, target site #16 has {ℛ, S } = {1, 1} which means that eastern Ash has the potential of being more active before pre–perihelion than #3 in Ma’at with {ℛ, S } = {0.58, 0.75}, and much more active than #24 in western Ash ({ℛ, S } = {0.11, 0.31}) and #31 in Imhotep ({ℛ, S } = {0.22, 0.07} ). 3.3. Ice loss from airfall deposits To perform the illumination calculations (including shadowing and self heating effects from the irregular nucleus) and the thermophysical simulations for the entire orbit, with rotation resolved (solving the thermophysical equations, Eq. 7–9 requires a time step of ~10 s) is computationally demanding. Therefore, we have not performed these calculations for all 31 target areas, but for a subset of four. These areas were selected to represent each of the major regions under study. Specifically, we consider #3 in Ma’at on the small lobe, #16 in eastern Ash, #24 in western Ash, and #31 in Imhotep on the large lobe. The results of these simulations are shown in Fig. 10. The black curves in those plots show the surface temperature, with the daytime peak as solid curves and the nighttime minimum as dashed curves. The gray areas mark the part of the orbit between the inbound and outbound equinoxes when the level of southern hemisphere activity is highest and most airfall takes place. Target areas #3, #16, and #24 are all located on the northern hemisphere. Because of the spin axis orientation they are illuminated at aphelion and the surface temperatures increases as the comet travels toward perihelion. The temperature peaks around 200 K at all three locations, and the day/night amplitude is typically 10–30 K. However, when the comet approaches the inbound equinox, the importance of the spin axis orientation becomes evident and the temperature plummets, as the regions enter into polar night. During the perihelion passage, the temperatures are as low as 60–70 K (with little to no dependence on nucleus rotational phase) and they are completely inactive. Therefore, airfall material originating from the strongly sublimating southern hemisphere has no problem landing in these target areas and they are flash–frozen after their coma flight. The generally good correspondence between inactivity and the deposition period in gray is the reason why we ignore self–cleaning for these areas. When target areas #3, #16, and #24 emerge from polar night after the outbound equinox at 2.6 AU, the surface temperatures are lower and (because of their strongly non–linear relation) the gas production rates are much lower than on the inbound branch, again due to the spin obliquity effect. Therefore, it is expected that much of the airfall ice collected at perihelion will survive the journey towards aphelion. Consistent with observations (e.g., Gulkis et al., 2015; Hässig et al., 2015; Fink et al., 2016), the northern hemisphere produces water vapor vigorously within ~3 AU pre–perihelion because the airfall material in these regions reach T≳190–200 K for the first time since it was expelled from the southern hemisphere six years earlier. The situation for target area #31, that is located on the southern hemisphere fairly close to the equator, is quite different. It is poorly illuminated at aphelion, and its peak surface temperature is ~40 K below the others. However, as the region approaches the Sun it remains fully illuminated and the temperature peaks at ~230 K, which is the highest among the four by far. Although it is bombarded with airfall material from even more active regions closer to the subsolar point at perihelion, it is capable of re–ejecting this material or preventing it from landing in the first place. That is why we considered self–cleaning for target areas #30–31. Fig. 10 also shows solid and dashed red curves. Those are the maximum and minimum temperatures 0.2 m below the surface. Because of the low conductivity of the porous comet material, the diurnal surface temperature amplitude is rapidly damped with depth, and at this level of magnification the solid and dashed curves are inseparable. In brief, 0.2 m below the surface, a thermometer would not be capable of telling whether the Sun is above or below the local horizon. We note that these calculations were performed without accounting for airfall. We ran one simulation for target area #16 where a 0.5 m thick layer at temperature 200 K was deposited instantaneously onto the nucleus at perihelion when this location had a surface temperature of ~70 K. This procedure was intended to simulate the arrival of warm airfall (originating from the fully illuminated southern hemisphere) to the northern hemisphere that had polar night. After ~150 days the 4. Discussion A number of authors have discussed the transfer of material from one location on 67P to another, the details of the airfall deposition process, and its effect on observable properties, using different methods. We here recapitulate their work and compare their results with our own. Pajola et al. (2017a) calculated the maximum daily lift pressure (the product between production rate and expansion velocity that is proportional to the maximum size of liftable chunks) for different regions on 67P and used that information to offer an explanation for the differing size distributions at Agilkia (coarse and fine chunks) and Sais (coarse chunks with D≳3 cm only). Their discussion provides valuable insight into the dynamics of dust deposition and self–cleaning. The region Bes on the southern hemisphere reaches the highest lift pressure at perihelion and is expected to send a wide range of chunk sizes toward the northern hemisphere. Hapi experiences an extensive polar night period and can collect a substantial amount of material. Once illuminated, its own lift pressure is many orders of magnitude smaller than that of Bes, which, according to Pajola et al. (2017a), may explain why the coma in August 2014 and in the following months (primarily fed by Hapi contributions) was dominated by rather small (0.1–1 mm) chunks. Hapi is unable to eject most of the large chunks it receives and may experience an unusually large accumulation of airfall material over time, unless fragmentation of dm–m chunks into 0.1–1 mm chunks on the surface is efficient. As a contrast, the lift pressure at Sais (once it emerges from polar night) is just an order of magnitude lower than for Bes, suggesting that all but the largest airfall chunks will be ejected after a brief residence on the nucleus. This may explain why the size distribution at that location (measured at the end of September 2016, 3.8 AU post–perihelion) was skewed towards coarse chunks. Agilkia, on the other hand, does not experience polar night at perihelion, which prevents it from accumulating chunks smaller than ~0.1 m at that time. However, it was observed in November 2014 (3.0 AU from the Sun, inbound), at a time when its lift pressure had been lower than that of Hapi since aphelion. The fine airfall seen at Agilkia but missing in Sais therefore originated from Hapi while its coarse airfall came from Bes and other regions on the southern hemisphere (Pajola et al., 2017a). In our work we only consider self–cleaning from regions #30–31 that are strongly illuminated around perihelion. This means that we ignore the self–cleaning observed by Pajola et al. (2017a) in, e.g., Sais at Ma’at, 18 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 because the level of activity after the outbound equinox is still very low (Fig. 10). Also, on the inbound leg we ignore the contribution arriving to our target areas from Hapi. In both cases, our assumptions concern very small particles, and we have already demonstrated that the mass carried by chunks smaller than a few centimeters is a relatively small fraction of the total mass (see Sec. 2.6). Therefore, our results are only slightly affected by this simplification (because we here only consider the total amount of mass, and do not attempt to accurately describe the size distribution of accumulated airfall at the target sites). Kramer and Noack (2015) performed numerical integration of dust chunk trajectories under the influence of realistic nucleus gravity, centrifugal and Coriolis forces (as they were working in a non–inertial frame), and gas drag. Their purpose was to explain the orientation of wind tails observed by Mottola et al. (2015) at Agilkia. They assumed that each point on the nucleus surface emitted the same fixed amount of gas flux. Thus, there was no dependence of local outgassing rate on solar location, the model nucleus had no day/night side, and the model coma was static and highly isotropic compared to the comet coma. However, the irregular nucleus shape created local concentrations of gas in strongly concave areas, like the Hapi neck region. The resulting dust flow had two major characteristics; 1) a flow out of Hapi that diverged into one flow over Seth and Ash on the large lobe and another flow over Hathor, Ma’at, and Hatmehit on the small lobe; and 2) a clockwise flow in the equatorial plane, as seen from the cometary north pole in a body–fixed frame, caused by the counter–clockwise nucleus rotation. Using this model, Kramer and Noack (2015) obtained dust flow vectors in the Agilkia region that were consistent with the wind tail directions observed by Mottola et al. (2015). Lai et al. (2017) performed an investigation of airfall deposition that also considered realistic nucleus gravity, centrifugal and Coriolis forces, and gas drag. In their work, a fixed outgassing rate was applied to each point on the surface. These outgassing rates were taken to be proportional to the average cosine of the local day–time solar zenith angle during nucleus rotation. This anisotropic gas flow roughly captures the dependence of coma number density with nucleus latitude but smears its local time dependence, except for a strong day/night asymmetry. This anisotropic coma was merged at a 0.85 weight with an isotropic component with weight 0.15, and the total gas production rate was normalized to match the observed one. This was done once per month from January–December 2015 in order to account for the changing subsolar latitude and evolving total gas production rate with heliocentric distance. For each of the dozen static outgassing patterns, Direct Simulation Monte Carlo simulations were performed to calculate the gas number density, translational temperature, and expansion velocity as functions of position within the coma, needed for gas drag evaluations. Lai et al. (2017) used 14 dust classes with sizes in the 2.8 μm–4.5 cm range, each represented by ~4 ⋅ 105 particles launched from random locations on the surface that had different weights depending on local gas production conditions and chunk cross sections. In this manner, the local airfall and self–cleaning rates could be calculated at a given time, and the results were integrated over the considered orbital arc to estimate the local net deposition during one orbit. We note that Thomas et al. (2015a) applied a similar model, however, limited the study to ejection of dust from Hapi to other parts of the northern hemisphere. They found that ejection speeds of vd ≳0.5 m s−1 are necessary to leave the Hapi valley and that ≳50% of the particles escape the nucleus if vd ≳1 m s−1 , thus providing tight speed constraints for that transfer route. Lai et al. (2017) found that the largest liftable chunk on the southern hemisphere has D = 1 cm and that a net mass loss occurs everywhere except in Hapi where a net airfall deposit of 0.4 m depth is formed. We note that the D ≤ 1 cm limit is small compared to the observed numerous coma chunks in the D = 0.1–1 m class (e.g., Rotundi et al., 2015; Davidsson et al., 2015; Agarwal et al., 2016), as well as compared to the size of chunks observed in resolved smooth terrains. Furthermore, their 0.5–1 m net loss on most of the northern hemisphere, combined with their 1–1.8 m net loss on most of the southern hemisphere, translates to a total loss of ~2.5 ⋅ 1010 kg (if applying an average loss of 1 m, a nucleus surface area of 46.9 km2 and assuming that the density is that of the bulk nucleus, 530 kg m−3, see Jorda et al., 2016). This is 2.4 times higher than the observed nucleus net mass loss of 1.05(±0.34) ⋅ 1010 kg according to the Rosetta/RSI radio science experiment (Pätzold et al., 2019). We suspect that the maximum liftable dust size at the perihelion subsolar point becomes significantly underestimated when calculating production rates based on the average cosine of the solar zenith angle. As a consequence, the amount of airfall is underestimated because the size distribution of ejected particles is truncated at a size below which return as airfall is highly uncommon. Simultaneously, the outgassing is rather strong in polar night regions (15% of the daytime production), which leads to a high self–cleaning rate. Except for an initial acceleration phase aimed at determining a relation between local gas production rates, dust size, and their velocities, the current work does not consider gas drag effects on the dust orbits. This is admittedly a weakness. However, we do not think it is meaningful to engage in computationally demanding gas coma modeling that enables continuous gas drag, unless the nucleus model source that feeds the coma model is reproducing the strong temporal, latitudinal, and longitudinal variations in gas production rate known to occur on the real nucleus (during nucleus rotation and throughout the orbit). The summary above reminds us that such an elaborate treatment is beyond the current state–of–the–art of the field. It is unlikely that attempts to account for the orbital perturbations due to outgassing from the chunks would be realistic (due to their small size, chunks are likely having an excited spin state, and when combined with thermal inertia effects the time–evolution of the net non–gravitational force vector is unpredictable). Therefore, this “rocket effect” is ignored as well. Because of the point–mass assumption we also have not applied a realistic gravity field, but based on the comparison with a more appropriate approach described in Sec. 2.2 we think that this may not have drastic consequences concerning our overall conclusions. Furthermore, we have ignored the effect of solar radiation pressure on the dynamics of coma chunks. Radiation pressure does have an affect on micron–sized chunks (e.g., Tenishev et al., 2011), but considering the small contribution of such particles to airfall deposits in terms of mass, this simplification has no important effect on our conclusions. With these limitations in mind, our calculations support the hypothesis (Keller et al., 2015a, 2017) of a significant transport route from the southern to the northern hemisphere and demonstrate that both Ash and Ma’at are plausible destinations of airfall, as previously has been suggested (Thomas et al., 2015b; Thomas et al., 2015a). We find that the amount of airfall varies systematically between the four main regions of study (Ma’at, eastern and western Ash, and Imhotep) and within those regions. About 40% of the considered target areas receive at most 0.5 m airfall and about 30% receive at least 1 m according to our study. Considering the global coverage OSIRIS resolution of 0.2–3 m px−1 (Preusker et al., 2017) it is therefore not surprising that morphological changes because of airfall are difficult or impossible to verify observationally in many locations, but that such changes definitively are detected in others (Hu et al., 2017). Our second goal consisted of estimating the amount of water ice loss during the transfer of chunks through the coma. We also wanted to explore the survivability of the substantially more volatile CO2 because of the chemical north/south dichotomy measured at 67P (Hässig et al., 2015; Fougere et al., 2016a; Fink et al., 2016) and interpreted in the context of coma transport and airfall (Keller et al., 2017). For a nominal set of model parameters our NIMBUS simulations showed that a dm–sized chunk would retain more than 90% of its water ice when suspended for 12 h and that a cm–sized chunk would keep almost half its water ice in the same time. This lends substantial credibility to the claim by Keller et al. (2017) that the airfall is “wet” and naturally explains the substantial level of water–driven activity on the northern hemisphere observed by Rosetta. We also found that a dm–sized chunk 19 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Fig. 11. Angle between the pole facing sunward at perihelion and the comet–Sun vector versus time. Spacecraft flyby dates are indicated. Perihelion is marked with a vertical line. The airfall deposit process is least efficient on the horizontal line (no polar night). containing 5% CO2 relative to water would lose all of it in 2 h. It also means that a meter–sized chunk might be able to transfer supervolatiles from the south to the north. If so, then the low levels of CO2 outgassing from the northern hemisphere is a consequence of the small number of such bodies (and it may indicate that most CO2 ice could have been lost already prior to departure). Alternatively, some CO2 may have been transported from the south to the north, not as an independently condensed form of ice but trapped within water ice. Laboratory experiments (e.g., Edridge et al., 2013) demonstrate that a fraction of CO2 trapped in amorphous ice survives the release during crystallization as well as the cubic–hexagonal transformation and is released during water sublimation. We did not model that process explicitly. We note that Philae detected some highly volatile species at Agilika, such as CO and CO2 with Ptolemy (Wright et al., 2015), and CH4 with COSAC (Goesmann et al., 2015). Those species may have been trapped in water ice within airfall material and released during water sublimation. Alternatively, they may have originated from underneath the airfall layer. NIMBUS simulations by Davidsson et al. (2020) show that the measured pre–perihelion CO2 production rate curve is reproduced when the depth of the CO2 sublimation front is located 3.8 m below the surface on the northern hemisphere, and at 1.9 m below the surface on the southern hemisphere at aphelion (these depths are modified in the model over time because of dust mantle erosion as well as CO2 sublimation, see Sec. 3.2). Supervolatiles like CO and CH4 may be released during sublimation of CO2 at the front, or from even larger depths by segregation processes, if they are trapped in CO2, as suggested by Gasc et al. (2017). Although presence of CO2 in the airfall material itself remains a possibility, it cannot be the major source of CO2 release from the northern hemisphere before the inbound equinox, because the CO2 production does not correlate with that of H2O (Luspay-Kuti et al., 2015; Gasc et al., 2017). Our third goal was to quantify the relative capacity of different airfall deposition sites to retain their water ice content after deposition. To this end we introduced a survivability index S and combined it with a reachability index ℛ that measures the relative capability to collect airfall material. We performed an investigation limited to four sites, where #16 in eastern Ash came out on top, with #3 in Ma’at second, and #24 (western Ash) and #31 (Imhotep) on a shared third place. Comparing with the H2O potential activity map of Fougere et al. (2016a) derived from long measurement series by ROSINA and corrected for illumination biases, it is clear that #16 is substantially more active than #24 and #31, which is consistent with our work. However, #3 is the most active of them all. Interestingly, the high–activity belt shown by Fougere et al. (2016a), stretching from Ash, Seth, to Ma’at along the (anti)meridian with an epicenter in Hapi, has a strong resemblance with the simulated model flow pattern of airfall material originating from Hapi (Thomas et al., 2015a; Kramer and Noack, 2015). It is therefore possible that the ice variability introduced by airfall from the southern hemisphere at perihelion is mixed and masked by redistribution of airfall due to Hapi activity at the time ROSINA performed the measurements. It is encouraging that the activity map of Lara et al. (2015), based on the footprint of dust jets observed early during the Rosetta mission, suggest high levels of activity both at #3 and #16. To what extent are airfall deposits common features on comet nuclei? We expect that three criteria need to be fulfilled in order to create thick and wide–spread airfall deposits. First, the spin pole needs to have such an orientation that a large fraction of the nucleus surface experiences polar night at perihelion. In terms of the spin axis obliquity and argument (for definitions see, e.g., Sekanina, 1981; Davidsson and Gutiérrez, 2004) the most favorable configuration is an obliquity near 90∘ and an argument near 90∘ or 270∘ (i.e., the perihelion subsolar point coincides with one of the comet nucleus poles). Second, the post–perihelion equinox should take place as late as possible. The equinox passage means that the Sun is in the equatorial plane of the comet and that fresh airfall deposits no longer are protected from solar illumination. If airfall is illuminated sufficiently close to perihelion, strong activity and self–cleaning may take place that prevents long–term deposits from forming. Third, the comet needs to produce abundant chunks that are too large to reach escape velocity. The presence of a detectable debris trail is a potential indicator that a specific comet is prone to eject particularly large chunks. A Spitzer survey of 34 comets showed that 27 objects, or 79%, had trails consisting of millimeter–sized particles and larger (Reach et al., 2007), suggesting that 67P is not unusual. Fig. 11 shows the time evolution of the angle between the solar direction and the comet nucleus pole that is sunlit at perihelion (here referred to as the polar angle), for four comets that were targets of spacecraft missions. The smaller the polar angle, the larger the fraction 20 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 5. Conclusions of the comet surface has polar night at perihelion. If it reaches 90∘ (marked by the horizontal line), equinox takes place and the Sun is in the equatorial plane of the comet (i.e., no region has polar night). Among the four, Comets 67P and 19P/Borrelly seem to have the best conditions for airfall deposition. The perihelion polar angles are 42∘ for 67P and 29∘ for 19P/Borrelly, i.e., the polar night regions are large. Comet 67P reaches equinox later (221 days post–perihelion, outside the figure) than 19P/Borrelly (97 days post–perihelion), but self–cleaning is probably limited in both cases. 67P evidently produces large chunks (as observed by Rosetta, but also inferred from the existence of a trail; Sykes and Walker, 1992; Kelley et al., 2009; Ishiguro et al., 2009), that enables airfall deposit formation, but this may not be the case for 19P/Borrelly. The reports on the Deep Space 1 flyby of 19P/Borrelly do not mention large coma chunks (Soderblom et al., 2002; Boice et al., 2002), and the comet does not have a trail (Ishiguro et al., 2009). Despite the favorable geometric conditions, Comet 19P/Borrelly may therefore not have large smooth terrains. Unfortunately, Deep Space 1 images cannot provide an answer because the flyby took place just 8 days after perihelion (Soderblom et al., 2002), i.e., a potential airfall coverage would have been hidden from view on the dark polar night side. Comets 81P/Wild 2, and particularly 9P/Tempel 1, are least suitable for vast airfall deposit formation. 9P/Tempel 1 has a perihelion polar angle of 86∘ and passes equinox just three weeks later. It is therefore surprising that 9P/Tempel 1 has prominent smooth areas (A’Hearn et al., 2005). Particularly the smooth area “S2” near the south pole (Veverka et al., 2013) would be most consistent with airfall onto the weakly illuminated and least active part of the nucleus, while the origin of localized smooth areas elsewhere is less obvious. It is also interesting to note that observations of a slowly expanding coma arc during the Deep Impact flyby suggests the presence of fairly large (millimeter–sized) slow–moving chunks (Farnham et al., 2007), and a trail associated with Comet 9P/Tempel 1 has been observed both with IRAS (Sykes and Walker, 1992) and with Spitzer (Reach et al., 2007). This suggests that the spin axis orientation plays a minor role for the existence of smooth airfall terrain, though a large coverage probably still requires a small perihelion polar angle. If true, this could mean that all comets with an appropriate dust size distribution are capable of building smooth airfall terrains. For that reason, the case of 81P/Wild 2 is somewhat puzzling. With a polar angle of 70∘, a rather large part of the nucleus has polar night at perihelion. The equinox passage is just 40 days post–perihelion, but judging from Comet 9P/Tempel 1, that may not prevent airfall debris from remaining on the surface. The Stardust spacecraft collided with a couple of mm–sized dust grains (Tuzzolino et al., 2004) and the comet has a trail according to Ishiguro et al. (2009). Furthermore, the Stardust flyby took place 98 days post–perihelion (Brownlee et al., 2004), at a time when the regions that experienced polar night at perihelion had become visible. Yet, the Stardust images (Brownlee et al., 2004) do not reveal smooth terrains of the type seen on 9P/Tempel 1 or like Imhotep or Hapi on 67P. Potentially, the airfall material on 81P/Wild 2 is hiding in plain view, i.e., as a thin coverage on top of the underlying rugged topography, similar to the Ash and Seth regions on 67P. The nature of those regions only became evident in Rosetta images that had significantly better resolution than the best Stardust images (14 m px−1; Brownlee et al., 2004). Alternatively, airfall deposits do not form on Comet 81P/Wild 2. Two spacecraft targets, Comets 1P/Halley and 103P/Hartley 2, were omitted from Fig. 11 because they are complex rotators (the polar angle is not easily calculated). Large chunks are common around 103P/ Hartley 2 (Kelley et al., 2013) and it has a trail (Reach et al., 2007). EPOXI images revealed a smooth neck region (A’Hearn et al., 2011) reminiscent of Hapi on 67P. With airfall confirmed on 67P, and suspected on 9P/Tempel 1 and 103P/Hartley 2, it indeed appear to be a common process in comets. We have modeled the transfer of material from the highly active southern hemisphere of comet 67P to its northern hemisphere during the perihelion passage. We find that such transport routes exist and that the Ash and Ma’at regions receive plenty of airfall, thereby supporting previous hypotheses (Thomas et al., 2015b; Thomas et al., 2015a; Keller et al., 2015a, 2017). The amount of airfall accumulated between the inbound and outbound equinoxes is typically a few times 0.1–1 m (some of which will be removed in other parts of the orbit, i.e., these numbers should not be interpreted as net orbital accumulations, but rather as an estimate of the material the northern hemisphere has to its disposal to drive activity pre–perihelion). The distribution of airfall is heterogeneous, with the most substantial accumulation in eastern Ash, followed by western Ash and Ma’at at similar levels, with the least airfall at central Imhotep because of efficient self–cleaning (only these four regions were considered, e.g., Hapi was not included in the study). We have also modeled the loss of both H2O and CO2 from 0.01–0.1 m diameter chunks in the coma, using the elaborate NIMBUS model. We find that a cm–sized chunk can retain roughly half its water ice during a 12h exposure to the Sun at perihelion, while a dm–sized chunk holds on to more than 90% of the water ice. If there is CO2 at a 5% level relative to water by number, a dm–sized chunk loses all carbon dioxide in two hours. We therefore support the scenario described by Keller et al. (2017) of a “wet” airfall deprived of supervolatiles that is responsible for the observed water–dominated comet activity prior to the inbound equinox. Finally, we studied the longterm loss of ice following the near–perihelion airfall accumulation. We demonstrated that the surface temperatures are kept comparably cool on the northern hemisphere during the outbound orbital branch, which helps explain the high level of activity in the north on the inbound orbital branch when the deposits finally are heated to high temperature. We introduce the reachability and survivability indices ℛ and S as a way to classify the potential of a region to have a high level of activity by simultaneously accumulating a high amount of wet airfall and to preserve it to the following perihelion passage. This work also serves as a guide for selecting the most promising volatile rich smooth terrain sites in the northern hemisphere for future comet nucleus sample return missions. Declaration of Competing Interest None. Acknowledgments Parts of the research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We gratefully acknowledge JPL funding. References A’Hearn, M.F., Belton, M.J.S., Delamere, W.A., Kissel, J., Klaasen, K.P., McFadden, L.A., Meech, K.J., Melosh, H.J., Schultz, P.H., Sunshine, J.M., Thomas, P.C., Veverka, J., Yeomans, D.K., Baca, M.W., Busko, I., Crockett, C.J., Collins, S.M., Desnoyer, M., Eberhardy, C.A., Ernst, C.M., Farnham, T.L., Feaga, L., Groussin, O., Hampton, D., Ipatov, S.I., Li, J.-Y., Lindler, D., Lisse, C.M., Mastrodemos, N., Owen Jr., W.M., Richardson, J.E., Wellnitz, D.D., White, R.L., 2005. Deep Impact: Excavating Comet Tempel 1. Science 258, 258–264. A’Hearn, M.F., Belton, M.J.S., Delamere, W.A., Feaga, L.M., Hampton, D., Kissel, J., Klaasen, K.P., McFadden, L.A., Meech, K.J., Melosh, H.J., Schultz, P.H., Sunshine, J. M., Thomas, P.C., Veverka, J., Wellnitz, D.D., Yeomans, D.K., Besse, S., Bodewits, D., Bowling, T.J., Carcich, B.T., Collins, S.M., Farnham, T.L., Groussin, O., Hermalyn, B., Kelley, M.S., Kelley, M.S., Li, J.-Y., Lindler, D.J., Lisse, C.M., McLaughlin, S.A., Merlin, F., Protopapa, S., Richardson, J.E., Williams, J.L., 2011. EPOXI at comet Harley 2. Science 332, 1396–1400. Agarwal, J., A’Hearn, M.F., Vincent, J.-B., Güttler, C., Höfner, S., Sierks, H., Tubiana, C., Barbieri, C., Lamy, P.L., Rodrigo, R., Koschny, D., Rickman, H., Barucci, M.A., Bertaux, J.-L., Bertini, I., Boudreault, S., Cremonese, G., Da Deppo, V., Davidsson, B., 21 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Davidsson, B.J.R., Samarasinha, N., Farnocchia, D., 2020. Modeling the water and carbon dioxide production rates of Comet 67P/Churyumov–Gerasimenko. In: Preparation. Edridge, J.L., Freimann, K., Burke, D.J., Brown, W.A., 2013. Surface science investigations of the role of CO2 in astrophysical ices. Phil. Trans. R. Soc. A 371, 20110578. El-Maarry, M.R., Thomas, N., Garcia Berná, A., Marschall, R., Auger, A.-T., Groussin, O., Massironi, M., Marchi, S., Preusker, F., Scholten, F., Jorda, L., Kührt, E., Hofmann, M., Hoefner, S., Deller, J., the OSIRIS team, 2015a. Fractures on comet 67P/Churyumov–Gerasimenko observed by the Rosetta/OSIRIS camera. Geophys. Res. Lett. 42 (13), 5170–5178. El-Maarry, M.R., Thomas, N., Giacomini, L., Massironi, M., Pajola, M., Sierks, H., Barbieri, C., Lamy, P., Rodrigo, R., Rickman, H., Koschny, D., Keller, H., Agarwal, J., A’Hearn, M.F., Auger, A.-T., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Deppo, V.D., Davidsson, B., Cecco, M.D., Debei, S., Güttler, C., Fornasier, S., Fulle, M., Groussin, O., Gutiérrez, P.J., Hviid, S., Ip, W.-H., Jorda, L., Knollenberg, J., Kramm, J.-R., Kührt, E., Küppers, M., Lara, L.M., Lazzarin, M., Moreno, J.J.L., Marchi, S., Marzari, F., Michalik, H., Naletto, G., Oklay, N., Pommerol, A., Preusker, F., Scholten, F., Tubiana, C., Vincent, J.-B., Wenzel, K.-P., 2015b. Regional surface morphology of comet 67P/Churyumov–Gerasimenko. Astron. Astrophys. 583, A26. El-Maarry, M.R., Thomas, N., Gracia-Berná, A., Pajola, M., Lee, J.-C., Massironi, M., Davidsson, B., Marchi, S., Keller, H.U., Hviid, S.F., Besse, S., Sierks, H., Barbieri, C., Lamy, P.L., Koschny, D., Rickman, H., Rodrigo, R., A’Hearn, M.F., Auger, A.-T., Barucci, M.A., Bertaux, J.-L., Bertini, I., Bodewits, D., Cremonese, G., Da Deppo, V., De Cecco, M., Debei, S., Güttler, C., Fornasier, S., Fulle, M., Giacomini, L., Groussin, O., Gutierrez, P.J., Ip, W.-H., Jorda, L., Knollenberg, J., Kovacs, G., Kramm, J.-R., Kührt, E., Küppers, M., Lara, L.M., Lazzarin, M., Lopez Moreno, J.J., Marschall, R., Marzari, F., Naletto, G., Oklay, N., Pommerol, A., Preusker, F., Scholten, F., Tubiana, C., Vincent, J.-B., 2016. Regional surface morphology of comet 67P/Churyumov–Gerasimenko from Rosetta/OSIRIS images: The southern hemisphere. Astron. Astrophys. 593, A110. El-Maarry, M.R., Groussin, O., Thomas, N., Pajola, M., Auger, A.-T., Davidsson, B., Hu, X., Hviid, S.F., Knollenberg, J., Güttler, C., Tubiana, C., Fornasier, S., Feller, C., Hasselmann, P., Vincent, J.-B., Sierks, H., Barbieri, C., Lamy, P., Rodrigo, R., Koschny, D., Keller, H.U., Rickman, H., A’Hearn, M.F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Besse, S., Bodewits, D., Cremonese, G., Da Deppo, V., Debei, S., De Cecco, M., Deller, J., Deshapriya, J.D.P., Fulle, M., Gutierrez, P.J., Hofmann, M., Ip, W.-H., Jorda, L., Kovacs, G., Kramm, J.-R., Kührt, E., Küppers, M., Lara, L.M., Lazzarin, M., Lin, Z.-Y., Lopez Moreno, J.J., Marchi, S., Marzari, F., Mottola, S., Naletto, G., Oklay, N., Pommerol, A., Preusker, F., Scholten, F., Shi, X., 2017. Surface changes on comet 67P/Churyumov–Gerasimenko suggest a more active past. Science 355 (6332), 1392–1395. Fanale, F.P., Salvail, J.R., 1984. An idealized short–period comet model: Surface insolation, H2O flux, dust flux and mantle evolution. Icarus 60, 476–511. Farnham, T.L., Wellnitz, D.D., Hampton, D.L., Li, J.-Y., Sunshine, J.M., Groussin, O., McFadden, L.A., Crockett, C.J., A’Hearn, M.F., Belton, M.J.S., Schultz, P., Lisse, C.M., 2007. Dust coma morphology in the Deep Impact images of Comet 9P/Tempel 1. Icarus 187, 26–40. Filacchione, G., Raponi, A., Capaccioni, F., Ciarniello, M., Tosi, F., Capria, M.T., Sanctis, M.C.D., Migliorini, A., Piccioni, G., Cerroni, P., Barucci, M.A., Fornasier, S., Schmitt, B., Quirico, E., Erard, S., Bockelee-Morvan, D., Leyrat, C., Arnold, G., Mennella, V., Ammannito, E., Bellucci, G., Benkhoff, J., Bibring, J.P., Blanco, A., Blecka, M.I., Carlson, R., Carsenty, U., Colangeli, L., Combes, M., Combi, M., Crovisier, J., Drossart, P., Encrenaz, T., Federico, C., Fink, U., Fonti, S., Fulchignoni, M., Ip, W.-H., Irwin, P., Jaumann, R., Kuehrt, E., Langevin, Y., Magni, G., McCord, T., Moroz, L., Mottola, S., Palomba, E., Schade, U., Stephan, K., Taylor, F., Tiphene, D., Tozzi, G.P., Beck, P., Biver, N., Bonal, L., Combe, J.-P., Despan, D., Flamini, E., Formisano, M., Frigeri, A., Grassi, D., Gudipati, M.S., Kappel, D., Longobardo, A., Mancarella, F., Markus, K., Merlin, F., Orosei, R., Rinaldi, G., Cartacci, M., Cicchetti, A., Hello, Y., Henry, F., Jacquinod, S., Reess, J. M., Noschese, R., Politi, R., Peter, G., 2016. Seasonal exposure of carbon dioxide ice on the nucleus of comet 67P/Churyumov–Gerasimenko. Science 354 (6319), 1563–1566. Fink, U., Doose, L., Rinaldi, G., Bieler, A., Capaccioni, F., Bockelée-Morvan, D., Filacchione, G., Erard, S., Leyrat, C., Blecka, M., Capria, M.T., Combi, M., Crovisier, J., De Sanctis, M.C., Fougere, N., Taylor, F., Migliorini, A., Piccioni, G., 2016. Investigation into the disparate origin of CO2 and H2O outgassing for Comet 67/P. Icarus 277, 78–97. Fornasier, S., Hasselmann, P.H., Barucci, M., Feller, C., Besse, S., Leyrat, C., Lara, L., Oklay, N., Tubiana, C., Scholten, F., Sierks, H., Barbieri, C., Lamy, P.L., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., A’Hearn, M.F., Bertaux, J.-L., Bertini, I., Cremonese, G., Deppo, V.D., Davidsson, B., Debei, S., Cecco, M.D., Fulle, M., Groussin, O., Güttler, C., Gutiérrez, P.J., Gutierrez-Marques, P., Hviid, S.F., Ip, W.-H., Jorda, L., Knollenberg, J., Kramm, R., Kührt, E., Küppers, M., La Forgia, F., Lazzarin, M., Lopez Moreno, J.J., Marzari, F., Massironi, M., Matz, K.-D., Michalik, H., Mottola, S., Naletto, G., Pajola, M., Pommerol, A., Preusker, F., Snodgrass, C., Thomas, N., Vincent, J.-B., 2015. Spectro–photometric properties of the 67P/Churyumov–Gerasimenko’s nucleus from the OSIRIS instrument onboard the Rosetta spacecraft. Astron. Astrophys. 583, A30. Fougere, N., Altwegg, K., Berthelier, J.-J., Bieler, A., Bockelée-Morvan, D., Calmonte, U., Capaccioni, F., Combi, M.R., De Keyser, J., Debout, V., Erard, S., Fiethe, B., Filacchione, G., Fink, U., Fuselier, S.A., Gombosi, T.I., Huang, Z., Le Roy, L., Leyrat, C., Migliorini, A., Piccioni, G., Rinaldi, G., Rubin, M., Shou, Y., Tenishev, V., Toth, G., Tzou, C.-Y., the VIRTIS team, and the ROSINA team, 2016a. Three–dimensional direct simulation Monte–Carlo modeling of the coma of comet Debei, S., De Cecco, M., Deller, J., Fornasier, S., Fulle, M., Gicquel, A., Groussin, O., Gutiérrez, P.J., Hofmann, M., Hviid, S.F., Ip, W.-H., Jorda, L., Keller, H.U., Knollenberg, J., Kramm, J.-R., Kührt, E., Küppers, M., Lara, L.M., Lazzarin, M., Lopez Moreno, J.J., Marzari, F., Naletto, G., Oklay, N., Shi, X., Thomas, N., 2016. Acceleration of individual, decimetre–sized aggregates in the lower coma of comet 67P/Churyumov-Gerasimenko. Mon. Not. R. Astron. Soc. 462 (Suppl_1), S78–S88. Anisimov, S.I., 1968. Vaporization of metal absorbing laser radiation. Soviet Physics JETP 27 (1), 182–183. Auger, A.-T., Groussin, O., Jorda, L., Bouley, S., Gaskell, R., Lamy, P.L., Capanna, C., Thomas, N., Pommerol, A., Sierks, H., Barbieri, C., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., A’Hearn, M.F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Da Deppo, V., Davidsson, B., Debei, S., De Cecco, M., ElMaarry, M.R., Fornasier, S., Fulle, M., Gutiérrez, P.J., Güttler, C., Hviid, S., Ip, W.-H., Knollenberg, J., Kramm, J.-R., Kührt, E., Küppers, M., La Forgia, F., Lara, L.M., Lazzarin, M., Lopez Moreno, J.J., Marchi, S., Marzari, F., Massironi, M., Michalik, H., Naletto, G., Oklay, N., Pajola, M., Sabau, L., Tubiana, C., Vincent, J.-B., Wenzel, K.P., 2015. Geomorphology of the Imhotep region on comet 67P/ Churyumov–Gerasimenko from OSIRIS observations. Astron. Astrophys. 583, A35. Auger, A.-T., Groussin, O., Jorda, L., El-Maarry, M.R., Bouley, S., Séjourné, A., Gaskell, R., Capanna, C., Davidsson, B., Marchi, S., Höfner, S., Lamy, P.L., Sierks, H., Barbieri, C., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., A’Hearn, M.F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Da Deppo, V., Debei, S., De Cecco, M., Fornasier, S., Fulle, M., Gutiérrez, P.J., Güttler, C., Hviid, S., Ip, W.-H., Knollenberg, J., Kramm, J.-R., Kührt, E., Küppers, M., Lara, L.M., Lazzarin, M., Lopez Moreno, J.J., Marzari, F., Massironi, M., Michalik, H., Naletto, G., Oklay, N., Pommerol, A., Sabau, L., Thomas, N., Tubiana, C., Vincent, J.B., Wenzel, K.-P., 2018. Meter–scale thermal contraction crack polygons on the nucleus of comet 67P/Churyumov–Gerasimenko. Icarus 301, 173–188. Berry, M.M., Healy, L.M., 2004. Implementation of Gauss–Jackson integration for orbit propagation. J. Astronaut. Sci. 52 (3), 331–357. Biele, J., Ulamec, S., Maibaum, M., Roll, R., Witte, L., Jurado, E., Noz, P.M., Arnold, W., Auster, H.-U., Casas, C., Faber, C., Fantinati, C., Finke, F., Fischer, H.-H., Geurts, K., Güttler, C., Heinisch, P., Herique, A., Hviid, S., Kargl, G., Knapmeyer, M., Knollenberg, J., Kofman, W., Kömle, N., Kührt, E., Lommatsch, V., Mottola, S., de Santayana, R. Pardo, Remetean, E., Scholten, F., Seidensticker, K.J., Sierks, H., Spohn, T., 2015. The landing(s) of Philae and inferences about comet surface mechanical properties. Science 349 (6247) (aaa9816–1). Birch, S.P.D., Tang, Y., Hayes, A.G., Kirk, R.L., Bodewits, D., Campins, H., Fernandez, Y., Kutsop, N.W., Sierks, H., Soderblom, J.M., Squyres, S.W., Vincent, J.-B., 2017. Geomorphology of comet 67P/Churyumov–Gerasimenko. Mon. Not. R. Astron. Soc. 469, S50–S67. Bird, G.A., 1994. Molecular gas dynamics and the direct simulation of gas flows. Oxford University Press Inc., New York. Boice, D.C., Soderblom, L.A., Britt, D.T., Brown, R.H., Sandel, B.R., Yelle, R.V., Buratti, B. J., Hicks, M.D., Nelson, R.M., Rayman, M.D., Oberst, J., Thomas, N., 2002. The Deep Space 1 encounter with Comet 19P/Borrelly. Earth, Moon, Planets 89, 301–324. Boulet, D., 1991. Methods of Orbit Determination for the Micro Computer. Willmann–Bell, Inc., Richmond, VA. Brownlee, D.E., Horz, F., Newburn, R.L., Zolensky, M., Duxbury, T.C., Sandford, S., Sekanina, Z., Tsou, P., Hanner, M.S., Clark, B.C., Green, S.F., Kissel, J., 2004. Surface of young Jupiter family comet 81P/Wild 2: View from the Stardust spacecraft. Science 304, 1764–1769. Burden, R.L., Faires, J.D., 1993. Numerical analysis. PWS Publishing Company. Cercignani, C., 2000. Rarefied gas dynamics. From basic concepts to actual calculations. Cambridge University Press, Cambridge. Crifo, J.F., Loukianov, G.A., Rodionov, A.V., Khanlarov, G.O., Zakharov, V.V., 2002. Comparison between Navier–Stokes and Direct Monte–Carlo simulations of the circumnuclear coma. I. Homogeneous, spherical source. Icarus 156, 249–268. Danby, J.M.A., 1989. Fundamentals of celestial mechanics. Willmann–Bell, Inc., Richmond. Davidsson, B.J.R., 2008. Comet Knudsen layers. Space Sci. Rev. 138, 207–223. Davidsson, B.J.R., 2020. Thermophysical evolution of planetesimals in the Primordial Disk. In: preparation. Davidsson, B.J.R., Gutiérrez, P.J., 2004. Estimating the nucleus density of Comet 19P/ Borrelly. Icarus 168 (2), 392–408. Davidsson, B.J.R., Rickman, H., 2014. Surface roughness and three–dimensional heat conduction in thermophysical models. Icarus 243, 58–77. Davidsson, B.J.R., Skorov, Y.V., 2002. On the light–absorbing surface layer of cometary nuclei.II. Thermal modeling. Icarus 159, 239–258. Davidsson, B.J.R., Skorov, Y.V., 2004. A practical tool for simulating the presence of gas comae in thermophysical modeling of cometary nuclei. Icarus 168 (1), 163–185. Davidsson, B.J.R., Gulkis, S., Alexander, C., von Allmen, P., Kamp, L., Lee, S., Warell, J., 2010. Gas kinetics and dust dynamics in low–density comet comae. Icarus 210, 455–471. Davidsson, B.J.R., Gutiérrez, P.J., Sierks, H., Barbieri, C., Lamy, P.L., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., A’Hearn, M.F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Bodewits, D., Cremonese, G., Da Deppo, V., Debei, S., De Cecco, M., Fornasier, S., Fulle, M., Groussin, O., Güttler, C., Hviid, S.F., Ip, W.-H., Jorda, L., Knollenberg, J., Kovacs, G., Kramm, J.-R., Kührt, E., Küppers, M., La Forgia, F., Lara, L.M., Lazzarin, M., Lopez Moreno, J.J., Lowry, S., Magrin, S., Marzari, F., Michalik, H., Moissl-Fraund, R., Naletto, G., Oklay, N., Pajola, M., Snodgrass, C., Thomas, N., Tubiana, C., Vincent, J.-B., 2015. Orbital elements of the material surrounding comet 67P/Churyumov–Gerasimenko. Astron. Astrophys. 583, A16. 22 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 67P/Churyumov–Gerasimenko observed by the VIRTIS and ROSINA instruments on board Rosetta. Astron. Astrophys. 588, A134. Fougere, N., Altwegg, K., Berthelier, J.-J., Bieler, A., Bockelée-Morvan, D., Calmonte, U., Capaccioni, F., Combi, M.R., De Keyser, J., Debout, V., Erard, S., Fiethe, B., Filacchione, G., Fink, U., Fuselier, S.A., Gombosi, T.I., Huang, L. Le Roy, Leyrat, C., Migliorini, A., Piccioni, G., Rinaldi, G., Rubin, M., Shou, Y., Tenishev, V., Toth, G., Tzou, C.-Y., the VIRTIS team, the ROSINA team, 2016b. Direct Simulation Monte Carlo modelling of the major species in the coma of comet 67P/ Churyumov–Gerasimenko. Mon. Not. R. Astron. Soc. 462, S156–S169. Fox, K., 1984. Numerical integration of the equations of motion of celestrial mechanics. Cel. Mech. 33, 127–142. Gasc, S., Altwegg, K., Balsiger, H., Berthelier, J.-J., Bieler, A., Calmonte, U., Fiethe, B., Fuselier, S., Galli, A., Gombosi, T., Hoang, M., De Keyser, J., Korth, A., Le Roy, L., Mall, U., Réme, H., Rubin, M., Sémon, T., Tzou, C.-Y., Hunter Waite, J., Wurz, P., 2017. Change of outgassing pattern of 67P/Churyumov-Gerasimenko during the March 2016 equinox as seen by ROSINA. Mon. Not. R. Astron. Soc. 469, S108–S117. Giauque, W.F., Egan, C.J., 1937. Carbon dioxide. The heat capacity and vapor pressure of the solid. The heat of sublimation. Thermodynamic and spectroscopic values of the entropy. J. Chem. Phys. 5, 45–54. Glassmeier, K.-H., Boehnhardt, H., Koschny, D., Kührt, E., Richter, I., 2007. The Rosetta mission: Flying towards the origin of the solar system. Space Sci. Rev. 128, 1–21. Goesmann, F., Rosenbauer, H., Bredehöft, J.H., Cabane, M., Ehrenfreund, P., Gautier, T., Giri, C., Krüger, H., Le Roy, L., MacDermott, A.J., McKenna-Lawlor, S., Meierhenrich, U.J., Caro, G.M. Muñoz, Raulin, F., Roll, R., Steele, A., Steininger, H., Sternberg, R., Szopa, C., Thiemann, W., Ulamec, S., 2015. Organic compounds on comet 67P/Churyumov–Gerasimenko revealed by COSAC mass spectrometry. Science 349 (6247), aab0689. Gombosi, T.I., 1994. Gaskinetic Theory. Cambridge University Press, Cambridge. Groussin, O., Sierks, H., Barbieri, C., Lamy, P., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., A’Hearn, M.F., Auger, A.-T., Barucci, M.A., Bertaux, J.-L., Bertini, I., Besse, S., Cremonese, G., Da Deppo, V., Davidsson, B., Debei, S., De Cecco, M., ElMaarry, M.R., Fornasier, S., Fulle, M., Gutiérrez, P.J., Güttler, C., Hviid, S., Ip, W.-H., Jorda, L., Knollenberg, J., Kovacs, G., Kramm, J.R., Kührt, E., Küppers, M., Lara, L. M., Lazzarin, M., Lopez Moreno, J.J., Lowry, S., Marchi, S., Marzari, F., Massironi, M., Mottola, S., Naletto, G., Oklay, N., Pajola, M., Pommerol, A., Thomas, N., Toth, I., Tubiana, C., Vincent, J.-B., 2015. Temporal morphological changes in the Imhotep region of comet 67P/Churyumov–Gerasimenko. Astron. Astrophys. 583, A36. Gulkis, S., Allen, M., von Allmen, P., Beaudin, G., Biver, N., Bockelée-Morvan, D., Choukroun, M., Crovisier, J., Davidsson, B.J.R., Encrenaz, P., Encrenaz, T., Frerking, M., Hartogh, P., Hofstadter, M., Ip, W.-H., Janssen, M., Jarchow, C., Keihm, S., Lee, S., Lellouch, E., Leyrat, C., Rezac, L., Schloerb, F.P., Spilker, T., 2015. Subsurface properteis and early activity of comet 67P/Churyumov–Gerasimenko. Science 347, aaa0709. Hässig, M., Altwegg, K., Balsiger, H., Bar-Nun, A., Berthelier, J.J., Bieler, A., Bochsler, P., Briois, C., Calmonte, U., Combi, M., De Keyser, J., Eberhardt, P., Fiethe, B., Fuselier, S.A., Galand, M., Gasc, S., Gombosi, T.I., Hansen, K.C., Jackel, A., Keller, H. U., Kopp, E., Korth, A., Kührt, E., Le Roy, L., Mall, U., Marty, B., Mousis, O., Neefs, E., Owen, T., Réme, H., Rubin, M., Sémon, T., Tornow, C., Tzou, C.-Y., Waite, J.H., Wurz, P., 2015. Time variability and heterogeneity in the coma of 67P/ Churyumov–Gerasimenko. Science 347 (6220), aaa0276. Hu, X., Shi, X., Sierks, H., Fulle, M., Blum, J., Keller, H.U., Kührt, E., Davidsson, B., Güttler, C., Gundlach, B., Pajola, M., Bodewits, D., Vincent, J.-B., Oklay, N., Massironi, M., Fornasier, S., Tubiana, C., Groussin, O., Boudreault, S., Höfner, S., Mottola, S., Barbieri, C., Lamy, P.L., Rodrigo, R., Koschny, D., Rickman, H., A’Hearn, M., Agarwal, J., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Da Deppo, V., Debei, S., De Cecco, M., Deller, J., El-Maarry, M.R., Gicquel, A., GutierrezMarques, P., Gutiérrez, P.J., Hofmann, M., Hviid, S.F., Ip, W.-H., Jorda, L., Knollenberg, J., Kovacs, G., Kramm, J.-R., Küppers, M., Lara, L.M., Lazzarin, M., Lopez-Moreno, J.J., Marzari, F., Naletto, G., Thomas, N., 2017. Astron. Astrophys. 604, A114. Huebner, W.F., Markiewicz, W.J., 2000. The temperature and bulk flow speed of a gas effusing or evaporating from a surface into a void after reestablishment of collisional equilibrium. Icarus 148, 594–596. Huebner, W.F., Benkhoff, J., Capria, M.-T., Coradini, A., De Sanctis, C., Orosei, R., Prialnik, D., 2006. Heat and gas diffusion in comet nuclei. ESA Publications Division, Noordwijk, The Netherlands. Ishiguro, M., Sarugaku, Y., Nishihara, S., Nakada, Y., Nishiura, S., Soyano, T., Tarusawa, K., Mukai, T., Kwon, S.M., Hasegawa, S., Usui, F., Ueno, M., 2009. Report on the Kiso cometary dust trail survey. Adv. Space Res. 43, 875–879. Ivanovski, S.L., Della Corte, V., Rotundi, A., Fulle, M., Fougere, N., Bieler, A., Rubin, M., Ivanovska, S., Liuzzi, V., 2017a. Dynamics of non–spherical dust in the coma of 67P/ Churyumov- Gerasimenko constrained by GIADA and ROSINA data. Mon. Not. R. Astron. Soc. 469, S774–S786. Ivanovski, S.L., Zakharov, V.V., Della Corte, V., Crifo, J.-F., Rotundi, A., Fulle, M., 2017b. Dynamics of aspherical dust grains in a cometary atmosphere: I. axially symmetric grains in a spherically symmetric atmosphere. Icarus 282, 333–350. Jorda, L., Gaskell, R., Capanna, C., Hviid, S., Lamy, P., Ďurech, J., Faury, G., Groussin, O., Gutiérrez, P., Jackman, C., Keihm, S.J., Keller, H.U., Knollenberg, J., Kührt, E., Marchi, S., Mottola, S., Palmer, E., Schloerb, F.P., Sierks, H., Vincent, J.-B., A’Hearn, M.F., Barbieri, C., Rodrigo, R., Koschny, D., Rickman, H., Barucci, M.A., Bertaux, J.L., Bertini, I., Cremonese, G., Deppo, V.D., Davidsson, B., Debei, S., De Cecco, M., Fornasier, S., Güttler, C., Kramm, J.R., Lara, L.M., Lazzarin, M., Moreno, J.J. Lopez, Marzari, F., Naletto, G., Oklay, N., Thomas, N., Tubiana, C., Wenzel, K.-P., 2016. The global shape, density and rotation of comet 67P/ Churyumov–Gerasimenko from preperihelion Rosetta/OSIRIS observations. Icarus 277, 257–278. Keller, H.U., Barbieri, C., Lamy, P., Rickman, H., Rodrigo, R., Wenzel, K.-P., Sierks, H., A’Hearn, M.F., Angrilli, F., Angulo, M., Bailey, M.E., Barthol, P., Barucci, M.A., Bertaux, J.-L., Bianchini, G., Boit, J.-L., Brown, V., Burns, J.A., Büttner, I., Castro, J. M., Cremonese, G., Curdt, W., da Deppo, V., Debei, S., de Cecco, M., Dohlen, K., Fornasier, S., Fulle, M., Germerott, D., Gliem, F., Guizzo, G.P., Hviid, S.F., Ip, W.-H., Jorda, L., Koschny, D., Kramm, J.R., Kührt, E., Küppers, M., Lara, L.M., Llebaria, A., López, A., López-Jimenez, A., López-Moreno, J., Meller, R., Michalik, H., Michelena, M.D., Müller, R., Naletto, G., Origné, A., Parzianello, G., Pertile, M., Quintana, C., Ragazzoni, R., Ramous, P., Reiche, K.-U., Reina, M., Rodríguez, J., Rousset, G., Sabau, L., Sanz, A., Sivan, J.-P., Stöckner, K., Tabero, J., Telljohann, U., Thomas, N., Timon, V., Tomasch, G., Wittrock, T., Zaccariotto, M., 2007. OSIRIS – The scientific camera system onboard Rosetta. Space Sci. Rev. 128, 433–506. Keller, H.U., Mottola, S., Davidsson, B., Schröder, S.E., Skorov, Y., Kührt, E., Groussin, O., Pajola, M., Hviid, S.F., Preusker, F., Scholten, F., A’Hearn, M.F., Angrilli, F., Barbieri, C., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Deppo, V.D., Debei, S., Cecco, M.D., Fornasier, S., Fulle, M., Gutiérrez, P.J., Ip, W.-H., Jorda, L., Knollenberg, J., Koschny, D., Kramm, J.R., Küppers, M., Lamy, P., Lara, L.-M., Lazzarin, M., Moreno, J.J.L., Marzari, F., Michalik, H., Naletto, G., Rickman, H., Rodrigo, R., Sabau, L., Sierks, H., Thomas, N., Vincent, J.-B., Wenzel, K.-P., 2015a. Insolation, erosion, and morphology of comet 67P/Churyumov–Gerasimenko. Astron. Astrophys. 583, A34. Keller, H.U., Mottola, S., Skorov, Y., Jorda, L., 2015b. The changing rotation period of comet 67P/Churyumov–Gerasimenko controlled by its activity. Astron. Astrophys. 579, L5. Keller, H.U., Mottola, S., Hviid, S.F., Agarwal, J., Kührt, E., Skorov, Y., Otto, K., Vincent, J.-B., Oklay, N., Schröder, S.E., Davidsson, B., Pajola, M., Shi, X., Bodewits, D., Toth, I., Preusker, F., Scholten, F., Sierks, H., Barbieri, C., Lamy, P., Rodrigo, R., Koschny, D., Rickman, H., A’Hearn, M.F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Deppo, V.D., Debei, S., Cecco, M.D., Deller, J., Fornasier, S., Fulle, M., Groussin, O., Gutiérrez, P.J., Güttler, C., Hofmann, M., Ip, W.-H., Jorda, L., Knollenberg, J., Kramm, J.R., Küppers, M., Lara, L.-M., Lazzarin, M., Lopez-Moreno, J.J., Marzari, F., Naletto, G., Tubiana, C., Thomas, N., 2017. Seasonal mass transfer on the nucleus of comet 67P/ Churyumov–Gerasimenko. Mon. Not. R. Astron. Soc. 469, S357–S371. Kelley, M.S., Wooden, D.H., Tubiana, C., Boehnhardt, H., Woodward, C.E., Harker, D.E., 2009. Spitzer observations of comet 67P/Churyumov–Gerasimenko at 5.5–4.3 AU from the sun. Astron. J. 137, 4633–4642. Kelley, M.S., Lindler, D.J., Bodewits, D., A’Hearn, M.F., Lisse, C.M., Kolokolova, L., Kissel, J., Hermalyn, B., 2013. A distribution of large particles in the coma of Comet 103P/Hartley 2. Icarus 222, 634–652. Klinger, J., 1980. Influence of a phase transition of ice on the heat and mass balance of comets. Science 209, 271–272. Klinger, J., 1981. Some consequences of a phase transition of water ice on the heat balance of comet nuclei. Icarus 47, 320–324. Kossacki, K.J., Markiewicz, W.J., Skorov, Y., Kömle, N.I., 1999. Sublimation coefficient of water ice under simulated cometary–like conditions. Planet. Space Sci. 47 (12), 1521–1530. Kramer, T., Noack, M., 2015. Prevailing dust–transport directions on Comet 67P/ Churyumov–Gerasimenko. Astrophys. J. Lett. 813, L33. Kührt, E., 1984. Temperature profiles and thermal stresses in cometary nuclei. Icarus 60, 512–521. Lai, I.-L., Ip, W.-H., Su, C.-C., Wu, J.-S., Lee, J.-C., Lin, Z.-Y., Liao, Y., Thomas, N., Sierks, H., Barbieri, C., Lamy, P., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., A’Hearn, M.F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Boudreault, S., Cremonese, G., Da Deppo, V., Davidsson, B., Debei, S., De Cecco, M., Deller, J., Fornasier, S., Fulle, M., Groussin, O., Gutiérrez, P.J., Güttler, C., Hofmann, M., Hviid, S.F., Jorda, L., Knollenberg, J., Kovacs, G., Kramm, J.-R., Kührt, E., Küppers, M., Lara, L.M., Lazzarin, M., Lopez Moreno, J.J., Marzari, F., Naletto, G., Oklay, N., Shi, X., Tubiana, C., Vincent, J.-B., 2017. Gas outflow and dust transport of comet 67P/Churyumov–Gerasimenko. Mon. Not. R. Ast. Soc. 462, S533–S546. Lara, L.M., Lowry, S., Vincent, J.B., Gutiérrez, P.J., Rozek, A., La Forgia, F., Oklay, N., Sierks, H., Barbieri, C., Lamy, P.L., Rodrigo, R., Koschny, D., Rickman, H., Keller, H. U., Agarwal, J., Auger, A.-T., A’Hearn, M.F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Besse, S., Bodewits, D., Cremonese, G., Davidsson, B., Da Deppo, V., Debei, S., De Cecco, M., El-Maarry, M.R., Ferri, F., Fornasier, S., Fulle, M., Groussin, O., GutiérrezMarques, P., Güttler, C., Hviid, S.F., Ip, W.-H., Jorda, L., Knollenberg, J., Kovacs, G., Kramm, J.-R., Kührt, E., Küppers, M., Lazzarin, M., Lin, Z.-Y., López-Moreno, J.J., Magrin, S., Marzari, F., Michalik, H., Moissl-Fraund, R., Moreno, F., Mottola, S., Naletto, G., Pajola, M., Pommerol, A., Thomas, N., Sabau, M.D., Tubiana, C., 2015. Large–scale dust jets in the coma of 67P/Churyumov–Gerasimenko as seen by the OSIRIS instrument onboard Rosetta. Astron. Astrophys. 583, A9. Lee, J.-C., Massironi, M., Ip, W.-H., Giacomini, L., Ferrari, S., Penasa, L., El-Maarry, M.R., Pajola, M., Lai, I.-L., Lin, Z.-Y., Ferri, F., Sierks, H., Barbieri, C., Lamy, P., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., A’Hearn, M.F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Da Deppo, V., Davidsson, B., Debei, S., De Cecco, M., Deller, J., Fornasier, S., Fulle, M., Groussin, O., Gutiérrez, P. J., Güttler, C., Hofmann, M., Hviid, S.F., Jorda, L., Knollenberg, J., Kovacs, G., Kramm, J.-R., Kührt, E., Küppers, M., Lara, L.M., Lazzarin, M., Marzari, F., Lopez Moreno, J.J., Naletto, G., Oklay, N., Shi, X., Thomas, N., Tubiana, C., Vincent, J.-B., 2016. Geomorphological mapping of comet 67P/Churyumov–Gerasimenko’s Southern hemisphere. Mon. Not. R. Astron. Soc 462 (Suppl_1), S573–S592. Levison, H.F., Duncan, M.J., 1994. The long–term dynamical behavior of long–period comets. Icarus 108, 18–36. 23 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Luspay-Kuti, A., Hässig, M., Fuselier, S.A., Mandt, K.E., Altwegg, K., Balsiger, H., Gasc, S., Jäckel, A., Le Roy, L., Rubin, M., Tzou, C.-Y., Wurz, P., Mousis, O., Dhooghe, F., Berthelier, J.J., Fiethe, B., Gombosi, T.I., Mall, U., 2015. Composition–dependent outgassing of comet 67P/Churyumov–Gerasimenko from ROSINA/DFMS. Astron. Astrophys. 583, A4. Mekler, Y., Prialnik, D., Podolak, M., 1990. Evaporation from a porous cometary nucleus. Astrophys. J. 356, 682–686. Mottola, S., Arnold, G., Grothues, H.-G., Jaumann, R., Michaelis, H., Neukum, G., Bibring, J.-P., Schröder, S.E., Hamm, M., Otto, K.A., Pelivan, I., Proffe, G., Scholten, F., Tirsch, D., Kreslavsky, M., Remetean, E., Souvannavong, F., Dolives, B., 2015. The structure of the regolith on 67P/Churyumov–Gerasimenko from ROLIS descent imaging. Science 349 (6247), aab0232–1. Özişik, M.N., 1985. Heat transfer. A basic approach. McGraw–Hill, Inc., New York. Pajola, M., Mottola, S., Hamm, M., Fulle, M., Davidsson, B., Güttler, C., Sierks, H., Naletto, G., Arnold, G., Grothues, H.-G., Jaumann, R., Michaelis, H., Bibring, J.P., Barbieri, C., Lamy, P.L., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., A’Hearn, M.F., Barucci, M.A., Bertaux, J.L., Bertini, I., Boudreault, S., Cremonese, G., Da Deppo, V., Debei, S., De Cecco, M., Deller, J., El Maarry, M.R., Feller, C., Fornasier, S., Gicquel, A., Groussin, O., Gutierrez, P.J., Hofmann, M., Hviid, S.F., Ip, W.H., Jorda, L., Knollenberg, J., Kramm, J.R., Kührt, E., Küppers, M., La Forgia, F., Lara, L.M., Lin, Z.Y., Lazzarin, M., Lopez Moreno, J.J., Lucchetti, A., Marzari, F., Massironi, M., Michalik, H., Oklay, N., Pommerol, A., Preusker, F., Scholten, F., Thomas, N., Tubiana, C., Vincent, J.B., 2016. The Agilkia boulders/ pebbles size-frequency distributions: OSIRIS and ROLIS joint observations of 67P surface. Mon. Not. R. Astron. Soc. 462, S242–S252. Pajola, M., Höfner, S., Vincent, J.B., Oklay, N., Scholten, F., Preusker, F., Mottola, S., Naletto, G., Fornasier, S., Lowry, S., Feller, C., Hasselmann, P.H., Güttler, C., Tubiana, C., Sierks, H., Barbieri, C., Lamy, P., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., A’Hearn, M.F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Besse, S., Boudreault, S., Cremonese, G., Da Deppo, V., Davidsson, B., Debei, S., De Cecco, M., Deller, J., Deshapriya, J.D.P., El-Maarry, M.R., Ferrari, S., Ferri, F., Fulle, M., Groussin, O., Gutierrez, P., Hofmann, M., Hviid, S.F., Ip, W.-H., Jorda, L., Knollenberg, J., Kovacs, G., Kramm, J.R., Kührt, E., Küppers, M., Lara, L.M., Lin, Z.Y., Lazzarin, M., Lucchetti, A., Lopez Moreno, J.J., Marzari, F., Massironi, M., Michalik, H., Penasa, L., Pommerol, A., Simioni, E., Thomas, N., Toth, I., Baratti, E., 2017a. The pristine interior of comet 67P revealed by the combined Aswan outburst and cliff collapse. Nature Astron. 1, 0092. Pajola, M., Lucchetti, A., Fulle, M., Mottola, S., Hamm, M., Da Deppo, V., Penasa, L., Kovacs, G., Massironi, M., Shi, X., Tubiana, C., Güttler, C., Oklay, N., Vincent, J.B., Toth, I., Davidsson, B., Naletto, G., Sierks, H., Barbieri, C., Lamy, P.L., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., A’Hearn, M.F., Barucci, M.A., Bertaux, J.L., Bertini, I., Cremonese, G., Debei, S., De Cecco, M., Deller, J., El Maarry, M.R., Fornasier, S., Frattin, E., Gicquel, A., Groussin, O., Gutierrez, P.J., Höfner, S., Hofmann, M., Hviid, S.F., Ip, W.H., Jorda, L., Knollenberg, J., Kramm, J. R., Kührt, E., Küppers, M., Lara, L.M., Lazzarin, M., Lopez Moreno, J.J., Marzari, F., Michalik, H., Preusker, F., Scholten, F., Thomas, N., 2017b. The pebbles/boulders size distributions on Sais: Rosetta’s final landing site on comet 67P/ChuryumovGerasimenko. Mon. Not. R. Astron. Soc. 471, 680–689. Pajola, M., Lee, J.-C., Oklay, N., Hviid, S.F., Penasa, L., Mottola, S., Shi, X., Fornasier, S., Davidsson, B., Giacomini, L., Lucchetti, A., Massironi, M., Vincent, J.B., Bertini, I., Naletto, G., Ip, W.H., Sierks, H., Lamy, P.L., Rodrigo, R., Koschny, D., Keller, H.U., Agarwal, J., Barucci, M.A., Bertaux, J.L., Bodewits, D., Cambianica, P., Cremonese, G., Da Deppo, V., Debei, S., De Cecco, M., Deller, J., El Maarry, M.R., Feller, C., Ferrari, S., Fulle, M., Gutierrez, P.J., Güttler, C., Lara, F., La Forgia, L.M., Lazzarin, M., Lin, Z.-Y., Moreno, J.J. Lopez, Marzari, F., Preusker, F., Scholten, F., Toth, I., Tubiana, C., 2019. Multidisciplinary analysis of the Hapi region located on Comet 67P/Churyumov-Gerasimenko. Mon. Not. R. Astron. Soc. 485 (2), 2139–2154. Pätzold, M., Andert, T., Hahn, M., Asmar, S.W., Barriot, J.-P., Bird, M.K., Häusler, B., Peter, K., Tellmann, S., Grün, E., Weissman, P.R., Sierks, H., Jorda, L., Gaskell, R., Preusker, F., Scholten, F., 2016. A homogeneous nucleus for comet 67P/ChuryumovGerasimenko from its gravity field. Nature 530, 63–65. Pätzold, M., Andert, T.P., Hahn, M., Barriot, J.-P., Asmar, S.W., Häusler, B., Bird, M.K., Tellmann, S., Oschlisniok, J., Peter, K., 2019. The nucleus of comet 67P/ChuryumovGerasimenko – Part I: The global view – nucleus mass, mass–loss, porosity, and implications. Mon. Not. R. Astron. Soc. 483, 2337–2346. Poulet, F., Lucchetti, A., Bibring, J.-P., Carter, J., Gondet, B., Jorda, L., Langevin, Y., Pilorget, C., Capanna, C., Cremonese, G., 2017. Origin of the local structures at the Philae landing site and possible implications on the formation and evolution of 67P/ Churyumov–Gerasimenko. Mon. Not. R. Astron. Soc. 462, S23–S32. Preusker, F., Scholten, F., Matz, K.-D., Roatsch, T., Willner, K., Hviid, S.F., Knollenberg, J., Jorda, L., Gutiérrez, P.J., Kührt, E., Mottola, S., A’Hearn, M.F., Thomas, N., Sierks, H., Barbieri, C., Lamy, P., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Da Deppo, V., Davidsson, B., Debei, S., De Cecco, M., Fornasier, S., Fulle, M., Groussin, O., Güttler, C., Ip, W.-H., Kramm, J.R., Küppers, M., Lara, L.M., Lazzarin, M., Lopez Moreno, J.J., Marzari, F., Michalik, H., Naletto, G., Oklay, N., Tubiana, C., Vincent, J.-B., 2015. Shape model, reference system definition, and cartographic mapping standards for comet 67P/Churyumov-Gerasimenko – Stereo–photogrammetric analysis of Rosetta/OSIRIS image data. Astron. Astrophys. 583, A33. Preusker, F., Scholten, F., Matz, K.-D., Roatsch, T., Hviid, S.F., Mottola, S., Knollenberg, J., Kührt, E., Pajola, M., Oklay, N., Vincent, J.-B., Davidsson, B., A’Hearn, M.F., Agarwal, J., Barbieri, C., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Da Deppo, V., Debei, S., De Cecco, M., Fornasier, S., Fulle, M., Groussin, O., Gutiérrez, P.J., Güttler, C., Ip, W.-H., Jorda, L., Keller, H.U., Koschny, D., Kramm, J.R., Küppers, M., Lamy, P., Lara, L.M., Lazzarin, M., Lopez Moreno, J.J., Marzari, F., Massironi, M., Naletto, G., Rickman, H., Rodrigo, R., Sierks, H., Thomas, N., Tubiana, C., 2017. The global meter–level shape model of comet 67P/Churyumov–Gerasimenko. Astron. Astrophys. 607, L1. Prialnik, D., 1992. Crystallization, sublimation, and gas release in the interior of a porous comet nucleus. Astrophys. J. 388, 196–202. Ratke, L., Kochan, H., 1989. Fracture mechanical aspects of dust emission processes from a model comet surface. In: Physics and mechanics of cometary materials, pp. 121–128 (ESA SP–302). Reach, W.T., Kelley, M.S., Sykes, M.V., 2007. A survey of debris trails from short–period comets. Icarus 191, 298–322. Robie, R.A., Hemingway, B.S., Takei, H., 1982. Heat capacities and entropies of Mg2SiO4, Mn2SiO4, and CO2SiO4 between 5 and 380 K. American Mineralogist 67, 470–482. Rotundi, A., Sierks, H., Corte, V.D., Fulle, M., Gutiérrez, P.J., Lara, L., Barbieri, C., Lamy, P.L., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., López-Moreno, J.J., Accolla, M., Agarwal, J., A’Hearn, M.F., Altobelli, N., Angrilli, F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Bodewits, D., Bussoletti, E., Colangeli, L., Cosi, M., Cremonese, G., Crifo, J.-F., Deppo, V.D., Davidsson, B., Debei, S., Cecco, M.D., Esposito, F., Ferrari, M., Fornasier, S., Giovane, F., Gustafson, B., Green, S.F., Groussin, O., Grün, E., Güttler, C., Herranz, M.L., Hviid, S.F., Ip, W., Ivanovski, S., Jerónimo, J.M., Jorda, L., Knollenberg, J., Kramm, R., Kührt, E., Küppers, M., Lazzarin, M., Leese, M.R., López-Jiménez, A.C., Lucarelli, F., Lowry, S.C., Marzari, F., Epifani, E.M., McDonnell, J.A.M., Mennella, V., Michalik, H., Molina, A., Morales, R., Moreno, F., Mottola, S., Naletto, G., Oklay, N., Ortiz, J.L., Palomba, E., Palumbo, P., Perrin, J.-M., Rodríguez, J., Sabau, L., Snodgrass, C., Sordini, R., Thomas, N., Tubiana, C., Vincent, J.-B., Weissman, P., Wenzel, K.-P., Zakharov, V., Zarnecki, J.C., 2015. Dust measurements in the coma of comet 67P/Churyumov–Gerasimenko inbound to the Sun between 3.7 and 3.4 AU. Science 347 (6220), aaa3905. Sekanina, Z., 1981. Rotation and precession of cometary nuclei. Ann. Rev. Earth Planet. Sci. 9, 113–145. Shi, X., Hu, X., Sierks, H., Güttler, C., A’Hearn, M., Blum, J., El-Maarry, M.R., Kührt, E., Mottola, S., Pajola, M., Oklay, N., Fornasier, S., Tubiana, C., Keller, H.U., Vincent, J.B., Bodewits, D., Höfner, S., Lin, Z.-Y., Gicquel, A., Hofmann, M., Barbieri, C., Lamy, P.L., Rodrigo, R., Koschny, D., Rickman, H., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Da Deppo, V., Davidsson, B., Debei, S., De Cecco, M., Fulle, M., Groussin, O., Gutiérrez, P.J., Hviid, S.F., Ip, W.-H., Jorda, L., Knollenberg, J., Kovacs, G., Kramm, J.-R., Küppers, M., Lara, L.M., Lazzarin, M., Lopez-Moreno, J.J., Marzari, F., Naletto, G., Thomas, N., 2016. Sunset jets observed on comet 67P/Churyumov–Gerasimenko sustained by subsurface thermal lag. Astron. Astrophys. 586, A7. Shoshany, Y., Prialnik, D., Podolak, M., 2002. Monte Carlo modeling of the thermal conductivity of porous cometary ice. Icarus 157, 219–227. Sierks, H., Barbieri, C., Lamy, P.L., Rodrigo, R., Koschny, D., Rickman, H., Keller, H.U., Agarwal, J., A’Hearn, M.F., Angrilli, F., Auger, A.-T., Barucci, M.A., Bertaux, J.-L., Bertini, I., Besse, S., Bodewits, D., Capanna, C., Cremonese, G., Deppo, V.D., Davidsson, B., Debei, S., Cecco, M.D., Ferri, F., Fornasier, S., Fulle, M., Gaskell, R., Giacomini, L., Groussin, O., Gutierrez-Marques, P., Gutiérrez, P.J., Güttler, C., Hoekzema, N., Hviid, S.F., Ip, W.-H., Jorda, L., Knollenberg, J., Kovacs, G., Kramm, J.-R., Kührt, E., Küppers, M., Forgia, F.L., Lara, L.M., Lazzarin, M., Leyrat, C., Moreno, J.J.L., Magrin, S., Marchi, S., Marzari, F., Massironi, M., Michalik, H., Moissl, R., Mottola, S., Naletto, G., Oklay, N., Pajola, M., Pertile, M., Preusker, F., Sabau, L., Scholten, F., Snodgrass, C., Thomas, N., Tubiana, C., Vincent, J.-B., Wenzel, K.-P., Zaccariotto, M., Pätzold, M., 2015. On the nucleus structure and activity of comet 67P/Churyumov-Gerasimenko. Science 347 (6220), aaa1044. Skorov, Y.V., Rickman, H., 1998. Simulation of gas flow in a cometary Knudsen layer. Planet. Space. Sci. 46 (8), 975–996. Skorov, Y.V., Rickman, H., 1999. Gas flow and dust acceleration in a cometary knudsen layer. Planet. Space. Sci. 47 (8), 935–949. Soderblom, L.A., Becker, T.L., Bennett, G., Boice, D.C., Britt, D.T., Brown, R.H., Buratti, B.J., Isbell, C., Giese, B., Hare, T., Hicks, M.D., Howington-Kraus, E., Kirk, R. L., Lee, M., Nelson, R.M., Oberst, J., Owen, T.C., Rayman, M.D., Sandel, B.R., Stern, S.A., Thomas, N., Yelle, R.V., 2002. Observations of Comet 19P/Borrelly by the Miniature Integrated Camera and Spectrometer aboard Deep Space 1. Science 296, 1087–1091. Sykes, M.V., Walker, R.G., 1992. Cometary dust trails. I - Survey. Icarus 95, 180–210. Takigawa, A., Furukawa, Y., Kimura, Y., Davidsson, B., Nakamura, T., 2019. Exposure experiments of amorphous silicates and organcis to cometary ice and vapor analogs. Astrophys. J. 881, 27. Tancredi, G., Rickman, H., Greenberg, J.M., 1994. Thermochemistry of cometary nuclei. I: The Jupiter family case. Astron. Astrophys. 286, 659–682. Tenishev, V., Combi, M.R., Rubin, M., 2011. Numerical simulations of dust in a cometary coma: Application to comet 67P/Churyumov–Gerasimenko. Astrophys. J. 732, 104. Thomas, N., Davidsson, B., El-Maarry, M., Fornasier, S., Giacomini, L., Berna, A.G., Hviid, S., Ip, W.-H., Jorda, L., Keller, H., Knollenberg, J., Kührt, E., Forgia, F.L., Lai, I., Liao, Y., Marschall, R., Massironi, M., Mottola, S., Pajola, M., Poch, O., Pommerol, A., Preusker, F., Scholten, F., Su, C., Wu, J., Vincent, J.-B., Sierks, H., Barbieri, C., Lamy, P., Rodrigo, R., Rickman, H., Koschny, D., A’Hearn, M.F., Barucci, M.A., Bertaux, J.-L., Bertini, I., Cremonese, G., Deppo, V.D., Debei, S., Fulle, M., Groussin, O., Gutiérrez, P., Kramm, J.-R., Küppers, M., Lara, L.M., Lazzarin, M., Moreno, J.J.L., Marzari, F., Michalik, H., Naletto, G., Güttler, C., 2015a. Re–distribution of particles across the nucleus of comet 67P/ Churyumov–Gerasimenko. Astron. Astrophys. 583, A17. Thomas, N., Sierks, H., Barbieri, C., Lamy, P.L., Rodrigo, R., Rickman, H., Koschny, D., Keller, H.U., Agarwal, J., A’Hearn, M.F., Angrilli, F., Auger, A.-T., Barucci, M.A., Bertaux, J.-L., Bertini, I., Besse, S., Bodewits, D., Cremonese, G., Deppo, V.D., 24 B.J.R. Davidsson et al. Icarus 354 (2021) 114004 Davidsson, B., Cecco, M.D., Debei, S., El-Maarry, M.R., Ferri, F., Fornasier, S., Fulle, M., Giacomini, L., Groussin, O., Gutiérrez, P.J., Güttler, C., Hviid, S.F., Ip, W.H., Jorda, L., Knollenberg, J., Kramm, J.-R., Kührt, E., Küppers, M., Forgia, F.L., Lara, L.M., Lazzarin, M., Moreno, J.J.L., Magrin, S., Marchi, S., Marzari, F., Massironi, M., Michalik, H., Moissl, R., Mottola, S., Naletto, G., Oklay, N., Pajola, M., Pommerol, A., Preusker, F., Sabau, L., Scholten, F., Snodgrass, C., Tubiana, C., Vincent, J.-B., Wenzel, K.-P., 2015b. The morphological diversity of comet 67P/ Churyumov–Gerasimenko. Science 347 (6220), aaa0440. Tuzzolino, A.J., Economou, T.E., Clark, B.C., Tsou, P., Brownlee, D.E., Green, S.F., McDonnell, J.A.M., McBride, N., Colwell, M.T.S.H., 2004. Dust measurements in the coma of Comet 81P/Wild 2 by the Dust Flux Monitor Instrument. Science 304, 1776–1780. Veverka, J., Klaasen, K., A’Hearn, M., Belton, M., Brownlee, D., Chesley, S., Clark, B., Economou, T., Farquhar, R., Green, S.F., Groussin, O., Harris, A., Kissel, J., Li, J.-Y., Meech, K., Melosh, J., Richardson, J., Schultz, P., Silen, J., Sunshine, J., Thomas, P., Bhaskaran, S., Bodewits, D., Carcich, B., Cheuvront, A., Farnham, T., Sackett, S., Wellnitz, D., Wolf, A., 2013. Return to Comet Tempel 1: Overview of Stardust–NExT results. Icarus 222, 424–435. Werner, R.A., Scheeres, D.J., 1997. Exterior gravitation of a polyhedron derived and compared with harmonic and Mascon gravitation representations of Asteroid 4769 Castalia. Cel. Mech. Dyn. Astron. 65 (3), 313–344. Wright, I.P., Sheridan, S., Barber, S.J., Morgan, G.H., Andrews, D.J., Morse, A.D., 2015. CHO–bearing organic compounds at the surface of 67P/Churyumov–Gerasimenko revealed by Ptolemy. Science 349 (6247), aab0673. Yomogida, K., Matsui, T., 1983. Physical properties of ordinary chondrites. J. Geophys. Res. 88 (B11), 9513–9533. Ytrehus, T., 1977. Theory and experiments on gas kinetics in evaporation. In: Potter, J.L. (Ed.), Rarefied gas dynamics, Progress in Astronautics and Aeronautics, vol. 51. The American Institute of Aeronautics and Astronautics, Inc., Washington, pp. 1197–1212. 25