High-Resolution Thermophysical Analysis of the OSIRIS-REx Sample Site and Three Other Regions of Interest on Bennu

,

understanding of carbonaceous asteroids and to identify and characterize the most suitable sites on the surface for sample collection.
The first global imaging of Bennu by OSIRIS-REx revealed a surface covered in boulders, with no obvious fine-grained regolith covering any substantial portion of the surface Lauretta et al., 2019;Walsh et al., 2019). Spacecraft spectroscopic observations  confirmed predictions (Clark et al., 2011) that Bennu is a B-type asteroid composed of material similar to CM or CI carbonaceous chondrites, with abundant hydrated phyllosilicates. Shape and density measurements Scheeres et al., 2019) confirmed expectations that Bennu is a rubble-pile asteroid-that is, a microgravity aggregate of unconsolidated material (Walsh, 2018).
The boulder-rich surface, however, was unexpected. Ground-based radar polarization observations of Bennu had suggested a smooth surface at the scale of the radar wavelength (∼3 to 13 cm; Nolan et al., 2013), and thermal-infrared observations had indicated a global average thermal inertia of only 310 ± 70 J m −2 K −1 s −1/2 , the most likely explanation of which was mm-to-cm-sized regolith particles (Emery et al., 2014). To explain this discrepancy, Nolan et al. (2019) speculated that the ground-based radar scattering was dominated by the surface textures of Bennu's boulders rather than the grain size of pebbles.
Thermal inertia is the resistance of a material to changes in temperature and is given by Γ = (ρkC P ) 1/2 , where ρ is density, k is thermal conductivity, and C P is the heat capacity. Values of thermal inertia smaller than the solid rock value-about 1030 J m −2 K −1 s −1/2 for CM meteorites analogous to Bennu (Opeil et al., 2010(Opeil et al., , 2020)indicate disruption of the thermal wave at scales smaller than the thermal skin depth (d S = [kP/πρC P ] 1/2 , where P is the rotation period). That disruption is usually caused by particle sizes being smaller than d S (Gundlach & Blum, 2013;Ryan et al., 2020;Sakatani et al., 2017). Despite the boulder-covered surface apparent in images, global thermal observations by the OSIRIS-REx Thermal Emission Spectrometer (OTES; 6-50 μm; Christensen et al., 2018) and the OSIRIS-REx Visible and near-InfraRed Spectrometer (OVIRS; 0.4-4 μm; Reuter et al., 2018) during the spacecraft's approach confirmed the low global average thermal inertia of Bennu (350 ± 20 J m −2 K −1 s −1/2 ), with minimal variation with rotation , suggesting that the boulders themselves are thermally distinct.
To systematically narrow in on the best potential sites on Bennu to collect a sample, observations by the OSIRIS-REx spacecraft were divided into several phases that mapped the surface (or parts of the surface) at increasing spatial resolution (Lauretta et al., 2021). The Detailed Survey phase provided global imaging at ∼5 to 25 cm/pixel, including stereo, and global spectroscopy and thermal observations at ∼20 m (OVIRS) and ∼40 m (OTES) spatial resolution at seven 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.). The global coverage of thermal observations with such broad sampling of the diurnal temperature curve provided global maps of thermal inertia and thermal roughness ( Figure 1; Rozitis, Ryan et al., 2020). Global topographic data were collected by the OSIRIS-REx Laser Altimeter (OLA; Daly et al., 2017) during orbital phases. From these data sets, four sites situated within craters were selected as the most promising for sample collection, dubbed Nightingale, Osprey, Kingfisher, and Sandpiper (Lauretta et al., 2021). These sites were observed more closely during the Reconnaissance (Recon) phase of the mission. In the Recon A subphase, all four sites were observed from an altitude of ∼1 km. Based on analysis of those data, Nightingale and Osprey were selected as the primary and backup sample sites, respectively, and were therefore targeted for higher-spatial-resolution observations in the Recon B and C subphases (∼0.6 and ∼0.3 km altitude, respectively). After two rehearsals, OSIRIS-REx collected its sample of regolith from Nightingale . Measured diurnal temperature curves for the four sites from Detailed Survey observations, along with best-fit model curves, are shown in Figure 2, and Table 1 summarizes the properties of each site.
The global mapping from the Detailed Survey mission phase showed heterogeneity among the boulders on Bennu's surface: Lower-reflectance boulders tend to have rougher, more hummocky textures and lower thermal inertias (∼180 to 250 J m −2 K −1 s −1/2 ), and higher-reflectance boulders tend to have smoother, more angular textures and higher thermal inertias (∼400 to 700 J m −2 K −1 s −1/2 ) Rozitis, Ryan et al., 2020). The thermal inertias of both boulder types are much lower than the thermal inertia of analog CM carbonaceous chondrite meteorites (∼1030 J m −2 K −1 s −1/2 ; Opeil et al., 2010Opeil et al., , 2020. The measured diurnal thermal curves are not consistent with the presence of enough dust to significantly lower the thermal inertias 10.1029/2021JE007153 3 of 25 of either boulder type . Spectroscopic observations show very little variation in the phyllosilicate composition across the surface Praet et al., 2021), indicating that the boulders all have similar bulk composition. The low thermal inertias are inferred to be due to high porosity, and therefore low density and strength, of the boulders (∼49% to 55% for the low-reflectance boulders and ∼24% to 38% for the high-reflectance boulders; Rozitis, Ryan et al., 2020). Preliminary analysis of Recon A data of the Nightingale and Osprey sites by Rozitis, Ryan et al. (2020) confirmed these trends at higher spatial resolution. Cambioni et al. (2021) combined machine learning and thermal modeling to separate the thermal contributions of unresolved fine-grained material from those of the boulders in Detailed Survey OTES observations (spatial resolution ∼40 m). They found a similar range of rock thermal inertias and argued that the higher-thermal-inertia boulders produce more regolith from both impact and thermal fracturing than the lower-thermal-inertia boulders, consistent with a greater strength for the higher-thermal-inertia boulders.
The inferred high porosity of Bennu's boulders, and low bulk macroporosity of Bennu as a whole, were independently confirmed by Tricarico et al. (2021) by extrapolating the boulders' size-frequency distribution and comparing their total volume and mass with that of global Bennu. Additionally, Ballouz et al. (2020) showed that Bennu's boulders also have low strength from analysis of mini-craters found on their surfaces.
The JAXA Hayabusa2 mission performed thermal analyses of boulders on the surface of the C-type asteroid (162173) Ryugu. The average thermal inertia of Ryugu is similarly low to that of Bennu (225 ± 45 J m −2 K −1 s −1/2 ; Okada et al., 2020;Shimaki et al., 2020;Sugita et al., 2019). Hayabusa2 deployed the MASCOT lander to the surface of Ryugu; Grott et al. (2019) and Hamm et al. (2020Hamm et al. ( , 2022 report a thermal inertia of the boulder in the field of view of the thermal imager of MASCOT of 256 ± 4 J m −2 K −1 s −1/2 . Imaging with a resolution of ∼0.2 mm/pixel shows no evidence for dust cover, and the diurnal curve is not consistent with dust cover. They conclude that the low thermal inertia is intrinsic to the boulder and is indicative of high porosity (and low density and strength). Grott et al. (2020) also confirmed this inference of high boulder porosity, and low bulk macroporosity of Ryugu, through extrapolation of the size-frequency distributions for boulders observed on Ryugu's surface. Using high-resolution optical and thermal imaging of Ryugu from the Hayabusa2 spacecraft, Sakatani et al. (2021) found several boulders that are darker than the average surface and lower in thermal inertia than the MASCOT boulder. The thermal inertias are as low as ∼50 to 70 J m −2 K −1 s −1/2 , corresponding to porosities of ∼70% to 90%. These boulders were observed in the centers of craters; Sakatani et al. (2021) inferred that they are sub-surface primordial boulders exposed by impacts that have not been compacted at all since their formation.
Most boulders on Bennu, as well as possible patches of fine-grained regolith (particles smaller than the skin depth of ∼1 to 5 cm; Rozitis, Ryan et al., 2020), were not resolvable by the spectrometers during the Detailed  (101955) Bennu, derived from OSIRIS-REx Thermal Emission Spectrometer (OTES) observations collected during the Origins, Spectral Interpretation, Resource Identification, and Security-Regolith Explorer Detailed Survey (adapted from Rozitis, Ryan et al. (2020)). The locations of the four candidate sample sites investigated in this study are indicated by the labeled datapoints.
Survey. No resolved region of fine-grained regolith on either Bennu or Ryugu has yet been analyzed thermally. It therefore remains unclear how the thermal signature of fine-grained regolith produced from the very porous boulders compares to that of the boulders themselves. To investigate whether the distinct properties of the boulders persist at smaller scales, and to search for the thermal signature of fine-grained regolith, here we analyze Figure 2. OSIRIS-REx Thermal Emission Spectrometer (OTES; red) and OSIRIS-REx Visible and near-InfraRed Spectrometer (OVIRS; blue) diurnal temperature curves of the four candidate sample sites obtained from the Origins, Spectral Interpretation, Resource Identification, and Security-Regolith Explorer Detailed Survey (data extracted from Rozitis, Ryan et al. (2020)): (a) Nightingale, (b) Osprey, (c) Sandpiper, and (d) Kingfisher. The datapoints are the 9.24-and 3.8-μm brightness temperatures measured by OTES and OVIRS, respectively, and the lines give the Advanced Thermophysical Model fits to the data using the thermophysical properties given in Table 1 Rozitis, Ryan et al. (2020)) the thermal observations of the four candidate sample sites from the three Reconnaissance subphases (Recon A, Recon B, and Recon C).

Diurnal and Seasonal Temperature Predictions
To model surface and sub-surface temperature distributions at high spatial and temporal resolution, we used the variant of the Advanced Thermophysical Model (ATPM; Rozitis & Green, 2011, 2012, 2013 that has previously been adapted for use on Bennu Rozitis, Ryan et al., 2020). In particular, the ATPM solves the 1-D heat conduction equation with a surface boundary condition to compute temperature T as a function of time t and depth z for each triangular facet of a user-specified digital terrain model (DTM). Lateral heat conduction is ignored because only facets larger than Bennu's diurnal thermal skin depth (∼1 to 5 cm) are considered. The ATPM also assumes that the three components of thermal inertia are constant with temperature and depth, such that the heat conduction and surface boundary condition equations that it solves for each facet are given by (1) respectively. Here, A B is the bolometric Bond albedo, F SUN is the integrated solar flux at 1 AU (i.e., 1,367 W m −2 ), ε B is the bolometric emissivity, and σ is the Stefan-Boltzmann constant. S(t), ψ(t), r H (t), F SCAT (t), and F RAD (t) are functions that evaluate shadowing, cosine of the Sun illumination angle, heliocentric distance, total multiple-scattered sunlight, and total re-absorbed thermal emission, respectively, for each facet at time t. For a given facet i, F SCAT ,i(t) and respectively, where f i,j is the view factor between facets i and j. Informatively, the fractional amount of facet i's sky that is obscured by other facets is given by the total view factor, T VIEW,i , which is calculated by and gives a qualitative indication of how much facet i is exposed to shadowing and self-heating (Rozitis & Green, 2011 is then used to extrapolate the quantities for the continuous orbit. In both cases, the time and depth domains are initialized at the time-averaged temperature, calculated assuming instantaneous equilibrium with absorbed sunlight, and the ATPM is iterated with user-specified thermal inertia(s) until surface temperatures converge to a difference level of 10 −3 K between successive iterations.
For modeling diurnal surface temperature variations, the ATPM adopted Bennu input parameters that were previously used and/or derived in Rozitis, Ryan et al. (2020). These included a bolometric Bond albedo of 0.02 (Li et al., 2021), a bolometric emissivity of 0.95, a rotation period of 4.296 hr, and a pole orientation of λ H = 69.92° and β H = −83.45° . Thermal inertia was allowed to vary between 0 and 800 J m −2 K −1 s −1/2 to fully capture the range in values previously resolved on Bennu . For modeling seasonal sub-surface temperature variations, the adopted Keplerian orbital parameters included a period of 1.196 years, semimajor axis of 1.126 AU, eccentricity of 0.204, inclination of 6.035°, argument of periapsis of 66.223°, and an ascending node of 2.061° (Farnocchia et al., 2021).
The ATPM was applied to the detailed Poisson-reconstructed DTMs of the four candidate sample sites derived by OLA ( Figure 3; Daly et al., 2020). In particular, the v20 and v21 5-cm OLA DTMs were reduced to 20-cm spatial resolution using MeshLab's quadric edge collapse decimation routine (Cignoni et al., 2008). This reduction was necessary because our previous study demonstrated that thermal roughness most correlated with OLA roughness when measured at 20-cm spatial scales . The ATPM also has the ability to add extra unresolved surface roughness, in the form of spherical-section craters, but this option was not utilized because modeling was performed directly at the spatial scale relevant for thermal roughness. However, any residual beaming effects were accounted for by employing an empirical beaming function in the infrared flux calculations, as described in Section 2.2. For the adopted 20-cm DTMs, view factors and roughness statistics were computed using the methods outlined in Rozitis and Green (2011). Figure 3 and Table 2 summarize the resulting DTM properties for the four candidate sample sites. As shown in Figure 3, regions with many boulders and overhangs had the highest total view factors, indicating that they were the most susceptible to the effects of shadowing and self-heating.
As an example application, Figure 4 demonstrates the diurnal and seasonal variations in illumination and temperature that occurred at the TAG contact point, as computed from the Nightingale OLA DTM assuming different thermal inertia values. The rugged terrain within Nightingale produces strong non-isothermal, shadowing, and self-heating effects, which are further exemplified in Movies S1 and S2. These movies also show that Nightingale's high latitude leads to a strong north-south facing asymmetry in surface temperature that is dependent on the azimuthal tilt directions of the DTM facets.

OTES Observations and Model Fitting
Thermal-infrared observations of the four candidate sample sites were collected by OTES during several low-altitude sorties, as summarized in Table 3. OTES is a Fourier transform infrared spectrometer that measures infrared radiance over the wavelength range of 5.71-50 μm (1,750-200 cm −1 ) using ∼2-s integration times (Christensen et al., 2018). It is a point spectrometer with an 8 mrad field of view, and therefore observations of the candidate sample sites were collected during a series of spacecraft slews that were primarily optimized for imaging. Several hundred spots were collected during each sortie, and the geometry and local solar time of observation were computed for each OTES spot. Following the observations, the raw data were run through the instrument pipeline to convert from instrument signal to physical spectral radiance units. As OTES operates by measuring differences in radiance between the observed scene and its detector, the calibration was not robust when the scene radiance crossed the detector radiance, which created a phase inversion issue. OTES spots that suffered from this issue were filtered out, which typically produced a gap in the collected brightness temperature data between ∼290 and ∼310 K.
For comparisons against the ATPM predictions, synthetic radiance spectra for each OTES spot were generated by combining the model radiance from all DTM facets that fell within the OTES field of view at the instances of the individual spots. In particular, the model observed radiance, R MOD,n (Γ,λ,t,α,τ), for spot n at wavelength λ is given by where B(λ,T) is the Planck function that was evaluated using the facet temperatures T i (Γ,t) predicted by the ATPM for a given thermal inertia. Additionally, a i , d i (t), v i (t), and φ i (t) are the facet area, distance to OTES, visibility to OTES (i.e., returns 1 or 0 depending on the presence or absence, respectively, of a clear line-ofsight), and observation angle, respectively, which were calculated from the positional and pointing geometry of OSIRIS-REx at the time of the OTES spot observations. ε(λ) and Ω(t) are the spectral emissivity and spot solid angle, respectively. Finally, bf(λ,α,τ) is an empirical beaming function that accounts for residual beaming effects (i.e., in addition to those implicitly induced by the 20-cm-scale topography), and its functional form is given by  where a(λ), b(λ), c(λ), and d(λ) are wavelength-dependent coefficients that are optimized to remove residual trends that vary with observation phase angle, α, and/or local solar time, τ.
To assess the quality of fits, correlation and root-mean-square-error analyses were performed on the measured and modeled radiance spectra by converting them to brightness temperatures, T OBS and T MOD , using the inverse of the Planck function. In particular, the root-mean-square-errors between the data and model for given thermal inertia values, RMSE(Γ), were calculated using where N SPOT and N λ were the number of OTES spots and wavelengths included in the comparisons, respectively. Furthermore, this quality of fit metric ensured that the infrared spectra measured by OTES were uniformly fit across its entire wavelength range (e.g., see Supporting Information of Rozitis, Ryan et al. (2020)).

Thermal Inertia Interpretation
Thermal inertia features were interpreted in terms of porosity for boulders (diameters > d S ) and/or particle size for fine-grained regolith (diameters < d S ). Following Grott et al. (2019) and Rozitis, Ryan et al. (2020), boulder porosity ϕ B is related to thermal conductivity and thermal inertia via as determined from laboratory studies of meteorites by Henke et al. (2016) and Flynn et al. (2018), respectively. For these calculations, a heat capacity of 750 J kg −1 K −1 at 270 K and a grain density of 2,960 kg m −3 were assumed based on average measurements for CM meteorites by Opeil et al. (2020) and Macke et al. (2011), respectively. Figure S1 in Supporting Information S1 demonstrates the conversions between thermal inertia and boulder porosity for these two relationships. Rozitis, Ryan et al. (2020) averaged the results from the two relationships, but Equation 10 gives better fits to more recent data of very weak CM-like materials (Avdellidou et al., 2020). Therefore, we produced results from both relationships separately, and quoted those from Equation 10 for feature comparison in Sections 4 and 7.
Similarly, for the regolith grain size interpretation, particle diameter d P is related to thermal conductivity and thermal inertia via  where k SOLID (d P ,T,ϕ R ) is the solid thermal conductivity (i.e., heat conduction through solids at the contacts of particles), k RAD (d P ,T,ϕ R ) is the radiative thermal conductivity (i.e., radiation exchange between particles), f K (d P ,T) is a non-isothermal correction factor (i.e., accounts for temperature gradients within particles), and ϕ R is the regolith macroporosity. Gundlach and Blum (2013) and Sakatani et al. (2017) give similar semi-analytical models for calculating k SOLID (d P ,T,ϕ R ) and k RAD (d P ,T,ϕ R ) from multiple regolith properties, and both models were implemented in this work using the non-isothermal correction factor derived by Ryan et al. (2020). To limit the explored parameter space, we assumed CM-like particles and a fixed regolith macroporosity of 50%. Tables S1 and S2 in Supporting Information S1 summarize the input parameter values adopted in the Gundlach and Blum (2013) and Sakatani et al. (2017) models, respectively, and Figure S2 in Supporting Information S1 demonstrates the resulting conversions between thermal inertia and particle diameter for a regolith temperature of 270 K. As shown in Figure S2 in Supporting Information S1, the particle size interpretation was only valid for thermal inertia values less than ∼360 and ∼230 J m −2 K −1 s −1/2 for the Gundlach and Blum (2013) and Sakatani et al. (2017) models, respectively, as the derived particle diameters exceeded the thermal skin depth otherwise. For suitable regions, we produced results from both models separately, and quoted those from the Sakatani et al. (2017) model for feature comparison in Sections 4 and 7.

Nightingale
As Nightingale was the primary sample site for OSIRIS-REx, it was subject to the most scientific scrutiny and characterization. Its three Reconnaissance passes occurred at a spacecraft range of 0.2-1.1 km, and the collected OTES spots ranged from 2.0 to 8.7 m in diameter (Table 3). To begin our analysis, we filtered out OTES spots that did not fully project onto the OLA DTM for Nightingale (Table S3 in Supporting Information S1), as some of the spacecraft slews extended beyond the terrain modeled during the three passes. The remaining OTES spots were modeled using the diurnal-only method outlined in Section 2 with the Bennu geometry given in Table 3, and the resulting brightness temperature predictions were then compared with the OTES measurements. Figure 5 shows the comparisons between the OTES data and the detrended model (i.e., with residual beaming corrections) for the three Nightingale passes, assuming a constant spot thermal inertia of 260 J m −2 K −1 s −1/2 (the value derived from the Detailed Survey data, given in Table 1). Large variations were observed in measured surface temperature across Nightingale's rugged terrain, which strongly correlated with the DTM-based model predictions (Table 4). However, as shown in Figure 6a, detrending the model predictions with Equation 7 was necessary because a residual phase angle trend was apparent in the original and unadjusted comparisons. In particular, the low-phase-angle observations from Recon A were hotter than predicted, and the high-phase-angle observations from Recon B and C were cooler than predicted. This residual trend mostly depended on the observation phase angle but a small dependence on local solar time was apparent too. If the entire region was treated as a single geological unit in the unadjusted model, then an increasing thermal inertia value was derived with increasing observation phase angle (Figure 6c). This was resolved by detrending the model with Equation 7 (Figure 6b), which enabled a repeatable thermal inertia value to be obtained from each of the three passes ( Figure 6d, Table 4). Possible causes of this residual phase angle trend are discussed in Section 6.
Despite accounting for the systematic phase angle trend, there were still some large residuals for a majority of the OTES spots (Figure 6b), which indicated that the spots were sampling heterogeneous thermal inertia at Nightingale. Therefore, thermal inertia variations were quantified by optimizing the model fits to the individual spots and their overlaps. The derived spot thermal inertias were then projected back onto the OLA DTM facets, which were averaged for facets that had multiple spots projected onto them. Figure 7a shows the resulting thermal inertia map of Nightingale produced after combining the Recon A and B data, and it displays more detail than the preliminary Recon A version presented in Rozitis, Ryan et al. (2020). Recon C data were not included in this analysis because the collected spots had an incomplete spatial coverage over the OLA DTM.
Hokioi Crater (described in Barnouin et al., 2022), in which the Nightingale sample site is located, has relatively low thermal inertia that distinguishes it from adjacent surrounding boulders (Figure 7a)

Osprey
Osprey was the backup sample site for OSIRIS-REx, and therefore it was also subject to detailed scientific scrutiny and characterization. Like Nightingale, its three Reconnaissance passes occurred at a spacecraft range of 0.2-1.1 km, and the collected OTES spots ranged from 1.9 to 8.7 m in diameter (Table 3). The same thermophysical analysis performed for Nightingale was applied to the three Osprey datasets. In particular, OTES spots that did not fully project onto the OLA DTM for Osprey were filtered out (Table S3 in Supporting Information S1), and the remaining spots were modeled using the diurnal-only method with the geometry given in Table 3.     (Table 1). There were strong correlations between the data and DTM-based model (Table 4), but there were smaller spatial variations in measured surface temperature when compared with the Nightingale datasets. The smaller spatial variations in surface temperature result from Osprey's near-equatorial location, which minimizes the north/south-facing temperature asymmetry that occurs very strongly at Nightingale's high latitude. Like for Nightingale, it was necessary to detrend the model predictions with Equation 7 because a residual phase angle trend was also apparent in the original and unadjusted comparisons (Figure 9a). This phase angle trend also impacted the thermal inertia values derived separately from each of the three passes (Figure 9c), but the differences were not as extreme as those obtained for Nightingale because of the narrower phase angle range of the Osprey observations. Detrending the residuals with Equation 7 demonstrated that the Osprey phase angle trend was very similar to that previously characterized for Nightingale, and it too enabled a repeatable thermal inertia value to be obtained from each of the three passes ( Figure 9d, Table 4).
Like at Nightingale, there were still some large residuals for the OTES spots after detrending that indicated the presence of heterogeneous thermal inertia at Osprey (Figure 9b). Therefore, a map of thermal inertia was derived from the Recon A and B data by optimizing the model fits to the individual spots and their overlaps, and then by  , and the large white circles in panels (b, c, d) indicate the areas formerly proposed for sampling within the other candidate sample sites (Walsh et al., 2022). High-resolution images of the features labeled with letters in panels (a and b) are provided in Figures 12 and 13, respectively. Boulder porosity and regolith grain size interpretations of these features are also provided in Tables 5 and 6. Additionally, corresponding maps of thermal inertia uncertainty, RMS brightness temperature error, and OSIRIS-REx Thermal Emission Spectrometer (OTES) spot coverage are provided in Figures S3, S4, and S5 in Supporting Information S1, respectively. Figure S6 in Supporting Information S1 illustrates the frequency distributions of the mapped thermal inertia values.

10.1029/2021JE007153
14 of 25 projecting the results back onto the OLA DTM of Osprey. Figure 7b shows the resulting thermal inertia map of Osprey, which displays more detail than the preliminary Recon A version presented in Rozitis, Ryan et al. (2020).
Wuchowsen Crater, in which the Osprey backup site is located, generally has lower thermal inertia than its surroundings, and its center has a distinct low-thermal-inertia feature (i.e., 220 ± 20 J m −2 K −1 s −1/2 ) that is co-located with a patch of dark-toned material (Figure 7b). Surrounding the crater are mostly boulders with relatively high thermal inertia, but a distinct low-thermal-inertia boulder is also present in the upper left.  Figure 5, but for Osprey.

Sandpiper and Kingfisher
Sandpiper and Kingfisher were candidate sample sites for OSIRIS-REx that were not selected as primary or backup, and as such had only one Reconnaissance pass each (during Recon A) that occurred at a spacecraft range of 0.9-1.1 km (Table 3). Therefore, the collected OTES spots ranged from 7.5 to 8.9 m in diameter for these two sites. As in the analyses previously performed for Nightingale and Osprey, the OTES spots that did not fully project onto the OLA DTMs for Sandpiper and Kingfisher were filtered out (Table S3 in Supporting Information S1), and the remaining spots were modeled using the diurnal-only method with the geometries given in Table 3. Figure 10 shows the comparisons between the OTES data and the detrended models assuming constant spot thermal inertia values of 270 and 320 J m −2 K −1 s −1/2 for Sandpiper and Kingfisher, respectively (Table 1). For both candidate sample sites, there were strong correlations between the data and DTM-based model (Table 4), and both exhibited a residual phase angle trend (Figure 11a). However, in these two cases, the residual phase angle trends were much less apparent than those for Nightingale and Osprey because of the narrow range of the observation phase angle near where the offsets are zero. As a result, the residual phase angle trends had minimal impacts on the thermal inertia values derived for Sandpiper and Kingfisher (Figures 11c and 11d, Table 4). However, the models were still detrended to ensure that the analyses were consistent with those previously performed for Nightingale and Osprey.
It was unclear, from the small model residuals for the majority of the OTES spots, whether thermal inertia is substantively heterogeneous at Sandpiper and/or Kingfisher (Figure 11b). To investigate, maps of thermal inertia were derived from their Recon A data by optimizing the model fits to the individual spots and their overlaps, and then by projecting the results back onto the OLA DTMs for Sandpiper and Kingfisher. No strong thermal inertia features were observed at Sandpiper and Kingfisher; only very subtle variations were evident (Figures 7c and 7d). At Sandpiper, there is a small patch of low thermal inertia in the upper left ( Figure 7c) that coincides with a region that is relatively devoid of meter-scale boulders (Figures 3e and 3f), which suggests a local difference in the beaming function and/or a local abundance of unresolved fine-grained material (Cambioni et al., 2021). At Kingfisher, there were no obvious correlations with its surface features (Figure 7d), but this does not necessarily imply a homogenous thermal inertia because the OTES spot size was much larger than the features visible at this candidate sample site.

Local Spatial Variations in Thermal Inertia
The high-spatial-resolution results for the Nightingale and Osprey sites (Figures 7a and 7b) show several variations in thermal inertia that correlate with geologic features (Figures 12 and 13, Tables 5 and 6). Most observations are in good agreement with previous findings from the Detailed Survey and initial Recon A analyses . For example, for Osprey, a low-reflectance boulder to the northwest of Wuchowsen Crater (Figure 13c) has a low thermal inertia of 190 ± 50 J m −2 K −1 s −1/2 (62%-76% porosity) compared with the surrounding terrain. The shape and surface textures of this boulder, and its thermal inertia, are very similar to other low-reflectance, low-thermal-inertia boulders on Bennu, such as Roc Saxum (194 ± 7 J m −2 K −1 s −1/2 with 67%-69% porosity; Rozitis, Ryan et al., 2020). Conversely, the lighter-toned boulder to the north of Wuchowsen Crater (Figure 13a) has an elevated thermal inertia value of 400 ± 30 J m −2 K −1 s −1/2 (43%-48% porosity), which is consistent with the previously described high-reflectance, high-thermal-inertia boulders . Furthermore, this boulder also contains abundant light-toned veins up to ∼1 cm in thickness (Figure 13b), like the carbonate-bearing veins described by Kaplan et al. (2020). One vein is jutting outwards from the rock surface, suggesting that it is more resistant to erosion than the host rock. Although the link between the high-reflectance boulders, thermal inertia, and the occurrence of these veins was inferred by Rozitis, Ryan et al. (2020), this is the first direct thermal inertia measurement that confirms that vein-bearing boulders indeed have thermal inertia values higher than the global average. The presence of these veins may indicate a history of alteration and cementation that could have acted to more strongly bind the grains within the boulder together, acting to decrease porosity and increase density, strength, thermal conductivity, and, as a result, thermal inertia.
At the Nightingale sample site, just to the northwest of Hokioi Crater , we observe a boulder that morphologically and texturally resembles the high-reflectance, high-thermal-inertia boulders observed elsewhere on Bennu (Figure 12a). However, the boulder's thermal inertia value of 330 ± 30 J m −2 K −1 s −1/2 (49%-55% porosity) is very similar to the global average (300 ± 30 J m −2 K −1 s −1/2 ; Rozitis, Ryan et al., 2020), that is, lower than expected for a high-reflectance boulder. This may be due to the presence of visible patches of mm-to-cm scale regolith particles, which could work to suppress the apparent thermal inertia to a value lower than that of the bare boulder. Particles are observed on another boulder with a very low thermal inertia just to the east of Hokioi Crater (Figure 12b). This boulder's thermal inertia of 180 ± 50 J m −2 K −1 s −1/2 (63%-77% porosity) is consistent with its morphology and texture, which resemble those of the typical low-reflectance, low-thermal-inertia boulders. In this case, it is not clear whether the presence of particles on the boulder surface affected its apparent thermal inertia.
The large boulder that is closest to Hokioi Crater (Figure 12c) has a variable but generally intermediate thermal inertia of 250 ± 50 J m −2 K −1 s −1/2 (55%-67% porosity). In lower-resolution images, this boulder has the appearance of a typical low-reflectance, low-thermal-inertia boulder type with a very hummocky surface texture. However, in higher-resolution images, it appears to be composed of clasts of many different sizes (<0.5 m) and could be classified as a breccia. Another rock with a similar texture is found to the north of Hokioi Crater (Figure 12d) that has a relatively high thermal inertia value of 400 ± 30 J m −2 K −1 s −1/2 (43%-48% porosity). Rozitis, Ryan et al. (2020) and DellaGiustina et al. (2020) noted that the low-reflectance, low-thermal-inertia boulders have hummocky or sometimes brecciated surface textures. However, these two breccia-like boulders have thermal inertia values that are somewhat higher than the low-thermal-inertia boulders identified from the global thermal analysis (∼180 to 250 J m −2 K −1 s −1/2 with 61%-70% porosity). This suggests that these boulders should be classified in a third category based on their brecciated appearance and moderate-to-high thermal inertia values. The elevated thermal inertia could be a function of the clasts (e.g., if they are composed of material like the high-reflectance boulders) and/or be indicative of a strong, dense, and/or thermally conductive cement that is holding the clasts together.
The breccia-like boulders have important implications for how the returned sample analysis can be used to constrain the origin and history of Bennu's different boulder types. Although it is already expected that some lithologies in the returned sample may represent relatively rare, exogenic material (DellaGiustina et al., 2021), we must also be prepared for the possibility that some lithologies may be an isolated component of a multi-component (brecciated) boulder type. Furthermore, such lithologies may even be unique to that brecciated setting. For example, the clasts within the brecciated boulders may be composed exclusively of near-surface materials from Bennu's parent body that were gardened by meteoroid impacts. If the impact gardening was extensive, then those near-surface materials may only exist as breccia clasts, with few or no homogenous meter-scale boulders composed of the near-surface lithology remaining on the surface of present-day Bennu.
Both the Nightingale and Osprey sites have low thermal inertia signatures within their respective craters (Figures 7a and 7b). Within Hokioi Crater at the Nightingale site (Figures 12e and 12f), the low thermal inertia signature (190 ± 30 J m −2 K −1 s −1/2 equivalent to 64%-73% porosity and/or 0.8-1.5 cm particle diameter) appears to be more clearly correlated with the relative absence of rocks and boulders, based on particle counts (Burke et al., 2021), a sampleability analysis that found evidence of abundant particles <2 cm in diameter (Walsh et al., 2022), and the success of the OSIRIS-REx TAG sample collection . However, the low thermal inertia values found within the crater in our analyses may not yet represent those of pure fine-grained regolith, given that many rocks larger than the diurnal skin depth (∼1 to 5 cm) are present that could still elevate the thermal inertia. Nonetheless, the thermal data and corresponding thermal inertia values of the Nightingale sample site are the ripest for future quantitative regolith thermal analyses, compared to anywhere else on Bennu.
The low thermal inertia signature in Wuchowsen Crater, where Osprey is located, is more concentrated toward the center of the crater and is well correlated with the presence of a dark-toned patch of material (Figures 13d  and 13e). It has been debated throughout the course of the mission whether this material is fine-grained regolith or the exposed top of a low-reflectance, hummocky boulder. The low thermal inertia of 220 ± 20 J m −2 K −1 s −1/2 is consistent with both scenarios (62%-67% porosity and/or 1.3-1.8 cm particle diameter). However, high-resolution images of this patch (Figures 13d and 13e) show a surface texture that is quite different from the fine-grained regolith observed at the Nightingale site (Figures 12e and 12f), and its appearance is not unlike the surface of Roc Saxum. If this is indeed the top of a buried rock, its occurrence at the center of the crater is similar to that of the buried primordial rocks with very low thermal inertia (∼50 to 70 J m −2 K −1 s −1/2 ) on Ryugu, described by Sakatani et al. (2021), that are apparently exhumed and exposed by impacts. However, the thermal inertia of the dark-toned patch in Wuchowsen Crater is much higher than the values found by Sakatani et al. (2021), suggesting that it may be no different from the other low-reflectance boulders on Bennu.  Table S4 in Supporting Information S1). The average pixel scale is 4 mm/pixel and the average phase angle is 57.4°.

Temperature Conditions at the TAG Contact Point
For the study of primitive asteroidal material in the laboratory-the culmination of the OSIRIS-REx mission-it is important to assess how well the collected samples have been preserved from active space weathering processes. One important process is thermal alteration from solar heating, as several organic chemical compounds of high interest have relatively low breakup temperatures (e.g., ∼300 to 670 K from Table 1 of Delbó & Michel, 2011), and thermal cycling can also lead to rock cracking and breakdown (Delbó et al., 2014;Molaro et al., 2015Molaro et al., , 2017. Therefore, to provide contextual information for the future laboratory analysis of the collected samples, we evaluated the surface and sub-surface temperature conditions of the TAG contact point using the high-resolution thermophysical model constructed for Nightingale. To perform this analysis, it was necessary to run a full seasonal thermophysical model for the sample site because OSIRIS-REx likely collected material from depths up to ∼50 cm within Bennu's sub-surface, in addition to surface material . Temperature variations at this sub-surface depth on Bennu are primarily driven by the seasonal thermal wave rather than by its diurnal cycles, which makes the diurnal-only modeling approximation unsuitable for use in this case. Therefore, the seasonal methodology outlined in Section 2 was applied to the Nightingale OLA DTM and run for 50-800 J m −2 K −1 s −1/2 thermal inertia to cover the full range of values resolved on Bennu. The adopted thermal inertia values were assumed to remain constant with depth to limit the explored parameter space, but it is a possibility that Bennu and other asteroids have some degree of near-surface layering that could also influence the model temperatures (Harris & Drube, 2016;Rozitis et al., 2018). However, thermophysical analysis of the pre-and post-sampling spectra acquired by OTES and OVIRS of the Hokioi Crater finds no notable change in thermal inertia following the TAG sample collection despite excavation to a depth of ∼70 cm . Additionally, our previous work demonstrated that Bennu's observed diurnal temperature variations were well reproduced by thermophysical models that assumed no or very minimal layering (Figure 2; Rozitis, Ryan et al., 2020). Therefore, it was a suitable approximation to assume that Nightingale's thermal inertia was constant with depth. Figure 14 shows the resulting temperature predictions for the TAG contact point computed as a function of thermal inertia and depth. In order to convert normalized skin depth to physical depth in these predictions, we adopted a heat capacity of 750 J kg −1 K −1 (Opeil et al., 2020) and a sub-surface bulk density of 600 kg m −3  in the skin depth size calculations. As indicated, the surface conditions of the TAG contact point are described by minimum and maximum temperatures of 192 ± 4 and 357 ± 3 K, respectively, when assuming a thermal inertia value of 190 ± 30 J m −2 K −1 s −1/2 (i.e., that measured within Hokioi Crater). However, these quantities change to 232 ± 1 and 261 ± 3 K at a depth of 50 cm, respectively, with most of this change occurring within the top few centimeters of the sub-surface. Additionally, the orbit-averaged temperature for all depths was estimated to be 245 ± 2 K, and the maximum diurnal amplitude experienced by the surface was 146 ± 9 K. Finally, the surface temperature at the date and time of TAG was estimated to be 272 ± 2 K.
In comparison to Table 1 of Delbó and Michel (2011), most of the listed organic chemical compounds are stable at the TAG contact point  Table S5 in Supporting Information S1). The average pixel scale is 3 mm/pixel and the average phase angle is 42.0°.  temperatures, particularly at depth, for Bennu's current orbit. For the potential presence of volatile materials, the surface and sub-surface temperatures are much too warm for water ice to be stable for any length of time (Jewitt & Guilbert-Lepoutre, 2012;Schorghofer, 2008), but it is possible that some high-atomic-mass volatile materials would be stable at these temperatures instead, for example, such as the polycyclic aromatic hydrocarbons . In terms of thermal fracturing, the diurnal amplitude of the TAG contact point is comparable to that experienced at Bennu's equator (∼140 K), but this is not surprising given Bennu's rugged shape . The amount of thermal fracturing experienced by the collected samples would ultimately depend on their shape, size, porosity, and strength (Cambioni et al., 2021;Molaro, Hergenrother, et al., 2020;. However, the diurnal amplitude of the TAG contact point is as high as it can be on Bennu, which implies that it could have been a significantly active process .

Sample site Feature
An important caveat is that these temperature predictions are applicable to Bennu's current orbit and the state of Nightingale's surface and sub-surface immediately before sampling. There are several factors that may have produced higher surface and sub-surface temperatures in the past. The influence of close planetary encounters precludes precise determination of the past orbits of near-Earth asteroids by numerical integration, preventing Sample site Feature Thermal inertia (J m −2 K −1 s −1/2 ) G13 particle diameter (cm) S17 particle diameter (cm) Nightingale e 190 ± 30 0.5-0.9 0.8-1.5 Osprey d 220 ± 20 0.8-1.1 1.3-1.8 Note. The G13 and S17 particle diameter models are those by Gundlach and Blum (2013) and Sakatani et al. (2017), respectively. Values from S17 are quoted in Sections 4 and 7.  a definitive assessment of whether Bennu's orbit has brought it significantly closer to the Sun, and therefore exposed it to higher temperatures. However, Delbó and Michel (2011) used the dynamical history of objects originating in the main-belt that could subsequently reach Bennu's orbit to perform a statistical study. They determined the probabilities of different perihelion distances, and used a thermophysical model to estimate a 50% probability that Bennu's surface has been heated to a temperature greater than 500 K. Although the probability of reaching such temperatures falls rapidly just a few centimeters below the surface, it is also possible that surface material exposed to high temperatures could now reside below the surface because of past impact gardening (Hörz & Cintala, 1997;Housen et al., 1979). The YORP effect (i.e., a torque resulting from asymmetric thermal emission and scattering of sunlight from an irregularly shaped body; Rubincam, 2000) can significantly change the spin-rate and pole orientation of km-sized near-Earth asteroids on timescales of ∼10 6 years, comparable with their orbital evolution (Bottke et al., 1996), resulting in changes in solar illumination conditions. The YORP effect can also potentially cause re-surfacing if the spin-rate approaches the spin-stability limit (Scheeres, 2015).
Although it is not possible to quantify the influence of these effects on the temperature history of the Nightingale sample site on Bennu, the reference temperatures determined in this study provide a suitable starting point to guide the future laboratory analysis of the collected samples.

Thermophysical Model Limitations and Possible Improvements
As exemplified in Figures 5, 8 and 10, Bennu's surface temperature at small spatial scales depends very strongly on its rugged terrain, which made it challenging to extract reliable information on thermal inertia via thermophysical modeling. The most diagnostic measure of thermal inertia is observing the change in surface temperature between night and day (i.e., the amplitude of the diurnal temperature variation), but with only daytime temperature data, the interpretation is very sensitive to the exact details of the topography and observation geometry. Fortunately, the OLA instrument onboard the OSIRIS-REx spacecraft allowed construction of incredibly detailed DTMs (Daly et al., 2020), which enabled the thermophysical models that we implemented to be sufficiently accurate for thermal inertia determination. However, there are drawbacks with our approach that limit the ultimate accuracy of our derived thermal inertia values and maps, the most notable being the residual phase angle trends of the modeled data and the non-uniform spatial coverage of the OTES spots.
In our previous study , we found that thermal roughness most correlated with OLA roughness when measured at 20-cm spatial scales, and therefore we directly modeled surface temperatures on 20-cm-scale topography in this work. Indeed, we found that the 20-cm-scale terrain captured the topographic variations in surface temperature very well, but it was insufficient to directly account for the residual phase angle trends demonstrated in Figures 6a, 9a and 11a. Unmodeled topography at sub-20-cm scales is likely to be the cause of these residual phase angle trends, as Bennu's diurnal skin depth of ∼1 to 5 cm would introduce some non-isothermal temperatures at spatial scales of ∼1 to 20 cm. Therefore, it would be logical to model the surface temperatures directly at the cm-scale instead, but this would be very computationally demanding. For instance, a 3-D thermophysical model that accounts for lateral heat conduction would be required (e.g., Davidsson & Rickman, 2014), along with an enormous view factor table, because the 1-D heat conduction approximation would not necessarily still apply at the cm-scale. Without going to such extreme detail, we were able to account for the residual phase angle trends by simply adopting an empirical beaming function (Equation 7). The measured beaming functions were very similar between the four sites investigated ( Figure S7 in Supporting Information S1), but it was unclear whether they varied within the sites. It would require spacecraft observations of specific geologic units acquired at multiple different observation geometries and local solar times of day to fully disentangle local beaming effects from thermal inertia. Alternatively, thermal emission phase function observations of the returned samples could be performed in the laboratory using a goniometer (e.g., Warren et al., 2017) to inform any necessary improvements to the model implemented here.
Interestingly, a similar residual phase angle trend was reported in Sakatani et al. (2021) from infrared images acquired by Hayabusa2 during descents to Ryugu's surface, but no attempts were made to correct for it. However, the Hayabusa2 observations were acquired at a phase angle of ∼30°, which is close to the ∼37° given by the OSIRIS-REx Reconnaissance observations, where the phase angle offsets are effectively zero. Therefore, it is unlikely that the thermophysical analysis performed by Sakatani et al. (2021) was strongly affected by not correcting for the trend.
Alternatively, we considered whether the phase angle trends could be attributed to temperature-dependent changes in thermal inertia caused by changes in the heliocentric distance of Bennu between the different Reconnaissance passes. However, the thermal inertia results from the unadjusted models do not monotonically increase with decreasing heliocentric distance (Tables 3 and 4), as would be expected from such a temperature dependence Rozitis et al., 2018). This affirms our finding that the unadjusted model residuals are most likely attributable to a phase angle-induced roughness effect, as discussed above.
Finally, the other limitation in our study was the non-uniform spatial coverage of the OTES spots. In particular, the maps shown in Figure 7 imply that we had continuous coverage over the sample sites, but some of the distinct features that appeared in imaging did not necessarily stand out in their mapped thermal inertia. An obvious example is Strix Saxum on the northern rim of Wuchowsen Crater at the Osprey site, visible in Figure 7b. Strix Saxum is a high-reflectance boulder and therefore expected to have relatively high thermal inertia, but its mapped value was not different from that of its immediate surroundings. However, as shown in the OTES spot coverage map in Figure S5b in Supporting Information S1, no OTES spots landed directly on that boulder; they only partially overlapped it. Consequently, its mapped thermal inertia value was primarily derived from its immediate surrounding area. A similar situation arises for the Kingfisher candidate sampling site, as shown in Figure 7d and Figure S5d in Supporting Information S1. Therefore, for future studies, accurate interpretations of the thermal inertia maps shown in Figure 7 should consider the spatial resolution and coverage of the OTES datasets that went into creating them ( Figure S5 in Supporting Information S1). As a mitigation, it might be possible to apply a two-component model to the OTES data to extract sub-spot detail, such as previously done for the Detailed Survey data by Cambioni et al. (2021), but this would require accurate modeling of the residual phase angle trends to prevent degeneracies in the results from occurring.

Summary and Conclusions
In our study, we find that brightness temperatures acquired of Bennu's surface at spatial scales of 2-9 m depend strongly on small-scale topography, local variations in thermal inertia, and the observation phase angle. By combining the ATPM with OLA-derived DTMs and an empirical beaming function, we were able to map spatial variations in thermal inertia for the site sampled by OSIRIS-REx and three other candidate sample sites. These local maps of thermal inertia identify boulders with both low (∼180-250 J m −2 K −1 s −1/2 with 61%-70% porosity) and high (≥400 J m −2 K −1 s −1/2 with ≤45% porosity) thermal inertia, consistent with the initial findings of Rozitis, Ryan et al. (2020). The maps also confirm that vein-bearing boulders (Kaplan et al., 2020) have high thermal inertia, as previously inferred by Rozitis, Ryan et al. (2020). However, the observations of some boulders with moderate thermal inertia (∼250 to 400 J m −2 K −1 s −1/2 with 45%-61% porosity) suggests that Bennu's boulders could be classified into three or more types, consistent with the finding of a boulder thermal inertia continuum by Cambioni et al. (2021).
Additionally, the thermal inertia of the site sampled by OSIRIS-REx, Nightingale, is 190 ± 30 J m −2 K −1 s −1/2 (equivalent to 64%-73% porosity and/or 0.8-1.5 cm particle diameter), which represents an average value for the mixed material within Hokioi Crater prior to sampling. This value appears consistent with sampleability metrics derived from in situ imaging (Walsh et al., 2022) and the interaction of the OSIRIS-REx spacecraft with Bennu's surface during sample collection . In particular, ∼47% of the sample site area was covered with particles <2 cm in diameter (Walsh et al., 2022), and the OSIRIS-REx spacecraft made contact with multiple rocks and mobilized a small quantity of dust . Therefore, the sample site likely contained a complex mixture of materials with both very low (<150 J m −2 K −1 s −1/2 ) and low-to-high (200-700 J m −2 K −1 s −1/2 ) thermal inertia, which is consistent with the derived average value. Disentangling the thermal properties of the different components is scope for future work.
Finally, given that the material at the Nightingale sample site within Hokioi Crater is likely derived from and/ or related to the boulders immediately surrounding it , we predict that the sample being returned by OSIRIS-REx will have unique thermophysical properties in comparison to the existing meteorite collection. In particular, the Nightingale boulders all have thermal inertia (∼200 to 400 J m −2 K −1 s −1/2 ) lower than that measured for analog CM carbonaceous chondrite meteorites (∼1030 J m −2 K −1 s −1/2 ; Opeil et al., 2010Opeil et al., , 2020, which implies that they would not survive atmospheric entry due to their inferred weak tensile strength . Furthermore, the maximum temperatures of the TAG contact point were