Isostatic Control of Axial Rivers and Large Drainage Basins on Passive Margins

More than half of the world's large rivers flow towards the ocean crossing passive continental margins. Here using an analytical solution and numerical models, we demonstrate that on passive margins, river basins may be integrated by major margin‐parallel channels, which form as a flexural isostatic response of the lithosphere to mechanical/erosional unloading along the margin. We analyzed the downstream courses of large rivers flowing across the passive margins and find that the majority of them (31 of 36) have major margin‐parallel channels. Occurrences of these channels are generally consistent with the model predictions, although the exact locations and geometry of these rivers may also be controlled/changed by other factors. Our results suggest that the lithosphere strength has an important control on the geometry of large river systems on passive margins, linking the evolution and routing of the Earth's freshwater systems to its deep interior dynamics.


Introduction
The geometry of river networks determines how water runoff is routed across continents and therefore has a direct influence on freshwater input into the ocean, carbon cycling, and the pattern of freshwater ecosystems. In steep mountainous areas, stream orientation is mostly controlled by local-scale topographic relief and active structures (e.g., Clark et al., 2004). However, many continents are dominated by low-relief interiors, where river networks often do not follow the topographic gradient (Black et al., 2017). It has been argued that low-amplitude topography caused by deep mantle processes may control the geometry of some drainage systems at the Earth's surface (e.g., Cox, 1989), but no universal pattern has been observed to support this hypothesis.
In low-relief cratonic areas, rivers connect to the ocean across passive margins. Drainage networks along passive margins are typically characterized by a dichotomy of river basin sizes, with small coastal rivers that drain directly to the ocean separated from large river basins that drain the continental interior surrounded by a rim of high topography. This topographic bulge is thought to form during continental breakup as a flexural isostatic response to continental stretching (Braun & Beaumont, 1989) and/or progressive erosion of an escarpment or topographic step (Gilchrist & Summerfield, 1990) that separates the coastal plain from the relatively flat, elevated continental interior (e.g., Figures 1a and 1b).
Along most passive margins, this escarpment or bulge is a drainage divide. Evolution of this drainage divide is determined by the competition between denudation and isostatic uplift (Braun, 2018;  is the angle between the downstream flow direction and the direction towards the margin. d a is the distance to the margin, that is the shortest distance from the river pixel to the transitional crust. Green line marks the location of the transitional crust with the thickness of 20 km. Box indicates the location of the topographic profile in b. (b) Topography along a profile within a swath. Black line indicates the median topography. (c) Sketch illustrating the flexural responses to erosion along a passive margin and the theoretical position of the topographic depression where an axial river may form. Slingerland, 1994). On the one hand, river incision and hillslope erosion reduce the elevation of the escarpment/bulge. On the other hand, the lowering of the surface elevation is partially compensated by flexural isostatic rebound triggered by erosional unloading (Figure 1c). These processes act to maintain the escarpment/bulge as a drainage divide and preserve it for tens of millions of years or longer after rifting, even if it has migrated several tens to hundreds of kilometers landward from its original position (Gilchrist & Summerfield, 1990).
Flexural isostasy which works to create and maintain the escarpment also results in low-amplitude subsidence of the continental area behind it ( Figure 1c) (Sacek et al., 2012). We postulate that this topographic low leads to the formation of a valley that parallels the passive margin and truncates the inland river network, diverting flow to a large, margin-parallel river. This "axial river" collects water and sediment from the truncated tributaries that drain expansive portions of the continental interior. Unable to breach the topographic barrier of the escarpment, it flows parallel to the coast until it reaches a coastal river that has incised through the topographic bulge to capture its drainage ( Figure 2). By combining multiple river flows into a large one, the occurrence of an axial river significantly increases the size of the continental interior drainage basins. Therefore, we propose that this flexure-controlled continent-scale drainage network reorganization is responsible for the first-order geometry of major river basins on large, tectonically stable continents, such that the geometry of most rivers draining continental interiors can be predicted from the strength of the underlying lithosphere which ultimately determines the scale and width of its flexural isostatic response.
In this paper, we test this hypothesis by comparing the locations of the axial rivers (1) predicted by an analytical solution of the flexure of the lithosphere (Section 2), (2) simulated by numerical landscape evolution models (Section 3), and (3) observed in the largest river basins on Earth (Section 4). The results support our hypothesis and indicate that the strength of the lithosphere is a fundamental component to the drainage network geometry on passive margins, although the detailed patterns are also controlled by other parameters.

Topographic Depression Due to Lithospheric Flexure
On a continental margin, erosion of a narrow column of the lithosphere can be considered as a line load: in which l is the lithospheric density, g is the gravitational acceleration, h e is the eroded thickness, and Δx is the width of the column. At the incipient stage of continental rifting, we assume that the erosion of the margin was concentrated at a position x = 0, where the initial escarpment developed. Assuming that the continental lithosphere acts as thin elastic plate of effective elastic plate thickness T e , vertical surface displacement of the lithosphere caused by erosion is given by (Turcotte & Schubert, 1982) u in which is the flexural wavelength given by where a is asthenospheric rock density and D the flexural rigidity of the lithosphere defined by where E and are elastic parameters, namely Young's modulus and Poisson's ratio.
As we consider the load is concentrated at x = 0, which makes Δx → 0, it is reasonable to assume that Δx ≪ . In this case, the solution expressed in Equation 2 has the shape as shown in Figure 1c. It not only explains the uplift of the coastal region and the formation of the bulge/escarpment observed along many passive margins but also predicts the formation of a depression at a distance of the region of maximum erosion. As suggested by Sacek et al. (2012), a small amplitude, secondary bulge may also form, which will further contribute to the formation of an axial river on the landward side of and running parallel to the escarpment.
Assuming that the initial escarpment was developed at a tectonic hinge line that parallels the margin, this line should mark a significant transition in the crustal thickness and topography (e.g., Karner & Driscoll, 1999). The position of the axial river should be located at a distance from the crustal transition line. This simple model predicts a power-law relationship between the locations of axial rivers along passive margins and the effective elastic thickness of the underlying lithosphere: Therefore, in the case where isostatic adjustment from sediment loading can be neglected, the position of axial rivers is solely dependent on the flexural strength of the lithosphere, that is its distance to margin, d a , is proportional to the effective elastic thickness raised to a power, T 4∕3 e . Note that for simplification in the following text, we term the location of the hinge line of transitional crust as the margin of the continental plate.

Landscape Evolution Modeling
The above prediction is dependent on important assumptions that the erosional unloading at the initial escarpment is significantly more than anywhere else on the margin and sediment deposition in the fluvial systems can be neglected. To test these assumptions and improve the predictions, we use a landscape evolution model to predict how the landscape responds to the distributed erosion and sedimentation along a passive margin following continental rifting and breakup. We use the approach developed by Braun and Willett (2013) and Yuan et al. (2019) to solve the basic equations governing the rate of fluvial incision, sediment deposition, and hillslope processes (Ahnert, 1977;Davy & Lague, 2009;Whipple & Tucker, 1999): where h is topographic elevation, t is time, A is upstream drainage area, S is local topographic slope, G is a dimensionless number controlling the balance between erosion and sediment deposition in rivers (see Yuan et al., 2019, for details), ∇ 2 h is surface curvature, and K f and K d are fluvial erosion and diffusion coefficients, respectively. The values of these rate coefficients are variable as they are likely dependent on lithology, rock erodibility, and/or climate. In this equation, U is the uplift rate due to erosional unloading, which we obtain by solving the biharmonic equation governing the flexure of a thin elastic plate (Turcotte & Schubert, 1982): We solved these coupled equations in a 1,200 × 1,200 km 2 domain at 5-km resolution and started each model with a 500-to 700-m-high flat plateau to which we added an infinitesimal tilt to force the development of an initial drainage network flowing towards the left-hand side of the model representing a rifted/passive margin. We ran the models for 50 Myr using 10-kyr long time steps and different randomly selected values of K f (10 −7 − 10 −5 m 0.2 /yr), K d (0-0.1 m 2 /yr), and T e (20-120 km). We tested landscapes dominated by detachment-limited rivers and that with sediment deposition by prescribing G in Equation 7 at 0 and 1, respectively; m was fixed at 0.4. For the flexure model (Equations 4 and 8), E, , l , and a were prescribed at 10 11 Pa, 0.25, 2,800 km/m 3 and 3,200 kg/m 3 , respectively. We performed >50,000 model experiments for the settings without (G = 0) and with (G = 1) sediment deposition, respectively. For each model, the shape of river networks is recorded at every 10-Myr interval (e.g., Figure 2).
In order to determine the position of the axial channels, we estimated the river flow direction relative to the orientation of the margin. For each pixel of the river channel, we introduce a measure, where is the angle between the observed downstream direction and the direction normal to the margin ( Figure 1a); the former is calculated over a critical downstream length of 100 km. A value of 1 indicates that the river channel is perpendicular to the margin, whereas a value of 0 corresponds to channels flowing parallel to the margin. Assuming that the river network branching is purely topological and not controlled Gray lines depict plate boundaries (Bird, 2003). Green lines depict locations of transitional crust (thickness of 20 km) with no significant tectonic deformation. Estimates of the effective elastic thickness (T e ) of the continental lithosphere are from Audet and Bürgmann (2011). by local topography, the bifurcation angle between any two branches should be around 72 • (Devauchelle et al., 2012), which should lead to the formation of branches with values around 0.8. We therefore used a low value (<0.3) to define the occurrence of axial rivers on the modeled landscape (e.g., Figure 3).
Our modeling results predict that axial rivers form on >90% of the modeled landscapes, which are located in a wide range of locations relative to the initial position of the escarpment. By summarizing the distances from the margin to the first occurrences of the axial rivers and plotting these distances versus the model T e , we find that the highest densities of points representing the model outputs are distributed along the theoretical line predicted by the 1-D analytical solution ( Figure 5), albeit that the model results with continental deposition are more dispersed.

Axial Channels of the Largest Rivers
In this section, we further test the control of lithosphere flexure on the formation of axial rivers by examining the downstream courses of the 56 rivers with the largest drainage areas that are located between 40 • N and 40 • S (Figure 4). High-latitude drainage systems are excluded from our analysis, as we wish to focus on basins dominated by fluvial processes. Of these 56 rivers 36 flow towards passive margins, where we consider that the perturbation of topography by tectonic activities is minimal. As described above, we use the value of a pixel on the river channel to determine its flow direction relative to the orientation of the margin. To perform the analysis, we used an upscaled (1/8 • resolution) river network model (Wu et al., 2012) and calculated using critical downstream lengths of 50, 100, and 150 km. For most rivers, this choice of the critical segment length had little influence on our estimate of the location of low-segments (see Figure S1 for full details), so here we present only the results using a critical segment length of 100 km. We restricted our analysis to river segments with upstream drainage areas over 10,000 km 2 and within 1,000-km distance from the continental margin. To define the trace of the initial position of the escarpment, that is the transitional zone between continental and oceanic crusts, we used a simplification and considered it at the location where the crustal thickness is 20 ± 5 km (according to the model of Laske et al., 2013).
The results show that for all but one (Webi Jubba in Somalia) of the 36 rivers, high values occur near their outlets, indicating that these channels flow in a direction perpendicular to the margin. Upstream and further inland, rivers start to develop more dendritic patterns. The majority (31 out of 36) of the rivers we analyzed start to diverge or branch into channels with low values (<0.3), implying that the developments of their networks were not purely topological. To determine the positions of axial rivers on a margin, we binned the pixels of a river according to their distances to the margin and estimated a location where at least 40% channel pixels have low-values (<0.3). The threshold (40%) was determined by testing a range of thresholds to capture the average flow direction of the river channels, but for most rivers, using a threshold between 30% and 50% leads to no difference for determining the position of the axial river. The full results are shown in Figure S1, and summarized in Figure 5. To estimate the mechanical strength of the lithosphere, we calculated the mean T e of the lithosphere within 500 km from the margin in a river's drainage area, using the results from spectral coherence analysis between topography and Bouguer gravity anomalies (Audet & Bürgmann, 2011).  (Audet & Bürgmann, 2011). The location of a margin is assumed at the position where the crust thickness is 20 ± 5 km; the crust thickness is according to Laske et al. (2013). Data points are color coded according to the sediment volume within a 100-km-wide swath from the river outlet to the margin, which is calculated using the model of Straume et al. (2019). Shading shows density contours of the individual model predictions from landscape evolution modeling.
Along the Atlantic margin of Africa, the determined distances between the margin and the axial rivers increase with the rigidity of the lithosphere (Figure 5a). The general consistency between theoretical prediction, modeling results, and the data of natural rivers in this region supports our hypothesis that flexural isostasy in response to erosion of the margin has a control on the drainage pattern. In contrast, along other passive margins such a relationship is not clearly expressed, and the data points show a considerable scattering around the theoretical prediction (Figure 5b). We suggest that the scattering is in part due to the non-deterministic nature of surface processes, but more importantly, the developments of axial rivers may be also controlled by other factors besides the isostatic response to the erosional unloading, such as active tectonics and sediment deposition offshore. To examine the potential link between the offshore sediment loading and the location of an axial river, we compiled the sediment volume in a 100-km-wide swath from the river outlet to the margin based on the model of Straume et al. (2019), which will be discussed in the following section.

Effects of Sediment Loading and Other Factors
From the numerical modeling results, we observe that contours of the highest density of models are systematically above the theoretical curve, suggesting that locations of axial rivers may be further controlled by isostatic response to erosion over a region wider than the initial escarpment. This may reflect the erosion of the continent interior when the axial river started to form, or that in some cases, the continuous retreats of the escarpment and the locus of the highest erosion rate displaced the axial river landwards.
In addition to the flexural uplift induced by erosion of the escarpment, we consider sediment deposition to be another important process that affects the formation and location of the axial rivers. First, sediment deposition in the drainage basin, especially on the alluvial plains in the downstream part of the large rivers, complicates the distribution of mass unloading and thus the pattern isostatic rebound on the margin. Compared to the modeling results with detachment-limited rivers, models with onshore deposition predict a more dispersed locations of the axial rivers (Figure 5b), suggesting that redistribution of sediments contributes significantly to the variability of the data. Second, the accumulation of marine sediments on the margin has a flexural effect that is opposite in sign to the erosional unloading, and thus, the location of the resulting peripheral bulge associated with sediment loading can also be predicted using Equation 5. Therefore, the combined isostatic effect due to both the erosional unloading and the sediment loading on the margin will drive the locations of axial rivers further inland from the distance from the margin.
Compilation of the sediment accumulation on the shelf show that on a margin with high sediment volume (e.g., the Amazon; Figure 5b), the axial river usually occurs further inland.
The scattering of data from natural rivers is also contributed by the efficiency of the surface processes. We notice, for example, that due to the low magnitude of flexural depression, the drainage network reorganization in a substantial amount (<10%) of model simulations yield no axial rivers, which is consistent with our observations that some large rivers on passive margins (e.g., the Colorado River in Texas) do not have major axial channels. Furthermore, contrasting rock erodibility between different lithological units can also lead to the formation of axial rivers. For example, in the drainage basins of the Ogooue and Mahanadi rivers, the axial river channels developed in the relatively weak metasedimentary belts that are parallel to the margin, which are the West Congo belt (Thiéblemont et al., 2018) and Khondalite zone (Dobmeier & Raith, 2003), respectively.
The drainage patterns of some rivers may also be perturbed by regional tectonics, even though the rivers do not flow towards the convergent margins. A clear outlier in our analysis is the Rio Grande. We suspect that, despite having its outlet along a passive margin on the Gulf of Mexico, geometry of the river network is likely controlled by tectonic uplift as the majority of its upstream course lies in or near the southern Rocky Mountains. Along the east coast of Africa, development of the East African Rift lifted the topography and created grabens, which have influenced the drainage patterns of rivers in the area (e.g., the Zambezi and Limpopo).

Long-term Stability of the Axial Rivers
The modeling results also suggest that the impact of lithospheric flexure on river networks can be preserved tens of millions of years after rifting (Figure 2). Once an axial river develops, it captures the drainage areas of several upstream inland rivers, therefore aggregating their discharge/drainage. This results in high erosion rates and the entrenchment of its position in the landscape, a process amplified by the isostatic uplift of surrounding areas. Consequently, the position of the axial river is likely to remain at its initial position set during continental break-up even in situations where the initial topographic relief has been strongly degraded or valley growth has extended further into the plateau.
A clear example of the long-term stability of axial river channels is expressed in the present-day geometry of the Cuanza River. The river had lost a significant portion of its drainage area to the present-day neighboring Congo River during a period of drainage reorganization , which was potentially driven by sub-lithosphere mantle upwelling (Lithgow-Bertelloni & Silver, 1998) or late marginal flexure . Yet, segments of its axial channels have been preserved at locations consistent with our predictions ( Figure 5), perhaps set during the late Cretaceous following rifting along the Atlantic margin of Africa.

Potential Impacts
The formation of axial rivers causes continental-scale river drainage reorganization, which has led to the development of most of the largest river systems on the Earth. This, in turn, has had strong consequences on many aspects of the Earth's system. For example, the routing of continental freshwater export to the ocean affects the relative salinity of the large ocean basins, which controls planetary-scale ocean circulation (Maffre et al., 2018;Sinha et al., 2012). Without lithospheric flexure, the freshwater flux to the ocean would be more evenly distributed along continental coastlines, which would likely decrease the salinity gradient between ocean basins and lead to more sluggish global ocean circulation.
The geometry and evolution of continental drainage systems are also known to exert a strong control on the biodiversity and species richness within freshwater ecosystems. For example, it has been documented across many river systems that the number of freshwater species per unit area increases with the river discharge (McGarvey & Ward, 2008;Xenopoulos & Lodge, 2006). Connecting smaller rivers into a larger drainage basin can also potentially increase fish species richness. For example, higher fish species richness has been observed downstream of confluences (Fernandes et al., 2004;Osborne & Wiley, 1992), and tributaries are also known to support more species than same-sized rivers directly connected to the ocean (Oberdorff et al., 1997). The stability of large river networks is also conducive to higher freshwater fish species richness (Griffiths, 2018). Therefore, our findings that continental river network geometry is preserved long after rifting suggests that lithospheric flexure has facilitated species diversification.