Limits on Runoff Episode Duration for Early Mars: Integrating Lake Hydrology and Climate Models

Fluvio‐lacustrine features on the martian surface attest to a climate that was radically different in the past. Since climate models have difficulty sustaining a liquid hydrosphere at the surface, multiple cycles of runoff episodes may have characterized the ancient Mars climate. A fundamental question thus remains: what was the duration of these runoff‐producing episodes? Here we use morphometric measurements from newly identified coupled lake systems, containing both an open‐ and a closed‐basin lake (n = 7). We combined hydrological balances with precipitation outputs from climate models, and found that breaching runoff episodes likely lasted 102−105 yr; other episodes may have been shorter but could not be longer. Runoff episode durations are model‐dependent and spatially variable, and no climate model scenario can satisfy a unique duration for all coupled systems. In the near future, these quantitative constraints on early Mars lake persistence may be tested through in situ observations from Perseverance rover.

• Coupled lake system morphologies were combined with climate model outputs to quantify upper and lower limits on runoff episode durations • Breaching runoff episodes lasted 10 2 -10 5 years depending on models and are spatially variable; other episodes could be shorter but not longer • Our constraints on lake persistence may be tested through in situ observations made by the Mars 2020 Perseverance rover in Jezero crater

Supporting Information:
Supporting Information may be found in the online version of this article.
an improved understanding of the spatio-temporal complexity of the hydroclimate, but also a shift toward viewing different timescale-deducing methods as complementary rather than contrasting.
Most geomorphic evidence suggests a wet climate lasting at least 10 kyr, and perhaps 100 Myr, up to the Noachian-Hesperian boundary (3.7 Ga), whether inferred from erosional (e.g., Barnhart et al., 2009;Craddock & Howard, 2002) or depositional (e.g., Armitage et al., 2011;Grotzinger et al., 2015;Schon et al., 2012) systems. In order to align these geomorphic constraints with hydroclimatic limits set by geochemical records and climate models, it is often proposed that surface liquid water was episodic, although the mechanism behind this episodicity remains uncertain (Kite, 2019;Wordsworth, 2016). An important point that goes into this hypothesis is that geomorphic timescales are cumulative, that is, they record the total amount of time taken to create a landform, whether a bedrock valley or fan deposit. As such, these studies typically rely on assumptions for runoff intermittency to calculate total durations (e.g., Buhler et al., 2014). If no hiatuses are considered, local depositional timescales could be substantially shorter (1-100 years; e.g., Jerolmack et al., 2004;Kleinhans et al., 2010). Further, calculated water availability rates, such as precipitation, runoff, and melting rates, vary over several orders of magnitude ( Figure 1); these estimates are not only sensitive to the methodology used, but also the spatial and temporal resolution. Clearly, there remains a need to shed further light on the uninterrupted availability of liquid water on the martian surface. Here we aim to provide new insights into this problem by addressing the following question: How can we constrain the duration of a runoff episode?
Our approach focuses on valley network-fed paleolakes, which provides an opportunity to assess individual surface runoff episodes: a discrete episode of time with net average positive runoff (i.e., runoff rate exceeds any losses; Supplementary Figure S1). Lake systems fall into one of two broad hydrological categories: openor closed-basins (Cabrol & Grin, 1999Fassett & Head, 2008;Goudge et al., 2015Goudge et al., , 2016. Open-basin lakes accumulated enough water to overflow and erode an outlet canyon (e.g., Goudge et al., 2019), whereas closed-basin lakes did not. As such, the presence or absence of an outlet canyon directly records whether a threshold event-lake overspill-was achieved in any single runoff episode ( Figure S1; Supplementary Text S1). Water input cannot be considered cumulative if separated by periods of water loss. We capitalize on this threshold relationship by investigating newly identified coupled lake systems, which contain both an open-and a closed-basin lake that are hydrologically connected (Figure 2a;  Note that (b) and (c) are not to scale and do not represent the large variability in basin sizes (see Figure S2). (d) Simple model for lake volume changes over time following Equations 2 and 3. Gray shaded region = permitted runoff episode duration between breaching ( B T ) and the maximum duration set by the volume of the closed-basin lake ( max T ). Here,

Paleolake Geometry Mapping
We investigated a subset of valley network-fed coupled paleolakes from Stucky de Quay et al. (2020), which, for both open-and closed-basin lakes contained within a coupled system, provide lake basin area, L A , lake basin volume, L V , and watershed area, W A ( Figure 2b). We used 7 coupled lake systems from this morphologic database (Supplementary Figure S2; Table S2) and measured an additional fourth parameter for all systems: lake volume remaining after the open-basin lake breached and drained, R V . This was done by identifying the highest closed contour in the basin before it spilled into the outlet canyon (e.g., Fassett & Head, 2008). For this we used the 100 m/pixel global daytime infrared mosaic (Edwards et al., 2011) from the Thermal Emission Imaging System (THEMIS; Christensen et al., 2004), and 463 m/pixel Mars Orbiter Laser Altimeter global gridded elevation data (MOLA; Smith et al., 2001). Subsequently, coupled systems were classified as either embedded (where the open-basin lake is contained within the watershed of the closed-basin lake watershed; n = 6) or adjacent (where the open-and closed-basin lake watersheds share drainage divides; n = 1). Importantly, sediment deposition is unlikely to significantly affect measured basin morphologies (see Supplementary Text S2; Mangold et al., 2009).

Derivation of Lake Hydrological Balance
Paleolake hydrology can be expressed using standard water balance equations (Horton, 1943;Fassett & Head, 2008;Matsubara et al., 2011). In a simplified system, the total lake volume, L V , accumulated over a runoff episode of duration, T , is given by where P is rainfall and/or snowmelt rates and E is evaporation rate, where both are averaged in space (across the total system area) and in time (across T ). Since lakes are fed by integrated valley networks ( Figure 1; Figure S2), this implies they were predominantly fed by surface runoff, and we suggest groundwater infiltration is likely to have limited effect on hydrology (Supplementary Text S3; Matsubara & Howard, 2009;Stephenson & Zuzel, 1981). In an embedded coupled system where an upstream open-basin lake (O) breaches at a time, B T , and overflows into a downstream closed-basin lake (C), we express both lake volumes, O v and C v , as a function of time, t: , where the full derivation is provided in Supplementary Text S4. This assumes a uniform P, although we find that variability is unlikely to significantly affect our results (Supplementary Text S5). Here, X denotes the aridity of the system, and can be expressed as 1 ( 1) AI  , where the aridity index, AI, is the ratio of runoff to evaporation ( / P E; Matsubara et al., 2011;Stucky de Quay et al., 2020). Because closed-basin lakes did not overflow, the observed basin volume provides an upper limit for C v , that is, C v in Equation 3 cannot exceed , L C V . These expressions allow us to assess the permitted timescales for runoff generation, which must be greater than the breaching timescale, B T , but less than the maximum timescale, max T (where Figure S1). closed-basin lake. After the breach, the open-basin lake volume remains constant and the closed-basin lake has a higher filling rate equivalent to the combined catchments of both lakes (since inflow and outflow in the upstream open-basin lake are now balanced). Without independent constraints, these runoff episode limits can only be constrained relative to the breaching timescale.

Temporal Constraints Using Climate Models
To provide absolute constraints on the runoff episode duration, we derived the lower, B T , and upper, max T , limits by building on the expressions provided above (Equations 1-3; Supplementary Text S4; Figure S1). For a given embedded coupled system, the runoff episode limits are given as To explicitly solve for these durations, we use precipitation rates from existing climate model outputs as a proxy for P during the runoff episode. Out of the existing constraints outlined in Figure 1, we selected precipitation rate outputs from four global climate models based on data availability and their full coverage of the planet (Figures 3a-3d; Guzewich et al., 2021;Kamada et al., 2020;Steakley et al., 2019;Wordsworth et al., 2015). For each coupled lake system, we extracted the P value averaged within the total lake watershed separately from each model output (e.g., Figure 3e). Note that each model considered various different scenarios, resulting in a total of 16 model outputs (Supplementary Table S3). This provides, for each coupled system, a range of durations that are permitted both by its morphology and the model-dependent runoff rate.
STUCKY DE QUAY ET AL.
10.1029/2021GL093523 5 of 10  Table S3 lists all 16 climate scenarios. Black polygons = total watershed of lake systems. Latitudes and longitudes are the same for (a)-(d). (e) Example calculation of regional runoff rate, P, for each system, which is averaged across the combined watershed and lake area. Hillshade from MOLA topography.

Relative Runoff Episode Limits
The geometries of paleolakes allow us to assess a range of timescales for which a These results suggest that, for the most tightly constrained system (i.e., the lowest / max B T T value; black arrow in Figures 2e and 2f), runoff generation can only continue for ∼1.6 2.6  longer than the time it took to initially breach the open-basin lake. For this system (Basin ID 187/9; Table S2), the open-basin lake spends a minimum of ∼40% of the runoff episode duration filling up before it breaches. Systems with larger / max B T T values may not require runoff cessation shortly after breaching, but do not explicitly preclude it. As a result, open lake systems on Mars may spend a large portion of their evolution filling up as closed lakes, as opposed to as stable open lakes.

Distribution of Absolute Runoff Episode Durations
Using Equations 4 and 5, we solved for B T and max T for each of the seven coupled systems using P values extracted from global climate models (Figure 3; Table S3). The resulting distributions for runoff episode durations are shown in Figure 4: each panel illustrates the number of coupled systems with B T and max T values that bound the episode duration, T , specified on the x-axis (i.e., systems where B max T T T   is satisfied) for that climate scenario. As a reference, Figure 4a shows the permitted temporal distributions if we assume globally constant runoff rates from geological observations by Irwin et al. (2005). Here, a 10,000-year runoff episode duration only satisfies one coupled system if the runoff rate was 60 mm/d, but would satisfy three systems if it was 1 mm/d. Each remaining panel (b-f) shows the timescale distributions for four different global climate studies, including different scenarios within each to explore how various parameters affect timescale distributions (Supplementary Table S5). Notably, none of the models satisfy all 7 coupled systems for a single duration bin (see Section 4.1). Figure 4b compares timescales using two end-member precipitation scenarios for rainfall and snowfall (wet, warm at 1.0 bar vs. cold, icy at 0.6 bar, respectively, from Wordsworth et al., 2015). Figure 4c shows the effect of increasing surface pressure (from 0.5 to 2.0 bar; Kamada et al., 2020). Climate scenarios with higher pressures, and consequently greater rainfall, generally result in shorter timescales, except for the highest surface pressure of 2.0 bar, where timescales increase. Using model outputs from Guzewich et al. (2021), we find that lower obliquities result in reduced runoff episode durations; however, global equivalent water inventory size has negligible effects (Figures 4d and 4e). Finally, durations required for an impact-induced atmosphere in Figure 4f show that impactor size does not significantly affect the distributions (50 and 100 km-diameter; Steakley et al., 2019). The distributions in Figure 4 assume a semiarid regime with AI  0.26, however results with no evaporation show less than an order of magnitude difference (Supplementary Figure S3).

Implications for Ancient Hydroclimate Scenarios
For all coupled systems, max T always exceeds B T , consistent with the observation that the closed-basin lakes tend to have larger basin volumes (relative to total contributing area; Figure S2). The mere existence of coupled lake systems implies that runoff generation was sufficiently intense and/or prolonged such that the contained open-basin lake breaches, but not enough to breach the downstream closed-basin lake (Figure 2). Unlike most terrestrial systems, these lakes would be less capable of achieving steady-state because STUCKY DE QUAY ET AL.
10.1029/2021GL093523 6 of 10 the flat-floored, steep-walled crater basins would not allow lake area to continually increase until enhanced evaporation could offset runoff input. These systems thus point to a climate regime that comprises lake filling and overflow, followed-sometimes shortly-by runoff cessation, all within the timeframes indicated by Figure 4.
Overall, the shortest runoff episode durations (1 year) are observed for the 100 km-diameter impact heating scenario from Steakley et al. (2019), whereas the longest durations correspond to Wordsworth et al. (2015)'s snowfall scenario (up to 6 10 yr; see following section for discussion on snowfall vs. snowmelt). In addition to the relative position of each distribution in Figure 4, both the width and height of each individual distribution provide further information. The wider the distribution, the larger the total range of timescales that are permitted by coupled lake hydrology. The taller the distribution (where maximum = 7), the greater the number of lakes that are satisfied by any given runoff episode duration bin. As such, the peak of each distribution denotes the T bin that satisfies the most coupled lake systems. If all systems were formed in runoff episodes of similar durations, then this peak T bin corresponds to the most probable duration for that climate scenario's distribution.
Distributions of T tend to span 3-4 orders of magnitude, suggesting a wide range of permissible timescales for all scenarios. When assessing the distribution peaks, we find that no single T bin can satisfy all 7 coupled systems for any scenario. Most scenarios can satisfy up to 5 systems, with the two exceptions being the 2.0 bar scenario from Kamada et al. (2020), which satisfies 6 systems for T  1000 yr, and the Wordsworth et al. (2015) snowfall scenario which only satisfies 4 at most (Figure 4). In general, most distribution peaks lie between 100 and 10,000 years, suggesting this range of breaching episode duration satisfies the greatest number of coupled system hydrologic constraints across the planet. Importantly, other runoff episodes (either before or after the breaching episode; Figure S1) could have been shorter, but no episode could exceed the maximum durations at any point in early Mars' history, since closed-basin lakes did not overflow.   Table S5); GE = global equivalent; precip. = precipitation; atm. = atmospheric; dashed line in (f) = upper limit of 10 Mars years for impact heating scenario.

Runoff Rates and the Geological Record
Global climate models provide valuable quantitative inputs for assessing runoff rates required to fill our lake systems as a function of space. Although geological estimates are important for understanding reach-scale, channel-forming hydrology, they are more challenging to extend over regional-to global-scales precisely due to being both localized and related to peak hydrologic conditions (Figure 1; Table S1). For example, at first glance, runoff rates modeled by Steakley et al. (2019) most closely match those estimated from geomorphic observations (Figure 1). However, such an impact-heated atmosphere can only persist for a maximum 10 years (Steakley et al., 2019;Turbet et al., 2020), and systems requiring longer time periods are not permissible. Despite this cut-off, the runoff rates provided are sufficient to fill the lakes, and are able to satisfy 5 systems for a runoff episode duration 10 years (Figure 4f). Importantly, though, for discrete events such as an impact, it seems most likely that T was spatially homogeneous, as it would reflect a global heating event.
Aside from the precipitation models used in this study (Figure 4), others have also explored snowfall or ice accumulation rates as a function of space (Figure 1). These snow/ice accumulation distributions are commonly compared to locations of fluvio-lacustrine features such as valley networks (e.g., Wordsworth et al., 2015). However, the relationship between snow accumulation and subsequent runoff rates is not well understood, and variability in processes such as snow ablation could result in spatial discrepancies between snowfall and resulting runoff. Other studies have derived snowmelt production rates (Figure 1), however, when calculated in global models they have yet to generate sufficient runoff in the required mid-latitude regions (e.g., . This could explain why the snowfall distribution in Figure 4b satisfies less systems than all other models for a single time bin, since the duration of liquid water availability (e.g., snowmelt) may not be directly related to snow/ice accumulation at any given location.

Lake Persistence and Episodic Climate Forcing
The sparsity of well-preserved fluvial deposits in older, valley network-fed lake basins results in limited independent geological constraints on lake persistence to test our results. Nonetheless, estimates for total delta building timescales have been previously calculated for two valley network-fedbasins: Eberswalde (a potential closed-basin lake) and Jezero (an open-basin lake). For Eberswalde, total delta-building runoff duration could have lasted 4 6 10 10  yr (Irwin et al., 2015;Moore et al., 2003), with perhaps a lake persisting for > 5 10 yr (Bhattacharya et al., 2005). For Jezero crater, total delta-building duration estimates have similarly ranged from 2 6 10 10  yr (Fassett & Head, 2005;Lapôtre & Ielpi, 2020). Thus, although these two systems are not included in our coupled database, their depositional fluvio-lacustrine temporal constraints are consistent with our runoff episode duration results, despite the wide range of climate scenarios invoked. Future insights into lake level persistence and variability, as captured by delta aggradation and progradation, can hopefully be obtained as the Mars 2020 Perseverance rover explores Jezero crater, providing key independent constraints to test our results.
How do our episode duration estimates fit into the bigger picture of Mars' paleoclimate? Previous studies of valley network evolution suggest that early Mars was characterized by a long-lived runoff-producing climate lasting > 5 10 yr, assuming some intermittency frequency (Barnhart et al., 2009;Hoke et al., 2011).
Our estimated durations suggest runoff production occurred in individual episodes lasting 2 5 10 10  yr, separated by periods of water loss sufficiently long to reset the lake systems. As such, numerous individual runoff episodes likely comprised the total runoff-producing climate, interspersed by periods of negligible runoff. This is broadly consistent with mineralogical records suggesting wet climates were punctuated by hyper-arid intervals (Ehlmann & Edwards, 2014), as well as the presence of deeply incised inlet valleys that fed our paleolakes ( Figure S2), implying multiple cycles of runoff/erosion (e.g., Rosenberg & Head, 2015;Luo et al., 2017). Ultimately, our lake-filling runoff episodes likely occurred during favorable climatic conditions associated with extremes of quasi-cyclical climate changes. Whether this cyclicity was modulated through astronomical variability (e.g., obliquity; Toon et al., 1980), geologically derived fluctuations (e.g., redox oscillations; Wordsworth et al., 2021), or other driving forces, an intermittent climate generating repeated runoff episodes that lasted hundreds to thousands of years would be needed in order to reconcile with the Late Noachian/Early Hesperian hydrological record on Mars.