The Timing of Alluvial Fan Formation on Mars

River deposits on Mars record past climatic conditions and thus serve as a proxy for Mars's climate evolution toward 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, 2019; Hauber et al. 2013; Kite et al. 2017; 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. 2021). The apparently recent ages of these fluvial features are at odds with the belief that intense atmospheric loss early in Mars's history would have rapidly rendered Mars's 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 therefore a cold climate) in order to overpressurize 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 (Kite et al. 2021; Wordsworth et al. 2021), and obliquity-shift-induced snowmelt/precipitation (e.g., Irwin et al. 2015; Mansfield et al. 2018). 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₂, 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 (Kite et al. 2019), at least some formed over an extended time span (>20–200 Myr) (Kite et al. 2017), and some are up to ∼1 km thick, recording a minimum of 100 yr–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 (Wilson et al. 2021), 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 that are, by virtue of their size, more susceptible to survey complications like atmospheric filtering (e.g., Kite et al. 2014; Warren et al. 2019) 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 (Wilson et al. 2021), where erosional backwasting of crater rims provided alluvium that was transported downslope 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 order-of-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 (Wilson et al. 2021), 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 data sets will affect fan-timing inferences.

To investigate the formation of alluvial fans, we used a recently compiled global database of alluvial fans (including potential deltas) in craters (Wilson et al. 2021) that incorporates both previous surveys (Moore & Howard 2005; Kraal et al. 2008; Morgan & Wilson 2019), and also fans newly identified 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 fan data set 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 (Figure A1 in the Appendix). 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.

Zoom In Zoom Out Reset image size

To obtain crater densities for each AHi subunit'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 (S. Robbins 2021, personal communication), we restricted the data set 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 predate 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 owing 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 predate crater ejecta based on cross-cutting relationships. In general, we were conservative and did not discard ambiguous craters, preventing us from erroneously underestimating the age of the subunits. This procedure, along with removing craters within each AHi subunit'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 subunit-forming-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 two 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., predating 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.

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, ${ \mathcal L }$, given by the standard equation:

$$\begin{eqnarray}&&{ \mathcal L }(t| k,A)=P(k| t,A)=\displaystyle \frac{{\left(A{\lambda }_{t}\right)}^{k}}{k!}{e}^{-A{\lambda }_{t}},\end{eqnarray} \tag{ 1 }$$

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; Neukum et al. 2001):

$$\begin{eqnarray}\begin{array}{rcl}{\lambda }_{t}(D) & = & C(D)\left[(3.79\times {10}^{-14})({e}^{6.93t}-1)\right.\\ & & \left.+(5.84\times {10}^{-4})t\right],\end{array}\end{eqnarray} \tag{ 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 25th and 75th age percentiles. These corresponded to their final age estimates and error bounds, respectively.

In our framework, we had strong reason to believe a nonuniform 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's 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 Section 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):

$$\begin{eqnarray}\begin{array}{rcl}P(t) & \propto & \displaystyle \frac{d}{{dt}}\left[(3.79\times {10}^{-14})({e}^{6.93t}-1)\right.\\ & & \left.+(5.84\times {10}^{-4})t\right].\end{array}\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

Using the likelihood function and prior distribution, we computed posterior age pdf's (Figure 2) by multiplying and normalizing:

$$\begin{eqnarray}&&P(t| k,A)\propto { \mathcal L }(t| k,A)P(t).\end{eqnarray} \tag{ 4 }$$

Our estimates are consistent with the method of Michael et al. (2016) in ∼85% 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).

Zoom In Zoom Out Reset image size

Using the posterior age CDFs for each AHi subunit, we derived the CDF for the age of the youngest fan-bearing AHi subunit. To do so, we assumed that the ages of each AHi subunit are independent random variables:

$$\begin{eqnarray}&&P({t}_{\min }\geqslant {t}_{\mathrm{ref}}| {t}_{\mathrm{ref}})=\prod _{i}P({t}_{i}\geqslant {t}_{\mathrm{ref}}| {k}_{i},{A}_{i},{t}_{\mathrm{ref}}),\end{eqnarray} \tag{ 5 }$$

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 }\lesssim 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 }\lesssim 2.5$ Ga.

Zoom In Zoom Out Reset image size

3.2. 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 (pulse-formation model), we assumed that alluvial fans formed at a single moment in history, τ_form, and that all AHi subunits with age t_i ≥ τ_form had an equal probability, p_max, of acquiring alluvial fans. In the second model (prolonged-formation model), each AHi subunit had a constant instantaneous probability, ${p}_{\max }/{\rm{\Delta }}\tau$, 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 τ_form 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.

Zoom In Zoom Out Reset image size

For both statistical models, the formation scenarios can be viewed both in an instantaneous framework (left panels of Figure 5) and in a cumulative framework (right panels of 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:

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal L }({\tau }_{\mathrm{form}},{p}_{\max }| {F}_{i}=1,{t}_{i})\\ \ \ =\ {\displaystyle \int }_{0}^{{t}_{i}}P(F\to 1| {\tau }_{\mathrm{form}},{p}_{\max },\tau )d\tau \end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal L }({\tau }_{\mathrm{end}},{\rm{\Delta }}\tau ,{p}_{\max }| {F}_{i}=1,{t}_{i})\\ \ \ =\ {\displaystyle \int }_{0}^{{t}_{i}}P(F\to 1| {\tau }_{\mathrm{end}},{\rm{\Delta }}\tau ,{p}_{\max },\tau )d\tau ,\end{array}\end{eqnarray} \tag{ 7 }$$

where F_i = 1 when the ith AHi subunit contains fans (0 otherwise) and P(F → 1) is the probability of acquiring fans, which varies with time (Figure 5). Note that, because F_i is a binary variable, P(F_i = 0) = 1 − P(F_i = 1), which defines ${ \mathcal L }$ in the case in which 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:

$$\begin{eqnarray}&&\begin{array}{l}E({ \mathcal L }({\tau }_{\mathrm{form}},{p}_{\max }| {F}_{i},{k}_{i},{A}_{i}))\\ \ \ =\ {\displaystyle \int }_{0}^{4.5}{ \mathcal L }({\tau }_{\mathrm{form}},{p}_{\max }| {F}_{i},t)P(t| {k}_{i},{A}_{i}){dt}\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\begin{array}{l}E({ \mathcal L }({\tau }_{\mathrm{end}},{\rm{\Delta }}\tau ,{p}_{\max }| {F}_{i},{k}_{i},{A}_{i}))\\ \ \ =\ {\displaystyle \int }_{0}^{4.5}{ \mathcal L }({\tau }_{\mathrm{end}},{\rm{\Delta }}\tau ,{p}_{\max }| {F}_{i},t)P(t| {k}_{i},{A}_{i}){dt},\end{array}\end{eqnarray} \tag{ 9 }$$

where k_i is the observed crater count and A_i the count region area for the ith (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:

$$\begin{eqnarray}&&\begin{array}{l}E\left({ \mathcal L }({\tau }_{\mathrm{form}},{p}_{\max }| {\boldsymbol{F}},{\boldsymbol{k}},{\boldsymbol{A}})\right)\\ \ \ =\ \prod _{i}E\left({ \mathcal L }({\tau }_{\mathrm{form}},{p}_{\max }| {F}_{i},{k}_{i},{A}_{i})\right)\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\begin{array}{l}E\left({ \mathcal L }({\tau }_{\mathrm{end}},{\rm{\Delta }}\tau ,{p}_{\max }| {\boldsymbol{F}},{\boldsymbol{k}},{\boldsymbol{A}})\right)\\ \ \ =\ \prod _{i}E\left({ \mathcal L }({\tau }_{\mathrm{end}},{\rm{\Delta }}\tau ,{p}_{\max }| {F}_{i},{k}_{i},{A}_{i})\right).\end{array}\end{eqnarray} \tag{ 11 }$$

We computed the expected likelihood for the pulse-formation model as a function of τ_form and p_max, allowing them to vary in [0 Ga, 4 Ga] and [0, 1], respectively (Figure 6). The model achieved maximum likelihood at τ_form = 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 τ_form (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 rightmost 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 require fan formation to occur within the last ∼2.8 Gyr.

Zoom In Zoom Out Reset image size

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 data set requires fan formation to have persisted into the last ∼2.3 Gyr.

Zoom In Zoom Out Reset image size

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 crater counts 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 τ_form 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 may 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 τ_form 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 owing 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 }\lt 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 data set. 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, 2019) and are also consistent with Amazonian age estimates made for some nonfan 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 do 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 τ_form and τ_end, which had large uncertainties associated with them. Indeed, our crater counts were sparse owing 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 owing 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 τ_form = 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).

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 pulse-formation 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 that 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).

Zoom In Zoom Out Reset image size

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; Farley et al. 2014; Robbins 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 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 reperformed 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.

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 (Wilson et al. 2021), 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 data set 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's climate.

We thank Stuart Robbins and Marisa Palucis for insightful discussions that helped conceive this study. We thank two reviewers for useful comments. We acknowledge funding from NASA grants NNX16AJ38G, NNX15AM49G, and 80NSSC20K0144. The scripts used to make this paper can be obtained for unrestricted further use by contacting the senior author (E.S.K.).

Figure A1 shows the map of fans and craters.

Zoom In Zoom Out Reset image size