Thermophysical Properties and Surface Heterogeneity of Landing Sites on Mars From Overlapping Thermal Emission Imaging System (THEMIS) Observations

Orbitally derived thermal inertia (TI) values of surfaces allow for remote interpretation of rock and sediment physical characteristics. The evolving local times of the Mars Odyssey Thermal Emission Imaging System (THEMIS) mission have enabled surface temperature data collection over multiple seasons and local times over ∼9 Mars years (MY). We utilize this higher temporal resolution data set to separate TI values of individual materials within THEMIS pixels (100 m sampling). In this study, we focus on geologic units within Gusev, Gale, and Jezero crater landing sites to determine their respective TI and ground‐truth our methods. We use the KRC model to predict temperatures for a range of homogeneous and two‐component thermophysical mixing scenarios (laterally and vertically heterogeneous), and compare those to THEMIS brightness temperatures. The Gusev and Gale crater results are consistent with rover‐based evaluations, indurated or clast‐covered sands. The best‐fit scenarios along the Jezero crater volcanic floor unit show low to moderate TI values representative of coarse sediment, volcaniclastic, or pyroclastic rocks. The light toned unit and western fan deposits both indicate sands and a moderate TI component; we interpret this as heavy fracturing within the rock. Difficulties in modeling some THEMIS temperatures are attributed to daily weather effects, local atmospheric dust variability, and real surface changes over time (e.g., dust deposition and removal); we also observe that some temperature observations lead to non‐unique modeled physical solutions. Nevertheless, these methods can still rule out a significant range of material properties and provide meaningful geologic information about Mars’ surface.

which include but are not limited to: fluvial, lacustrine and deltaic sedimentation, volcanism, aeolian transport and erosion, and climatic and atmospheric evolution. The physical properties of sediment (i.e., sorting, grain size, induration) and rock (i.e., porosity, vesicularity, cementation quality) reflect the processes that formed them; thus, constraining those properties is an important goal of remote geological characterization. Thermal inertia (TI) is a material property that describes a material's resistance to changing temperature, is highly sensitive to physical property variations (Section 2.1, Jakosky, 1986;Kieffer et al., 1977;Presley & Christensen, 1997a, 1997bPiqueux & Christensen, 2009a, 2009bWechsler & Glaser, 1965;Wechsler et al., 1972). TI is derived from orbital or lander/rover temperature observations in combination with a thermal model and has been used in many previous studies to provide insight into the physical nature of surface materials (Section 2.1, 2.3). However, subpixel mixing of materials (e.g., sand and outcrops) and general surface complexity can greatly complicate the interpretation of TI values from orbit, particularly when those values are derived using temperature observations from a single time of day, or during the daytime when surface albedo is not well constrained and introduces larger error in modeled temperatures (Putzig & Mellon, 2007a). Temperature observations of the same surface at multiple local times, however, can help to isolate the likely TIs of individual surface components and search for evidence of near surface layering, such as duricrusts, desert pavements, or sediment cover over bedrock (Audouard et al., 2014;Edwards et al., 2018;Mellon, 2007a, 2007b;Vasavada et al., 2017).
The 2001 Mars Odyssey mission carries a nadir-viewing thermal multispectral imager, the Thermal Emission Imaging System (THEMIS), that has been collecting temperature observations of Mars for over 9 Mars years (MY) across multiple seasons and local times. Here, we take advantage of the longevity of the Mars Odyssey mission and use years of overlapping THEMIS images to observe surface temperature changes over diurnal and seasonal cycles. We then model surface temperatures using a series of homogeneous, lateral mixing, and vertical layering scenarios to better constrain the thermophysical properties of sites of interest (Kieffer, 2013). Specifically, we focus on characterizing the physical nature of geologic units in past, present and future landing sites where there is an abundance of data and a means to ground-truth the methods presented herein using rover observations. Key questions addressed include: (a)

Thermal Inertia
Both diurnal and seasonal surface temperature cycles are mostly governed by the physical properties of the materials that make up the surface and near subsurface and their response to solar heating. Physical properties are quantitatively related to surface temperature through the thermal inertia (TI) parameter. TI (measured in J m −2 K −1 s −1/2 ), is an intrinsic property of a material that describes how efficiently that material can store, conduct, and reradiate heat. It is given by: (k: thermal conductivity, W/mK; ρ: bulk density, kg/m 3 ; and c: specific heat, J/kgK). In other words, TI describes the ability of the upper several centimeters of a planetary surface to resist change in kinetic temperature and is not dependent on latitude, time of day, or season. It is dependent solely on the physical nature of the materials themselves (Kieffer et al., 1977;Presley & Christensen, 1997a, 2009a, 2009bWechsler & Glaser, 1965;Wechsler et al., 1972). Temperature fluctuation in the near subsurface decreases exponentially with depth until it reaches a value of 1/e at a depth known as the thermal skin depth (Turcotte & Schubert, 2002). Assuming an idealized sinusoidal forcing function, thermal skin depth is given by: Under Mars atmospheric surface pressure ranges, thermal conductivity (k) of particulate regolith is expected to vary by 3-4 orders of magnitude, whereas bulk density and specific heat vary by a factor of ∼2-3 or less, making variations in TI on Mars most sensitive to variations in k (Jakosky, 1986;Presley & Christensen, 1997a;Wechsler et al., 1972). Particle size determines solid and gas conductivity in loose sediment, while cementation, porosity, vesicularity, and/or fracturing controls conductivity in indurated sediment or rock (Presley & Christensen, 1997a, 1997b, 2009a, 2009bWechsler et al., 1972).
The relationship between sediment grain size and thermal conductivity under Martian pressures has been experimentally determined by Presley and Christensen (1997b), but quantitative relationships between TI and properties in rock under Martian atmospheric pressures are not as well constrained. Finally, orbitally derived TI values can be interpreted in multiple ways when using only a single local time, particularly for values >∼350 J m −2 K −1 s −1/2 , where it could indicate very coarse grained, homogeneous sediment or mixtures of finer-grained sediment and higher-TI materials (Piqueux & Christensen, 2011).
TI values for a planetary surface are derived through the use of thermal infrared observations and thermophysical models. Thermal models (e.g., Haberle & Jakosky, 1991;Kieffer, 2013;Mellon et al., 2000;Vasavada, Bandfield, et al., 2012, Vasavada, Chen, et al., 2012 are used to generate lookup tables of predicted surface temperatures at different times of the year, times of day, geographic locations, surface pressures, albedos, atmospheric dust opacity values, and atmospheric layering conditions for a range of TI values. The model temperatures are then compared with brightness temperature observations taken from orbit (Audouard et al., 2014;Fergason, Christensen, & Kieffer, 2006;Mellon et al., 2000;Putzig & Mellon, 2007a) or from a rover (Fergason, Christensen, Bell, et al., 2006;Edwards et al., 2018;Vasavada et al., 2017). Interpolation between modeled and measured temperatures is then used to determine the best-fit TI value for the surface. This process is described in further detail in Section 3.2.

Challenges to Interpretation of TI Values
Interpreting the TI of a surface with a single temperature measurement presents a significant challenge. A diurnal temperature curve for any surface has an asymmetrical bell shape that peaks in the early afternoon and flattens throughout the night. The shape and amplitude of the curve changes with TI and also with lateral or vertical mixtures of materials of differing TI (Figure 1). Mixtures can increase or decrease surface temperatures differently depending on the individual material TIs and the proportions and order of those materials present. For example, a rocky surface topped by a thin layer of dust will have a lower temperature in morning observations and a higher temperature in afternoon observations than that of an uncovered, homogenous rocky surface. Low TI materials, such as dust, cool much faster through the night and warm to a greater degree during the day than do high TI materials, such as rocks. On a diurnal temperature curve, a dust-covered rocky surface appears to be anomalously cool in mornings and warm in evenings as compared to a dust-free rocky surface (Figure 1a). In the cases of rover-derived TI values, it is possible to measure surface temperature fluctuations throughout an entire diurnal cycle and have more points of comparison on a temperature curve to those generated by models (e.g., Fergason, Christensen, Bell, et al., 2006;Vasavada et al., 2017;Edwards et al., 2018). However, many orbital-based TI studies are forced to rely on a limited AHERN ET AL.

10.1029/2020JE006713
3 of 30 Figure 1. Modeled diurnal temperature curves for comparison. (a) A homogenous surface of dust (TI = 100 J m −2 K −1 s −1/2 ; blue line) and exposed basaltic bedrock (TI = 1800 J m −2 K −1 s −1/2 ; brown line), and a 1 mm covering of dust over rock (orange line) at 0° latitude/longitude and Ls = 1 (∼northern spring/southern autumn). (b) A homogenous surface of intermediate thermal inertia (TI) (TI = 600 Jm −2 K −1 s −1/2 ; blue line) as compared with a bedrock surface (TI = 1200 Jm −2 K −1 s −1/2 ) overlain by 2 mm of dust (TI = 100 Jm −2 K −1 s −1/2 ; orange line) at 0° latitude/longitude and Ls = 270 (∼northern winter/southern summer). The gray boxes in both (a) and (b) indicate the local time windows of Mars Odyssey Thermal Emission Imaging System (THEMIS) observations used in this study. (b) Illustrates the relative difficulty of separating different thermophysical scenarios that produce similar temperatures in the pre-sunrise window as compared with the postsunset window. number of observations due to orbital schedule constraints or limited mission lifetimes, which can make it difficult to distinguish between different thermophysical scenarios ( Figure 1b).
Because specific temperatures at a single time of day/season can arise from multiple thermophysical scenarios (a homogeneous surface, a checkerboard-type mixture, or near-surface vertical layering), a TI value derived from a single temperature is considered an "apparent" TI (e.g., Audouard et al., 2014;Putzig & Mellon, 2007a, 2007b. (We note that "apparent TI" has a slightly different meaning in terrestrial thermal remote sensing studies, where it is defined as one minus albedo, divided by the difference between daytime and nighttime temperatures, e.g., Price, 1977). Using single apparent TI values has allowed for mapping of large-scale global units based on relative TI differences, as well as general characterization of the physical properties/grain size of those units (e.g., Fergason, Christensen, & Kieffer, 2006;Mellon et al., 2000;Putzig et al., 2005). However, for local scales, particularly where there is a significant portion of bedrock exposed, apparent TI values from a single time of day do not provide an accurate or complete view of the thermophysical properties of that surface (McCarty & Moersch, 2020). Therefore, using multiple overlapping temperature measurements taken at different times of day and times of year can be used to better distinguish between a range of possible thermophysical scenarios for Martian surfaces.
Several studies have seen evidence for thermophysical surface heterogeneity and some have attempted to characterize the nature of the heterogeneity, recognizing the importance of using temperatures from multiple times of the Martian diurnal cycle. "Anomalous afternoon cooling" of about 32% of the Martian surface was observed by the Viking IRTM instrument (Ditteon & Kieffer, 1980;Jakosky, 1979;Kieffer et al., 1976;Palluconi & Kieffer, 1981), which Ditteon (1982) attributed to surface layering (∼1-4 cm) of fine sediments over possible duricrusts. Mellon (2007a, 2007b) recognized temporal variations in apparent TI using both daytime (∼1400 local time) and nighttime (∼0200 local time) observations with TES and attributed those variations to near-surface vertical layering, with some lateral heterogeneity. They also created maps of surface heterogeneity at 5 deg/pixel based on seasonal changes in apparent TI. Audouard et al. (2014) used overlapping "end-of-dawn" and "beginning-of-dusk" temperature measurements ranging 4 MY using the Observatoire pour la Minéralogie, l'Eau, les Glaces et l'Activité (OMEGA) instrument onboard Mars Express (MEX) to show changes in apparent TI that can be attributed to heterogeneous surfaces. Their work mapped surfaces of interest at 20-32 ppd resolution, which represented a significant advance in resolution compared with TES TI studies. Mellon, (2007a, 2007b) and Audouard et al. (2014) reported that the majority of the Martian surface displays moderate TI values, which could often be attributed to mixed surfaces. In some of their local-case studies, these authors reported the presence of possible dust or dune cover or indurated top layers as reasons for changing apparent TI values over time. Simurda et al. (2019) noted that individual flows of Daedalia Planum, although compositionally similar, exhibited different diurnal heating patterns caused by vertical near-surface layering and/or horizontal mixing of lava, dust, and sand. Additionally, McCarty and Moersch (2020) showed that diurnal variations in orbitally derived apparent TI values are sufficient to detect lateral surface heterogeneity and have ground-truthed these results with rover imagery.
These studies have laid the foundation for our work, in which we use multiple overlapping THEMIS brightness temperature observations taken in the presunrise (0300-0600 local time) and postsunset (1800-2000 local time) windows that correspond with over 9 MY worth of data. Our goal is to better constrain the TI, and thus the physical properties, of geologic components within certain regions of interest on Mars. Specifically, for purposes of this study, we focus on two previously visited landing sites (Gusev and Gale craters) and one upcoming landing site (Jezero crater), where we have larger volumes of data and the opportunity to groundtruth our orbital-based methods. The techniques presented in this study could be applied to anywhere with clear diurnal temperature cycles and multiple spatially overlapping temperature observations. Information gained through these methods can provide insight into locations that lack ground truthing opportunities, while also providing useful tools in future mission planning and rover trafficability studies.

Temperature Derivation and Data Selection
Surface temperatures for this study were observed in two nighttime windows: presunrise (0300-0600) and postsunset , when surface temperature contrasts are highest in the diurnal cycle and uncertainty from albedo effects on temperature is lower. We considered using daytime observations in this study but found that albedo uncertainties at this time of day increased modeled temperature uncertainties by ∼0.2-0.5 K compared to nighttime observations (Section S2). The evolving local times of the Mars Odyssey mission have allowed for data collection in the 0300-0600/1500-1800 window earlier in the mission and the 0600-0800/1800-2000 windows later in the mission, where it is currently operating. In total, these data span 9 MY and have captured the surface in all seasons.
Surface brightness temperatures were estimated by fitting Planck curves to THEMIS band 3 (centered at 7.93 μm) or band 9 (centered at 12.57 μm) calibrated radiance-at-sensor data; no atmospheric correction was applied, and this likely contributes in part to the ∼2 K accuracy assessment for the THEMIS instrument that was suggested by Fergason, Christensen, and Kieffer (2006) based on examination of polar cap surfaces. We note that (a) of the available THEMIS bands, the Martian atmosphere is most transparent in the wavelength regions corresponding to bands 3 and 9 (Text S1), excluding cloudy periods, and (b) many nighttime images only have three bands available, greatly complicating an atmospheric correction. Among the surface-sensing THEMIS bands, the signal-to-noise ratio (SNR) is highest in band 9, which is centered closer to the peak of the planetary radiance curve. Thus band 9-derived temperatures were used in the pre-sunrise (0300-0600) window, when surface temperatures are colder (Fergason, Christensen, & Kieffer, 2006;Edwards et al., 2009). However, for the postsunset observations, we chose to use band 3 to avoid the strong CO 2 transmission feature that is present in band 9 during daytime and postsunset  observations, which can lead to overestimates of surface temperature. To further corroborate these band choices, we compared temperatures derived from TES radiance spectra with those derived from TES radiance spectra convolved to THEMIS resolution, from bands 3 and 9 (Text S1, Figure S1). Finally, in this study, we assumed an emissivity of 1.0 at the chosen band and thus brightness temperatures are considered to be kinetic temperatures. This assumption is valid for band 3, which overlaps with the Christiansen frequency of most silicates. For band 9, surface emissivity is likely to be slightly less than 1.0 (e.g., Fergason, Christensen, & Kieffer, 2006), but for presunrise THEMIS observations, this emissivity overestimate is effectively offset by the atmospheric CO 2 emission feature that is present (Text S1, Figure S1).
The regions of interest in this study are Gusev crater, Gale crater, and Jezero crater: past, present, and future landing sites. We selected these regions because of the relatively large amount of THEMIS IR coverage over these areas, as well as the ability to ground-truth our methods with rover-based observations. Within these regions of interest, we identified geologic units with distinct geomorphological or thermophysical characteristics. These units are identified in Table 1 and Figure 2, and are further discussed in Section 4. We studied each unit in THEMIS, Mars Reconnaissance Orbiter's (MRO) Context Camera (CTX), and MRO's High Resolution Imaging Science Experiment (HiRISE) imagery to determine the best surface of each unit for temperature sampling. Our criteria for each sampling location was: (a) an apparently horizontal surface to limit slope effects on temperature derivation, (b) a surface with the least amount of visible surface complexity (e.g., craters, bedforms), (c) the best representative surface of physical characteristics of the unit (i.e., spectral character, average sediment cover or rock exposure, geomorphological expression), and (d) the maximum surface area, while maintaining the other three criteria to reduce THEMIS temperature uncertainty (Fergason, Christensen, & Kieffer, 2006). The polygons from which we have sampled temperatures are found in Figure 2 and are further discussed in Section 4.
Once locations within each unit of interest were determined, we averaged the brightness temperatures from all pixels covered in a sample polygonal area for each available THEMIS image overlapping that polygon. By doing so, we obtained averaged temperatures from across the polygon from each local time and solar longitude (Ls) that a THEMIS observation was acquired. Images that contained artifacts, distortions or visual haziness were discarded to reduce uncertainty in temperature derivation. We also removed observations AHERN ET AL. that yielded temperatures <160 K due to the possible presence of CO 2 frosts and the inherent thermophysical complexities that frosts would introduce.

Thermal Model (KRC), Inputs, and Uncertainties
To investigate the possibility and nature of surface heterogeneity at these sites, we first used the THEMIS-derived temperature measurements to predict TI values, using the KRC model (Kieffer, 2013). The KRC model (k: conductivity, ρ: density, c: specific heat; the terms in the TI equation) has been widely used to study the surface of Mars from orbit and with rovers (Fergason, Christensen, Bell, et al., 2006, Fergason, Christensen, & Kieffer, 2006Edwards et al., 2009;Kieffer, 2013). Among its many features, this model can be used to predict planetary surface temperatures or surface apparent TI values at all latitudes, times of day, and seasons up to any depth using temperature-dependent thermal conductivity and specific heat. KRC assumes a one-layer, two-stream (solar and thermal wavelengths) atmosphere that is spectrally gray at solar wavelengths and accounts for condensation and global pressure variations. The surface and atmosphere are radiatively coupled, with thermal radiation assumed to be isotropic and spectrally gray. Surface and subsurface temperatures are calculated by KRC through the use of explicit forward finite-differences in solving the heat conduction equation, while taking into account surface upward emission, downwelling radiation, an insulating bottom layer, both direct and diffuse insolation, and the latent heat of CO 2 , if it reaches saturation temperatures (Kieffer, 2013 Table 3)/interpolated by model from surface T

INERTIA2
Layer 2 TI Set by user (see Table 3)

Thick
Thickness of layer 1 Set by user (see Table 3) Abbreviations: THEMIS, thermal emission imaging system; TI, thermal inertia.

Table 2 Input Parameters Used in This Study and Their Sources
(e.g., Audouard et al., 2014;Putzig & Mellon, 2007a, 2007bPutzig et al., 2005). The normal mode of forward-computing temperature for a set of input TIs, including vertical layering, was then further used to evaluate the nature of that heterogeneity. In one-point mode, a lookup table of surface temperatures as a function of TI values is calculated for a given set of inputs, after a 2-year modeled spin-up period.
In this work, we specified latitude, longitude, Ls, local time, visible dust opacity, visible albedo, and a broadband infrared/visible dust extinction ratio (see Table 2). Some inputs, such as elevation and associated surface atmospheric pressures, are set by the model (see Table 2), which defaults to values from the 2 ppd Mars Orbiter Laser Altimeter (MOLA) elevation map (D. E. Smith et al., 2001). Albedo and dust opacity inputs, and associated uncertainties, are described below.
Bolometric lambert albedo values (∼0.3-3.0 μm) are available from the Mars Global Surveyor TES instrument, however, the spatial resolution (∼3 × 8 km) is too coarse to resolve our regions of interest. Thus, surface albedo values were estimated from THEMIS VIS-ALB images using a previously determined empirical relationship between TES albedo and THEMIS VIS-ALB lambert albedo values, as reported in the THEMIS Geometric Processing User's Guide and given in Table 2 (Murray et al., 2020). We adopt an albedo uncertainty of 0.05, based on the reported "albedo sigma" threshold reported for THEMIS VIS-ALB processing by (Murray et al., 2020) (interpreted as the maximum deviation from a globally determined relationship between TES albedo and THEMIS-derived albedo). Our regions of study are expected to be on the higher end of the uncertainty range, given that they are relatively close to the equator, in these latitudes, the Sun is above the horizon for a longer time of day and gives less time for albedo effects to disperse through the night (Fergason, Christensen, & Kieffer, 2006). Finally, we note that we use a single albedo for each polygon; the assumption of an unchanging albedo likely leads to some incorrect modeled temperatures for some observations. This is unavoidable due to lack of corresponding albedo observation for each temperature observation, and is discussed in Section 5.  (Montabone et al., 2020). The climatology dust opacities are at infrared (IR) wavelengths; to convert to visible dust opacity (TAUD), we multiplied the IR values by 2.6, after Montabone et al. (2020). The ability to retrieve TAUD values from the same Ls, Mars Year, and geographic location as the THEMIS temperature data greatly reduce the uncertainty due to atmospheric dust. However, the dust opacity values of Montabone et al. (2020) have associated uncertainties that must be accounted for in our analysis. We obtained those uncertainties from the climatology database for each THEMIS observation used in this study, and found median values of 0.10 (absolute) for Jezero and Gale craters, and 0.08 for Gusev crater. In this study, we adopt an IR opacity uncertainty of 0.10.
We estimated the modeled temperature error associated with the albedo and dust opacity uncertainties described above. For one location at each of the three general study regions, modeled temperatures were calculated for the latitude, longitude, Ls, and local time for each THEMIS observation used in this study, for three albedo values (measured, and measured ± 0.05), and for three IR dust opacity values (measured, and measured ± 0.10), for a range of TI values (Figure 3). Because albedo and dust opacity uncertainties are uncorrelated, we calculated a total modeled temperature uncertainty by taking the root-mean-square value of the two sources of error, and then finally calculated the average root-mean-square value for all modeled TIs at each site. Thus, the modeled temperature errors associated with albedo and dust opacity are 1.28 K for Gale, 1.66 K for Gusev, and 1.52 K for Jezero. Finally, we add an observation uncertainty of 2K, based on the nighttime THEMIS temperature accuracy range of 1.8-2.8 K reported by Fergason, Christensen, and Kieffer (2006); because we are averaging several THEMIS pixels per polygon and are working in low latitude, warmer regions, we used an observation uncertainty of 2 K. Thus, the total model plus observation temperature error is estimated to be 3.28, 3.66, and 3.52 K for Gale, Gusev, and Jezero craters, respectively.

Initial Assessment of Surface Heterogeneity From Apparent TI
With apparent TI values for each temperature observatio,n we assessed each sampling site for potential surface heterogeneity. If surface mixing is present, apparent TI values appear to change with local time or season, whereas for a thermophysically homogenous surface TI values remain constant ( Figure 4). We used the diurnal apparent TI trends to make initial assessments and predictions about surface heterogeneity for each of the sampling polygons. Although TI is temperature-dependent, the apparent TI variations that we observe in this work are beyond expected diurnal temperature variation-induced changes (<∼30 J m −2 K −1 s −1/2 for dust-sized particles and up to ∼100 J m −2 K −1 s −1/2 for mm-sized grains, Figure 6 of Piqueux & Christensen, 2011). Therefore, our modeling assumes temperature-independent properties.

Modeling Surface Temperatures for Vertical and Lateral Mixing Scenarios
Examining the changes in TI values over days and seasons can indicate surface heterogeneity; however, the goal of this work is not only to identify the presence of heterogeneity, but also to narrow down the potential ranges and configurations of materials in these heterogeneous surfaces. To this end, we forward-modeled the temperatures for a range of two-component scenarios of subpixel lateral and vertical mixtures to represent the possibilities of lateral surface mixing and near-surface layering of materials with different thermophysical properties (Table 3)   represent the span of materials expected for the surface of Mars, those of sediment grains ranging from dust to gravel (Presley & Christensen, 1997b;Wentworth, 1922) and those of rocks spanning from the most weakly consolidated or highly porous to the most dense and competent. The values used in this study are listed in Table 3, along with their corresponding grain sizes and/or verified measured samples. Modeled temperatures were predicted for the same Ls and local times as the temperatures derived from the actual THEMIS observations.
For lateral binary surface mixtures, we started by modeling the surface temperatures from potential end-member components listed in Table 3, using KRC. We then converted those derived temperatures to blackbody curves using the Planck equation and mathematically mixed every combination of two blackbody curves in proportions of 25/75, 50/50, and 75/25 for all TI values to create new blackbody curves. This is possible because radiance curves of spatial mixtures of materials mix linearly (Audouard et al., 2014;Bandfield, 2009;Nowicki & Christensen, 2007). Finally, the new curves were convolved to THEMIS spectral band passes and used to determine surface brightness temperatures in the same manner, as was done for the THEMIS data (described in Section 3.1). These brightness temperatures represent the observed temperature when there is subpixel lateral mixing of two types of materials. The total number of lateral mixing scenarios used in this study is 474. Albedo for the individual components within laterally mixed surfaces was kept at one constant value for simplicity.

Table 3
The

Range of Thermophysical Values and Their Corresponding Physical Properties Used in This Study
The KRC model also has the ability to calculate brightness temperatures for two-layer scenarios, in which the user specifies the TI values for the upper and lower layers and the thickness of the top layer. We input combinations of all of the values in Table 3 as upper and lower layers, with the thickness of each top layer determined by the annual skin depth. The thicknesses used in each combination of upper and lower layers ranged from 1/512 to 3 full annual skin depths for each top layer with each step increasing by a factor of 2 ( Figure 5). With these combinations, the total number of vertically layered thermophysical scenarios considered in this study is 2076. After temperatures were retrieved from the 2076 total vertical layering scenarios, we then discarded scenarios in which increasing the thickness of the top layer no longer changed the predicted temperatures. This can happen short of a full annual skin depth because even with the increased temporal coverage of THEMIS IR observations we are unable to measure over full annual cycles. Therefore, we cannot always thermally "see" through to the full annual skin depth. In these cases of redundant temperature measurements, we down-selected the vertical models to include the maximum thickness at which temperatures change and then removed those that gave the exact same temperature measurements with thicker top layers. This conservative approach to determining the maximum depth of sensitivity within each set of observations results in the total number of vertical layering models being different for each sampling polygon, depending on the number and temporal scope of available THEMIS observations. AHERN ET AL.
10.1029/2020JE006713 11 of 30 For completeness, we also considered scenarios in which the surface is laterally and vertically homogenous at each of the TI values used in our mixing and layering scenarios described above. A surface would appear thermophysically homogeneous if there was no lateral surface mixing present and if the top layer of material was thick enough to exceed the annual skin depth.

Evaluating Model Fits and Assessment of Nonunique Solutions
For each sampling surface, we ran the KRC model with all of the lateral surface mixtures, vertical layering, and thermophysically homogenous scenarios as described in Section 3.2.2 to obtain predicted surface temperatures at the same Ls, local time, and MY as our orbitally observed temperatures (see Table 4, Tables S1, S2, S3, S4). The best fit models for a given polygon were selected based on the maximum number of THEMIS observations fit by model within the total temperature uncertainty, which ranges from 3 K to 4 K for the polygons in this study (similar to error is given in Fergason, Christensen, & Kieffer, 2006, Section 3.2). In most cases, no models could match all of the THEMIS observations; potential reasons for this are discussed in Section 5. A root mean square (RMS) value was also calculated for each model based on differences between measured and modeled temperatures. The number of observations fitting within error and the best-fit models are given throughout Section 4 in Tables 5-9.
Nonunique solutions are possible and expected. In some cases, the scenarios that cluster in the lowest RMS groupings may mean roughly the same thing geologically, with minor differences in layer thicknesses or sediment grain size, for example. There are other instances, however, where there are geologically different solutions that could account for the observed temperatures. This is discussed in detail for each polygon in Section 4.

Rationale and Background
Gale crater (d∼154 km) is the result of a late Noachian/early Hesperian impact that occurred along the dichotomy boundary between the northern lowlands and the southern highlands. It is a partially filled crater with a ∼5 km layered mound at its center, known as Mount Sharp or Aeolis Mons (Thomson et al., 2011). The Curiosity rover landed in Aeolis Palus, a plain along the northern foothills of Mount Sharp. Its traverse up the flanks of Mount Sharp has taken it through basaltic conglomerates, cross-bedded sandstones, mudstones, and dune fields . The physical and chemical nature of these rocks and sediments, along with fans and deltas present within the crater, suggest that Gale crater housed a lake shortly after its formation (Grotzinger et al., , 2015. Since that time aeolian processes have scoured the floor of the crater. We selected two regions of study for this work ( Figure 6). The first is from the Bagnold dune field, which has predawn apparent TI values of 250-315 J m −2 K −1 s −1/2 from THEMIS (Fergason et al., 2012). The Curiosity rover traversed the Bagnold dunes in a narrow portion of the dune field that is several km southwest of our study area. Our study area choice was driven by the need to minimize complexity, thus we selected an area that appeared to consist entirely of sand, with no bedrock exposure, at CTX scale. The Gale observations are included here as an example; Gusev and Jezero observations are shown in supplementary Tables S1-S4. THEMIS, thermal emission imaging system.

Table 4 Thermal Emission Imaging System (THEMIS) Observational Data Retrieved From the Gale Crater Sampling polygons
The second study region lies on the "bedded fractured unit," as defined by Grotzinger et al. (2014), which exhibits the highest pre-dawn apparent TI values in the area of the landing ellipse, in THEMIS imagery, ranging from ∼450 to 555 J m −2 K −1 s −1/2 (Fergason et al., 2012;Grotzinger et al., 2014). The southwestern margin of this unit was traversed by the Curiosity rover in a region informally named Yellowknife Bay.
Ground-based imagery of the Yellowknife Bay region revealed finely laminated mudstones with patches of dark-toned unconsolidated material (Figure 5b). Vasavada et al. (2017) presented ground-based TI estimates of this unit derived from the Rover Environmental Monitoring System (REMS) Ground Temperature Sensor (GTS) diurnal temperature measurements, using both "simple models" (homogeneous surfaces) as well as "complex models" (lateral/vertical mixing scenarios). Using a simple model, the best fit TIs ranged between ∼370 and 540 J m −2 K −1 s −1/2 . However, complex models of a cm-thick high TI material (∼1700 J m −2 K −1 s −1/2 ), overlaying a low TI surface of 160 J m −2 K −1 s −1/2 , plus subpixel lateral mixing of fines (TI ∼160 J m −2 K −1 s −1/2 ), were found to be slightly better fits to these data .

Gale Crater Results (Figure 6, Tables 4 and 5)
Best fit model results of the Bagnold Dunes are most consistent with either thick fine-to-medium sands overlaying a high-TI material (1400-2000 J m −2 K −1 s −1/2 ) or fine-to-medium sands laterally mixed with more moderate TI materials (265-600 Jm −2 K −1 s −1/2) The latter seems more likely based on previous observations and studies (Presley & Christensen, 1997b;Ehlmann et al., 2017;Ewing et al., 2017;Lapotre et al., 2017;Edwards et al., 2018), in which medium-sized sands were reported. The sampling area selected visually appears to be mostly homogeneous sand (Figures 6c, 6d, 6f, and 6g), but the higher-TI components selected in some of the best-fit models could indicate patchy armoring of fine-sand dunes by coarser clasts, as also reported in the study by Edwards   The rocky surface is best modeled by a layering of ∼1.5 mm of fine sand (182 J m −2 K −1 s −1/2 , 100 μm) with a moderate TI unit (500 J m −2 K −1 s −1/2 ) below. The moderate component is consistent with expected values of poorly compacted sandstones or mudstones, volcaniclastic rocks, or moderately fractured rocks, via comparison with Mini-TES-derived TI values from clastic rocks in the Columbia Hills (600 J m −2 K −1 s −1/2 , Fergason, Christensen, Bell, et al., 2006). The moderate component is consistent with GTS-derived TI values using a homogeneous-surface model , see Section 4.1.1), but not with the complex best fit model described by Vasavada et al. (2017). The thermophysical differences between the dunes and adjacent rocks are visible both in diurnal apparent trends and in the best fit thermophysical scenarios selected.

Rationale and Background
The Spirit rover traversed the Gusev volcanic plains from its landing until it reached the base of the Columbia Hills at sol 156 of the mission. Along the way, Spirit confirmed variability in rock size and abundance as had been predicted by Golombek et al. (2003) from THEMIS observations. Both rock size and rock abundance in the plains along the traverse seem to be related to impact ejecta from nearby craters larger than 50 m in diameter, with larger and more frequent clasts occurring closer to the their respective craters Greeley & Batson, 1990;Maki et al., 2003). The plains are olivine-rich basaltic materials that have been reworked by impacts and surface deflation and have experienced minor chemical alteration. Pancam survey images of the plains surface reveal sediment that is poorly sorted, made of subangular and subrounded grains indicative of emplacement as ejecta, rather than by fluvial or aeolian means. The surface rocks observed were either small (a few cm) and vesicular or larger and not vesicular (Crumpler et al., 2005;Grant et al., 2004). The nonvesicular rock clasts have high TI values (∼1200 J m −2 K −1 s −1/2 ), as measured from Spirit's Mini-TES instrument, suggesting emplacement as effusive basaltic flows (Fergason, Christensen, Bell, et al., 2006). A thin dust coating is observed on most rock and plains surfaces (Crumpler et al., 2005;Grant et al., 2004). The sampling area we selected in this study is located along the traverse at Spirit's location at sols 113-122. We selected this location based on the proximity to the traverse, while avoiding the craters (Figure 7a-7c and 7e).
The contact between the volcanic plains and the Columbia Hills is well defined from orbit and as observed by Spirit. Older Columbia Hills rocks were observed to be more friable and more chemically and physically weathered than the younger volcanic plains that embay them . Indeed, Mini-TES measurements indicate that Columbia Hills rocks have lower TI (∼620 J m −2 K −1 s −1/2 ) than some of the individual clasts found within the volcanic plains (Fergason, Christensen, Bell, et al., 2006), indicating a more porous, vesicular, or weakly consolidated material. However, in THEMIS nighttime IR imagery, the Columbia Hills seem relatively warm, indicating higher TI. We sampled a region with relatively high nighttime radiance, ∼1.5 km south of Home Plate (Figures 7a, 7d, and 7f).
The differences in physical properties as measured and observed in imagery between the volcanic plains and the Columbia Hills suggest that they have different geologic origins. Fergason, Christensen, Bell, et al. (2006) interpret the volcanic plains as effusive, dense basaltic flows and the Columbia Hills rock as volcaniclastic sediment or ash flow deposits. We selected these two distinct units as geologic units of interest at Gusev crater to evaluate the aforementioned depositional hypotheses and shed further light on processes that have further modified these units (Figure 7; Table 6). A ground-based image from the Spirit navigation camera, acquired on Sol 113, is shown in Figure 7b. The Columbia Hills region selected was not directly traversed by the Spirit rover because we wanted to characterize the region with the largest area of likely rock exposure, indicated by the higher nighttime radiance.
AHERN ET AL.

10.1029/2020JE006713
16 of 30  Table 6 (the x-axes in these plots are arranged by the local time of the THEMIS observations). (j) Provides a visual representation of the best fit models plotted in (h-i) and listed in Table 6 Table 5 (the x-axes in these plots are arranged by the local time of the THEMIS observations). (j) Provides a visual representation of the best fit models plotted in (h-i) and listed in  (Figure 7, Table S1, Table 6)

Gusev Crater Results
The plains unit sampling area (PU) is most consistent with moderate TI values for both the lower and higher TI materials (Table 6). The best fit models indicate a dominance of medium to coarse sediment with possibly larger clasts mixed in (265-500 J m −2 K −1 s −1/2 ). This is consistent with the sediment and clasts seen strewn along the plains as observed by the Spirit rover. Additionally, Spirit encountered slightly indurated sands on the order of ∼6-10 cm depths as it trenched the plains unit, so the higher TI values between 500 and 700 J m −2 K −1 s −1/2 could also be reflecting the slight increase in cementation of the sands (Sullivan et al., 2008).
The model fits at CH predict moderate TI values overlying low TI values ( Table 6). The moderate values of the upper layer could indicate some induration or a block contribution. However, CH values are also consistent with very porous or weakly consolidated rock. The hypotheses by Fergason, Christensen, Bell, et al. (2006) of this unit consisting of volcaniclastic sediment or being emplaced as an ash flow would be consistent with our observed TI values of 400-1000 J m −2 K −1 s −1/2 for the higher-TI component. It is possible that mechanical breakdown of the uppermost sections of the rocks here and surface-atmosphere exchange has allowed for patches of induration to occur. HiRISE imagery of this surface shows clear evidence for lithified material protruding through unconsolidated surface materials, which may be reflected by the higher TI components

Rationale and Background
The volcanic floor unit (VFU), as mapped by Goudge et al. (2015), is considered to be the youngest and uppermost stratigraphic unit within Jezero crater. Goudge et al. (2015) interpreted it to be effusive volcanic in origin because of its moderate-toned and crater-retaining appearance and because it exhibits infrared spectral properties consistent with unaltered mafic material. VFU also appears to be relatively thin as compared with other units in the crater (<10 m thick at its margins) (Goudge et al. 2015). Its formation age is estimated to be no earlier than ∼3.45 Ga and as late as ∼1.4 Ga, when the light toned floor unit and the fan deposits were already in place ( Note: Thicknesses given in mm.

Table 6
Best Fit Model Results for the Gusev Crater Sampling Polygons been mapped as one distinct layer (Goudge et al., 2015), it has notable variation in its surface expression across the crater; closer to the western fan, the volcanic unit appears smoother than in the east/southeast, where it looks rougher (Warner et al., 2020. Note that Warner et al. refer to this unit as the mafic floor unit, or MFU). The VFU differs from other Martian lava plains in that its ∼10 m scale craters have rocky ejecta/ rims and low spatial density (Warner et al., 2020). Combined, these observations suggest a surface that lacks a thick regolith, and/or has been stripped of the fines portion of the regolith, where the rough terrain represents a surface in a more advanced state of deflation or exhumation (Warner et al., 2020). We sampled three locations within the volcanic floor unit to further evaluate its properties and thermophysical variation (Table 1, Table S2, Table 7). VFU1 samples the smoother section proximal to the fan (Figure 8 a,b), whereas VFU2 and VFU3 sample two locations progressively further from the fan (Figures 8a, 8c, and 8d). Visually, VFU3 appears to have the greatest bedrock exposure, whereas VFU1 appears mantled.
The second geologic unit of interest selected for this study in Jezero was the light toned floor unit (hereafter, LTU) as mapped by Goudge et al. (2015), which lies stratigraphically below the volcanic floor unit and the fan deposits. Bedforms are prevalent across LTU; and in places where it is exposed below the dunes or in erosional windows below VFU, it appears to be heavily fractured (Figures 8a-8g). The mineralogy of the light toned unit is difficult to ascertain given the frequent cover of olivine-rich dunes. However, in select locations where it is better exposed, Mg-rich carbonates have been identified (Goudge et al., 2015;Horgan et al., 2020;Mustard et al., 2009). The physical nature of the light toned unit, although often obscured by dunes, is important to understanding the early sedimentary processes that formed and modified it over time. Our two sampling areas, LTU1 and LTU2 (Table S3), show the variation of sediment cover over the fractured unit. LTU1 samples a surface with more of the light toned rock exposed, whereas LTU2 covers a section with more sediment cover (Figure 9a-9e).
The third and final Jezero unit examined in this study was the western fan deposit (WFD), also mapped by Goudge et al. (2015). The northern and western fans have been interpreted as delta deposits that resulted from fluvial materials being transported and deposited into the Jezero paleolake (Ehlmann et al., 2008;Fassett & Head, 2005;Schon et al., 2012). Thin sedimentary layers and preserved channel meanders are visible in the western fan from orbit (Fassett & Head, 2005;Schon et al., 2012), and spectral studies have indicated that it is likely composed of olivines, Mg-rich carbonates, and Fe/Mg smectites (Goudge et al., 2015). The thermophysical properties of the western fan deposit may provide insight into the types of rocks available within the watershed that contributed to its formation, as well as the erosional byproducts that either supply sediment throughout the crater or are no longer present. We have selected two sampling polygons within the western fan (Table S4). (Figure 8, Table S2)

Volcanic Floor Unit Results and Discussion
The best fit models for the three volcanic floor unit locations are summarized in Table 7. VFU1 temperatures are consistent with having a dust to fine sand component (∼5-100 μm; 91-182 J m −2 K −1 s −1/2 ) and another low-to-moderate component ranging from 251-400 J m −2 K −1 s −1/2 (Figure 8n). The values of the second component are consistent with medium to coarse sand (400 μm-2 mm, respectively; Presley & Christensen, 1997b). Best fit models of VFU2 are largely the same as those for VFU1. The low-to-moderate TI component in these two polygons could be consistent with indurated material, very coarse sediment (mm to cm-sized clasts), or porous/fractured rock. Best fit models of VFU3 suggest either a thin covering of dust over coarse sands (365 J m −2 K −1 s −1/2 ), or lateral mixing of silt (60 μm; 162 J m −2 K −1 s −1/2 ) and a moderate-to-high TI material (1000 J m −2 K −1 s −1/2 ), where the higher TI component could represent very coarse sediment (e.g., cm-sized clasts), indurated sediment, or rock ( Figure 8n).
The apparent TI values at VFU3 are slightly higher than those at VFU1 and VFU2. This is also seen in the possible fit to a lateral mixture with a 1000 J m −2 K −1 s −1/2 component. Orbital imagery seems to show greater bedrock exposure at VFU3. There could be a mechanism eroding or reworking the VFU more efficiently closer to the fan and causing it to break down, thereby lowering the TI in this area relative to the center of the crater. Warner et al. (2020) suggest that the smooth portion of the VFU may have been less stripped of fine sediments by wind deflation compared to the rough portions, or that the material covering VFU1 and VFU2 has been derived from erosion of the adjacent WFD materials. Ruff (2017) suggested the possibility that the WFD unit could be stratigraphically younger than the VFU, and that the deposits may have been    Table 7 (the x-axes in these plots are arranged by the local time of the THEMIS observations). (n) Provides a visual representation of the best fit models plotted in (k-m) and listed in partially eroded away over time. Thus the smooth, lower TI materials could represent a lag deposit from eroded delta materials. Alternatively, the VFU itself may have spatial physical variations based on lava chemistry, lava gas content, emplacement stage, or its flow behavior as related to rocks below. (Figure 9, Table S2, Table 8)

Light Toned Unit Results and Discussion
The best model fits for LTU1 indicate either vertical layering involving low-to-moderate TI components, or lateral surface mixtures of low TI, fine grained sediments, and a higher TI component (500-1200 J m −2 K −1 s −1/2 ). A value of 365-500 J m −2 K −1 s −1/2 corresponds to very coarse sand (2 mm), but could be representative of the highly fractured rock that has been observed in orbital imagery. When compared with a nonporous, competent bedrock, highly fractured rocks have much lower TI values due to the lessened heat storage ability and increased surface area for heat to escape. The corresponding imagery (Figure 9d) shows both bedforms and lithified material exposed at the surface. The best fit models seem to be capturing either the bedforms and the rock, or the bedforms and a light dust layer. It is possible that a three component scenario actually exists, where exposures of rock between the bedforms are also covered in a thin layer of dust.
The best fit models at LTU2 are consistent with subpixel lateral mixing scenarios, with differing proportions of higher and lower TI materials based on the pairings of values (Figures 9j, Table 8). Overall, at the light toned unit sampling sites, we observed mixtures of what appear to be fine sediments with low to moderate TI materials that could be identified as larger-grained sediments or low TI rocks. We note that, with the exception of one model fit at LTU1, no models indicate a well-cemented, low-porosity rock at either of the LTU sites; however, a well-cemented rock TI would be lowered by fracturing, and additionally, our models may be impacted by the surface complexity mentioned above. (Figure 10, Table S2, Table 9)

Western Fan Deposit Results and Discussion
Both the WFD1 and WFD2 polygons yield nearly identical results geologically. Best fit models predict a lateral surface mixture of 75% 100 μm (182 J m −2 K −1 s −1/2 ) sands and 25% 700-800 J m −2 K −1 s −1/2 . The mantled appearance of the unit supports the larger proportion of sand, likely derived from the underlying intact fan material. The values of 700-800 J m −2 K −1 s −1/2 are consistent with rocks of sedimentary origin that are lightly or moderately cemented, as would be expected for rocks of delatic origin.

Discussion and Conclusions
The previous sections detailed the results and discussion for each region of study. Here, we summarize some of the overall challenges with the available data and methods for interpreting surface heterogeneity, and conclude by listing the major results for our regions of study.
One of the universal findings from this work is that there were few cases in which all of the THEMIS observations could be fit with a single lateral, vertical or homogeneous modeled scenario, despite having over a thousand potential combinations to choose from. Though the surface is certainly more complex than can be captured in our two-endmember scenarios, the magnitude of measured-modeled temperature differences for some of our THEMIS observations require explanations beyond just surface complexity. We suggest the following possible reasons for this, including: use of dust opacity values that are not representative of the local conditions at the time of the THEMIS observation, local-scale and/or daily variability in near-surface atmospheric or ground conditions that are unaccounted for in the model (e.g., winds, haze, H 2 O frosts) (e.g., Fergason, Christensen, & Kieffer, 2006;Rangarajan & Ghosh, 2020), seasonal variability in surface warming from equatorial clouds, and real surface changes (e.g., dust deposition and removal) that have occurred during the lifetime of the THEMIS mission. We describe these further below.
For dust opacity, we utilized the best-available dust opacity information as a function of latitude/longitude (to the nearest degree), and Ls (Montabone et al., 2015(Montabone et al., , 2020. However, it is conceivable that local conditions could occasionally differ from the gridded dust opacity values, both spatially and temporally.  Table 8 (the x-axes in these plots are arranged by the local time of the THEMIS observations). (j) Provides a visual representation of the best fit models plotted in (h-i) and listed in Similarly, temperature variations due to daily weather conditions (Fergason, Christensen, & Kieffer, 2006), and possibly to frost formation (Martínez et al., 2016;Rangarajan & Ghosh, 2020) are expected and could result in inaccurate modeled temperatures. Next, tropical clouds can result in nighttime surface warming by preventing efficient radiative cooling, during the spring and summer. Though this effect is most prominent in the low TI regions of Tharsis and Arabia (Wilson and Guzewich, 2014), it could have a potential effect on our measured temperatures (up to ∼5 K) during these time periods. Such strongly affected observations would be discarded in our process of downselecting to observations that fit models within ∼3-4 K. Finally, surface dust deposition and removal could have occurred multiple times over the ∼9 MY of THEMIS investigation; this has been observed both from orbit and ground and can be either seasonal or long-term (e.g., Edgett & Newsom, 2018;Fenton et al., 2007;Rice et al., 2018;Robbins, 2020). This can lead to changes in albedo relative to our assumed single albedo values, but also represents a change in surface component TI, if dust abundances are changing significantly enough (e.g., tens of microns, Edwards et al., 2011). Because THEMIS predawn and postsunset temperature observations and THEMIS VIS images are not usually acquired within the same time periods, it is not possible to account for changing albedo or dust deposition in our models. Some of the problems described above could be alleviated by simultaneous dust opacity measurements with surface temperature measurements, knowledge of near-surface wind activity, and better long-term constraints on albedo changes at the scale of our observations.
Another issue is the problem of nonunique solutions, which would likely persist even without the above-mentioned sources of error. A major advantage of this work is the range of temporal coverage that THEMIS has provided during the course of its extended missions. However, THEMIS observations do not cover an entire annual temperature cycle due to orbital time constraints or image quality issues caused by weather and dusty time periods. A number of available images in this work had to be discarded (see Table 1 Table 9 (the x-axes in these plots are arranged by the local time of the THEMIS observations). (j) Provides a visual representation of the best fit models plotted in (h-i) and listed in Table 9, with gray boxes representing the higher TI component and orange boxes representing the lower TI component. Values in (j) are TI values in J m −2 K −1 s −1/2 . months, CO 2 frosts can accumulate on the surface or in the near subsurface; therefore, when we had images with temperatures below 160 K, we discarded those images to minimize surficial complexity. Ideally, temperature observations would span the diurnal cycle for every Ls over an entire MY, which would allow for better interpolation of TI from more complete temperature curves; but this type of coverage is not currently available.
Despite the problem of nonunique results, we are able to rule out many thermophysical scenarios and in particular, narrow down the TI of any bedrock component that is exposed. Based on our analyses of multiple overlapping THEMIS observations from 0300 to 0600 and from 1800 to 2000 local times, we can conclude the following: 1. The best fit models for the Bagnold Dunes ROI in Gale crater indicate a surface of medium-sized sands with patches of higher TI, possibly due to small induration surfaces or patchy armoring by coarser clasts 2. The bedded fractured unit  ROI in Gale crater is best modeled as a vertical layering of fine sands and a 500 J m −2 K −1 s −1/2 component below. This value is consistent with weak or poorly compacted clastic rocks, based on comparisons with Mini-TES observations at Gusev crater (Fergason, Christensen, & Kieffer, 2006;Thomson et al., 2013) 3. The best fit models for the Gusev plains unit (PU) indicate a dominance of medium-to-coarse sediments.
Some models indicated the presence of a low-to-moderate-TI component, which could indicate patchy induration, block contribution, or weakly consolidated rock. The ubiquitous thin dust layer observed by Spirit cameras is not reflected in the best-fit models, suggesting it is thermally insignificant 4. The Columbia Hills polygon (CH) polygon results are consistent with a weak rock or indurated layer overlying fine sands. This is consistent with the assessment by Fergason, Christensen, and Kieffer (2006) in which they proposed a surface composed of volcaniclastic rocks 5. The Jezero crater volcanic floor unit (VFU) best fit models suggest a competent rock component (∼1000 J m −2 K −1 s −1/2 ) mixed with fine sands in our VFU3 polygon, far from the WFD, or a thin layer of dust over coarse sediment/very weak rock (∼365 J m −2 K −1 s −1/2 ). In the smooth portion of the VFU, closer to the WFD, our models are consistent with mixtures of dust or fine sands with moderate TI materials (moderate to coarse sediments or weak rocks) 6. The two ROIs covering the light toned unit (LTU) in Jezero crater indicate bedform surface cover and possibly highly fractured rocks. The best fit models yield low to moderate TI values for both components at both of these polygons; these values are consistent with both sand and highly fractured/very weak rock 7. The polygons over the western fan deposit (WFD) have temperatures consistent with fine sands laterally mixed with low to moderately cemented sedimentary rocks The longevity and changing orbital schedule of the THEMIS mission has provided unprecedented temporal coverage. This has been vital in our efforts to determine individual unit TI values and physical properties. The workflow developed and discussed in this study can be applied to any surface on Mars that has multiple THEMIS IR observations for a more complete local thermophysical description of that surface. Detailed knowledge of thermophysical properties will provide greater understanding of the diversity of bedrock-forming and subsequent modification processes that have occurred in Mars' past. This can aid in both reconstruction of the geologic history of Mars as well as in providing detailed physical characterization for future Mars surface rover and human-based exploration.
AHERN ET AL.

Data Availability Statement
The full set and best-fit sets of modeled temperatures for each polygon, polygon shapefiles, and modeled apparent TI values are included in a data repository (Ahern et al., 2021;Zenodo). We note that an earlier version of the data exists in which there was an error in the TAUD values for the sampling polygons (Ahern et al., 2020; Zenodo). The KRC model (Kieffer, 2013) (Montabone et al., 2015(Montabone et al., , 2020 can be accessed from http://www-mars.lmd.jussieu.fr/mars/dust_climatology/.