[8] We first present a geometrical statistic for quantifying the shape of drainage networks, and demonstrate its usefulness for evaluating the extent of erosional modification of a surface. We calibrate our statistic with a simple landscape evolution model, and we validate it through comparisons with terrestrial drainage networks with constrained exhumation histories. Finally, we map fluvial networks on Titan and apply our method to those networks. Our results suggest that some of Titan's most prominent drainage networks have caused spatially averaged erosion of only a small fraction (<∼10%) of the initial topographic relief in the regions where they have formed.

[7] In this study, we analyze the planform geometries of Titan's drainage networks in order to quantify erosional exhumation, and thereby describe the evolution of the landscape through time. Metrics of drainage network shape—such as drainage density, the width function, and channel sinuosity—are commonly applied on Earth to understand the way network structure relates to environmental variables such as precipitation or bedrock lithology [e.g., Rodriguez ‐ Iturbe and Rinaldo , 1997 ]. Our measurements of purely planform geometry are unusual by terrestrial standards, however, because analyses of fluvial incision on Earth typically also rely on elevation data to measure landscape characteristics such as slope, drainage area, and valley relief. In the case of Titan the limited available topographic data currently inhibit such measurements. Even on Earth, it is challenging to assess the cumulative erosion that a landscape has experienced without independent knowledge of the initial conditions. The goal of this paper is to develop and apply a method to assess the extent of erosional development of a fluvial landscape, using only features that can be measured in plan view.

[2] Prior to the arrival of the Cassini‐Huygens spacecraft in orbit around Saturn in 2004, the surface of Titan was the subject of wide‐ranging interest and speculation. Gerard Kuiper was the first scientist to detect methane in Titan's atmosphere [ Kuiper , 1944 ]. Methane is near its triple‐point (91 K and 0.2 bars) at surface conditions on Titan of ∼94 K and 1.47 bars [ Fulchignoni et al. , 2005 ]. This observation led Sagan and Dermott [1982] to envision methane oceans on Titan hundreds of meters deep. Lunine et al. [1983] suggested that methane photolysis might instead result in an ocean dominated by ethane, of even greater depth. Other researchers, however, concluded that global methane oceans seemed unlikely [ Eshleman et al. , 1983 ].

[14] In summary, ΔW is a measure of how much the distribution of segments within a fluvial network deviates from a uniform distribution. In order to test our hypothesis that ΔW should increase as drainage networks dissect a surface, and to quantify any such relationship between ΔW and cumulative erosion, we turn to a numerical model of landscape evolution.

[13] We predict that the cumulative width functions of less mature drainages (such as Figure 3a ) will reflect a uniform distribution of channel segments, because the drainage will not have developed an equant shape with dendritic tributaries. Basins like those in Figure 3b–3c , which are further along in their development, will have sigmoidal or concave cumulative width functions. In order to compare the cumulative width functions of various networks, we calculate the absolute value of the area between each network's cumulative width function and the straight line that represents a uniform width function ( Figures 3g–3i ), and we define this dimensionless quantity as ΔW. This statistic, which measures the deviation of the measured width function from a uniform width function, is conceptually similar to the Cramér‐von Mises criterion, which measures the difference between a theoretical cumulative distribution function and a measured distribution. The ΔW statistic is distinct, however, from the Kolmogorov‐Smirnov statistic, which only considers the maximum difference between a measured cumulative distribution and a theoretical one.

[12] Given these general expectations for how a drainage network should evolve over time, especially the development of dendritic tributaries and the trend toward more equant and integrated basins, we can define a morphologic statistic to quantify the relative maturity of a network. Rodriguez ‐ Iturbe and Rinaldo [1997] define the width function as the number of links in a drainage network at a given distance from the drainage basin outlet. The cumulative width function, as demonstrated in Figures 3g–3i , is the integral of the width function: the fraction of the network that lies within a given flow distance from the outlet.

[11] As the geometric sequence described above proceeds, the drainage networks incise through the initial relief (the initial range of elevations in a landscape over a given horizontal area), gradually overprinting the initial topography ( Figures 3d–3f ). We can define the relative “age” of a fluvial landscape as the extent to which it has modified an initial surface. All else being equal, this sequence will proceed more slowly at larger spatial scales with higher relief. By normalizing the total vertical erosion E (spatially averaged over the landscape) by the initial relief H 0 , we obtain a dimensionless measure of the relative erosional “age” of the landscape. Because landscapes can develop at different rates, the morphologic age does not correspond to an absolute calendar age. Henceforth, we will use the qualitative terms “immature” and “mature” to refer to networks that are in the early stages of fluvial dissection (small E / H 0 ) and networks that have caused extensive erosion of an initial surface ( E / H 0 approaching unity, or exceeding unity if uplift or base level lowering generate additional relief), respectively.

[10] By devising a method to quantify these trends, we can use planform drainage network morphology to gauge the extent of erosional modification. Although this assessment sounds like a simple task, it has not yet been completed for Earth—in part because elevation data and various lines of geological evidence can often be used to constrain the development of terrestrial networks.

(a–c) Schematic drawings illustrating the expected evolution of drainage networks in a landscape with sinks at the top and bottom of the frame. Initially, headward advance dominates. Basins develop roughly uniform width (Figure 3a). Once upslope area is exhausted, lateral tributary growth causes widening (Figure 3b). Interaction among basins concentrates area at intermediate distances from outlet (Figure 3c). (d–f) Three frames from a landscape evolution model run illustrating the evolution predicted in Figures 3a–3c. Time increases by a multiple of five in each successive frame. Colors denote elevation, with blue to red showing lowest to highest elevations, respectively. Maximum initial elevations are ∼400 m, yielding average slopes of ∼0.008. (g–i) Idealized illustration of the relationship between the cumulative width function (red line) and the ΔW statistic (shaded area) for the drainage networks shown in Figures 3a–3c. As shown in the figure, ΔW is proportional to the deviation of the cumulative width function from a uniform distribution of network links. Figures 3g–3i are schematic; actual width functions may not be symmetric about their midpoints.

[9] With few impact craters and little ground‐based or historical data available, Titan's drainage networks provide one of the best sources of information about past erosion and resurfacing patterns. Over time, we expect the geometry of these networks to evolve. Topographic data that resolve the cross‐sectional and longitudinal forms of fluvial valleys on Titan are currently limited. In the absence of widespread topographic data, as discussed above, we must focus our analysis on the planform geometry of drainage networks. Our conceptual model for drainage network development is similar to that outlined by Horton [1945] . The basic stages in this model are shown schematically in Figure 3 , starting with initial perturbations in a random surface ( Figure 3a ). These initial perturbations develop into channels that compete for water and therefore drainage area. Because the erosion rate in channels is proportional to flow discharge (see section 3 ), incipient channels that have smaller drainage areas and therefore smaller discharge than their neighbors will grow less rapidly than their faster‐eroding neighbors. The more rapidly eroding networks will gradually integrate more drainage area and develop simple tributary channels ( Figure 3b ). At first, the channel networks will expand mainly through headward growth of the trunk channel, but once they reach the drainage divide, tributaries will grow laterally and capture additional drainage area. As the surviving networks approach a steady form, they develop dendritic drainage networks with evenly spaced trunk channels and less‐elongated drainage basins ( Figure 3c ).

3. Landscape Evolution Model

[15] Numerical modeling provides a clear theoretical framework for both testing and calibrating our predictions for how the ΔW statistic should behave as drainage networks evolve. We identify several processes that are likely to be most important in shaping fluvial landscapes on Titan and include them in our evolution model.

[16] Drainage basins on Titan, like those on Earth, have developed through a combination of fluvial incision, which creates networks of incised channels, and hillslope processes that govern the response of the surrounding topography to channel incision. The limited stereo topography currently available for Titan confirms that many of the networks visible in SAR are incised valleys with steep valley walls, and therefore bedrock erosion must have occurred during their formation [Tomasko et al., 2005; Soderblom et al., 2007]. Images from the Huygens probe [Tomasko et al., 2005] show branched networks dissecting the landscape down to a certain scale, but with channels that appear to truncate at a scale coarser than the image resolution. This observation implies that at the finest scales, non‐fluvial erosion mechanisms are more important than fluvial incision. However, because we are interested in the shapes of networks which span 20 to 500 km, the relative contribution of hillslope processes to network geometry is likely to be small. Cross sections through valleys near the Huygens landing site show nearly straight hillslope profiles with slopes of roughly 30 degrees [Tomasko et al., 2005], suggesting that hillslope erosion may be modeled with a critical angle, where any hillslopes that steepen beyond this angle experience rapid mass wasting.

[18] We modeled the evolution of the landscape by solving equation (1)forward in time on a regular grid using a finite difference method. We began with random, self‐affine initial surfaces defined by the variance in elevation and the slope of the surface's power spectrum. We assume there are no structures that would introduce strongly oriented initial topography. The lowest 10% of elevations in the initial surface were set to a fixed elevation of zero to simulate a base level, as might be provided for high‐latitude networks by Titan's polar lakes. Initial relief for our model landscapes was approximately 400 m. The average relief of equatorial mountains reported byRadebaugh et al. [2007] was ∼400 m over horizontal profiles tens of kilometers long. The grid size is 400 by 400 points, with Δx = Δy = 125 m. We selected K = 0.005 yr−1 with m = 0.5, and ran the model for 25 Myr. Because the magnitude of K is not known for Titan, these model results are not tuned to be specifically applicable to the rate of erosion on Titan. Rather, they form a test set for our general predictions of landscape evolution and the ΔW statistic. There is a tradeoff between initial relief, K, horizontal scale, and the time interval over which the landscape erodes. We are primarily interested in how drainage network form changes as a function of how much of the initial relief has eroded. Because this evolution is independent of the absolute rate of erosion, our results are valid without reference to the value of K or elapsed time.

[19] Our model runs employed a drainage area exponent of m = 0.5 and a slope exponent of one, which implies that erosion rate is linearly proportional to stream power [Whipple and Tucker, 1999]. Varying these exponents will shift the concavity and response time of the drainage profiles [Whipple and Tucker, 1999], but exploratory model runs not presented here indicate that this has only a minor effect on ΔW. The 50 km by 50 km scale of the model runs is similar to that of most of the networks we analyzed on Titan, though some of the networks are several times larger. In Figure 4, the imprint of initial topography is clearly visible even after 25 Myr of landscape evolution. The self‐affine random surface that we use to initiate the models is qualitatively consistent with available data from Titan, which generally show gentle slopes over long baselines [Zebker et al., 2009], though higher resolution data from the Huygens landing site contain valley side slopes up to 30 degrees over short distances [Tomasko et al., 2005].

Figure 4 Open in figure viewer PowerPoint Example 25 Myr model run, with one drainage network highlighted in red. We only analyze the portion of the network that has eroded through >20% of the initial relief. The ΔW increases over time, eventually approaching a steady value for a given drainage network, even though the topography does not reach a steady state.

[20] Figures 3d–3f show three stages in a typical run of the landscape evolution model. To investigate our hypothesized relationship between fluvial network morphology and the extent of erosion, we require a criterion for defining fluvial networks in our model topography that is comparable to the visibility of fluvial networks in coarsely resolved Cassini SAR images. As we discuss in section 5, deeper valleys should be much easier to detect in SAR imagery of Titan, and very shallow valleys may be invisible. Therefore, in our model runs, we calculated the ΔW only for those portions of the network that have incised a depth greater than 20% of the initial relief of the landscape. 10–25% thresholds all yield similar trends, but lower thresholds correspond with lower confidence in our upper limits on erosion. We discuss this sensitivity in greater depth in section 6.2. In order to avoid fragmented networks, we included all portions of a network downstream from a segment that has eroded to this depth. This approach is comparable to our procedure for mapping discontinuous but clearly related fluvial features.

[21] Four timesteps from a 25 Myr model run are shown in Figures 4a–4d. In Figures 4e–4h, we isolate a single network and show the portion that we analyzed after applying the “visibility” criterion described above. These frames show how, for the most deeply eroded fluvial features, the evolution of the network generally matches our predictions as shown in Figure 3. Because stream incision rate is proportional to drainage area, A, in equation (1), there is a positive feedback between stream incision and drainage area. The more drainage area a given network claims, the faster it incises, which allows it to lengthen or widen, thereby capturing more drainage area. Incision rate is also proportional to land‐surface slope, so networks grow more quickly in the upslope direction. As networks begin to form on an initially undissected surface, they will first grow headward, because headward growth is unimpeded by competition with neighboring networks. This growth leads them to become elongated, as inFigures 4a, 4b, 4e, and 4f. Once networks begin to impinge on one another on opposite sides of a drainage divide, however, they widen and become more equant, as in Figure 4c, 4d, 4g, and 4h.

[22] These trends are reflected in ΔW, as shown in Figures 4i–4l. On the x axis, instead of time, we express the progressive development of the drainages in terms of E/H 0 , the normalized measure of average erosion defined in section 2.1. Very early in the development of the networks, the ΔW is low. It rises quickly to a nearly steady value that varies from network to network. The precise final configuration of each drainage network depends on initial conditions [Perron and Fagherazzi, 2012]; the overall shape of the basin will control the shape of the cumulative width function (of which Figures 4i–4l show examples), and will therefore also influence the ultimate value of ΔW achieved by each network. In order to better estimate the distribution of possible outcomes, we employed a Monte Carlo procedure in which we performed twenty model runs with twenty different initial surfaces that shared similar statistical characteristics. The largest networks (drainage area >10% of the total grid) in these runs were used to construct a calibration data set for ΔW as a function of E/H 0 . We binned these data according to E/H 0 . The resulting distributions of ΔW for various values of E/H 0 are plotted in Figure 5, with the mean trend with increasing E/H 0 shown in Figure 6. Because the distributions for low and high ΔW statistics are associated with different E/H 0 , we can use our results to calculate the statistical likelihood of a given erosional exhumation for drainage networks on Titan or Earth.

Figure 5 Open in figure viewer PowerPoint Six snapshots of the distribution of ΔW from an ensemble of 20 model runs. These histograms show that mean ΔW increases with increasing erosion relative to initial relief, eventually approaching a steady value of ΔW.