Drainage integration and sediment dispersal in active continental rifts: A numerical modelling study of the central Italian Apennines

Progressive integration of drainage networks during active crustal extension is observed in continental areas around the globe. This phenomenon is often explained in terms of headward erosion, controlled by the distance to an external base‐level (e.g. the coast). However, conclusive field evidence for the mechanism(s) driving integration is commonly absent as drainage integration events are generally followed by strong erosion. Based on a numerical modelling study of the actively extending central Italian Apennines, we show that overspill mechanisms (basin overfilling and lake overspill) are more likely mechanisms for driving drainage integration in extensional settings and that the balance between sediment supply vs. accommodation creation in fault‐bounded basins is of key importance. In this area drainage integration is evidenced by lake disappearance since the early Pleistocene and the transition from internal (endorheic) to external drainage, i.e. connected to the coast. Using field observations from the central Apennines, we constrain normal faulting and regional surface uplift within the surface process model CASCADE (Braun & Sambridge, 1997, Basin Research, 9, 27) and demonstrate the phenomenon of drainage integration, showing how it leads to the gradual disappearance of lakes and the transition to an interconnected fluvial transport system over time. Our model results show that, in the central Apennines, the relief generated through both regional uplift and fault‐block uplift produces sufficient sediment to fill the extensional basins, enabling overspill and individual basins to eventually become fluvially connected. We discuss field observations that support our findings and throw new light upon previously published interpretations of landscape evolution in this area. We also evaluate the implications of drainage integration for topographic development, regional sediment dispersal and offshore sediment supply. Finally, we discuss the applicability of our results to other continental rifts (including those where regional uplift is absent) and the importance of drainage integration for transient landscape evolution.


| INTRODUCTION
In many continental settings undergoing active extension river network geometries change considerably over time (e.g. Jackson & Leeder, 1994;Leeder & Jackson, 1993). An often-observed trend is the progressive development of fluvial connections between initially isolated, endorheic, drainage basins and the eventual formation of a regional drainage network (e.g. Connell, Hawley, & Love, 2005;D'Agostino, Jackson, Dramis, & Funiciello, 2001;Dickinson, 2015;Duffy, Brocklehurst, Gawthorpe, Leeder, & Finch, 2015;Menges, 2008;Smith, 2013). This phenomenon, so-called drainage integration, explains why lake sediments often characterise older parts of the stratigraphy of fault-bounded extensional basins, whereas fluvial sediments are observed higher up in the record (e.g. Cavinato & DeCelles, 1999;Connell et al., 2005;D'Agostino et al., 2001;Miccadei, Piacentini, & Barberi, 2002). An area where drainage integration has clearly occurred is the central part of the Italian Apennines (Figure 1), which has been affected by active extension since approximately the beginning of the Pleistocene (Cavinato & DeCelles, 1999;Roberts & Michetti, 2004). While lakes were widespread during the Early-Middle Pleistocene in this area, most of them disappeared in the course of the Middle-Late Pleistocene as tectonic basins became progressively fluvially connected (e.g. D' Agostino et al., 2001;Piacentini & Miccadei, 2014). Understanding the mechanisms that control drainage integration is clearly important for interpreting the stratigraphic record preserved in such extensional settings.
In the central Apennines, drainage integration has previously been explained in terms of headward erosion from the coast, i.e. the capturing of basins at higher elevations by major streams that enlarge their catchments in an upstream direction (D'Agostino et al., 2001). However, there are other mechanisms that can lead to drainage integration between adjacent extensional basins. Drainage integration may partly be explained by the structural evolution of normal fault systems as adjacent fault segments propagate and link (Cowie, 1998). This leads to the structural lowering of topographic thresholds between these basins so they can become fluvially connected in an along-strike direction (Connell et al., 2005;Gawthorpe & Leeder, 2000;House, Pearthree, & Perkins, 2008;Menges, 2008). Another structural mechanism allowing integration to occur can be a reduction of fault slip rates over time (Connell et al., 2005). However, for explaining drainage integration across-strike and at a regional scale, as observed in the central Italian Apennines, additional mechanisms based on the dynamics of the fluvial system itself are required. Besides headward erosion (e.g. D' Agostino et al., 2001;Dickinson, 2015), other important mechanisms, proposed mainly for other areas, are the spilling over of lakes (e.g. Bishop, 1995;Douglass, Meek, Dorn, & Schmeekle, 2009;Garcia-Castellanos, Verges, Gaspar-Escribano, & Cloetingh, 2003;Smith, 2013) and the complete infilling of tectonic basins with sediment (e.g. Bishop, 1995;D'Agostino et al., 2001;Douglass et al., 2009). Although we have a fairly good understanding of these different mechanisms at a local scale, i.e. for individual basins, many fundamental questions remain regarding the conditions under which the different mechanisms may dominate and the impact of drainage integration on landscape evolution, sediment dispersal and, ultimately, basin stratigraphy in continental rifts (Smith, 2013).
There are additional reasons why improving our understanding of drainage integration is important. First of all, it forms a key aspect of transient landscape development in extensional settings but has, in contrast to the evolution of normal fault systems, received surprisingly little attention (e.g. Bishop, 1995;Stokes, Mather, & Harvey, 2002). Secondly, drainage integration has a profound impact on the volumes and characteristics of sediment supplied to tectonic basins (e.g. Smith, 2013). Thirdly, through its impact on sediment dispersal and hence mass redistribution, it is of great relevance for studies on the feedback between surface processes and tectonics in extensional settings (e.g. Buiter, Huismans, & Beaumont, 2008;Maniatis, Kurfeb, Hampel, & Heidbach, 2009). However, studying drainage integration in the field is complicated due to poor preservation of evidence. This is because drainage integration generally produces a wave of erosion in response to base-level changes. To overcome the problem of limited field evidence, we investigate the processes of drainage integration by means of numerical modelling. We use a simple model setup that includes the main features of tectonic deformation in the central Apennines to drive surface processes through time. By applying our modelling approach to this area, we make use of a wealth of field observations for calibrating our model and for evaluating our results. While previous modelling studies have demonstrated aspects of drainage reorganisation in rifts at a local scale (Cowie et al., 2006;Douglass & Schmeekle, 2007;Garcia-Castellanos et al., 2003;Smith, 2013), this approach allows us to address the problem at a regional scale (>100 km), involving a large number of extensional basins and fault-blocks.
The elevated topography in the central part of the Apennines cannot be explained by crustal or lithospheric isostasy (Faccenna, Becker, Miller, Serpelloni, & Willett, 2014). However, a clear correlation exists between topography, surface uplift and regional extension rates, suggesting that uplift and extension are driven by the same underlying mechanism (Faure Walker et al., 2012). Although the exact mechanism is debated (see review provided by Faccenna et al., 2014) uplift and extension are likely related to either flow or buoyancy variations in the uppermost mantle and removal of mantle lithosphere (e.g. Bartolini, D'Agostino, & Dramis, 2003;Cowie, Scholz, Roberts, Faure Walker, & Steer, 2013;D'Agostino & McKenzie, 1999;D'Agostino et al., 2001;Faccenna et al., 2014).
As first discussed by D' Agostino et al. (2001), many field observations demonstrate the combined impact of uplift and faulting on the geomorphologic development of the central Apennines and on the evolution of the drainage network (see also Ascione et al., 2008;D'Alessandro, Miccadei, & Piacentini, 2003, 2008. A key observation is that most of the major fault-bounded basins contain lake sediments in the older parts of their stratigraphy (Cavinato, 1993;Cavinato & DeCelles, 1999;Cavinato et al., 1994Cavinato et al., , 2002Miccadei et al., 2002). Based mainly on these sediments, it has been concluded that many large lakes coexisted during the Lower-Middle Pleistocene suggesting that endorheic drainage was prevalent at that time (D'Agostino et al., 2001;Piacentini & Miccadei, 2014). Today, most of these basins are fluvially dissected and connected to one another and to the coast. In other words, a temporal transition is inferred to have occurred from internal to external drainage leading to the integration of previous isolated basins with the regional river network (Bartolini et al., 2003;D'Agostino et al., 2001;Piacentini & Miccadei, 2014). Developing a better understanding of this transition via numerical modelling is the focus of this study.

| METHODOLOGY
For simulating regional landscape evolution for the setting of the central Apennines, we use the surface process model CASCADE developed by Braun and Sambridge (1997). Its suitability has been demonstrated for modelling landscape development in extensional settings, where both fluvial erosion and deposition occur and where lakes are common features in the landscape (Cowie et al., 2006). There is a one-way coupling in our model in that we allow surface processes to respond to surface deformation due to tectonics, but there is no feedback of surface processes on the tectonics. Besides extensional faulting, our model also includes regional uplift, and both are simulated by means of simple surface deformation functions (see below).
The model domain covers all land area between the modern coastlines in central Italy (Figure 1b). The region is rotated 45°clockwise relative to true North so that the dominant SW-NE direction of extension coincides with the x-direction in the model domain (Figure 1a,b). The model domain is 170 9 170 km and has a spatial resolution of 1 km in both directions (Table 1). The left and right boundaries of the model domain represent the Tyrrhenian and Adriatic coastlines respectively. These coastal boundaries are fixed in order to keep base-level constant, as climatically induced sea-level oscillations are small compared to the tectonic deformation we impose. We return to this assumption in the Discussion section. The other two boundaries of the model domain delimit our study area in the along-strike direction, i.e. along the Apennines (y-direction in the model), and are free to slip vertically. All four boundaries are open in the sense that water and sediment Scaling factor for channel width W = c 9 Q / 1 -can cross them. There is a free surface above for enabling topography to develop and vertical surface displacements are imposed from below. We run all our experiments for 3 million years, i.e. the estimated duration of extension (Roberts & Michetti, 2004). Although some authors suggest that regional uplift may have commenced more recently (e.g. Pizzi, 2003), we impose regional uplift from the beginning of the model runs for the sake of simplicity. The calculation time step in the model is adjusted dynamically but is ca. 100 years on average. We do not assume any pre-existing topography, except for 1-m-scale random noise to initiate flow, even though the central Apennines were likely characterised by some relief at the time extension commenced (e.g. D' Alessandro et al., 2003). This means that there is no inheritance effect on drainage network development. We evaluate the potential implications of our zero pre-existing topography assumption in the Discussion section.

| Normal faulting surface deformation and regional uplift function
For our calculations of vertical surface deformation in response to normal faulting, we use the elastic dislocation model Coulomb 3.4 (Lin & Stein, 2004;Toda, Stein, Richards-Dinger, & Bozkurt, 2005), which is based on linear elasticity laws and a half-space assumption (Okada, 1992; for more details see Appendix S1). The main input to the elastic dislocation model is a fault map that includes all normal faults thought to have accommodated extension in the central Apennines since the Early Pleistocene (principally based on Michetti, 2004 andWedmore et al., 2017). Except for some faults located in the southwestern part of the area, they are all considered as active today and throughout the modelling period ( Figure 1b). In order to focus on the main topographic features only, the fault map was simplified by removing faults shorter than 5 km and by straightening the fault traces (compare Figures 1b and 2a). The simplified fault map comprises 50 faults with lengths between 5 and 40 km. Nearly all faults dip to the southwest (towards the left in the model domain) and we assume pure dip-slip for all of them. Rarely observed minor strike-slip motions do not contribute to relief and are thus ignored. Table S1 (Appendix S1) shows the parameter values used in the elastic dislocation model. Parameters for which no field area-specific data exist are assigned published values and are kept constant in all our calculations (Poisson's ratio, Young's modulus and coefficient of friction). The three fault-related parameters dip angle ("dip"), fault root depth ("root"), and a linear fault length-displacement scaling factor (c or "gamma") are most important for our study. The latter scales maximum fault displacement D (experienced by the central part of the fault) linearly to fault length L as given by: D = c 9 L (Cowie & Scholz, 1992). For each of these parameters, we test the impact on the vertical surface displacement field for three different values (two extremes and one intermediate value) based on published data from the central Apennines (Table S1; Appendix S1). The parameter c has the greatest impact on the vertical surface displacement field. Our intermediate value for c (0.07) produces total throws which correspond best to those estimated in the field (Roberts & Michetti, 2004), and we use this surface deformation field as our standard faulting scenario in all of our experiments (Figure 2a). For transforming the fault map into surface deformation rates used in the landscape evolution model, we divide the total uplift and subsidence values by 3 million years (see Figure 2a for the resulting uplift and subsidence rates). This implies that fault offsets accumulate linearly over time (with c = 0.07, maximum uplift and subsidence rates are ca. 0.24 and À0.54 mm/year, respectively). Because field evidence suggests that some faults in the central Apennines experienced an increase in slip rate around 0.5-1.0 Ma (Cowie & Roberts, 2001;Roberts & Michetti, 2004;, we address the potential implications of changes in fault slip rates in the Discussion section. We simulate long-wavelength regional uplift across the mountain belt and its forelands using a Gaussian function (coast-to-coast transect; Figure 2b). In the direction parallel to the mountain range, i.e. parallel to the y-axis in our model domain, we assume regional uplift to be uniform. We scaled our Gaussian function based on published field observations and some new estimates of regional uplift for the mountain range interiors, in order to obtain the right order of magnitude of total Pleistocene plus Holocene uplift ( Figure 2b). Because of the limited number of welldated regional uplift estimates and their considerable spatial variability across our study area, we emphasise that our regional uplift function is only a first-order approximation. However, most important for our modelling study is that it accounts for the strongest uplift in the mountain range interiors and a gradual decline when moving across the foreland areas towards the Tyrrhenian and Adriatic coastlines. The published data that we used for constraining our uplift function in the foreland areas are paleoshorelines and exposed shoreface deposits (e.g. D' Agostino et al., 2001;Pizzi, 2003). The data from the different sites are provided in Figure 2b and described in more detail in Appendix S2 (Table S2). Besides these published observations, we use structural data from four normal faults to provide some additional constrains on our regional uplift function in the interior part of the central Apennines. For these four normal faults (see Figures 1b and 2a for their locations), we estimate the amount of uplift of their fault planes by assuming typical long-term ratios of footwall uplift to hanging-wall subsidence and by assuming that the land surface was close to sea-level before regional uplift started. A detailed description of our method, the data, and our regional uplift estimates are provided in Appendix S3 (Table S3). These new uplift estimates suggest a total amount of regional uplift of around 1,000 m in the innermost part of the central Apennines that corresponds well with reconstructions made by others (Ascione et al., 2008;Pizzi, 2003). It is important to note that we use a F I G U R E 2 (a) Vertical surface displacement map produced by the elastic dislocation model Coulomb 3.4 based on a simplified fault map from the central Apennines (modified from Wedmore et al., 2017). It shows our "standard faulting scenario" using dip = 60°, root = 15 km and c = 0.07 (see Methodology section and Appendix S1). This displacement field is assumed to represent the accumulated impact of normal faulting after 3 Myr. Uplift and subsidence rates (mm/year) are the total uplift and subsidence values divided by 3 million years. (b) Regional uplift curve showing the total amount of long-wavelength surface uplift along a coast-to-coast transect projected on top of a 20-km-wide topographic swath (in grey) across the central Apennines (see Figure 1b for swath location). Regional uplift rates (mm/year) are the total regional uplift values divided by 3 million years (see vertical axis on the right). Also shown are localities and elevations of field observations that were used to constrain the amplitude and shape of the regional uplift function. These observations comprise four different sea-level markers (M1-M4; see also Figure 1b and Appendix S2) and the localities of four faults (FiF, FuF, BaF and SuF) where the amount of regional uplift was estimated (see also Figure 1b and Appendix S3) symmetrical uplift function in most of our numerical experiments even though some studies suggest the Adriatic flank of the mountain range may have experienced more uplift than its Tyrrhenian counterpart (e.g. Pizzi, 2003). We assume a symmetric function for simplicity and because there seems to be no general agreement about the exact pattern of regional uplift. The potential implications of this assumption are addressed in the Discussion section. Regional uplift rates are kept constant through time in our model (see Figure 2b for regional uplift rates).

| Surface process model
We use CASCADE (Braun & Sambridge, 1997) for simulating fluvial erosion and sediment deposition in lakes ( Table 1). The fluvial erosion algorithm follows the "under-capacity model" and can generate both erosion and deposition (Kooi & Beaumont, 1996;van der Beek & Bishop, 2003): where dh dt is elevation change. Transport capacity Q c is the volume of sediment that is theoretically possible to be carried by the flowing water and its magnitude depends on discharge Q w and local channels slope S: This linear dependency is scaled by the dimensionless transport capacity constant K f . The sediment volume Q s in Equation (1) is determined by integrating all the elevation changes that are occurring upstream and represents the sediment passed to every node in each time step: where A is the total upstream drainage area and da is the downstream increment of upstream area. According to Equation (1), the rate of erosion or deposition dh dt is primarily a function of the disequilibrium between the transport capacity Q c of the river and the volumetric sediment flux Q s . If Q c > Q s there is erosion, if Q c < Q s there is deposition, and the difference between them controls the rate of erosion or deposition. However, erosion and deposition rates are additionally controlled by the width of the channel W and the fluvial length scale parameter L f , which both reduce erosion rates as their values increase. Because of the large dimensions of our study area, we assume channel width to vary as a function of discharge (W ¼ ffiffiffiffiffiffi Q w p ). Both parameters K f and L f affect the erosive conditions in our model. Simply stated, higher values for K f generate higher erosion rates and vice versa, whereas lower values for L f generate higher erosion rates and vice versa. However, as discussed in detail by Cowie et al. (2006), L f additionally controls the way in which rivers respond to changes in base level, either in a more transport-limited or in a more detachment-limited manner. We systematically varied K f and L f between 0.08-0.12 and 30-70 km, respectively, in order to test the sensitivity of our model (see Appendix S4). We do not consider spatial lithological differences and temporal changes in climate in this study, and K f and L f are consequently kept constant in space and time. We address the potential implications of assuming a uniform lithology in the Discussion section. Climate variability is out of the scope of our study as it is not possible to resolve its crucial aspects (e.g. storm intensity) on geological timescales (e.g. Whittaker, 2012). Land-sliding is locally important for landscape evolution in the central Apennines (Whittaker, Attal, & Allen, 2010) but we do not include it because the spatial resolution (1,000 m) of our regional-scale model means that no slopes exceed the critical angle for landslide initiation (typically ≥21°). The fluvial algorithm in CASCADE does not distinguish fluvial channels from the interfluve areas and thus erosion occurs across the entire landscape not only along channels.
Important for this study is the treatment of water and sediment when a stream enters a local minimum in an extensional basin. First of all, the model calculates the lowest point on the rim of the basin (i.e. the spill-point) and defines all nodes in the basin at lower elevation as lake nodes. All sediment entering a basin is trapped as long as the basin is underfilled and supports a lake. The sediment is deposited in nodes closest to the river mouth, causing basins to become progressively filled from their edges. With regard to water conservation we simulate truly endorheic drainage, i.e. closed basins where water loss through evaporation or seepage (including karst) exceeds water supply. This is chosen because at least two large lakes in central Italy, i.e. the historical Fucino lake (which is now artificially drained) and the Trasimeno lake (Umbria; Ludovisi, Gaino, Bellezza, & Casadei, 2013), demonstrate the occurrence of truly endorheic drainage under modern-day (interglacial) climatic conditions. Additionally, some studies on Italian lakes have demonstrated the important role of evaporation in controlling their hydrological balance also in glacial times (e.g. Zanchetta, Borghini, Fallick, Bonadonna, & Leone, 2007). Finally, by comparing model experiments in which we implemented either endorheic or nonendorheic (water 100 % conserved) drainage, we found that characteristic topographic features of the central Apennines and important aspects of its evolution are only reproduced by means of the endorheic type of drainage (see Appendix S5).

| Topographic development
Here we present results mainly from our "reference model" (using K f = 0.10 and L f = 50 km) as it shows the general behaviour of the system we model. Varying erosional conditions produces slightly different patterns and rates of landscape development but does not change the main trend of landscape evolution (see Appendix S4). The surface displacement field (faulting [ Figure 2a] and regional uplift [ Figure 2b] together) produces +1,600 m and À900 m of maximum uplift and subsidence, respectively, corresponding to maximum rock uplift and subsidence rates in between À0.3 and +0.6 mm/year over 3 Myr. The steadystate concavity of major river systems crossing both the faulted domain and the foreland area lies between ca. 0.35 and 0.6 for the K f and L f values used in our reference model (Appendix S6). This range encompasses concavity values that are typical for steady-state river profiles in general and also corresponds well with those observed in the central Apennines .
Figure 3a-d illustrate over four time steps how the topography evolves through time in the reference model. Initially, elevations remain low everywhere (<ca. 500 m during the first 1.5 Myr of run-time), since we do not assume any preexisting topography. However, with time, mean elevations in the central part of the model domain increase as a consequence of regional uplift (Figure 3a-d). Our reference model produces just over 1,000 m of topography after 3 Myr runtime (Figures 3d and 4). A local-scale morphology of longitudinal ridges and basins develops due to the normal faulting superimposed on the regional topography (Figures 3d and  4). This gradual increase in relief at two different spatial scales (regional vs. local scale) is characteristic of the topographic development in our model and is consistent with the topography of central Italy today (see Discussion section). While the final regional relief is approximately 1,000 m, the local-scale (ca. 10-20 km), fault-related relief is of the order of hundreds of metres, but varies greatly throughout the model (top panel of Figure 4). This large spatial variation in fault-related relief in our model is caused by variations in fault length, fault spacing, the orientation of faults relative to one another, and the position of faults relative to the regional uplift field. This is because surface deformation at any location in our model is the sum of all surface deformation fields produced by the individual faults plus the regional uplift field (Figure 2). Local relief is additionally affected by the degree of basin infilling. Because most basins experience, successively, sedimentation and incision, the degree of infilling is strongly time-dependent. Another striking feature of the final topography is its asymmetry (higher topography on the Adriatic side) even though our regional uplift function is symmetrical (bottom panel of Figure 4). This asymmetry is partly due to the SW preferential fault dip in combination with the relative small fault spacing (so that the uplift-subsidence fields of individual faults overlap), generating higher faultrelated topography on the Adriatic side ( Figure 2a). However, as discussed below in "Regional-scale sediment F I G U R E 3 Time evolution maps from our reference model (see main text) showing the main landscape features after 1.5, 2, 2.5 and 3 Myr.  Figure 3d) together with normal faults (schematic) and basin deposits. Vertical arrows demonstrate regional-and local-scale relief (see main text). The regional uplift function is shown as reference in the bottom panel (see also Figure 2b) dispersal", the asymmetry in topography additionally results from different rates of erosion and overall landscape evolution between the Adriatic and Tyrrhenian domains.
Both spatially averaged mean and maximum elevations continue to increase even at the end of each model run ( Figure 5a). In other words, the landscape does not reach a topographic steady state within the 3 Myr time period we consider here. This is consistent with the transient landscapes observed today in the central Apennines (e.g. Whittaker et al., 2008) and ongoing surface uplift (D'Anastasio et al., 2006;Serpelloni et al., 2013). In our reference model, steady state is reached approximately after 6 Myr, i.e. after twice the normal model run-time. In the central Apennines today elevations can exceed 2,000 m, whereas in the model the highest elevations are around 1,000 m. This difference can be attributed to pre-existing topography, something we come back to in the Discussion section.

| Drainage network evolution
At the beginning of the experiments, small stream networks initiate over the entire model domain. A large number of lakes form particularly in the faulted domain where local topographic minima develop in the hanging-wall basins. Each of the drainage basins that support lakes are endorheic, i.e. internally drained (see "Surface process model"). The lakes act as local base-levels and trap all the sediment delivered from upstream. Initially, the whole area affected by normal faulting is internally drained, i.e. ca. 40%-50% of the total model domain (Figures 3e and 5b). However, through time we observe a consistent trend of progressive integration of the drainage network, resulting in the disappearance of lakes and shrinkage of the total endorheic area (Figures 3e-h and 5b). Although both lake and endorheic area show a progressive change over time, it is important to note that the total surface area occupied by lakes ("total-lake-area") declines in a different way compared to the total area that is internally drained (compare Figure 5b,c). The total-lake-area shrinks from the beginning of the model run, with the most drastic decline occurring during the first 1.5 Myr of the experiment (from ca. 24% down to ca. 7% of the total model domain, see Figure 5c). On the other hand, the total endorheic area remains fairly constant until 1.5 Myr and successively shrinks in a step-wise manner (Figure 5b). The reason why the total-lake-area decline is so different from that of the endorheic area (Figure 5b,c) is because the extent of the endorheic area is determined by the presence of lakes most proximal to the coast. For instance, the westernmost basin (e.g. Figure 3e) keeps the Tyrrhenian flank internally drained until ca. 1.6 Myr although many lakes upstream have already disappeared. The transition from internal to external drainage means that sediment produced in the upland area is henceforth transported out of the faulted domain, and thus exported to the coast, at localities that we define as fluvial "exit points" (Figure 3e-h).
Characteristic of the drainage network in general is the strong contrast in drainage network geometry, within and outside the central area affected by normal faulting. Outside the faulted domain the network has a parallel to slightly dendritic appearance, formed by channels that follow the regional slope of the land surface towards both coastlines (Figure 3). Within the central area, however, many streams or stream segments flow axially, parallel to fault strike, forming a trellis-like drainage pattern (Twidale, 2004). The planview geometry of the drainage network and the position of the central drainage divide are established early on and remain fairly stable over time (Figure 3). The position of this drainage divide is controlled by the regional uplift field (see also Appendix S7).

| Drainage integration
The dominant mechanism that causes drainage integration in our model is what we call "basin overfilling". By this we mean the filling of basins with sediment up to the elevation of their spill-point, i.e. the lowest point on their morphological boundaries. When a basin becomes overfilled with sediment the water can spill over and the lake environment is replaced by a through-going river system (Figure 6a). From this moment onwards, some sediment is still deposited within the basins to balance newly created accommodation due to fault-controlled basin subsidence, but most sediment is now transported downstream towards other basins (Figure 6a,c) or all the way to the coast. The long-term regional-scale tendency of basins to become overfilled demonstrates that sedimentation rates gradually start to outpace the rate at which accommodation is created through basin subsidence. This happens over time as mean erosion rates increase owing to an increase in both fault-related and regional relief (Figure 5e). This increase in mean erosion rates, in turn, causes a gradual shift in the balance between F I G U R E 6 (a) Maps and cross-sections of two fault-bounded basins in the model illustrating the gradual filling of basins through time and the mechanism of basin overfilling (for location see Figure 7a). Basin I and basin II become overfilled around 1.3 and 2.2-2.8 Myr respectively. I.G. = Interior gorge, located in between basin I and II. (b) Lake volume curves of basin I and II shown in (a). These curves show that basins initially tend to become more undersupplied but later on less undersupplied due to constant fault slip rates but increasing erosion rates. (c) Sediment supply to basin I and II, showing that overfilling of basin I generates a sudden increase in sediment supply to basin II sediment supply and accommodation creation within the basins. Using lake volume as a proxy for how undersupplied a basin is, we can demonstrate this shift (Figure 6b). Lake volumes firstly tend to increase, meaning that the basins become increasingly undersupplied. However, this trend reverses as soon as sediment supply outpaces accommodation creation causing the lake to shrink and the supporting basin to become progressively less undersupplied. It is important to note that each individual basin/lake follows its own curve (Figure 6b). In our reference model, the total volume of all lakes together increases until circa 1.2 Myr and successively decreases thereafter (Figure 5d). The order in which the individual basins become overfilled does not follow any clear spatio-temporal pattern. For instance within the Tyrrhenian part of the chain interior, lakes with either a proximal or distal location relative to the coast disappear early on in time (e.g. Figure 3f). Moreover, the longest surviving endorheic basin on the Tyrrhenian side has an intermediate position and is not located closest to the central drainage divide (Figure 3g). A clear spatio-temporal pattern is lacking because basin overfilling is a function of a large number of local factors that affect the balance between sediment supply and accommodation creation. The rate at which accommodation is created is not only a function of fault length and slip rate but is also affected by the position of faults relative to one another. Sediment supply on the other hand is controlled by the size of the source area and its internal relief, which are also strongly controlled by the pattern of faulting. Furthermore, sediment supply to individual basins depends on the infilling histories of basins located upstream.

| Regional-scale sediment dispersal
Erosion rate maps (Figure 7a-c) show the general pattern of erosion and deposition in our reference model: Sediment F I G U R E 7 (a-c) Erosion-deposition maps showing the total amount of erosion and deposition that occurred in the model during 100 kyr periods, namely 0.9-1 Myr, 2.4-2.5 Myr and 2.9-3 Myr. Yellow stars correspond to the spill-point of Basin II in Figure 6a. "low e" = reduced erosion within the endorheic area. (d) Cumulative erosion (for the total 3 Myr time period) summed along-strike and projected on a coast-to-coast transect. The transect has been divided into a Tyrrhenian flank domain (dark-shaded zone on left-hand side), a faulted domain (white zone in the middle) and an Adriatic flank domain (dark shaded zone on right-hand side) based on the extent of the area affected by normal faulting. The two light-shaded zones are transition zones owing to the 3D geometry of the fault array that is projected on a 2D cross-section. "low e" = reduced erosion because of endorheic drainage (see also Figure 7c). Percentages show the relative contribution of the total Tyrrhenian and Adriatic domains (both including a faulted domain and mountain flank part) to the overall amount of eroded material delivered to the coastlines is mainly produced at the footwall highs and along the major river valleys and is deposited in the fault-controlled basins or offshore (outside the model boundaries). However, the three different time windows (Figure 7a-c) also show how sediment dispersal changes over time. The most important regional trends are: (i) The gradual increase in erosion rates and hence sediment production due to increasing relief, and (ii) the progressive decline in deposition in tectonic basins due to drainage integration (e.g. compare Figure 7a,c). Because less sediment is trapped within the basins over time, progressively more sediment becomes removed from the faulted domain as the landscape evolves. In other words, there is a delayed export of sediment out of the mountain range towards the offshore.
At a local scale, on the other hand, we observe abrupt shifts between erosion and deposition. These shifts are again related to drainage integration that acts as a threshold phenomenon. While most basins firstly experience a relative stable phase of lake sedimentation (e.g. Figure 7a), they abruptly switch to a fluvial environment with strong incision as soon as they become overfilled. Incision initiates in the area of the spill-point as a new base-level is established at a lower level and there is an abrupt increase in discharge (Figures 7b and 8a) as the fluvial system becomes connected. Lowering of the spill-point, in turn, generates a wave of erosion that starts to propagate upstream and deeply dissects the basin fill (Figures 7c and  8b). In other words, in our model, sedimentary basins themselves and their spill-point areas are most prone to abrupt local changes in erosion or deposition. All the surrounding terrain successively adapts in a more gradual manner. However, it is important to note that these local developments, due to drainage integration, strongly affect the downstream parts of the catchment. For instance when a basin becomes overfilled, the sediment is no longer trapped and is henceforth transported to another basin downstream (Figure 6a,c). Additionally, strong incision commences in the valley in between the basins causing the sediment supply to the downstream basin to become enhanced even more. As such, the infilling history of each basin is a function of the infilling histories of all the other basins located upstream.
After 3 Myr of landscape evolution, the Tyrrhenian and Adriatic domains (i.e. measured from the central divide) have experienced approximately the same amount of erosion (respectively 49% and 51%; Figure 7d), implying that both offshore areas have received similar sediment volumes overall. However, Figure 7d clearly shows that erosion is most intense on the Adriatic flank of the central Apennines in our model. This is because the SW preferential fault dip is opposite to the regional slope in the Adriatic domain, generating higher relief along the flank of the mountain range and thus higher erosion rates (Figure 9). The reason why this does not produce a higher total sediment output to the Adriatic offshore compared to the Tyrrhenian offshore, is that a large part of Adriatic faulted domain is still internally drained after 3 Myr in our model (Figure 7c). This latter effect can also be attributed to the structural setting of the Adriatic domain (fault dip opposite to regional slope) as it slows down basin overfilling and therefore drainage integration. In other words, within the Adriatic faulted domain, local relief and therefore sediment delivery to

Basin I Basin II
Spill-point of Basin II in Fig. 6 basin fill is overfilled I.G.
Interior gorge that formed due to overfilling of Basin I I.G. hanging-wall basins are relatively low, whereas the rate of accommodation creation is relatively high compared to its Tyrrhenian counterpart (Figure 9). The offshore as a whole, i.e. the Tyrrhenian and Adriatic coasts together, experiences a long-term progressive increase in sediment supply in our model (Figure 5f). At around 1.9 Myr, an abrupt increase in the offshore sediment flux is observed when most of the internally drained area on the Tyrrhenian flank becomes fluvially connected to the coast (Figure 3f). Every time a significant part of the faulted domain becomes externally drained, a step-wise increase is observed in the offshore sediment flux (Figure 5f).

| DISCUSSION
In this study, surface process modelling is used to investigate the impact of regional uplift and normal faulting on long-term landscape evolution across the central Apennines. Our model results enable us to improve our general understanding of drainage integration in extensional continental areas and allow field observations from the Apennines to be evaluated in a temporal perspective. The benefit of our study lies in the simplicity of our model set-up. However, it may not explain detailed field observations on a local scale.

| Model vs. observations
The drainage integration trend seen in our model explains the commonly observed transition from lacustrine to fluvial sedimentation in basin stratigraphy in the central Apennines, followed by strong incision of the basin fill (Cavinato, 1993;Miccadei et al., 2002;Pucci et al., 2014). While widespread lacustrine deposition characterised the Lower-Middle Pleistocene, progressively more basins became externally drained post late Middle Pleistocene (Piacentini For comparison, we also show the localities of the two fluvial "exit points" produced by our reference model ("S" and "T" in (d)) by means of pink dots. Also shown are the Pescara ("PES") and Fucino ("FUC") catchments. (c) and (d) show the same type of data as shown in (a) and (b) but from our reference model after 3 Myr. Yellow dots show the two fluvial "exit points" ("S" and "T") on the Adriatic side and in grey the area that is still endorheic after 3 Myr. The model catchments marked "X" and "Y" show alternative geometries for the real Sulmona and l'Aquila-Campo Imperatore catchments, and emerge in the absence of an influence from pre-existing topography (see also (e)). These model catchments are connected to the Adriatic foreland area through fluvial exit points "S" and "T". (e) Residual between the DEM and the surface displacement field used in our model. The residual is derived by subtracting the 3 Myr surface deformation field (including both normal faulting and regional uplift; Figure 2a & Miccadei, 2014). Our modelling results are therefore in general agreement with D' Agostino et al. (2001) in concluding that drainage integration in the central Apennines is related to the development of a topographic bulge in combination with normal faulting along its crest. However, our results allow us to investigate the processes controlling this transition in more detail. We compare the final topography (after 3 Myr) from our model with the Digital Elevation Model (Tarquini et al., 2007) from the central Apennines (Figure 10a,c). Both show a combination of long-wavelength and more local-scale fault-related topography, demonstrating the importance of both normal faulting and regional uplift for landscape evolution in the central Apennines. In addition, our model reproduces the observed strong tectonic imprint on the stream network (Figure 10b,d). The river network has a dominantly along-strike orientation within the faulted domain and exhibits a parallel drainage pattern in the tilted foreland areas. Although the modelled and observed stream networks overlap to great extent, the exact catchment geometry differs in detail (compare Figure 10b,d). This is because catchment geometry in extensional settings is strongly controlled by the localities where streams find their way across fault-related topography. Because these transverse reaches are sensitive to many factors that are not included in our model (e.g. pre-existing topography, lithological differences, fault propagation and rock damage, karst drainage, etc.), some are not exactly reproduced. An important example is the Popoli gorge that receives water from the large Pescara catchment, and is the locality where most surface water exits the faulted domain on the Adriatic side (Figures 10b, 11, and 12f). In our model, there are instead two smaller catchments, one supplying the Sulmona basin (catchment "Y" in Figure 10d) and the other one in the area around l'Aquila and Campo Imperatore (catchment "X" in Figure 10d). This is because the model predicts the presence of two main "exit points" instead of one near Popoli (Figures 10d and 11). Although this is an obvious mismatch, we believe it provides some interesting insights. First of all, our simple model setup always produces high topography in the area around Popoli instead of producing a relative low area that can become an exit point. This suggests that active tectonics alone probably cannot explain the Popoli gorge and another factor is needed to explain it, e.g. pre-existing topography (Figure 10e, see below) and possibly karstification processes (Boni, 2000). Secondly, the localities of the two exit points produced by the model actually do coincide with a deeply incised valley that receives water from Campo Imperatore and a large windgap in between Maiella and Sulmona ( Figure 11). Although this valley and windgap may have other explanations, our results clearly demonstrate that these two localities are favoured as potential exit points based on faulting and regional uplift only. Evidence for pre-existing topography is clear from the difference in maximum elevation (ca. 2,000 vs. 1,000 m) between our model and observations (note different colour bar scaling in Figure 10a,c, see also Figure 10f). In addition, there is an intermediate-scale morphology in the central Apennines consisting of 20-30 km wide ridges that cannot be explained by normal faulting alone (Figure 10a). These differences are clearly visible in our calculated residual topography (present-day topography minus our tectonic uplift function) shown in Figure 10e. Because the landscape morphology in Figure 10e lines up with mapped thrust faults it confirms that the central Apennines were likely characterised by significant thrust-related topography and deformation structures from the earlier phase of compression prior to Quaternary extension. In other words, our model results support the idea that inherited thrust-related topography has also contributed to the modern-day landscape, (e.g. D' Alessandro et al., 2003) and possibly influenced the extensional fault pattern (D'Agostino, Chamot-Rooke, Funiciello, Jolivet, & Speranza, 1998;Scisciani, Tavarnelli, & Calamita, 2002). However, here we show that inherited topography is not a necessary ingredient to produce drainage integration.
Even though local peak elevations >1,000 m are not reproduced in our model, the hypsometric distributions show a striking similarity between model and reality marked by a local maximum around 600 m (Figure 10f). In the model this local maximum cannot be explained only by the tectonic uplift function, as demonstrated by the hypsometric distribution produced by normal faulting and regional uplift only (pink line in Figure 10f). It can be explained, however, by the prevalence of internal drainage for a considerable part of the 3 Myr model time and the existence of local (perched) base-levels. As long as there is internal drainage, rivers transport material towards the altitude of their local base-level, leading to the development of a local maximum in the hypsometric distribution. The real local maximum in the central Apennines corresponds to the elevation of the internally drained Fucino basin, at circa 650 m. The Fucino basin remains internally drained today because there is insufficient sediment supply compared to the high rate of accommodation creation (see "Overspill vs. headward erosion from the coast" below). In our model run, the local hypsometric maximum corresponds to the local base-level elevation of the l'Aquila-Campo Imperatore area that is still internally drained after 3 Myr (Figure 10d). The primary reason why drainage integration in this area of the model is slowed down is the structural setting of the Adriatic part of the faulted domain where the dominant fault dip direction is opposite to the regional slope (Figure 9; see "Regional-scale sediment dispersal" in the Model results section).

| Overspill vs. headward erosion from the coast
Our model results demonstrate that an important mechanism driving drainage integration is basin overfilling, i.e. the progressive filling of basins with sediment up to the level of their spill-point enabling water to spill over (Figure 6a). We note, however, that our endorheic model setup (see Methodology section) does not allow us to distinguish basin overfilling from lake overspill, i.e. the spill over of water when the lake surface (and not the sediment surface) reaches the altitude of the spill-point. In theory the potential for lake overspill is mainly climate-dependent, likely making lake overspill more important under wetter climatic conditions (e.g. Garcia-Castellanos et al., 2003;Heidarzadeh, Ballato, Hassanzadeh, Ghassemi, & Strecker, 2017;House et al., 2008). Even though our model cannot distinguish between basin overfilling and lake overspill, we can consider them both as "overspill mechanisms" (Bishop, 1995;Smith, 2013), as they both act in a downstream or "top-down" direction and are mainly controlled by sediment and water supply from upstream. Therefore, in turn, we believe our model suggests that overspill mechanisms mainly drive drainage integration in the central Apennines (Figure 12a,c). This finding contradicts previous field-based studies on the central Apennines that suggest headward erosion from the coast to be the dominant driving mechanism (Figure 12b,d;D'Agostino et al., 2001;Bartolini et al., 2003), i.e. "bottom-up" fluvial integration (Bishop, 1995;Smith, 2013). We do observe headward erosion from the coast in our model, but its contribution to drainage integration is negligible, and this result is irrespective of the erosional parameters that we use (see Appendix S4).
It is important to note that there is no reason to expect the contribution of overspill in our model to be over-estimated relative to headward erosion. First of all, increased sediment supply to hanging-wall basins can only be generated under more erosive conditions, which in turn also increases headward erosion. Their relative importance thus remains the same and explains why overspill remains the dominant process driving drainage integration when varying erosional parameters K f and L f (Appendix S4). F I G U R E 1 1 Top figures: Google Earth images of the Popoli gorge and the two modelled fluvial "exit points" on the Adriatic side produced by our reference model (see also Figure 10b,d). Main fault systems (according to Roberts & Michetti, 2004) shown by means of red lines. Bottom figure: Google Earth image of the Adriatic foreland area and the central Apennines in the background, also showing the main river systems. Over a distance of 120 km only one river system penetrates far into the high topographic area (the Pescara river system), whereas most rivers drain the mountain range flank and foreland area only Secondly, we do not expect the dominant role of overspill to be related to major assumptions underlying our model setup. If lithology is not uniform, as we assume, overspill would most likely become even more important as the main lithological contrast between basin alluvium and the more resistant bedrock ridges would lead to more rapid excavation of sediment from the basins and more rapid incision at the spill-point directly following a drainage integration event (e.g. Cowie et al., 2008). Initiating the model with pre-existing topography, on the other hand, is likely to increase the rate of basin filling and hence overspill. Increasing the model resolution also would not affect our main results because of the strong control on the scale of the local relief exerted by the fault pattern. Climatically induced sea-level low stands could theoretically enhance headward erosion but are small compared to the tectonic uplift. Finally, using an asymmetric regional uplift function instead of a symmetric function does not affect the dominant role of overspill in long-term drainage integration, even though it produces a significantly different landscape after 3 Myr (Appendix S7).
There are a number of field observations from the central Apennines that point towards basin overfilling, or lake overspill, being an important process. The most direct evidence comes from the Terni basin in the northwest corner of our study area (Figure 1b). Here continental deposits are preserved on top of the adjacent Narnese-Amarina ridge, marking the location of the former outlet of the Terni basin. The high elevation of this former outlet relative to the present-day basin surface shows that the basin has been totally filled up to its spill-point and that basin overfilling caused it to become interconnected with the Tyrrhenian foreland area (D'Agostino et al., 2001).
Based on the position of the closed Fucino basin on the central drainage divide and at greatest distance to the Tyrrhenian and Adriatic coasts (Figure 10a), previous work has argued for headward erosion from the coast to be the main mechanism driving drainage integration (D'Agostino et al., 2001). D' Agostino et al. (2001) hypothesised that over time all other major and initially endorheic basins have been captured except for the Fucino basin, which has "survived" and remained internally drained because of its distal position relative to regional (marine) base-level. However, based on our results we suggest that the Fucino basin is internally drained today due to an insufficient sediment supply that has been outpaced by fast accommodation creation (see also Whittaker et al., 2008). Its stratigraphy does not support the "survivalconcept" as it shows a transition from overfilled to underfilled conditions over time (Cavinato et al., 2002) that can be explained by a 93 to 95 increase in slip rate at around 1-0.5 Ma along the main basin-bounding fault (Roberts & Michetti, 2004;Cowie & Roberts, 2001;Whittaker et al., 2008;Appendix S8). Our study suggests that the main reason why the Fucino basin is endorheic today is simply because its central position within the fault array caused the Fucino fault to become the largest and most active fault in the area (Cowie & Roberts, 2001;Roberts & Michetti, 2004). In other words, the preservation of the Fucino basin confirms the importance of accommodation creation vs. sediment supply and hence the major role of overspill mechanisms rather than headward erosion from the coast in controlling drainage integration.
Another reason why we do not expect headward erosion from the coast to be important for drainage integration in the central Apennines is the small number of fluvial connections of significant size between the foreland area and the interior of the mountain range. For instance, in the Adriatic domain only one such connection, i.e. the Popoli gorge, exists over a total along-strike distance of ca. 120 km, i.e. between the Sangro and Tronto river valleys ( Figure 11). Moreover, the young age of the Popoli gorge (ca. 400-350 ka according to Miccadei et al., 2002) implies that for most of the Pleistocene no fluvial connections existed at all between the mountain range interior and the Adriatic foreland area. Although fluvial incision in the foreland areas is clearly significant, these field observations suggest that most foreland draining streams have not been successful in enlarging their catchments into the faulted domain. Moreover, our modelling results support the idea that the Popoli gorge is more likely controlled by other local factors like pre-existing topography (Figure 10e), perhaps in combination with the collapse of underground drainage (Boni, 2000;Piacentini & Miccadei, 2014), rather than by efficient headward erosion from the coast.
A type of field observation that we also consider as indicative of overspill is what we call "interior gorges", i.e. deeply incised river valleys located in the interior part of the faulted domain that are not related to an erosional wave propagating upstream from the coast (Figures 7b and 12). Theoretically, this kind of gorge could be produced by headward erosion at a more local scale, e.g. by a first-order stream draining an individual hanging-wall basin margins. F I G U R E 1 2 Main features in extensional systems where drainage integration is either dominated by overspill mechanisms (a) and (c) or headward erosion from the coast (b) and (d) based on our model results. (e) Picture of the Sagittario gorge (see Figure 1b for locality), example of an "interior gorge" located in between two fault-bounded basins. It cannot be explained by headward erosion from the coast but might have formed as a consequence of basin overfilling. (f) Picture of the Popoli gorge (see Figure 1b for locality), located in between the Sulmona basin and the Adriatic foreland area. Our model results suggest that it cannot be explained by faulting, regional uplift and fluvial incision only, but it might be explained by pre-existing topography perhaps in combination with karst For example, Smith (2013) suggested that inter-basin headward erosion is favoured when a lower elevated faultbounded basin (the one containing the headward eroding stream) subsides at a faster rate than an adjacent higher elevated basin (the one becoming captured). Even where these tectonic conditions occur (e.g. Figure 6a), overspill can still dominate and can lead to local incision and gorge formation between adjacent basins (e.g. Figure 8a). Therefore, based on our model results, we expect interior gorges in the central Apennines to be mainly produced through overspill-driven drainage integration. One example of such an interior gorge is the San Venanzio gorge located between the Lower Aterno valley and Sulmona basin (Figure 1b). The fact that alluvial fan deposits at the outlet of the gorge interfinger with lacustrine deposits in the Sulmona basin (Cavinato & Miccadei, 2000) implies that this gorge was formed before the Sulmona basin was captured by headward eroding rivers that drain to the coast. Another example is the Sagittario river gorge, also located upstream of the Sulmona basin, but downstream of Lake Scanno (Figures 1b and 12e). The dimensions of this gorge suggest that it cannot be explained by an upstream propagating wave of erosion considering its position in the hangingwall of a large normal fault (Figure 1b) and the much more limited amount of incision in the downstream Sulmona basin. Therefore, based on our model results, we suggest both gorges most likely formed due to overspill from basins located directly upstream, leading to the formation of fluvial connections with the downstream located Sulmona basin followed by rapid local incision. Finally, our model results are consistent with an increasing number of studies that call into question headward erosion as being an important drainage integrating process (e.g. Bishop, 1995;Connell et al., 2005;Douglass et al., 2009;Heidarzadeh et al., 2017;Spencer & Pearthree, 2001). Theoretically, true headward erosion, i.e. the uphill lengthening of first-order streams, is expected to be a relatively inefficient process as discharge and consequently stream power are low close to the water divide (Bishop, 1995;Connell et al., 2005;Spencer & Pearthree, 2001). It is potentially relevant at the scale of gully systems for which headward erosion has been mainly described, where erosion is strongly associated with high runoff events and therefore relatively large amounts of water entering the gully heads due to sheet flow (e.g. Bishop, 1995;Bocco, 1991). We think that the relative inefficiency of headward erosion is clearly demonstrated by the drastic increase in incision rate that is generally observed directly following a drainage integration event (Figure 8a; see also e.g. Stokes et al., 2002). As long as the basin is still internally drained, erosion affecting the basin margins proceeds typically at a low rate and can only be explained by headward erosion by first-order streams. However, as soon as a fluvial connection becomes established, the spill-point area experiences an increase in discharge and slope causing a rapid increase in erosion rates (e.g. Garcia-Castellanos et al., 2003;Smith, 2013;Stokes et al., 2002). For instance, for the model river analysed in Figures 6 and 8, drainage integration results in a >6 times increase in incision rate (Figure 8a). Our conclusion is that headward erosion may have been invoked too often in regional or catchment-scale landscape evolution studies because drainage integration events are usually followed by intense erosion so that field evidence necessary for distinguishing between bottom-up and top-down integration mechanisms tends to become lost (Douglass et al., 2009).

| Impact of drainage integration on sediment dispersal
Our model results also have important implications for studying regional-scale sediment dispersal in the central Apennines and comparable settings. Top-down (basin overfilling and lake overspill) and bottom-up (headward erosion) mechanisms clearly produce different spatio-temporal patterns of sediment dispersal (Figure 12). In the case of headward erosion, a systematic pattern emerges which is a function of distance to the coast (Figure 12b,d). The more proximal to the coast the earlier lake sedimentation ceases, fluvial activity starts and incision of the basin fill commences. In the case of overspill, in contrast, the pattern is complex as local conditions become more important (Figure 12a,c): Lacustrine sedimentation ceases first in those basins that either; (i) have relatively high sediment input due to a large source area with high relief, (ii) have relatively low rates of basin subsidence or (iii) have a relatively low spill-point (e.g. due to pre-existing topography). However, in the case of overspill, sediment dispersal also strongly depends on the geometry of the drainage network and modifications to it over time. For instance, basins experience a significant increase in sediment supply when an upstream basin becomes externally drained and its sediment-fill becomes excavated. In other words, the top-down pattern of drainage integration is more difficult to predict because overfilling of a single basin is the integrated effect of all landscape developments occurring upstream and depends strongly on the regional-scale geometry and temporal evolution of the upstream drainage network. The temporal evolution of the drainage network, in turn, depends strongly on the growth of the extensional fault population (Cowie et al., 2006).
For the offshore area, our numerical experiments suggest a long-term increase in sediment supply due to the progressive increase in regional relief. This corresponds to field observations from the Adriatic where strong progradation started ca. 1.8 Ma (e.g. Artoni, 2013). On top of this gradual trend, however, our model predicts more step-wise increases in sediment supply due to drainage integration events. Considering the age of the Popoli gorge (ca. 0.4-0.35 Ma; Miccadei et al., 2002), which is the main sediment exit point on the Adriatic side of the central Apennines, we would expect a sudden increase in sediment to the Adriatic around this time. Based on the limited data available, no clear evidence exists that could confirm this but it is possible that the increase in sediment supply due to formation of the Popoli gorge is overprinted by effects due to a possible acceleration of regional uplift in the Adriatic foreland area around circa 0.7 Ma (e.g. Pizzi, 2003). Another implication of our model results is that the Adriatic mountain range flank has likely experienced more intense erosion than its Tyrrhenian counterpart (Figure 7d). This is because the dominant SW dip direction of the normal faults produces enhanced uplift and high relief driving erosion (Figure 9), even though the long-wavelength uplift is symmetric in our model.

| Transient landscape evolution as a function of regional uplift and normal faulting
Our model clearly demonstrates that landscape development in the central Apennines is transient even after 3 Myr. Even though the tectonic forcing is constant and climatic oscillations are not considered, we show that the landscape adapts continuously to modifications to the connectivity of the drainage network. This has an important implication because drainage integration represents a transient development that forms the background to other transient responses related to changes in allogenic forcing such as fault slip rate variations  or climate (Wegmann & Pazzaglia, 2009). Therefore, we consider drainage integration as an autogenic process inherent to many continental extensional systems and recommend it to be considered as an important element in future transient landscape studies in such settings. Furthermore, our model results suggest that in the central Apennines, drainage integration can be explained by the unique combination of normal faulting and differential regional uplift. On their own, these individual tectonic processes do not lead to drainage integration, either because no closed basins develop (in the case of regional uplift only) or because they do not become interconnected over time (in the case of normal faulting only). Besides fault development (controlling accommodation creation), we believe that the availability of sediment is a crucial factor in driving drainage integration and is potentially more important than external base-level fall. This means that in settings like the central Apennines where sediment originates only from the extensional domain itself, it is of key importance that there is enough relief to produce enough erosion and thereby sufficient sediment to fill the basins (favouring both basin overfilling and lake overspill). This relief can either be produced by active regional uplift or be inherited from pre-extensional times. Because of the high amplitude of regional uplift (up to ca. 1,000 m) across relative short distance (ca. 150 km), this requirement is fulfilled in the central Apennines, whereas the exact pattern of regional uplift is less relevant (see "Asymmetric uplift experiment" in Appendix S7).
In the Basin and Range Province, in contrast, the lack of sufficient relief and hence sediment supply may explain why drainage integration at a regional scale (including across-strike integration) is, in several areas, less advanced. The lack of relief can be overcome if an external sediment source is available (external to the extensional domain), e.g. the Gila river system (Arizona) that transports sediment from the southern edge of the Colorado Plateau to basins in the southern Basin and Range (Dickinson, 2015). Although the Gila river and its tributaries drain most of the fault-bounded basins in this region, there are also a few basins that remain internally drained (Dickinson, 2015). Importantly, these endorheic basins all have a distal position relative to the Colorado Plateau (the main sediment source), supporting the idea that sediment supply and overspill play a key role in controlling drainage integration.
Finally our study shows that drainage integration occurs even if both faulting and regional uplift accumulate uniformly over time. Although changes in tectonic deformation, for example due to fault propagation and interaction (Cowie et al., 2006), likely affected the evolution of the central Apennines river network, our model shows they are not needed to explain drainage integration. In other words, our simple model setup demonstrates that landscape evolution is highly dynamic even if the tectonic forcing is not. We expect changes in tectonic conditions over time to have made long-term drainage integration even more dynamic and to have enabled some basins to go through multiple cycles of internal and external drainage (e.g. Galli, Giaccio, & Messina, 2010; see also Appendix S8). Therefore, we expect the trend of lake disappearance seen in our model to be even more complex in reality.

| CONCLUSIONS
We have used a surface process model to investigate the phenomenon of drainage integration in the actively extending central Italian Apennines. By using a simple model setup that accounts for the main aspects of tectonic deformation in this area, i.e. regional uplift and normal faulting, we investigated the evolution of drainage integration, the roles of the main controlling mechanisms, and its impact on supporting her PhD research and Open Access publication, RLG acknowledges support from VISTA, and AHG and PAC acknowledge financial support from Statoil -University of Bergen Academia Agreement. LW was partly funded by NERC Urgency grant NE/P018858/1.