Global Drivers and Transport Mechanisms of Lunar Rockfalls

The long‐ and short‐term drivers and transport mechanisms of lunar rockfalls are currently not well understood, but could provide valuable information about the geologic processes that still shape the surface of the Moon today. Here, we compare the global distribution of rockfalls with relevant geophysical data, such as seismic, topographic, thermal, gravity anomaly, and tidal displacement data sets. Rockfalls appear to predominantly occur (a) on equator‐facing slopes and thus in regions with large thermal amplitudes, (b) on slope angles well above‐average (Δ ∼ 10°), and (c) in regions with above‐average rock abundance. We do not observe a qualitatively or statistically relevant relation between rockfall abundance, monitored Apollo‐era shallow seismic activity, and the distribution of visible tectogenetic features. Informed by our global analysis, we conduct a targeted, in‐depth study of 687 rockfall boulders and trajectories in 13 sites across the Moon, including 7 craters, 2 volcanic vents, 2 tectonic structures, and 2 unclassified geomorphic regions. We identify four different source region types, where the type appears to control the occurrence of rockfalls. The source region type in turn is controlled by surface age rather than geomorphic context. We find that rockfall trajectories are mainly controlled by the trigger energy and the geometry of the slope. Our results suggest that erratic small‐scale impacts (mainly in old, Imbrian‐Nectarian, shallow terranes), aided by solar‐induced thermal fatigue of fractured bedrock (mainly in young, Copernican‐Eratosthenian steep terranes), were the dominant, global‐scale long‐ and short‐term drivers of rockfalls in the Moon's recent geologic past.

In the Apollo-era (1960s and1970s) there were many potential mechanisms hypothesized as the drivers of lunar rockfall detachment, including erosion through micro-impacts, thermal breakdown, moonquakes, and impact events (e.g., Hovland & Mitchell, 1973;Titley, 1966). More recently, Xiao et al. (2013) used LRO NAC imagery to map ∼50 rockfalls across the surface of the Moon, and argue that seismic activity and meteoritic impacts are the main long-and short-term causal drivers of mass wasting features, including rockfalls. Here, a "long-term driver/causal factor" or "preparatory causal factor" describes ground conditions and processes that reduce (slowly or rapidly) the stability of potentially unstable rock compartments until failure occurs (weathering), while a "short-term driver/causal factor" or "trigger" describes a process that ultimately initiates the release of a boulder (erosion; see, e.g., Popescu, 1994). The usage of the terms "long-and short-term driver" is required, as the available data cannot be used to distinguish clear triggers for the vast majority of cases. Kumar et al. (2013Kumar et al. ( , 2016 identified a number of rockfalls in Schrödinger Basin and suggested that these events have likely been triggered by seismicity and meteoritic impacts, due to their proximity to thrust faults and small impact craters, respectively. More recently, Kumar et al. (2019) and Mohanty et al. (2020) looked at rockfalls in Laue crater and in the proximity of Mare Orientale which are located next to two extrapolated Apollo-era shallow moonquake epicenters from January 3, 1975 and December 9, 1972, respectively. They argue that the presence of rockfalls in these regions is connected to seismicity along lobate scarps, wrinkle ridges, and ring faults in the area, i.e., are predominantly seismically triggered. Similarly, Watters et al. (2010) and Watters et al. (2019) use the abundance of rockfalls on or close to tectogenetic features, such as lobate scarps as indicators for present-day seismic activity in the respective regions. In contrast, Valantinas and Schultz (2020) do not observe any rockfalls on or in the direct vicinity of presumably seismically active wrinkle ridges (extensional faults), adding that the slope angles of these tectonic features might be too low for boulders to start moving, even when experiencing ground motion.
Rockfalls have been studied on other planetary bodies as well, such as Mars. Roberts et al. (2012) noticed that rockfall and boulder abundance and size vary across the Cerberus Fossae fault system, leading them to suggest that these events have been triggered by local seismic shaking. Due to an apparent lack of uniformly distributed boulder sizes, they specifically exclude temperature-and climate-related triggers for rockfalls in the Cerberus Fossae fault system. In strong contrast, Tesson et al. (2020) analyzed the spatial distribution of a series of rockfalls as a function of latitude and found that events mainly occur on N-and S-facing, as well as equator-facing slopes, indicating that solar-induced thermal breakdown of bedrock plays a crucial role in rockfall occurrence on Mars. Martian and lunar rockfalls could potentially share their physical drivers, as in contrast to Earth, both planetary bodies lack precipitation events in their most recent geologic past, and feature substantial evidence of impact-related processes as well as strong temperature amplitudes. Previous studies observed that solar-induced thermal breakdown can drive rockfall occurrence on Earth as well (Collins et al., 2018;Collins & Stock, 2016).
Once detachment has occurred, lunar rockfalls rapidly move downslope, leaving distinctive tracks across the lunar surface. On Earth, rockfall transport generally initiates with a period of free fall, followed by bouncing, rolling, and sliding, depending on boulder shape and topographic gradients (Evans & Hungr, 1993). On the Moon, energetic triggering (through meteorite impacts) and generally low topographic gradients may influence rockfall transport mechanisms, although little is known about this at present. The tracks of displaced boulders were used to estimate some of the geomechanical properties of the lunar regolith in preparation of the Apollo missions (e.g., Hovland & Mitchell, 1973;Moore, 1970). More recently, lunar rockfall tracks mapped in NAC imagery have been used to estimate geomechanical properties of regions of increased exploration interest, such as large pyroclastic deposits (Bickel et al., 2019), sunlit south polar regions (Bickel & Kring, 2020), and south polar permanently shadowed regions (Sargeant et al., 2020).
It is important to note that all conclusions concerning rockfall drivers, on the Moon as well as on other celestial bodies, are based on observations and data sets with significant spatial limitations, and thus suffer from a strong observational bias. In particular, the majority of studies do not distinguish between longand short-term rockfall initiation processes, and do not present detailed analyses of the characteristics of rockfall source regions. Additionally, studies of lunar rockfall transport mechanisms have mostly focused on specific geographic regions, and an understanding of the prevalence of the various kinematic modes, typical runout lengths, boulder shapes, and slope angles at deposition across the surface of the Moon is largely missing.
Recently, Bickel et al. (2018) and Bickel, Conway, et al. (2020) developed a method to automate rockfall mapping in satellite imagery using a convolutional neural network (CNN), and used that tool to map the global distribution of lunar rockfalls in more than 250,000 NAC images. For the first time, a global study of lunar rockfalls and their drivers was possible (Bickel, Jordan, et al., 2020). Based on the observation that the total number of rockfalls as well as the spatial density of rockfalls is significantly increased in impact craters, Bickel, Jordan, et al. (2020) argue that impact processes are the main driver of lunar rockfall occurrence. However, the majority of lunar rockfalls did not have an obvious trigger, and other preparatory factors (long-term drivers) may contribute to the occurrence of rockfalls.
This work aims at analyzing the 136,610 rockfalls mapped by Bickel, Jordan, et al. (2020), in order to better understand the boulder source regions, the short-and long-term initiation, and the transport mechanisms that govern lunar rockfall occurrence on a global scale. This is achieved through combining the rockfall catalog from Bickel, Jordan, et al. (2020) with global information about the geologic, seismic, thermophysical, topographic, gravitational, and tidal properties of the Moon. Moreover, we perform a detailed local analysis of 687 rockfalls across 13 selected focus regions to study details of rockfall source regions and transport 10.1029/2021JE006824 4 of 26 processes which are not resolvable on large spatial scales. The results of this work have important implications for understanding erosion processes on the Moon and other airless planetary bodies.
The article is structured as follows. Section 2 introduces the data sets that have been used for the global and local scale analyses. Section 3 describes the applied methodology and the selected focus regions. Section 4 reports the results, which are discussed in Section 5 and concluded in Section 6.

Data Sets
We use a variety of existing, global geophysical data sets for this study. Some of these data sets are only used for the global analysis, some only for the local analysis and some for both. In order to avoid feature distortion issues at very high latitudes (e.g., affecting slope aspect determination), all direction-sensitive, non-latitude dependent data sets, are re-projected to a conformal Mercator projection before they were used for the global analysis. All other data sets used for the global analysis are in a default equirectangular projection (Moon 2000). All data sets used for the local study are either in a regular equirectangular projection (Moon 2000) or a polar stereographic projection (Moon 2000), depending on the geographic location of the area of interest, as further discussed in Section 3.1 and 3.2. Bickel, Jordan, et al. (2020) created a global catalog of lunar rockfalls using a CNN, which mapped 136,610 rockfalls. This data set consists of point locations of each mapped rockfall, scattered across the surface of the Moon from ∼70°N to ∼70°S, as well as CNN-derived meta annotations, such as estimated boulder diameter. In this work, we only consider rockfalls between 60°N and S, as well as with CNN confidence values larger than 0.6, derived with images with solar incidence angles between 10° and 60°.

Auxiliary Data Sets
We used subnadir and oblique high-resolution imagery from LRO's NAC (e.g., Robinson et al., 2010) for the detailed local geologic site analyses. Subnadir images were processed and map-projected using ISIS3, either to an equirectangular or a polar stereographic projection, depending on the location of the area of interest. Processing steps included the default steps: lronac2isis, spiceinit, lronaccal, lronacecho, maptemplate, cam2map, and gdal_translate. The spatial resolution of the projected subnadir images range from 0.49 to 1.53 m/pixel. Oblique images were downloaded as PTIF versions and manually stitched to form NAC leftright mosaics. The spatial resolution of the oblique images is difficult to define due to the high emission angle (ground sampling distance changes significantly across the image), but is lower than the resolution of the sub-nadir NAC imagery due to the increased distance between spacecraft and lunar surface.
Furthermore, we use the 118 m/pixel LOLA DEM (Lunar Orbiter Laser Altimeter Digital Elevation Model; "Moon LRO LOLA DEM 118m"; USGS, 2014) for the global and local analyses related to topography. We include a DIVINER-derived rock abundance map retrieved from the PDS with a spatial resolution of ∼250 m/pixel (e.g., Bandfield et al., 2011). For all analyses related to lunar surface temperature, we made use of global thermal maps created by Moseley et al. (2020) from LRO DIVINER data, specifically maps of daily minimum temperature, daily maximum temperature, and the daily temperature amplitude. Here, the daily temperature amplitude describes the difference between the daily minimum (midnight) and maximum (noon) temperature value at each point on the lunar surface. We further use three thermophysical (anomaly) maps generated by a Variational AutoEncoder (Lord VAEder, Moseley et al., 2020). These maps entangle (and thus represent) different physical processes: (a) solar thermal onset, (b) thermal conductivity, and (c) solar effective albedo. All six temperature-related maps have been constructed from ∼425 billion level 1 DIVINER measurements taken over 9 yr that have been resampled to ∼250 m/pixel (modified from Moseley et al., 2020). The thermophysical maps are described in more detail below: Solar thermal onset describes the time of local sunrise, that is, this map visualizes how east-facing slopes are illuminated first, and west-facing slopes are illuminated last during a lunar day. Thus, the solar thermal onset map distinguishes E-and W-facing slopes, as well as non-inclined regions. We point out that this data set is not sensitive to N-and S-facing slopes. The thermal conductivity map describes how the surface maintains its temperature over the course of a lunar night, that is, it describes the thermal inertia (and thermal conductivity) of the surface, where young, rocky regions usually have high thermal conductivity. Therefore, the thermal conductivity map distinguishes rocky and non-rocky regions. The solar effective albedo map describes how hot the surface becomes during lunar noon (peak temperature), which is a function of the albedo of the surface and the influence of local topography, which affects the amount of received solar illumination. In the context of the present work, the solar effective albedo map distinguishes latitudinal location, albedo (e.g., highland vs. mare), and N-as well as S-facing slopes (see Figure S4 in Supporting Information S1 for a visualization of these maps).
We include maps of the free-air (unprocessed gravity anomalies) and Bouguer gravity (gravity anomalies with topography removed) disturbances created using the GRGM1200A gravity model (e.g., Lemoine et al., 2014). These two maps have a spatial resolution of ∼2,000 m/pixel. Additionally, we utilize maps of estimated lunar radial tidal displacements generated by Thor et al. (2021), which is based on LOLA altimetry data that has been processed to retrieve vertical surface displacements (tides), and of nightside surface albedo (Lyman-α wavelength bands), which is based on Lyman Alpha Mapping Project (LAMP) far-ultraviolet data (e.g., Gladstone et al., 2012). The latter has a spatial resolution of ∼1,000 m/pixel.  Table S1 in Supporting Information S1.

Global Geographic Information System Study
In order to understand potential long-as well as short-term rockfall drivers across the surface of the Moon, we first imported all above-listed global, auxiliary data sets into QGIS (https://qgis.org/de/site/). Using the LOLA DEM (USGS, 2014), we created two additional data sets for the analysis, specifically a global slope aspect and slope angle map, using QGIS' built-in functions "Raster/Analysis/DEM/Slope" and "Raster/Analysis/DEM/Aspect." We subdivide each global map into three individual, equal-sized, geographic regions (latitude bands): (a) equatorial (20°N-20°S), (b) northern sub-equatorial (20°N-60°N), and (c) southern sub-equatorial (20°S-60°S). We limit the global-scale analysis to latitudes between 60°N and S, as the used rockfall catalog has been derived using images taken with incidence angles between 10° and 60°, that is, topographic depressions beyond 60°N and S could be partially shadowed.
To analyze global drivers of lunar rockfall, we compare the distribution of the attributes of a given geophysical data set at the locations of deposited rockfalls with the distribution of attributes of the entire data set, per latitude band (as defined above). If both distributions match, then the attributes derived at the rockfall locations are sampled randomly, meaning that this particular geophysical property does not influence the distribution of rockfalls in a given latitude band (with the spatial resolution and quality available). In contrast, a mismatch of both distributions could indicate that this particular geophysical property does have an influence on the rockfall distribution. To do so, we used the rockfall catalog derived by Bickel, Jordan, et al. (2020) to extract the values (attributes) of each global map at each rockfall location. We refer to the distributions of these extracted attributes (values derived at all rockfall locations) as "distribution of data set <type> at-target," or short "at-target." Similarly, we refer to the attributes derived from the entire data set as "background distribution of data set <type>," or short "background." Next, we investigated the potential role of lunar seismicity as a short-and long-term driver. To do this, we analyzed our global database of lunar rockfall using similar methodologies that previous researchers have used to argue that seismicity is an important driver locally (e.g., Kumar et al., 2019;Mohanty et al., 2020). As summarized in the introduction, this previous work has been partly based on spatial correlations between the distribution of rockfall and known tectonic features and shallow moonquake epicenters (based on Apollo-era seismic data). In the present work, we perform two separate spatial correlations: (a) a visual, qualitative comparison of the locations of rockfall hotspots with vector data that describes the geographic position and extent of visible lobate scarps (digitized from Watters et al., 2019), wrinkle ridges (Thompson et al., 2017;Valantinas and Schultz, 2020), graben (Nahm et al., 2018), multi-ring basins larger than 200 km in diameter (Kadish et al., 2011), and the extrapolated Apollo-era shallow moonquake epicenters (Nakamura et al., 1979;Watters et al., 2019); (b) a statistical, quantitative analysis of the relation of spatial rockfall frequency and seismic event magnitude or stress drop.
For (b), we first derived the distances between the 28 original Apollo shallow moonquake epicenter locations (Nakamura et al., 1979), the 12 recently confirmed, re-located Apollo-era epicenter locations (Watters et al., 2019) and all rockfalls mapped by Bickel, Jordan, et al. (2020). We then used three different spatial buffers (0.5°, 1°, and 5° radius around each epicenter, i.e., ∼15, ∼30, and ∼150 km at the equator, respectively), and determined the number of rockfalls within each buffer. We specifically highlighted the strongest recorded shallow moonquake (magnitude 3.2) from January 3, 1975, as well as the December 9, 1972 event close to Mare Orientale (magnitude 1.2). We then assessed whether the rockfall density around epicenter locations is anomalously high, by comparing rockfall densities at epicenter locations to rockfall density at 1,000 randomly selected locations, using the same three spatial buffers. Finally, we compare the relation between seismic event magnitude and the number of rockfalls around each epicenter, within the three spatial buffers. Due to the nature of the available data sets, distance has been measured in degrees (°).

Local Geologic and Geomorphic Analysis
In order to further study the characteristics of rockfalls on a local scale we selected and performed measurements in 13 areas of interest (AoIs) across the Moon (Bickel, 2021). These specific locations have been selected because they are rockfall hotspots and belong to one of the four relevant lunar geomorphic (sub) classes as outlined by Bickel, Jordan, et al. (2020): crater (impact-caused, -induced, and -ejected), volcanic vent, tectonic structure, or other, unclassified geomorphic context. We included an increased number of impact-related AoIs because the vast majority of rockfalls are located in impact craters (Bickel, Jordan, et al., 2020). We manually mapped 687 rockfalls in these AoIs, where we expect that an estimated ∼60%-70% of these manually mapped features are also present in the global rockfall catalog derived by Bickel, Jordan, et al. (2020) used for the global study, based on the recall (percentage of rockfalls found in a testset) of the deployed CNN (Bickel et al., 2018). All AoIs are characterized in more detail in Table 1 and their locations are shown in Figure 2. Local overview maps of all AoIs are showcased in Figure S6 in Supporting Information S1.
The local geologic analysis consists of a series of detailed observations and measurements made in high-resolution, subnadir NAC imagery. Where available, oblique NAC imagery has been used to gain a qualitative understanding of the overall topography of an area of interest as well as of the structure and geology of the source and deposition regions of rockfalls, thus, complementing the main analysis. Using subnadir NAC imagery, we measured a large number of rockfall and track characteristics, specifically: boulder short axis, boulder long axis, boulder shadow length, track length, slope angle at source and at deposition, and the elevation drop from source to deposition (using the LOLA DEM).
We also recorded several qualitative aspects of the AoIs and rockfalls, such as the texture and roughness of the slope substrate along each rockfall trajectory (regolith and debris), the rockfall kinematic (bouncing, rolling, and sliding), the geologic characteristics of the source region (cliff, outcrop, boulder field, and plain slope), if a rockfall deviated off the steepest descent path during its displacement (yes or no), the geologic unit, and the geologic age of the host surface (Pre-Nectarian, Nectarian, Imbrian, Eratosthenian, and Copernican). Using these observations, we calculated a number of additional values, specifically: boulder height (using the shadow length and solar incidence angle), boulder volume (assuming an ellipsoidal shape), approximated boulder (macro-) sphericity (assuming an ellipsoidal shape, using the three orthogonal boulder axes), and the reach angle or Fahrböschung (the relation of vertical drop height and horizontal runout, e.g., Heim, 1932) of each rockfall. We also derived accurate trajectories of two individual rockfalls, one in an impact-caused and one in an impact-ejected setting. For these two events, we created map-and longitudinal-views of their trajectories, which describe the spatial variation of jump heights, widths, and kinematics.
We observe that the vast majority of rockfalls begin their displacement by sliding or rolling. In order to understand at what slope angle boulder movement initiation can occur without additional external forces (assuming exclusively gravitational acceleration) we apply a simplified limit-equilibrium planar slope failure analysis with: assuming that sliding occurs when F < 1, that is, the driving forces (governed by the local slope angle  E ) are larger than the resisting forces (governed by the regolith's angle of internal friction  E ), neglecting sinkage of the boulder. The literature reports a wide range of mare and highland regolith friction angles, ranging from 13° to 50° (e.g., Carrier et al., 1991). We use an estimate of    29 E as derived by Bickel et al. (2019) who estimated the average angle of repose of lunar regolith by measuring the slope angle of inclines that are partially covered by deposits of regolith that flowed downslope, using high-resolution, LRO NAC-based elevation models. Boulders located on slopes that are steeper than the critical slope angle for sliding will start to move unless restrained by other factors, such as cohesion, micro-topography, or other obstacles, which may degrade through time.

Global Analysis
The global, GIS-based analysis provides a series of important findings, summarized in Figure 3- Note. This study follows the typology for impact-driven rockfalls as recently suggested by Bickel, Jordan, et al. (2020), that is, impact-caused-A rockfall occurs in a crater as consequence of substantial, impact-caused fracturing of the upper crater rim but without obvious trigger-and impact-induced-A rockfall is triggered by a close, small-or large-scale impact event-, and impact-ejected-A rockfall is directly ejected from an impacting, small-scale meteoritic body. Rockfall density per Bickel, Jordan, et al. (2020). Bickel, Jordan, et al. (2020) background distribution of slope aspect is uniform, that is, equally distributed over north, east, south, and west. In contrast, the overall at-target (rockfall) slope aspect distribution is not uniform: around the equator more rockfalls appear to be located on N-and S-facing (equator-facing) slopes, that is, on slopes that receive more solar insolation and thus experience larger thermal amplitudes. Interestingly, E-and W-facing slopes appear to host more rockfalls in the sub-equatorial regions ( Figure 3). The slope angle at-target distribution reveals that more rockfalls are located on steeper-than-average slopes, over all studied latitudes. Interestingly, the majority of rockfalls are located on slopes with angles less than ∼20°. It is important to note that all at-target values are derived at rockfall deposition locations, that is, these slope angles do not represent the global source region slope angle distribution. This will be addressed further when we present the results of the local analysis of the areas of interest (Section 4.2). The absolute topographic elevation background and at-target distributions match relatively well, indicating that the absolute elevation does not significantly influence the abundance of rockfalls. However, we note that there is a lack of rockfalls in equatorial regions at elevations of ∼−500 to ∼−2,000 and ∼3,000 to ∼7,000 m; these elevations are predominantly occupied by the nearside maria and the highest regions of the farside highlands, respectively. In contrast, there is a relative abundance of equatorial rockfalls between ∼0 and ∼2,000 m, that is, terrain occupied by highlands and a relatively large number of impact basin walls. In the northern sub-equatorial regions we observe an over-abundance of rockfalls at elevations from ∼−1,800 to ∼−5,000 m; we attribute this to the presence of a large number of young, intra-mare impact craters which are among the locations with the highest spatial 10.1029/2021JE006824 9 of 26 density of rockfalls on the lunar surface (e.g., craters Buerg, Hercules, Atlas, and Aristoteles; Bickel, Jordan, et al., 2020; Figure 3). As in equatorial regions, there is an apparent lack of northern sub-equatorial rockfalls at elevations between ∼1,000 and ∼5,000 m, that is, in the flat portions of the farside highlands. In southern sub-equatorial regions, there appears to be a lack of rockfalls at elevations between ∼−3,000 to ∼−5,000 m, in the lower sections of the South Pole Aitken Basin (SPA).

Table 1 Overview of the 13 Areas of Interest Used for Local Geomorphic Mapping, Using Four Distinct Geomorphic Classes: (a) Impact Structures, (b) Volcanic Features, (c) Tectonic Features, and (d) Unclassified Geomorphic Regions, Per
The free-air and Bouguer gravity anomaly distributions generally match well, suggesting that gravity anomalies do not significantly influence rockfall abundance. Interestingly, there appears to be a lack of rockfalls in equatorial regions with Bouguer anomaly values below ∼−500 mGal, that is, in the highest regions of the equatorial farside highlands, agreeing with observations made in the elevation data. Similarly, there appear to be fewer northern sub-equatorial rockfalls in locations with Bouguer values between ∼−100 and ∼−500 mGal, which correspond to the highest sections of the farside highlands. We note a slight overabundance of rockfalls in southern sub-equatorial regions with Bouguer gravity values between ∼−500 and ∼−50 mGal, which could be caused by the two distinct rockfall hotspots located in Mare Orientale and Tsiolkovsky crater. There is a lack of rockfalls in the southern sub-equatorial regions at Bouguer gravity values between ∼250 and ∼500, which corresponds to the lower regions of SPA, as also suggested by observations in the elevation data ( Figure 4). We further observe a slight difference between the tidal radial displacement background and the at-target distributions at various latitudes, but are not able to identify a potential physical reason ( Figure 4).
The at-target distributions derived with the surface minimum and maximum temperature as well as the temperature amplitude generally match well, across all studied latitudes ( Figure 5). We point out that we are only able to sample rockfall deposition locations, not their source regions, meaning that the extracted values are not entirely descriptive of the acting of thermal fluctuations at and around source regions. Figure 5 additionally shows that rockfalls occur in regions with thermal amplitudes ranging between ∼110 and ∼270 K or more, to a maximum of ∼290 K.
The at-target distribution of the solar thermal onset is slightly shifted from the background distribution over all studied latitudes, indicating that more rockfalls occur on slopes which are illuminated in the morning and afternoon of a lunar day versus during noon. In other words, more rockfalls are located on inclined, E-and W-facing slopes than in flat regions. We note that the solar onset map is not sensitive to illumination on N-and S-facing slopes (Moseley et al., 2020). The surface thermal conductivity at-target distribution suggests that more rockfalls tend to be located in regions with higher thermal conductivity, that is, usually young and/or rocky regions ( Figure 6). The subtle, absolute shift of the effective albedo at-target from the background distribution to higher values in equatorial and northern sub-equatorial regions might suggest that (a) in equatorial regions more rockfalls tend to be located in the highlands, which feature slightly higher effective albedo values as the mare regions in the used data set (Moseley et al., 2020, see Figure S4 in Supporting Information S1) and (b) in northern sub-equatorial regions more rockfalls are located on E-W-facing slopes, potentially supporting observations in other auxiliary data sets, such as slope aspect ( Figure 6).
The at-target rock abundance distribution shows that more rockfalls tend to occur in regions with increased rock abundance, that is, in generally rougher terranes. This observation agrees with the observations made in the thermal conductivity data set, mentioned earlier ( Figure 6). The at-target and background distributions of nightside surface albedo generally agree well, while there is a subtle shift in the northern sub-equatorial region (Figure 7), potentially indicating that the northern hemisphere features an increased number of rockfalls in regions with high surface (UV) albedo, that is, in rocky and/or young terranes, which usually are UV-bright. Potentially, this could be caused by the globally distinct rockfall clusters in the northern mare regions, for example, craters Buerg, Hercules, Atlas, and Aristoteles (Bickel, Jordan, et al., 2020). This observation might further support findings in the rock abundance and thermal conductivity data sets ( Figure 6).
For context, all auxiliary maps used for the global analysis are shown in Figures S1-S5 in Supporting Information S1.
In order to better understand the role of tectonic activity and moonquakes we visually analyzed the distribution of rockfall hotspots and visible tectonic structures, including lobate scarps, wrinkle ridges, graben, and basins >200 km, as well as Apollo-era shallow moonquake epicenters (Figure 8). This qualitative analysis shows that the geographic locations of rockfall hotspots, quake epicenters, and visible tectonic structures only partially overlap. For example, the accumulation of rockfall hotspots in the Lacus Mortis and Mare Frigoris regions are relatively well aligned with five Apollo-era epicenters as well as a number of lobate scarps and presumably active wrinkle ridges. However, a significant number of epicenter and   tectonic feature locations do not overlap with rockfall hotspots. Notably, the moonquake with the highest magnitude (3.2) in the Laue crater region does not correlate with a rockfall hotspot.
We performed a statistical test to further improve our understanding of the influence of lunar seismicity on rockfall occurrence, analyzing the number and distance distributions of rockfalls from Apollo-era quake epicenters as well as randomly selected locations. The results, showcased in Figure 9, reveal that the average number of rockfalls and mean distance, standard deviation, and variance for all events, real quakes and random locations, and all spatial buffers are almost identical. Similarly, the test performed using only the Laue crater epicenter data (the Apollo event with the highest magnitude value, 3.2, January 3, 1975) does not show any statistical relation, just as the test with the Mare Orientale event (1.2, December 9, 1972). The rockfall frequency and seismic event magnitude relation for all original and re-located Apollo-era events ( Figure 9) does not reveal any clear trend either.

Local Analysis
To overcome some of the limitations of the global analysis, such as the lack of information about source region types and slope angles, as well as boulder track trajectories, we further study 687 individual rockfalls that originate from 4 different types of source regions in the 13 AoIs considered, using subnadir and oblique NAC imagery. We label these source regions as: (a) cliffs, that is, large-scale, continuous (intact), inclined to (near)vertical bedrock outcrops that extent over tens of meters, (b) outcrops, that is, fractured bedrock outcrops with spatially limited extent, intermixed with closely spaced boulders, (c) boulder fields, that is, accumulations of individual blocks that are not as closely spaced (that are, e.g., part of ejecta blankets), and (d) plain slopes, that is, smooth slopes without signs of bedrock outcrops or larger boulder fields and only a Figure 9. Statistical analysis of moonquake epicenters and randomly selected locations versus the global rockfall distribution as derived by Bickel, Jordan, et al. (2020). Scatter plots show number of rockfalls as a function of event magnitude for the three spatial buffers (0.5° red, 1° blue, and 5° black), the Laue crater, and Mare Orientale events are marked; the table shows statistical descriptors for distance distributions (epicenter-rockfall) for the three spatial buffers.
few scattered boulders. We note that there is a certain overlap between classes, as this classification is only based on qualitative geologic observations made in satellite images with limited spatial resolution. Figure 10 shows examples of the different source region types; in addition, Figure S7 in Supporting Information S1 shows an Apollo 17 ground-based photograph of a boulder field-type source region on the Taurus Littrow South Massif summit. It is important to note that all source region types occur in all geomorphic contexts, but cliff and outcrop source types predominantly occur in geologically young terranes (Copernican).
The analysis of the collected measurements across the AoIs (Figure 11) shows that source region slope angles are very similar across all analyzed geomorphic contexts, ranging between ∼15° and ∼35° with a mean of ∼24.9° across all regions (excluding impact-ejected rockfalls), which is close to the approximated angle of repose of regolith on lunar slopes (29°) as estimated by Bickel et al. (2019), among others. Only impact-ejected rockfalls feature a wider range of source region slope angles, from ∼0 to ∼35°, as they occur wherever an impactor strikes-they do not require a topographic gradient (mean source region slope angles are shown in Figure 11). A comparison with the globally derived slope angle values is shown in Figure S8 in Supporting Information S1. Using these measurements, we evaluated Equation 1 using the average source  (4) plain slope (images of a pre-Imbrian-aged crater, Gibbs, 18.4°S-84.3°E) source region type. Sketches feature bedrock (dark gray), exposed bedrock and boulders (black), surface regolith (gray), and subsurface fracture networks (light gray). Relevant features are indicated by red arrows. Typical rockfall deposition regions in sub-nadir NAC imagery (bottom), from left to right: sliding-type (Rima Ariadaeus, 7.2°N-11.1°E), rolling-type (Archytas, 58.9°N-4.2°E), and bouncing-type (Giordano Bruno, 36°N-102.8°E). Image credits: LROC/ ASU/GSFC. slope angle for each geomorphic context (and an internal friction angle of 29°), which resulted in F values of between 1.11 and 1.29. This basic, numerical experiment indicates that mean source region slope angles in all geomorphic AoIs are close to or slightly below the angle of repose (∼29°). We note that there are source regions with F < 1 across our AoIs, that is, some source regions are steep enough for boulders to start moving without external force applied, considering the assumptions made.
Displaced boulders across all geomorphic regions feature similar shapes (sphericity between 0.8 and 1 assuming an ellipsoidal shape) and volumes (up to ∼1,000 m 3 assuming an ellipsoidal shape) except for volcanic regions, which feature slightly more angular (sphericity between 0.6 and 1) and smaller boulders (<<1,000 m 3 ; average shape values are displayed in Figure 11). This could point to differences in the fracture networks (extent and spacing) that produce the boulders that ultimately get displaced. Boulders have very similar average dimensions across the AoIs as well, ranging from 5.95 m × 3.5 m × 3.0 m (volcanic-related) to 8.56 m × 5.82 m × 6.36 m (impact-related; see Figure 11), agreeing with the global distribution of boulder diameters (∼7-10 m) as derived by Bickel, Jordan, et al. (2020).
The vast majority of rockfalls across the 13 AoIs do not have an obvious trigger, independent of the geomorphic context. Only rockfalls that belong to the impact-induced and impact-ejected class feature distinct triggers, specifically small-scale, meter to hundreds of meter-sized, meteoritic impactors that either eject rockfalls directly or cause strong ground motion which in turn activates boulder detachment and/or displacement downslope (see Figure S9 in Supporting Information S1).
The track length of rockfalls varies across geomorphic regions: in impact craters and unclassified geomorphic contexts runout can reach up to ∼4,000 m and more (average is ∼1,070 m); impact-ejected rockfall runout can reach up to ∼2,000 m (average is ∼430 m); in tectonic regions runout does not exceed ∼2,000 m (average is ∼670 m), and in volcanic regions runout does not exceed ∼1,000 m (average is ∼270 m; summarized in Figure 11). It is important to note that runout lengths could be influenced by differences in topographic scale across the different regions, as volcanic and tectonic regions do not feature as drastic topographic gradients and absolute elevation differences as impact craters and unclassified regions (e.g., valleys). The longest track lengths are achieved by boulders that either have high sphericity values, predominantly bounce, and/or originate from steep source regions. Interestingly, boulders with larger volumes do not appear to travel further than boulders with smaller volumes (see Figure S10 in Supporting Information S1).
We observe that boulders only bounce on slopes steeper than ∼20° and move by rolling or sliding as soon as the slope gradient declines ( Figure S11 in Supporting Information S1). Across our AoIs we note that rolling is the most common kinematic mode, where 79.6% of all mapped rockfalls predominantly roll, 20.1% bounce, and 0.3% slide. Most rockfalls across the AoIs feature similar movement types at the beginning of their trajectories: We observe that usually boulders start to roll, gain speed, potentially start bouncing, eventually slow down due to a slope angle decline, and then roll or slide until they finally stop moving. There is no indication that the approximated boulder (macro-) sphericity (as derived using only the long and short boulder diameters as well as the height) influences the displacement mode. The only exception to the usual sequence of movement types as described above are impact-ejected rockfalls: usually, impact-ejected boulders are initially airborne, then bounce over the surface in slowly declining rebounds, eventually slow down and finally roll or slide to a halt. Figure S11 in Supporting Information S1 showcases the trajectories of one impact-ejected and one impact-caused rockfall, visualizing the described observations and differences.
We observe that the majority of rockfall trajectories follow the path of steepest descent along their complete trajectory. However, we note a series of rockfalls that deviate off steepest descent, specifically in the Rowland, Ingalls, and Vavilov crater AoIs. The rockfalls in these AoIs either belong to the impact-induced or impact-ejected class, that is, these rockfalls have evident, highly energetic triggers that could influence their trajectories: in the Rowland crater AoI, the trigger impact ejected boulders in all directions, irrespective of the local topography and topographic gradient; in the Ingalls and Vavilov crater AoIs the trigger impact appears strong enough to not only have initiated the observed rockfalls, but to have given the respective boulders a lateral "push," enough to influence their initial descent path (see Figure S9 in Supporting Information S1).
The rockfalls in the different geomorphic contexts show minor differences in the relation between boulder drop height (vertical displacement component) and boulder track length (horizontal displacement component), that is, the reach angle or Fahrböschung. The overall mean Fahrböschung angle is 18.4°, while rockfalls in volcanic and tectonic regions feature angles of 17.41° and 15.08°, respectively, below the overall mean Fahrböschung angle. The rockfalls in the other geomorphic contexts feature mean Fahrböschung angles slightly above the overall mean value.

Long-and Short-Term Drivers of Lunar Rockfall
The combination of a global rockfall catalog, detailed mapping of select areas of interest, and a large number of geophysical auxiliary data sets has revealed a number of interesting findings regarding the long-and short-term factors that control lunar rockfall occurrence. First, detailed mapping of the areas of interest has revealed that rockfalls originate from four main source region types: (1) cliff-type, (2) outcrop-type, (3) boulder field-type, and (4) plain slope type, where all source region types are present in all different geomorphic contexts. We note that the occurrence of source region types is controlled by the age of the host terrane rather than the geomorphic class: cliffs and outcrops are common on the steep slopes of, for example, fresh (Copernican) terranes, while boulder fields and plain slopes are common on the shallower inclines of, for example, old (Imbrian and pre-Imbrian) terranes. Here, regolith-dominated source regions (type 3 and 4) can either be the result of continuous weathering of cliff-and outcrop-type source regions, or can be created through impact processes that exhume and deposit boulders on the surface of the Moon (ejecta blankets). We observe that cliff-and outcrop-type source regions (young source regions) are more productive, that is, they appear to release more rockfalls per given source region area. This agrees with observations made in previous work that younger terranes feature higher spatial concentrations of rockfalls (Bickel, Jordan, et al., 2020).
The global analysis shows that rockfalls are generally located in rockier than average terranes and on steeper than average slopes. The measured mean source region slope angles are very similar for all analyzed geomorphic contexts and source region types and range from ∼24 to ∼27° across all regions. These angles are close to the slope angle needed (∼29°) to initiate passive shearing and sliding in the more regolith-supported source region types (3) and (4), based on the limit equilibrium analysis presented above. Across all AoIs, ∼10% of source region slope angles exceed ∼29°. We note that impact-ejected rockfalls are an exception, as these occur regardless of local topography-these rockfalls have been excluded from the planar sliding analysis. The simplified limit equilibrium analysis indicates that a limited number of the analyzed source regions features conditions favorable for passive rockfall initiation, in the absence of any cohesion. In rockier source regions, such as type (1) and (2), cohesion could be provided by rock bridges, and degraded through time by mechanisms such as (micrometeorite) impacts, thermal breakdown, tidal displacements, and local seismic shaking (e.g., Hörz et al., 2020). We note that all source region slope angles have been derived using a DEM with a spatial resolution of ∼118 m/pixel, whereas the average boulder has a diameter of about ∼7-∼8 m; the local, small-scale source region slopes angles could therefore be even steeper than reported here, facilitating rockfall initiation in all source region types.
The globally apparent relation between rockfall occurrence and slope aspect suggests that solar-driven thermal fatigue may be an important global-scale long-term driver of lunar erosion and rockfall. In equatorial regions, Figure 3 shows there are more rockfalls located on N-S-facing slopes. These regions are generally expected to have higher noon temperatures and larger thermal amplitudes than the global average, however we do not observe this ( Figure 5), potentially due to the limited spatial resolution of the available thermal data sets and sampling these distributions at deposition locations, rather than source locations. We additionally derive the temperature amplitude distribution for all source regions as listed in the local rockfall catalog (687 measurements) and compare it to the at-target (deposition) distribution, but are not able to observe any significant difference either, likely caused by the same limitations of the data as discussed above, assuming that the spatial scale of source regions (<50 m) is much smaller than the spatial resolution of the used data sets (>>∼250 m/pixel).
Figure 5 further shows that rockfall sites at higher latitudes feature smaller thermal amplitudes than equatorial sites. The same figure indicates that: (a) a daily temperature amplitude of ∼120 K might be sufficient to drive fracture growth and rockfall occurrence and (b) larger temperature amplitudes might generate a larger rockfall population (in the equatorial regions). Solar-induced thermal fatigue might also act as shortterm causal factor and true trigger of lunar rockfalls in steep (>∼29°), young and exposed types (1) and (2) source regions, where thermal fatigue has been shown to be most effective (e.g., Molaro et al., 2017). Past work has shown that damage accumulation rates appear to decrease in older, more fractured rocks, meaning that thermal fatigue would become a less important driver in older terranes (e.g., Molaro et al., 2017).
The reason for the rockfall at-target slope aspect rotation from N-S to E-W ( Figure 4) at higher latitudes is currently unknown; on Mars, a similar, weaker trend has been observed by Tesson et al. (2020) for rockfalls and by Dundas et al. (2015) for gullies, but has not been specifically recognized or addressed yet. Dundas et al. (2015) mention that martian slope aspect analyses could generally be affected by the sun-synchronous, mid-afternoon orbit of the Mars Reconnaissance Orbiter (mapping bias), but this observation cannot be directly translated to the Moon, as the LRO does not follow a sun-synchronous orbit. Still, the more complicated illumination situation at higher latitudes (>∼40°) could potentially impact the CNN-driven rockfall detection and slope aspect analysis, as discussed in more detail below. Alternatively, the rockfall E-W orientation at higher latitudes might be controlled by the peak interior and surface thermal stresses that have been shown to occur in the local morning and afternoon, that is, on E-and W-facing boulder faces and slopes (Molaro et al., 2017). The important role of directional insolation in driving mechanical weathering of boulders has been observed on Earth (Eppes et al., 2010;McFadden et al., 2005) and Mars (Eppes et al., 2015) and has additionally been used to explain asteroid activity (Molaro et al., 2019). The decreased influence of the effective albedo with increasing latitude could result in a more pronounced effect of the sunrise/sunset-related thermal effects. Another potential explanation for the observed E-W rotation could be that existing impact-generated fracture networks are reworked by the global, near-surface stress field of the Moon as proposed by, for example, Watters et al. (2019) and Matsuyama et al. (2021), which features an extensive E-W trend at latitudes beyond ∼50°N and S. High resolution, outcrop-scale data returned from future missions, such as Lunar Trailblazer (Ehlmann et al., 2021), will help to further address the potential relation between the thermophysical properties of the lunar surface and the weathering/erosion of the Moon.
In contrast to thermal fatigue, the importance of impact events in driving lunar rockfalls is more evident. Earlier work found that rockfalls predominantly cluster in impact features, indicating that the impact process itself facilitates or initiates topographic erosion through rockfalls, among other mass wasting processes (Bickel, Jordan, et al., 2020). Across all geomorphic contexts, small-scale impacts trigger rockfalls, either directly, that is, by ejecting boulders in all directions (impact-ejected type), or indirectly, that is, by inducing the displacement of pre-existing boulders or by triggering the failure of slope compartments on topographic inclines (impact-induced type). Similar processes have been observed by Kumar et al. (2013), Xiao et al. (2013), and Hovland and Mitchell (1973). We note that the vast majority of rockfalls in old terranes (Imbrian and pre-Imbrian) as well as in polar regions appear to be triggered by impact processes (impact-induced and impact-ejected), indicating that terrane age and thus source region type influence rockfall occurrence and triggers. Generally, rocky terranes host an increased number of rockfalls, as indicated by the relation of rockfall distribution and rock abundance, thermal conductivity, nightside surface albedo, and slope angle.
Past work extensively highlighted the seemingly decisive role of lunar tectonic activity in contributing to or controlling mass wasting and rockfall occurrence (e.g., Kumar et al., 2019;Mohanty et al., 2020;Titley, 1966;Xiao et al., 2013). This assumption is mainly based on the apparent spatial proximity of rockfalls to visible topographic expressions of tectonic activity (e.g., Kumar et al., 2016;Watters et al., 2019) and/or to the epicenters of Apollo-era shallow moonquakes (e.g., Kumar et al., 2019;Mohanty et al., 2020). However, all past work is based on locally constrained data, thus, is not representative of larger regions or the entire lunar surface. In order to test whether endogenic activity has been an important, global-scale driver of rockfall, we use the first consistent and global map of lunar rockfalls (Bickel, Jordan, et al., 2020) in combination with relevant auxiliary data and find that a large fraction of Apollo-era shallow moonquake epicenters and geomorphically visible lobate scarp, wrinkle ridge, graben, and basin-associated ring fault features do not spatially align with rockfall hotspots (Figure 8). Some regions show some level of overlap, such as Mare Frigoris and the Lacus Mortis region, but those regions appear to be rare. We further observe that a number of rockfall hotspots are associated with impact basins, but subsequent, detailed observations suggest this to be a consequence of: (a) the large size of impact basins, meaning that they have a larger chance to contain basin-unrelated rockfall hotspots, such as small Copernican/Eratosthenian impact craters that are located within the basin and (b) the impact process itself, that is, they feature impact-generated, subsurface fracture networks that likely drive boulder production and displacement over long timescales (Bickel, Jordan, et al., 2020), rather than seismic activity along potential faults.
Similar to these qualitative observations, our quantitative analysis shows no statistically significant difference between the spatial distribution of rockfalls around known moonquake epicenters and randomly selected locations on the lunar surface. Furthermore, there seems to be no global relation between moonquake magnitude or stress drop and spatial rockfall frequency. This is also the case for the January 3, 1975 and December 9, 1972 events, that is, the moonquake with the largest recorded magnitude as discussed by Kumar et al. (2019) and a quake close to Mare Orientale as discussed by Mohanty et al. (2020). On Earth, the minimal, local earthquake magnitude expected to trigger rockfall is 4 (Keefer, 1984), which suggests that none of the globally mapped rockfalls occurred during the Apollo-era. However, we note that magnitude is not a good measure for local rockfall susceptibility, as the key factors controlling seismicity-driven rockfall are local peak-ground acceleration, earthquake duration, local topography, and local geology, such as degree of fracturing and weathering (Keefer, 1993;Massey et al., 2014;Newmark, 1965) meaning that the correlation of quake epicenters location and magnitude with rockfall hotspots is not ultimately conclusive. It is important to note that the majority of observable rockfalls and tracks (which likely survive erosion over millions of years, e.g., Hurwitz & Kring, 2016) have likely been produced outside of the small Apollo-era seismic monitoring time window (∼9 yr), meaning that a comparison of rockfall and moonquake locations is not particularly conclusive. In addition, the recorded shallow seismic activity probably does not represent the true lunar endogenic activity as it occurred over the past few millions of years, which likely features events with return periods much larger than 9 yr, that is, events with potentially significantly increased magnitudes. In addition, endogenic activity might originate from blind thrust faults that do not have a geomorphic expression on the surface.
Acknowledging the limitations of this analysis, mostly governed by the available data, our results indicate that lunar tectonic activity as recorded during the Apollo era has not been a main, global-scale driver of rockfalls in the Moon's recent geologic past. Our findings are reflected in the observation that the number and spatial concentration of rockfalls in the proximity of geomorphically visible tectonic features is very 10.1029/2021JE006824 20 of 26 low: only 0.8% of all mapped lunar rockfalls are close to tectonic features such as graben (Bickel, Jordan, et al., 2020). Thus, while we cannot definitively exclude the possibility that lunar seismic activity was a dominant, global driver of rockfalls in the lunar geologic past and that lunar seismicity still triggers rockfalls on local or even regional scales (see, e.g., Bickel, Jordan, et al., 2020;Mohanty et al., 2020), this possibility is not supported by the analysis presented here. Future, ground-based, long-term geologic investigations and global geophysical monitoring networks, such as for example, promoted by Jawin et al. (2018) and Nunn et al., (2021), will help to address many of the remaining questions. Our results suggest that observations of mass wasting features, such as rockfalls, on planetary surfaces may not always be a direct proxy for recent seismic activity or even seismic hazard, but can have other, potentially more likely drivers.
We observe that regions with large radial tidal displacement values can feature an increased number of rockfalls ( Figure 4). Currently, tidal forces can only be mapped on large spatial scales (e.g., Thor et al., 2021) and we are not able to observe the same trend in all studied geographic regions, meaning that our observations are not entirely conclusive. Tidally driven, cyclic, and vertical displacements of the lunar crust on the order of tens of centimeters could, in theory, favor fracture growth, weathering, and erosion in general, resulting in more rockfall events on average. On Earth, tides have been observed to trigger shallow thrust fault earthquakes (Cochran et al., 2004;Metivier et al., 2009). Tidal displacement could also interact with the near-surface stress state of the Moon (Watters et al., 2019), potentially helping to explain the observed at-target rockfall slope aspect rotation from N-S to E-W at higher latitudes ( Figure 3).
All observations combined indicate that short-term meteoritic impacts, aided by long-term temperature cycles, could have been the dominant drivers of rockfall occurrence in the Moon's recent geologic past. The distinct daily temperature cycle-predominantly controlled by the Sun-might interact with impact-generated fracture networks to drive fracture growth over long timescales (e.g., Molaro et al., 2017), which could eventually lead to kinematically free boulders and their subsequent detachment and failure in rocky and steep source regions (types (1) and (2)), as hypothesized by Hovland and Mitchell (1973) in the Apollo-era. This means that thermal fatigue might be a major long-term, global-scale causal driver of lunar erosion and rockfall. Short-term thermomechanical triggering of rockfalls, as observed on Earth during extremely hot days (e.g., Collins & Stock, 2016), is less probable on the Moon, however, but might occur in steep, exposed types (1) and (2) source regions. Our observations further indicate that in more regolith-supported, debris-type source regions (granular soil slopes, types (3) and (4)) with shallower slopes angles, predominantly present in older terranes, meteoritic impacts act as the dominant short-term causal drivers. Small-scale impacts also appear to be the dominant short-term driver of rockfalls in the polar regions (>85°N and S), as previously observed by e.g., Bickel and Kring (2020). Rockmass degradation and boulder initiation in all source regions could potentially be facilitated by continuous, abrasive meteoritic micro-bombardment. We summarize our observations in a conceptual model, depicted in Figure 12, illustrating how a new impact crater is formed, how thermal fatigue drives wall erosion on equator-facing slopes in the time period after Figure 12. Conceptual model of long-and short term lunar rockfall drivers, modified from Bickel, Jordan, et al. (2020). In young terranes with cliff-and outcrop-type source regions (type (1) and (2)), thermal fatigue can act as dominant long-and possibly short-term driver of rockfall, in addition to impacts and potentially moonquakes. In old terranes with boulder field-and plain slope-type source regions (type (3) and (4)) thermal fatigue appears to become a less important driver, while small-scale impacts appear to play an increasingly important role in driving rockfall occurrence. We note that in reality all drivers can be active in young and old terranes, but we only show the dominant drivers per panel to maintain clarity.
the impact, how impact-induced and -ejected rockfalls become the dominant short-term drivers of rockfall as crater age increases, and how long-term seismicity might contribute both to long-term damage and shortterm rockfall triggering in the vicinity of visible or invisible tectonic features.
On Mars, an atmospheric body unlike the Moon, Tesson et al. (2020) also observed that rockfall occurrence is likely connected to thermal stress induced by solar insolation, as the rockfalls in their study area tend to cluster on N-S-facing slopes as well. Using their observations, Tesson et al. (2020) argue that thermal fatigue plays an important role in causing rockfall events on Mars. Studies by Eppes et al. (2015) provides an additional evidence that sun-induced thermal stress is a dominant driver of fracture propagation in Martian boulders. Similarly, Molaro et al. (2017Molaro et al. ( , 2019 observe that directional insulation controls rock breakdown on airless planetary bodies such as asteroid Bennu. Despite the geologic and environmental differences between the Moon, Mars, and other (airless) planetary bodies, the agreement between the different and independent sets of observations increases our confidence in the presented results.
We note that this study has several important limitations. The performance of the CNN used by Bickel, Jordan, et al. (2020) to map the global population of lunar rockfalls might be limited; for example, a bias could be introduced by a non-direction invariant detector that might miss rockfalls of a certain orientation, even if such biases were not observed in the test set (Bickel et al., 2018). Further, the Klee-based (Klee, 1977) image selection algorithm (Bickel, Jordan, et al., 2020) might select inappropriate imagery; for example, an image with a solar incidence angle of 60° (sunrise) would result in a shadowed west-facing slope of an equatorial crater with walls steeper than 30°. We applied the following logic when selecting the incidence angle range from 10 to 60° for this work: (a) the vast majority of the lunar surface, including craters, features slope angles below 30° (Kreslavsky & Head, 2016); (b) an incidence angle of 60° allows for the analysis of the surface up to 60°N and S; (c) the number of images available for the analysis is generally larger, allowing for a nearly seamless coverage (i.e., there are no coverage gaps).
Locally, there are slopes that are steeper than 30°, meaning that there is a small chance that rockfalls would be missed, but only if: (a) the image selection algorithm happens to select an image with a very high incidence angle and (b) if the respective slope is facing away from the Sun. Also, we point out that rockfalls are usually deposited at or beyond the bottom of slopes, meaning that they might not be shadowed even if the steepest portion of the slope is slightly steeper than 30°. We further note that the rockfalls that were mapped manually for this study (i.e., should be illumination bias free; Bickel, 2021) deposited on a ∼19° slope angle, on average, with a maximum of ∼32°-only 1% of all manually mapped rockfalls deposited on slopes steeper than 30°. We point out that the detection of rockfalls at higher latitudes is more likely to be biased by the illumination conditions, as the solar incidence steadily increases and thus the portion of shadowed regions and images with higher incidence angles. As more NAC images become available over time, future work might be able to address some of the remaining, illumination-related uncertainties, for example, by applying an incidence angle range from 10 to 50°.
Another limitation of our global analysis is that the used CNN only detects deposited rockfalls. Thus, the globally derived at-target distributions represent the physical properties of the rockfall deposition regions, not necessarily the source regions. Considering that the measured (across 13 AoIs), average track length of a lunar rockfall is ∼600 m, the effect of this on the comparison with low-resolution data sets (>2 km/ pixel, e.g., GRAIL gravity data) is probably minor, but could affect the comparison with high-resolution data sets (∼250 m/pixel), specifically the temperature-related data sets. Data sets like slope aspect are more robust versus that source-deposition uncertainty, as very low slope angles (<10°) still produce the same aspect value as high slope angles (>30°), meaning that the attributes derived from rockfall source and deposition regions are more likely to be very similar or identical. In some cases the separation of deposition and source region properties is desired, however, for example, when calculating the elevation difference or Fahrböschung; here, the relatively low spatial resolution of some auxiliary data sets can also introduce uncertainties, for example, in cases where the rockfall track length is shorter than 1 pixel of an auxiliary map (note that the spatial resolution of the used DEM is ∼118 m/pixel). We point out that the generally low spatial resolution of all available auxiliary data sets complicates the interpretation of our global-scale results: rockfall-deposited boulders and their source regions are meters in size, not hundreds of meters or kilometers (see e.g., Figure 10). We additionally point out that substantial portions of the used LOLA DEM are derived by interpolation, potentially introducing inaccuracies or errors that could affect all slope angle and 10.1029/2021JE006824 22 of 26 aspect-related analyses, particularly for very shallow slopes. We note that slope aspect values derived from pixels/samples with slope angles below 5° have been omitted from the figures and analysis to minimize a potential bias of the used LOLA DEM. The current set of observations and interpretations will significantly benefit from future, high resolution, and high-quality data sets.
Finally, the qualitative and quantitative assessment of the relation between rockfall abundance and seismic activity is potentially biased, too, as the epicenters of all Apollo-era moonquakes are only spatial approximations due to the small spatial extent of the Apollo-era geophysical network, meaning that Nakamura et al.'s (1979)'s and Watters et al. (2019)'s catalogs of events contain a localization error that could bias our analysis in unforeseen ways. Our analysis is additionally impacted by the short Apollo-era monitoring period that likely introduces a spatial, temporal, and event scale-related bias. As for the rest of the analysis, future data will help to address these uncertainties and limitations.

Transport Mechanisms of Lunar Rockfall
After detachment, boulders usually start to roll and speed up, before they potentially start to bounce. Once they reach the toe of the slope, they slow down and stop. Our (locally constrained) measurements show that rockfalls originate from source regions with average angles of ∼25° (in our AoIs) and are deposited on slopes of around ∼13° (globally, Bickel, Jordan, et al., 2020), on average, across all geomorphic contexts (see Figure S8 in Supporting Information S1). Rockfalls that originate on steeper slopes achieve longer runout lengths, as one would expect. Interestingly, there are rockfalls in all AoIs that stop on inclines as steep as ∼30°. These boulders are not anomalous with regard to their shape or trigger; apparently, a very small reduction of the local slope angle and/or collisions with micro-topography or obstacles which are not resolved by the NAC camera are sufficient to force the respective boulders to a halt. In comparison, on Earth rockfall trajectories are strongly influenced by the roughness of the slope and the size and character of obstacles present (e.g., Dorren, 2003). Again, impact-ejected rockfalls are excluded from this comparison, as they occur on the full range of slope angles present on the lunar surface. Impact-ejected rockfalls are initially airborne and keep bouncing until they lose enough energy to enter a rolling or sliding displacement mode. Most impact-ejected rockfall trajectories have a distinct first bounce mark, usually a particularly deep and wide circular depression that likely is the result of an initial impact from high altitudes (impact ejecta).
We additionally observe that: (a) bouncing boulders with (b) higher approximated sphericity values achieve the longest runout distances. As non-impact-ejected boulders only bounce on slopes steeper than ∼20°, the slope geometry, combined with a higher likelihood of rolling kinematics (indicated by higher sphericity values), could be the dominant factor influencing lunar rockfall transport kinematics. Interestingly, and in contrast to Earth, we are not able to identify any relation between boulder volume and runout length. We note that we use a simplified characterization of boulder shape and volume, based on only three orthogonal axes (length, width, and height) that assumes an ellipsoidal shape, which may omit potential relations between runout and boulder properties. On Earth, boulder micro-topography has been identified as one of the main aspects that control rockfall trajectories (e.g., Dorren, 2003); unfortunately, the data available for the Moon (space-borne sensors) does not match the quality and amount of data available for Earth (groundbased and air-borne sensors) which means that a more detailed reconstruction of lunar boulder shapes will only become possible in the future.
We observe a number of rockfalls that do not follow a steepest descent path, specifically during the first tens of meters of their displacement: these boulders initially deviate from the direction of steepest descent before eventually returning to a steepest descent path further down the slope (see Figure S9 in Supporting Information S1). This means that these boulders have been triggered by a (highly) energetic trigger which provided a lateral push. The observations in our AoIs indicate that rockfalls that deviate from the steepest decent direction are exclusively associated with small-scale impact features, suggesting that these boulders either are impact-induced but with a lateral component, or impact-ejected but with very low ejection velocities (given their spatial proximity to the impact itself). Other than impact-ejected rockfalls, we do not observe any boulders that start bouncing directly after being triggered, indicating that neither impact processes (impact-induced type) or seismic shaking appear to have a sufficient amount of energy to get lunar boulders off the ground directly.
We note that rockfalls with longer runouts tend to have relatively more vertical drop than rockfalls with shorter runouts, apparently independent of the geomorphic context (see Figure 11). In contrast, the rockfalls in AoIs with a tectonic geomorphic context, Vallis Alpes and Rima Ariadaeus, feature a very constant, quasi-linear distribution of Fahrböschung angles (Figure 11 and highlighted in Figure S12 in Supporting Information S1). In other words, rockfalls in tectonic regions travel farther per meter dropped when compared to the other geomorphic regions. This is surprising, as source region slope angles, approximated boulder sphericity, and boulder volumes are very similar across all AoIs. The anomalous mean Fahrböschung angles could be the result of a sampling bias (too few AoIs per context were studied) or could represent differences in the absolute relief across the AoIs, where tectonic features would feature smaller absolute elevation differences than, for example, impact craters. The difference in mean Fahrböschung could also be related to regional or local variations of slope characteristics (e.g., roughness and restitution coefficients) that are not obvious in the available imagery. However, additional mapping efforts in various AoIs are required in order to exclude observational bias and to substantiate these initial observations.
Ultimately, many of the open questions could be addressed by closely monitoring a presumably active lunar slope with a ground-based geophysical monitoring system which not only includes seismometers but also optical and thermal instruments that allow for the cross-correlation of rockfall occurrence, visual changes, and thermal fluctuations with high spatial and temporal resolution.

Conclusions
We performed the first comprehensive and global study of lunar rockfall long-and short-term causal drivers and transport mechanisms using data about more than 130,000 rockfalls in combination with highly detailed local data about 687 rockfalls in 13 focus regions, covering volcanic, tectonic, and impact geomorphic contexts. Lunar rockfalls appear to predominantly occur on the slopes of N-S-facing, rocky impact structures that experience large thermal amplitudes over the course of the lunar day. We do not observe a qualitatively or statistically relevant relation between rockfall abundance, recorded Apollo-era shallow moonquake activity, and the distribution of visible tectogenetic features, indicating that moonquakes have not been a main, global driver of lunar rockfalls in the Moon's recent geologic past.
Our observations suggest that rockfalls originate from four source region types, specifically: cliff-, outcrop-(rock-supported), boulder field-, and plain slope-type (regolith-supported) source regions, where cliff and outcrop type source regions are more common in young, Copernican terranes and produce larger numbers of rockfalls. The transportation length of boulders appears to depend on their sphericity and the source region slope angle, but not on their volume. We note that the displacement kinematic of regular lunar rockfalls appears to be controlled by the inclination of the slope, specifically, where boulders potentially start to bounce above ∼20°.
Our results suggest that, in addition to impacts, continuous, solar-induced fracture propagation (thermal fatigue) is an important and global long-term driver of lunar rockfall mainly in young terranes with steep, exposed bedrock, potentially supported by continuous and abrasive micrometeorite bombardment and long-term tectonic and seismic activity. Erratic, small-scale impacts appear to be an important and global short-term driver of lunar rockfall (more pronounced in old terranes with shallow slopes and fewer bedrock outcrops). We note that shallow, seismic events could act as short-term drivers as well, particularly in the vicinity of tectonic features such as grabens, lobate scarps, and wrinkle ridges. We find that the trajectories of lunar rockfalls are mainly controlled by the trigger energy and the geometry of the slope.
Our findings can be applied to study the endo-and exogenic activity of the Moon, while complementing existing and ongoing studies of the weathering and erosion of atmospheric and airless planetary bodies in general, including Earth, Mars, and small bodies such as asteroid Bennu. Our findings can further inform and support the planning of future, geophysical surface exploration missions to the Moon.