The Timing of Alluvial Fan Formation on Mars

The history of rivers on Mars is an important constraint on Martian climate evolution. The timing of relatively young, alluvial fan-forming rivers is especially important, as Mars' Amazonian atmosphere is thought to have been too thin to consistently support surface liquid water. Previous regional studies suggested that alluvial fans formed primarily between the Early Hesperian and the Early Amazonian. In this study, we describe how a combination of a global impact crater database, a global geologic map, a global alluvial fan database, and statistical models can be used to estimate the timing of alluvial fan formation across Mars. Using our global approach and improved statistical modeling, we find that alluvial fan formation likely persisted into the last ~2.5 Gyr, well into the Amazonian period. However, the data we analyzed was insufficient to place constraints on the duration of alluvial fan formation. Going forward, more crater data will enable tighter constraints on the parameters estimated in our models and thus further inform our understanding of Mars' climate evolution.


INTRODUCTION
River deposits on Mars record past climatic conditions and thus serve as a proxy for Mars' climate evolution towards its current cold and dry state (e.g. Kite 2019). Since the end of the era of regionally-integrated valley networks in the early Hesperian ∼ 3.4 Ga (Fassett & Head 2008), fluvial activity on Mars primarily produced shorter landforms (Goudge et al. 2016), including small, 'pollywog', exit breach craters (Wilson et al. 2016;Warren et al. 2021), alluvial fans (Moore & Howard 2005;Kraal et al. 2008;Grant & Wilson 2011;Hauber et al. 2013;Kite et al. 2017;Grant & Wilson 2019;Wilson et al. 2021), and gullies (which may or may not be fluvial in origin; see Conway et al. (2019)). Some of the alluvial fans are associated with aqueous minerals (hydrated silica; e.g. Pan et al. (2020)). The apparently recent ages of these fluvial features are at odds with the belief that intense atmospheric loss early in Mars' history would have rapidly rendered Mars' surface cold and dry (e.g. Kite et al. 2019). Further, Amazonian/Hesperian catastrophic floods (e.g. Baker & Milton 1974) are commonly (but not exclusively) believed to require a ≥ 1 km thick cryosphere (and thereby a cold climate) in order to over-pressurize the aquifers that fed these floods (e.g. Carr 1979). Multiple mechanisms have been proposed to explain late-forming rivers, including impact-induced snowmelt/precipitation (e.g. Williams & Malin 2008;Mangold et al. 2012b), transient greenhouse atmospheres (e.g. Kite et al. 2021;Wordsworth et al. 2021), and obliquity-shift induced snowmelt/precipitation (e.g. Irwin et al. 2015). Each of these scenarios predicts transient excursions from an otherwise cold, dry climate state. Different proposed mechanisms have different dependencies on volcanism and/or pCO 2 , both of which declined over time (Fassett & Head 2011). Thus, the timing of young fluvial activity is an important constraint on Martian climate evolution.
In this study, we focus on the alluvial fans (for the purposes of this study, we also include deltas, which are fans formed in standing bodies of water) and how improved statistical methodologies shed light on their formation ages. The alluvial fans are of particular interest because many are precipitation-fed , at least some formed over an extended time span   (Kite et al. 2017), and some are up to ∼ 1 km thick, recording a minimum of 100 years -1 Myr of intermittent river flow (Irwin et al. 2015;Stucky de Quay et al. 2019). These observations, in addition to the widespread distribution of alluvial fans within observed latitude bands , are consistent with a synoptic water source rather than localized impact trigger and therefore provide true constraints on the Martian climate through time.
The primary method for estimating the absolute age of geologic surfaces on Mars is through counting the craters in an area and comparing the observed density to a chronology function, which relates crater density to absolute age (e.g. Michael 2013;Fassett 2016, and references therein). Dating small areas like fan surfaces requires counting small craters which are, by virtue of their size, more susceptible to survey complications like atmospheric filtering (e.g. Kite et al. 2014) and geologic modification (e.g. erosion, burial; Palucis et al. (2020)). Indeed, derived ages for small areas such as individual deltas vary over a huge time range (e.g. Hauber et al. 2013;Warner et al. 2015). To reduce these uncertainties, previous workers have averaged over numerous close-proximity fans (e.g. Moore & Howard 2005;Grant & Wilson 2011;Morgan & Wilson 2019). However, this methodology (1) implicitly assumes that fan formation occurred simultaneously for each fan (such that their surfaces reflect a single age), (2) requires that there has not been sufficient, potentially spatially-varying, resurfacing of the fan surfaces to obscure any craters in the considered size range (Palucis et al. 2020), and (3) sets aside the likely possibility that alluvial fans formed intermittently over tens or hundreds of Myr, rather than instantaneously (Kite et al. 2017;Stucky de Quay et al. 2019). We propose that these shortcomings can be vaulted by a new approach involving careful consideration of the fans' cross-cutting relationships and choice of statistical modeling methodology.
The majority of alluvial fans are hosted on the interior walls of impact craters , where erosional backwasting of crater rims provided alluvium that was transported down-slope and deposited on crater floors. These fans unambiguously post-date their host-craters, whereas determination of alluvial fan stratigraphic relationships within host crater interiors requires more detailed investigation on an individual fan basis. Further, the ages of the impact craters can be estimated via crater counts on their ejecta deposits (where preserved), which provide us orderof-magnitude greater count areas (and thus the ability to count larger craters) than the fan surfaces themselves. In the following sections, we demonstrate how combining a global database of impact craters (Robbins & Hynek 2012), a global geologic map (Tanaka et al. 2014), a global database of alluvial fans , and a novel statistical model can constrain alluvial fan formation timing to have persisted until < 3 Ga. Further, we will discuss limitations of our approach and how updated datasets will affect fan-timing inferences.

DATASETS AND METHODS
To investigate the formation of alluvial fans, we used a recently compiled global database of alluvial fans (including potential deltas) in craters ) that incorporates previously published surveys (Moore & Howard 2005;Kraal et al. 2008) and abstracts (e.g. Morgan & Wilson 2019) in addition to newly identified fans using images from the Context Camera (CTX) (Malin et al. 2007). Although the use of CTX imagery enabled surveying smaller fans than previous studies, we expect that the smallest features are missing from the survey, due to the finite resolution of the images used. For this study, we restricted the fans dataset to those contained within the 212 Amazonian-Hesperian Impact Unit (AHi) subunits from the Tanaka et al. (2014) global geologic map of Mars that are fully contained within 40 • latitude of the equator (Appendix Figure 1). These subunits are impact craters that have well-preserved rims and ejecta (Tanaka et al. 2014).These AHi subunits provide us a sample of areas that are globally distributed spatially, share similar steep slopes and topography, and record their ages via crater-counts on their ejecta. Because fans on the crater walls must post-date the crater formation, ages obtained from the ejecta provide robust upper bounds on each fan's age. Further, not all AHi subunits possess fans (70 out of 212 subunits contain fans), which enables comparison between the formation times of fan-bearing and non-fan-bearing surfaces.
To obtain crater densities for each AHi subunits's ejecta, we used a global database of impact craters > 1 km in diameter (Robbins & Hynek 2012). However, due to concerns over survey completeness for small craters (Stuart Robbins, personal communication), we restricted the dataset to craters with diameters of at least 4 km. Manual inspection of craters on the AHi subunits revealed three additional potential data quality issues. First, ejecta blankets often contain easily-identified chains of secondary craters. Second, some craters mapped on ejecta actually pre-date ejecta, as evidenced by cross-cutting relationships with the central crater's rim, impact ejecta, or secondary craters ( Figure 1). Finally, craters in the interior of AHi subunits' central craters are often obscured from surveys due to resurfacing by burial/erosion.
To remedy these crater-counting issues, we manually inspected each crater in Thermal Emission Imaging System (THEMIS) (Christensen et al. 2004) and CTX imagery and removed obvious secondary craters (elongated craters often in chains extending from the center of large nearby impact craters) and craters that obviously pre-date crater ejecta based on cross-cutting relationships. In general, we were conservative and did not discard ambiguous craters, preventing us from erroneously under-estimating the age of the subunits. This procedure, along with removing craters within each AHi subunits's central crater's interior, reduced our total crater count from 1768 to 1308. Further, we restricted our count areas on each AHi subunit to the annulus extending from the crater rim to 2× the subunitforming-crater's radius (i.e., we restricted our count areas to 1-2 crater radii from the crater center). In the < 5% of cases where the AHi subunit is defined by the ejecta blankets of 2 large craters, we chose the larger crater as our nominal center. In many cases, the ejecta blankets do not extend to 2× the crater radius, in which case we included in our count area only the intersection of the ejecta blanket and our defined annulus. We imposed the outer limit on our count area to accommodate the fact that craters "poking through" (i.e. pre-dating ejecta) near the ejecta edges cause crater densities to increase by a factor of ∼ 2-3 but that unambiguous cross-cutting relationships are increasingly difficult to identify as the ejecta thins radially outward. This indicates that the fraction of craters "poking-through" the ejecta blankets increases near the edges, and thus, the outer limit prevents erroneously biasing our AHi subunit ages upward. Imposing the outer annulus limit further reduced our total crater count to 498.  (Dickson et al. 2018) (illumination from north-west to south-east) image of a large crater-hosted alluvial fan deposit, formed by erosional backwasting of the southwestern crater rim (333 • E, 28 • S). Right: Colorized elevation, ranging from ∼ −1.8 km (indigo) to ∼ 0.7 km (yellow), over THEMIS daytime IR image of an AHi subunit (0 • E, 15 • S) (illumination from west to east). In this example two craters (C, D) pre-date the central crater (A) based on cross-cutting relationships with the central crater's northern rim and ejecta deposits. The crater in the top right (E) has well-preserved ejecta that cross-cuts the central crater's ejecta deposits, indicating that it post-dated the AHi subunit formation. The crater in the bottom left (B) appears fresh and has ambiguous relationships with the terrain around it, so it was included in the count. In both images, blended High Resolution Stereo Camera (HRSC) and Mars Orbiter Laser Altimeter (MOLA) 200m elevation data (Fergason et al. 2018) was used to produce colorized elevation. In both images, north is towards the top of the page.

STATISTICAL MODELING AND RESULTS
Our goal was to take crater counts on the AHi subunits' ejecta blankets, note which subunits host alluvial fans, and make quantitative statements about the timing and duration of alluvial fan formation. We achieved this in two steps. First, we developed a probabilistic understanding of the age of each AHi subunit by estimating their posterior age distributions in a Bayesian framework (similar to Michael et al. 2016). Second, we implemented two statistical models that represent hypothetical formation scenarios for the alluvial fans. Parameter estimation for these models enabled inferences about the timing of alluvial fan formation.
3.1. Dating the AHi subunits Michael et al. (2016) described the use of Bayesian posteriors for estimating planetary surface ages, under the name "Poisson Timing Analysis." In this framework, crater counts were assumed to be generated as a Poisson process, with underlying parameter likelihoods, L, given by the standard equation: where k is the observed number of craters, A is the area of the count region, and λ t is the expected density of craters on a surface with age t. The correspondence between age and expected density is known as the chronology function, and throughout this study, we adopt the chronology function of Michael (2013) (based on Hartmann (2005) and Neukum et al. (2001)): ( 2) where λ t (D) is the expected density of craters > D km in diameter, and C(D) is a function that accounts for the size-frequency distribution of craters (Table 1 in Michael (2013)). If, a priori, each candidate age is assumed to be equally likely, the likelihood function can be normalized (scaled such that its integral equals 1) and treated as a posterior probability distribution for subunit age. From there, Michael et al. (2016) computed median ages, as well as the 25% and 75% age percentiles. These corresponded to their final age estimates and error bounds, respectively. Example normalized likelihood functions, as well as prior and posterior distributions, for the ages of two AHi subunits. For the subunits in the left panel, no craters were counted on the AHi unit ejecta (none > 4 km in diameter were found within the count annulus), hence the maximum likelihood is achieved at t = 0 Ga. For the subunits in the right panel, two craters were counted on the AHi subunit ejecta, hence the maximum likelihood is achieved at t ∼ 3.5 Ga. Prior distributions are assumed to be the same for each AHi subunit. Despite different likelihood functions, the posterior distributions are similar.
In our framework, we had strong reason to believe a non-uniform prior distribution in AHi subunit formation age. In particular, impact fluxes on Mars were elevated > 3 Ga, relative to today (see among many others, Michael (2013)). This implies that the AHi subunits, which are themselves formed by impact craters, are more likely to have formed early in Mars' history. Thus, for our application, we took the prior cumulative distribution function (CDF) of AHi subunit ages to be proportional to the chronology function, truncated at 4.5 Ga. A sensitivity test truncating at 3.54 Ga is discussed in §4.). We then computed the prior probability density function (PDF) on AHi subunit ages as the normalized derivative of the Michael (2013) chronology function ( Figure 2): Using the likelihood function and prior distribution, we computed posterior age PDF's ( Figure 2) by multiplying and normalizing: Our estimates are consistent with the method of Michael et al. (2016) in ∼ 90% of cases ( Figure 3). In every case, our strong prior forces our central estimate to be older than the central estimate from the Michael et al. (2016) method ( Figure 3). Almost all AHi subunits are consistent at the 95% confidence level with being Hesperian or Amazonian in age, however (considering only crater data and setting aside cross-cutting relationships) most of the AHi subunits are also consistent with being Noachian in age (e.g. Fassett 2016). Using the posterior age CDF's for each AHi subunit, we derived the CDF for the age of the youngest fan-bearing AHi subunit. To do so, we assumed the ages of each AHi subunit are independent random variables: where i is an arbitrary index over the fan-bearing AHi subunits. We computed this distribution for two populations: (a) all fan-bearing AHi subunits (N = 70), and (b) those fan-bearing AHi subunits where localized impact trigger is excluded as an explanation for alluvial fan formation (N = 13). Population (b) comprises those fan-bearing subunits in which pre-or syn-fluvial impact craters were identified (Kite et al. 2017;Kite 2021). The results are summarized in Figure 4. We found at the 95% confidence level that t min 0.9 Ga if we include all fan-bearing AHi subunits.
Considering only fan-bearing AHi subunits where localized impact trigger is excluded, we found at the 95% confidence level that t min 2.5 Ga. The blue curve is the CDF for all fan-bearing subunits, while the red curve is the CDF for those fan-bearing subunits where localized impact trigger is excluded. The horizontal, black dashed line indicates the P = 0.05 level. Thus, the red and blue vertical dashed lines (defined by the intersection of the CDF's with the P = 0.05 level) represent our 95% confidence level upper estimates for the age of the youngest fan in each population. Our results indicate that at least one fan likely formed in the last ∼ 0.9 Ga and that at least one fan where localized impact trigger is excluded formed in the last ∼ 2.5 Ga.

Modeling Fan Formation Scenarios
To investigate the formation timing of alluvial fans in a way that incorporates information from crater-counts on non-fan-bearing AHi subunits, we developed two simple statistical models ( Figure 5). In the simplest model (pulseformation model), we assumed that alluvial fans formed at a single moment in history, τ f orm , and that all AHi subunits with age t i ≥ τ f orm , had an equal probability, p max , of acquiring alluvial fans. In the second model (prolongedformation model), each AHi subunit had a constant instantaneous probability, p max /∆τ , of acquiring alluvial fans during the ∆τ -long time period ending at τ end . This prolonged-formation model is approximately equivalent to a model with many short bursts of activity that are evenly spaced or randomly spaced throughout the ∆τ -long time period ending at τ end . In both models, the probability of an AHi subunit younger than τ end or τ f orm acquiring fans was assumed to be 0. Both models are motivated by searching for changes in the probability of acquiring fans as a function of time, for example due to atmospheric evolution. Neither model is set up to identify preferred spatial patches for fan formation, for example latitude dependent climate variations. Any such variations are averaged out in the globally uniform (but time-varying) probability of acquiring alluvial fans.
For both statistical models, the formation scenarios can be viewed both in an instantaneous framework (left panel, Figure 5), and in a cumulative framework (right panel, Figure 5). Both viewpoints are instructive, as the instantaneous view relates to a hypothetical climate scenario (e.g. short-lived, global transient greenhouse), while the cumulative view relates fan-acquisition probabilities to AHi subunit ages, t i . The relationship between these two views can be formalized mathematically via simple integrals for both models: where F i = 1 when the i th AHi subunit contains fans (0 otherwise), and P (F → 1) is the probability of acquiring fans, which varies with time ( Figure 5). Note, because F i is a binary variable, P (F i = 0) = 1 − P (F i = 1), which defines L in the case that F i = 0. Because we only know the AHi subunit ages, t i , probabilistically, we computed expected likelihoods (instead of the likelihoods themselves) as a function of model parameters and crater density data by marginalizing over the age posterior distributions: where k i is the observed crater count and A i the count region area for the i'th (out of 212) AHi subunits. Taking the data from every AHi subunit together, we can multiply the expected likelihoods together to estimate the overall expected likelihood as a function of model parameters and data: We computed the expected likelihood for the pulse-formation model as a function of τ f orm and p max , allowing them to vary in [0 Ga, 4 Ga] and [0,1], respectively ( Figure 6). The model achieved maximum likelihood at τ f orm = 1.5 Ga and p max = 0.36. Assuming a uniform prior over p max , we marginalized the likelihoods over p max ("summed over p max ") to estimate a posterior PDF for τ f orm (right panel of Figure 6). The 95% confidence level for the estimated posterior PDF in the right panel of Figure 6 is slightly more restrictive than the right-most extent of the middle red contour in the right panel of Figure 6 -this is due to the marginalization over p max . At the 95% confidence level, we found that in the pulse-formation scenario, the data requires fan formation to occur within the last ∼ 2.8 Gyr. We computed the expected likelihood for the prolonged-formation model, as a function of τ end , ∆τ , and p max , allowing them to vary in [0 Ga, 4 Ga], [0, 4 Ga −τ end ], and [0,1], respectively (Figure 7). The model achieves maximum likelihood at τ end = 0 Ga, ∆τ = 3 Ga, p max = 0.36 (Figure 7, top left panel). We compared the maximum likelihoods from the two models (this comparison does not require assuming that either model is a good model). Despite a maximum likelihood solution with ∆τ > 0, the difference in maximum likelihoods between the prolonged-formation model and the pulse-formation model (equivalent to the prolonged-formation model when ∆τ = 0) was not sufficiently large to reject the pulse-formation model via a likelihood-ratio test (p < 0.81, where our computed test statistic, λ LR ≈ 0.06, is given by twice the difference in the maximum log-likelihoods for each model and is asymptotically chi-square distributed with 1 degree of freedom). Assuming uniform priors, we marginalized the likelihoods over ∆τ and p max to estimate a posterior distribution for τ end (Figure 7). At the 95% confidence level, we found that in the prolonged-formation scenario, our dataset requires fan formation to have persisted into the last ∼ 2.3 Gyr. Red contours show 1-, 2-, and 3-σ confidence regions in the marginalized, 2-D parameter spaces. Note, this is not equivalent to the projection of confidence regions in 3-D parameter space down to 2-D. The apparent simultaneous maxima of τ end = 0 and ∆τ = 0 in the bottom panels arise from the fact that we restricted the domains of τ end and ∆τ (greyed out region in top-left panel) and assumed uniform prior distributions over those domains. The bottom panels do, however, show that pmax is strongly determined by the data and can safely be marginalized out of the inferred posterior distribution (i.e. one can interpret the top-left panel without too much concern over the effect of prior distribution choices). Expected likelihoods are calculated using Equation 11.

DISCUSSION
Using the crater counts on AHi subunits' ejecta blankets, as well as some simple statistical models, we were able to place constraints on the timing of the formation of alluvial fans. Our estimate of the CDF for the age of the youngest fan indicated that at least one fan likely formed by a process other than localized impact trigger in the last ∼ 2.5 Gyr. Further, our CDF estimate indicates that at least one fan likely formed in the last ∼ 0.9 Gyr. This is broadly consistent with the results from the pulse-formation model (which incorporates information from cratercounts on non-fan-bearing AHi subunits), in which the best-fitting time of fan formation was ∼ 1.5 Gyr. This indicates that fan formation likely persisted well into the Amazonian epoch, consistent with Palucis et al. (2020) and Grant & Wilson (2019). Although the pulse-formation model does not favor formation occurring in the last ∼ 1 Gyr (due to the fact that very young values of τ f orm force non-fan-bearing AHi subunits to take on likelihood 1 − p max instead of 1 with too high a probability), it does not explicitly reject it. Further, at least some of the young fans can be explained by localized-impact trigger (e.g. Williams & Malin 2008;Goddard et al. 2014), which is not permitted in the pulse-formation model scenario.
We estimated parameters for a more complex model, the prolonged-formation model, and compared the results with the pulse-formation model. The introduction of a third parameter introduced a large degeneracy between τ end and ∆τ in our parameter retrieval (Figure 7). Although the model selected τ end = 0 Ga with fan-formation lasting several Gyr, a likelihood ratio test revealed that the maximum likelihood is not sufficiently higher in the prolonged-formation model than in the pulse-formation model to definitively justify the addition of the third parameter. Thus we cannot reject either the pulse-formation or prolonged-formation model. Further, we do not have tight constraints on τ end +∆τ , due to the degeneracy in our prolonged-formation model parameter retrieval. As a result, we cannot conclusively say, using this crater data alone, whether the alluvial fans formed in an era entirely separate from the valley networks (i.e. τ end + ∆τ 3.4 Ga with high confidence), or if the alluvial fans represent a long tail of continuous fluvial activity (i.e. τ end + ∆τ 3.4 Ga with high confidence).
For both models, we showed marginalized distributions for τ f orm and τ end , assuming uniform priors over the other parameters, p max and ∆τ . Unlike the prior distributions for the ages of the AHi subunits, we had no particular reason to choose these parameter prior distributions (although we note that marginalizing over p max is likely okay due to its strong determination by the data, see Figure 7). The marginalized distributions should thus be viewed only as condensed data visualizations, rather than true posterior distributions. Despite this, maximum likelihood estimates from both models indicate fan formation continued into the last ∼ 1.5 Gyr. This global age is substantially (∼ 2 Gyr) younger than existing quantitative estimates of the age of fan formation from a study of two regions (Mangold et al. 2012a). This can be explained by the fact that our model allows for spatially patchy fan formation (i.e. p max < 1) while the method of Mangold et al. (2012a) assumes uniform fan formation and thus would overestimate the young limit on alluvial fan formation if applied to our dataset. We note that the assumption that p max = 1 may be perfectly acceptable within the Mangold et al. (2012a) study areas. Our results are consistent with young age estimates obtained using smaller count areas for some fans (e.g. Grant & Wilson 2011;Grant & Wilson 2019), and are also consistent with Amazonian age estimates made for some non-fan features that indicate flowing surface liquid water (e.g., Dickson et al. 2009;Fassett et al. 2010;Howard & Moore 2011;Weitz et al. 2013;Ehlmann & Buz 2015;Brossier et al. 2021;Butcher et al. 2021). A K-Ar age of (2.12 ± 0.36) Ga obtained for aqueous minerals at Gale crater by the Curiosity rover (Martin et al. 2017) is consistent with our scenario but does not require it. Our strategy in this study does not incorporate the prior indications of young activity from these previous studies, and our crater data does not overlap with that used by these previous studies, therefore our results are independent.
By themselves, our data do not allow us to determine whether fan formation was intermittent (either synchronous between fan-forming regions, or regionally diachronous), or alternatively if water flow occurred in more-or-less every year in all regions where fans are present. Consensus on this important topic may have to await a future landed mission.
In general, data availability limited our parameter retrievals for τ f orm and τ end , which had large uncertainties associated with them. Indeed, our crater counts were sparse due to our lower diameter cutoff of 4 km, causing our AHi subunit age posteriors to be largely influenced by our choice of prior distribution and have high variance (Figures 2  and 3). We performed a series of sensitivity tests to illustrate how the model performs both with more data and with idealized data. First, we performed the parameter estimation with the same data, but with the outer annulus cutoff removed ( Figure 8). This had the effect of biasing our retrievals to older ages due to the larger influence of old craters "poking-through" the ejecta blankets near the edges (Figure 8). Indeed, the pulse-formation model achieved maximum likelihood at τ f orm = 1.7 Ga and p max = 0.36. Further, the increased quantity of data decreased the variance of our estimates, which is evident in the decreased degeneracy in the prolonged-formation model ( Figure 8).
Next, in order to evaluate how effectively the mathematical model can retrieve parameter values, we carried out a test using synthetic data. Specifically, we performed the parameter retrievals using the original count areas (including outer annulus limit), using synthetic data where the absence or presence of alluvial fans and ejecta crater counts were randomly generated assuming that the pulse-formation model was true, that the fans formed at 1.5 Ga with probability 0.5, and that AHi subunits have ages drawn from the prior distribution, truncated at 4 Ga. In general, the model performs well and achieves maximum likelihood within 1σ of the correct value ( Figure 9). However, there is still substantial variance in the results, and the model is unable to substantially break the degeneracy between the end of fan formation and the duration of the fan formation era, although it does achieve maximum likelihood at τ end = 2.1 Ga, p max = 0.47 and ∆τ = 100 Myr (Figure 9). This suggests that more craters are needed to constrain the formation timing.  Figure 4.6, as expected given increased contamination by "poke-through" craters. Right: expected likelihood for the prolonged-formation model, as a function of τ end and ∆τ , marginalizing over pmax. Greyed-out region was not modeled. Red contours show 1-, 2-, and 3-σ confidence regions in both panels (constructed in marginalized, 2-D space in the right panel). Figure 9. Model retrievals where synthetic data was created assuming the pulse-formation model is correct (i.e. ∆τ = 0). Left: expected likelihood for the pulse-formation model, as a function of τ f orm and pmax. Dashed red lines show true parameters from the synthetic data generation process. Right: expected likelihood for the prolonged-formation model, as a function of τ end and ∆τ , marginalizing over pmax. Greyed-out region was not modeled. Red contours show 1-, 2-, and 3-σ confidence regions in both panels (constructed in marginalized, 2-D space in the right panel).
We repeated the experiment with synthetic data, artificially increasing the AHi count areas by a factor of 10. This had the effect of placing roughly 10 times as many craters on each subunit, thereby decreasing the uncertainty in the posterior age distributions. As a result, the variance in our parameter retrievals shrunk (compare Figure 10 with Figure  9). Further, the parameter retrieval has a slightly increased ability to discriminate between the prolonged-and pulseformation models (relative to the synthetic data experiment with uninflated AHi subunit areas) and preferentially selects a fan formation era duration of 0 ( Figure 10). This indicates that in order to obtain better estimates for the timing and duration of alluvial fan formation, one needs substantially more data and must consider craters at smaller diameters (recall our lower cutoff was 4 km). Although we ignored the possibility that surface processes efficiently obliterated 4 km diameter craters in the post-Noachian, one must correct for crater obliteration to avoid underestimating age when smaller craters are considered (e.g. Palucis et al. 2020).
In our methodology, we did not account for uncertainty in the chronology function, which is continually being reevaluated from a combination of crater counting, geochronology, and solar system dynamics perspectives (e.g. Werner & Tanaka 2011;Robbins 2014;Farley et al. 2014;Marchi 2021). It is possible that our age estimates are systematically younger or older than the true ages, which would have biased our fan formation timing estimates. This potential effect is especially important in the data-limited case, where the prior strongly informs posterior age distributions (e.g. if Figure 10. Model retrievals where synthetic data was created assuming the pulse-formation model is correct (i.e. ∆τ = 0) and inflating the count area by a factor of 10. Left: expected likelihood for the pulse-formation model, as a function of τ f orm and pmax. Dashed red lines show true parameters from the synthetic data generation process. Right: expected likelihood for the prolonged-formation model, as a function of τ end and ∆τ , marginalizing over pmax. Greyed-out region was not modeled. Inflating the count area had the effect of decreasing the variance in our estimates. Red contours show 1-, 2-, and 3-σ confidence regions in both panels (constructed in marginalized, 2-D space in the right panel).
there were ample data, the posterior curves in Figure 2 would coincide with the normalized likelihood curves). A related point is that if we had imposed an exponentially decaying prior likelihood on fan formation time, then our fan age estimates would have been older. Further, we assumed a maximum age of 4.5 Ga in our prior age distribution. In principle, one could impose a younger maximum age based on stratigraphy to obtain tighter constraints. We re-performed our parameter retrieval assuming a maximum age of 3.54 Ga (approximately the Noachian/Hesperian boundary). However, roughly half of the AHi subunits achieved maximum posterior probability at 3.54 Ga, indicating that this prior was too restrictive. It is possible that with improved crater surveys, the data would produce maximum a-posteriori ages that are well into the Amazonian.
Finally, we did not account for the possibility that the fan catalog is incomplete. It is possible that some fans were wrongly identified and that other fans were missed in certain craters in the survey. To explore this effect, we performed a sensitivity test where the probability of an AHi subunit acquiring fans after the end of fan formation (or after the pulse of fan formation) was 0.001, instead of 0. This had a negligible effect on the outcome, and we concluded that survey incompleteness needs to be more severe (e.g. a 10% chance of missing all the fans in a crater) to substantially alter the results.

CONCLUSIONS
We combined a global geologic map (Tanaka et al. 2014), a global database of Martian impact craters (Robbins & Hynek 2012), a global database of alluvial fans , and statistical models to make inferences about the timing and duration of alluvial fan formation. Although differing in detail due to differing assumptions, the results of all our models suggest that fan formation had to have persisted into the last ∼ 2.5 Gyr (based on our most robust and conservative estimate, see Figure 4), well into the Amazonian time period. This estimate is ∼ 1 Gyr younger than previous estimates (e.g. Mangold et al. 2012a). Overall, our modeling procedure was not able to constrain the duration of fan formation. As a result, we were unable to determine from this dataset whether the alluvial fans formed in a long tail of continuous fluvial activity that overlapped with valley-network formation, or if they formed in an era entirely separate from the valley networks.
In general, our inferences were limited by a lack of crater count data, and our tests with synthetic data suggest that detailed surveys of smaller craters on AHi subunit ejecta blankets will enable much tighter constraints on fan formation. Going forward, we hypothesize that more data and more complex statistical models (i.e. ones with more parameters) for the population of alluvial fans (or other geomorphic features) will enable determination of both (a) their joint temporal and spatial evolution and (b) the relative importance of climate-triggered and localized impact-triggered fan formation. The method we have developed can also be applied to determining the ages of other crater-hosted features on Mars. Unraveling the evolution of alluvial fan (and/or other crater-hosted fluvial features') formation would thereby provide even stronger constraints on the evolution of Mars' climate. Figure 11. Appendix Figure. Map showing the distribution of AHi subunits used in this study (outlined in yellow, from Tanaka et al. (2014)), as well as the distribution of alluvial fans/deltas (red, from Wilson et al. 2021). Latitude/longitude grid spacing is 30 • . Study region extends from 40 • S -40 • N and from 0-360 • E. Background is Mars Orbiter Laser Altimeter shaded relief.