Asteroid (101955) Bennu’s weak boulders and thermally anomalous equator

See allHide authors and affiliations

Science Advances  08 Oct 2020:
Vol. 6, no. 41, eabc3699
DOI: 10.1126/sciadv.abc3699


Thermal inertia and surface roughness are proxies for the physical characteristics of planetary surfaces. Global maps of these two properties distinguish the boulder population on near-Earth asteroid (NEA) (101955) Bennu into two types that differ in strength, and both have lower thermal inertia than expected for boulders and meteorites. Neither has strongly temperature-dependent thermal properties. The weaker boulder type probably would not survive atmospheric entry and thus may not be represented in the meteorite collection. The maps also show a high–thermal inertia band at Bennu’s equator, which might be explained by processes such as compaction or strength sorting during mass movement, but these explanations are not wholly consistent with other data. Our findings imply that other C-complex NEAs likely have boulders similar to those on Bennu rather than finer-particulate regoliths. A tentative correlation between albedo and thermal inertia of C-complex NEAs may be due to relative abundances of boulder types.


Temperature changes in response to time-varying illumination conditions are used to determine the physical properties of the surfaces of planetary bodies. For example, thermal emission measurements of the lunar surface (1) predicted the presence of a powder-like surface material years before the fine-particulate lunar regolith was confirmed by robotic landers and studied in situ by the Apollo missions. Global and local thermal emission mapping of the Moon and Mars has since facilitated regional and local characterization of regolith properties (24), boulder abundance (57), and a variety of geologic processes such as impact cratering, volcanic resurfacing, and rock breakdown (813).

The two thermophysical properties that dictate the temperature distribution and thermal emission of an airless planetary surface are the thermal inertia and surface roughness. Thermal inertia (Γ) is a measure of a material’s resistance to temperature change with time. It is defined by thermal conductivity (k), specific heat capacity (CP), and density (ρ) as Γ = (kρCP)1/2. Roughness is the irregularity of a surface measured over spatial scales that are comparable to or larger than the diurnal thermal skin depth but smaller than the spatial resolution of the known topography included in the thermal model. The diurnal thermal skin depth, dS = (kP/πρCP)1/2, where P is the rotation period, is the depth at which the amplitude of temperature variations decreases by 1/e (a few centimeters on Bennu). Surface roughness causes the thermal-infrared beaming effect, which is an enhancement of the surface thermal emission in the sunward direction (14, 15). Both properties can be quantified by fitting remote sensing observations of infrared thermal emission with a thermophysical model (16).

Spitzer Space Telescope measurements of the thermal emission of C-complex near-Earth asteroid (NEA) (101955) Bennu, the target of NASA’s OSIRIS-REx (Origins, Spectral Interpretation, Resource Identification, and Security–Regolith Explorer) asteroid sample return mission, indicated a low thermal inertia (310 ± 70 J m−2 K−1 s−1/2), which was interpreted as evidence that Bennu would be covered by millimeter- to centimeter-scale regolith (17). A spatially variable thermal inertia was expected that could be used to map boulder abundance and the particle size of fine regolith across the surface of the asteroid, as has been done on the Moon and Mars (4, 7, 18), to investigate the geologic history of the asteroid surface and aid in sample site selection and characterization (19, 20). However, when the OSIRIS-REx spacecraft arrived at Bennu in December 2018, images revealed that the surface is densely covered with boulders and that centimeter-scale or smaller regolith is apparently rare (21), in contrast to the prearrival predictions.

The OSIRIS-REx spacecraft is equipped with two instruments capable of measuring infrared thermal emission: the OSIRIS-REx Thermal Emission Spectrometer (OTES) (22) and the OSIRIS-REx Visible and InfraRed Spectrometer (OVIRS) (23). Disk-integrated data returned by OTES and OVIRS during the Approach phase of the mission (21, 24) confirmed the previous Spitzer-based thermal inertia measurement of Bennu: A global-average thermal inertia of 350 ± 20 J m−2 K−1 s−1/2 and a surface roughness RMS (root mean square) slope of 43° ± 1° (21) were derived from the OSIRIS-REx Approach-phase thermal emission light curves and the encounter-based shape model of Bennu (25). This thermal inertia value is much lower than the typical boulder and bedrock values of ~1500 to 2500 J m−2 K−1 s−1/2 (2, 7) and those measured for the presumed Bennu analogs, the CM and CI carbonaceous chondrite meteorites (26). These rare meteorites have thermal inertia values of ~1030 J m−2 K−1 s−1/2 (27, 28). The disconnect between the ground-based predictions and the encounter observations of Bennu led to the hypotheses that the boulders on Bennu must be either coated in thin layers of dust and/or highly porous (21).

The Hayabusa2 mission recently explored (162173) Ryugu, a similar boulder-covered C-complex NEA (29). Measurements by the onboard Thermal Infrared Imager (TIR) indicated a global-average thermal inertia of ~300 J m−2 K−1 s−1/2, with some small (meter-scale), apparently rare boulders with higher thermal inertias (600 to 1000 J m−2 K−1 s−1/2) identified in higher-resolution descent thermal images (29, 30). The MARA infrared radiometer on the MASCOT (Mobile Asteroid Surface Scout) lander measured a full diurnal temperature cycle of one boulder on the surface. Analysis of those data indicates a thermal inertia of ~282 (+95/−35) J m−2 K−1 s−1/2 at around 320 K (31), consistent with the global-average thermal inertia determined by TIR and indicating that the boulders on Bennu and Ryugu might share similar properties. The analysis furthermore showed that the MASCOT thermal data are not consistent with dust coatings, suggesting that the low thermal inertia values are characteristic of the rocks themselves (31).

The boulders on Bennu and Ryugu apparently have distinct thermophysical properties from those so far observed on the Moon and Mars (57). These findings have forced a reexamination of how thermal data and thermophysical model results should be interpreted in terms of the relative abundance and physical properties of boulders and regolith (defined here as loose particles smaller than the thermal skin depth). In addition, early images from the OSIRIS-REx encounter with Bennu showed a range of boulder reflectances (21, 24), a finding that persists in more detailed observations that indicate a bimodal reflectance distribution (32). This result raises the question of whether the differences in reflectance correspond to differences in thermophysical properties.

The thermophysical properties of spacecraft-explored asteroids are key to understanding their origins and interpreting telescopic observations of the broader asteroid population. Here, we present detailed global maps of Bennu’s thermal inertia and roughness and discuss the implications for boulder properties and histories, surface processes, other NEAs, and the sample ultimately returned by OSIRIS-REx.


Global thermal inertia and surface roughness maps

Remote observations of Bennu during OSIRIS-REx proximity operations were divided into several phases (20). Detailed Survey was the main global mapping phase, and it was split into two subphases: Baseball Diamond (March and April 2019), optimized for imaging, and Equatorial Stations (late April through early June 2019), optimized for spectroscopy. Targeted observations of potential sample sites at higher spatial resolution took place during the Reconnaissance (Recon) phase of the mission. In the Recon A subphase (October 2019), the spacecraft observed potential sample sites with a spatial resolution about five times higher than in Detailed Survey.

The primary datasets presented here were collected in Equatorial Stations, during which OTES and OVIRS measured the thermal emission from the surface of Bennu at ~40- and ~20-m spot sizes and wavelengths of 6 to 50 μm and 3.5 to 4.0 μm, respectively, at seven different local times of day: 3:20 a.m., 6:00 a.m., 10:00 a.m., 12:30 p.m., 3:00 p.m., 6:00 p.m., and 8:40 p.m. (fig. S1). The two instrument datasets complement each other well. The OTES wavelength range captures the bulk of the thermal spectral flux from the asteroid and is the primary thermal instrument. OVIRS measurements at shorter wavelengths are more sensitive to hotter facets and have a higher spatial resolution. We carried out thermophysical analysis of these data using a variant of the Advanced Thermophysical Model (ATPM) (3335) that has been adapted for use on Bennu (see Materials and Methods). On the basis of the apparent absence of fine-grained regolith on the surface (21) and the weak temperature dependence of thermal properties of analog meteorite materials (27, 28), we began with the assumption that Bennu’s thermal properties are constant with both temperature and depth. Analyses described below validate this assumption in the case of Bennu. We adopted a fractional coverage of hemispherical section craters to represent unresolved surface roughness. The OTES and OVIRS data were then analyzed separately to produce two sets of independently derived maps of thermal inertia and surface roughness for Bennu (Fig. 1).

Fig. 1 Thermal inertia and surface roughness maps of (101955) Bennu derived from OSIRIS-REx infrared observations.

(A and B) Thermal inertia maps from OTES (A) and OVIRS (B). (C and D) Surface roughness maps from OTES (C) and OVIRS (D). The color scales are offset by 20 U between the two thermal inertia maps but are equivalent for the two roughness maps. The black boxes indicate the locations of the regions with the lowest and highest thermal inertia that are shown in detail in Fig. 2 (C and D). Colorized values are overlaid on the global basemap of Bennu (99).

We find spatial variations in both thermal inertia and surface roughness on Bennu, and the OTES- and OVIRS-derived maps are in good spatial agreement (Fig. 1). Mean thermal inertia values of 300 ± 30 and 320 ± 30 J m−2 K−1 s−1/2 were derived from the OTES and OVIRS observations, respectively. For surface roughness, mean values from OTES and OVIRS were 40° ± 2° and 40° ± 3° RMS slope, respectively. The offset in mean thermal inertia between OTES and OVIRS appears to be systematic in nature, as the distributions of derived values are identical when the OVIRS-derived thermal inertia values are shifted by −20 J m−2 K−1 s−1/2 (fig. S2). Such a systematic offset suggests a subtle difference in instrument calibration and/or a small difference in spectral emissivity at near- and mid-infrared wavelengths. The difference of the results in this work from the previously reported global-average thermal inertia (350 ± 20 J m−2 K−1 s−1/2) (21) was caused by updating Bennu’s bolometric emissivity from 0.9 to 0.95 in this study, as the latter gave a better fit to the OTES spectra. The global-average measured thermal inertia and rotation period of Bennu (4.296 hours) correspond to a thermal skin depth of ~1 to 5 cm.

Our global maps are derived from both daytime and nighttime OSIRIS-REx measurements. Daytime-only measurements tend to result in a strong degeneracy between thermal inertia and surface roughness in model fit quality (i.e., equally good fits are obtained when thermal inertia is increased with increasing roughness; fig. S3), which has generally been difficult to break in previous thermophysical studies of asteroids. This daytime effect is caused by surface elements tilted toward the Sun being hotter than they would be on a smooth surface, mimicking a lower thermal inertia (33). However, the addition of nighttime data resolves this because the degeneracy has an opposite manifestation on the nightside (i.e., equally good fits are obtained when thermal inertia is decreased with increasing roughness; fig. S3). This nighttime effect is a result of surface facets facing each other. When a facet’s view to space is partially obscured by another part of the surface, it is not able to cool as efficiently, so a rough surface stays warmer at night than a smooth surface, mimicking a higher thermal inertia (33). We obtained mean uncertainties of 8 and 12 J m−2 K−1 s−1/2 for thermal inertia values derived from OTES and OVIRS, respectively. For surface roughness, these mean uncertainties were 0.8° and 1.4° RMS slope, respectively. Although these uncertainties may appear small, those for thermal inertia are consistent with the distribution of model result differences between OTES and OVIRS after accounting for their systematic offset in thermal inertia (fig. S4). The thermal inertia and roughness uncertainties were generally larger for OVIRS because it had less sensitivity than OTES for the measurement of nighttime thermal emission. Animations of Bennu’s daytime and nighttime brightness temperature variations are provided in movies S1 to S4.

Spatial variations in surface roughness

Spatial variations in surface roughness correlate with the spatial number density of boulders on the surface of Bennu (Fig. 1, C and D). This correlation is consistent with the interpretation that the infrared data are sensitive to all roughness spatial scales below the shape model resolution but larger than the diurnal thermal skin depth of a few centimeters (36). In particular, images of the areas of high surface roughness on Bennu reveal boulder piles that were not resolved in the 6-m-facet shape model used in the ATPM analysis. Such boulder piles also dominate the surface roughness measured on Eros and Itokawa by laser altimetry (37, 38). Our finding from Bennu thermal data is therefore consistent with independent methods of measuring surface roughness.

To determine the smallest spatial scale of roughness to which the infrared data are sensitive, we compared our surface roughness maps with the global 20-cm-facet shape model of Bennu derived from data acquired by the OSIRIS-REx Laser Altimeter (OLA) (3941). For a direct comparison, we first calculated the RMS slopes of the global 20-cm shape model with respect to the equivalent 6-m shape model. The results were then degraded to the spatial resolution of the OTES thermal roughness map for direct comparison. The global mean roughness value from the OLA map is 39° ± 7°, which is in excellent agreement with the OTES-derived global roughness (40° ± 2°). The two maps correlate spatially (figs. S5 and S6), and we therefore conclude that ~20 cm is indeed at or near the smallest spatial scale to which the infrared data are sensitive. This scale is 4 to 20 times as large as Bennu’s estimated diurnal thermal skin depth (~1 to 5 cm), so negligible lateral heat conduction occurs at the length scales at which roughness contributes to the thermal signature. We conclude that thermal surface roughness is a measure of the number and shape of rocks larger than ~20 cm resting on the surface of Bennu.

Spatial variations in thermal inertia

Variations in thermal inertia across Bennu’s surface appear to be primarily related to variations in the relative proportions of different boulder types. Reflectance and color data acquired by the MapCam imager in the OSIRIS-REx Camera Suite (42, 43) show that the majority of the surface is represented by a bimodal population of boulders (32). The average reflectance and color of Bennu can be reproduced by mixing the color spectra of these two boulder categories: “low reflectance” (~3.5 to 4.9% normal reflectance, MapCam 0.55-μm v-band) and “high reflectance” (~4.9 to 7.4% normal reflectance). Other, less common boulder types on the surface, such as very high-reflectance, pyroxene-bearing boulders (44), are typically small (less than a few meters), isolated, constitute a very small fraction of the global surface area, and would have a negligible effect on global thermal inertia studies; they are thus excluded from this analysis.

The areas with the lowest thermal inertia (~180 to 250 J m−2 K−1 s−1/2) on Bennu correspond most commonly with areas populated by the large, low-reflectance boulders. The areas with higher thermal inertia (~350 to 400 J m−2 K−1 s−1/2) tend to consist of mixtures of low- and high-reflectance boulders (Fig. 2). All high-reflectance boulders are smaller than an OTES spot in the Equatorial Stations observations, but several low-reflectance boulders are large enough to be resolved by OTES. We find that thermal inertia values correlate well with average surface reflectance values (Fig. 3). As noted in (32), the global variability in reflectance can be explained as a mixture of the two most dominant boulder types (low reflectance and high reflectance). Thus, we conclude that a substantial proportion of the global thermal inertia variation can be attributed to boulder-type distribution, where the low-reflectance boulders have lower mean thermal inertia values and the high-reflectance boulders have higher mean thermal inertia values. Furthermore, the average Bennu reflectance and color appear to be dominated by low-reflectance dark boulders or particles formed from them (32), indicating that low-reflectance boulders are more abundant than the high-reflectance boulders and thus contribute more strongly to the global average and local thermal inertia values.

Fig. 2 Regions with the lowest and highest thermal inertia on the surface of (101955) Bennu.

(A and B) Diurnal temperature curves. TI, thermal inertia. (C and D) Reference images. The OTES 9.24-μm brightness temperatures have been scaled to a common heliocentric distance and adjusted to a common spot geometry. The best model fits for these two end-member regions are shown for comparison. The images are photometrically corrected to units of normal reflectance (MapCam v-band global mosaic) (32).

Fig. 3 Two-dimensional histogram of thermal inertia (OTES) and MapCam v-band (550 nm) normal reflectance map from −70° to +70° latitude on (101955) Bennu.

Thermal inertia and reflectance from a total of 46,638 facets (~6-m facet size) are sorted into 50 bins.

The low-reflectance boulders tend to be the largest—i.e., all boulders with ≳20 m in diameter are of the low-reflectance type (32). This trend explains the observation that the low–thermal inertia areas on Bennu tend to correspond with areas covered by large boulders. The largest of these boulders, the Roc, is about 95 m in the longest dimension. Numerous OTES and OVIRS spots were acquired directly on the surface of this boulder and yield thermal inertia values of ~180 to 220 J m−2 K−1 s−1/2. Other large, low-reflectance boulders that were well resolved in the global thermal inertia map (Fig. 1, A and B) yield similar values. The high-reflectance boulders tend to be smaller (≲10 m in diameter) (21), which prevents the direct determination of their thermal inertia in our global-scale OTES and OVIRS datasets collected during Equatorial Stations. However, higher-resolution data were acquired during Recon A flybys of the OSIRIS-REx mission’s primary and backup sample sites, known as Nightingale and Osprey, respectively. These data provide an opportunity to test at the local scale whether the global spatial trends are linked to boulder type and to better estimate the thermal inertia of the high-reflectance boulders.

The Recon A flybys were conducted at an altitude of ~1 km, leading to an OTES spot diameter of ~10 m. The data were acquired during a narrow range of local times of day (~9:00 a.m. to ~2:00 p.m.). As a result, it is not possible to reliably separate thermal inertia and surface roughness in the ATPM model fits of the Recon A data owing to the strong degeneracy between these two parameters. To circumvent this degeneracy, we measured local variations in surface roughness from the OLA-derived digital terrain models of the candidate sample sites and converted them to OTES-equivalent roughness variations using the relationship derived from the Equatorial Stations global analysis. Therefore, in the ATPM fits to the Recon A data, surface roughness was kept fixed at the OLA-based values, and only thermal inertia was varied to find the best model fit. Figure 4 shows the resulting local thermal inertia maps for the primary (Nightingale) and back-up (Osprey) sample sites.

Fig. 4 Local thermal inertia maps of the Nightingale and Osprey candidate sample sites.

These maps were produced from Recon A OTES data using surface roughness measured by OLA. Values are overlaid on the global basemap (99).

Local thermal inertia variations in the Recon A maps generally reflect the trends seen in the global thermal inertia map. In particular, high-reflectance boulders have higher thermal inertia values than low-reflectance boulders at both the global and local scales. The uncertainties in the absolute thermal inertia values obtained from these local maps are higher, given the limited time-of-day coverage. Nevertheless, we find that the high-reflectance boulders have thermal inertia values in the range of ~400 to 700 J m−2 K−1 s−1/2, which, although higher than that of the low-reflectance boulders, is still lower than measured values for carbonaceous chondrite meteorites (~1030 J m−2 K−1 s−1/2) (27, 28). The modal thermal inertia values of these local maps are also in good agreement with their corresponding global map values (fig. S7). This agreement indicates that the global thermal inertia variations and trends indeed represent the relative local abundance of the different rock types and that Bennu’s surface is generally well mixed.

Data acquired at the Nightingale sample site afford the opportunity to analyze the thermophysical properties of what appears to be a fine-particulate regolith (unresolved material at pixel scales of ~1 to 2 cm per pixel) with relatively few meter-scale or larger boulders. The Nightingale crater has a lower thermal inertia compared to its surroundings, indicating the possible presence of fine particles smaller than the skin depth. However, the thermal inertia values at Nightingale are ~150 to 200 J m−2 K−1 s−1/2, comparable to the values for the low-reflectance boulders. Although advances have been made in models that used to predict regolith particle size from thermal inertia, a precise determination is still complicated by assumptions of the material properties of possible regolith particles in this crater (45, 46). If the particles are composed of material derived from the low-reflectance boulders, then this thermal inertia value could indicate that the particle diameters are predominantly greater than the diurnal skin depth of a few centimeters. Conversely, if the particles are composed of a material with higher density and higher thermal conductivity, then this thermal inertia value would correspond to particle diameters of ~0.7 to 1.8 cm. Very-high-resolution images (~0.4 cm per pixel) obtained during the Recon C phase from an altitude of ~250 m show that Nightingale crater hosts an abundance of centimeter-scale regolith uncharacteristic of most of the rest of Bennu.

Thermal inertia interpretations: Dust cover or structural properties?

We evaluated the inferences in (21) that the low overall thermal inertia of Bennu’s boulder-dominated surface could be due to (i) high porosity in the boulders from pores, voids, and fractures and/or (ii) dust coatings that would lower the apparent thermal inertia. Studies of Ryugu (31, 47) eliminated the hypothesis of dust coatings at the MASCOT site by showing that, in a layered thermal model, dust coatings on boulders would cause the diurnal temperature curves to have distinct shapes that were not observed in the data. We performed a similar study in which we fit a dust-layer model to OTES data of the regions identified in Fig. 2 with lowest and highest thermal inertia. We found that thin (>50 μm) laterally continuous layers of dust fail to reproduce the measured diurnal curves, particularly the afternoon and evening cooling portions (fig. S8). The model fit quality is acceptable for very thin layers of dust (<50 μm), so we cannot rule out such a very thin dust cover. Nevertheless, the best model fits are obtained for models with no dust present. Furthermore, a 50-μm dust layer has a minimal effect on the apparent thermal inertia: It is lowered by only ~10 J m−2 K−1 s−1/2, so even if a very thin dust layer is present, low–thermal inertia rock properties would still need to be invoked.

In addition to laterally continuous layers of dust, we considered an alternative, checkerboard mixing scenario where discrete pockets of thick dust layers are present next to bare exposures of rock. To determine whether a fractional mixture of dust with very low thermal inertia could explain the properties of Bennu’s boulders, we fitted a two-component model to OVIRS data of the regions of lowest and highest thermal inertia identified earlier (fig. S9). OVIRS data were used because the shorter wavelengths are more sensitive to the hotter daytime temperatures that dust would cause. As expected, the model compensated for higher dust fractions by increasing the rock thermal inertia values. However, the quality of model fits to the OVIRS data worsened progressively as the dust fraction increased. For example, very large quantities of dust (>20 and >80% surface coverage for the regions identified in Fig. 2 with the highest and lowest thermal inertia, respectively) led to an increase in the required rock substrate thermal inertia up to ~600 J m−2 K−1 s−1/2 but increased the RMS error of the fit by a factor of ~3 to 4. We cannot completely rule out a small areal coverage of dust (5 to 10%) due to the broad minima near the no-dust case in the model fit quality (fig. S9). Nevertheless, dust-free models achieve the best fits to the data. As with thin dust layers, even if small patches of dust were present, the underlying rock thermal inertias would be very similar to the dust-free values described above—i.e., the rock thermal inertia need only be raised by about 10 J m−2 K−1 s−1/2 to compensate for the presence of the upper limit of 10% dust.

Given that the diurnal thermal data are inconsistent with the presence of more than minor quantities of dust, we conclude that the thermal inertia of the boulders on Bennu represents the material properties of the boulders themselves. A negligible or minor contribution of dust is consistent with high-resolution images of boulder surfaces (~1 cm per pixel) obtained during the Recon A phase of the mission, which do not show evidence of dust ponding in local depressions. Spectral and color variations on Bennu (32, 48) furthermore indicate that the compositional diversity among the low- and high-reflectance boulders is not sufficient to drive the observed differences in thermal inertia. Thus, the thermal inertia of the boulders on Bennu is likely a consequence of their internal and near-surface physical structure.

Laboratory measurements of meteorite samples have shown that increased porosity leads to decreased thermal conductivity and thermal inertia; microcracks within measured meteorites act as barriers to the heat transfer (49, 50). The presence of void spaces, in general, also decreases thermal conductivity by decreasing the cross-sectional area and increasing the effective length of the pathways for heat conduction (51). The meteorite with the lowest thermal conductivity so far reported, the CM2 Cold Bokkeveld (27, 28, 49), is a mineralogical analog for Bennu and has a porosity of ~15% (52). However, its thermal inertia of ~1030 J m−2 K−1 s−1/2 at ~260 K is far above our values for low-reflectance boulders and moderately higher than our estimates for the high-reflectance boulders. The overall lower thermal inertias of boulders on Bennu therefore suggest porosities much higher than 15%. Models for thermal conductivity as a function of porosity (31, 50, 51) predict that the low-reflectance boulder thermal inertia of ~180 to 250 J m−2 K−1 s−1/2 would correspond to a porosity value of ~49 to 55%, whereas our high-reflectance boulder thermal inertia value of ~400 to 700 J m−2 K−1 s−1/2 would correspond to a porosity of ~24 to 38%, at least within the upper few diurnal skin depths of the exposed surface (fig. S10).

The strength of these boulders can be estimated from the derived values of thermal conductivity and porosity. According to the model in (31), the tensile strength scales linearly with the boulder thermal conductivity and gives tensile strengths of ~0.10 to 0.15 MPa and ~0.31 to 0.78 MPa for the low- and high-reflectance boulders, respectively (fig. S10). Therefore, the tensile strength of the high-reflectance boulders is a factor of ~2 to 8 higher than that of the low-reflectance boulders. The inference of different strengths is qualitatively supported by the appearance of the boulders: The high-reflectance boulders generally have smoother surface textures, are more angular, and sometimes have pronounced linear fracture planes, giving the impression that they are denser and more competent. The low-reflectance boulders have hummocky surface textures and, in some cases, appear to be breccias, giving the impression that they might be more friable and less dense (32). Many of the high-reflectance boulders contain bright putative veins proposed to consist of carbonate minerals (48). These features suggest that the high-reflectance boulders have decreased porosity due to infill from hydrothermal mineralization.

Elevated thermal inertia at the equator

The global thermal inertia maps exhibit a distinct latitudinal trend that is not fully explained by variations in local relative abundance of low- and high-reflectance boulders, as inferred from surface reflectance data. The latitudinally binned thermal inertia is 40 to 60 J m−2 K−1 s−1/2 higher at the equator than at latitudes of >60o and decreases gradually with increasing latitude (Fig. 5). This trend is apparent in the 3:20 a.m. brightness temperature data (fig. S1). Furthermore, the southern mid-latitudes have a thermal inertia that is ~20 to 25 J m−2 K−1 s−1/2 lower than that in the northern mid-latitudes. We evaluated whether this latitudinal thermal inertia trend is due to a temperature-dependent thermal inertia by comparing results from the first flyby of the Baseball Diamond subphase, where the heliocentric distance was ~1.0 astronomical unit (AU) versus ~1.2 AU during the Equatorial Stations subphase. Despite an estimated difference in the diurnal mean temperature of ~25 K at these distances, which is comparable to the difference in diurnal mean temperature between equatorial and mid- to high latitudes, we find that the thermal inertia values from the two subphases are approximately equivalent (fig. S11). Any offset between the two derived sets of values is well below the thermal inertia uncertainty of ~10 J m−2 K−1 s−1/2. Thus, we find that the thermal inertia change cannot be more than 10 J m−2 K−1 s−1/2 for a 25-K temperature change. Particulate regoliths on airless bodies are expected to have strongly temperature-dependent thermal inertias (4), whereas meteorites have been found to be less temperature-dependent (27, 28, 49). The minimal temperature dependence of the thermal inertia is consistent with the observation that Bennu’s surface is dominated by boulders rather than regolith. We therefore conclude that the latitudinal trend in thermal inertia is likely to reflect real differences in the properties of the surface materials and that at most a small part of the trend (≲10 J m−2 K−1 s−1/2) may be attributed to temperature dependence.

Fig. 5 Variation of mean thermal inertia with latitude on (101955) Bennu.

This plot was produced by averaging the global map values of Fig. 1 within latitude bins of 5° in width. The OVIRS-derived values are offset from the OTES-derived values by +50 J m−2 K−1 s−1/2, and the error bars represent the variation of the mean thermal inertia uncertainty with latitude.

Several dynamical quantities of Bennu’s surface—e.g., gravitational slope, magnitude, and potential—vary with latitude (25, 53, 54). These variations arise because of Bennu’s “spinning-top” shape, low bulk density, and relatively fast spin rate. Therefore, it is plausible that the relationship between thermal inertia and latitude could be a result of mass movement and/or wasting induced by the dynamical environment of Bennu’s surface. Dozens of mass movement features have been identified on Bennu, with the majority pointing toward the equator (55). Thus, it is possible that mass movement toward the equator has led to a sorting of materials, leading to a greater abundance of higher–thermal inertia boulders at the equator. The relationship between the equatorial bulge and the impact craters on Bennu’s surface suggests that the bulge itself is a comparatively old feature (56). However, global mass movement might cause a continual resurfacing of the equatorial bulge but yet deposit a thin enough layer to preserve the underlying crater topography (55). In such a case, the materials deposited on the equator would still have longer exposure ages compared to the mid- and high latitudes, where mass movement has carried the older surface layers equatorward to expose fresher, underlying material.

Boulders with longer exposure ages are likely to have undergone more cumulative weathering due to processes such as meteoroid bombardment (57, 58) and thermal fatigue from temperature cycling (59, 60). The inferred strength difference between the low- and high-reflectance boulders suggests that these two boulder types may respond differently to these two weathering mechanisms, which could lead to changes in their relative size frequency distribution over time. Thus, the increased thermal inertia at the equator may be indicative of a boulder maturation process that favors the long-term survival of boulders with higher thermal inertia. However, the MapCam reflectance data do not show a clear indication that the equator has a significantly elevated proportion of high-reflectance boulders. Alternatively, compaction due to meteoroid bombardment may cause both boulder types to increase in density over time, which would lead to an increase in thermal inertia without necessarily changing the reflectance values.


Possible formation and evolution of the two boulder types

Spectral analysis of the OVIRS and OTES data reveals similar phyllosilicate mineralogy across the surface of Bennu, indicating that the two boulder types are also similar in this respect (48, 61). Nevertheless, MapCam images show color and reflectance diversity among the boulders (32). In addition, some high-reflectance boulders have been found to contain centimeter-thick veins, which may be composed of carbonates (48). The difference in the morphology between the low- and high-reflectance boulders also indicates differences in their structural properties. We now consider the stages in the formation and evolution of these boulders where porosity would have been introduced or reduced, including how and when the histories of these two boulder types may have diverged.

Bennu’s spectra most closely resemble those of moderately to heavily altered CI and CM chondrites (26, 32, 48). The minerals on Bennu’s surface have therefore passed through the earlier stages of alteration consisting of the transition from cronstedtite (Fe-rich phyllosilicate) to serpentine (Mg-rich). The extent of alteration of cronstedtite to serpentine is used as an indicator of the extent of alteration in CM chondrites (62). The alteration of cronstedtite to serpentine in Bennu’s parent body could introduce microporosity in the mineral fabric; the insertion of Mg into the cronstedtite structure causes a warping of the relatively planar cronstedtite sheets into more curved serpentine sheets (63). However, any porosity induced by this stage of alteration cannot explain the differences between the two boulder types, given that spectra show that they have both reached this stage.

In the model in (64), the aqueous alteration process occurs before lithification in a convecting mudball system. It is possible that the physical differences in Bennu’s boulders arose during lithification into coherent rock. The process of lithification likely involved the cementation of grains by minerals precipitated from the pore fluid. The local temperature, pressure, and availability and mobility of water would affect this process. Thus, vertical (i.e., with depth) or lateral variability in the hydrothermal systems on the parent body could have led to varying degrees of lithification in Bennu’s boulders. This variability in lithification may be responsible for the veins and their proposed pore-filling carbonate minerals observed exclusively in the high-reflectance boulders (48). This infilling could reduce porosity and increase the coherence of the rock, which could, in turn, increase the thermal inertia.

Impact brecciation on the parent body may have introduced considerable porosity. Nearly all carbonaceous chondrite meteorites are interpreted to be breccias (65). Materials closer to the surface of the parent body would be more substantially brecciated, whereas materials at depth would experience little or no brecciation. If the boulders have high initial porosity, then it is difficult to determine whether brecciation would increase or decrease the porosity and, therefore, the thermal inertia. As noted above, our inferred porosities of boulders on Bennu’s surface are substantially higher than the porosities of carbonaceous chondrite meteorite breccias.

Another possibility for the difference in porosity of the two boulder types is a process of porosity reduction. Modeling and experimental studies suggest that chondritic parent bodies accreted with an initially very high porosity, on the order of ~60 to 70% (6668). Porosity reduction would then occur through aqueous alteration, burial compaction, and/or impact. Porosity reduction from aqueous alteration does not explain the thermal inertia differences among Bennu’s boulders, as relative degrees of alteration would be readily detected spectrally (6971). Burial compaction would require a much larger body—a ~200-km-diameter asteroid is required to generate ~1 MPa of lithostatic pressure at its center, equivalent to ~10 m in depth on Earth (72)—but this could have occurred before the disruption and reaccumulation of Bennu (25, 73, 74). Considerable porosity within chondrites can also be efficiently removed via impacts with minimal dehydration or compositional changes (7577). Therefore, impact compaction of an initially highly porous chondritic material (the low-reflectance, low–thermal inertia boulders) would result in compositionally similar, but less porous material (high-reflectance, high–thermal inertia boulders). The highest porosity measured in a carbonaceous chondrite meteorite to date is 40% for the ungrouped C chondrite Tagish Lake (78), which also shares some spectral similarities to Bennu (61, 79).

The respective boulder populations may evolve differently over time in response to weathering mechanisms such as meteoroid bombardment and thermal fracturing. For example, meteoroid bombardment may more thoroughly fracture the comparatively denser, high-reflectance boulders, leading to large-scale fracturing and boulder subdivision, whereas the low-reflectance boulders would not dissipate energy as efficiently and might instead experience fragmentation that is more local to impact points. The two boulder types may experience different thermally induced stresses fields due to potential differences in their effective elastic moduli and may respond differently to said stresses. Although the lower-strength boulders might initially be expected to be more susceptible to thermal fatigue, they may be more elastic than the high-reflectance boulders, which would reduce the magnitude of thermally induced stresses. Conversely, they would, on average, experience larger diurnal temperature amplitudes due to their lower thermal inertia, which would act to increase thermal stresses. It is possible that the larger temperature amplitudes could be sufficient to compensate for stress reduction due to higher elasticity.

Nevertheless, we find it unlikely that thermal fracturing alone would be sufficient to produce the differences observed between the two boulder types. If all the boulders started out the same, then we would expect thermal fracturing to lead to a continuum of boulder types, controlled by degree of fracturing, but we do not observe such a continuum. In addition, thermal fracturing is only operative in the near surface of the boulders (within the upper few skins depths). Thus, a fractured rind would develop on the surface of the boulders, which then must be removed if the boulders are to be reduced in size. Such a fractured rind would act to further reduce the thermally induced stresses, leading to a self-limiting process. This may explain the low thermal inertia of the low-reflectance boulders and would in fact suggest that they have a longer lifetime.

Conclusions and implications

Analysis of thermal observations of Bennu by the OSIRIS-REx mission reveals that the surface is dominated by boulders with lower thermal inertias than their meteorite analogs. Furthermore, they fall into two distinct populations: (i) boulders with lower thermal inertias (i.e., ~180 to 250 J m−2 K−1 s−1/2) and lower inferred densities and strengths and (ii) boulders with higher thermal inertias (i.e., ~400 to 700 J m−2 K−1 s−1/2) and higher inferred densities and strengths, relative to one another. Photometric and textural analyses (32) also identify a bimodality that corresponds to the one we describe, such that the lower–thermal inertia (weaker) boulders have low reflectances and a hummocky texture, whereas the higher–thermal inertia (stronger) boulders have higher reflectances, are more angular, and contain evidence for hydrothermal mineral deposition in preexisting cracks and void spaces (48). The first-order thermal inertia variations on the surface of Bennu are explained by mixtures of these two boulder types. The spatial distribution of these boulders, as inferred from the thermal inertia maps, may point to a dynamic history of preferential preservation based on different weathering rates of the two boulder types.

The inferred thermal conductivity of the boulders on Bennu falls below any measured meteorite thermal conductivity value. However, the high-reflectance boulders have thermal inertias that approach the values for some CM chondrites (27, 28). Thus, CM and similar meteorites may be similar to the high-reflectance boulders in their phyllosilicate mineralogy and porosity. The low-reflectance boulders, on the other hand, have thermal conductivities far below those of our meteorite collections (28) and likely are distinct from any recognized meteorite types in terms of their thermomechanical properties. Recent laboratory experiments (80) have, however, shown that is possible to create simulants that have the composition of CMs but thermal conductivities intermediate of those of known meteorites and those of the boulders of Bennu and Ryugu. Regardless of the mechanism that caused the high porosity, the inferred weakness of this material is consistent with the implication that materials of this type would rarely or never survive atmospheric entry. Given the abundance of the low-reflectance boulders on Bennu (32), it is likely that some of this material will be present in the returned sample. We therefore expect OSIRIS-REx to return for analysis material that is not currently in Earth’s meteoritic collections; see also (48).

Peak stresses in meter-scale boulders on different planetary bodies are expected to increase with decreasing perihelion distance [(81); see also (59, 82, 83)]. Thus, if thermal fracturing more rapidly reduces the surface coverage of lower–thermal inertia, low-reflectance boulders, then this could explain the observed absence of low-albedo asteroids at low perihelia (59, 84) and, reciprocally, the comparatively high thermal inertia (~600 J m−2 K−1 s−1/2) and high albedo of (3200) Phaethon (85), which has a perihelion distance of only 0.14 AU. Whatever mechanism dominates, preferential degradation of weaker, lower-reflectance boulders would lead to an increase in albedo and thermal inertia of C-complex NEAs.

Thermal inertia has been determined for seven C-complex NEAs, including Bennu and Ryugu. These seven NEAs show an apparent correlation between albedo and thermal inertia (fig. S12) that is reminiscent of the trend with boulder type observed on Bennu. Three of these NEAs (Betulia, 1996 FG3, and 2002 CE26) (86, 87) have global thermal inertias that are comparable to those of Bennu’s low-reflectance boulders, and the other two (Phaethon and 2008 EV5) (85, 88) have thermal inertias similar to those of Bennu’s high-reflectance boulders. Given that the thermal inertias are all in the range of those of Bennu’s boulders, we infer that the surfaces of these C-complex NEAs are also dominated by boulders rather than by finer-particulate regolith. The apparent correlation suggests that mixing of boulder types of different thermal properties may be a common feature of small C-complex asteroids.


OSIRIS-REx global infrared observations

OTES is a Fourier transform infrared spectrometer that operates over the 5.71- to 100-μm (1750 to 100 cm−1) spectral range (22). OVIRS uses wedge (linearly varying) filters mounted just above the detector to measure spectra over the 0.4- to 4.0-μm spectral range (23). In daytime data, thermal radiance from Bennu is greater than reflected radiance longward of about 3.0 μm (depending on the illumination and surface temperature). Both of these instruments are point spectrometers with fields of view of 8 and 4 mrad for OTES and OVIRS, respectively. During the Baseball Diamond Flyby 1 (BBD 1) and Equatorial Stations, the spacecraft was at a range of ~5 km above the equator, leading to spot sizes of ~40 and ~20 m on the surface of Bennu for OTES and OVIRS, respectively.

At each of the Equatorial Stations, the spacecraft was positioned above the equator at the targeted local solar time (3:00 p.m., 3:20 a.m., 12:30 p.m., 10:00 a.m., 6:00 a.m., 8:40 p.m., and 6:00 p.m. for stations 1, 2, 3, 4, 5, 6, and 7, respectively). Spacecraft pointing was then slewed north-south (N-S) across the disk of Bennu, extending off to space for calibration at the ends of each slew, with slew rates set to provide some overlap between OTES and OVIRS measurement spots. These N-S slews were continued without interruption for a bit more than a full rotation of Bennu to ensure complete surface coverage at each station. For the 12:30 p.m. station, a few small gaps in OVIRS coverage occurred near the equator, leading to the uncolored facets in Fig. 1. The observations were therefore nadir at the equator and increased in emission angle with latitude. The Equatorial Stations were executed in the order given above, approximately a week apart from each other: 25 April 2019; 2, 9, 16, 23, and 30 May 2019; and 6 June 2019. The spectral observation sequences started between 17:19 and 18:02 UT on these days. Integration times were ~2 s for OTES and ~1 s for OVIRS, leading to ~3000 OTES spot spectra and ~5000 to ~7000 OVIRS spot spectra per station. The observing sequence for the BBD 1 station was slightly different. That station was optimized for imaging and included N-S slews pointed 45° to the east and west along with the nadir-relative slew. These slews did not extend fully to the poles and include large gaps in the OTES and OVIRS data between slews. BBD 1 was executed on 7 March 2019 and spectral observations started at ~17:00 UT.

We began our analysis with data that were run through the instrument pipeline to convert from instrument signal to physical spectral radiance units, performing all necessary corrections and resampling. The spectra we fit with the thermal model (see subsequent sections) are therefore measurements of total radiance from the surface over the wavelength ranges listed above. We restricted thermal analysis of OVIRS data to the 3.5- to 4.0-μm spectra range, where thermal radiance exceeds reflected by factors of ~5 to ≥20. For daytime observations, it was necessary to remove the contributions from reflected solar radiance at these wavelengths. We fit a solar spectrum, slightly adjusted for the mean near-infrared spectral slope of Bennu, to each OVIRS spectrum at ~2.0 to 2.1 μm. This scaled solar spectrum was subtracted from the total measured radiance to isolate the thermal component in each OVIRS spectrum.

OTES operates by measuring the difference in radiance between the observed spot and the detector. The calibration is not robust when the average surface temperature within the spot is within a few degrees of the detector temperature, which creates a phase inversion issue. These spots are flagged and omitted from our analysis. The measured spots cover surface facets that have a distribution of temperatures, and the spot radiances are therefore not generally fit well by a single temperature blackbody. As described in subsequent sections, we therefore model the measured radiances rather than derived temperatures. Nevertheless, fig. S1 shows maps of brightness temperature at 9.24 μm for the 12:30 p.m. and 3:20 a.m. stations to illustrate day- and nighttime measured temperature distributions. The temperatures are mapped to the stereophotoclinometric (SPC) v20 shape model (25) degraded to ~200,000 triangular facets with a resolution of ~3 m. Each 40-m OTES spot therefore includes many facets, and because of the overlap of the spots, each facet is included in multiple OTES spots. To create the map, we identify all spots that cover a given facet on the shape model and average the temperatures of those spots, weighted by T4.

Thermophysical model

To model and fit the OSIRIS-REx infrared observations of Bennu, we use the ATPM developed in (3335). For instance, the ATPM has previously been applied to Bennu in (21, 53, 87, 89, 90). To briefly describe how it works, the ATPM solves the one-dimensional (1D) heat conduction equation with a suitable surface boundary condition for each triangular facet of a user-specified input global shape or local digital terrain model. In particular, 1D heat conduction is given bydTdt=1ρCPddz(kdTdz)(1)where T is temperature, t is time, and z is depth. In addition, k, ρ, and CP are the thermal conductivity, density, and specific heat capacity of the asteroid surface material, respectively, which are combinable into the single thermal inertia parameter, Γ, viaΓ=kρCP(2)

To ensure conservation of energy between incoming and outgoing radiation, the surface boundary condition is given by(1AB)([1S(t)]ψ(t)FSUNrH2(t)+FSCAT(t))+FRAD(t)+k(dTdz)z=0εBσTz=04=0(3)where AB is the Bond albedo, FSUN is the integrated solar flux at 1 AU (i.e., 1367 W m−2), rH(t) is the heliocentric distance of the asteroid in AU at time t, εB is the bolometric emissivity, and σ is the Stefan-Boltzmann constant. S(t), ψ(t), FSCAT(t), and FRAD(t) are functions that evaluate shadowing, the cosine of the Sun illumination angle, the total amount of multiple-scattered sunlight, and the total amount of self-heating for each facet, respectively. Detailed descriptions of how these functions are evaluated in the ATPM are provided in (3335).

Equations 1 and 3 are solved in the ATPM with a finite difference method and a series of Newton-Raphson iterations, respectively, after normalizing them to the length scalel=kP2πρCP(4)where P is the asteroid rotation period (91). In particular, subsurface temperatures are computed over 56 depth steps to a maximum depth of eight diurnal thermal skin depths and over 650 time steps per asteroid rotation (92). After initializing the ATPM with rotationally averaged instantaneous equilibrium temperatures, the model is then run for several tens of asteroid rotations until the surface temperatures converge to a tolerance level of 10−3 K.

Unresolved surface roughness is also accounted for within the ATPM by adding a fractional coverage, fR, of spherical section craters to each facet of the input shape/topography model. These spherical section craters are of the form specified in (36) and include full heat conduction for each crater subfacet. For our analysis in this work, we adopt hemispherical craters (i.e., with crater opening angle of 180°) that are constructed from 100 equally spaced subfacets, upon which Eqs. 1 and 3 are also applied. Conveniently, a fractional coverage of hemispherical craters can be specified in terms of surface roughness RMS slope, θ, viaθ=49fR(5)

For various shape/topography models, the ATPM was run for the midpoint positional geometry of Bennu for each of the observational stations for which we analyzed infrared data. In particular, table S1 summarizes the positional geometry of Bennu that we used in the ATPM temperature modeling for the seven Equatorial Stations, BBD 1, and the two Recon A flybys analyzed in this work. In addition, we adopted within the ATPM Bennu’s measured pole orientation, λH = 69.92° and βH = −83.45°, and rotation period, 4.296 hours. (24). In terms of thermophysical properties, surface temperatures were generated for a Bond albedo of 0.02, bolometric emissivity values of 0.90, 0.95, and 1.00, and thermal inertia ranging from 0 to 800 J m−2 K−1 s−1/2 in steps of 10 J m−2 K−1 s−1/2. The modeled surface temperatures were then stored in lookup tables to enable the production of synthetic radiance spectra for each OTES and OVIRS spot.

To generate synthetic radiance spectra for a given OTES or OVIRS spot, the model flux was determined for each facet and crater subfacet that fell within the spot’s field of view. In particular, the model observed flux, FMOD,i(Γ,λ,t), from facet i at wavelength λ is given byFMOD,i(Γ,λ,t)=B(λ,Ti(Γ,t))ε(λ)aidi2(t)cosφi(t)(6)where B(λ,Ti(Γ,t)) is the Planck function, ε(λ) is the spectral emissivity, and ai is the area of the facet. In addition, Ti(Γ,t), di(t), and φi(t) are the temperature, distance to OSIRIS-REx, and observation angle of the facet, respectively, at time t. In particular, the Ti(Γ,t) were retrieved from the previously built lookup tables, and the di(t) and φi(t) were calculated from the positional and pointing geometry of OSIRIS-REx at the time of the observations. The total model observed flux for a given spot n, FMOD,n(Γ,fR,λ,t), was then determined by summing across all facets and crater subfacets that fell within its field of view, as given byFMOD,n(Γ,fR,λ,t)=Σi=1NFACETvi(t)((1fR)FMOD,i(Γ,λ,t)+fRΣj=1100vj(t)FMOD,j(Γ,λ,t))(7)where vi(t) and vj(t) indicate whether facet i and crater subfacet j were visible to OSIRIS-REx at the time of the observations, respectively and fR is the previously defined fractional coverage of hemispherical craters. Last, the total model observed flux was converted to radiance, RMOD,n(Γ,fR,λ,t), by dividing by the total solid angle of the observed facets, Ωn(t)RMOD,n(Γ,fR,λ,t)=FMOD,n(Γ,fR,λ,t)Ωn(t)(8)where Ωn(t) was calculated byΩn(t)=Σi=1NFACETvi(t)aicosφi(t)di2(t)(9)

Equations 6 to 9 were only applied to facets and crater subfacets that fell within an angular radius of 4 or 2 mrad from the OSIRIS-REx pointing vector depending on whether an OTES or OVIRS spot, respectively, was being simulated (22, 23).

Fitting a single OTES or OVIRS spot

The OSIRIS-REx infrared observations were fitted in terms of brightness temperature because the OTES measurements were primarily reported in these units. As in previous work that fitted similar datasets for other planetary bodies, e.g., (7, 93), we calculated the RMS error, RMSE(Γ,fR), between the measured and modeled brightness temperatures, TOBS,n(λ) and TMOD,n(sfn,Γ,fR,λ), respectively. For fitting a single OTES or OVIRS spot, RMSE(Γ,fR) was calculated usingRMSE(Γ,fR)=1NλΣNλ(TOBS,n(λ)TMOD,n(sfn,Γ,fR,λ))2(10)where Nλ is the number of wavelengths that were observed by OTES or OVIRS. In addition, the TMOD,n(sfn,Γ,fR,λ) were calculated from the modeled radiance spectra using the inverse of the Planck function given byTMOD,n(sfn,Γ,fR,λ)=hc/λkBln(2hc2λ5sfnRMOD,n(Γ,fR,λ,t)+1)(11)where h is the Planck constant, c is the speed of light, and kB is the Boltzmann constant. sfn was an optional scale factor for each spot that optimized the model fits to the spectral shapes (or apparent color temperatures) of the spot radiance spectra. Including this optional scale factor allowed the model fits to account for small modeling errors caused by incorrect assumptions of emissivity, inaccuracies in the shape/topography model used for thermophysical modeling, and/or deficiencies in the hemispherical crater model used to represent the unresolved surface roughness. The scale factor, when used, was varied to optimize the model fit for a given spot for each thermal inertia and roughness value tested. Otherwise, the scale factor had a fixed value of 1. The use of this scale factor was limited to the fitting of the OTES spots because the OVIRS observations had insufficient spectral range to perform model fitting based only on spectral shape.

In general, fitting a single OTES or OVIRS spot could not simultaneously constrain both thermal inertia and roughness due to the strong degeneracy between these two parameters. To circumvent this in single-spot analyses, we fixed the amount of roughness at a specified value and varied just the thermal inertia to find the best model fit using Eqs. 10 and 11. Therefore, the uncertainties of the single-spot thermal inertia values were primarily driven by the uncertainties of the specified input roughness values. Examples of this single-spot fitting are shown in fig. S13 for the lowest and highest thermal inertia regions identified on Bennu. This single-spot fitting was performed in the analysis of the Recon A data described below.

Fitting multiple overlapping OTES or OVIRS spots

For the simultaneous fitting of NSPOT overlapping spots, Eq. 10 was modified toRMSE(Γ,fR)=1NSPOTNλΣn=1NSPOTΣNλ(TOBS,n(λ)TMOD,n(sfn,Γ,fR,λ))2(12)

In particular, suitable combinations of spots to simultaneously fit were sought by finding all OTES or OVIRS spots that projected onto a specified location on Bennu, as indicated by the OSIRIS-REx positional and pointing information. Similar to the single-spot fitting, the simultaneous fitting of multiple spots could be performed with or without the use of optional scale factors to optimize the fit quality. If they were included, then a unique best-fit scale factor was determined for each spot for each thermal inertia and roughness value tested. Therefore, the amount of variation in the derived scale factors for a multiple-spot fit gave an alternative measure of the absolute model fit quality. As in the single-spot fitting, this multiple-spot fitting with multiple scale factors was limited to the analysis of the OTES observations.

Unlike the single-spot fitting, the simultaneous fitting of multiple spots could uniquely constrain both thermal inertia and roughness if the included spots had a large spread in local solar time. This was especially so for day- and nighttime spots that captured the local amplitude of the diurnal temperature variation. Examples of this multiple-spot fitting are shown in fig. S14 for the regions of the lowest and highest thermal inertia identified on Bennu.

The fitting of multiple spots also allowed the uncertainties of the derived thermal inertia and roughness values to be independently quantified through a Monte Carlo method. It was not possible to quantify the uncertainties through a χ2 analysis because measurement errors were not reported by OTES for all the wavelengths that it observed. However, this was not usual because previous work using similar datasets for other planetary bodies also used Monte Carlo methods to quantify uncertainties for similar reasons, e.g., (7). For our Monte Carlo uncertainty analysis, we use the bootstrap method described in (94). In particular, if given a suitable number of independent observations, then the bootstrap method estimates the uncertainties by performing multiple fits of different randomly chosen combinations of those observations. Each different fit gives a slightly different set of best-fit properties, and the resultant property distributions are characteristic of the model uncertainties. For the analysis of the OTES and OVIRS data, this translated to fitting different randomly chosen combinations of the overlapping spots. Therefore, for any given location on Bennu, we produced 1000 different spot combinations from all the spots that projected onto that location and found the best-fit thermal inertia and roughness for each spot combination using Eq. 12. In the bootstrap random spot choosing, some of the spots would be repeated multiple times, while others would be neglected. After the best-fit thermal inertia and roughness were found for each spot combination, we calculated their average value and SD to characterize the overall fit uncertainties. Examples of this Monte Carlo bootstrap uncertainty analysis are shown in fig. S15 for Bennu’s regions of the lowest and highest thermal inertia. These multiple-spot fitting and uncertainty quantifications were performed in the analysis of the Equatorial Stations and BBD 1 data described below.

Global analysis of Equatorial Stations data

To produce the global maps of thermal inertia and surface roughness for Bennu, we analyzed the OTES and OVIRS data collected by OSIRIS-REx from the seven Equatorial Stations using the multiple-spot fitting method described above. In particular, we chose spots with local solar times that fell within 1 hour of each of the respective station times for best spatial resolution (e.g., spots with local solar time between 2 p.m. and 4 p.m. for the 3 p.m. Equatorial Station) and filtered out OTES spots that suffered from a phase inversion issue (i.e., spots where the scene temperatures were too similar to the OTES detector temperature). Table S1 summarizes the number of OTES and OVIRS spots that were selected for use in this analysis.

Surface temperature lookup tables were generated by running the ATPM on three different versions of the 6-m-resolution Bennu shape model to capture any shape sensitivity in our results. These shape model versions were the SPC v20 model derived from imaging only (25), the SPC/OLA v34 model derived from a combination of stereophotoclinometry and laser ranging (40), and the OLA v16 model derived from laser ranging only (41). As described earlier, these lookup tables were built using a Bond albedo of 0.02, bolometric emissivity values of 0.90, 0.95, and 1.00, and thermal inertia ranging from 0 to 800 J m−2 K−1 s−1/2 in steps of 10 J m−2 K−1 s−1/2. Synthetic model radiance spectra were then produced for each OTES and OVIRS spot using these surface temperature lookup tables.

For the multiple-spot fitting, we created lists of OTES and OVIRS spots that projected onto each facet of the 6-m shape model being tested and filtered out facets that had less than 10 spots projected onto them. In particular, a minimum of 10 spots was chosen to ensure that the Monte Carlo bootstrap analysis gave reliable uncertainties. This mainly filtered out polar facets due to the lack of completely filled spots at extremely high latitude, and some equatorial facets in the OVIRS-derived maps due to the occurrence of fewer overlaps between its smaller spots. However, this criterion of more than 10 spots still allowed us to derive thermal properties for 99.6 and 98.2% of Bennu’s surface from the OTES and OVIRS data, respectively. Last, the multiple-spot fitting and uncertainty analysis was performed for each 6-m facet using the lists of overlapping spots, which was separately performed for the OTES and OVIRS data to give two sets of independently derived thermal property maps. This included analyses that kept the scale factors fixed at a value of 1 and analyses that optimized the fits with variable scale factors. Tables S2 and S3 summarize the results of the different global analyses performed on the OTES and OVIRS data, respectively.

As shown in table S2, a mean scale factor of ~1 was only derived from the OTES data when a bolometric emissivity of 0.95 was assumed in the thermophysical modeling. For assumed bolometric emissivity values of 0.90 and 1.00, the mean scale factors were ~1.05 and ~0.95, respectively, which indicated that a bolometric emissivity of 0.95 was the true value for Bennu. However, the same thermal inertia of Bennu was derived despite using these wrong bolometric emissivity values, which demonstrated the usefulness of having variable scale factors. In addition, very similar thermal inertia and roughness values were also derived for Bennu from the three different shape model versions tested. However, the OLA/SPC v34 and OLA v16 shape models did give slightly better fits to the data than the SPC v20 shape model. This was particularly evident on and near large boulders that were not adequately captured in the SPC v20 shape model.

As shown in table S3, the OVIRS analysis gave consistent results between the three different shape models tested, but the derived thermal inertia was consistently 20 J m−2 K−1 s−1/2 higher than that derived by OTES. Figure S2 demonstrates that this was a systematic offset, which can simply be accounted for by subtracting 20 J m−2 K−1 s−1/2 from the OVIRS-derived values. Regardless of this systematic offset, the OVIRS-derived maps are in good spatial agreement with those derived from OTES data (Fig. 1).

Verifying global map uncertainties

We checked that the uncertainties derived by the Monte Carlo bootstrap method were realistic by confirming that the uncertainty-normalized dispersions between the OTES- and OVIRS-derived global maps were normally distributed. For a quantity that was independently derived by both OTES and OVIRS, the uncertainty-normalized dispersion for a given facet, σDIS, was calculated byσDIS=fOVIRSfOTESσOVIRS2+σOTES2(13)where fOTES, fOVIRS and σOTES, σOVIRS were the model derived values and uncertainties, respectively. Figure S4 shows the σDIS distributions determined from the OTES- and OVIRS-derived maps of thermal inertia and surface roughness for Bennu. For thermal inertia, its σDIS distribution had a best-fit Gaussian function of N(0,1) after accounting for the systematic offset in thermal inertia between OTES and OVIRS, which indicated that the estimated thermal inertia uncertainties were realistic. For surface roughness, its σDIS distribution had a best-fit Gaussian function of N(0,1.2), which indicated that the roughness uncertainties were perhaps underestimated by ~20%. However, this was primarily a consequence of the OVIRS-derived roughness range being ~50% larger than that derived by OTES due to its higher spatial resolution. Increasing the OTES-derived roughness uncertainties by ~50% to account for its smaller range in roughness ensured that a N(0,1) Gaussian function could be fit to the roughness σDIS distribution. Therefore, the derived roughness uncertainties were also realistic despite the initial comparison, suggesting a ~20% underestimation.

Comparison of thermal roughness with OLA 20-cm shape model

The surface roughness RMS slope, θ, of the OLA v16 20-cm-facet shape model can be quantified in the same way as thermal roughness usingθ=Σi=1Nϑi2aicosϑiΣi=1Naicosϑi(14)where ϑi is the angular slope (tilt) of facet i measured against a local horizontal surface (33, 36). Here, we have modified the definition in (36) to include overhangs using the absolute value returned by cos ϑi. Otherwise, the many overhangs captured in the OLA v16 20-cm shape model would have reduced the calculated RMS slope rather than increased it.

To compute spatial variations in OLA-derived RMS slope, we initially registered the facets of the 20-cm shape model to the closest facet of the 6-m shape model. This registering was performed using a nearest-neighbor approach between the central points of the 20-cm and 6-m facets. After all the 20-cm facets had been registered to a suitable 6-m facet, their angular slopes were calculated by comparing their surface normal against the surface normal of the 6-m facet with which they were registered. Equation 14 was then repeated for each 6-m facet using the calculated angular slopes and areas of the 20-cm facets to build up a 6-m resolution map of RMS slope measured at 20-cm scales.

For a direct comparison with the OTES thermal roughness, the 6-m resolution map of OLA RMS slope had to be degraded to match the spatial resolution of OTES, which was performed in two averaging steps. First, an “observed” RMS slope value was calculated for each OTES spot by averaging the RMS slope of all 6-m facets onto which the spots projected. Second, the OTES RMS slopes were then mapped back to the 6-m shape model by averaging the RMS slopes of all OTES spots that projected onto each 6-m facet. As a result, this method produced the degraded-resolution OLA surface roughness map shown in fig. S5, which was in good spatial agreement with the OTES thermal roughness map shown in Fig. 1.

Figure S6 shows a comparison between the surface roughness RMS slopes derived from OTES and OLA. As shown, there was a good correlation between the two datasets, but the relationship was not one to one. In particular, a polynomial fit gave a best-fit cubic function ofθOTES=3.062×104θOLA34.526×102θOLA2+2.332θOLA(15)which can be used to predict the OTES RMS slope from an OLA-derived RMS slope with a typical RMS error of 1.3°. This cubic function was later used in the analysis of the Recon A data described below.

Searching for temperature dependency in OTES BBD 1 data

In our thermophysical analysis of Bennu, we considered that its thermal inertia was constant with both temperature, T, and depth. For the Moon and Mars, although, their fine-particulate regolith gives rise to thermal inertia values that are temperature dependent. This temperature dependency is chiefly driven by thermal conductivity, which is proportional to T3 on an airless body due to the dominance of radiative heat transfer between regolith particles (4, 16). A T3 dependence of thermal conductivity translates to a T3/2 dependence for thermal inertia. Heat capacity is also temperature dependent but has less variation with temperature than radiative thermal conductivity. If a strongly temperature-dependent thermal inertia existed on Bennu, then it could result in the observed trend with latitude because the average temperature on Bennu also decreases with increasing latitude. Temperature-dependent thermal inertia would also result in overpredicted model temperatures for both the day- and nighttime observations (16), but we did not see any evidence for this in our model diurnal temperature curves (Fig. 2). However, these model overpredictions are expected to be subtle, and they could have been masked by surface roughness effects.

Bennu’s relatively high orbital eccentricity provides a more robust way of detecting a temperature-dependent thermal inertia, as its distance to the Sun, rH, varies from 0.90 to 1.36 AU during its orbit. Given that surface temperature is proportional to rH−1/2, it follows that thermal inertia is proportional to rH−3/4 for a T3/2 dependence of thermal inertia. Heliocentric distance changes in thermal inertia have previously been observed for the main-belt asteroid (6) Hebe (95), a few NEAs (87), and groups of trans-Neptunian objects and Centaurs (96). To determine whether Bennu exhibited changes in thermal inertia with heliocentric distance, we analyzed OTES data collected from the mission’s BBD 1. Because these data were collected when Bennu was at a heliocentric distance of ~1 AU rather than at ~1.2 AU during the Equatorial Stations, Bennu’s thermal inertia should have systematically increased by ~50 J m−2 K−1 s−1/2 if a T3/2 temperature dependency was present.

For the analysis of the BBD 1 data, we selected OTES spots that had local solar times that fell within 1 hour of 12:30 p.m. to ensure that they had the best possible overlap with spots collected from the 12:30 p.m. Equatorial Station 3. Temperature lookup tables and model radiance spectra were then produced by the ATPM using the positional geometry of Bennu given in table S1. For the OTES spot fitting, the single local solar time of day of the BBD 1 data meant that we could not independently constrain thermal inertia and surface roughness using it. Therefore, to allow us to constrain thermal inertia, we used the surface roughness map determined from the global analysis of the Equatorial Stations data as an additional input in our multiple-spot fitting procedure. As shown in fig. S16, the resulting BBD 1 thermal inertia map only covered 60.4% of Bennu’s shape model facets, so we reevaluated the thermal inertia trend with latitude by only averaging facets that had values reported in both maps. The resulting latitude trends are shown in fig. S11, and they demonstrate that changes in Bennu’s thermal inertia between the heliocentric distances of ~1.0 and ~1.2 AU were minimal; the data are well fit by a one-to-one line. Although a small offset may exist between the two sets of thermal inertia values (we compute a mean difference of ~2 J m−2 K−1 s−1/2 for the data in fig. S11), it is far below the average derived thermal inertia uncertainty value of ~10 J m−2 K−1 s−1/2. We thus use this uncertainty value to conclude that the thermal inertia change between the two stations is ≲10 J m−2 K−1 s−1/2 or ≲3.7%. Correspondingly, a mean diurnal surface temperature of 260 K at 1.2 AU would be expected to increase to ~285 K, an increase of 25 K. We can compare this increase in thermal inertia with temperature to the well-studied Cold Bokkeveld CM meteorite (27, 28, 51), which would have thermal inertia values of ~1030 and 1090 J m−2 K−1 s−1/2 at 260 and 285 K, respectively (a 6% increase). Some other CM meteorites measured in (28) show an even weaker temperature dependence of thermal inertia in this temperature range (20 to 30 J m−2 K−1 s−1/2). Although these meteorites are perhaps not ideal specimens for comparison given that their thermal inertia values are much higher than Bennu’s average, they serve to show that the estimated upper limit for temperature dependence on Bennu is comparable to experimental values for chondritic meteorites.

To make a contrasting comparison, we can estimate the temperature dependence of a regolith with thermal inertia that is equivalent to Bennu’s global average 300 J m−2 K−1 s−1/2 at 260 K. We use the model in (45) with the addition of the nonisothermality correction factor developed in (46) and find that a thermal inertia of 300 J m−2 K−1 s−1/2 corresponds to a particle size of approximately 3 cm in diameter. See “Particle size estimates” section below for a list of model parameters, with the only change being that, here, we use the temperature-dependent heat capacity for Cold Bokkeveld from (28). We can then calculate the thermal inertia of such a regolith at 285 K and find that it increases to ~350 J m−2 K−1 s−1/2, a ~17% increase. Given that we know that Bennu is not predominantly covered in fine-grained regolith, it is expected that this level of temperature dependence is not observed.

The average diurnal surface temperatures of Bennu, from model fits to Equatorial Stations data, decrease by up to ~25 K from the equator to the high latitudes. Given that this is roughly equivalent to the temperature difference between the Equatorial Stations and BBD 1 mission phases, we should expect that thermal inertia decreases due to temperature dependence cannot be more than ~10 J m−2 K−1 s−1/2. The reported thermal inertia latitudinal trend (Fig. 5) shows a decrease of up to ~40 to 60 J m−2 K−1 s−1/2 from the equator to high latitudes. Therefore, we conclude that the latitude trend is a real variation in the physical properties of Bennu’s surface; less than ~25% of the trend could be attributed to temperature-dependent material properties.

Local analysis of OTES Recon A data

To produce local maps of thermal inertia for the Nightingale and Osprey candidate sample sites, we analyzed the OTES data collected by OSIRIS-REx during the low-altitude Recon A passes using the single-spot fitting method described above. For generating suitable temperature lookup tables, we carved out subsections from the 1.6-m-resolution shape model that were 80 m in radii about the central points of the Recon A scan patterns. Temperature lookup tables were then built by the ATPM using these shape model subsections with the positional geometry of Bennu given in table S1. To avoid boundary issues with the edge of these shape model subsections, we only analyzed and produced model radiance spectra for OTES spots that fell within 75-m radii of the subsection central points. Most facets were only observed a few times at a single local solar time of day despite there being a ~7-hour range for the collected OTES spots. Therefore, it was not possible to independently constrain thermal inertia and surface roughness from the Recon A OTES data alone. To circumvent this, we measured the surface roughness of the OLA-derived v16 20-cm shape model for these regions by comparing its facets against those of the 1.6-m shape model subsections using Eq. 14. This allowed us to calculate θOLA for each OTES spot, which could then be converted to θOTES using the previously derived cubic function given in Eq. 15. Last, these calculated roughness values were given as an input to the single-spot fitting procedure to determine a best-fit thermal inertia for each OTES spot analyzed. This resulted in an average thermal inertia uncertainty of ~20 J m−2 K−1 s−1/2 for each OTES spot due to the uncertainty in the measurement and conversion of the OLA roughness.

To visualize the OTES spot thermal inertias, we projected their values back onto the shape model subsection facets and averaged the values for facets that had multiple spots projected onto them. The resulting local thermal inertia maps are shown in Fig. 4. For comparison with the global analysis, we characterized the local thermal inertia distributions by binning the values that were within 20-m radii of the candidate sample site coordinates (i.e., the spatial resolution of OTES during the Equatorial Stations). These local thermal inertia distributions are shown in fig. S7, and they agree well with the corresponding thermal inertia values derived from the global analysis.

Dust layering analysis

To investigate whether thin dust layers that could lower thermal inertia are present on Bennu, we implemented the dust layering methodology outlined in (47) into the ATPM analysis. In modeling the temperatures of dust layers, Eqs. 1 and 3 were solved in terms of actual physical depth, rather than normalized depth, to allow k, ρ, and CP to take different quantities at different depths. This modification was applied to both the shape model facets and the roughness crater subfacets in the ATPM temperature modeling. For the dust layer, we assumed a lunar-like thermal inertia of 50 J m−2 K−1 s−1/2 (1), which was specified in the ATPM by setting k, ρ, and CP to 0.0028 W m−1 K−1, 1190 kg m−3, and 750 J kg−1 K−1, respectively. A variable thickness of this dust layer was then placed on top of rocks with assumed thermal inertia values of 200, 300, 400, and 800 J m−2 K−1 s−1/2, which were specified in the ATPM by setting k, ρ, and CP to the quantities summarized in table S4. The thickness of the dust layer was then varied from 50 μm to 1 mm in 50-μm steps, and the subsurface temperatures of the underlying rock were computed down to a maximum physical depth of 1 m.

For these thin dust layers, solving Eqs. 1 and 3 with the finite difference method required much smaller time steps to obtain a stable solution (10,000 to 20,000 time steps per revolution were needed), which increased the ATPM memory requirements. Therefore, to reduce the ATPM memory cost, we investigated the potential presence of dust layers for just the regions with the lowest and highest thermal inertia, rather than the whole of Bennu, by carving out their subsections from the 6-m shape model. Temperature lookup tables were built using the methodology and physical properties outlined in this section, and model radiance spectra were generated for the OTES spots that projected onto these end-member thermal inertia regions. The resulting diurnal temperature curves are shown in fig. S8.

Lateral dust mixing analysis

To determine whether a fractional mixture of dust with very low thermal inertia on boulders with high thermal inertia could explain the low thermal inertia of Bennu’s boulders, we fitted a two-component model to OVIRS data of the regions with the lowest and highest thermal inertia identified earlier. OVIRS data were used for this investigation, rather than OTES data, because the short wavelengths observed by OVIRS were more sensitive to the potential hot day-side temperatures induced by the presence of any very low thermal inertia dust.

In the two-component model fits, the total radiance, RTOT,nROCKDUST, fDUST, fR, λ, t), was calculated from the single-component model radiances byRTOT,n(ΓROCK,ΓDUST,fDUST,fR,λ,t)=(1fDUST)RMOD,n(ΓROCK,fR,λ,t)+fDUSTRMOD,n(ΓDUST,fR,λ,t)(16)where fDUST was the fractional surface coverage of dust and ΓROCK and ΓDUST were the rock and dust thermal inertia values, respectively. As before, this total radiance could be converted to brightness temperatures using Eq. 11, which could be fit using the multiple-spot fitting method. For the two end-member thermal inertia regions, we adopted a lunar-like dust thermal inertia of 50 J m−2 K−1 s−1/2 (1) and calculated the quality of model fit and the best-fit rock thermal inertia for the full allowable range of dust surface cover (fig. S9). To avoid lengthy searches through a 4D parameter domain, we fixed the surface roughness at the values determined from the single-component fits.

Albedo trend analysis

We began with a photometrically corrected MapCam v-band equirectangular mosaic that has been calibrated to units of normal reflectance (97). A Gaussian convolution kernel is applied to the map to smear it to approximate OTES spectrometer resolution during the Equatorial Stations. The size of the kernel was set to approximate the OTES spot sizes. The resulting smeared maps were then facetized to the SPC/OLA v34 6-m shape model of Bennu by mapping each pixel centroid to the nearest facet and averaging the binned values. These facetized albedo maps were plotted against the corresponding OTES thermal inertia map from Fig. 1 to generate the histogram data shown in Fig. 3.

Particle size estimates

To estimate the effective regolith particle size in Nightingale crater, we use the model in (45) with the addition of the nonisothermality correction factor developed in (46). The model in (45) for the thermal conductivity of regolith is first used as described in the paper with the following input parameters used in (46): Macroporosity varies from 0.40 to 0.70; specific heat capacity of 750 J kg−1 K−1 (28); particle density is a function of thermal conductivity, using the average of two models (50, 51) for meteorite thermal conductivity versus microporosity and assuming an average CM grain density of 2960 kg m−3 (52); surface energy of 0.032 J m−2 (98); bolometric emissivity of 0.95; Bennu mid-latitudes surface acceleration of 40 μm s−2 (53); Poisson’s ratio of 0.27; Young’s modulus of 5.63 GPa; mean burial depth of 10 mm; Sakatani model tunable solid contact conductivity parameter ξ = 0.63 (98); and tunable radiative conductivity parameter ζ = 0.85 (98).

The model is then augmented from its form as described in (45) by multiplying the calculated regolith radiative thermal conductivity by a nonisothermality correction factor, fk, following the methods described in (46). A lookup table of regolith thermal conductivity and thermal inertia is calculated as a function of regolith particle size, macroporosity, and the thermal conductivity of the individual particles. This lookup table is then used to determine the range of possible particle sizes for a desired value of thermal inertia.

The thermal inertia of nightingale crater is approximately 150 to 200 J m−2 K−1 s−1/2 in the OTES Recon A thermal inertia map (Fig. 4). When the thermal conductivity of the particle material is set to 0.60 W m−1 K−1 (similar to CM meteorite Cold Bokkeveld) (27), the predicted regolith particle diameter for a thermal inertia of 150 J m−2 K−1 s−1/2 is approximately 0.68 to 1.0 cm. When the thermal inertia is 200 J m−2 K−1 s−1/2, the predicted particle diameter is approximately 1.2 to 1.8 cm.

If the assumed thermal conductivity of the individual particles is lowered to a value below ~0.1 W m−1 K−1 for the case of 200 J m−2 K−1 s−1/2, then the predicted particle diameter exceeds the length scale of the diurnal skin depth. In these cases where the particle size prediction exceeds the diurnal skin depth, it is assumed that the model is unconstrained. Such a regolith should have the appearance of bedrock in thermal measurements.

Calculation of rock porosity and tensile strength

We follow the methods described in (31) for the estimation of rock porosity and tensile strength for a given thermal conductivity. We assume a specific heat capacity of 750 J kg−1 K−1, based on measurements of CM meteorites (28), a grain density of 2960 kg m−3, based on average density values for CM meteorites (52). Grott et al. (31) used two models to relate rock thermal conductivity to porosity (50, 51). We use the average predictions from the two models to give approximate estimates of the expected rock porosity values (fig. S10). For the calculation of rock tensile strength, we adopt the same parameters as (31), i.e., a Young’s modulus of 10 GPa and a solid thermal conductivity of 2.95 W m−1 K−1. As before, we use the average predictions from the two porosity models to give approximate estimates of the expected rock tensile strengths (fig. S10).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We are grateful to the entire OSIRIS-REx Team for making the encounter with Bennu possible, C. Wolner for editorial help, and H. Roper and D. Worden for help with figure formatting. We thank the J-Asteroid/JMARS development team for assistance with the production of map figures. Funding: This material is based on work supported by NASA under contract NNM10AA11C issued through the New Frontiers Program. Some of the research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. B.Ro. acknowledges funding support from the Royal Astronomical Society (RAS) and the U.K. Science and Technology Facilities Council (STFC). A.J.R., M.B., M.D., and D.P.M. acknowledge support from the Academies of Excellence on Complex Systems and Space, Environment, Risk and Resilience of the Initiative d’EXcellence (IDEX) Joint, Excellent, and Dynamic Initiative (JEDI) of the Université Côte d’Azur and from the Centre National d’Études Spatiales (CNES). The work of M.D. was partially supported by ORIGINS of the Agence National de la Recherche (ANR-18-CE31-0014). J.L.B., C.M.E., R.D.H., J.L.M., and M.A.S. were supported by the OSIRIS-REx Participating Scientist Program. Author contributions: B.Ro., A.J.R., and J.P.E. led the thermophysical analysis and interpretation of the OSIRIS-REx infrared data of Bennu. P.R.C. and V.E.H. are the lead instrument scientists for OTES supported by C.W.H. A.A.S. and D.C.R. are the lead instrument scientists for OVIRS. M.G.D. is the lead instrument scientist for OLA supported by M.A.A., O.S.B., and H.C.M.S. B.Ri. is the lead instrument scientist for the OSIRIS-REx Camera Suite. R.-L.B., J.L.B., C.A.B., M.B., K.N.B., S.C., B.E.C., M.D., D.N.D., C.M.E., R.D.H., E.S.H., D.R.G., E.R.J., H.H.K., L.F.L., J.L.M., D.P.M., M.A.S., and K.J.W. contributed to the thermophysical analysis and interpretation of results. M.C.N. and D.S.L. provided scientific leadership and oversight. All authors contributed to the paper writing. Competing interests: The authors declare that they have no competing interests. Data and materials availability: OTES and OVIRS data from the Detailed Survey (Baseball Diamond and Equatorial Stations) and Recon phases of the OSIRIS-REx mission are available via the Planetary Data System (PDS) at and, respectively. OCAMS (PolyCam) images from the Recon phase are available on the PDS at OLA data used to derive roughness for the Recon phase thermal analysis are available via the PDS at Data are delivered to the PDS according to the schedule in the OSIRIS-REx Data Management Plan, available in the OSIRIS-REx mission bundle at

Stay Connected to Science Advances

Navigate This Article