www.sciencemag.org/content/356/6339/727/suppl/DC1

Supplementary Materials for
Global drainage patterns and the origins of topographic relief on Earth, Mars,
and Titan
Benjamin A. Black,* J. Taylor Perron, Douglas Hemingway, Elizabeth Bailey,Francis Nimmo, Howard Zebker*Corresponding author. Email: bblack@ccny.cuny.edu
Published 19 May 2017, Science 356, 727 (2017)
DOI: 10.1126/science.aag0171

This PDF file includes:
Materials and Methods
Figs. S1 to S8Tables S1 and S2
Captions for Databases S1 and S2References
Other Supplementary Materials for this manuscript includes the following:(available at www.sciencemag.org/content/356/6339/727/suppl/DC1)
Database S1 (Excel)
Database S2 (Zipped archive)

Materials and Methods
Our work relies on three principal datasets: maps of fluvial features on Titan, Earth, and Mars;
spherical harmonic models for the topography of each of these bodies; and results from anumerical model of landscape evolution under a range of uplift conditions. We outline ourmethods for each dataset below.
Mapping fluvial features
Earth. To identify major fluvial features on Earth, we used 1/8th degree grids of global flow
accumulation, basins, and flow length derived from orbital altimetry and corrected manually(17). We identified each sink as a point of maximum flow accumulation within each basin, andthe main trunk as the flow path linking each sink to the most distant point (along flow paths)
within each basin. For rivers that spanned more than ~50 km, we identified sampling points at aninterval of ~50 km, or every 5th grid point (the results are relatively insensitive to this selection;
we chose this interval because it was sufficiently large to avoid sensitivity to the 1/8th degree
resolution of the drainage dataset). We analyzed only rivers longer than this 50 km interval.
Mars. We relied on a global database of Martian river networks (18). Sinks were identified as the
topographically lowest extremities within each network; the validity of these sinks was spotchecked visually. The main trunk of each network was taken to be the path linking the most
distal point within a network to the sink. We identified sampling points along this main trunk at
an interval that matches the interval we used on Earth, but scaled by the ratio
RadiusMars/RadiusEarth (this interval equates to ≈30 km on Mars). We analyzed only networks
longer than this interval.
Titan. We manually selected 71 drainage networks from a global database compiled from theCassini spacecraft’s Synthetic Aperture Radar (SAR) swaths Ta to T71 (12). These drainages
were delineated on the basis of: i) linear geometry cross-cutting other SAR features; ii) light-darkpairings indicative of narrow topographic features; iii) branching geometries; iv) drainage intofeatures interpreted as lakes (11, 12). We selected drainages where we could confidently identify
a sink location on the basis of junction angles, progressive downstream widening of valley
features, orientation relative to features interpreted as lakes, and/or locally available stereo
topography (34). The source and sink coordinates we identified for all 71 networks we analyzed
on Titan are given in Database S1. The main trunk of each network was taken to be the pathlinking the most distant point within a network to the sink. We identified sampling points along
this main trunk at an interval that matches the interval we used on Earth, but scaled by the ratio
RadiusTitan/RadiusEarth (this interval equates to ≈21 km on Titan). We analyzed only networks
longer than this interval. We defined a subset of North Polar networks as all eligible drainages
with sinks located north of 60 °N latitude.
Spherical harmonic topography
We constructed spherical harmonic models of topography for Titan from RADAR-based
topography (21, 31, 35) referenced to the geoid (36); for Earth from the Earth2012
topography/bathymetry model (referenced to sea level) (22); and for Mars from the Mars Orbiter2

Laser Altimeter (MOLA)-derived, aeroid-referenced topography model (37, 38). A robust
spherical harmonic expansion for Titan’s topography is available only up to degree (ℓ) 6 (ref.
(21)). Titan has less long-wavelength relief than Earth or Mars (Fig. 1), indicating that themagnitude of topography alone cannot explain conformity. Our spherical harmonic models ofterrestrial topography include bathymetry (Fig. S6). Oceanic basins form an integral part ofEarth’s plate tectonically derived topography (though the weight of the oceans deepens ocean
basins, and fluvial erosion does not operate on the seafloor). Ancient oceans may also once haveexisted on Mars (3).
Data availability
Coordinates for all drainage networks analyzed on Titan are available in spreadsheet form in
Database S1. The data required to generate spherical harmonic topography for Earth, Titan, andMars are available from the references provided in the text.
Landscape evolution model
Our landscape evolution model considers the effects of fluvial incision, uplift, and hill slope
erosion (parameterized with a critical slope) on surface topography (11, 28, 39). The model isavailable at https://github.com/MITGeomorph/Tadpole.
Variable vs. uniform uplift. The landscape evolution model solves the stream-power equation for
the time evolution of fluvially eroded topography(28, 39):
𝑑𝑑𝑑𝑑
= 𝑈𝑈 − 𝐾𝐾𝐴𝐴𝑚𝑚 |∇𝑧𝑧|𝑛𝑛
𝑑𝑑𝑑𝑑

(S1)

with elevation z, contributing area A, stream-power coefficient K, rate of uplift relative to aboundary U, m a constant taken to be 0.5(40), and n a constant taken to be 1 (41). Hill slope
processes, which occur at finer scales than are of interest for comparison with our global
topographic models, are approximated with a critical slope of 0.6, which prevents slopes frombecoming unrealistically steep at drainage divides (11).
The power-law relationship between incision rate and channel slope and contributing
drainage area in Eq. (1) can be derived from the assumption that channel erosion rates scale with
bed shear stress due to flow in a channel (28, 40), and is motivated by the observed negative
correlation between channel slope and contributing drainage area on Earth (42). In our
simulations, we assume a spatially constant stream-power coefficient, though in practice climaticand lithologic variations can lead to spatial variability in K (41). Plate tectonics produces apatchwork of continental lithologies, and the orogenic feedback loop links mountain-building,exhumation, and climate (43). The first effect might be expected to prolong the effects of plate
tectonics in suppressing topographic conformity even after active deformation has ceased; thesecond effect might be expected to complicate the relationship between rock uplift, topography,
and drainage patterns.
We derive a non-dimensional form of the governing equation following (28):3𝑑𝑑𝑑𝑑′ 𝐾𝐾𝐾𝐾2𝑚𝑚 ′𝑚𝑚=
𝐴𝐴 |∇′𝑧𝑧′| + 1
𝑑𝑑𝑑𝑑′
𝑈𝑈

(S2)

with 𝑧𝑧 = 𝑧𝑧 ′ 𝐿𝐿 , 𝐴𝐴 = 𝐴𝐴′ 𝐿𝐿2 , 𝑡𝑡 = 𝑡𝑡 ′ 𝐿𝐿/𝑈𝑈 , ∇= ∇′ /𝐿𝐿 . The lengthscale L is chosen to represent the
distance from the drainage divide to the sink.
We initialize our simulations with randomly generated, autocorrelated initial topography;
the lowest 70% of the initial surface is set to be a topographic sink as an analog for Earth’s
oceans (Fig. S7). We tabulate model parameter values in Table S1. In the simulations weanalyzed, 13 ± 12% (1σ) of the land area is internally drained. For comparison, 14 ± 10% (1σ;the uncertainty reflects variability across continents) of Earth’s unglaciated land area is internallydrained (44). To approximate the effects of plate tectonics, which localizes crustal thickening,
the variable uplift cases have a pseudo-random pattern of autocorrelated, spatially non-uniform
uplift superposed on the background uniform uplift field.
In place of spherical harmonic decompositions to characterize the model topography, weused 2-D Fourier transforms as described below.
Impact cratering. Impact cratering is one of the key mechanisms for relief generation on Mars.
Crater counting of fluvial landscapes suggests that river activity reached a climax around theNoachian-Hesperian boundary, followed by a decline in activity that roughly coincided withwaning cratering activity after the Late Heavy Bombardment (10, 15). Landscape evolution
modeling and crater counting further suggests that the terminal period of relatively intense
fluvial activity lasted at least 103-104 years, and more likely spanned 105-107 years of episodic
activity (45, 46). Ancient large-scale topography on Mars, including the hemispheric dichotomy,
predates fluvial incision (1, 5). This topography may have been generated through basin-scale
impacts (4) or through other unknown processes. Geomorphic mapping suggests that valley
networks were strongly influenced by ancient topographic gradients, but that younger NoachianHesperian cratering did modify and disrupt some river paths (2, 23).
To investigate and quantify the influence of impact cratering on topographic conformity,
we ran three landscape evolution model ensembles, each with ten simulations (Fig. S8), in which
we considered (i) an initial surface with randomly generated, autocorrelated initial topography(our control ensemble, which represents ancient topography without any younger, fresh craters),(ii) an initial surface with randomly generated, autocorrelated topography, with superposed
impact craters (to represent erosion of ancient topography that has experienced more recent
cratering), (iii) an initial surface with randomly generated, autocorrelated topography, withsuperposed impact craters, and with additional impact cratering occurring in tandem with fluvial
erosion (to represent cratered, ancient topography that undergoes fluvial erosion in tandem withcratering, for example during the Late Heavy Bombardment).
To account for the greater occurrence of smaller craters, we assumed the size-frequency
distribution for Martian impact craters as adapted from the lunar production function for craters
1-16 km in diameter (47). The autocorrelated initial topography represents ancient topography,
for example due to basin-forming impacts. However, we did not generate fresh craters larger
than 16 km diameter, because our goal was to study the interaction of impacts with pre-existing
topography and valley networks (and larger impacts obliterate both for our chosen grid size of 50
km by 50 km). We used scaling relationships for Mars highland craters (48) to calculate crater
depths for the strength and gravity regimes, and we employed polynomial fits from (48) to
calculate the shapes of axisymmetric cavities, rims, and uplift zones. Following (49), we
4

calculated the final topography outside the crater rim as a weighted average of the pre-impact
topography and the impact-generated topography, where the weighting declines linearly from therim to a distance of three radii from the crater center.
Our simulations to investigate impact cratering differ from those designed to investigate
the effects of plate tectonic-style deformation in that for the purposes of accurate impact crater
depth-diameter scaling, we have chosen to make these simulations dimensional. We used model
grids that represented lateral dimensions of 50 km by 50 km, with 125 meter horizontal
resolution. Crater interiors in the model are assumed to be flooded and therefore do not erode.
The complete list of parameters and parameter values used in the impact cratering simulations isgiven in Table S2. The true duration and rate of valley network incision on Mars are unknown(10, 50), and the rate and duration trade off in the model to determine the total amount ofincision. Trunk channels of Martian valley networks incised ~50 to 350 m into older terrains(10), creating drainage densities of ~10-2 km-1 (51). We selected values for K (see Table S2) andsimulation duration (60 My) that generated cumulative erosion that qualitatively matched typicalvalues of drainage density and trunk channel fluvial incision observed on Mars.
Analysis
The metrics we applied to compare drainage orientations with topography are illustrated
graphically in Fig. S1. To avoid the need to weight fluvial features by size, and to integratechanges in flow direction and downslope direction along river courses, we performed both the%d and Λ analyses at intervals along the main trunks of major fluvial features (see Fig. S1).
Synthetic networks. The maximum size of observed drainage networks differs on Titan, Earth,and Mars, which might a priori influence the scale of topography reflected in drainage
orientations. To avoid scale dependence in our comparison, we used a point spacing that is
uniform relative to planetary radius on each body. To test whether our algorithm displays any
bias related to drainage network scale, and to investigate the likelihood of false positive results,
we repeated our analyses on synthetic datasets. These datasets comprised line segments of
uniform length distributed and oriented randomly across the simulated landscape (Fig. S2). We
generated 1,000-10,000 such synthetic segments for each test, and we analyzed the portion of
those synthetic segments that crossed the landmasses in our simulations. In total, we conductedfour tests on representative simulations from the ensemble analyzed in the manuscript: with line
segments spanning either 1/2 or 1/20 the domain of our numerical simulations, and on
landscapes in which our simulations included either uniform or variable uplift. In all cases, ouranalysis yielded statistically indistinguishable results with a topographic conformity of zero (Fig.
S2, panels E-F). In other words, the random synthetic networks showed no preferred orientation
relative to the topography, and our algorithm displayed no measureable bias due to systematicdifferences in drainage size. These tests support the robustness of our results for Titan, Earth, andMars.
Downhill percentage (%d). At each upstream-downstream pair of sampling points along eachdrainage path, we determined whether the upstream point is at a higher elevation according to the
topographic model at a given maximum degree. For each body (and each model run) %d is thepercentage of all upstream-downstream pairs that pass (i.e. the upstream point has a higher
elevation according to the model topography). At infinite resolution (or as wavenumber k→200
5

for the 400×400 model grids), %d should approach 100%, because liquid flows downhill withrespect to the geoid.
Topographic conformity (Λ). At each sampling point along each drainage path, we determined δ,
the angle between the steepest descent direction and the flow direction. For the planetary bodies,steepest descent was determined using MATLAB’s gradientm function for the gradient on
spheroidal bodies; for the model runs, steepest descent was calculated according to the D-infinity
flow routing algorithm(52). We calculated flow direction as the azimuth of the vector linking
each sampling point to the next downstream point. We defined the topographic conformity at agiven maximum spherical harmonic degree or wavenumber as the median value of cosine(δ)
calculated at all sampling points: Λℓ = median(cos(δ)).
Uncertainties. The uncertainties indicated in Figure 1, 2, S2, and S3 represent the 95%
confidence interval for the median. The 95% confidence interval for the median is bounded by
the jth and kth observation in a ranked list of n observations, where (53, 54):
𝑗𝑗 = 𝑛𝑛 × 𝑞𝑞 − 1.96�𝑛𝑛 × 𝑞𝑞(1 − 𝑞𝑞) , and 𝑘𝑘 = 𝑛𝑛 × 𝑞𝑞 + 1.96�𝑛𝑛 × 𝑞𝑞(1 − 𝑞𝑞) .
Here q = 0.5, because by definition the median divides the dataset into two quantiles.
The uncertainties in Figure 3 represent two standard errors of the mean %d (Figure 3c)and Λ (Figure 3d) values across the ensemble of ten simulations with spatially variable uplift andten simulations with spatially uniform uplift.

6Fig. S1. Schematic illustration of two proxy metrics for the agreement between riverorientations and topography at a given scale. The upper panel illustrates the definition of thepercent downhill metric (%d) and the lower panel the definition of the conformity factor (Λ).

7

Fig. S2. Synthetic tests to identify any bias related to the scale of measured drainages. (A)Shaded relief map of model topography produced with spatially variable uplift. (B) As in A, but
for spatially uniform uplift. (C) Short synthetic drainages, superposed on outline of topography
from A. (D) As in C, but with spatially uniform uplift from B. (E) As in C, but with long
synthetic drainages (F) As in E, but for uniform uplift. (G) Topographic conformity Λ forsynthetic dataset with spatially variable uplift. (H) Topographic conformity Λ for synthetic
dataset with spatially uniform uplift. Error bars correspond to the 95% confidence interval for the
median.
8

Fig. S3. Long-wavelength topographic conformity is lower in Titan’s north polar region. As
in Fig. 2, but including Titan’s north polar region (defined here as the region northwards of60°N). Error bars in (B) correspond to the 95% confidence interval for the median. The median
conformity factor values for Titan’s north polar region at degrees 4-6 are lower than those for
Titan as a whole. Given the small sample size, the 95% confidence intervals overlap, but theoffset in median values supports differences in the geologic history of Titan’s north polar region
relative to the rest of Titan.

9Fig. S4. Drainage networks and deformation interact through time. To investigate how
temporal variations in the pattern of deformation influence topographic conformity, we
conducted simulations in which variable uplift (which represents dominantly short-wavelength
deformation associated with plate tectonics) gave way to either (A) uniform uplift, which
represents dominantly long-wavelength deformation, or (B) zero uplift, which represents tectonic
quiescence. Curves represent individual simulations. We find that low Λ is a signature of actively
generated variable-uplift plate tectonics. If variable uplift is followed by uniform uplift, then thesignature of that variable uplift will gradually be erased. If variable uplift is followed by zero
uplift, then Λ increases at first as drainages conform with topography. Ultimately, landscapes
where virtually all relief has been erased can also display poor topographic conformity.

10

Fig. S5. Differences in topographic conformity on Mars versus Earth extend to shorter
wavelengths. Solid lines (and filled circles in right panel) show the percent of points that are
uphill (according to the topography at a given spherical harmonic degree) of the nearest
downstream point. Dashed lines (and empty circles in right panel) show the percent ofheadwaters that are uphill (according to the topography at a given spherical harmonic degree) of
the drainage network outlet. Small panel at right shows the same metrics computed for the
ETOPO2 gridded 2-minute global relief dataset for Earth (55) and four pixel-per-degree MOLA
topography for Mars (37). We expect %d to reach 100% when topography is perfectly resolved,
drainage networks are perfectly delineated, and drainages have not experienced deformation after
the era of fluvial activity. On Earth, the delineation of fluvial features was performed at 1/8
degree resolution from HydroSHEDS and Hydro1k data (17), distinct from the higher resolution
ETOPO2 dataset, resulting in slight misalignments at the scale at which fluvial valleys areresolved. These misalignments do not affect long wavelength topographic conformity, because
fluvial valleys are not resolved at these wavelengths. We attribute the plateau in %d on Mars to
deformation after incision of valley networks (24, 25). For Earth and Mars, calculation of %d forthe headwaters versus outlet of each drainage shows that %d values on Earth converge to 100%across the scale of the drainage, whereas %d values on Mars reach ~95%.

11

Fig. S6. Spherical harmonic models for the topography of Titan, Earth, and Mars. For
Earth and Mars, we show maximum spherical harmonic degrees (ℓ) of 3,6, and 20. At present,
topography is not well constrained for ℓ>6 for Titan (21). Coloring reflects elevation referenced
to the geoid.

12

Fig. S7. Example model topography, filtered to increasing maximum spectral
wavenumbers. We show the model state at the conclusion of variable uplift (left column) anduniform uplift (right column) runs with identical initial conditions. We filter the topography tomaximum spectral wavenumbers of k=1, k=5, and k=200 wavelengths across the domain, as
indicated on each panel.
13

Fig. S8. Shaded relief maps of landscape evolution simulations with cratering. We ran threeensembles (each with ten simulations spanning 60 My) to examine the influence of impact
cratering on topographic conformity. (A) In the control ensemble, fluvial erosion of an initiallyuncratered surface proceeded without interference. (B) We considered fluvial erosion of thesame initial surfaces in the control ensemble, but with impact topography superposed on theinitial surface as described in the Materials and Methods. (C) We also considered fluvial erosionof the same cratered initial surfaces, but with further impacts that occurred during the course of
our simulations, disrupting the topography (note the presence of truncated valley networks). This
simulation is also shown in Fig. 4A. Each snapshot shows the state of the simulation after 60 My.

14

Table S1. Landscape evolution modeling parameters.
Parameter

Value

Notes

Lateral grid dimensions

400 × 400

Enhanced uplift relative tobackground uplift

10

Steam power coefficient K

5 × 10-6 m(1-2m) yr-1

Drainage area exponent m

0.5

ref. (40)

Slope exponent n

1.0

ref. (41)

Critical slope

0.6

Slope of the power
spectrum of initial red noise
topography

2.0

Fraction of initial
topography assigned to be a
fixed base level (torepresent oceans)

0.7

Slope of the power
spectrum of red noise
surface used to define zones
of enhanced uplift invariable uplift simulations

1.3

Only in variable uplift
simulations. Less positive
values translate to morevariance at shorter
wavelengths

Fraction of this surface
assigned to experienceenhanced uplift in variable
uplift simulations

0.35

Only in variable uplift
simulations.

Simulation duration

10 My

Only in variable uplift
simulations

ref. (27)15

Table S2. Landscape evolution modeling parameters for impact cratering simulations.
Parameter

Value

Notes

Lateral grid dimensions

400 × 400

Lateral grid spacing

125 m × 125 m

Steam power coefficient K

1 × 10-8 m(1-2m) yr-1

Drainage area exponent m

0.5

ref. (40)

Slope exponent n

1.0

ref. (41)

Critical slope

0.6

Slope of the power
spectrum of initial red noise
topography

2.0

Fraction of initial
topography assigned to be a
fixed base level (torepresent lakes or oceans)

0.1

Simulation duration

60 My

ref. (27)16Database S1. Source and sink coordinates for analyzed drainage networks on Titan (see Excel
file with tabulated coordinates).
Database S2. Spherical harmonic expansion coefficients representing Titan’s topography withrespect to its geoid.

17

References and Notes
1. R. J. Phillips, M. T. Zuber, S. C. Solomon, M. P. Golombek, B. M. Jakosky, W. B. Banerdt,
D. E. Smith, R. M. Williams, B. M. Hynek, O. Aharonson, S. A. Hauck II, Ancientgeodynamics and global-scale hydrology on Mars. Science 291, 2587–2591 (2001).
doi:10.1126/science.1058701 Medline
2. R. P. Irwin, A. D. Howard, Drainage basin evolution in Noachian Terra Cimmeria, Mars. J.
Geophys. Res. Planets 107, 10-1–10-23 (2002). doi:10.1029/2001JE001818
3. J. T. Perron, J. X. Mitrovica, M. Manga, I. Matsuyama, M. A. Richards, Evidence for anancient martian ocean in the topography of deformed shorelines. Nature 447, 840–843(2007). doi:10.1038/nature05873 Medline
4. J. C. Andrews-Hanna, M. T. Zuber, W. B. Banerdt, The Borealis basin and the origin of themartian crustal dichotomy. Nature 453, 1212–1215 (2008). doi:10.1038/nature07011
Medline5. S. Bouley, D. Baratoux, I. Matsuyama, F. Forget, A. Séjourné, M. Turbet, F. Costard, Late
Tharsis formation and implications for early Mars. Nature 531, 344–347 (2016).
doi:10.1038/nature17171 Medline
6. G. C. Collins, W. B. McKinnon, J. M. Moore, F. Nimmo, R. T. Pappalardo, L. M. Prockter, P.
M. Schenk, “Tectonics of the outer planet satellites,” in Planetary Tectonics, T. R.
Watters, R. A. Schultz, Eds. (Cambridge Univ. Press, 2009), pp. 264–350.
7. F. Nimmo, B. Bills, Shell thickness variations and the long-wavelength topography of Titan.
Icarus 208, 896–904 (2010). doi:10.1016/j.icarus.2010.02.020
8. J. M. Moore, W. B. McKinnon, J. R. Spencer, A. D. Howard, P. M. Schenk, R. A. Beyer, F.
Nimmo, K. N. Singer, O. M. Umurhan, O. L. White, S. A. Stern, K. Ennico, C. B. Olkin,
H. A. Weaver, L. A. Young, R. P. Binzel, M. W. Buie, B. J. Buratti, A. F. Cheng, D. P.
Cruikshank, W. M. Grundy, I. R. Linscott, H. J. Reitsema, D. C. Reuter, M. R.
Showalter, V. J. Bray, C. L. Chavez, C. J. A. Howett, T. R. Lauer, C. M. Lisse, A. H.
Parker, S. B. Porter, S. J. Robbins, K. Runyon, T. Stryk, H. B. Throop, C. C. C. Tsang, A.
J. Verbiscer, A. M. Zangari, A. L. Chaikin, D. E. Wilhelms, New Horizons Science
Team, The geology of Pluto and Charon through the eyes of New Horizons. Science 351,
1284–1293 (2016). doi:10.1126/science.aad7055 Medline
9. S. D. Willett, S. W. McCoy, J. T. Perron, L. Goren, C. Y. Chen, Dynamic reorganization ofriver basins. Science 343, 1248765 (2014). doi:10.1126/science.1248765 Medline
10. A. D. Howard, J. M. Moore, R. P. Irwin, An intense terminal epoch of widespread fluvial
activity on early Mars: 1. Valley network incision and associated deposits. J. Geophys.
Res. Planets 110, E12S14 (2005). doi:10.1029/2005JE002459
11. B. A. Black, J. T. Perron, D. M. Burr, S. A. Drummond, Estimating erosional exhumation on
Titan from drainage network morphology. J. Geophys. Res. Planets 117, E08006 (2012).
doi:10.1029/2012JE004085
12. D. M. Burr, S. A. Drummond, R. Cartwright, B. A. Black, J. T. Perron, Morphology offluvial networks on Titan: Evidence for structural control. Icarus 226, 742–759 (2013).
doi:10.1016/j.icarus.2013.06.016
18

13. C. Hoorn, F. P. Wesselingh, H. ter Steege, M. A. Bermudez, A. Mora, J. Sevink, I.
Sanmartín, A. Sanchez-Meseguer, C. L. Anderson, J. P. Figueiredo, C. Jaramillo, D. Riff,
F. R. Negri, H. Hooghiemstra, J. Lundberg, T. Stadler, T. Särkinen, A. Antonelli,Amazonia through time: Andean uplift, climate change, landscape evolution, andbiodiversity. Science 330, 927–931 (2010). doi:10.1126/science.1194585 Medline
14. H. V. Frey, J. H. Roark, K. M. Shockey, E. L. Frey, S. E. Sakimoto, Ancient lowlands on
Mars. Geophys. Res. Lett. 29, 22-1–22-4 (2002). doi:10.1029/2001GL013832
15. C. I. Fassett, J. W. Head III, The timing of martian valley network activity: Constraints from
buffered crater counting. Icarus 195, 61–89 (2008). doi:10.1016/j.icarus.2007.12.009
16. C. D. Neish, R. D. Lorenz, Titan’s global crater population: A new assessment. Planet. Space
Sci. 60, 26–33 (2012). doi:10.1016/j.pss.2011.02.016
17. H. Wu, J. S. Kimball, H. Li, M. Huang, L. R. Leung, R. F. Adler, A new global river network
database for macroscale hydrologic modeling. Water Resour. Res. 48, W09701 (2012).
doi:10.1029/2012WR012313
18. B. M. Hynek, M. Beach, M. R. Hoke, Updated global map of Martian valley networks andimplications for climate and hydrologic processes. J. Geophys. Res. Planets 115, E09008(2010). doi:10.1029/2009JE003548
19. Materials and methods are available as supplementary materials.
20. M. A. Wieczorek, “Gravity and topography of the terrestrial planets,” in Treatise on
Geophysics, G. Schubert, Ed. (Elsevier, ed. 2, 2015), chap. 10.05, pp. 153–193.
21. H. A. Zebker, B. Stiles, S. Hensley, R. Lorenz, R. L. Kirk, J. Lunine, Size and shape ofSaturn’s moon Titan. Science 324, 921–923 (2009). doi:10.1126/science.1168905
Medline22. C. Hirt, M. Kuhn, W. Featherstone, F. Göttl, Topographic/isostatic evaluation of newgeneration GOCE gravity field models. J. Geophys. Res. Solid Earth 117, B05407(2012). doi:10.1029/2011JB008878
23. R. P. Irwin III, R. A. Craddock, A. D. Howard, H. L. Flemming, Topographic influences ondevelopment of Martian valley networks. J. Geophys. Res. Planets 116, E02005 (2011).
doi:10.1029/2010JE003620
24. W. Luo, T. Stepinski, Orientation of valley networks on Mars: The role of impact cratering.
Geophys. Res. Lett. 39, L24201 (2012). doi:10.1029/2012GL054087
25. A. Lefort, D. M. Burr, F. Nimmo, R. E. Jacobsen, Channel slope reversal near the Martiandichotomy boundary: Testing tectonic hypotheses. Geomorphology 240, 121–136 (2014).
doi:10.1016/j.geomorph.2014.09.028
26. G. Tobie, J. I. Lunine, C. Sotin, Episodic outgassing as the origin of atmospheric methane on
Titan. Nature 440, 61–64 (2006). doi:10.1038/nature04497 Medline
27. D. L. Turcotte, A fractal interpretation of topography and geoid spectra on the Earth, Moon,
Venus, and Mars. J. Geophys. Res. Solid Earth 92, E597–E601 (1987).
doi:10.1029/JB092iB04p0E597

1928. J. T. Perron, W. E. Dietrich, J. W. Kirchner, Controls on the spacing of first-order valleys. J.
Geophys. Res. Earth Surf. 113, F04016 (2008). doi:10.1029/2007JF000977
29. F. Nimmo, D. Stevenson, Influence of early plate tectonics on the thermal evolution andmagnetic field of Mars. J. Geophys. Res. Planets 105, 11969–11979 (2000).
doi:10.1029/1999JE001216
30. V. Sautter, M. J. Toplis, R. C. Wiens, A. Cousin, C. Fabre, O. Gasnault, S. Maurice, O.
Forni, J. Lasue, A. Ollila, J. C. Bridges, N. Mangold, S. Le Mouélic, M. Fisk, P.-Y.
Meslin, P. Beck, P. Pinet, L. Le Deit, W. Rapin, E. M. Stolper, H. Newsom, D. Dyar, N.
Lanza, D. Vaniman, S. Clegg, J. J. Wray, In situ evidence for continental crust on earlyMars. Nat. Geosci. 8, 605–609 (2015). doi:10.1038/ngeo2474
31. D. Hemingway, F. Nimmo, H. Zebker, L. Iess, A rigid and weathered ice shell on Titan.
Nature 500, 550–552 (2013). doi:10.1038/nature12400 Medline
32. T. Schneider, S. D. B. Graves, E. L. Schaller, M. E. Brown, Polar methane accumulation andrainstorms on Titan from simulations of the methane cycle. Nature 481, 58–61 (2012).
doi:10.1038/nature10666 Medline
33. R. D. Lorenz, E. P. Turtle, B. Stiles, A. Le Gall, A. Hayes, O. Aharonson, C. A. Wood, E.
Stofan, R. Kirk, Hypsometry of Titan. Icarus 211, 699–706 (2011).
doi:10.1016/j.icarus.2010.10.002
34. B. W. Stiles, S. Hensley, Y. Gim, D. M. Bates, R. L. Kirk, A. Hayes, J. Radebaugh, R. D.
Lorenz, K. L. Mitchell, P. S. Callahan, H. Zebker, W. T. K. Johnson, S. D. Wall, J. I.
Lunine, C. A. Wood, M. Janssen, F. Pelletier, R. D. West, C. Veeramacheneni,
Determining Titan surface topography from Cassini SAR data. Icarus 202, 584–598(2009). doi:10.1016/j.icarus.2009.03.032
35. H.A. Zebker, L. Iess, S. D. Wall, R. D. Lorenz, J. I. Lunine, B. W. Stiles, “Titan’s Figure
Fatter, Flatter Than Its Gravity Field,” American Geophysical Union Fall Meeting,abstract P23F-01 (2012).
36. L. Iess, R. A. Jacobson, M. Ducci, D. J. Stevenson, J. I. Lunine, J. W. Armstrong, S. W.
Asmar, P. Racioppa, N. J. Rappaport, P. Tortora, The tides of Titan. Science 337, 457–
459 (2012). doi:10.1126/science.1219631 Medline
37. D. E. Smith, M. T. Zuber, S. C. Solomon, R. J. Phillips, J. W. Head, J. B. Garvin, W. B.
Banerdt, D. O. Muhleman, G. H. Pettengill, G. A. Neumann, F. G. Lemoine, J. B.
Abshire, O. Aharonson, C. D. Brown, S. A. Hauck, A. B. Ivanov, P. J. McGovern, H. J.
Zwally, T. C. Duxbury, The global topography of Mars and implications for surface
evolution. Science 284, 1495–1503 (1999). doi:10.1126/science.284.5419.1495 Medline
38. D. Smith, G. Neumann, R. Arvidson, E. Guinness, S. Slavney, Mars Global Surveyor laser
altimeter mission experiment gridded data record (NASA Planetary Data System, 2003).
39. J. T. Perron, J. W. Kirchner, W. E. Dietrich, Formation of evenly spaced ridges and valleys.
Nature 460, 502–505 (2009). doi:10.1038/nature08174
40. A. D. Howard, G. Kerby, Channel changes in badlands. Geol. Soc. Am. Bull. 94, 739–752(1983). doi:10.1130/0016-7606(1983)94<739:CCIB>2.0.CO;220

41. K. L. Ferrier, K. L. Huppert, J. T. Perron, Climatic control of bedrock river incision. Nature
496, 206–209 (2013). doi:10.1038/nature11982 Medline
42. K. X. Whipple, G. E. Tucker, Dynamics of the stream-power river incision model:
Implications for height limits of mountain ranges, landscape response timescales, andresearch needs. J. Geophys. Res. Solid Earth 104, 17661–17674 (1999).
doi:10.1029/1999JB900120
43. S. D. Willett, Orogeny and orography: The effects of erosion on the structure of mountainbelts. J. Geophys. Res. Solid Earth 104, 28957–28981 (1999).
doi:10.1029/1999JB900248
44. C. Vörösmarty, B. Fekete, M. Meybeck, R. Lammers, Global system of rivers: Its role inorganizing continental land mass and defining land‐to‐ocean linkages. Global
Biogeochem. Cycles 14, 599–621 (2000). doi:10.1029/1999GB900092
45. C. J. Barnhart, A. D. Howard, J. M. Moore, Long‐term precipitation and late‐stage valley
network formation: Landform simulations of Parana Basin, Mars. J. Geophys. Res.
Planets 114, E01003 (2009). doi:10.1029/2008JE003122
46. E. S. Kite, A. Lucas, C. I. Fassett, Pacing early Mars river activity: Embedded craters in the
Aeolis Dorsa region imply river activity spanned ≳(1–20) Myr. Icarus 225, 850–855(2013). doi:10.1016/j.icarus.2013.03.029
47. B. A. Ivanov, Mars/Moon cratering rate ratio estimates. Space Sci. Rev. 96, 87–104 (2001).
doi:10.1023/A:1011941121102
48. S. T. Stewart, G. J. Valiant, Martian subsurface properties and crater formation processes
inferred from fresh impact crater geometries. Meteorit. Planet. Sci. 41, 1509–1537(2006). doi:10.1111/j.1945-5100.2006.tb00433.x
49. N. K. Forsberg‐Taylor, A. D. Howard, R. A. Craddock, Crater degradation in the Martian
highlands: Morphometric analysis of the Sinus Sabaeus region and simulation modelingsuggest fluvial processes. J. Geophys. Res. Planets 109, E05002 (2004).
doi:10.1029/2004JE002242
50. N. Mangold, E. S. Kite, M. G. Kleinhans, H. Newsom, V. Ansan, E. Hauber, E. Kraal, C.
Quantin, K. Tanaka, The origin and timing of fluvial activity at Eberswalde crater, Mars.
Icarus 220, 530–551 (2012). doi:10.1016/j.icarus.2012.05.026
51. B. M. Hynek, R. J. Phillips, New data reveal mature, integrated drainage systems on Mars
indicative of past precipitation. Geology 31, 757–760 (2003). doi:10.1130/G19607.1
52. D. G. Tarboton, A new method for the determination of flow directions and upslope areas ingrid digital elevation models. Water Resour. Res. 33, 309–319 (1997).
doi:10.1029/96WR03137
53. W. J. Conover, Practical Nonparametric Statistics (Wiley, ed. 2, 1980).
54. M. Longnecker, R. Ott, An Introduction to Statistical Methods and Data Analysis (Nelson
Education, Toronto, ed. 5, 2001).
55. W. H. Smith, D. T. Sandwell, Global sea floor topography from satellite altimetry and ship
depth soundings. Science 277, 1956–1962 (1997). doi:10.1126/science.277.5334.1956
21