The Effects of Terrain Properties Upon the Small Crater Population Distribution at Giordano Bruno: Implications for Lunar Chronology.

The distribution of impact craters on the ejecta of Giordano Bruno, a recent (<10 Ma) 22-km diameter crater within the lunar highlands, exhibits substantial variations. We surveyed craters D ≥ 10 m across a 1,323 km2 area of Giordano Bruno's ejecta and compared the distribution of craters with variations in thermophysical properties derived from the Lunar Reconnaissance Orbiter Diviner instrument. We used Diviner-derived rock abundance and nighttime regolith temperatures along with thermal model-predicted surface temperatures for a diversity of terrains to identify and isolate areas of the ejecta based on thermophysical properties such as bulk density and thermal conductivity. We found that thermophysical properties of the ejecta vary considerably both laterally and vertically, and consistently differ from typical regolith, indicating the presence of higher thermal inertia materials. Crater-size frequencies are significantly lower in areas with terrain properties exhibiting higher: rock abundance, nighttime temperatures, and/or modeled thermal inertia. This discrepancy in crater distribution increases for craters smaller than ∼25 m. These thermophysical variations indicate changes in the mechanical properties of the target materials. We suggest that these variations-specifically, terrain-dependent crater scaling variations and impactor-scale heterogeneities in material properties such as the presence or absence of large boulders-may influence crater diameters or inhibit crater production altogether in Giordano Bruno's ejecta; furthermore, these factors are size-dependent.

Introduction

Inconsistencies in crater distributions in and around lunar craters (e.g., different superposed crater distributions observed in counts on ejecta vs. counts on crater floors or melt ponds for the same large crater; and additionally, spatial variations within ejecta counts) were observed in early work (Hartmann, 1968; Strom & Fielder, 1968, 1970) and more recently by investigations using images from the Lunar Reconnaissance Orbiter Camera (LROC) (e.g., Ashley et al., 2012; Hiesinger et al., 2012; Plescia & Robinson, 2019; van der Bogert et al., 2017, 2010; Zanetti et al., 2017). These differences manifest as differences in cumulative numbers and power law slopes of crater size‐frequency distributions (CSFDs). It has long been known that target properties influence crater scaling, particularly at smaller crater sizes in the strength scaling regime (e.g., Holsapple, 1993; Melosh, 1989; Schultz & Spencer, 1979) and Shoemaker et al. (1968) suggested that late‐arriving fragments ejected at high angles create self‐secondary craters—both of which may contribute to the observed discrepancies in crater populations. This is of particular concern as crater counts conducted on the proximal ejecta of lunar craters are used to anchor the crater chronology of the inner solar system (e.g., Hiesinger et al., 2012; Robbins, 2014; Stöffler & Ryder, 2001).

Giordano Bruno (GB), a Copernican‐age crater within the lunar highlands (36°N, 103°E), is one of the youngest sizable impact craters (D = 22 km) observed on the Moon and exhibits substantial spatial variations in crater frequency and slope of the CSFD of craters superposed on its ejecta (Plescia & Robinson, 2019; Williams, Paige, et al., 2014; Williams et al., 2016, 2018). Recent work by Plescia and Robinson (2019) concluded that a substantial fraction of these craters were likely to be self‐secondary craters. This is supported by the observation of partial burial of some craters by impact melt, indicating the craters formed prior to the final emplacement of the melt (Plescia & Robinson, 2019; Williams et al., 2018; Xiao, 2018). Brightness temperatures measured using the Lunar Reconnaissance Orbiter (LRO) Diviner Lunar Radiometer Experiment (Paige, Foote, et al., 2010) reveal substantial spatial heterogeneity in thermophysical properties of GB's ejecta that show a correlation with crater densities, suggesting that variations in terrain properties have also likely exerted an influence on the production of craters and the resulting CSFDs (Williams et al., 2016, 2021).

Because of GB's young age of ∼1–10 Ma (Basilevsky & Head, 2012; Morota et al., 2009; Shkuratov et al., 2012), relatively fresh morphology, and minimally disrupted ejecta materials, it provides an optimal location to explore the influence of impact target effects on CSFDs. In this study, we systematically map the crater population superposed on the proximal ejecta of GB (Figure 1) and compare the distribution of craters with Diviner‐derived thermophysical properties. In addition to using published rock abundance and rock‐free regolith nighttime temperatures maps (Bandfield et al., 2011, 2017), we develop a one‐dimensional heat transfer model capable of modeling temperatures resulting from mixtures of regolith and higher thermal inertia materials. Nighttime temperatures are controlled by the thermophysical properties of the materials. By fitting the nighttime surface temperatures observed by the Diviner, we estimate vertical and lateral variations of the thermophysical properties within the diurnal skin depth of the regolith (upper ∼10 cm), which relate to variation in mechanical properties (e.g., Grott et al., 2019). From this we characterize how crater populations vary with terrain properties.

Figure 1

(a) LROC NAC images (Table 1) of Giordano Bruno crater. (b) Crater count survey boundary (white) with crater diameters marked (red). (c) Crater point density heat map: color scale ranges from 0 to 60 craters/km2. (d) Rock abundance (Bandfield et al., 2011, 2017) of the crater count area: color scale ranges from 0.0 to 0.2 surface rock fraction with contours at 0.05 and 0.10. (e) Regolith nighttime temperatures of the crater count area normalized for local time and latitude (Bandfield et al., 2011, 2017): color scale ranges from ‐20–60 K with contours at 10/20/30/40 K.

Table 1

LROC NAC Images Used in the Study in Order of Superposition as Presented in Figure 1 (Top to Bottom)

Product ID	Incidence angle (°)	Sub solar azimuth (°)	Resolution (m pix−1)
M1153560774L	64.76	162.91	1.25
M1168832870R	65.40	157.52	1.19
M1153560774 R	64.99	163.19	1.25
M1214775929L	67.55	165.91	1.36
M1214775929R	67.76	166.10	1.35
M1288827855L	58.19	152.17	1.20
M1288827855R	58.35	149.30	1.20
M1122929850L	65.12	163.38	1.51
M1122929850R	65.35	163.65	1.51
M1122936950L	65.03	163.69	1.50
M1122936950R	65.26	163.95	1.50
M1168832870L	65.20	157.06	1.19
M103838997R	65.70	198.54	1.53
M103831840R	65.13	200.62	1.53

(NASA/GSFC/Arizona State University).

Data and Methods

Crater counting was conducted using LROC Narrow Angle Camera (NAC) images (Robinson et al., 2010). Images were calibrated and map‐projected using the Integrated Software for Imagers and Spectrometers (ISIS) software (Anderson et al., 2004) and constrained to eastward‐only illumination with incidence angles 56°–75° optimal for crater measurements (e.g., Antonenko et al., 2013; Ostrach et al., 2011) with the exception of two images with westward illumination to fill in a 6.8 km2 area within the count area (Figure 1a; Table 1). Prior to map projecting the images, a bundle adjustment was performed to reduce image seam offsets (e.g., Klem et al., 2014). Counts were conducted using the Cratertools plugin for Arcmap (Kneissl et al., 2011). Crater diameters were measured down to 10 m in a region around the rim of Giordano Bruno representing a total area of 1,323 km2 (Figure 1b). CSFDs were fit with modeled crater‐age isochrons of Neukum et al. (2001) using the Craterstats2 program (Michael & Neukum, 2010) with absolute model ages (AMAs) derived using the Poisson timing analysis as described in Michael et al. (2016).

The crater distributions were compared to the thermophysical properties of the ejecta derived from Diviner observations. The crater counts were converted to a point density heat map for the count area using a circular search area with a radius of 1 km (Figure 1c). The crater heat map was compared to the 128 pixels per degree rock abundance and rock‐free nighttime regolith temperature maps available at the PDS (https://pds-geosciences.wustl.edu/missions/lro/diviner.htm) (Bandfield et al., 2011, 2017). The rock abundance and rock‐free regolith temperatures are derived by fitting anisothermality in the gridded nighttime observations of Diviner channels 6–8 (Paige et al., 2011) with modeled mixtures of rocks and rock‐free regolith where rock temperatures are estimated using a one‐dimensional thermal model (Bandfield et al., 2011). Anisothermality in the Diviner channels occurs as a result of the nonlinear nature of Planck radiance with respect to wavelength. Surfaces radiating at warmer temperatures within the instrument's field of view have an increased proportional influence on brightness temperatures at shorter wavelengths resulting in differing brightness temperatures in the Diviner channels when the surface footprint contains sub‐pixel variations in temperature (e.g., Williams et al., 2016, 2017). The anisothermality in the nighttime temperatures is assumed to result from the slower cooling rocks relative to the surrounding regolith and the modeled rock abundance represents the surface fraction of rocks that are large enough to be thermally isolated from the regolith (≳50–100 cm) (Bandfield et al., 2011; Hayne et al., 2017).

The rock abundance (RA) map was used to define areas of high rock abundance (RA > 0.10), moderate rock abundance (0.10 ≥ RA > 0.05), and low rock abundance (RA ≤ 0.05) (Figure 1d). Similarly, the regolith temperature map was used to define areas of warmer and cooler regolith temperatures in five temperature ranges: <10 K, 10–20 K, 20–30 K, 30–40 K, and >40 K (Figure 1e) where the temperatures are shown relative to average temperatures normalized by latitude. CSFDs were then determined for each of these temperature regimes.

We additionally developed a one‐dimensional thermal model to further quantify variations in thermophysical properties of the ejecta (Appendix A). The model was fit to bolometric temperatures, the equivalent blackbody temperature of the wavelength‐integrated radiance from all seven channels, derived from Diviner observations across a 285‐km2 portion of the ejecta south of the crater rim (102.25°–103.50°E and 35.25°–35.55°N) selected to avoid large slopes and shadowing from the rim while providing the greatest diversity of ejecta materials. Radiance from Diviner Reduced Data Records (RDRs) from 5 July 2009 to 1 May 2020 for channels 3–9 were binned at 128 pixels per degree and 0.1 hr of local time to create a gridded data set for the modeled area. Emission angles were constrained to <20° to minimize phase angle‐dependent radiative behavior of the lunar surface (Bandfield et al., 2015; Warren et al., 2019), while including targeted observations with small emission angles made as part of the Diviner “twilight campaign” to capture post‐sunset observations (Russell et al., 2017, 2019). The effective field of view of each RDR record was modeled and projected onto a digital elevation model during the binning process as described in Williams et al. (2016) and Sefton‐Nash et al. (2017). The binned radiances were then converted to bolometric temperatures as described in Paige, Siegler, et al. (2010).

Thermal diffusion models developed to derive the thermophysical properties of the lunar regolith from Diviner observations have found that the observed nighttime surface temperatures can be well characterized by regolith with density and thermal conductivity profiles that increase exponentially with depth such that where ρ and ρ are the densities at the surface and at depth z ≫ H, where the H‐parameter is the exponential scale height governing the rate of increase in density with depth and with a solid conductivity that is coupled to density by where k and k are the conductivities at the surface and at depth z ≫ H, (Hayne et al., 2017; Vasavada et al., 2012).

The thermophysical properties can be collectively characterized by thermal inertia, I = , where ρ is density, k, is thermal conductivity, and c is the specific heat capacity, which describes the resistance to changes in temperatures of materials. The regolith is generally comprised of low thermal inertia material typically within ∼10–100 J m−2 K−1 s−1/2 at 273 K within the diurnal thermal skin depth (Hayne et al., 2017). Thermal inertia depends on temperature due to the temperature dependence of specific heat capacity and radiative conductivity, which contributes to the overall effective thermal conductivity (Appendix A).

By varying the H‐parameter while holding other values constant, Hayne et al. (2017) mapped the global variations in regolith thermophysical properties by fitting the rock‐free nighttime regolith temperatures of Bandfield et al. (2017). However, at GB and many other young Copernican‐age impact craters, varying the single parameter H is not adequate for fitting the observed temperatures as the ejecta is comprised of materials that substantially exceed the nominal density and conductivity values of the model assumed for regolith. We therefore employ a model that can accommodate higher thermal inertia materials and the diversity of observed terrain types such as impact melt, sizable blocks, granular material with variable clast sizes, and any mixtures thereof.

Analytic models have been developed that attempt to account for the relevant physics of heat flow through regolith (e.g., Sakatani et al., 2017; Wood, 2020). Such models determine heat conduction pathways within the material, which are largely determined by the nature of the thermal contacts between grains. However, the parameters within such models are rarely known a priori and will vary with changes in terrain properties. Therefore, we take a more generalized approach that is agnostic toward the configuration of materials in order to capture the range of diverse thermal inertias encountered in the ejecta without requiring knowledge of poorly constrained model parameters.

For simplicity, we have preserved the overall approach of the Hayne et al. (2017) model, but allow the lower boundary thermophysical parameters ρ and k to vary in addition to H (Appendix A). This extends the range of permissible thermal inertias while retaining the ability of the vertical thermal inertia profile to vary. This also allows the model to still reproduce an increase in density and conductivity with depth, while allowing higher density materials to be at the surface when H = 0. The lower boundary is assumed to be a volumetric mixture of rock and regolith while the nominal surface density and conductivity are that of surface regolith, remaining unchanged from Hayne et al. (2017). Using the weighted arithmetic mean of rock and regolith is such a way represents a parallel model for a two‐component mixture (see Appendix A). However, the choice of mixing model will influence how density and conductivity are coupled and other valid mixing models can be employed (e.g., Carson et al., 2006, 2005; Wood, 2020). Absolute values of density and conductivity derived from model fits should therefore be interpreted with caution. However, this approach does provide meaningful estimates of relative differences in thermal inertia across the ejecta. A discussion on the modeling approach chosen is provided in Appendix A along with model details and the parameter values used are listed in Table A1.

Table A1

Model Parameters

Parameter	Symbol	Value	Reference
Surface layer conductivity	k s	7.4 × 10−4 W m−1 K−1	Hayne et al. (2017)
Deep layer conductivity, regolith	k d,reg	3.4 × 10−3 W m−1 K−1	Hayne et al. (2017)
Deep layer conductivity, rock	k d,rock	1.492 W m−1 K−1	Bandfield et al. (2011)
Surface layer density	ρ s	1,100 kg m−3	Hayne et al. (2013)
Deep layer density, regolith	ρ d,reg	1,800 kg m−3	Hayne et al. (2017)
Deep layer density, rock	ρ d,rock	2,940 kg m−3	Bandfield et al. (2011)
Specific heat capacity	c p	c 0  + c 1 T + c 2 T 2 + c 3 T 3 + c 4 T 4	Hemingway et al. (1981); Ledlow et al. (1992)
Coefficients for specific heat capacity, regolith	c 0	−3.6125 J kg−1 K−1	Hayne et al. (2017)
	c 1	+2.7431 J kg−1 K−2
	c 2	+2.3616 × 10−3 J kg−1 K−3
	c 3	−1.2340 × 10−5 J kg−1 K−4
	c 4	+8.9093 × 10−9 J kg−1 K−5
Coefficients for specific heat capacity, rock	c 0	−154.9 J kg−1 K−1	Horai and Simmons (1972)
	c 1	+4.983 J kg−1 K−2
	c 2	−8.207 × 10−3 J kg−1 K−3
	c 3	+5.192 × 10−6 J kg−1 K−4
	c 4	0 J kg−1 K−5
Surface layer radiative conductivity parameter	χ s	2.7	Vasavada et al. (2012); Hayne et al. (2017)
Deep layer radiative conductivity parameter	χ d	χ d = (1– v rock ) χ s	This study

Results

Crater Populations, Rock Abundances, and Nighttime Temperatures

A total of 31,685 craters were identified and measured on the GB ejecta blanket with diameters D ≥ 10 m with the largest crater D = 224 m. All derived model ages discussed below are summarized in Table 2. The CSFD of the total population yields an AMA of 6.3 Ma (Figure 2) consistent with Morota et al. (2009) who estimated an age of 1–10 Ma using images from the SELENE (Kaguya) Terrain Camera. However, we find that AMA values vary with variations in thermophysical properties. For instance, we divided the crater population into surface areas of high, intermediate, and low rock abundances, where high and low rock abundance areas were defined as having values above 0.10 and below 0.05, respectively (Figure 2b). AMAs for these two areas are 4.5 Ma and 8.4 Ma for the high and low rock abundance areas, respectively, and 6.1 Ma for the areas with intermediate rock abundances. The AMAs for the intermediate and high rock abundance areas used craters D ≥ 25 m as the slopes of the CSFDs for these areas become increasingly shallow at smaller diameters. This downturn in the CSFDs is not seen in the low rock abundance area and is not attributed to the downturn commonly observed as CSFDs approach the image resolution as sub‐10 m craters are easily identified in the images. Discrepancies in ages are greater if craters D ≥ 10 m are used (see Table 2).

Table 2

Summary of AMA Results

	Count Area (km2)	Number of craters	AMA (Ma) a	N(1) (km−2) b	N(10 m) (km−2) c	N(10 m) (km−2) obs. d	Chronology model e	Figure
Total survey area	1,322.62	31,685	6.3−0.04+0.05	5.26 × 10−6	23.92	23.96	NPF	2a
Rock abundance
Low (<0.05)	722.08	23,082	8.4−0.05+0.07	7.02 × 10−6	31.89	31.97	NPF	2b
Moderate (0.05–0.10)	341.67	3,263	4.5−0.05+0.05/6.2+0.3−0.3 f	4.65 × 10−6	17.09/23.54 f	17.35	NPF	2b
High (>0.10)	258.87	1,148	2.5−0.05+0.05/4.5+0.3−0.3 f	3.30 × 10−6	9.49/17.09 f	9.68	NPF	2b
Regolith temperature
<10 K	373.47	13,087	9.2−0.08+0.08	7.70 × 10−6	34.93	35.04	NPF	2c
10–20 K	327.4	9,182	7.3−0.08+0.08	6.16 × 10−6	27.72	28.04	NPF	2c
20–30 K	377.14	3,367	4.3−0.06+0.06/5.6+0.3−0.3 f	4.35 × 10−6	16.33/21.26 f	16.57	NPF	2c
30–40 K	216.52	1,609	3.3−0.06+0.06/5.3+0.3−0.3 f	3.62 × 10−6	12.53/20.12 f	12.64	NPF	2c
>40 K	28.08	201	2.4−0.1+0.1/4.5+0.8−0.8 f	3.49 × 10−6	9.11/17.09 f	9.26	NPF	2c
Terrain types
Clastic	197.41	5,218	6.9−0.01+0.01	5.80 × 10−6	26.2	26.43	NPF	6b
			6.4−0.09+0.09	9.38 × 10−6	26.32	26.43	WPF (regolith)	7b
Impact melts	43.72	310	–	–	–	7.09	–	6b
Blocky	35.14	194	–	–	–	5.52	–	6b
Melt flow deposit	2.51	8	0.90−0.03+0.03	7.58 × 10−7	3.42	3.19	NPF	–
			0.84−0.3+0.3	1.23 × 10−6	3.45	3.19	WPF (regolith)	7b
			3.2−0.1+0.1	9.51 × 10−6	3.44	3.19	WPF (rock)	7b

Uses the Poisson timing analysis of Michael et al. (2016) with uncertainty representing the ±34 percentile of the probability density function for craters D ≥ 10 m excepted where noted.

The extrapolated cumulative number of craters ≥1 km per km2 of the crater chronology model.

The cumulative number of craters ≥10 m per km2 of the crater chronology model.

The cumulative number of observed craters ≥10 m per km2.

NPF (Neukum et al., 2001); WPF (Williams, Paige, et al., 2014 and Williams, Pathare & Aharonson 2014).

AMA determined using craters .

Figure 2

(a) CSFD and AMA for the total crater population measured in our survey of a 1323‐km2 area on the GB continuous ejecta. (b) CSFDs and AMAs for the crater population for surface areas of low rock abundance (surface fraction <0.05; blue), high rock abundance (surface fraction >0.10; red), and intermediate rock abundances (yellow). (c) CSFDs and AMAs for the crater population for surface areas with differing nighttime regolith temperatures (<10, 10–20, 20–30, 30–40, and >40 K; blue to red). The boundaries in the rock abundances and regolith temperatures correspond to the contours in Figures 1d and 1e. The model ages use the Neukum et al. (2001) chronology system. For the total survey, the low rock abundance, and the regolith temperatures <10 K and 10–20 K, AMAs use D ≥ 10 m; all others use D ≥ 25 m.

The count area was similarly divided by nighttime regolith temperatures into five categories from cooler to warmer temperatures in 10 K increments (Figure 2c). AMAs for these areas decrease with increasing temperature: 9.2 Ma, 7.3 Ma, 5.6 Ma, 5.3 Ma, and 4.5 Ma. A similar downturn at the smallest diameters is also observed for the warmer areas. Thus, craters smaller than 25 m are excluded from the AMA calculations for the three warmest areas with regolith temperatures >20 K.

This represents a factor of ∼2 difference in AMAs between the high and low rock abundance areas and the areas with the warmest and coolest regolith temperatures. Consistent with this trend in AMAs, we observe a trend in crater densities with thermophysical properties (Figure 3). Crater densities are observed to generally decline with increasing rock abundance and nighttime regolith temperatures across the ejecta.

Figure 3

The density of craters D ≥ 10 m on the ejecta of GB versus (a) rock abundance and (b) nighttime regolith temperatures relative to average temperatures normalized by latitude. Crater densities are binned with error bars representing the standard deviation.

Crater Populations and Thermophysical Modeling of Terrain Types

The southern portion of the ejecta of GB is particularly diverse in surface morphologies with sizable areas of impact melt deposits and clastic materials varying from dense concentrations of meter‐scale and larger blocks to finer‐scale granular materials as seen at NAC resolution (1.2–1.5 cm/pixel; Table 1). Variations in these terrain properties are also apparent in the thermophysical expression observed in Diviner temperatures. We apply our thermal model to this area to characterize variations in these materials based on nighttime temperatures to further compare crater populations to variations in terrain properties (Figure 4).

Figure 4

(a) Southern portion of the GB ejecta blanket with the RGB composite map of model results. Portions of NAC images m1153560774L and m1122936950L show examples of terrains classified as melt, blocky, clastic stratified, and clastic uniform with 100 m scale bars. (b) The individual maps used in the RGB composite showing: (red) magnitude of nighttime cooling, ΔT ; (green) thermal inertia at 273 K, I ; (blue) H parameter.

Our modeling approach (described in Section 2 and Appendix A) is used to identify a best‐fit H parameter for each 128 ppd bin of data in the area, along with ρ and k to define a best‐fit constant‐temperature (T = 273 K) thermal inertia at depth, I . These quantities have been mapped along with the magnitude of the temperature change, ΔT , between dusk and dawn, which will be sensitive to differences in the vertical stratification, and defined here as the difference in mean temperature during local times 18.5–19.0 hr and 5.0–5.5 hr, respectively. We have mapped these values into an RGB composite image where ΔT has been mapped into the red channel, I into the green channel, and H into the blue channel to characterize lateral and vertical variation in the ejecta properties (Figure 4).

From this mapping, we can broadly classify the materials into four end‐member categories. Materials that are mapped as yellow have large ΔT and I values with no vertical stratification in thermophysical properties in the thermally active layer, that is, H = 0. These areas contain high densities of visible blocks with many blocks on the order of tens of meters in size. Cyan colors correspond to materials that also have high thermal inertias. However, they have high H values, indicating that the high thermal inertia material is not at the surface, and lower ΔT values, resulting from rapid early cooling at dusk followed by slower cooling later in the night, suggesting an upper layer of low‐thermal inertia material overlying higher‐thermal inertia material. These surfaces contain visible impact melt deposits in NAC images. Materials that are mapped as red‐orange have very low, or zero, H values and I values that are lower than both the impact melt and blocky surfaces. They consist of a mixture of granular materials that include meter‐scale and smaller rocks. The zero H values indicate a uniform, homogenous vertical mixture of thermophysical properties for a given temperature. Areas that map as magenta are similar, but have modest, nonzero H values, suggesting some vertical stratification of properties. These areas appear generally devoid of visible blocks at the surface in NAC images. Data from individual 128 ppd bins representative of each of these terrain types are plotted in Figure 5 along with the fit models showing the nighttime surface cooling and corresponding vertical I profiles for these terrains. Note in this figure that approximately 3 hr after sunset local time, the surface temperatures of the melt and the uniform clastic materials cross due to differences in vertical structure (i.e., stratified vs. uniform).

Figure 5

(a) Nighttime bolometric temperatures observed by the Diviner for four terrain types with model curves derived using a parallel mixing model of regolith and rock (Section 2 and Appendix A) and nominal regolith (black curve) for latitude 35.5°. (b) The resulting best‐fit exponential vertical thermal inertia profiles for the models for T = 273 K, I .

Crater count areas were then selected by the I and H values. Mapped thermophysical properties were defined as having low or high thermal inertias, and vertically uniform or stratified thermophysical properties based on cutoff values of I  = 400 J m−2 K−1 s−1/2 and H = 2 cm. Using these criteria, we mask out three areas that correspond to regions with impact melt deposits (I  > 400 J m−2 K−1 s−1/2 and H > 2 cm), blocky materials (I  > 400 J m−2 K−1 s−1/2 and H < 2 cm), and submeter, granular clasts (both vertically stratified and uniform: I  < 400 J m−2 K−1 s−1/2) (Figure 6, Table 2). The CSFDs of these areas generally overlap within ∼5–10 Ma isochrons for D ≳ 25 m. However, at smaller diameters the CSFDs diverge, with a deficit of sub‐25 m diameter craters in the melt deposits and blocky materials. This results in lower crater densities on these terrains, and similar to the rock abundance and regolith temperatures, the crater density also decreases with increasing thermal inertia of materials (Figure 7).

Figure 6

Regions used to derive CSFDs for ejecta materials characterized as (a) stratified and homogenous clastic materials (submeter clasts/granular), (b) impact melts, and (c) blocky materials. (d) The CSFDs for the three areas with 1, 10, and 100 Ma model age isochrons from Neukum et al. (2001).

Figure 7

The density of craters D ≥ 10 m versus the modeled I on a portion of the ejecta of GB (Figure 4). Crater densities are binned with error bars representing the standard deviation.

Discussion

Thermophysical Properties of Terrain Types

The impact melt deposits and blocky materials within our modeled regions were found to have the highest I values (≳600 J m−2 K−1 s−1/2) (Figure 8). However, these areas differ in rock abundance values. The blocky materials have rock abundance values exceeding 0.2 while the impact melts are ∼0.1. The low ΔT values observed in the melt region resulting from the early relatively rapid cooling followed by the relatively slow cooling later, require a low thermal inertia material overlying the high thermal inertia materials. This configuration appears to limit the rock abundance values retrieved in the model of Bandfield et al. (2011). The much slower early evening cooling of the blocky terrains, by comparison, requires the high thermal inertia materials at the surface to have H = 0. Surfaces with rock abundances ≳0.12 are best‐fit by models with H close to, or equal to, zero (Figure 8b). Many of the observed blocks are meter‐scale and larger, exceeding the diurnal thermal skin depth. This suggests that above ∼0.12 surface fraction of rocks, the rocks on the surface dominate the bolometric temperatures, and the model is no longer sensitive to any interspersed or underlying lower thermal inertia materials. At lower rock abundance values, our model is able to distinguish between vertically stratified and homogenized materials within the diurnal skin depth of the ejecta.

Figure 8

Rock abundance versus modeled (a) thermal inertia, I , and (b) H parameter. Symbols representing mapped pixels are colourized by the same RGB model results as Figure 4a with the terrain types labeled.

CSFDs on Terrain Types

Crater density is observed to vary with ejecta materials of differing thermophysical properties. Where materials have a higher modeled thermal inertia, rock abundance, or regolith temperature, the number of craters decrease (Figures 3 and 7). The CSFDs of these terrains (Figures 2 and 6) show a downturn in the power law slope at smaller crater diameters (≲25 m), which is not seen in areas of lower thermal inertia, rock abundances, or nighttime temperatures. This suggests that the observed variations in crater densities across the different ejecta materials are largely due to variations occurring in the smaller, sub‐25 m diameter craters. This can be seen in Figure 9, where we compare the crater frequencies at D ≥ 10 m [N(10)], D ≥ 20 m [N(20)], and D ≥ 25 m [N(25)], for the CSFDs from Figures 2b, 2c and 6b. The crater frequencies are normalized to the isochron for the AMA (6.3 Ma) of the entire survey, N AMA = 6.3, from Figure 2a. The N(10) values for high rock abundance, warm regolith temperatures, and melt and blocky ejecta material type show a deficit in craters relative to the total survey (factor of ∼0.2–0.4). The low rock abundance, cooler regolith temperatures, and smaller‐clast materials all have higher N(10) values than the total survey (i.e., N(10)/N(10)AMA = 6.3 Ma > 1). The deviations in crater frequency from the total survey decrease for N(20) and N(25), confirming that the discrepancies in crater densities largely occur at smaller diameters. Taking only craters D ≥ 25 m, we recalculated the crater density for the survey as before, but now the trend of decreasing crater density with increasing rock abundance, regolith temperatures, and thermal inertia is significantly reduced (Figure 10).

Figure 9

The crater frequencies N(10), N(20) and N(25), representing the cumulative number of craters D ≥ 10, 20, and 25 m, respectively, relative to the isochron for an AMA of 6.3 Ma, which was derived for the entire survey (Figure 2a) for the CSFDs in (a) Figure 2b, (b) Figure 2c, and (c) Figure 6b. Points above the horizontal dashed line are cumulative frequencies at that crater diameter exceeding the 6.3 Ma isochron, and points below the line represent a deficit in craters relative to the isochron.

Figure 10

The density of craters D ≥ 25 m on the ejecta of GB versus (a) rock abundance, (b) nighttime regolith temperatures relative to average temperatures normalized by latitude, and (c) modeled I . Crater densities are binned with error bars representing the standard deviation. The increase in crater density between I ∼ 400–500 is due to a cluster of D > 25 m self‐secondary craters (see Figure 12).

Figure 12

A portion of LROC NAC image m1153560774R centered on 35.50°N and 102.64°E showing a cluster of decameter‐size craters. Several of these craters contain impact melt deposits (example shown with inset) suggesting these are self‐secondary craters formed prior to emplacement of the impact melt. Craters are not observed within the blocky materials.

While much of the heterogeneity in the crater distribution appears to occur for craters D ≲ 25 m, systematic differences in AMAs do persist even when only D ≥ 25 m craters are used as seen in the AMAs in Figure 2, which still vary by as much as a factor of ∼2. This suggests that (a) terrain properties exert an influence on the crater population at all sizes in our survey, and (b) influence increases at smaller sizes due to some size‐dependent factor. We suggest target scaling effects and impactor‐scale heterogeneities in material properties such as the presence or absence of boulders could be influencing the crater production on the ejecta.

Target‐Scaling Effects

Target properties such as density, porosity, and yield strength can influence resulting crater diameters to the extent that geologic interpretations using derived model ages could be affected (e.g., Dundas et al., 2010; van der Bogert et al., 2017). The characteristic yield strength of the target materials will begin to influence the resulting crater volume where it becomes large relative to the lithostatic pressure, ρ gR , where ρ is the target density, g is surface gravity, and R is the projectile radius, taken to be the characteristic depth of the relevant lithostatic pressure (e.g., Holsapple, 1993; Housen & Holsapple, 2011; Melosh, 1989). For the Moon, the transition to this strength‐controlled scaling occurs at crater diameters of several hundreds of meters and will increase for targets with higher yield strengths (Melosh, 1989; Schultz & Spencer, 1979). Since the crater diameters of our survey (10 m ≤ D ≤ 224 m) are well within this strength‐scaling regime, target property contrasts could be influencing the CSFDs in the GB ejecta field.

We demonstrate this further by isolating crater counts from a distinct melt flow on the southern flank of GB (Figure 11), within the area we applied our thermal model, and comparing the CSFD with the craters in the area comprised of granular, clastic materials (Figure 6a). The melt flow was selected to provide optimal contrasting target properties. The flow has been estimated to be several meters thick (Bray et al., 2010), and given the youthful nature of the GB impact, should provide a relatively unaltered, rock target for comparison with the surrounding granular, clastic materials. The largest crater on the melt flow, D = 31 m, suggests this is the case as the crater interior and surroundings have a blocky texture, indicating that the impact generated numerous meter‐scale clasts and fragments.

Figure 11

(a) Impact melt flow on the southern flank of the GB ejecta. White outline is boundary of count area. Inset shows the largest crater found in the count area with coarse, meter‐scale blocky materials within the crater interior and immediate surroundings suggesting a coherent, competent target consistent with a relatively young melt deposit. (b) CSFDs of the melt deposit and the clastic materials from Figure 6a. AMAs are from Williams, Pathare, and Aharonson (2014). The clastic materials have a model age of 6.4 Ma and the melt deposit of 840 ka assuming target properties of nominal regolith. However, assuming hard rock target properties, the melt deposit has a model age of 3.2 Ma.

We derive AMAs for the melt flow and the clastic materials using the model from Williams, Pathare, and Aharonson (2014) based on the observed flux of terrestrial fireballs (Brown et al., 2002). This model can accommodate different target materials, including effective strengths and scaling parameters that have been shown to relate to the coefficient of friction and porosity (Prieur et al., 2017), to account for scaling effects in the production function. Assuming parameters appropriate for lunar regolith (see Table 2 of Williams, Pathare, & Aharonson, 2014) yields similar AMAs as the Neukum et al. (2001) model. Applying this to CSFD of the clastic terrain gives a model age of 6.4 Ma (Figure 11; Table 2). However, the melt flow has an AMA of only 840 ka. Impacts into rock targets yield smaller diameter craters for a given impact event relative to regolith. If we account for this by assuming crater‐scaling parameters for solid rock, we derive an AMA of 3.2 Ma which reduces the discrepancy in age from a factor of 7.6 to a factor of 2. The remaining discrepancy may be due to statistical uncertainty of using a small count area (2.5 km2) (Pasckert et al., 2015; van der Bogert et al., 2015; Warner et al., 2015) or the presence of self‐secondaries generated on the ejecta prior to the emplacement of the melt flow (Plescia & Robinson, 2019; Williams et al., 2018).

Impacts Into Large Clasts and Boulders

The apparent heterogeneities in the target materials as observed in LROC NAC images and Diviner‐derived thermophysical properties likely influence the cratering process. Our thermal modeling confirms that target characteristics vary both laterally and vertically; also, the ejecta materials differ from typical regolith, requiring higher thermal inertia materials to be present throughout the area we modeled. The ejecta surfaces can therefore be interpreted to represent differing concentrations of coarse‐grained materials, larger clasts, and boulders which have experienced relatively minor physical disruption, abrasion, overturn, and gardening by impact bombardment based on GB's young age (Basilevsky et al., 2015; Costello et al., 2018, 2021; Ghent et al., 2014, 2016; Hörz et al., 2020; Ruesch et al., 2020) or fracturing and fragmentation by thermal fatigue (Delbo et al., 2014; Molaro et al., 2017). Meter‐scale boulders are thought to comminute into regolith over about 200 million years (e.g., Basilevsky et al., 2013, 2015); thus, after less than 10 million years, Giordano Bruno's ejecta are in a relatively early stage of their evolution into fine‐grained regolith. The presence of coarse grains, coherent boulders, and heterogeneity in their distribution influences the craters which form on these surfaces.

Laboratory and numerical studies exploring the effects of coarse‐grained targets and target heterogeneities on impact cratering have demonstrated that as target grain sizes or heterogeneities approach length scales comparable to the impactor size, complex interactions develop that can influence crater growth and efficiency and the final crater morphology (e.g., Barnouin et al., 2019; Daly et al., 2019; Durda et al., 2011; Güttler et al., 2012; Susorney et al., 2017; Tatsumi & Sugita, 2018). In coarse‐grained targets, lower energy impacts expend their energy fragmenting the first‐contact grain. With increasing impact energy, surrounding grains will begin to be disrupted and at sufficiently large impact energies, enough high‐speed fragments from the first‐contact grain transfer momentum and kinetic energy to the surrounding grains to develop full excavation flow (Tatsumi & Sugita, 2018). This would suggest that self‐secondary impacts will be more sensitive to target particle sizes as their impact velocities will be below the lunar escape velocity, 2.38 km s−1. If a substantial fraction of the craters on the ejecta are self‐secondaries, as suggested by Plescia and Robinson (2019), the terrain effects may be enhanced relative to a primary‐only crater population.

For the larger clasts and boulders, impactors of similar length scale will expend considerable energy on fracturing or fragmenting the boulders with no recognizable impact crater being formed. This “armoring effect” has been suggested to explain the significant depletion of small craters observed on boulder‐rich asteroids such as Eros, Itokawa, and Bennu (Bierhaus et al., 2022; Chapman et al., 2002; Hirata et al., 2009; Walsh et al., 2019). Fragmented blocks were observed on the ejecta of GB during our crater survey providing evidence that collisional‐fragmentation of blocks has occurred, which is consistent with the findings of Ruesch et al. (2020),who surveyed blocks at several lunar craters, including GB, and found the abundance of fragmented blocks increased with crater age.

Implications for Lunar Chronology

Our AMA estimates varied from 0.84 ka to 9.2 Ma (Table 2). If we assume that the lower AMA values are the result of target scaling effects and the inhibition of crater formation by blocks and large clasts, this would suggest that the older age of this range is more representative of the age of GB. However, a substantial fraction of the craters in our survey could be self‐secondary craters, which complicates the age estimate and would suggest that the higher age estimates are an overestimate of the age of GB (Plescia & Robinson, 2019). If 50% of the craters counted are self‐secondaries, this would suggest an age of ∼4–5 Ma.

These complicating factors are highlighted in Figure 12 showing a cluster of decameter‐size craters on the southwest portion of the GB ejecta (see also Figure 10). Several of these craters contain impact melt materials within their interiors supporting their origin as self‐secondary craters as the craters must have formed prior to the emplacement of the impact melt. However, craters are not observed within the area containing a high density of boulders, possibly due to the armoring effect discussed above.

The influence of terrain properties on the crater production along with the uncertain contribution of self‐secondary craters has implications for chronology modeling. Radiometric and cosmic‐ray exposure ages of Apollo and Luna samples (e.g., Arvidson et al., 1975; Stöffler & Ryder, 2001, 2006), correlated with crater populations, anchor the lunar crater chronology and enable systems of crater retention age isochrons to be developed for other solar system bodies (e.g., Hartmann, 1999, 2005; Neukum et al., 2001; Neukum & Ivanov, 1994; Neukum & Wise, 1976). Much of the uncertainty in this technique results from a complete absence of lunar samples of known provenance from locations with surface ages of ∼1–∼3 Ga, representing a substantial observational gap in the chronology (e.g., Robbins, 2014; van der Bogert & Heisinger, 2020). This has been recently alleviated by samples returned by the Chang'e−5 mission of basalts from the Oceanus Procellarum region which have yielded ages of ∼2 Ga (Che et al., 2021; Li et al., 2021). However, all absolute ages for Copernican age (<1.1 Ga)‐sampled locations are derived from the comparison of exposure ages with crater counts on the ejecta of impact craters. Similar processes that resulted in the crater population being nonhomogenous at GB most likely also have influenced craters around these Copernican‐age craters, further adding to the uncertainty of chronology models. For example, counts conducted by Hiesinger et al. (2015) in 9 different areas around Cone crater provide ages ranging from ∼16 to ∼82 Ma; counts in a limited area by Williams, Pathare, and Aharonson (2014) produced an age of ∼25 Ma and Plescia & Robinson (2011) derived an age of ∼73 Ma. This suggests a critical need to better understand crater production on the ejecta of Copernican‐age craters and for additional locations to be sampled to provide greater confidence in crater chronology modeling, which is the primary method used for age‐date planetary surfaces.

Conclusions

Giordano Bruno possesses a relatively young and well‐preserved surface morphology that has undergone minimal disruption from subsequent impacts, providing an early stage view into the formation of craters on all types of crater ejecta. We used Diviner‐derived rock abundances and nighttime regolith temperatures along with thermal inertias to identify and isolate areas of the ejecta based on thermophysical properties for comparison with our crater counts. Our thermal model demonstrated the existence of lateral and vertical variations in thermal inertia. We found that the number of craters decreases in areas of higher thermal inertia, rock abundance, and nighttime temperatures with CSFDs that yield younger absolute model ages. This discrepancy increases at smaller crater sizes and we suggest that target material scaling effects and heterogeneities play a role. The areas with the highest rock abundance values correspond to areas with high numbers of boulders visible in NAC images and the lowest crater densities, suggesting boulders are inhibiting the formation of craters.

Many of the craters in our survey are likely to be self‐secondary craters. However, regardless of whether craters formed as primary or secondary craters, the cratering process will be subject to the influence of target effects — indeed, secondary crater production will be more strongly inhibited by boulders, large clasts, and other heterogeneities as their impact energies will be smaller than primary impacts. Therefore, regardless of the crater origin, terrain properties will exert a strong influence on the size‐frequency distribution of craters.

We conclude that terrain properties at GB have greatly affected the distribution of craters on the ejecta, as we observe almost order‐of‐magnitude discrepancies in derived model ages. This adds to the growing body of evidence (e.g., Ashley et al., 2012; Hiesinger et al., 2015, 2012; Plescia & Robinson, 2019; van der Bogert et al., 2017; Williams et al., 2018; Xiao, 2018; Zanetti et al., 2017) that crater populations superposed on the ejecta of craters can exhibit substantial heterogeneities. This has implications for chronology models as the young ages are tied to crater counts on the ejecta of Copernican‐age craters. It is likely that such heterogeneities as observed at GB also exist to varying degrees at these other locations used to calibrate the lunar chronology.