P-A Tesson

S J Conway
email: susan.conway@univ-nantes.fr

N Mangold

J Ciazela

S R Lewis

D Mège

Evidence for thermal-stress-induced rockfalls on Mars Impact Crater Slopes

Keywords: Mars, surface, Thermal stress, Ices, Solar radiation, Weathering, Highlights

Here we study rocks falling from exposed outcrops of bedrock, which have left tracks on the slope over which they have bounced and/or rolled, in fresh impact craters (1-10 km in diameter) on Mars. The presence of these tracks shows that these rocks have fallen relatively recently because aeolian processes are known to infill topographic lows over time. Mapping of rockfall tracks indicate trends in frequency with orientation, which in turn depend on the latitudinal position of the crater. Craters in the equatorial belt (between 15°N and 15°S) exhibit higher frequencies of rockfall on their N-S oriented slopes compared to their E-W ones. Craters >15° N/S have notably higher frequencies on their equator-facing slopes as opposed to the other orientations. We computed solar radiation on the surface of crater slopes to compare insolation patterns and rockfall spatial distribution, and find statistically significant correlations between maximum diurnal insolation and rockfall frequency. Our results indicate that solarinduced thermal stress plays a more important role under relatively recent climate conditions in rock breakdown and preconditioning slopes for rockfalls than phase transitions of H 2 O or CO 2, at mid and equatorial latitudes. Thermal stress should thus be considered as an important factor in promoting masswasting process on impact crater walls and other steep slopes on Mars.

1-Introduction

Geomorphological processes are active at the surface of Mars at present day. Repeat-coverage and high resolution images (better than 50 cm/pix - [START_REF] Mcewen | Mars Reconnaissance Orbiter's High Resolution Imaging Science Experiment (HiRISE)[END_REF] have revealed terrestrial-like gravitational mass movements [START_REF] Tsige | Modeling of Landslides in Valles Marineris, Mars, and Implications for Initiation Mechanism[END_REF] ranging in size from landslides [START_REF] Lucchittta | A large landslide on Mars[END_REF] to single rockfalls [START_REF] Roberts | Possible evidence of paleomarsquakes from fallen boulder populations, Cerberus Fossae, Mars[END_REF]Senthil Kumar et al., 2019). The spatial distribution of tracks left on slope-materials by falling rocks has been used as a passive seismometer on Mars and the Moon [START_REF] Roberts | Possible evidence of paleomarsquakes from fallen boulder populations, Cerberus Fossae, Mars[END_REF]Senthil Kumar et al., 2016;Senthil Kumar et al., 2019). On Mars, few studies have been carried out investigating the factors that control rockfall activity. From single images alone, rockfalls appear to occur in a similar way as they do on Earth: a clast detaches from a cliff and rolls or bounces downslope, while leaving a track in the surface (Fig. 1 &2A). However, weathering mechanisms leading to rockfall could differ substantially given the environment and climatic differences between Mars and Earth. For instance, liquid water plays a major-role in terrestrial slope stability (e.g., [START_REF] Terzaghi | Stability of Steep Slopes on Hard Unweathered Rock[END_REF]. The involvement of water in active surface processes on Mars is widely debated (e.g., [START_REF] Schmidt | Formation of recurring slope lineae on Mars by rarefied gas-triggered granular flows[END_REF][START_REF] Ojha | Seasonal Slumps in Juventae Chasma, Mars[END_REF]. Since liquid water is generally not stable today on the martian surface [START_REF] Farmer | Liquid water on Mars[END_REF][START_REF] Haberle | On the possibility of liquid water on present-day Mars[END_REF][START_REF] Hecht | Metastability of liquid water on Mars[END_REF], it is unlikely, however, that water plays such an important role in slope instability on Mars as it does on Earth. [START_REF] Tsige | Modeling of Landslides in Valles Marineris, Mars, and Implications for Initiation Mechanism[END_REF]. On the other hand, [START_REF] Eppes | Cracks in Martian boulders exhibit preferred orientations that point to solar-induced thermal stress[END_REF] have linked boulder fracture patterns from Mars Exploration Rover (MER) Spirit images to directional solar-induced thermal stress on Mars. Further, on Earth, [START_REF] Collins | Rockfall triggering by cyclic thermal stressing of exfoliation fractures[END_REF] showed that, on Earth, rockfalls can be linked to thermal stress. Therefore, we consider thermal stress as a potential weathering mechanism linked to rockfall activity on Mars at present day. In order to evaluate its relative importance, we have catalogued recent rockfalls on the slopes of fresh impact crater walls across a range of latitudes to highlight different patterns of frequency, block size and orientation. Craters walls are receiving different amount of insolation depending on their orientation and latitude. Therefore, they are exposed to relative differences of potential thermal stress intensity. If a link exists between thermal stress and rockfall activity, then there should be an aspect relationship between them.

2-Theoretical Background

2.1-Slope Stability

Slope stability is a well-documented topic on Earth, especially in terms of geohazards related to landslides and rockfalls. Mechanically, a rockfall can be modelled as small-scale mass movement. It is a result of the breakup of an individual rocky object from the top of a slope [START_REF] Selby | Hillslope materials and processes[END_REF]. Figure 1 outlines the factors that can contribute to a rockfall event. Exposed rock (i.e. an outcrop) at the top of a topographic slope will be subject to gravity, and cohesion counteracts gravity's pull resisting material failure and rockfall. A steep slope is the primary factor in controlling rockfall as the tangential component of the gravitational acceleration is a function of the slope angle. For failure to occur, the material needs to undergo a loss of cohesion as well. In a natural setting, environmental stresses lead to the growth of fractures via subcritical fracture that occurs due to stress magnitudes much lower than the critical strength of the rock (e.g. [START_REF] Eppes | Mechanical weathering and rock erosion by climate-dependent subcritical cracking[END_REF]. Over time, subcritical fracture growth results in loss of cohesion and reduction in material strength overall.

2.2-Longevity of rockfall tracks

Aeolian processes are active at the surface of Mars today [START_REF] Bourke | Recent aeolian dune change on Mars[END_REF][START_REF] Hansen | Seasonal erosion and restoration of Mars' northern polar dunes[END_REF] and study of the wheel tracks left by MER Opportunity and Spirit show that they persist for more than one martian year [START_REF] Geissler | Gone with the wind: Eolian erasure of the Mars Rover tracks[END_REF][START_REF] Sullivan | Wind-driven particle mobility on Mars: Insights from Mars Exploration Rover observations at "El Dorado" and surroundings at Gusev Crater[END_REF]. Rover tracks are susceptible to infill because they are located in areas where aeolian grains are easy to transport [START_REF] Sullivan | Wind-driven particle mobility on Mars: Insights from Mars Exploration Rover observations at "El Dorado" and surroundings at Gusev Crater[END_REF][START_REF] Sullivan | Aeolian processes at the Mars Exploration Rover Meridiani Planum landing site[END_REF]. On the other hand, rockfall tracks are located on sloping terrains, and may have slower infill rates because they are not exposed to the same aeolian transport mechanisms. To be conservative therefore we estimate that a meter-scale rockfall track could persist for more than a thousand years, and we acknowledge that this rate will vary from one context to another (aeolian activity, track depth and width, and substrate). Although this time period is short over geological timescales, orbital forcing related climate cycles are relevant [START_REF] Laskar | Long term evolution and chaotic diffusion of the insolation quantities of Mars[END_REF] and should be considered. We assume that if a track can still be observed, the rockfall should be recent and the weathering mechanism involved should still be ongoing at the surface today.

Thermal stress

High temperature contrasts experienced by rock surfaces in diurnal cycles lead to thermal expansion of the material at the surface, but less so at depth, and hence rock breakdown. This mechanism is widely studied on Earth. Rock breakdown linked to solar-induced thermal stresses is known to occur in very arid environments in both cold (e.g. Antarctica, [START_REF] Lamp | Thermal stress weathering and the spalling of Antarctic rocks[END_REF] and hot deserts [START_REF] Hall | The role of thermal stress fatigue in the breakdown of rock in cold regions[END_REF][START_REF] Mckay | High-frequency rock temperature data from hyperarid desert environments in the Atacama and the Antarctic Dry Valleys and implications for rock weathering[END_REF][START_REF] Eppes | Cracks in desert pavement rocks: Further insights into mechanical weathering by directional insolation[END_REF]. Recent studies demonstrated that high stress can be correlated to high absolute maximum temperature and diurnal temperature range. These temperature parameters could thus be used as proxies of thermal stress [START_REF] Boelhouwers | Critical assessment of the 2°cmin-1 threshold for thermal stress weathering[END_REF][START_REF] Collins | Thermal influences on spontaneous rock dome exfoliation[END_REF][START_REF] Eppes | Deciphering the role of solar-induced thermal stresses in rock weathering[END_REF][START_REF] Molaro | Rates of temperature change of airless landscapes and implications for thermal stress weathering[END_REF]. Thermal stress could cause rock breakdown on airless bodies (e.g. The Moon) [START_REF] Molaro | Grain-scale thermoelastic stresses and spatiotemporal temperature gradients on airless bodies, implications for rock breakdown[END_REF][START_REF] Molaro | Thermally induced stresses in boulders on airless body surfaces, and implications for rock breakdown[END_REF] as well as on Earth [START_REF] Eppes | Cracks in desert pavement rocks: Further insights into mechanical weathering by directional insolation[END_REF][START_REF] Warren | Automated field detection of rock fracturing, microclimate, and diurnal rock temperature and strain fields[END_REF]. Therefore on Mars, arid conditions and high surface temperature gradients (e.g. [START_REF] Spanovich | Surface and near-surface atmospheric temperatures for the Mars Exploration Rover landing sites[END_REF] suggest that thermoelastic stress might also occur on the surface at present day. [START_REF] Viles | Simulating weathering of basalt on Mars and Earth by thermal cycling[END_REF] have been able to lower the strength of pre-stressed basaltic rocks under martian atmospheric conditions by exposing them to Mars-like diurnal cycles temperature variations. Based on boulder fracture patterns from Mars Exploration Rover (MER) Spirit images, [START_REF] Eppes | Cracks in Martian boulders exhibit preferred orientations that point to solar-induced thermal stress[END_REF] advocate active thermal stress-related rock breakdown. This directional relationship between thermal stress on Mars and the fracture patterns could exist also on rock walls. Growth of fractures caused by repetitive thermal stress cycles could lead to a fall, as known from Earth [START_REF] Collins | Rockfall triggering by cyclic thermal stressing of exfoliation fractures[END_REF][START_REF] Vargas | On the effect of thermally induced stresses in failures of some rock slopes in Rio de Janeiro, Brazil[END_REF][START_REF] Gischig | Natural hazards: Cracking cliffs feel the heat[END_REF]. If such a relation exists on Mars, there should be a correlation between number of rockfalls and maximum temperatures or/and maximum temperature ranges on slopes.

3-Material and methods

3.1-Site Selection

Impact craters are widely distributed over the martian surface, and for that reason, can be used as sample locations to test the potential factors controlling rockfall distribution. Conveniently, they are circular allowing a relatively unbiased assessment of slope-orientation influence. Here, we focus on relatively fresh impact craters (Fig. 2B), from 1 to 10 km in diameter (Fig. 3) to reduce the influence of slope-inheritance from other long-term processes. Impact craters with a diameter <10 km tend to have a "simple" bowl-shaped morphology, whereas craters with diameters >10 km tend to have a more complex morphology [START_REF] Melosh | Impact cratering: A geologic process[END_REF], including multiple wall terraces, which could complicate our analysis. Site selection was made by cross-referencing the global database of Mars impact craters from [START_REF] Robbins | A new global database of Mars impact craters ≥1 km: 2. Global crater properties and regional variations of the simple-to-complex transition diameter[END_REF], and image data availability. Selected craters are located between 50°S and 40°N with most of them being located between 30°S and 30°N (Fig. 3) to avoid the latitude-dependent mantle (LDM). LDM is a meter-thick layer of ice and dust covering at least 23% of the surface that extends from the mid-latitudes to the poles [START_REF] Conway | Decameter thick remnant glacial ice deposits on Mars[END_REF][START_REF] Kreslavsky | Mars: Nature and evolution of young latitude-dependent waterice-rich mantle[END_REF]. It is a draping unit which likely formed during the many orbital variation-driven climate excursions that have occurred during the Amazonian period [START_REF] Kreslavsky | Mars: Nature and evolution of young latitude-dependent waterice-rich mantle[END_REF]. Fresh craters covered by LDM would introduce a bias in the results since they would tend to cover the slopes making rockfall tracks harder to observe, we decided to discard craters where it can be observed.

3.2-Dataset and mapping

Throughout this paper, we use the clast nomenclature suggested by [START_REF] Bruno | Something more than boulders: A geological comment on the nomenclature of megaclasts on extraterrestrial bodies[END_REF]. Because most of the clasts mapped in this study range from 1 to 10 m, we will refer to them as "blocks" as opposed to "boulders", the term usually used, which refers to clasts smaller than 1 m in diameter according to this nomenclature. In order to map rockfalls in impact craters, we have used images from High Resolution Imaging Science Experiment (MRO) instrument aboard the Mars Reconnaissance Orbiter (MRO), which provides a spatial resolution up to 0.25 m/pixel [START_REF] Mcewen | Mars Reconnaissance Orbiter's High Resolution Imaging Science Experiment (HiRISE)[END_REF] as listed in supplementary material, Table 1. Map-projected images were integrated in ArcGIS© 10.4 to identify and map recent rockfall tracks. We digitized the tracks left by the clasts as they fell (Fig. 2A) and where possible the long axis of block at the end of the track. From these polylines, we calculated their angle with respect to geographic north. For length measurements, we used a sinusoidal projection centered on the crater to avoid the distortion linked to map projection. We used two different ways of evaluating frequency, using a normalized and a non-normalized representation for the distribution of rockfall track orientations. For both methods, we first calculated the number of rockfall tracks in 20° azimuth bins for each crater in a given latitude range. For the normalized distribution, we calculated for each crater the percentage of the total number of tracks in that crater for each azimuth bin. For each crater, the percentage of rockfalls in each bin is then relative to the number of rockfalls in the crater. We then calculated the mean percentage of tracks for each bin in all craters in this specific latitude range (Fig. 4,right). For the non-normalized distribution, the number of tracks in each bin was summed for every azimuth bin in each crater in each latitude set (Fig. 4 left), then the percentage is calculated relative to the total number of rockfall in the latitude range.

3.3-Topographic measurements

In order to assess variations in slope angle from crater-to-crater, or for different slope-orientations, for a sub-sample of our craters we examined the slope angle at the bedrock outcrops. Where it was possible, we generated Digital Terrain Models (DTM) at 24 m/pix from stereo-pairs of MRO Context camera (CTX) images at ~6 m/pix [START_REF] Malin | Context Camera Investigation on board the Mars Reconnaissance Orbiter[END_REF] using the Ames Stereo Pipeline [START_REF] Broxton | The Ames Stereo Pipeline: Automated 3D surface reconstruction from orbital imagery[END_REF]. Before gridding, the generated point clouds were vertically controlled to Mars Orbiter Laser Altimeter (MOLA) elevation point data (PEDR) from Mars Global Surveyor. DTM accuracy was estimated by calculating the root mean square (RMS) between the elevation of the MOLA points and the mean elevation of the CTX DTM in a circle of 168 m in diameter surrounding those points. This diameter corresponds to the pulse diameter estimated from the point spread function (PSF) of MOLA instrument [START_REF] Neumann | Mars Orbiter Laser Altimeter pulse width measurements and footprint-scale roughness[END_REF]. RMS results are reported in supplementary material, Table 2. In order to estimate the slope of the crater wall outcrops, we split each crater into orientation bins, totaling 18 arc-segments each covering 20° of azimuth. In each bin, we extracted the elevation of every point between 50 and 200 m from the crater rim from where the clasts should typically fall. This elevation range is corresponds to where rock outcrops are generally observed [START_REF] Conway | Glacial and gully erosion on Mars: A terrestrial perspective[END_REF], and may vary from a crater to another. We extract the slope value for each bin by taking the slope of a linear fit between the elevation and the distance from the rim for every DTM pixel.

4-Results

We recorded 2040 recent rockfall tracks in 39 impact craters among which, 1584 tracks had an associated clast. Figure 4 displays the frequency of rockfalls in craters with respect to the orientation of the crater wall, for different latitude ranges. Both the normalized and non-normalized plots show similar overall trends in each latitude range. The mid-latitude craters have the most rockfall tracks on the equator-facing slopes, in the northern and southern hemisphere. Northern mid-latitude craters have the highest number of rockfall tracks on the N-NE slopes (> 40% combined between N0° and N80°) in both normalized and non-normalized plots. S slopes have very few recorded rockfall tracks (< 20% in total from N100° to N260°) in both normalized and non-normalized plots. The mean vector inferred from these distributions is respectively N13.3°±19.1 and N15.1°±19.6 (95% confidence interval) for normalized and non-normalized distribution. Rayleigh and Rao's test for the null hypothesis of uniformity yield P-value <0.01 for normalized distribution and 0.05 for non-normalized distribution. In the southern mid-latitudes, equator-facing slopes of craters have the most rockfall tracks. The normalized results show that up to 70% of the total number of rockfall tracks occur on the S slopes (i.e. north-facing, from N100° to N260°), whereas the non-normalized results only ~60% of total rockfall tracks occur in the same range. The non-normalized plot has larger percentage of rockfall tracks on the pole-facing slope (20% in total between N340° and N40°) than the normalized results (2-3% on average in the same bins). The mean vector inferred from these distributions is respectively N178.3°±27.6 and N141.6°±15.0 (95% confidence interval) for normalized and non-normalized distribution. Rayleigh and Rao's test for the null hypothesis of uniformity yield P-value <0.01 in both cases. At equatorial latitudes, there are more rockfall tracks on N-S-oriented slopes compared to E-W-oriented slopes. Specifically, the non-normalized distribution shows two peaks of 8% between N300°-N320° and N0°-N20°, and another one from N120° to N160° (15% combined). In the normalized distribution, the fraction of rockfall tracks is greater between N340° and N20° (18% combined) and from N180° to N200° (9%) compared to the non-normalized distribution. A peak in the N140°-N160° bin (7%) is also apparent in the normalized distribution. Mean vector from N270° to N90° (northern slopes) is respectively N0.0°±12.9 and N322.5°±4.0 for normalized and non-normalized distribution, while it is N176.6°±13.6 and N165.2°±4.0 for southern slopes, with P-value < 0.01 in each case.

5-Rockfall clast size

5.1-Magnitude-Frequency of rock volume

Figure 5 shows the cumulative frequency-volume distribution of martian rockfalls mapped in this study. Magnitude-cumulative frequency (MCF) distribution is commonly modelled by a power law in the middle section, with a deviation at the low and high sections [START_REF] Corominas | Magnitude and frequency relations: are there geological constraints to the rockfall size? Landslides 15[END_REF]. The power-law exponent (referred to as "scaling parameter") is thought to represent the fracture susceptibility of the rock mass under consideration. The scaling parameter of rockfalls in impact craters from our catalog is -1.23. On Earth, the scaling parameter of the power law usually varies from -0.90 to -0.40, depending on geological, morphological and climatic conditions [START_REF] Corominas | Magnitude and frequency relations: are there geological constraints to the rockfall size? Landslides 15[END_REF], but also varies based on the technique used to measure the rock-dimensions. Using Terrestrial Laser Scanner (TLS)-generated 1 m/pix DTM to measure in-situ detachable rock volumes on a rock cliff (chute of Forat Negre, Andorra), [START_REF] Mavrouli | Size Distribution for Potentially Unstable Rock Masses and In Situ Rock Blocks Using LIDAR-Generated Digital Elevation Models[END_REF] found an exponent of -1.3 while previous study yielded a value of -0.9 [START_REF] Santana | Magnitude-frequency relation for rockfall scars using a Terrestrial Laser Scanner[END_REF] for the same investigated area when using rockfall scar measure with TLS. The higher value of scaling factor for Mars compared to Earth could be explained by a difference in mapping method but could also represent a difference in the rock mechanics between the two planets.

5.2-Median length of clasts

Figure 6 shows the median long-axis length of rockfall clasts in each impact crater with respect to latitude. The median size of recent rockfalls in the southern hemisphere is higher at the mid-latitudes than near the equator in Hesperia Planum and other locations on Mars (Fig. 6). The median size of blocks is greater than 2.5 m at latitudes >35°S and decreases down to 2 m close to the equator. A weaker similar trend is observed in the northern hemisphere in Syrtis Major Planum and in the martian northern lowlands. To confirm whether this trend exists in the northern hemisphere more data would be needed >25°N.

6-Influence of slope on rockfalls

The sub-sample of ten craters with DTMs allows us to investigate the relationship between rockfall frequency and slope angle at the source. The aim of this analysis is to determine whether a systematic variation in slope angle with orientation and latitude could explain the frequency distribution of rockfall tracks reported in Figure 4. These craters are equally distributed within the latitude range of our total sample set of 39 craters, and we find that >88% of rockfalls originate from slopes steeper than 32° (Fig. 7A). Mars gravity is only 38% of Earth, therefore the question arises whether internal friction angle differs on Mars and Earth, where it is ~30°. Early experiments carried out by Viking landers have shown that angle of internal friction on Mars appears to be similar to typical Earth values, ranging from 27 to 39° [START_REF] Moore | Physical Properties of the Surface Materials at the Viking Landing Sites on Mars[END_REF]. Therefore, the 32° value we found (Fig. 7A) is in line with expectations. The DTM resolution (24 m/pix) means we are measuring the overall slope value within the assumed rockfall source area, rather than the meter-scale slope from where the rocks detach. The distribution of slope angle values follows a normal distribution for slopes with and without rockfalls (Fig 7B). However, the distribution of slopes with rockfalls is shifted towards higher values of topographic gradient, in accordance with the expectation that increasing slope angle should increase rockfall activity (Fig. 1). Mid-latitude equator-facing slopes on Mars are known to be steeper than pole-facing slopes [START_REF] Kreslavsky | North-south topographic slope asymmetry on Mars: Evidence for insolation-related erosion at high obliquity[END_REF]2018) and this could potentially explain the higher frequency of rockfalls on equatorfacing slopes at our mid-latitude sites. The proposed explanation for this asymmetry is deposition and removal of the LDM [START_REF] Kreslavsky | Mars Climate History: Insights From Impact Crater Wall Slope Statistics[END_REF], although recent work has shown it could be related to enhanced glacial erosion on pole-facing slopes [START_REF] Conway | Glacial and gully erosion on Mars: A terrestrial perspective[END_REF]. Low-angle impact craters can display specific ejecta and interior morphologies [START_REF] Herrick | The planforms of low-angle impact craters in the northern hemisphere of Mars[END_REF], with a latitude-dependent frequency [START_REF] Barlow | Martian impact craters: Correlations of ejecta and interior morphologies with diameter, latitude, and terrain[END_REF]. Non-circular rim morphologies could create azimuthal variations in slope value in walls and influence rockfall activity.

Although rockfalls preferentially occur on steeper slopes, slope steepness alone does not condition rockfalls. Figure 7C shows that a significant proportion of steep slopes do not have any rockfall tracks. For instance, 30% of 36-38° slopes have no rockfall tracks (Fig. 7C). Slopes angles of 42-44° slopes have no recorded rockfall tracks, but they only represent less than 2% of the population, meaning that this low proportion of rockfalls at steep slopes could be a statistical bias. Slopes steeper than 32° (excluding 42-44°) represent 27% of the measured population and have a proportion with rockfall tracks ranging from 52% to 10% (Fig. 7C). Therefore, if rockfall activity was dictated by topographic gradient only, the proportion of steep slopes with rockfalls should become closer to 100% the steeper the slope becomes. Such a correlation is not found in our results. Moreover, our impact craters are relatively fresh with a circular bowl-shaped morphology (Fig. 2B) and have been selected so to avoid influence of slopeinheritance from other long-term processes and LDM (see Section 3.1). Therefore, an anisotropic process must be involved in rock breakdown and rockfall activity on Mars at present-day to explain the observed frequency in occurrence of rockfall tracks in orientation with latitude (Fig. 4).

7-Other sources of rockfalls in impact craters

Ejecta blankets contain abundant clasts (ranging in size from silt to several ten-of-meter-large blocks) and can extent to the upper interior walls of impact craters (Krishna and Senthil Kumar, 2016;[START_REF] Kumar | Impact fragmentation of Lonar Crater, India: Implications for impact cratering processes in basalt[END_REF]. In addition, ejecta layer is underlain by bedrock highly fractured during the impact (Senthil [START_REF] Kumar | Structural effects of meteorite impact on basalt: Evidence from Lonar crater, India[END_REF]Senthil Kumar and Kring, 2008). Both of those impact-related clasts are potential source of rockfalls in craters as they can simply fall due to gravity. However, the distribution of these rockfalls would be random, independent on slope orientation, which is in contrast with our highly oriented distribution. (Fig. 4). As impact craters in this study are relatively young, they also could be exposed to ongoing modification process such as impact crater collapse [START_REF] Melosh | Impact Crater Collapse[END_REF] which could be the source of rockfalls. This factor, however, should also act randomly and is not expected to cause any latitude-dependent orientation pattern of rockfall on crater walls. At mid-latitudes (>30°), crater walls may display gullies that could be a source of rockfalls [START_REF] Harrison | Global documentation of gullies with the Mars Reconnaissance Orbiter Context Camera and implications for their formation[END_REF]. However, gullies between 30° and 40° prevail on the pole-facing slopes [START_REF] Conway | New slope-normalized global gully density and orientation maps for Mars[END_REF] while rockfalls seem to occur preferentially on equator-facing slopes at those latitudes (Fig. 4). One could also argue that gullies have the opposite effect and tend to mitigate rockfall activity, or simply introduce a bias in the mapping itself by reducing their visibility. However, this would need to cause lower detection of rockfalls on craters at >30° latitudes, which is not the case. Our data indicate relatively constant rockfall frequency from 15° to 40/50° in both hemispheres (Fig. 4).

8-Marsquakes

Although most of the rockfall energy is provided by gravity, loss of material cohesion may not be sufficient to trigger rockfall. Other local sources of energy could play this role. Rockfall activity has been used to infer present-day seismic activity on Mars and on the Moon [START_REF] Roberts | Possible evidence of paleomarsquakes from fallen boulder populations, Cerberus Fossae, Mars[END_REF]Senthil Kumar et al., 2016;Brown and Roberts, 2019). During earthquakes, energy decreases with distance from epicenter. [START_REF] Keefer | Landslides caused by earthquakes[END_REF] noted that spatial frequency and intensity of slope instabilities increase closer to earthquake epicenter. For instance, [START_REF] Roberts | Possible evidence of paleomarsquakes from fallen boulder populations, Cerberus Fossae, Mars[END_REF] reported an increase in clast size and spatial distribution of recent block falls on along Cerberus Fossae floor around a specific location and exclude other triggering factors than recent marsquakes. Marsquakes could thus also affect the rockfall distribution. In addition, rockfalls be caused by neighbor impacts or wind but only in rare cases. All in all, however, this phenomenon may occur only locally and could not explain the specific pattern observed here (Fig. 4).

9-Weathering mechanism 9.1-Phase changes of ice

On Earth, phase change of water from liquid to solid is usually responsible for rock-breakdown and can result in a rockfall. Freezing and thawing are rare on Mars as liquid water is unstable under current atmospheric conditions (average pressure of 6 mbar [START_REF] Farmer | Liquid water on Mars[END_REF][START_REF] Haberle | On the possibility of liquid water on present-day Mars[END_REF][START_REF] Hecht | Metastability of liquid water on Mars[END_REF] although special regions have been identified where this could happen (e.g. [START_REF] Chevrier | Stability of perchlorate hydrates and their liquid solutions at the Phoenix landing site, mars[END_REF] and availability of metastable liquid water can affect weathering rate, even during the late Amazonian period [START_REF] De Haas | Local late Amazonian boulder breakdown and denudation rate on Mars[END_REF].

Compared to liquid water, water ice is abundant on Mars and the ice itself can change volume and can produce stresses of >5 MPa [START_REF] Mellon | Small-scale polygonal features on Mars: Seasonal thermal contraction cracks in permafrost[END_REF]. Hence, rock may be broken down by the seasonal and/or diurnal thermal contraction. Ground ice is thought to exist from the mid to high latitudes (>45°) on modern Mars, and has been documented in multiple locations using both orbital [START_REF] Boynton | Distribution of hydrogen in the near surface of Mars: Evidence for subsurface ice deposits[END_REF][START_REF] Byrne | Distribution of mid-latitude ground ice on mars from new impact craters[END_REF][START_REF] Dundas | Exposed subsurface ice sheets in the Martian mid-latitudes[END_REF][START_REF] Mouginot | The 3-5MHz global reflectivity map of Mars by MARSIS/Mars Express: Implications for the current inventory of subsurface H2O[END_REF] and in-situ data [START_REF] Mellon | Ground ice at the Phoenix Landing Site: Stability state and origin[END_REF]. Also, ground H 2 O ice is inferred to exist from the observed distribution of CO 2 seasonal ices on polefacing crater slopes at latitudes as low as 25° in the southern hemisphere (Vincendon et al., 2010a). In addition to water, CO 2 seasonal frost resulting from condensation of atmospheric CO 2 in winter is known to form a continuous >10 g/cm 2 layer which extends from polar caps to 60° in latitude in both hemispheres [START_REF] James | The carbon dioxide cycle[END_REF][START_REF] Kelly | Seasonal polar carbon dioxide frost on Mars: CO2 mass and columnar thickness distribution[END_REF]. Small patches of seasonal frost can also be found in shadowed pole-facing slopes at latitudes as low as 33°S with modeled concentration reaching 8 g/cm 2 [START_REF] Schorghofer | Seasonal surface frost at low latitudes on Mars[END_REF]. Thin layers (<1 mm) of diurnal CO 2 frost also exists at low latitude, down to the equator, on low thermal inertia, dusty units [START_REF] Piqueux | Discovery of a widespread low-latitude diurnal CO 2 frost cycle on Mars[END_REF]. Small concentrations of water ice are associated with CO 2 frost, even at mid-latitudes [START_REF] Schorghofer | Seasonal surface frost at low latitudes on Mars[END_REF][START_REF] Carrozzo | Mapping of water frost and ice at low latitudes on Mars[END_REF]Vincendon et al., 2010b) and surface H 2 O ice frost deposits were observed by the Viking lander 2 at 47.64°N [START_REF] Farmer | Liquid water on Mars[END_REF]. Vincendon et al., (2010b) observed water ice deposits ranging from 2 to 200 μm at latitudes as low as 13°S and 32°N combining data from imaging spectrometers and a modeling approach. These thin deposits are derived from atmospheric humidity (recorded in TES data - [START_REF] Smith | The annual cycle of water vapor on Mars as observed by the Thermal Emission Spectrometer[END_REF] generated by sublimation of ground water ice and polar caps. Higher relative humidity itself can also increase crack growth velocity [START_REF] Nara | Influence of surrounding environment on subcritical crack growth in marble[END_REF]. Relative humidity is greater for a higher water vapor contain and a lower temperature [START_REF] Harri | Mars Science Laboratory relative humidity observations: Initial results[END_REF], meaning crack growth velocity increase could be the highest at the mid-latitude and above, where water content is greater [START_REF] Smith | The annual cycle of water vapor on Mars as observed by the Thermal Emission Spectrometer[END_REF] and temperatures lower. Martian rockfalls linked to phase changes of H 2 O and CO 2 should occur where they are expected to condense and/or be preserved from previous ice ages, namely on pole-facing slopes in the mid to highlatitudes and little to nowhere at the equator. Our results show that rockfalls occur on equator-facing slopes both in the mid and equatorial latitudes (Fig. 4), suggesting that such a trend does not exist at these latitudes, and that another mechanism is involved.

9.2-Solar radiation model

Latitude-dependence of rockfall orientation trends indicates that insolation plays a role on source-rock breakdown and preconditioning rockfall occurrence on impact craters slopes. For this reason, thermoelastic stress may be likely to play this role. To assess the variation of insolation on latitude, we computed solar insolation over a typical DTM of a crater placed at different latitude positions. We used the publicly available HiRISE DTM (DTEEC_002118_1510_003608_1510_A01) of Zumba crater (Fig. 2B) to represent a typical fresh crater. The DTM was reduced to 10 m per pixel and we only considered the crater walls to achieve a reasonable tradeoff between resolution and computation time. For each pixel, a "viewshed", which provides information on the sky visibility in every direction, is first derived from the DTM [START_REF] Rich | Using Viewshed Models to Calculate Intercepted Solar Radiation: Applications in Ecology[END_REF]. The model then computes solar incidence angle of each pixel with respect to slope angle and orientation, and sun position in the sky at a given moment. We ran the model for a martian day (sol) every 10° of solar longitude (Ls) for the whole martian year. For each sol, direct insolation in W/m² is computed every half an hour of local time. The daily mean is then calculated by averaging the solar flux at each step for the entire sol. Equations making up the mathematical framework of our model are all extracted from [START_REF] Appelbaum | Solar radiation on Mars[END_REF]. Details on the model are available in the supplementary materials (Appendix A1). We obtained a raster that represents the maximum diurnal average insolation received by the crater walls at a given latitude over a martian year.

9.3-Orbital parameters

The lifetime of the tracks left by rockfalls could be up to several tens of thousands of years (see section 2.5). This timescale requires that changes in orbital parameters be considered when studying latitudedependent processes depending on insolation. During the last 100,000 years, the eccentricity of Mars' orbit has ranged from 0.075 to 0.118, being 0.093 today [START_REF] Laskar | Long term evolution and chaotic diffusion of the insolation quantities of Mars[END_REF]. Eccentricity influences seasonal contrasts (see section 7.2), and therefore it is unlikely to play a role in the relative rates of weathering for different slope orientations at different latitudes. In the same period the obliquity has ranged from 22.5° to 26.8° (25.2° today - [START_REF] Laskar | Long term evolution and chaotic diffusion of the insolation quantities of Mars[END_REF], which is likely too small to have any noticeable effect on the insolation patterns per orientation with latitude. However, during the same period, Mars has experienced two full precession cycles, implying that the solar longitude at which perihelion occurs has changed significantly. For instance, 22,700 years ago, Mars' perihelion occurred at Ls = 90° [START_REF] Laskar | Long term evolution and chaotic diffusion of the insolation quantities of Mars[END_REF], meaning that maximum insolation was received during the northern summer (Fig. 8) and meaning solar insolation was at its yearly maximum on north-facing slopes at the equator. This is the opposite of the situation today. With perihelion at Ls = 90°, north-facing slopes at the equator would experience maximum diurnal and seasonal temperature contrasts, as it is currently the case for S-facing slopes. Precession cycles should therefore be considered when computing insolation received by crater slopes during the last 100,000 years. For each crater, we have then computed maximum diurnal average insolation for specific values of solar longitude of perihelion to estimate the average over a full precession cycle.

Comparison with rockfall distribution

The insolation model was run at the latitude of each crater, every 45° of solar longitude of perihelion, to cover a full precession cycle. To allow comparison with rockfall distribution, we have binned insolation data with same method used for rockfall tracks (Section 3 -Fig. 4). The bins are then averaged for a specific latitude range and compared to a corresponding rockfall distribution for the same range. Results are plotted in Fig. 9. Altogether, we can observe a significant linear correlation (Pearson r=0.69; R²=0.48; P-value <0.01) between maximum diurnal insolation averaged for a full precession cycle and rockfall frequency (Fig. 9).

Insolation and thermal stress

Thermal stress is intrinsically higher where contrasts of temperature are greater, thus, one can expect it to be stronger where solar flux is higher. [START_REF] Molaro | Grain-scale thermoelastic stresses and spatiotemporal temperature gradients on airless bodies, implications for rock breakdown[END_REF] claim that temporal gradient of temperature alone is a poor proxy for thermal stress at grain-scale and suggest using absolute temperatures and offset from diurnal means. The same conclusions are shared by [START_REF] Boelhouwers | Critical assessment of the 2°cmin-1 threshold for thermal stress weathering[END_REF], as well as [START_REF] Molaro | Rates of temperature change of airless landscapes and implications for thermal stress weathering[END_REF]. Solar flux and temperature are different physical quantities, but accurate temperature models are far more complex than insolation models and would make computation time for such a long timescale far greater. Also, temperature, including peak temperature, is controlled by solar flux and thermal proprieties of the material. Different rock compositions will have different thermal proprieties which also affect thermal stress intensity. However, one can expect thermal proprieties of rocks to be roughly the same within one crater, therefore only solar flux variations would control temperature disparities between walls. Higher the solar flux received by a surface within a sol, higher will be its maximum temperature, and potentially the thermal stresses experienced by the material underneath it. Our results indicate a spatial correlation between maximum diurnal insolation (and potentially average daily temperature range) and rockfall activity at the mid and equatorial latitudes (Fig. 9). Overall, insolation is higher close to the equator than at the mid-latitudes, suggesting that thermal-related weathering should be more efficient here. This could provide an explanation for why we observe a lower median size of clasts (yet more numerous) at lower latitudes (Fig. 6). These rockfalls could have occurred either during previous precession conditions and the track preserved or be recent rockfalls derived from bedrock weakened during previous precession conditions. The latter is more likely, as thermal stress is only responsible for preconditioning outcrops for rockfalls, not triggering them (Fig. 1). Equatorial north-facing slopes could have been weathered over longer-period but rockfalls would occur more recently. Hence, our results also indicate that the delay between weathering and rockfall is probably constrained to recent orbital conditions (last 100,000 years) as it does not seem to be affected by obliquity variations which have a longer timespan (e.g. >10° at >100 ka - [START_REF] Laskar | Long term evolution and chaotic diffusion of the insolation quantities of Mars[END_REF]. If rock breakdown occurred over similar timescales to obliquity cycles then, midlatitude craters should have a similar frequency of rockfalls to equatorial craters because obliquity has changed from 15 to 45° in the last 10 Myrs [START_REF] Laskar | Long term evolution and chaotic diffusion of the insolation quantities of Mars[END_REF], meaning the location of maximum insolation transfers between the equator and the mid-latitudes. Assuming a relatively constant obliquity in the last 100 kyrs, reversed perihelion would not change insolation conditions to the point where a switch of maximum solar flux between north and south facing slopes would occur in the mid-latitudes, as it does at the equator, so no change in the rockfall distribution would be expected.

10-Conclusion

We report on the first detailed study of individual recent rockfalls in impact craters on Mars inferred from the presence of tracks on the crater walls. We observe that the frequency of rockfalls and their orientation is dependent on latitude. The frequency of rockfall tracks is higher on the equator-facing slopes at the latitudes between 15°N/S and 40°N/50°S. In equatorial impact craters, the frequency of rockfall tracks is higher on both the north-and south-facing slopes compared to the east-and west-facing slopes (Fig. 4). Median clast size tends to decrease towards higher latitudes (Fig. 6). Topographic analysis shows that the signals observed are not a direct consequence of systematic variations in rock wall slopes (Fig. 7). Thus, these trends are more likely to be linked to the weathering mechanism responsible for rock breakdown prior to the rockfall rather than the slope inclination (Figs. 1). In addition, the observed patterns argue against a role of H 2 O or CO 2 phase changes in preconditioning slopes for rockfalls, considering that most rockfalls occur in the equatorial area where these volatiles are scarce or lacking. Instead, thermal stress-driven subcritical cracking [START_REF] Collins | Thermal influences on spontaneous rock dome exfoliation[END_REF][START_REF] Eppes | Mechanical weathering and rock erosion by climate-dependent subcritical cracking[END_REF] related to high contrasts in surface temperature is more likely to be responsible for rock breakdown on modern Mars. Comparison between our results and a solar flux model (Fig. 9) emphasizes the potential role of diurnal temperature cycles at preconditioning slopes for rockfalls. We suggest that thermal stress must be developed over timescales long enough for full cycles of precession to occur (Fig. 8) in order to explain the bimodal peaks in rockfall frequency at the equator.

Our study shows the key role of thermal stress in rock-breakdown on Mars. Thermal stress should thus be considered as an important factor in promoting mass-wasting process on impact crater walls.. Fig. 4. Distribution of rockfall tracks by orientation in impact craters within different latitude ranges, derived using two different methods. All directions are with respect to the crater centers (e.g. North in the roses refers to northern slopes, i.e. South facing, arrows indicate equator direction). Left: non-normalized distributions, where the number of rockfall tracks in each orientation bin is summed for every crater in the specific latitude range and expressed as a percentage of the total number of rockfall tracks in this latitude range. This method emphasizes the signal from craters which have a large rockfall population (>100). Right: normalized distribution where the percentage of rockfalls is calculated for each orientation bin in each crater from which the mean is then derived for all craters in each latitude range. This method emphasizes the signal from craters having a relatively low rockfall population (< 20). Lines correspond to vector means with 95% confidence interval depicted by brackets. Rao's spacing test and Rayleigh statistical test were performed for the null hypothesis of uniformity. Since equatorial plots display two distinct trends, we calculated two vector means for each half of the plots, assuming the other half is uniform. Volume is calculated from diameter assuming an elliptical-shaped object with an aspect ratio of 0.8 [START_REF] Kumar | Recent seismicity in Valles Marineris, Mars: Insights from young faults, landslides, boulder falls and possible mud volcanoes[END_REF] and is normalized to the maximum volume recorded Power-law fit is shown. . Rockfall track frequency compared with slope derived from eleven digital terrain models. Slope is calculated for a 150 m wide area starting 50 m below crater rim, which is the most probable boulder source. Craters are divided into 20° bins in which rockfalls are counted. A slope value is attributed to each rockfall track, corresponding to the slope in the bin where it is located. A: Frequency of rockfall tracks against slope value. Rockfall tracks occur more frequently on slopes steeper than 32°. B: Frequency of all slopes with and without rockfall tracks. In both, the distribution is normal, but with different medians. C: Frequency distribution of slope angles for slopes without any observed rockfall tracks. 100% of slopes between 18 and 20° are devoid of rockfalls because they are not steep enough (Fig. 1). 28% of slopes between 36 and 38° do not display any rockfall either. B and C highlight that a steep slope alone is not enough for a rockfall to occur, and that a stress-inducing mechanism is also necessary to allow cohesion loss of the material (Fig. 1). Fig. 8. Insolation received by a north-or south-facing slope located at the equator in two opposite cases. At 45.7 ka, the longitude of perihelion was very similar to the present day (Ls =251°), south-facing slopes received maximum insolation. At 22.7 ka at the opposite sense of perihelion, north-facing slopes received maximum insolation at the equator. Fig . 9. Year maximum diurnal insolation plotted with normalized rockfall frequency in northern, southern hemisphere, as well as equatorial region. Average diurnal insolation is computed for the latitude of each crater in a given range (from 50°S to 15°S and 15°N to 50°N). The model runs for a whole martian year every 45° of solar longitude of perihelion, to have an average over a full precession cycle. The solar data is then binned the with same method used for rockfall tracks (Section 3 -Fig. 4) and an average insolation value is obtained for the given latitude range. Solar data bins are compared to their corresponding normalized rockfall frequency bins (Fig. 4). P-value for the null-hypothesis is <<0.01 in both plots.

Fig. 5 .

 5 Fig. 5. Magnitude-cumulative frequency (MCF) relationship of clasts derived from martian rockfalls.

Fig. 6 .

 6 Fig. 6. Median long axis size of the rockfall clasts with respect of the latitude in each crater where recent rockfalls were mapped. The numbers next to each bar correspond to the number of clasts recorded in each crater. The error bars represent the standard errors. The craters with larger populations have lower standard errors. Craters are sorted by region to account for lithology variations, as different rock compositions would induce different thermal proprieties and thermal stress response.

Fig. 7

 7 Fig. 7. Rockfall track frequency compared with slope derived from eleven digital terrain models. Slope is

Fig

  Fig. S1. Example of slope extraction from the CTX DTMs. Each point corresponds to a pixel value from the DTM in a 20° bin of azimuth with their distance from crater slope. The linear fit used to calculate the slope is shown in red together with the derived slope angle.

Fig. S2

 S2 Fig. S2 Median clast size of rockfalls with respect to the host crater diameter

Acknowledgements

D. Mège, P-A. Tesson, and J. Ciazela are supported by the TEAM program of the Foundation for Polish Science (TEAM/16-3/20) co-financed by the European Union under the European Regional Development Fund. J. Ciazela is additionally supported within the START program of the Foundation for Polish Science. P-A Tesson, N. Mangold and S.J. Conway are also supported by CNES (French Space Agency). S.R. Lewis gratefully acknowledges the support of the STFC/UK Space Agency under grants ST/R001405/1 and ST/S00145X/1.

Appendix A1 -Solar Radiation Model

Our model computes daily mean solar insolation for every pixel of the input Digital Terrain Model (DTM), at a specific sol/Solar longitude (Ls). Each pixel has a slope and an aspect value as well as information about visibility in every direction ("viewshed"). Viewshed is computed for every pixel using Skyline_3d and SkylineGraph_3d function of ArcPy Python package for ArcGIS [START_REF] Rich | Using Viewshed Models to Calculate Intercepted Solar Radiation: Applications in Ecology[END_REF]. Equations 1-8 making up the mathematical framework of our model are all extracted from Appelbaum & Flood (1989).

Sun map calculation

Each sol is split in 48 equal timesteps, for which we calculate solar elevation angle (α s ) (Eq. 1) and Sun azimuth angle (Eq. 2). ϕ

(1) sin = cos cos cos + sin sin

(2) cos ϕ = sin cos -cos cos sin sin ( ) where is the hour angle (Eq. 3), the solar declination (Eq. 4), is the latitude, and is the solar zenith angle ( = 90° -). When < 0, Sun azimuth angle should be subtracted of 180°. Hour angle can be obtained with the following equation:

(3) = 15

-180 where LST is the local solar time in Martian hour. 15° is the rotation speed of Mars per Martian hour. Solar declination angle can be calculated using the following equation:

(4) sin = sin 0 sin where is Mars obliquity and Ls is the solar longitude. Comparison of Sun position in the sky at each 0 step and visibility provided by the viewshed assess whether solar flux should be calculated or not.

Direct irradiance calculation

Solar constant S (W/m²) is the solar flux received at the top of the atmosphere (Eq. 5).

(5)

where ecc is the eccentricity, Ls the solar longitude, solar longitude of perihelion and the mean solar constant (for Mars: 586 W/m²). Direct solar radiation D (W/m²) is the fraction of solar flux reaching the surface, and it is given by Eq. 6:

(6) = -( ) cos where is the optical depth of martian atmosphere, i the incidence angle (Eq. 8) and the airmass ( ) which can be approximated by Eq. 7:

(7)

Slope correction

Incidence angle of the solar flux with respect to a tilted surface is calculated (Eq. 8) and included in Eq. 6.

(8) cos = cos cos + sin sin cos ( -) where is the slope inclination and is the slope orientation.

Fig. 1. Schematic diagram illustrating how thermal stress could influence rockfall activity. Firstly, a topographic gradient is needed. The slope value and orientation control the amount of insolation received by the surface, and hence the potential solar-induced thermal stress intensity. Energy is mostly provided by gravity, although marsquakes may also contribute (e.g. [START_REF] Roberts | Possible evidence of paleomarsquakes from fallen boulder populations, Cerberus Fossae, Mars[END_REF]. Crack formation mechanism is required to weaken the exposed material and reduce its cohesion.