Studying the Global Spatial Randomness of Impact Craters on Mercury, Venus, and the Moon With Geodesic Neighborhood Relationships

Impact crater records on planetary surfaces are often analyzed for their spatial randomness. Generalized approaches such as the mean second closest neighbor distance (M2CND) and standard deviation of adjacent area (SDAA) are available via a software tool but do not take the influence of the planetary curvature into account in the current implementation. As a result, the measurements are affected by map distortion effects and can lead to wrong interpretations. This is particularly critical for investigations of global data sets as the level of distortion typically increases with increasing distance from the map projection center. Therefore, we present geodesic solutions to the M2CND and SDAA statistics that can be implemented in future software tools. We apply the improved methods to conduct spatial randomness analyses on global crater data sets on Mercury, Venus, and the Moon and compare the results to known crater population variations and surface evolution scenarios. On Mercury, we find that the emplacement of smooth plain deposits strongly contributed to a global clustering of craters and that a random distribution of Mercury's basins is not rejected. On Venus, the randomness analyses show that craters are largely randomly distributed across all sizes but where local nonrandom distributions due to lower crater densities in regions of recent volcanic activity may appear. On the Moon, the global clustering of craters is more pronounced than on Mercury due to mare volcanism and the Orientale impact event. Furthermore, a random distribution of lunar basins is not rejected.

Often, the spatial randomness of crater populations is quantified by measuring the spatial relationships between craters. Since craters are located on a curved planetary surface, the determination of such relationships requires the consideration of the planetary curvature, particularly when investigating large or global surface units. This is typically considered in individual approaches that are developed for global applications (e.g., Turcotte et al., 1999). However, two commonly applied methods that are implemented in a software tool do not consider such effects. To overcome these limitations, we use geodesic measurements that determine distances and point coordinates on a great circle and include these measurements in the referred methods. We apply the improved approaches to investigate the global spatial randomness of crater populations on Mercury, Venus, and the Moon and use the results to draw parallels to previous investigations of the surface evolution on the respective planetary bodies.

Quantifying the Spatial Randomness of Impact Craters
The spatial randomness of impact craters is often analyzed using Monte Carlo approaches (e.g., Hirata et al., 2020;Kreslavsky, 2007;Kirchoff, 2017;Michael et al., 2012). In such approaches, a measure that describes the population of a given set of craters is compared to the same measure of n randomly distributed crater populations, each with the same number of craters. Depending on how much the measure of the examined population deviates from the measures of the randomly distributed populations, a statement on how strongly a given set of craters differs from a random distribution can be made. The deviation from a random distribution can be quantified using statistical measures such as percentile or Z-score (e.g., Kirchoff, 2017;Michael et al., 2012;Squyres et al., 1997). The crater populations itself can in turn be described by measurements that quantify the neighborhood relationships between the craters (e.g., Kreslavsky, 2007;Michael et al., 2012;Squyres et al., 1997). This can involve the distance or the area between neighboring craters, for example.

M2CND and SDAA Statistics in the Craterstats Software
Due to the widespread use of the Craterstats software (Michael & Neukum, 2010), the mean second closest neighbor distance (M2CND) and the standard deviation of adjacent area (SDAA) (Michael et al., 2012) are widely used neighborhood relationships (e.g., Adeli et al., 2016;Hao et al., 2020;Iqbal et al., 2019;Neesemann et al., 2019) that are used to quantify the spatial randomness of impact craters. In the M2CND approach, the distance to the second closest crater centroid is determined for each crater centroid. The M2CND value results from the mean of these values. The SDAA approach is based on a Voronoi diagram that is constructed from the crater centroids. The Voronoi diagram consists of polygons in which the distance of any point to the associated crater centroid is less or the same as to any other crater centroid. The SDAA value results from the standard deviation of the area of all Voronoi polygons. In order to assess the spatial randomness of a given crater data set, the obtained measures are compared to those of randomly distributed crater data sets.
The randomness analysis in Craterstats uses two statistical measures to quantify the deviation from a random population: Percentile and Z-score (the latter is termed  n in the Craterstats software). The percentile marks how many percent of the randomly distributed data sets yielded a lower randomness measure than the given population. The Z-score shows the deviation from the histogram in standard deviations from the mean. A low measure in the M2CND approach indicates a clustered population (because the mean distance of second closest neighbors is smaller than for randomly distributed populations); a high measure indicates an ordered population (because the mean distance of second closest neighbors is larger than for randomly distributed populations). In the SDAA statistics, this is reversed. Here, a low measure marks an ordered population (because there is less variance in the size of Voronoi polygons compared to randomly distributed population) and a high measure marks a clustered population (because the variation in Voronoi polygon sizes is larger than for randomly distributed populations).
To illustrate this, we generated a clustered set of 40 same-sized craters and calculated M2CND and SDAA statistics using Craterstats. In this example, the randomness measure of the data set is compared to the randomness measures of 1,000 sets of randomly distributed craters of the same quantity. The results are shown in Figure 1. In the M2CND statistics (Figure 1a), the M2CND value of the data set is located outside the lower end of the histogram; at the 0th percentile, with a Z-score of less than −3 (more than  3 below the mean). In the SDAA statistics (Figure 1b), the SDAA value is located at the upper end of the histogram, with a Z-score of 4.814. 99.9% of the randomly distributed crater data sets yielded a lower SDAA value than our data set. The statistical measures in both the M2CND and the SDAA statistics strongly indicate the presence of crater clusters in the given data set.

M2CND and SDAA Statistics From Geodesic Measurements
In Craterstats, all measurements as well as the construction of the Voronoi diagram are carried out in a two-dimensional Cartesian coordinate system. Accordingly, the calculated distances and areas are susceptible to map distortion effects. This affects global measurements in particular, because the distortions typically intensify with increasing distance from the map projection center (e.g., Kneissl et al, 2011;Snyder, 1987). The randomness analysis in Craterstats for example, uses measurements in the Lambert Azimuthal Equal Area (LAEA) projection, which correctly maps area sizes on a sphere, but distorts distances and shapes with increasing distance from the projection center (Snyder, 1987). We therefore apply great circle measurements to account for the curvature of planetary surfaces directly and thus circumvent the limitations of projected measurements when determining M2CND and SDAA statistics. The measurements are applied in order to (1) measure distances between craters to identify second closest neighbors for the M2CND approach and (2) determine the geodesic boundaries of Voronoi polygons for the SDAA approach.

M2CND
We measure the geodesic distances to all other craters for each crater in the given population, determine the neighbor with the second lowest distance for each crater and determine the M2CND value. Because a crater must have at least two adjacent craters, a minimum of three craters is required in a data set for this measurement.

SDAA
The calculation of global SDAA statistics requires the generation of geodesic Voronoi diagrams. For this purpose, we use the SphericalVoronoi algorithm (version 0.18.0) which is implemented in Python's scipy library (version 1.2.1) (Virtanen et al., 2020). The algorithm returns the vertices of spherical Voronoi polygons but does not generate geodesic polygon edges in the given version. In order to produce geodesic polygon edges and to construct geodesic Voronoi polygons, we use the coordinates of the polygon vertices generated by the SphericalVoronoi algorithm (we refer to these vertices as preliminary polygon vertices) and add further polygon vertices between them along a geodesic line. This ensures that the planetary curvature RIEDEL ET AL.  is taken into account when constructing the edges of the geodesic Voronoi polygons. We calculate the geodesic distance s between two preliminary vertices n P and 1 n P . If s is larger than 15 km, we calculate the coordinates of additional vertices at fixed intervals (we use intervals of 15 km for our investigation) along a geodesic line between n P and P n+1 . The geodesic Voronoi polygon is eventually constructed from both preliminary and additional vertices (see Figures 2a-2c).
If a geodesic Voronoi polygon intersects the date line or one of the poles, the construction of the polygon is slightly modified to ensure that such polygons are correctly generated. If a Voronoi polygon intersects the date line, the intersection is used as an additional vertex. In such a case, the date line acts as a cutting line, which separates the polygon into two parts. Each part is constructed individually to avoid errors when creating geodesic Voronoi polygons that cross the date line ( Figure 2d). In the case of a polar intersection, we use the intersection with the date line and the pole as additional vertices and modify the order of the vertices to generate the polygon. The pole and the date line intersection form the first two and the last two vertices of the polar Voronoi polygon. All intermediate vertices are added from west to east with increasing longitude values. The polygon with its first and last vertex at the pole is thus drawn in counter-clockwise direction around the pole (Figure 2e).
When a geodesic Voronoi polygon is constructed, it is projected to the LAEA projection (centered at the corresponding crater centroid) to measure its area. The areas of all geodesic Voronoi polygons allow the determination of the SDAA value. However, a current limitation in the SphericalVoronoi algorithm allows us to only apply the approach when there is a minimum of 20 craters in the data set. A smaller number of craters can lead to incorrectly calculated polygon vertices (see Supporting Information S1 and S2).
RIEDEL ET AL.

Projected Measurements Falsify the Results of M2CND and SDAA Statistics
To illustrate how projected measurements could falsify the results of the randomness analysis, we calculate geodesic as well as Cartesian M2CND and SDAA statistics for a global data set of 100 randomly arranged craters on a sphere with lunar dimensions and compare the results. We applied the traditional crater counting technique in CSFD Tools  and conducted a randomness analysis in Craterstats (Michael & Neukum, 2010;Michael et al., 2012) using Cartesian M2CND and SDAA measurements. The results are shown in Figure 3. Although the differences between Cartesian and geodesic measurements in this example do not allow a reliable rejection of the random distribution, the results show that Cartesian measurements distort global M2CND and SDAA statistics. However, because it cannot be ruled out that Cartesian measurements can also lead to an incorrect assessment of the truly existing spatial randomness of impact craters, we apply geodesic measurements to calculate global M2CND and SDAA statistics for craters on Mercury, Venus, and the Moon. We apply both methods to consider the different sensitivities of the two approaches to particular crater distributions (e.g., Kreslavsky, 2007;Michael et al., 2012). In contrast to the Craterstats software, we do not apply a binning based on crater size but divide a given crater data set into overlapping bins with the same number of craters. This has several advantages: (1) the overlap allows randomness variations to be identified in a higher resolution, (2) the constant number of craters allows for a better comparison of crater populations on different planetary bodies, and (3) the constant number ensures that the randomness analysis is based on a sufficient number of craters. In our analysis, each bin contains 300 craters, with the 150 smallest craters in the next larger bin corresponding to the 150 largest craters in the next smaller bin. Due to the high number of histogram and map plots that show the results for each bin, we placed the detailed results to the Supplementary Information. A few selected histogram and map plots are included in the respective chapters.
In this context, Fassett et al. (2012) and Orgel et al. (2020) analyzed the magnitude of basin asymmetry between Mercury's hemispheres. They found that the distribution of basins (  300 km D ) on Mercury's hemispheres is nonuniform and that the erasure of basins by resurfacing events may have contributed to it.
RIEDEL ET AL.  However, the studies by Fassett et al. (2012) and Orgel et al. (2020) only considered the distribution of basins over longitude. To provide more detailed information about the basin spatial relationships, we reinvestigate the global spatial randomness of impact craters and basins on Mercury by applying geodesic M2CND and SDAA analyses. To this end, we use two data sets. (1) The crater catalog by Fassett et al. (2011) (v19) to analyze the crater populations (   20.36 km 300 km D ) and (2) the basin data set by Orgel et al. (2020) to analyze Mercury's basin populations (  300 km D ) ( Figure 4). The latter data set was created from recent investigations of Mercury's basin inventory and contains 94 basins that are labeled as "certain," "probable," and "tentative," based on their visual detectability.

Results
Figure 5 summarizes the results of the randomness analysis at normalized Z-scores for each bin of the crater data set by . The results show that the number of bins where randomness is rejected at a two-sigma confidence level is higher in the SDAA than in the M2CND analysis. Out of 47 bins, randomness is rejected in 37 (SDAA) and 14 (M2CND) populations at the given confidence level. This indicates that both approaches have a different sensitivity for the crater configurations in the given bins. Crater populations where randomness cannot be rejected at a two-sigma confidence level are more numerous in the M2CND approach and appear across all diameters in both methods. There are eight bins where a nonrandom distribution of impact craters cannot be rejected from M2CND and SDAA analyses. They occur in bins with mean diameters of 20.8, 21.2, 23.6, 24.1, 24.7, 27.1, 27.7, and 38.5 km ( Figure 5 and Figures S3-S5  tion of craters appears less grouped over the surface, which leads to a nonrejection of randomness at a two-sigma confidence level in the corresponding randomness analysis. However, because the results of the randomness analysis largely reject randomness and because crater size and the significance of randomness rejection does not correlate in our analysis, we consider this a stochastic effect rather than an indication that resurfacing did not affect the distribution of craters in the respective bins. The largest coherent areas where craters appear in a less dense configuration are located around Caloris basin and, to a lesser extent, in the northern polar region. These regions are dominated by extensive smooth plains deposits, which erased preexisting crater records. In fact, many craters in this region that are included in the crater catalog are covered by volcanic material. We therefore conclude that the emplacement of smooth plains had a significant influence on the global spatial randomness of Mercury's impact craters, where differences in surface age and crater densities cause the rejection of randomness at the given confidence level in various bins. Further small-scale resurfacing processes such as tectonics (e.g., Watters et al., 2016), the geometric overlap of preexisting craters by a new impact (sometimes referred to as cookie cutting-e.g., Kneissl et al., 2016;Richardson, 2009;Woronow, 1977), degradation from downslope diffusion (sometimes RIEDEL ET AL. 10.1029/2020JE006693 9 of 21 referred to as sandblasting-e.g., Minton et al., 2019;Richardson, 2009;Soderblom, 1970), and the deposition of young smooth plains (e.g., Byrne et al., 2016) or volcanic material (e.g., Herrick et al., 2018;Thomas et al., 2014) likely contributed to the erasure of preexisting craters and the appearance of clusters in certain regions on Mercury. However, such processes occurred at a regional scale and therefore had a much smaller impact on the global clustering of impact craters. Figure 6 shows the results from the basin data set by Orgel et al. (2020). Orgel et al. (2020) grouped basins into categories of "certain," "certain and probable," and "certain, probable, and tentative" basins and analyzed whether there is a basin asymmetry between the hemispheres on Mercury. We adopt their classification to re-analyze the distribution of Mercury's basin population. Note however, that the number of basins is much lower than the number of craters that are used to analyze the cratering record with  300 km D .
Therefore, the results are not directly comparable to the spatial randomness of Mercury's cratering record.
In all basin categories, randomness is not rejected at two-sigma confidence by either approach, with percentiles ranging from 17% to 96% (M2CND) as well as from 24% to 91.1% (SDAA) along with Z-scores between −0.95 and 1.71 (M2CND) as well as −0.51 and 1.39 (SDAA). Although the map plots in Figure 6 indicate that the number of basins is larger in the western than in the eastern hemisphere, as pointed out by  (2020), the results do not imply that such an asymmetry corresponds to a strictly nonrandom distribution of basins. Thus, the large basin distribution on Mercury could nevertheless represent a random population.

Background
Compared to other planetary bodies in the inner Solar System, the number of craters on Venus is low (Figure 7) with most craters showing pristine topographic features (e.g., Kreslavsky et al., 2015;Phillips et al., 1992;Schaber et al., 1992;Strom et al., 1994). Despite the presence of various geologic units that are stratigraphically different (e.g., Basilevsky & Head, 1998Basilevsky et al., 1997;Ivanov & Head, 2011;Ivanov et al., 2015;Kreslavsky et al., 2015;Phillips et al., 1992), crater size-frequency distribution analyses have shown that there is very little variation in the spatial densities of impact craters between them (e.g., Kreslavsky et al., 2015;Schaber et al., 1992;Strom et al., 1994). This indicates that the observable surface of Venus is very young and that the visible cratering record accumulated over a recent period of probably a few hundred million years (e.g., Herrick & Rumpf, 2011;Kreslavsky et al., 2015;Phillips et al., 1992;Strom et al., 1994). Accordingly, the geologic history of Venus has strongly been influenced by resurfacing events that erased preexisting craters, mainly due to volcanic and tectonic events (e.g., Schaber et al., 1992).
There are two scenarios that describe the resurfacing history of Venus-equilibrium resurfacing and global resurfacing. In the equilibrium scenario (e.g., Bjonnes et al., 2012;Guest & Stofan, 1999;Phillips et al., 1992), volcanic and tectonic resurfacing occurred at a somewhat constant rate throughout the geologic history of Venus. This implies that the observable cratering record is related to a global crater equilibrium, where craters accumulate at the same rate they are erased. In the global resurfacing scenario (e.g., Airey et al., 2017;Basilevsky & Head, 2000Basilevsky et al., 1997;Ivanov & Head, 2011;Ivanov et al., 2015;Kreslavsky et al., 2015;Phillips et al., 1992;Price & Suppe, 1995;Schaber et al., 1992;Strom et al., 1994), an epoch of intensive volcanism and tectonic resurfacing eliminated all craters that predate the current surface units. This was followed by an epoch of decreasing and spatially limited resurfacing activity (see e.g., Ivanov & Head, 2011) during which the observable geologic record was formed, and impact craters accumulated. Accordingly, the observable crater record on Venus largely resembles the production population and has only experienced limited volcanic and tectonic modification (e.g., Price & Suppe, 1995;Schaber et al., 1992; RIEDEL ET AL.  Strom et al., 1994). Due to the presence of many well-preserved impact craters (e.g., Schaber et al., 1992;Strom et al., 1994), similar crater size-frequency distributions and crater retention ages on stratigraphically different units (e.g., Kreslavsky et al., 2015;Schaber et al., 1992;Strom et al., 1994), and investigations on rift evolutions (e.g., Airey et al., 2017), most recent publications favor the global resurfacing scenario in order to describe the resurfacing history of Venus. However, this scenario is still under investigation. For example, the observations by Hauck et al. (1998) and Bjonnes et al. (2012) suggest that the equilibrium scenario cannot be excluded and that, according to Grindrod et al. (2010), Rumpf (2011), andO'Rourke et al. (2014), the influence of volcanic and tectonic activity on crater modification may be underestimated.
In this context, the spatial randomness of impact craters was also examined and it was found that the global cratering record on Venus is randomly distributed across all crater sizes (e.g., Hauck et al., 1998;O'Rourke et al., 2014;Phillips et al., 1992;Strom et al., 1994;Turcotte et al., 1999). While Strom et al. (1994) used a chisquare test on the sine of latitude for this purpose, Phillips et al. (1992), Hauck et al. (1998), Turcotte et al. (1999), andO'Rourke et al. (2014) used a more sensitive approach based on the nearest neighbor relationships from angular distances. Here, we reinvestigate the spatial randomness of Venus' observable impact crater record based on geodesic M2CND and SDAA measurements. To this end, we use the USGS Venus Crater Database (see e.g., Schaber et al., 1992Schaber et al., , 1995Strom et al., 1994) with a total of 967 craters from 1.3 to 270 km in diameter.

Results
Figure 8 summarizes the results from the geodesic M2CND and SDAA analyses for each bin of the crater data set (for details, see Figure S13). It shows that in most cases, the crater populations cannot be distinguished RIEDEL ET AL.
10.1029/2020JE006693 12 of 21 however, is an exception. While the M2CND analysis does not reject randomness even at a one-sigma confidence level, randomness is rejected at two-sigma confidence in the SDAA analysis. Accordingly, the SDAA approach reacts more sensitively to the different densities of impact craters in this bin than the M2CND approach. The map plot in Figure 9 indicates that areas where craters with   15.19 km 33.09 km D are somewhat less abundant are irregularly distributed across the surface. Extended areas with lower crater densities occur for example along 140°W, between 60°N and 40°S and along 20°S, between 60°E and 180°E. Although not exclusively, these areas include Aphrodite Terra area, as well as the rift zones and lobate plains of eastern Aphrodite and the Beta-Atla-Themis region. While the overall low number of small craters on Venus may be related to atmospheric filtering (e.g., Schaber et al., 1992;Strom et al., 1994) and the rejection of randomness may have a stochastic origin, the rejection of randomness in the SDAA analysis could as well be caused by local crater degradation from volcanic and tectonic events in those regions. Although recent volcanic activity on Venus is less pronounced than on Earth (e.g., Strom et al., 1994), the Beta-Atla-Themis region has been identified as an area of recent volcanic activity (e.g., Airey et al., 2017;Schaber et al., 1992). It is therefore possible that craters have not only been modified, but also erased in these regions (e.g., Grindrod et al., 2010;Herrick & Rumpf, 2011;Strom et al., 1994), which in turn could result in a nonrandom distribution of craters with   15.19 km 33.09 km D . However, given that randomness is not rejected at two-sigma confidence in the remaining bins, this effect would be very low compared to other planetary bodies in the inner Solar System. This implies that ever since the observable cratering record accumulated on Venus, no resurfacing event was strong enough to cause significant nonrandomness of craters on a global scale. However, this does not exclude the possibility of local crater erasure by volcanic or tectonic events. This alone however, would not necessarily reject the possibility of an equilibrium resurfacing scenario. In order to erase craters at a somewhat constant rate and to keep the global crater population spatially random, the local resurfacing events must have been efficient at degrading craters but weak at resurfacing on a global level. Still, such a scenario seems unlikely since the presence of extended areas of recent volcanic activity would contradict a random distribution of craters on Venus (e.g., Kreslavsky, 1996;Romeo, 2013;Romeo & Turcotte, 2009). Therefore, the results from the geodesic M2CND and SDAA analysis are most consistent with the global resurfacing scenario, where craters accumulated in an epoch of decreasing and spatially limited volcanic and tectonic activity (e.g., in the Beta-Atla-Themis region). In this scenario, we see a crater population that is spatially random for the most part, but where local variations in spatial crater densities due to volcanic and tectonic activity can cause a nonrandom cratering record to some extent.

Background
Due to the absence of an atmosphere and limited surface modification processes during its evolution, the Moon has the best preserved cratering record in the inner Solar System (Figure 10). The lunar cratering record has been used to understand the bombardment history and the planetary surface evolution on various RIEDEL ET AL.
Intensive bombardment and accretion of basins and craters as well as lunar volcanism have led to the formation of two major geologic units-low albedo lunar maria with smooth terrain and high albedo lunar highlands with rough terrain (e.g., Wilhelms, 1987). The lunar highlands contain the oldest lunar surface units and cover large parts of the lunar surface. They are heavily cratered, with a high density of large craters. Crater retention ages show that some lunar highland units may have surface ages of >4.3 Ga (e.g., Neukum & Ivanov, 1994;Orgel et al., 2018). The lunar maria, on the other hand, are sparsely cratered with a lower density of large craters. They are a result of flood volcanism and consist of younger surface units, which formed in topographic lows, often within old basins (e.g., Head & Wilson, 1992;Wilhelms, 1987). Lunar maria cover roughly 17% of the lunar surface (e.g., Head & Wilson, 1992), the vast majority of which is located on the lunar nearside around the Procellarum KREEP terrane (PKT). On the farside, the lunar maria are sparsely distributed and cover much smaller areas, mostly within large craters and basins. The asymmetric distribution of lunar maria has been attributed to increased volcanic activity on the lunar nearside due to a thinner crust and a higher abundance of heat producing elements (e.g., Head & Wilson, 1992;Joliff et al., 2000;Miljković et al., 2013;Wieczorek et al., 2012;Zhu et al., 2019). This in turn, resulted in higher subsurface temperatures compared to the lunar farside and the formation of larger craters and basins due to differences in target properties (Miljković et al., 2013). Crater retention ages revealed that most nearside mare were emplaced at the late stage of lunar basin formation, during the late Imbrian period, around 3.3-3.7 Ga (Hiesinger et al., 2000(Hiesinger et al., , 2003(Hiesinger et al., , 2011. Volcanic activity decreased after the late Imbrian and came to an end around 1.2 Ga, when the youngest mare units in the center of the Procellarum KREEP terrane were emplaced (Hiesinger et al., 2003). Due to the decrease in cratering rate (e.g., Neukum & Ivanov, 1994), the lunar maria are in a much more sparsely cratered state than the lunar highlands, despite their relatively long exposure to impact cratering.
With the exception of mare emplacement, geologic resurfacing and crater erasure on the Moon has largely been controlled by impact cratering processes. Cratering-related effects that contribute to the erasure of the lunar cratering record include "cookie cutting" which occurs when cratering becomes nonsparse (e.g., Kneissl et al., 2016;Minton et al., 2019;Orgel et al., 2018;Povilaitis et al., 2018;Richardson, 2009;Riedel et al., 2018;Woronow, 1977), sandblasting which contributes to topographic diffusion and simple crater equilibrium (e.g., Minton et al., 2019;Soderblom, 1970), and burial by proximal ejecta blankets (e.g.,  Richardson, 2009). The intensity to which the individual processes influence the modification of the topography likely depends on the size of the impact . Accordingly, the lunar landscape and the visible cratering record is highly influenced by crater formation and erasure. Here, we apply geodesic M2CND and SDAA statistics to investigate the global spatial randomness of lunar impact craters with D ≥ 20 08 . km. To this end, we use the lunar crater catalog by Robbins (2019) with a total of 6,973 digitized craters larger than 20 km in diameter. Smaller craters, although included in the catalog, are not taken into account in the investigation. Due to the low number of basins, we investigate the lunar cratering (  300 km D ) and basin record (  300 km D ) separately.

Results
The results from the geodesic M2CND and SDAA analyses are summarized in Figure 11 (for details, see Figures S14-S22). It shows that randomness is rejected at two-sigma confidence for nearly all binned crater populations. The significance at which randomness is rejected is stronger in the SDAA analysis, which implies that the approach reacts more sensitively to the given crater configurations than the M2CND. There are two of 45 bins where randomness is not rejected at two-sigma confidence. In bins with  mean 29.5 km D and  mean 87.0 km D , 13% and 2.75% of randomly distributed data sets yield a lower M2CND value than the given population, resulting in Z-scores of −1.13 and −1.96, respectively. However, given that the remaining binned populations as well as the SDAA analysis reject randomness at two-sigma confidence, we consider this a stochastic effect rather than an indication for a randomly distributed population.
RIEDEL ET AL.

(b)
The results from the randomness analysis largely confirm previous investigations on lunar surface evolution and can therefore be used as an indicator that our geodesic modification of the M2CND and SDAA statistics deliver valid results. For example, the map plots (Figures S14-S22) indicate that throughout all sizes, there is a global clustering of impact craters in the binned crater populations due to differences in crater densities between ancient lunar highlands and younger surface units. These younger units involve the extended mare deposits on the lunar nearside as well as the area that surrounds Orientale basin, the youngest large basin on the Moon (e.g., Orgel et al., 2018). Accordingly, the lower density of impact craters together with a rejection of randomness confirms that a large part of preexisting craters were erased by the Orientale impact event and mare volcanism (e.g., Povilaitis et al., 2018;Whitten et al., 2011).
Mare volcanism in the PKT occurred during the Imbrian, Eratosthenian, and Copernican periods, with the youngest volcanism occurring in the central PKT region (e.g., Hiesinger et al., 2003). In contrast, the mare units surrounding the center of the PKT as well as the mare areas east of it were emplaced during an earlier episode of mare volcanism, mainly during the late Imbrian period (e.g., Hiesinger et al., 2000Hiesinger et al., , 2003Hiesinger et al., , 2006Hiesinger et al., , 2011. As a result of the different timing of mare emplacement, smaller craters with  32 km D could largely accumulate in mare regions east of the PKT and around large basins, such as Crisium and Nectaris, where mare volcanism ended early, whereas ongoing mare volcanism in the central PKT contributed to a longer lasting erasure of craters, which can be seen in a low abundance of craters in all bins (Figures S14-S22). Furthermore, crater populations within the PKT are not necessarily affected by a superposing basin impact. This suggests that mare emplacement in this area could have caused the erasure of large preexisting craters larger than 100 km in diameter. Buried craters of such size have been identified by Evans et al. (2016) and Sood et al. (2017) beneath lunar mare deposits. Accordingly, lunar mare deposits in this region would have to reach thicknesses of several kilometers in order to completely cover such craters. Although the depth of lunar mare deposits is poorly constrained, various studies found evidence for a similar thickness of lunar nearside maria (e.g., Gong et al., 2016;Thomson et al., 2009;Williams & Zuber, 1998).
In contrast to the cratering record (  300 km D ), randomness is not rejected at two-sigma confidence for lunar basins with  300 km D (Figure 12). The obtained M2CND value is larger than in 87.6% of the randomly distributed data sets (  1.14 Z ) and the obtained SDAA value is larger than 29.4% of the randomly distributed data sets (  0.58 Z ). However, since all basins are summarized in one bin, a detailed distinction between basin sizes is not possible. This implies that the nonrejection of randomness at two-sigma confidence does not necessarily contradict an asymmetry in target properties at the time of basin accretion, which influenced the final size of lunar near side and far side basins (Miljković et al., 2013). Such an investigation using M2CND and SDAA statistics would require a further subdivision of the basin population, which would involve fewer craters in the M2CND and SDAA analyses. Since involving fewer craters would yield in less representative results, M2CND and SDAA are not suitable in making such distinction from the lunar basin record. Furthermore, as the number of basins per bin is far less than 300, the results from the randomness analysis of lunar basins are not directly comparable to the spatial randomness of the lunar cratering record (  300 km D ).

Conclusions
In this study, we presented improved techniques to quantify the global spatial randomness of impact craters with respect to planetary curvature. The methods are applicable on all planetary bodies that are well approximated by a sphere. As the presented approaches are sensitive to the number of input craters, the planetary cratering records are subdivided into populations of the same quantity. Thus, it is possible to identify size-dependent variations from spatial randomness and to directly compare the spatial randomness of cratering records on different planetary bodies. The size of the subdivisions is a balance between a minimum number of craters at which variations from spatial randomness due to resurfacing events are detected and a maximum number of craters at which size-dependent variations in spatial randomness can still be identified. In our investigation, we subdivided the cratering record of the investigated bodies into bins of 300 craters to compare the results of Mercury, the Moon, and the cratering record of Venus. This number can certainly be adjusted. However, we recommend using 300 as a minimum number to properly identify nonrandom populations and to minimize stochastic effects when quantifying the spatial randomness of impact craters. We also recommend using both (or further) approaches to identify deviations from nonrandom populations since different methods can react at specific sensitivities to certain nonrandom distributions. In general, the SDAA approach showed a higher sensitivity toward nonrandom distributions in the given crater populations.
On Mercury, we identified a global clustering of impact craters with   20.36 km 300 km D due to the emplacement of extensive smooth plains deposits and the Caloris impact event. However, since the timing between the emplacement of intercrater plains and smooth plains is relatively close, the deviations from random populations are generally less significant than on the Moon. Particularly with the M2CND analysis, a random distribution of craters was not rejected at two-sigma confidence in various populations. This implies that, although major resurfacing events occurred on Mercury that cause the global cratering record with   20.36 km 300 km D to be in a clustered distribution, their influence on the global spatial randomness is less intense than on the Moon. Furthermore, the SDAA approach is more sensitive in recognizing such a resurfacing scenario. In contrast, a random distribution of basin-sized craters, was not rejected at two-sigma confidence, although previous studies indicated a basin asymmetry between Mercury's eastern and western hemisphere Orgel et al., 2020).
The randomness analyses show that craters on Venus are largely randomly distributed across all sizes. However, randomness was rejected at two-sigma confidence for the crater population with  mean 22.4 km D using the SDAA approach. Although it is not clear whether the rejection of randomness is due to a stochastic effect, some of the areas where the crater density of this population is lower correspond to regions of recent volcanic activity. We therefore suggest that local volcanic activity on Venus may have contributed to crater erasure, which in turn can result in a certain nonrandom distribution of impact craters. Our results also confirm that the influence of recent volcanism on crater erasure was not strong enough to cause a significant clustering of impact craters on a global scale across all crater sizes. The given configuration could in principle be explained by an equilibrium resurfacing scenario in which the given crater configurations were caused by volcanic processes that are highly efficient at erasing craters at a local scale, but weak at resurfacing on a global level. However, such a scenario seems less likely since the geologic units in which recent volcanic modification occurred are too large to maintain a random distribution of craters (e.g., Kreslavsky, 1996;Romeo, 2013;Romeo & Turcotte, 2009). Therefore, our results are more consistent with a global resurfacing scenario, in which craters on Venus accumulated in an epoch of decreasing geologic activity and where local crater erasure in areas of recent volcanic and tectonic activity may cause the rejection of randomness in certain crater populations.
The investigation of the lunar cratering record with   20.08 km 300 km D shows that a random distribution of craters is rejected at two-sigma confidence in almost all populations, with normalized Z-scores that are clearly below the values of Mercury's crater populations in both approaches. Therefore, the global clustering of craters is more pronounced on the Moon than on Mercury. The stronger rejection of randomly distributed crater populations on the Moon is most likely caused by a greater difference in surface age between the oldest and the youngest geologic units, which leads to greater differences in crater densities. The results confirm that mare volcanism, as well as the Orientale impact event, had a major impact on the erasure of the preexisting crater record. The results also suggest that mare emplacement in the central PKT likely caused the erasure of preexisting craters larger than 100 km in diameter. This indicates that mare deposits in this region can reach a thickness of several kilometers. In the analysis of the lunar basin record, a random distribution was not rejected at two-sigma confidence by either approach. However, this does not contradict a possible asymmetry in target properties during the time of basin emplacement (e.g., Miljković et al., 2013).
In summary, our improved methods together with a binning of crater populations allow for an accurate randomness analysis of global crater populations. The results are consistent with known population variations on Mercury, Venus, and the Moon and can be used to support a number of surface evolution scenarios on the respective planetary bodies. Therefore, we consider the presented approach to be a robust improvement to measure the spatial randomness of impact craters with respect to the planetary curvature.

Data Availability Statement
The associated data are available via Mendeley Data: Riedel et al. (2020).