Controls on the Global Distribution of Martian Landslides

Recent acquisition of high‐resolution satellite imagery of the Martian surface has permitted landslides to be studied on a global scale on Mars for the first time. We apply the Scoops3D software package to compute slope stability for select regions of the Martian surface, combining calculations of slope stability with number of observed landslides (Crosta, Frattini, et al., 2018; Crosta, De Blasio, et al., 2018), as reported in a recently published inventory of Martian landslides, to understand controls on the global distribution of landslides on Mars. We find that the distribution of landslides does not simply follow the distribution of unstable slopes. In particular, there is an abundance of landslides around Tharsis, and especially in Valles Marineris and Noctis Labyrinthus, which is not explained by an abundance of unstable topography alone. We analyzed for but did not find a clear large‐scale lithologic or stratigraphic control on landslide occurrence from subsurface heterogeneities. Other possibilities to explain the increased occurrence of landslides in Tharsis include (1) thin weak unit(s) that is regionally widespread and at multiple stratigraphic levels, such as from interbedded ashes; (2) seismic activity related to the Tharsis's geological activity, and (3) possible groundwater near Valles Marineris into the Amazonian. Given the apparently young ages of many landslide deposits in Valles Marineris (Quantin et al., 2004), continued modern day analysis of lithologies in Valles Marineris and observations of Martian seismicity may act to strengthen or rebut the first two hypotheses.

More recently, the proliferation of high-resolution imagery of the Martian surface generated by the CTX camera (Malin et al., 2007) has enabled landslides to be studied on a global scale. In particular, Crosta, Frattini, et al. (2018b) and Crosta, De Blasio, et al. (2018a) have characterized Martian landslides at a global scale for the first time with a published inventory of mapped Martian landslides (Figure 1). They map 3,118 landslides with deposit areas greater than 0.1 km 2 between 60° N and 60°S. We apply this full inventory of landslides, excluding only those features mapped as rock glaciers, in this study. The highest density of landslides is found near Mars's equator, and particularly in the Valles Marineris region. Morphologically, rock avalanche landslides dominate across most of Mars's surface, with slump landslides more prevalent in the outflow channels in the eastern part of the Valles Marineris system. An area's landslide susceptibility is influenced by numerous factors. Some of these, such as hill slope angle and length, are observable from orbit. Others, such as rock strength, the presence of groundwater, and the history of seismicity in an area, are difficult or impossible to observe or infer from orbit. The study of spatial distributions of landslides can thus provide insight into these factors, if the different mechanisms which can trigger landslides can be distinguished uniquely.
We use the new (Crosta, Frattini, et al., 2018b) global landslide inventory in concert with global Martian datasets, including topography and faulting patterns, to undertake a comparative regional-scale analysis of factors driving the spatial distribution of landslides on the Martian surface. We first investigate the hypothesis that the distribution of landslides simply follows the distribution of unstable slopes on the Martian surface, which means it is controlled by topographic characteristics alone. We then consider the possible influence of spatially varying tectonic seismicity on Mars, the possible influence of variable rock strength, and the possible influence of groundwater on interpretation of our modeling results.

Digital Elevation Model Dataset
For input to slope stability models we used elevation data from High Resolution Stereo Camera (HRSC) DEMs, downsampled via bilinear interpolation to a common resolution of 200 m/pixel to remove any potential influence of DEM resolution on the stability analysis. The average height error of 11.6 m associated with HRSC data has little impact on our factor of safety (FOS) calculations, as we are considering the susceptibility to large-scale landslides, big enough to be mapped from orbit by Crosta, Frattini, et al. (2018b) on slopes hundreds of meters to kilometers in height. We chose this HRSC DEM dataset because of its availability, its large regional coverage, and its demonstrated greater accuracy in the steepest topography in comparison to CTX DEMs created by us. We attempted to rapidly produce CTX DEMs for some of our study areas via Ames ROBACK AND EHLMANN 10.1029/2020JE006675 2 of 14 Stereo Pipeline terrain-processing tools (Beyer et al., 2018), but experienced sporadic large errors (on the order of 50-100 m) in excess of those associated with HRSC data in high topography portions of Valles Marineris. Due to our need to acquire topographic data across large portions of Mars's surface quickly, we chose the already-available HRSC dataset rather than working to further improve the CTX DEM quality of each stereo pair.

Calculating Susceptibility to Landsliding With Scoops3D
To quantify the relative susceptibility to landsliding of various topographies on Mars, we apply the Scoops3D software package developed by the USGS. (Reid et al., 2015). Scoops3D calculates the forces driving and resisting motion on a large number of possible spherical trial landslide surfaces generated underneath an input digital elevation model (DEM). These trial surfaces are spheres drawn around nodes placed at intervals set by the user above an input DEM; we set this interval to be 200 m. The number of trial surfaces generated per DEM was typically on the order of 10 6 , varying with the spatial extent of each DEM.
Since the trial landslide surfaces are spherical, the landslides are simulated as rotational failures. Though we cannot know definitively that past landslides on Mars happened as rotational failures, arcuate planes of failure are commonplace in materials without large contrasts or discontinuities in strength (Hoek & Bray, 1981), and anticipated by numerical modeling (Katz et al., 2014). Additionally, many of the largest landslides on Mars (e.g., Lucchitta et al., 1979) show curved source scars. However, in materials with strong discontinuities or planes of weakness, such as weak rock layers, this assumption may lead to a relative underestimation of the actual landslide susceptibility. The possible implications of such heterogeneity in subsurface material properties are discussed more fully in Section 4.2.2.
The model splits a trial failure mass into columns, whose position and size are set by the horizontal resolution of the DEM. Over each of these, driving and resisting forces are calculated and summed to derive an estimate of the ratio between the total forces resisting motion divided by the total forces driving motion of the trial landslide, known as the FOS (Figure 2). Higher values of the FOS, output pixel by pixel for each DEM cell, indicate greater stability. Simpler slope stability modeling techniques commonly involve application of 2D infinite-slope modeling techniques (e.g., Cernica, 1995;Sharma, 1994) to individual cells of a DEM on a pixel-by-pixel basis; these approaches assume slopes are infinite in extent, and neglect effects such as added resistance at the toe of a slope, and the differing depths of failure planes by slopes of different lengths. Scoops3D simulates trial landslide failure surfaces in 3D, integrating the effects of topographic wavelength (the length of slopes) on the calculation of the stability of trial landslide surfaces, and the determination of the most-probable landslide size and depth associated with the least-stable trial failure surface. Commonly, slopes with averaged large-scale angles of ∼25° receive are assigned FOS values < 2, while slopes with ∼20° average gradients get FOS = 2-3.
To calculate the FOS, we input into Scoops3D a constant set of material parameters, including cohesion, angle of internal friction, and weight, for all study areas to efficiently isolate the impact of topography alone on the relative stability of different Martian landscapes. Parameters were set to values corresponding to a low-end intact rock strength, such as that it might correspond to a terrestrial shale. In particular, the cohesion was set to 1 MPa, the angle of internal friction to 30°, and the weight of material to 10 kN/m 3 . Typical cohesions for terrestrial rocks range from ∼1-100 MPa; angles of internal friction can vary from ∼15 to 50° (Goodman, 1980). For Martian rocks, strengths have been estimated both at the local scale from rover-based measurements (Okubo, 2007;Thomson et al., 2013) and at the large scale from strength modeling of Valles Marineris hillslopes (Schultz, 2002). At the rover scale, local rock outcrops analyzed at Eagle and Erebus craters by the Mars Exploration Rover (MER) Opportunity had unconfined compressive strengths of ∼10 MPa, corresponding to cohesions of ∼5 MPa (Okubo, 2007)  some weaker material (cohesions of 0.5-25 MPa). At the large scale, other slope stability analyses (Clow et al., 1988;Schultz, 2002) have estimated bulk strength of wallrock of 0.1-0.4 MPa, suggesting that the large-scale strength of rock may be set by fracture patterns.
Note that the weight of material is set lower than typical weights of terrestrial rocks and soils (20-30 kN/m 3 ) due to Mars's lower gravity (∼38% of Earth's). The permissible volume range of trial landslide failure surfaces was set to 10 8 -10 12 m 3 ; the high bound on possible landslide volume is set given the scale of landslide deposits observed in Valles Marineris (Lucchitta, 1979). A Bishop's simplified limit-equilibrium method was used (Hungr, 1987). This method assumes horizontal side forces, computed on the basis of the column's weight and the angle of the sliding surface below the center of each column, between the vertical columns comprising the modeled landslide failure mass. Earthquake loading and the presence of groundwater were not incorporated into the model setup; the point of this analysis is to consider the impact of topography alone (see "A role for groundwater?" in the Discussion section). However, as stated previously the FOS in this modeling scheme is a function of more than the local slope angle.

Statistical Analyses of Landslide Occurrence
The distribution of landslides was then compared to the distribution of unstable slopes, that is, slopes below a given FOS, via means of a window-based analysis. Analyzed parts of the Martian surface (see Figure 4, and link to our publicly archived slope stability data at the end of the paper) were divided into 20-km square windows in equal-area projections; the 20-km length-scale was chosen to balance minimizing the occurrences of large landslides spanning multiple terrain windows, along with generating a large population of windows to differentiate between terrains of different stability levels. Within each window, the number of landslides was tallied, along with the fraction of the topography with a calculated FOS value below a certain threshold deemed "unstable." For windows covered incompletely with slope stability data (i.e., on the edge of our HRSC DEM coverage), we calculated the percentage of unstable slopes with the available DEM data. We varied the threshold FOS to evaluate sensitivity to the choice, using values ranging from 3 to 6. In all plots presented here, we use 4 as the "unstable" threshold FOS value, unless otherwise noted, as the results were not very sensitive to the choice. The range of values we explored for the FOS threshold is informed partially by the range of FOS values computed for areas of observed, recent landslides on Mars; a histogram of FOS values acquired in the vicinity of mapped Martian landslides is provided in Figure 3. Landslides go from being extremely abundant in topographies at FOS < 3 to extremely rare at FOS = 6. Of the ∼2000 landslides that lie within areas covered by our FOS modeling, >95% of the landslides occur in areas with a modeled FOS < 4, hence the choice of 4 as the threshold value for instability value focused on in these analyses of % unstable topography. This pattern of FOS values holds whether we consider the FOS values computed at the centroid of the landslide's source area (  Windows were then binned by their fractional percentage of unstable slopes, and landslide populations of bins with the same percentage of unstable slopes were compared across different geographic regions of the Martian surface. Studied and compared regions include Valles Marineris, Elysium Mons, Olympus Mons, Kasei Valles, Nili Fossae, Libya Montes, Noctis Labyrinthus, and a general region for all other studied terrain, which is mostly mapped as Noachian Highlands (Tanaka et al., 2014) ( Figure 4). Across our studied areas, we delineated more than 50,000 terrain windows, covering an area of Mars's surface in excess of 20 million km 2 .
The statistical significance of differences in landslide occurrence between bins of different regions was evaluated via the use of a Wilcoxon rank-sum test on the population of windows within each bin. All windows were included in the Wilcoxon rank-sum comparison, whether they did or did not have landslides. The Wilcoxon rank-sum test was chosen because the populations of landslides in each bin were generally distributed in a non-Gaussian way, with most windows hosting 0-1 landslides and a typically small number of windows hosting multiple landslides. If the output p-value from the rank-sum test was < 0.05, indicating a >95% confidence that the difference in landslide abundance between the two bins was not a result of random sampling from the same distribution of landslides, we consider the difference in landslide abundance between the bins to be statistically significant. Further, we only consider differences in landslide abundance to be significant if p < 0.05 is satisfied for different choices of the "unstable" threshold FOS value.
The Wilcoxon rank sum test often, but not always, supports statistically significant differences between pairs of bins in which the average number of landslides per window is greatly different. Factors which can lead to p > 0.05 (meaning that the bins are not deemed significantly different) for pairs of bins with different overall abundances of landslides include small populations of windows in one or both bins as well as the presence of outlier windows, for example, when the number of windows in a one of the pairwise compared bins is very small (generally <5) while there exists one or a small number of windows with many (generally >5) landslides. These outlier windows can significantly increase the average landslide abundance in a particular bin but have less of an impact on the Wilcoxon rank sum test, which is more sensitive to the overall distribution and abundance of landslides across many windows.
In the analyses presented here, we group Martian terrain into 10 stability bins, corresponding to 10% increments of unstable slope abundance. In general, decreasing the number of bins acts to increase the number of windows in each bin, increasing the odds of differences being found statistically significant. However, larger bins group a greater range of terrains together, leaving more possibility for a difference in landslide abundance between regions to be explainable as a difference in topography. Sensitivity tests performed with different numbers of bins found that our results were not generally sensitive to the number of bins used, for ranges of ∼5-20 bins.

Examination of the Role for Faulting
To investigate the possible impact of tectonic seismicity on Martian landslides, we applied an inventory of surface faulting as mapped by Knapmeyer et al., (2006). They map faults ranging from 4 to 1445 km in length. The previously described method of window-based analysis was repeated, except that for the faulting case we compared populations of windows that were within 50 km of a mapped Martian fault to populations of windows that were more than 50 km from a mapped fault. Extremely large (M8-9) earthquakes on Earth can trigger landslides more than 200 km from the epicenter, but for M6-8 earthquakes 50-100 km ROBACK AND EHLMANN 10.1029/2020JE006675 5 of 14 is a more typical distance to the furthest landslides (Keefer, 1984). The distance of 50 km was also chosen in part because increasing the buffer around faults beyond this distance greatly decreases the amount of available terrain that is classified as far from a fault, hindering the ability of our statistical test to identify significant differences. Additionally, we do not consider the ages of either landslide scars or fault lines in our method; the ages of landslide scars are not reported by Crosta, Frattini (2018b), as many landslide scars are too small to be crater-counted, and the process would be prohibitively labor intensive for an inventory of >3,000 landslides.

Results
We observe that, in general, topography is a primary control on the distribution of Martian landslides. As expected, windows with abundant unstable terrain are much more likely to host a landslide, but even terrains that are >90% unstable have landslides in only ∼20-40% of the terrain area windows ( Figure 5). However, significant differences in landslide abundance not related to topography are observed between some regions of the Martian surface. In particular, regions around Tharsis show more landslides than regions elsewhere on Mars, even when controlling for percentage unstable topography and utilizing different factors of safety ( Figure 6).
These differences are most pronounced in comparing Noctis Labyrinthus (Figures 6a and 6c) and Valles Marineris (Figures 6b and 6d) to Noachian terrains. Landslides are more common in these two regions in the heart of Tharsis than other Martian regions, and comparison of most bins of similar topographies show differences in landslide abundance that are significant at >95% confidence level. These differences persist across different choices of the threshold FOS value used to define "stable" and "unstable" terrain. In other regions around Tharsis, such as Olympus Mons (Figure 6e) and Kasei Valles (Figure 6f) areas, landslides are more common relative to Noachian terrain, but the differences in landslides do not meet our criteria for statistical significance for most bins, in part due to the small sizes of these regions and correspondingly small populations of topographic "windows" with few landslides each. We also attempted to compare Olympus Mons to Elysium Mons, in an effort to compare large volcanic edifices near and far from Tharsis. However, there is a general lack of steep topography in the Elysium Mons area comparable to the aureoles surrounding the Olympus Mons edifice, which are the origin of most of Olympus's landslides. Thus, the only comparable bins between the two regions are low-instability bins with very few landslides; these bins do not host statistically different landslide abundances.
Comparisons between different parts of Tharsis (Figure 7) generally do not reveal statistically significant differences in landslide abundance, although landslides are more common in Valles Marineris and Noctis Labyrinthus relative to other Tharsis areas. Additionally, we compare the eastern and western parts of the Valles Marineris region to see if a contrast in the composition of subsurface materials exhumed by impact craters, which was reported by Quantin et al. (2012), has any apparent impact on landslides. There is no statistically significant variation in landslide frequency across this lithologic contrast (Figure 7c).
Mapped faulting also does not appear to be a major control on mapped landslides (Figure 8). A global comparison of mapped landslides for areas within 50 km of a mapped surface-exposed fault versus areas further away shows variation only in the most unstable terrain (Figure 8a). However, the terrain driving this variation turns out to be concentrated almost exclusively in the Valles Marineris and Noctis Labyrinthus areas (Figures 8b and 8c). When these areas are excluded, no significant differences associated with windows close to or far from faulting are observed elsewhere across the planet. In Valles Marineris, areas with >90% unstable terrain have statistically higher landslide occurrences if < 50 km from a mapped fault (Figure 8c

Global Landslide Prevalence
Our results comparing landslide occurrence in terrains of equivalent instability suggest that the abundance of landslides in regions around Tharsis, and particularly near Valles Marineris, is not solely related to the presence of abundant steep topography in this area. Landslide abundance correlates with topographic stability in Tharsis, but not in Noachian Highlands terrain, for the values of the "unstable" FOS threshold presented here. The cause of this lack of correlation may be due to the fact that most Noachian Highlands landslides are hosted in the walls of impact craters. The slopes of these impact craters are rather short in length, meaning that they occupy rather limited portions of the 20-km-sided windows considered in our analysis; this provides a population of landscape windows that can host landslides despite having comparatively small areal extents of unstable slopes. Larger craters with longer hillslopes filling larger portions of our 20-km windows do not exist in heavily cratered Noachian terrain, because craters above 3-8 km in diameter take on complex-crater morphologies (Pike, 1980) with wall slope angles set by the spontaneous gravitational collapse of the bowl-shaped "transient crater" during an impact event.
Large areas of unstable slopes in Noachian terrain are uncommon; the largest of these is found in the fretted terrain in Nilosyrtis Mensae, which contains few visible landslides. This region shows extensive exposures of hydrated minerals (  second row (c), (d) indicates curves calculated using an unstable FOS threshold = 4. Areas in white are statistically significant while shaded areas mark parts of the curve in which the differences in landslide abundance between the compared regions are not significant at a 95% confidence level. In the lower panels, the solid lines correspond to the left y-axis and the dashed lines to the right y-axis. also may have experienced extensive glacial modification (Levy et al., 2007;Davila et al., 2013), which may have erased traces of landslides which may have occurred in the region's early history.
Interestingly, the Aeolis Mensae region, a geomorphologically similar region with fretted terrain along the boundary of the Noachian highlands, contains far more landslides that Nilosyrtis Mensae, though unstable slopes are actually less abundant (Figure 9). The Elysium Mons volcanic edifice, located ∼500 km from landslide-affected regions in Aeolis Mensae, is farther away than the distances over which landslides are observed to be seismically triggered on Earth, even by extremely large earthquakes (Keefer, 1984). A possible factor is lithology, as portions of Aeolis Mensae are interpreted to be volcanic ash (Kerber et al., 2011), and interbedded ash layers might provide planes of weakness for sliding.
The main question arising from our work is the cause of the elevated landslide abundance in Valles Marineris, and the Tharsis rise region more broadly. Additional factors we have considered include possibly enhanced Tharsis seismicity, the possible presence of a "weak layer" or layers low in the stratigraphy in the Valles Marineris area, and a possible atypical abundance of groundwater in the vicinity of the Valles Marineris area.   (d) and (e). In panel C, the east/ west division (see Figure 4B) follows the contrast in subsurface material properties observed by . In panel D, an unstable FOS threshold of 3 is used to compute slope stability; in all other panels the unstable FOS threshold is assumed to be 4. Shading and axes correspondence to lines are as in Figure 6.

A Possible Role for Seismicity?
Seismicity associated with the extensive tectonic and volcanic activity surrounding Tharsis may be a key contributor to landslide abundance, as suggested previously in the literature (Quantin, Allemand, Mangold, & Delacourt, 2004b). There are far more faults around Tharsis than elsewhere on Mars (Golombek et al., 1992), implying more seismicity near Tharsis (Anderson et al., 2001(Anderson et al., , 2008. Although, as previously discussed, landslides do not follow patterns of mapped faulting in general across Mars, faults may have also been much more active in the vicinity of Tharsis, as would be expected given the larger scale of the chasmae and volcanic edifices present there. Indeed, there is a statistically significant increase in landslides in windows with the most unstable topography (>90%) when those windows are <50 km from a fault (Figure 8) in the Valles Marineris region. Additionally, Tharsis may harbor other deep faults obscured by comparatively recent volcanic activity. If seismicity is indeed the cause of most landslides in Valles Marineris, it may persist and be observable into the present day given the apparently youthful age of many landslide deposits (Quantin, Allemand, & Delacourt 2004a). Seismic data from InSight or future network missions may yet to fully map Martian global seismicity.

A Role for Heterogeneous Subsurface Lithologies?
Another idea that has been suggested in the literature (Montgomery et al., 2009) is a rheological "weak layer" at depth as a possible facilitator of large-scale gravity spreading in the Tharsis area. Such a layer could be dictated by mineralogy, for example clay minerals (e.g., Watkins et al., 2015; or by the physical properties of the rock, for example loosely consolidated sediments or pyroclastics (e.g., Bandfield et al., 2013). We mapped source areas of landslides throughout the Valles Marineris and Noctis Labyrinthus regions, searching for debris-free landslide scarps from which the elevation of the bottom of the landslide failure plane could be estimated (e.g., Figure 10, panel A). We mapped an area including 1,766 landslides mapped by Crosta, Frattini, et al. (2018b), but could only find clear basal scarps for 109 (∼6%) of these landslides. This in large part reflects modes of failure in Valles Marineris; most landslides in the canyon are rock avalanches found below steep cliffs which lack a clear source area (e.g., Figure 10  The 109 identified basal scarps are clustered at low elevations in Valles Marineris ( Figure 11) and are found preferentially in the eastern part of the canyon system. The preference for the eastern part of the canyon for this form of scarp is interesting as Quantin et al. (2012) previously reported a change in subsurface lithology from west to east. Most of the successfully identified basal scarps are relatively small (∼1 km 2 or less) and are located at or below the elevation (about −2000 m) at which low-calcium pyroxene (LCP)-rich, massive light-toned rocks are reported to typically occur Quantin et al., 2012) in eastern Valles Marineris. Relative to the distribution of topography in Valles Marineris, landslides are moderately more abundant than expected between −2000 and −4000 m elevation. It may be that these LCP lithologies or their boundaries with adjacent units are more conducive to the formation of rotational slumps of this scale; these rock units are described as "massive" Quantin et al., 2012), implying a lack of strong jointing or layering observable from orbit. Terrestrial studies have found that landslides in strongly jointed or layered rocks are much more likely to take other forms including rock falls, topples, and rock ROBACK AND EHLMANN 10.1029/2020JE006675 10 of 14 avalanches (Guzzetti et al., 1996;Hermanns & Strecker, 1999). However, slumping, especially translational slumping along a plane of weakness is still possible. Smooth rotational slumps can also still form if the scale of failure is much larger than the scale of layering or jointing, as intact rock is fractured between moving planes of weakness.
Despite this clustering of slumps at low elevations, the overall abundance of landslides remains similar to the lanslides from the lower-elevation slopes of Valles Marineris to Noctis Labyrinthus, where chasma walls host landslides at elevations of 3000-7000 m, above the levels at which much landslides are observed in Valles Marineris, and above the levels of previously suggested weak layers in the stratigraphy. If stratigraphy in the Tharsis area is assumed to be generally horizontal, this variation in elevations of landslides implies that landslides are hosted in a variety of rock units of different ages. Lower elevation outcrops in the Valles Marineris system have been interpreted as outcrops of primitive Noachian crust Quantin et al., 2012), while lava flows associated with Tharsis are thought to comprise the stratigraphy at higher elevations (Murchie et al., 2009), though lavas may also be present in the lower elevations of Tharsis (Golombek & Phillips, 2010). Additionally, landslides are similarly present across the differences in subsurface materials exhumed by impacts, observed by Quantin et al. (2012). This prevalence of landslides is an argument against characteristics of any particular lithologic unit controlling landslide distribution, although the morphological expression of landslides may vary depending on strength and homogeneity of the particular rock unit and lithologic heterogeneities at smaller scale than observable from orbit could still play a role in dictating the planes of failure.
The layer with abundant slump scarps may invoke comparisons to the weak layer suggested to exist by Montgomery et al. (2009) in the lower portions of Valles Marineris's stratigraphy, and to the low thermal inertia, putative volcaniclastic material inferred to exist by Bandfield et al. (2013). It is possible that material in the higher elevations of the Valles Marineris stratigraphy is indeed more well lithified than the lower-elevation massive layers, but the effective material strength at large scales may be set by systems of joints and fractures formed naturally from post-eruption cooling of the lava flows often assumed to comprise the higher parts of Valles Marineris stratigraphy. In summary, though slumps with clear basal scarps appear to be concentrated in a relatively low-elevation unit in the Valles Marineris system, these slumps represent only a minority of all landslides present near Valles Marineris; most are rock avalanches which lack a clear basal scarp. This makes it difficult to tie most landslides to particular stratigraphic units. However, given our coarse knowledge of lithologic units on Mars, lithologic controls could nonetheless be important.

A Role for Groundwater?
Groundwater presence, perhaps episodically related to periodic glacier formation on Tharsis, predicted by climate models (Fastook et al., 2008), and volcanic activity, is another possible contributor to enhanced landslide activity. Groundwater fills pore spaces in rock; this both makes a bulk volume of rock heavier and induces a pore pressure which reduces stress normal to the failure plane, but does not reduce the shear stress parallel to the failure plane, leading to easier failure of slopes. Numerous outcrops of aqueous minerals are found in the Valles Marineris region (e.g., Ehlmann & Edwards, 2014;Murchie et al., 2009), and many of these outcrops are interpreted to be relatively young (Milliken et al., 2008;Weitz et al., 2013), implying that groundwater and surface water may have persisted in this area longer than in other regions of the Martian surface. Given the inferred young age of many Valles Marineris landslides (Quantin et al., 2004b), these relatively young mineral deposits could advance a case for groundwater as a potentially important factor. Precise age dates for aqueous deposits on Valles Marineris, acquired from in situ or sample return missions, could confirm the hypothesis of persistent groundwater in Valles Marineris.

Conclusions
By applying a new global inventory of landslides and slope-stability modeling techniques, we analyzed spatial patterns of Martian landslides while controlling for the stability of topographies across different Martian terrains. The results indicate an abundance of landslides around Tharsis, and particularly around the Valles Marineris area. The effect is not solely due to the abundant steep topography in these areas. We analyzed for but did not find a clear local lithologic or stratigraphic control on landslide occurrence from heterogeneities in the crust. Other possibilities to explain the increased occurrence of landslides in Tharsis include (1) regionally widespread Tharsis weak unit(s), such as from interbedded ashes and lavas; (2) seismic activity related to the Tharsis's geological activity, and (3) possible groundwater near Valles Marineris into the Amazonian. Given the apparently young ages of many landslide deposits in Valles Marineris (Quantin et al., 2004b), continued modern day analysis of lithologies in Valles Marineris and observations of Martian seismicity may act to strengthen or rebut the first two hypotheses.

Data Availability Statement
Our dataset of slope-stability calculations generated for this work, as well as a list of all HRSC DEMs used and all window-by-window calculations of landslide count and unstable slope abundance, for different choices of the "unstable" slope threshold for all studied regions of Mars is publicly archived: Roback, K. P., and Ehlmann, B. L., 2020, Martian Slope Stability Data, version 1.0. CaltechDATA. https://doi.org/10.22002/ D1.1617.