Implications for Ice Stability and Particle Ejection From High‐Resolution Temperature Modeling of Asteroid (101955) Bennu

The finding by the OSIRIS‐REx (Origins, Spectral Interpretation, Resource Identification, and Security‐Regolith Explorer) mission that its target (101955) Bennu is an active asteroid has raised questions as to whether the observed particle ejection events are driven by temperature. To investigate sublimation of water ice and rock thermal fracture as possible temperature‐driven causes, we modeled the global temperatures of Bennu and searched for correlations with the identified ejection points on the asteroid surface. We computed temperatures with the Advanced Thermophysical Model and the 75‐cm‐resolution global shape model of Bennu derived by the OSIRIS‐REx mission. We find that ~1,856 m2 of Bennu's polar regions have orbit‐averaged temperatures that are sufficiently cold to enable water ice, if buried within the top few meters of the surface, to remain stable over geological timescales. Millimeter thick layers of surface water ice are also stable over ~103‐year timescales within polar centimeter‐scale cold traps. However, we do not find evidence of conditions enabling ice stability in the warmer equatorial regions, where ejection events have been observed, implying that sublimation of water ice is not the cause of particle ejection. Conversely, rock thermal fracture remains a possible mechanism of particle ejection. We find high amplitudes of diurnal temperature variation, a proxy for the efficacy of thermal fracturing, at all latitudes on Bennu due to its extreme ruggedness. Therefore, if rock thermal fracture is the mechanism, particles could be ejected from any latitude, which is consistent with the continued observations of particle ejection by OSIRIS‐REx.


Introduction
The OSIRIS-REx (Origins, Spectral Interpretation, Resource Identification, and Security-Regolith Explorer) spacecraft arrived at near-Earth asteroid (NEA) (101955) Bennu in December 2018 (Lauretta et al. 2019).Although the primary goal to return at least 60 g of material from the surface will not be completed until 2023, initial observations of Bennu have enriched our understanding of the surfaces of NEAs and carbonaceous asteroids.Most notable of these discoveries is that Bennu is an active asteroid (Lauretta and Hergenrother et al. 2019).
The navigation camera NavCam 1 (Bos et al. 2018) has observed hundreds of particles actively ejected from the surface and a smaller number in short-term orbits around Bennu.
Confidential manuscript submitted to JGR: Planets During each of the three largest detected events in January and February 2019, ≥60 particles were observed creating a spray-like pattern from the surface.From the particle trajectories and velocity distributions, Lauretta and Hergenrother et al. (2019) triangulated the ejection sites and concluded that the ejection events were impulsive (i.e.all particles ejected at approximately the same time).The first event was sourced from a high southern latitude (~60 to 75° S) and the other two from low-to-mid northern latitudes (~20° N), and all three occurred at late afternoon local solar times (between 15:22 and 18:05).The minor events (i.e.detections of one or a few particles at a time) are harder to trace back to their origins, but appear to occur at all local times of day (Chesley et al. submitted; Pelgrift et al. in press).Lauretta and Hergenrother et al. (2019) explored multiple hypotheses to explain the ejection events.They found that the particle size and velocity distributions were incompatible with rotational disruption of parts of the surface and with electrostatic lofting of particles, and that Bennu's high surface temperatures and lack of spectral evidence for H 2 O argue against sublimation of ice as a source.Other mechanisms that are tied to temperature cycling of the surface, namely thermal fracturing and volatile release by dehydration of phyllosilicate minerals in the surface rocks, are viable explanations.Meteoroid impacts could also eject particles from the surface, and secondary impacts of ejected particles could explain at least some of the smaller events.More detailed investigations into the electrostatic lofting and meteoroid impact mechanisms are addressed in Hartzell et al. (submitted) and Bottke et al. (submitted), respectively.
Although Lauretta and Hergenrother et al. (2019) did not favor volatile sublimation as a candidate mechanism based on global-scale evidence, volatile materials at and immediately below the surface have been detected in some unexpected places in the Solar System, including Mercury (e.g.Slade et al. 1992), the Moon (e.g.Siegler et al. 2016), Ceres (e.g.Ermakov et al. 2017), and Themis (e.g.Rivkin and Emery 2010;Campins et al. 2010).Rough topography on objects with low obliquity can create cold traps in shadows at high latitudes (e.g.Salvail and Fanale 1994;Hayne and Aharonson 2015), and unbound volatiles can migrate to these cold traps (e.g.Paige et al. 2013).Water molecules released from hydrated minerals by a combination of mechanical and thermal cycling or created by solar wind bombardment of the surface might accumulate in such cold traps.If heated suddenly, the sublimation and expansion of water ice could contribute to the particle ejection events observed on Bennu.The observation that the first detected event originated at high latitudes hinted at this explanation, but later lower-latitude events may be more difficult to explain this way.Determining the locations of potential cold traps on Bennu is therefore an important step to quantitative assessments of volatile sublimation as a potential source for any of the ejection events.
Thermally-induced fatigue in rock is driven by stress fields induced in response to diurnal cycles imposed by the rotations of planetary bodies (e.g.Viles et al. 2010;Eppes et al. 2015) and is thought to be an important driver of rock breakdown on some airless planetary surfaces (e.g.Dombard et al. 2010;Jewitt and Li 2010;Attree et al. 2018).Field and laboratory studies have demonstrated that terrestrial and chondritic meteorites subjected to thermal cycling experience crack propagation and, in some cases, disaggregation of material from rock (e.g.Delbó et al. 2014;Collins and Stock 2016;Eppes et al. 2016).Recent works simulating stress fields in boulders have shown that stresses at different times and locations can drive the development of fractures leading to effects such as exfoliation, surface disaggregation, and through-going fractures (Molaro et al. 2015(Molaro et al. , 2017)).The magnitude and timing of these stress fields depend on both the boulder composition and size, as well as the amplitude of diurnal temperature variation Confidential manuscript submitted to JGR: Planets (Molaro et al. 2017;El Mir et al. 2019;Hamm et al. 2019).The latter may be used as a proxy for the efficacy of thermal fatigue, which motivates the need to quantify temperature variations across Bennu's surface to better understand its potential role in particle ejections.
Through ground-based and OSIRIS-REx observations, the physical, orbital, and rotational properties of Bennu that are important for computing surface and sub-surface temperatures are better characterized than for most asteroids.Bennu has a spheroidal shape with an equatorial bulge (Nolan et al. 2013;Barnouin et al. 2019), very low obliquity (i.e.177.6 ± 0.1°), and extremely well-determined rotation period (Barnouin et al. 2019).The global (i.e.hemispherical-scale) thermal inertia of Bennu is relatively low (i.e.350 ± 20 J m -2 K -1 s -1/2 ; Emery et al. 2014;DellaGiustina and Emery et al. 2019), despite the dominance of boulders larger than the diurnal thermal skin depth on the surface (Walsh et al. 2019;DellaGiustina and Emery et al. 2019), and is fairly uniform with rotational phase.Hamilton et al. (2019) report the spectral detection of hydrated minerals on the surface of Bennu, suggesting that a reservoir of volatiles was present at some point in Bennu's past.
In this paper, we investigate sublimation of water ice and rock thermal fracture as possible temperature-driven causes of particle ejection by modeling the surface and sub-surface temperature distributions of Bennu, and searching for correlations with the identified particle radiant points (the locations of particle ejection).Section 2 describes the Advanced Thermophysical Model that we used to model temperatures on Bennu along with the temperature criteria that we used to assess the two potential mechanisms.Section 3 describes the temperature modeling results for the high resolution global shape model of Bennu, the particle radiant points identified in Lauretta and Hergenrother et al. (2019), and unresolved small-scale surface roughness.Section 4 provides a discussion of the temperature modeling results, and Section 5 provides a summary with conclusions.

The Advanced Thermophysical Model
To investigate sublimation of water ice and thermal fracturing as possible temperaturedriven causes of particle ejection from Bennu requires evaluation of Bennu's surface and subsurface temperature distribution.As described in sections 2.2 and 2.3, the stability of surface water ice is primarily governed by the maximum surface temperature experienced by a particular facet (surface element) at any given time during the orbit of Bennu, T MAX .Similarly, the stability of sub-surface water ice is primarily governed by the facet's orbital average surface temperature, T AVG .Additionally, the magnitude of thermally-induced stress is approximately proportional to the amplitude of the diurnal temperature variation that the facet experiences during successive Bennu rotations, ΔT.Calculation of these parameters therefore requires a suitable thermophysical model to predict temperatures at any given point and time on the surface and in the sub-surface of Bennu.For this purpose, we use the Advanced Thermophysical Model (ATPM) developed by Rozitis and Green (2011, 2012, 2013).
To model the surface and sub-surface temperatures of an asteroid as a function of time, the ATPM ingests a shape model in the triangular facet formalism and solves the 1D heat conduction equation with a suitable surface boundary condition for each triangular facet.Lateral heat conduction is ignored because only facets larger than the thermal skin depth (typically a few centimeters) are considered.For temperature T, time t, and depth z, 1D heat conduction is described by Confidential manuscript submitted to JGR: Planets where k is the thermal conductivity, ρ is the material density, and C P is the heat capacity.The latter three properties are assumed to be constant with temperature and depth, and can be combined into the single parameter known as thermal inertia, Γ, via For conservation of energy between incoming and outgoing radiation, the ATPM takes into account direct solar illumination, projected shadows, multiple scattering of sunlight, and self-heating effects within its surface boundary condition.This surface boundary condition is given by where A B is the Bond albedo, F  is the integrated solar flux at 1 AU (i.e.1367 W m -2 ), r H (t) is the heliocentric distance of the asteroid in AU at time t, ε is the bolometric emissivity, and σ is the Stefan-Boltzmann constant.S(t) is a function that determines whether the facet is shadowed at time t, and ψ(t) is a function that returns the cosine of the Sun illumination angle of the facet.
Finally, F SCAT (t) and F RAD (t) are functions that evaluate the total multiple-scattered sunlight and the total thermal emission that are imposed on the facet from neighboring interfacing facets, respectively.
For a given pole orientation, orbital position, and rotational phase of the asteroid, the ATPM computes the illumination geometry for each facet specified in the asteroid shape model.
Projected shadows are determined by ray-triangle intersection tests; they return S(t) = 1 if a particular facet is shadowed, and S(t) = 0 if it is not shadowed, at the specified time.If more precision is required (e.g. for high illumination angles), then the model splits the facet of interest into a set of smaller facets on which the ray-triangle intersection tests are repeated.A fractional value of S(t) that indicates how much of the original facet is covered in a projected shadow is then returned.
F SCAT (t) and F RAD (t) are calculated by using viewfactors.In particular, the viewfactor f i,j specifies the fractional amount of radiation that is reflected or emitted by facet i and is transferred to facet j when assuming Lambertian reflection and emission (Lagerros 1998).
F SCAT,i (t) for facet i is then calculated by and F RAD,i (t) by Viewfactors are pre-computed for each shape model investigated and are stored in a lookup table for use by the ATPM.In the model, F SCAT (t) is computed by using several Gauss-Seidel iterations to evaluate multiple bounces of reflected sunlight, and F RAD (t) is computed from the surface temperatures determined by the previous model revolution.
The asteroid shape models that are typically input into the ATPM have facets that are at least several meters in size and are therefore much larger than the diurnal thermal skin depth.
Unresolved surface roughness, which occurs at spatial scales between the diurnal thermal skin depth and the shape model facet size, has a tendency to re-radiate absorbed sunlight back towards the Sun, an effect known as thermal-infrared beaming (Lagerros 1998;Rozitis and Confidential manuscript submitted to JGR: Planets Green 2011).To model the thermal-infrared beaming effect, the ATPM also ingests an additional topography model to represent the unresolved surface roughness.For the primary temperature modeling in this work, we adopt spherical-section craters of the form specified in Spencer (1990) for simplicity.The ATPM is also capable of using more complex roughness models, and we explored a range of fractal roughness models when investigating the potential presence of smallscale cold traps.The ATPM places the topography model at the location and orientation of each of the shape model facets, and equations ( 1) to ( 4) given above are also applied to each of the topography model's facets.
After specifying suitable time and depth domains, a finite difference method is used to solve the 1D heat conduction equation given in (1), and a suitable number of Newton-Raphson iterations is used to solve the surface boundary condition given by equation ( 2).Typically, the ATPM only considers diurnal temperature variations for a specified fixed position in an asteroid's orbit.Following Spencer et al. (1989), equations ( 1) and ( 2) are normalized in the model to the diurnal thermal skin depth, l, given by where P is the asteroid rotation period.Sub-surface temperatures are computed down to a maximum depth of eight diurnal thermal skin depths, which are resolved by seven depth steps per diurnal thermal skin depth.Additionally, each asteroid rotation is resolved by 650 time steps, and zero temperature gradient is also assumed at the model's maximum depth to give a required internal boundary condition.These time and depth domain parameters were chosen to ensure that the ATPM accurately models the diurnal component of the Yarkovsky effect (Rozitis et al. 2013), which it was recently verified to do for Bennu by having correctly estimated its mass from orbital drift measurements (Chesley et al. 2014).For instance, the gravitational parameter of Bennu was determined to be 4.9 ± 0.1 and 4.892 ± 0.006 m 3 s −2 from the ATPM and OSIRIS-REx radio science analyses, respectively (Scheeres et al. 2019).To initialize the model, the initial facet temperatures are set to the rotational average temperature obtained when assuming instantaneous equilibrium with direct solar illumination.The model is then run for several tens of asteroid rotations until changes in surface temperature between successive revolutions diminish to less than 10 -3 K.
Seasonal variations in temperature are also important in assessing the stability of water ice and the rates of thermal fracturing on an asteroid.These seasonal variations arise from orbital variations in heliocentric distance due to the orbital eccentricity, from orbital variations of the sub-solar latitude due to the obliquity of the asteroid, and from the thermal lag induced by the seasonal thermal wave.Ideally, the model time and depth domains specified earlier should be set up in a way to allow capture of both the diurnal and seasonal thermal waves.However, this is computationally expensive, especially when unresolved surface roughness is also considered.To simplify modeling of the seasonal temperature variations, we only consider seasonal changes in illumination geometry and ignore the seasonal thermal wave by running a series of independent diurnal thermal models around the orbit of the asteroid.This is an acceptable approximation because T MAX depends primarily on the maximum irradiance imposed on a facet during a rotation; T AVG depends primarily on the total irradiance imposed on a facet during an orbit; and ΔT depends primarily on the difference between the maximum and minimum temperatures experienced by a facet during a rotation.All of these are adequately captured by running a series of diurnal thermal models around the orbit of an asteroid.
Confidential manuscript submitted to JGR: Planets The seasonal thermal wave does dictate the precise orbital timing at which deep subsurface layers reach their maximum temperature, but this information is not required for our investigation.For assessing the stability of sub-surface water ice, we only need to know T AVG , which can be calculated from the surface temperatures because it is a property that does not vary with depth.Tests with a seasonal thermal model (Vokrouhlický and Farinella 1998) demonstrate that the seasonal thermal wave only introduces a small orbital modulation (i.e.~1 K) to the surface temperatures derived by the series of diurnal thermal models for a Bennu-like thermal inertia and orbit.Therefore, T MAX , T AVG , and ΔT are calculated with sufficient accuracy in this approximation.The only disadvantage with this approach is that the series of diurnal thermal models tend to underestimate facet surface temperature during seasonal shadows.However, as such seasonal shadows only occur briefly during Bennu's orbit, they do not meaningfully affect our calculations of T MAX , T AVG , and ΔT.
Finally, in our thermophysical modeling of Bennu, we assume input parameters that were derived by DellaGiustina and Emery et al. ( 2019) from the OSIRIS-REx Approach-phase infrared observations.In particular, we adopt a thermal inertia of 350 J m -2 K -1 s -1/2 , a Bond albedo of 0.016, and a bolometric emissivity of 0.9 in the ATPM, and we run the ATPM for the derived shape model of Bennu (Barnouin et al. 2019) using its pole orientation, λ = 68.9°and β = -83.0°,and rotation period, 4.296 hr (Lauretta et al. 2019).We do not consider the relatively small uncertainties of the input properties within the temperature modeling of Bennu because the shape and topography have the greatest influence on the calculation of predicted temperature.
Similarly, we also do not consider spatial variations in the input properties, but any variations are expected to be small because of the lack of rotational variability seen in the Approach-phase infrared observations (DellaGiustina and Emery et al. 2019).From equation ( 5), this thermal inertia gives diurnal and seasonal thermal skin depths of ~2 cm and ~1 m, respectively, when a heat capacity of 750 J kg -1 K -1 and a material density of 1190 kg m -3 (the bulk density of Bennu; Scheeres et al. 2019) are assumed.Therefore, our computational domain for diurnal temperature modeling extends to a maximum physical depth of ~16 cm, and the resulting T AVG are applicable to sub-surface layers up to several meters deep.From equations (3) and (4), the low Bond albedo of Bennu causes F SCAT (t) to be a factor of ~40 smaller than F RAD (t) for midday illumination, so self-heating dominates the radiative exchange of energy between facets in our modeling.

Assessing the stability of water ice
Volatile stability is calculated based on two main criteria: the sublimation vapor pressure of a volatile molecule at a given temperature, and the ability for that molecule to diffuse through overlying material.For water ice, its sublimation rate when directly exposed to vacuum, E, can be calculated by the standard formula where p V (T) is the equilibrium vapor pressure of ice at temperature T, m is the molecular mass of water, and R is the universal gas constant (Estermann 1955).Therefore, water will leave the surface at a rate of ~10 -9 kg m -2 yr -1 at a temperature of approximately 100 K (e.g. Watson et al. 1961;Schorghofer and Taylor 2007), which makes it geologically stable.For 125 K, that rate rises to ~10 -5 kg m -2 yr -1 , and to about ~10 -2 kg m -2 yr -1 for 140 K.This exponential increase in the loss rate over a small range in temperatures makes water ice a precise marker of past temperature maxima.As we want to evaluate the possible presence of centimeter-scale cold traps, in addition to meter-scale cold traps, on Bennu, we use the T MAX < 131 K criterion of Confidential manuscript submitted to JGR: Planets Jewitt and Guilbert-Lepoutre (2012) for ~10 3 -year stability of millimeter-thick layers of surface water ice.Such a short stability period, in geological terms, would require recharge of the surface water ice from sub-surface reservoirs, and the potential presence of such reservoirs on Bennu are evaluated separately using the temperature criterion described next.
A small cover of regolith can also preserve volatiles at much higher temperatures (e.g.Schorghofer and Taylor 2007).This regolith layer can both insulate sub-surface volatiles from high temperatures and provide a tortuous diffusion pathway for molecules that do sublimate in the sub-surface.Using a simple estimation of residence time on regolith grain surfaces as a function of temperature, Schorghofer and Taylor (2007) show that just 10 cm of regolith cover can increase the 10 -9 kg m -2 yr -1 loss rate temperature of water ice to ~120 K.The maximum average temperature at which water ice has been modeled to be stable over geological timescales at any depth under vacuum is 145 K, above which the loss rates are too high for water ice to remain under any thickness of regolith cover (Schorghofer 2008).Therefore, to evaluate the potential presence of water ice buried within the top few meters of the surface of Bennu, we use this criterion of T AVG < 145 K. To calculate T AVG , we also follow Schorghofer (2008) by averaging the facet surface temperatures around the orbit of Bennu.

Assessing the efficacy of thermal fracturing
As described in section 1, the breakdown of rocks by thermal cycling is a complex process involving the propagation of cracks by thermal stresses induced by spatial and temporal temperature gradients.For an object with constant material properties, the magnitude of stress induced by thermal cycling is directly proportional to the amplitude of the temperature variation.
However, the resulting crack propagation varies non-linearly with stress (e. fracturing in a relative sense on NEA (162173) Ryugu.They assumed that regolith generation by thermal fracturing was proportional to the maximum ΔT experienced anywhere on the surface of Ryugu, and they subsequently determined the latitudinal distribution of ΔT by assuming a spherical shape for Ryugu.They find that the generation of regolith can take place within a surprisingly wide band around the equator of Ryugu.For our Bennu investigation, we also assume that if the particle ejection events were driven by thermal fractures, then they would have originated from areas of relatively high ΔT.However, we incorporate the measured shape of Bennu rather than assume a sphere.Therefore, to assess the relative rates of thermal fracturing, we compare the modeled ΔT for the particle ejection sites against the maximum ΔT we find anywhere on Bennu.We do not investigate the mechanism(s) by which thermal fractures eject particles, as this is addressed in Molaro et al. (submitted).In particular, that work demonstrates how particles may be ejected during exfoliation or other thermal cracking events via the release of stored thermal strain energy in boulders, and it estimates the resultant particle ejection speeds and sizes.Our work here informs where such processes may be active on Bennu due to the thermal cycle.
Confidential manuscript submitted to JGR: Planets 3 Temperature modeling results

High-resolution global shape model
To provide context for the individual particle ejection sites, we first assessed the global temperatures of Bennu by running the ATPM on the global shape model produced by stereophotoclinometry as described in Barnouin et al. (2019).In particular, we used global shape model v19 produced with ~3 million 75-cm-sized facets to ensure accurate assessments of water ice stability at meter-scale depths.The ATPM was run for 36 positions, equally spaced in true anomaly, around the orbit of Bennu, and the modeled surface temperatures were stored in lookup tables to allow comparison with the temperature criteria given in sections 2.2 and 2.3.
Although Bennu has an orbit with an eccentricity of 0.2 and a rotation pole with an obliquity of ~177.6°, the orbital variations in heliocentric distance and sub-solar latitude (see Figure 1) give rise to modest seasonal variations in temperature.For an equatorial surface element orientated with zero tilt, T MAX and ΔT vary respectively from 390 and 140 K at perihelion to 303 and 82 K at aphelion (Figure 2).Therefore, for the equatorial and mid-latitude regions on Bennu, the perihelion temperatures matter the most for assessing the stability of surface water ice, and for the peak efficacies of thermal fracturing.However, the small orbital variations in sub-solar latitude, which are offset in orbital phase to the heliocentric distance variations (Figure 1), make it necessary to evaluate T MAX and ΔT around the entire orbit of Bennu for the polar regions.Therefore, we sought the highest T MAX and ΔT that each facet attained at any point during Bennu's orbit, which was in addition to the orbital averaging of the facet surface temperatures needed to calculate T AVG .where millimeter-thick layers of surface water ice could persist at 75-cm spatial scales.However, in small pockets located near both poles, T AVG was less than 145 K and therefore would allow water ice buried within top few meters of the surface to remain stable over geological timescales.
The latitudinal distribution of these potential sub-surface water ice pockets is shown in Figure 4; they are more common at Bennu's south pole than at its north pole.This asymmetry arises because the south pole receives less solar illumination, on average, than the north pole owing to the timing of seasons for the two poles in relation to the orbital variation of heliocentric distance (Figure 1).For instance, as southern winter occurs near aphelion, it is deeper and longer-lasting than northern winter due to the greater heliocentric distance and slower orbital motion of Bennu at this part of the orbit.In total, ~1856 m 2 or ~0.2% of Bennu's surface has the potential to harbour sub-surface water ice, but all of these locations are at high latitude (above ~60° N and S).If the particle ejection events were driven by water ice sublimation, they should all originate from near Bennu's poles.Relatively high ΔT (~140 K) values were obtained at all latitudes due to the extreme ruggedness of Bennu's shape (Figures 3 and 4), which suggests that thermal fatigue is not limited to the equatorial region of Bennu.This is similar to the finding of Hamm et al. (2019) for Ryugu, but we find that the presence of boulders increases ΔT at latitudes beyond the wide latitudinal bands identified in that study.For instance, a boulder-dominated surface will always have boulder faces that point directly at the Sun, regardless of their latitude, during an asteroid rotation.Therefore, if the particle ejection events were driven by thermal fracturing, they could originate from any latitude on Bennu's surface.However, as indicated by the mean ΔT, particle ejection would occur more frequently on a per unit area basis from low latitudes because equatorial facets more frequently attain the maximum ΔT.Both of these findings-that ejection can occur from anywhere, and that it likely occurs preferentially at low latitudes-are consistent with the continued observations of particle ejection by OSIRIS-REx (Chesley et al. submitted; Pelgrift et al. in press) and the distribution of the events characterized in Lauretta and Hergenrother et al. (2019).

Particle ejection sites
To study the three particle ejection events described in Lauretta and Hergenrother et al.
(2019) more closely, we carved out regions surrounding their radiant points from the highresolution global shape model of Bennu.Figures 5, 6, and 7 show the modeled temperatures for the regions around the radiant points of the 6 January, 19 January, and 11 February 2019 particle ejection events, respectively.The 6 January 2019 event was located at a high southern latitude.
As this event was not as well captured by imaging data as the other two events, it has two possible radiant points, "near" and "far" (relative to the spacecraft; Lauretta and Hergenrother et al. 2019).The near radiant point is close to one of the potential pockets of sub-surface water ice that we identified earlier.Values of T AVG in this region were as low as 95 K, which, as mentioned previously, would enable sub-surface water ice to remain stable over geological timescales.Unlike the first event, the 19 January and 11 February 2019 events were located Confidential manuscript submitted to JGR: Planets much closer to the equator of Bennu, where no potential pockets of sub-surface water ice were identified.In these regions, T AVG was only as low as 202 K, far too warm to enable sub-surface water ice to remain stable for any amount of time.Therefore, if the particle ejection events are caused by a common mechanism then they cannot be caused by the sublimation of sub-surface water ice.However, this does not necessarily imply the absence of any sub-surface water ice near the poles of Bennu.
The regions surrounding the radiant points of the 19 January and 11 February 2019 particle events had relatively high ΔT, average of 122 K with a maximum of 142 K, due to their low latitudes on Bennu.For the 6 January 2019 particle event, both potential radiant points were within regions of relatively high ΔT, up to 140 K in this case, but the far radiant point was part of a much larger region of relatively high ΔT than the near radiant point.This consistent occurrence of relatively high ΔT implies that thermal fracturing could be the common mechanism of particle ejection for all three events, and that the far radiant point was the more likely source region for the 6 January 2019 event.

Unresolved small-scale surface roughness
From the imaging data and the analysis of the Approach-phase infrared observations (DellaGiustina and Emery et al. 2019), it is clear that Bennu exhibited small-scale topography that is not captured in the 75-cm-scale global shape model.In particular, an RMS slope of 43 ± 1° at diurnal thermal skin depth scales (~2 cm from equation 5) was derived from the Approachphase infrared observations in addition to the thermal inertia of Bennu.Therefore, it is possible that Bennu could host small-scale cold traps where surface water ice might be stable that were not resolved in the global shape model used here.To investigate the occurrence of centimeterscale cold traps, we considered several forms of unresolved small-scale roughness and added them to facets of the 12-m-scale global shape model (the shape model resolution used to analyze the Approach-phase infrared observations) within the ATPM.
First, we considered the ~77% fractional coverage of hemispherical craters (180° crater opening angle) used to fit the Approach-phase observations in DellaGiustina and Emery et al. (2019).However, the large amounts of self-heating that occurred within these hemispherical Confidential manuscript submitted to JGR: Planets craters prevented surface water ice from being stable, as the lowest T MAX obtained was 136 K.
We then considered a 100% coverage of craters with an opening angle of 120°, which gave a similar RMS slope but had less self-heating within the craters.In this case, surface water ice was stable near the poles of Bennu (Figure 8), as the lowest T MAX was 108 K.
where L is the horizontal baseline over which the surface roughness is measured, ν is the surface roughness RMS height deviation, and ν 0 is the RMS height deviation at the unit scale (Shepard et al. 2001).The surface roughness RMS height deviation over a specified baseline is given by ( ) where Δh is the change in topography over this baseline, and N is the total number of measurements of Δh.As we did not know whether the measured Hurst exponent from meterscale surface roughness extended to the centimeter-scale, we explored a range of Hurst exponents: 0.55, 0.75, 0.85 and 0.95.We also measured the RMS height deviation at the unit scale, i.e. ν 0 at L = 1 m, to be 0.1 m and 0.25 m from the global shape model and the infrared Confidential manuscript submitted to JGR: Planets observations, respectively.Both measurements of ν 0 and their averaged value of 0.175 m were used along with the four Hurst exponents to generate 12 different patches of artificial topography from a power spectral density function for fractal surfaces (Jacobs et al. 2017).In particular, we used the MATLAB code of Kanafi (2020) to generate these patches, and we subsequently added them to the south pole of Bennu (the colder of the two poles) within the ATPM to verify the spherical-section crater results.
Four examples of the artificial topography with RMS slopes of ~40° are shown in Figure 9.The coldest facets of these artificial topographies had values of T MAX below 131 K. Therefore, this finding with more realistic topography confirms that millimeter-thick layers of surface water ice could be stable at centimeter scales near the poles of Bennu.However, as the small-scale cold traps occur only at high latitudes, our conclusion is unchanged that thermal fracturing, as opposed to the sublimation of water ice, is a more likely temperature-driven mechanism for particle ejection on Bennu.Confidential manuscript submitted to JGR: Planets

Supply of material driving particle ejection
The continued observations of particle ejection by OSIRIS-REx (Pelgrift et al. submitted) raises questions as to whether this phenomenon is relatively new for Bennu or an ongoing process.If it is an ongoing process, then a steady source of the material driving particle ejection that has not become depleted over the dynamical age of Bennu is required.In terms of the more likely temperature-driven cause of rock thermal fracture, Bennu's many rocks and large boulders (Lauretta et al. 2019) would provide such a source.Furthermore, surface rocks that have been broken down to regolith by thermal fracturing may be replaced over time by rocks from Bennu's interior.Interior rocks deeper than a few meters are protected from the diurnal thermal cycling that leads to rock breakdown, and they will only start to breakdown once brought to the surface.
Granular matter processes, such as the Brazil Nut effect, can cause regolith to sink to the interior of a rubble-pile asteroid whilst bringing larger rocks up to its surface (Murdoch et al. 2015).
Therefore, if particle ejection is driven by rock thermal fracture then it may only stop once Bennu has been entirely converted to fine-particulate regolith.
In terms of maintaining the less likely temperature-driven cause of water ice sublimation, water ice present in Bennu's interior could be drawn towards the surface during colder points of the year.As the saturation vapor pressure is exponentially temperature-dependent, given enough free water molecules, an effective vapor density gradient will draw molecules towards colder areas.This driving force could re-supply near-surface cold traps from below.However, this process demands that a relatively large concentration, at least mono-layer coverage, of water is available and that the recharge occurs quickly (Schorghofer and Taylor 2007).This fast migration may require a high-porosity pathway to open, such as that believed to occur on comets, to move fast enough but still allow the initial ice reservoir to be seasonally stable.
To determine whether Bennu's interior could be a significant reservoir of water ice, we estimated the core temperature of Bennu by assuming that it has already reached its equilibrium temperature with its current orbit.We found an upper bound of 257 K for the core temperature of Bennu by spatially averaging the T AVG across its surface.This result is consistent with the upper bound of 267 K obtained from the analytical expression derived by Schorghofer and Hsieh (2018) for a spherical asteroid.With such a high core temperature, any interior water ice that was initially present within Bennu would have already migrated to its cold polar regions.This estimate breaks down if Bennu's core has had insufficient time to equilibrate after arriving at its current orbit.However, this is highly unlikely because the time it takes for heat to propagate from Bennu's surface to its core (~200,000 years given Bennu's radius and thermal inertia) is much less than the estimated time it takes for an asteroid to migrate from the main belt to a Bennu-like orbit (~1 to 10 million years; Bottke et al. 1996).
Additionally, the polar regions where water ice, if buried within the top few meters of the surface, could persist over geological timescales were found based on the assumption that Bennu's surface, pole orientation, and orbit are constant with time.In particular, the T AVG < 145 K criterion obtained by Schorghofer (2008) that we used in this work implies that sub-surface water ice could be stable for over a billion years if the surface illumination conditions are kept constant.Maintaining a fairly constant pole orientation and orbit is essentially impossible for a small near-Earth asteroid such as Bennu.The YORP effect is known to be accelerating the spin rate of Bennu (Nolan et al. 2019;Hergenrother et al. 2019), and the Yarkovsky effect is shrinking its orbit (Chesley et al. 2014;Scheeres et al. 2019).Therefore, the presence of any Confidential manuscript submitted to JGR: Planets water ice within Bennu depends strongly on its past dynamical history.However, previous work has shown that the obliquity component of the YORP effect rapidly moves the pole of Bennu to the stable orbit-perpendicular configuration it is in now (within ~10 5 years; Statler 2015), and that there is an ~85% probability that Bennu has not approached the Sun to within 0.4 AU (Delbó and Michel 2011).Re-surfacing events may also erase and create new polar cold traps over time, but from equatorial crater counts Bennu's surface is estimated to be rather old at between 100 million to 1 billion years in age (Walsh et al. 2019).If Bennu had water ice, these conditions may have ensured that it was retained during the asteroid's migration from the main belt.

Particle ejection mechanisms
In our temperature modeling of Bennu, we looked for relationships between modeled temperatures and radiant points, and we did not address precisely how particles are ejected from Bennu.For sublimation of water ice, particle ejection can be triggered when stable water ice is suddenly exposed to direct sunlight by impact excavation, or perhaps by a protecting rock that underwent thermal fracture.However, as noted, the particle ejection events should consistently originate from near the poles of Bennu if the sublimation of water ice is the driving mechanism.
For rock thermal fracture, simulations by rock thermomechanical models are required to explain how particles can be ejected via this mechanism.Such simulations are addressed in Molaro et al. (submitted) but here we can provide some evidence to help establish the processes involved.For instance, the late afternoon timing (between 15:22 and 18:05) of the three large perhaps the late afternoon timing is related to when sub-surface layers reach their diurnal maximum temperature.Figure 10 shows that layers located at one to two diurnal thermal skin depths below the surface of Bennu reach their maximum temperature at approximately this late afternoon local solar time.The surface is also cooling at this time, which makes it susceptible to tensile stresses induced by thermal contraction.As the diurnal thermal skin depth is ~2 cm for Bennu (from equation 5), the depths of these sub-surface layers are comparable to the size range of the ejected particles observed by OSIRIS-REx (<1 to 10 cm).Therefore, the rock thermomechanical failure that initiates particle ejection could be triggered at this depth.101955) Bennu.This is for the perihelion example given in Figure 2. From equation ( 5), the diurnal thermal skin depth is ~2 cm for Bennu.

Summary and conclusions
After modeling the global surface and sub-surface temperatures of Bennu, we find that ~1856 m 2 of Bennu's polar regions at latitudes above ~60° N and S have average temperatures that are below 145 K.These are sufficiently cold to enable water ice, if buried within the top few meters of the surface, to remain stable over geological timescales -up to a billion years if Bennu's surface, orbit, and pole orientation remain constant.Additionally, we find that millimeter-thick layers of surface water ice are stable over ~10 3 -year timescales within polar centimeter-scale cold traps with maximum surface temperatures below 131 K. Therefore, particle ejection would be limited to high latitudes if it were solely caused by water ice sublimation, as conditions enabling ice stability are only found near the poles.However, particle ejection has also been observed to occur from the much warmer equatorial region of Bennu (Lauretta and Hergenrother et al. 2019).We find relatively high amplitudes of diurnal temperature variation, a proxy for the efficacy of rock thermal fracture, at all latitudes on Bennu due to the extreme ruggedness of its shape (ΔT of ~140 K).Therefore, if rock thermal fracture is the driving mechanism behind particle ejection, it could occur from any latitude on Bennu's surface.
However, as the amount of surface area that reaches maximum ΔT is greatest near the equator, particle ejection is more likely at low latitudes.These findings are consistent with the continued observations of particle ejection by OSIRIS-REx (Chesley et al. submitted; Pelgrift et al. in press).
g. ElMir et al. 2019;Graves et al. 2019), and there are limited constraints on the rate at which this process may fracture or disaggregate rocks on planetary surfaces.A full treatment with rock thermomechanical models is beyond the scope of this work.Instead, we evaluate the spatial distribution of ΔT across the surface of Bennu, which has been shown to be the best and most simple proxy for first order estimation of the efficacy of rock breakdown on airless planetary bodies (e.g.Boelhouwers and Jonsson 2013;Molaro 2015;Molaro et al. 2015).Previous work byHamm et al. (2019) investigated the generation of regolith by thermal

Figure 2 .
Figure 2. Diurnal temperature variations for an equatorial surface element on (101955) Bennu at perihelion and aphelion (0.90 and 1.36 AU, respectively).In this example, maximum surface temperature (T MAX ) and change in temperature (ΔT) were 390 and 140 K, respectively, at perihelion, and 303 and 82 K at aphelion for a thermal inertia of 350 J m -2 K -1 s -1/2 .

Figure 3
Figure 3 shows the T MAX , T AVG , and ΔT that were modeled by the ATPM for the two hemispheres of Bennu.No location on the surface of Bennu had a T MAX of less than 131 K

Figure 3 .
Figure 3. Global distribution of (a,b) maximum surface temperature (T MAX ), (c,d) average temperature (T AVG ), and (e,f) maximum ΔT for (101955) Bennu.In the northern hemisphere (a,c,e), the circle and square indicate the locations of the radiant points for the 19 January and the 11 February 2019 particle ejection events, respectively.In the southern hemisphere (b,d,f), the circle and square indicate the locations of the near and far radiant points, respectively, for the 6 January 2019 particle ejection event.

Figure 4 .
Figure 4. Latitudinal distribution of ΔT and potential area of sub-surface water ice for (101955) Bennu.The solid and dashed red lines indicate the mean and maximum ΔT, respectively, for global shape model facets that have been binned into 2° latitudinal steps.The blue line indicates the total area of facets within these latitudinal bins that have average temperatures (T AVG ) below 145 K.

Figure 5 .Figure 6 .
Figure 5. Local distribution of (a) surface temperature at time of particle ejection (T), (b) maximum surface temperature (T MAX ), (c) average temperature (T AVG ), and (d) maximum ΔT for the two candidate sites (near and far radiant points) of the 6 January 2019 particle ejection event.The circle and square indicate the locations of the near (latitude = -75°, longitude = 325°) and far (latitude = -57°, longitude = 344°) radiant points, respectively.The global shape model segment shown here is from latitude -80° to -50° and longitude 310° to 360°.The positional uncertainties of the radiant points on these maps are approximately 20 m along the line that joins them.

Figure 7 .
Figure 7. Local distribution of (a) surface temperature at time of particle ejection (T), (b) maximum surface temperature (T MAX ), (c) average temperature (T AVG ), and (d) maximum ΔT for the 11 February 2019 particle ejection site.The circle indicates the location of the radiant point for this ejection event (latitude = 21°, longitude = 60°).The global shape model segment shown here is from latitude 10° to 30° and longitude 50° to 70°.The positional uncertainty of the radiant point on these maps is approximately 5 m in radius.

Figure 8 .
Figure 8. Distribution of minimum peak surface temperature within unresolved small-scale surface roughness for the (a) northern and (b) southern hemispheres of (101955) Bennu.In this example, the maximum temperatures of the coldest facets within spherical-section craters of 120° opening angle are shown.The circle and square indicate the same particle radiant points that were defined in Figure 3.

Figure 9 .
Figure 9. Example renderings of unresolved small-scale surface roughness with RMS slopes of ~40°.These were placed at the south pole of Bennu, the colder of its two poles, and the T MAX values were sought by running the ATPM around the orbit of Bennu.As demonstrated, hypothetical millimeter-thick layers of surface water ice were stable over ~10 3 -year timescales within these artificial topographies because the lowest T MAX values were <131 K.In particular, the lowest T MAX values obtained for H of (a) 0.55, (b) 0.75, (c) 0.85, and (d) 0.95 were 80, 94, 106, and 106 K, respectively.
particle ejection events characterized inLauretta and Hergenrother et al. (2019) is a curious feature.As shown in Figures5, 6, and 7, there are large spatial variations in surface temperature caused by shadows at the local solar times of the particle radiant points.If shadowing is an important part of the process, particle ejection should also occur shortly after local sunrise, when a similar amount of shadowing occurs.As suggested inLauretta and Hergenrother et al. (2019),

Figure 10 .
Figure 10.Diurnal temperature variations at different diurnal thermal skin depths below an equatorial surface element on (101955) Bennu.This is for the perihelion example given in Figure