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

Anneleen Geurts, Patience Cowie, Guillaume Duclaux, Rob Gawthorpe, Ritske Huismans, Vivi K. Pedersen, Luke Wedmore

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Anneleen Geurts, Patience Cowie, Guillaume Duclaux, Rob Gawthorpe, Ritske Huismans, et al. Drainage integration and sediment dispersal in active continental rifts: A numerical modelling study of the central Italian Apennines. Basin Research, 2018, 30 (5), pp.965-989. 10.1111/bre.12289. hal-01983446

HAL Id: hal-01983446

https://hal.science/hal-01983446v1

Submitted on 19 Jun 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Geurts, A. H., Cowie, P. A., Duclaux, G., Gawthorpe, R. L., Huismans, R. S., Pedersen, V. K., & Wedmore, L. N. J. (2018). Drainage integration and sediment dispersal in active continental rifts: A numerical modelling study of the central Italian Apennines. *Basin Research, 30*(5), 965-989. https://doi.org/10.1111/bre.12289

Publisher's PDF, also known as Version of record

License (if available):

CC BY

Link to published version (if available):

10.1111/bre.12289

Link to publication record on the Bristol Research Portal

PDF-document

This is the final published version of the article (version of record). It first appeared online via Wiley at https://onlinelibrary.wiley.com/doi/abs/10.1111/bre.12289. Please refer to any applicable terms of use of the publisher.

**University of Bristol – Bristol Research Portal**

**General rights**

This document is made available in accordance with publisher policies. Please cite only the published version using the reference above. Full terms of use are available: http://www.bristol.ac.uk/red/research-policy/pure/user-guides/brp-terms/
# Drainage integration and sediment dispersal in active continental rifts: A numerical modelling study of the central Italian Apennines

Anneleen H. Geurts<fnref n="1" />

<figure><img src="image_3.png" /></figure>

Patience A. Cowie<fnref n="1" />

Guillaume Duclaux<fnref n="1" />

<figure><img src="image_4.png" /></figure>

Rob L. Gawthorpe<fnref n="1" />

<figure><img src="image_5.png" /></figure>

Ritske S. Huismans<fnref n="1" />

Vivi K. Pedersen<fnref n="1" />

Luke N. J. Wedmore<fnref n="2" />

## Correspondence

Anneleen H. Geurts, Department of Earth Science, University of Bergen, Bergen, Norway.

Email: anneleen.geurts@uib.no

## Present addresses

Guillaume Duclaux, CNRS, OCA, IRD, Université Côte d'Azur, Géoazur, Valbonne, France.

Luke N. J. Wedmore, School of Earth Sciences, University of Bristol, Bristol, UK.

## Funding information

Universitetet i Bergen; VISTA; NERC
Urgency, Grant/Award Number: NE/
P018858/1

## Abstract

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.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2018 The Authors. Basin Research © 2018 John Wiley & Sons Ltd, European Association of Geoscientists & Engineers and International Association of Sedimentologists

<footnote n="1">¹Department of Earth Science, University of Bergen, Bergen, Norway</footnote>

<footnote n="2">²Institute of Geophysics and Tectonics, University of Leeds, Leeds, UK</footnote>
# 1 | 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

## Highlights

* Landscape evolution modelling study applied to the actively extending central Italian Apennines

* Transient landscape development even if tectonic and climatic forcing in the model are uniform over time

* Faulting and regional uplift explain drainage integration, i.e. the transition from internal to external drainage

* Drainage integration is driven by “top-down” overspill mechanisms rather than “bottom-up” headward erosion

* Drainage integration causes abrupt local changes in erosion/deposition and delayed sediment export offshore

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
<figure><img src="image_6.png" /><figcaption>**FIGURE 1** (a) Location map of the study area in the central Apennines, (b) the study area and model domain itself (DEM from Tarquini et al., 2007) and (c) a simplified geological map of the area showing the main lithological units (modified from Whittaker et al., 2008). On top of the topographic map in (b) we show active (solid pink lines) and inactive (dashed pink lines) normal faults (principally after Roberts & Michetti, 2004). We also show sea-level markers (*M1* = Early Pleistocene shoreline, *M2* = Last-interglacial shoreline, *M3* = Sicilian shoreface deposits, and *M4* = Last-interglacial floodplain, see also Appendix S2), fault sites (*FiF* = Fiamignano fault, *FuF* = Fucino fault, *BaF* = Barete fault and *SuF* = Sulmona fault; see also Figure 2b) and paleosurfaces which were used for estimating long-wavelength uplift (see also Appendix S3). Other abbreviations used on the map are: *TER* = Terni basin, *SAG* = Sagittario gorge, *VEN* = San Venanzio gorge and *POP* = Popoli gorge</figcaption></figure>

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.

## 2 | GEOLOGICAL SETTING

The central part of the Italian Apennines is the highest (up to ca. 2,900 m) and widest part of the Apennines mountain belt. After cessation of thrusting during the Pliocene (Patacca, Sartori, & Scandone, 1990), its Quaternary landscape evolution (ca. the last 3 million years) has been dominated by a combination of regional uplift and southwest-

northeast extension localised on dominantly southwest dipping normal faults (e.g. Roberts & Michetti, 2004; Figure 1b). Brackish marsh deposits at the base of some of the extensional basins suggest the area was close to sea-level when extension and uplift commenced (Gliozzi & Mazzini, 1998). Regional uplift has produced a topographic bulge as the mountain belt interior has undergone large (>800 m) differential uplift relative to the coastlines (Ascione, Cinque, Miccadei, Villani, & Berti, 2008; Centamore & Nisio, 2003; D'Agostino et al., 2001; D'Anastasio et al., 2006; Mancini, D'Anastasio, Barbieri, & de Martini, 2007; Pizzi, 2003; Serpelloni, Faccenna, Spada, Dong, & Williams, 2013). Most of the extensional deformation has occurred along the crest of this topographic bulge and is accommodated by a wide (>60 km) array of normal faults (Figure 1b; D'Agostino et al., 2001; Roberts & Michetti, 2004; Faure Walker, Roberts, Sammonds, & Cowie, 2010). Over time these faults have generated large footwall uplifts mainly consisting of Mesozoic limestone and have trapped thick sequences of continental deposits in their hanging-wall basins (Figure 1c; Cavinato, 1993; Cavinato,
Cosentino, de Rita, Funiciello, & Parotto, 1994; Cavinato & DeCelles, 1999; Miccadei et al., 2002; Cavinato, Carusi, Dall'Asta, Miccadei, & Piacentini, 2002). Total throw estimates along the (up to 40 km long) faults vary across the area but tend to be greatest (up to 2,200 m) across the more centrally located, higher elevation, faults (Cowie & Roberts, 2001; Roberts & Michetti, 2004).

The elevated topography in the central part of the Apennines cannot be explained by crustal or lithospheric instability (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., 1994, 2002; Miccadei 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.

## 3 | 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 × 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

<table><caption>**TABLE 1** Overview of parameter values used in the surface process model CASCADE</caption><thead><tr><td>Parameter</td><td>Description</td><td>Values</td><td>Units</td></tr></thead><tbody><tr><td>dx, dy</td><td>Grid resolution</td><td>1,000</td><td>m</td></tr><tr><td>dt</td><td>Calculation timestep (adjusted dynamically)</td><td>ca. 100</td><td>year</td></tr><tr><td>endtime</td><td>Length of model run</td><td>3 × 10<sup>6</sup></td><td>year</td></tr><tr><td>v</td><td>Effective precipitation rate</td><td>1</td><td>m/year</td></tr><tr><td>K<sub>f</sub></td><td>Dimensionless fluvial transport parameter</td><td>0.08–0.12</td><td>–</td></tr><tr><td>L<sub>f</sub></td><td>Fluvial erosion length scale</td><td>30–70 × 10<sup>3</sup></td><td>m</td></tr><tr><td>φ</td><td>Scaling exponent for channel width W = c × Q<sup>φ</sup></td><td>0.5</td><td>–</td></tr><tr><td>c</td><td>Scaling factor for channel width W = c × Q<sup>c</sup></td><td>1</td><td>–</td></tr></tbody></table>
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.

## 3.1 | 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 Roberts & Michetti, 2004 and Wedmore 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 (γ 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 = \gamma \times 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 γ has the greatest impact on the vertical surface displacement field. Our intermediate value for γ (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 γ = 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; Whittaker, Attal, Cowie, Tucker, & Roberts, 2008), 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 well-dated 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
<figure><img src="image_7.png" /><figcaption>FIGURE 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 *γ* = 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)</figcaption></figure>

assuming typical long-term ratios of footwall uplift to
hanging-wall subsidence and by assuming that the land sur-
face was close to sea-level before regional uplift started. A
detailed description of our method, the data, and our regio-
nal 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 inner-
most 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
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).

## 3.2 | 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):

$$ \frac{dh}{dt} = \frac{1}{W \cdot L_f} (Q_c - Q_s) \quad (1) $$

where $\frac{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$:

$$ Q_c = K_f \cdot Q_w \cdot S \quad (2) $$

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:

$$ Q_s = \int_{0}^{A} \frac{dh}{dt} da \quad (3) $$

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 $\frac{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 = \sqrt{Q_w}$). 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).
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>
# 4 | MODEL RESULTS

## 4.1 | 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 steady-state 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 (Whittaker et al., 2008).

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 pre-existing 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 run-time (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 fault-related topography on the Adriatic side (Figure 2a). However, as discussed below in “Regional-scale sediment

<figure><img src="image_16.png" /></figure>

**FIGURE 4** Topographic profile (top panel) and 100-km wide topographic swath profile (bottom panel) across the final (3 Myr) topography of the reference model (transect and swath shown in 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)

**FIGURE 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. (a–d) Topographic evolution. (e–h) Development of the drainage network. The maps show only lakes >10 km², associated with the fault-bounded basins. By fluvial “exit points” we mean localities on the edges of the faulted domain where the intermountain drainage network becomes integrated with streams draining towards the coast. At these localities, water and sediment exits the mountainous area affected by normal faulting. For model sensitivity tests, see Appendix S4
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. Whitaker 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.

## 4.2 | 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

<figure><img src="image_17.png" /></figure>

**FIGURE 5** Time evolution plots showing different aspects of modelled landscape evolution. (a) Maximum and mean elevation calculated for the total model domain. (b) Size of the total area that is internally drained (as a percentage of the total model domain), as indicated by the grey area in Figure 3e–h. (c) Total surface area occupied by lakes (as a percentage of the total model domain). (d) Total volume of all lakes together. The volume of each individual lake is determined by calculating the volume of water that is needed to fill a topographic depression up to its spill-point. (e) Mean erosion rate, based on the total area experiencing erosion (excluding depocentres). (f) Sediment flux crossing the coastal (left and right) boundaries of the model

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
<figure><img src="image_18.png" /><figcaption>**FIGURE 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</figcaption></figure>

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).

## 4.3 | 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
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

**FIGURE 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

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.

## 4.4 | Regional-scale sediment dispersal

Erosion rate maps (Figure 7a–c) show the general pattern of erosion and deposition in our reference model: Sediment
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

<figure><img src="image_23.png" /></figure>

**FIGURE 8** (a) Erosion/Deposition rate along the stream shown in Figure 6 (crossing two distinct basins) for the time period in between 1.7 and 2.7 Myr. This period encloses the event of basin overfilling for Basin II (Figure 6) for which its spill-point is marked by means of a yellow star (corresponding to yellow stars in Figure 7). Black lines show erosion rates before overfilling of Basin II takes place, and red lines do the same for after basin overfilling. (b) Longitudinal profile along the same stream as analysed in (a) and Figure 6 at different model time steps, showing the way in which the stream profile adapts to base-level change following basin overfilling. Although our standard model time is 3 Myr, we also show longitudinal profiles developed after 3.5 and 4 Myr when steady state is approximately reached

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
<figure><img src="image_24.png" /><figcaption>**FIGURE 9** Cartoon explaining the structural asymmetry between the Tyrrhenian (a) and Adriatic (b) sides of the central Apennines (due to the dominant southwest dip direction of the normal faults) and its implications for the rate of drainage integration and erosion on the mountain flanks</figcaption></figure>

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).

**5 | 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.

**5.1 | 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

**FIGURE 10** Comparison of the topography and stream network of the central Apennines with our reference model. (a) Ten-metre DEM of the central Apennines (Tarquini et al., 2007) interpolated at 1-km resolution, similar to the model resolution, together with normal faults, the stream network, the central water divide (for legend see also Figure 3) and the internally drained Fucino basin ("FUC"). (b) Stream network and catchment geometry of the central Apennines, derived from the DEM shown in (a). The yellow dot shows the Popoli gorge that is the main locality where surface water exits the Adriatic side of the faulted domain. 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,b) from the DEM (Figure 10a) and has been smoothed by means of a Gaussian kernel (sigma = 4 km) in order to reveal the main topographic features. Main thrust faults (modified from Miccadei, Piacentini, & Buccolini, 2017) are shown on top. The morphological pattern shown by the residual most likely reflects pre-existing, thrust-related, topography. (f) Hypometric distributions for both the DEM and the final (3 Myr) topography of our reference model shown in (a) and (c) respectively
& 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

Reality

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

Model

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

(e) Residual as indicator for the pattern of pre-existing topography

<figure><img src="image_29.png" /></figure>

(f) Reality

<figure><img src="image_30.png" /></figure>
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, litho-
logical 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 sug-
gests 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 wind-
gap in between Maiella and Sulmona (Figure 11).
Although this valley and windgap may have other explana-
tions, our results clearly demonstrate that these two locali-
ties 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 addi-
tion, there is an intermediate-scale morphology in the cen-
tral 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 resid-
ual topography (present-day topography minus our tectonic
uplift function) shown in Figure 10e. Because the land-
scape 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 com-
pression 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, Cha-
mot-Rooke, Funiciello, Jolivet, & Speranza, 1998; Scis-
ciani, Tavarnelli, & Calamita, 2002). However, here we
show that inherited topography is not a necessary ingredi-
ent 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 drai-
nage 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 Apen-
nines corresponds to the elevation of the internally
drained Fucino basin, at circa 650 m. The Fucino basin
remains internally drained today because there is insuffi-
cient sediment supply compared to the high rate of
accommodation creation (see “Overspill vs. headward ero-
sion 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 pri-
mary 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 (Fig-
ure 9; see “Regional-scale sediment dispersal” in the
Model results section).

## 5.2 | Overspill vs. headward erosion from the coast

Our model results demonstrate that an important mecha-
nism 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 (Fig-
ure 6a). We note, however, that our endorheic model set-
up (see Methodology section) does not allow us to distin-
guish 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,
<figure><img src="image_31.png" /><figcaption>**FIGURE 11** 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</figcaption></figure>

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).
(a) OVERSPILL

(basin overfilling / lake overspill)

Formation of 'interior gorges' after
fluvial connections between basins
are established

<figure><img src="image_32.png" /></figure>

No clear spatial pattern of
drainage integration

Longevity of closed basins is
controlled by local conditions
affecting the balance
between sediment supply
and accommodation creation

Headward erosion occurs but
is not the dominant process
causing drainage integration

"TOP-DOWN"

(b) HEADWARD EROSION

from the coast

Survival of closed basins most
likely furthest away from
regional baselevel (coast)

<figure><img src="image_33.png" /></figure>

Clear spatial-temporal pattern
as a function of proximity to the coast

Gorges clearly associated
with an upstream
migrating wave of erosion

"BOTTOM-UP"

Erosional wave
propagating upstream
from the coast

(c) IF overspill controls drainage integration

TIME

<figure><img src="image_34.png" /></figure>

s.l. = sea level

fluvial connection with the coast

local fluvial connections between
adjacent basins ('interior gorges')

Large endorheic area with many closed basins and
lakes and local sediment redistribution

<figure><img src="image_35.png" /></figure>

s.l.

Endorheic area shrinks because outermost basins are
captured by regional drainage network

<figure><img src="image_36.png" /></figure>

Closed basins survive
longer closer to the central
drainage divide

<figure><img src="image_37.png" /></figure>

s.l.

internal drainage and
sediment dispersal
external drainage and
sediment dispersal

(e)

<figure><img src="image_38.png" /></figure>

Sagittario gorge

N

(f)

<figure><img src="image_39.png" /></figure>

Popoli gorge

N ←
**FIGURE 12** 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

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 “survival-concept” as it shows a transition from overfilled to underfilled conditions over time (Cavinato et al., 2002) that can be explained by a ×3 to ×5 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.
For example, Smith (2013) suggested that inter-basin headward erosion is favoured when a lower elevated fault-bounded 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 hanging-wall 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).

## 5.3 | 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 prograda-tion 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.

**5.4 | 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 (Whittaker et al., 2008) 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.

**6 | 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
regional-scale sediment dispersal. Based on our modelling,
our specific conclusions are:

(1) Both regional uplift and extensional deformation are important for long-term landscape evolution in the central Apennines. Together they reproduce the main landscape features and essential transient aspects of its evolution, in particular, the cessation of lake sedimentation and drainage integration.

(2) Basin overfilling, and hence overspill and drainage integration occur in our model because of the increasing relief at both fault-block and regional scales that generates more erosion and in turn more sediment supply to basins over time. Even for the case of constant fault slip rates, this causes basins to become progressively overfilled and eventually a through-going river system to develop (e.g. Figure 6).

(3) Our model suggests overspill (basin overfilling and lake overspill) rather than headward erosion from the coast to be the dominant fluvial mechanism driving drainage integration in the central Apennines, i.e. "top-down" rather than "bottom-up" integration (Figure 12). These results are consistent with field observations from the central Apennines, in particular the formation of "interior gorges" (Figure 12), and with an increasing number of other studies that call into question headward erosion as being an important process for regional-scale drainage integration.

(4) Overspill depends largely on the balance between sediment supply and accommodation creation in individual fault-bounded basins. Because both of them depend on many local factors and developments along the drainage network upstream, basin overfilling does not produce a clear spatio-temporal pattern of drainage integration—unlike in the case of headward erosion where it depends primarily on distance to the coast in the case of Italy (Figure 12).

(5) We show that landscape evolution can be highly dynamic even if the tectonic forcing and climate are uniform over time. This is because drainage integration causes the landscape to adapt continuously to modifications in connectivity of the drainage network. Other processes like fault interaction are likely to make drainage integration even more dynamic in reality (e.g. enabling some basins to go through multiple cycles of internal and external drainage), although it is not needed in order to explain the phenomenon itself.

(6) Over long timescales of millions of years, drainage integration produces a delayed export of sediment out of the area affected by normal faulting and a step-wise increase in sediment supply offshore. At a local scale, it leads to abrupt changes in erosion/deposition patterns, marked variation in sediment supply to basins

and hence sedimentary environment (lacustrine vs. fluvial), and strong incision following drainage integration events (e.g. Figures 7 and 8).

(7) According to our model results, the dominant SW dip of the normal faults in the central Apennines favours overspill and therefore drainage integration within the Tyrrhenian part of the faulted domain, compared to its Adriatic counterpart. The Popoli gorge is an exception that is probably explained by local factors. Moreover, this structural asymmetry generates more intense erosion on the Adriatic flanks than on the Tyrrhenian flanks of the mountain range (e.g. Figure 9).

(8) We suggest that the most important factor for drainage integration to occur in continental extensional systems is the availability of sufficient sediment relative to the accommodation being created through normal faulting (more important than proximity to the coast, or other external base-levels). The important role that normal faulting plays both through the uplift of source areas and the accommodation creation in hanging-wall basins leads to the conclusion that better understanding of the underlying geodynamic mechanism(s) for fault growth is vital. In the case of the central Apennines, this is likely related to either flow or buoyancy variations in the uppermost mantle and associated with surface uplift at a regional scale (e.g. Cowie et al., 2013; Faccenna et al., 2014; Faure Walker et al., 2012). Sufficient sediment, on the other hand, can alternatively be provided through (pre-extensional) inherited relief, strong regional uplift (in case of the central Apennines) or an external sediment source (e.g. high topography adjacent to the continental rift).

(9) Finally, our results reveal abrupt and complex shifts in patterns of erosion/deposition at the fault-block scale, suggesting that feedbacks between surface processes and fault development may be enhanced, potentially contributing to temporal variations in fault slip rates and/or fault activity over time (e.g. Maniatis et al., 2009).

## ACKNOWLEDGEMENTS

We thank Jean Braun for providing the numerical code and for his help during the early stage of setting up the model experiments. We thank Gerald Roberts and Alessandro Michetti for providing the fault data and sharing their ideas at the time of constraining our model set-up. We also acknowledge Alberto Pizzi, Elsa Gliozzi and Tommaso Piacentini for providing more details on their data or publications. This study greatly benefitted from stimulating discussions with Alex Whittaker, Stefano Pucci and Fabio Villani, and from constructive comments from Peter van der Beek and an anonymous reviewer on an earlier draft of our manuscript. AHG acknowledges Bergen University for
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.

ORCID

Anneleen H. Geurts ID http://orcid.org/0000-0002-2973-5946

Rob L. Gawthorpe ID http://orcid.org/0000-0002-4352-6366

REFERENCES

Artoni, A. (2013). The Pliocene-Pleistocene stratigraphic and tectonic evolution of the Central sector of the Western Periadriatic Basin of Italy. *Marine and Petroleum Geology, **42*, 82–106.

Ascione, A., Cinque, A., Miccadei, E., Villani, F., & Berti, C. (2008). The Plio-Quaternary uplift of the Apennine chain: New data from the analysis of topography and river valleys in Central Italy. *Geomorphology, **102*, 105–118.

Bartolini, C., D'Agostino, N., & Dramis, F. (2003). Topography, exhumation, and drainage network evolution of the Apennines. *Episodes, **26*, 212–216.

Bishop, P. (1995). Drainage rearrangement by river capture, behead-
ing and diversion. *Progress in Physical Geography*, *19*, 449–473.

Bocco, G. (1991). Gully erosion: Processes and models. *Progress in Physical Geography, **15*, 392–406.

Boni, C. (2000). Karst aquifers of the Central Apennines. *Hydrogéologie, **4*, 49–62.

Braun, J., & Sambridge, M. (1997). Modelling landscape evolution on geological timescales: A new method based on irregular spatial discretization. *Basin Research, **9***, 27–52.

Buiter, S. J. H., Huismans, R. S., & Beaumont, C. (2008). Dissipation analysis as a guide to mode selection during crustal extension and implications for the styles of sedimentary basins. *Journal of Geophysical Research, **113*, B06406.

Cavinato, G. P. (1993). Recent tectonic evolution of the Quaternary deposits of the Rieti Basin (Central Apennines, Italy): Southern part. *Geologica Romana, **29*, 411–434.

Cavinato, G. P., Carusi, C., Dall'Asta, M., Miccadei, E., & Piacentini, T. (2002). Sedimentary and tectonic evolution of Plio-Pleistocene alluvial and lacustrine deposits of Fucino Basin (central Italy). *Sedimentary Geology, **148*, 29–59.

Cavinato, G. P., Cosentino, D., de Rita, D., Funiciello, R., & Par-
otto, M. (1994). Tectonic-sedimentary evolution of intrapenninic
basins and correlation with the volcano-tectonic activity in Cen-
tral Italy. *Memorie Descrittive Della Carta Geologica d'Italia*,
*49*, 63–76.

Cavinato, G. P., & DeCelles, P. G. (1999). Extensional basins in tectonically bi-modal central Apennines fold-thrust belt, Italy: Response to corner flow above a subducting slab in retrograde motion. *Geology, **27***, 955–958.

Cavinato, G. P., & Miccadei, E. (2000). Pleistocene carbonate lacustrine deposits: Sulmona Basin (Central Apennines, Italy). In E. H. Gierlowski-Kordesch & K. R. Kelts (Eds.), Lake basins through space and time (pp. 517–526). AAPG Studies in Geology, 46.

Centamore, E., & Nisio, S. (2003). Ejects of uplift and tilting in the Central-Northern Apennines, Italy. *Quaternary International*, **101-102**, 93–101.

Connell, S. D., Hawley, J. W., & Love, D. W. (2005). Late Cenozoic drainage development in the southeastern Basin and Range of New Mexico, southeasternmost Arizona, and western Texas. In S. G. Lucas, G. S. Morgan & K. Zeigler (Eds.), *New Mexico's ice ages* (Vol. 28, pp. 125-150). Albuquerque, NM: New Mexico Museum of Natural History and Science.

Cowie, P. A. (1998). Normal fault growth in three-dimensions in con-
tinental and oceanic crust. In W. Roger Buck, P. T. Delaney, J. A.
Karson & Y. Lagabrielle (Eds.), *Faulting and magmatism at mid-
ocean ridges* (Vol. 106, pp. 325–348). Washington, DC: Ameri-
can Geophysical Union Monograph.

Cowie, P. A., Attal, M., Tucker, G. E., Whittaker, A. C., Naylor, M., Ganas, A., & Roberts, G. P. (2006). Investigating the surface process response to fault interaction and linkage using a numerical modelling approach. *Basin Research, **18***, 231–266.

Cowie, P. A., & Roberts, G. P. (2001). Constraining slip rates and spacings for active normal faults. *Journal of Structural Geology*, **23**, 1901–1915.

Cowie, P. A., & Scholz, C. H. (1992). Displacement-length scaling relationship for faults: Data synthesis and discussion. *Journal of Structural Geology, **14*, 1149–1156.

Cowie, P. A., Scholz, C. H., Roberts, G. P., Faure Walker, J. P., & Steer, P. (2013). Viscous roots of active seismogenic faults revealed by geologic slip rate variations. *Nature Geoscience, **6***, 1036–1040.

Cowie, P. A., Whittaker, A. C., Attal, M., Roberts, G. P., Tucker, G. E., & Ganas, A. (2008). New constraints on sediment-flux-dependent river incision: Implications for extracting tectonic signals from river profiles. *Geology, **36*, 535–538.

D'Agostino, N., Chamot-Rooke, N., Funiciello, R., Jolivet, L., & Speranza, F. (1998). The role of pre-existing thrust faults and topography on the styles of extension in the Gran Sasso range (central Italy). *Tectonophysics, **292***, 229–254.

D'Agostino, N., Jackson, J. A., Dramis, F., & Funiciello, R. (2001). Interactions between mantle upwelling, drainage evolution and active normal faulting: An example from the central Apennines (Italy). *Geophysical Journal International*, **147**, 475–497.

D'Agostino, N., & McKenzie, D. (1999). Convective support of long-wavelength topography in the central Apennines (Italy). *Terra Nova, **11***, 234–238.

D'Alessandro, L., Miccadei, E., & Piacentini, T. (2003). Morphostructural elements of central-eastern Abruzzi: Contributions to the study of the role of tectonics on the morphogenesis of the Apennine chain. *Quaternary International*, **101–102**, 115–124.

D'Alessandro, L., Miccadei, E., & Piacentini, T. (2008). Morphotec-
tonic study of the lower Sangro River valley (Abruzzi, Central
Italy). *Geomorphology, **102*, 145–158.

D'Anastasio, E., de Martini, P. M., Selvaggi, G., Pantosti, D., Mar-
chioni, A., & Maseroli, R. (2006). Short-term vertical velocity
field in the Apennines (Italy) revealed by geodetic leveling data.
*Tectonophysics*, **418**, 219–234.

Dickinson, W. R. (2015). Integration of the Gila River drainage system through the Basin and Range province of southern Arizona and southwestern New Mexico (USA). *Geomorphology*, **236**, 1–24.

Douglass, J., Meek, N., Dorn, N. I., & Schmeekle, M. W. (2009). A criteria-based methodology for determining the mechanism of
transverse drainage development, with application to the south-
western United States. *Geological Society of America Bulletin*,
*121*, 586-598.

Douglass, J., & Schmeekle, M. W. (2007). Analogue modeling of transverse drainage mechanisms. *Geomorphology*, *84*, 22–43.

Duffy, O. B., Brocklehurst, S. H., Gawthorpe, R. L., Leeder, M. R., & Finch, E. (2015). Controls on landscape and drainage evolution in regions of distributed normal faulting: Perachora Peninsula, Corinth Rift, Central Greece. *Basin Research*, *27*, 473–494.

Faccenna, C., Becker, T. W., Miller, M. S., Serpelloni, E., & Willett, S. D. (2014). Isostasy, dynamic topography, and the elevation of the Apennines of Italy. *Earth and Planetary Science Letters*, *407*, 163–174.

Faure Walker, J. P., Roberts, G. P., Cowie, P. A., Papanikolaou, I., Michetti, A. M., Sammonds, P., ... Phillips, R. J. (2012). Relationship between topography and strain rate in the actively extending Italian Apennines. *Earth and Planetary Science Letters*, 325–326, 76–84.

Faure Walker, J. P., Roberts, G. P., Sammonds, P. R., & Cowie, P. A. (2010). Comparison of earthquake strains over 10<sub>2</sub> and 10<sub>4</sub> year timescales: Insights into variability in the seismic cycle in the central Apennines, Italy. *Journal of Geophysical Research*, 115, B10418. https://doi.org/10.1029/2009JB006462

Galli, P., Giaccio, B., & Messina, P. (2010). The 2009 central Italy earthquake seen through 0.5 Myr-long tectonic history of the L'Aquila faults system. *Quaternary Science Reviews*, *29*, 3768–3789.

Garcia-Castellanos, D., Verges, J., Gaspar-Escribano, J., & Cloetingh, S. (2003). Interplay between tectonics, climate, and fluvial transport during the Cenozoic evolution of the Ebro basin (NE Iberia). *Journal of Geophysical Research*, *108*, B7.

Gawthorpe, R. L., & Leeder, M. R. (2000). Tectono-sedimentary evolution of active extensional basins. *Basin Research*, *12*, 195–218.

Gliozzi, E., & Mazzini, I. (1998). Palaeoenvironmental analysis of Early Pleistocene brackish marshes in the Rieti and Tiberino intrapenninic basins (Latium and Umbria, Italy) using ostracods (Crustacea). *Palaeogeography, Palaeoclimatology, Palaeoecology*, *140*, 325–333.

Heidarzadeh, G., Ballato, P., Hassanzadeh, J., Ghassemi, M. R., & Strecker, M. R. (2017). Lake overspill and onset of fluvial incision in the Iranian Plateau: Insights from the Mianeh Basin. *Earth and Planetary Science Letters*, *469*, 135–147.

House, P. K., Pearthree, P. A., & Perkins, M. E. (2008). Stratigraphic evidence for the role of lake spillover in the inception of the lower Colorado River in southern Nevada and western Arizona. In M. C. Reheis, R. Hershler & D. M. Miller (Eds.), Late Cenozoic drainage history of the southwestern great basin and lower Colorado River region: Geological and biotic perspectives (pp. 335–353). Geological Society of America Special Paper, 439.

Jackson, J., & Leeder, M. (1994). Drainage systems and the develop-
ment of normal faults – An example from pleasant valley,
Nevada. *Journal of Structural Geology*, *16*, 1041–1059.

Kooi, H., & Beaumont, C. (1996). Large-scale geomorphology: Clas-
sical concepts reconciled and integrated with contemporary ideas
via a surface processes model. *Journal of Geophysical Research,*
*101*, 3361–3386.

Leeder, M. R., & Jackson, J. A. (1993). The Interaction between nor-
mal faulting and drainage in active extensional basins, with

examples from the Western United States and Central Greece.
*Basin Research*, 5, 79–102.

Lin, J., & Stein, R. S. (2004). Stress triggering in thrust and subduc-
tion earthquakes, and stress interaction between the southern San
Andreas and nearby thrust and strike-slip faults. *Journal of Geo-
physical Research*, 109, B02303.

Ludovisi, A., Gaino, E., Bellezza, M., & Casadei, S. (2013). Impact of climate change on the hydrology of shallow Lake Trasimeno (Umbria, Italy): History, forecasting and management. *Aquatic Ecosystem Health & Management*, 16(2), 190–197.

Mancini, M., D'Anastasio, E., Barbieri, M., & de Martini, P. M. (2007). Geomorphological, paleontological and ⁸⁷Sr/⁸⁶Sr isotope analyses of early Pleistocene paleoshorelines to define the uplift of Central Apennines (Italy). *Quaternary Research*, 67, 487–501.

Maniatis, G., Kurfeb, D., Hampel, A., & Heidbach, O. (2009). Slip acceleration on normal faults due to erosion and sedimentation – Results from a new three-dimensional numerical model coupling tectonics and landscape evolution. *Earth and Planetary Science Letters*, 284, 570–582.

Menges, C. M. (2008). Multistage late Cenozoic evolution of the Amargosa River drainage, southwestern Nevada and eastern California. In M. C. Reheis, R. Hershler & D. M. Miller (Eds.), Late Cenozoic drainage history of the southwestern Great Basin and lower Colorado River region: Geological and biotic perspectives (pp. 39–90). Geological Society of America Special Paper, 439.

Miccadei, E., Piacentini, T., & Barberi, R. (2002). Uplift and local tectonic subsidence in the evolution of intramontane basins: The example of the Sulmona basin (central Apennines, Italy). *Studi Geologici Camerti*, 2002, 119-133.

Miccadei, E., Piacentini, T., & Buccolini, M. (2017). Long-term geomorphological evolution in the Abruzzo area, Central Italy: Twenty years of research. *Geologica Carpathica*, 68(1), 19-28.

Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space. *Bulletin of the Seismological Society of America*, 82(2), 1018–1040.

Patacca, E., Sartori, R., & Scandone, P. (1990). Tyrrhenian Basin and Apenninic Arcs: Kinematic relations since late Tortonian times. *Memorie Della Societa' Geologica Italiana*, 45, 425-451.

Piacentini, T., & Miccadei, E. (2014). The role of drainage systems and intermontane basins in the Quaternary landscape of the Central Apennines chain (Italy). *Rendiconti Lincei*, 25, S139-S150.

Pizzi, A. (2003). Plio-Quaternary uplift rates in the outer zone of the central Apennines fold-and-thrust belt, Italy. *Quaternary International*, 101–102, 229–237.

Pucci, S., Villani, F., Civico, R., Pantosti, D., del Carlo, P., Smedile, A., ... Gueli, A. (2014). Quaternary geology of the Middle Aterno Valley, 2009 L'Aquila earthquake area (Abruzzi Apennines, Italy). *Journal of Maps*, *11*, 689–697.

Roberts, G. P., & Michetti, A. M. (2004). Spatial and temporal varia-
tions in growth rates along active normal fault systems: An exam-
ple from Lazio-Abruzzo, central Italy. *Journal of Structural
Geology*, *26*, 339–376.

Scisciani, V., Tavarnelli, E., & Calamita, F. (2002). The interaction of extensional and contractional deformations in the outer zones of the central Apennines, Italy. *Journal of Structural Geology*, 24, 1647–1658.

Serpelloni, E., Faccenna, C., Spada, G., Dong, D., & Williams, S. D. P. (2013). Vertical GPS ground motion rates in the Euro-
<figure><img src="image_40.png" /></figure>

Mediterranean region: New evidence of velocity gradients at dif-
ferent spatial scales along the Nubia-Eurasia plate boundary. *Jour-
nal of Geophysical Research: Solid Earth*, **118**, 1–22.

Smith, J. (2013). *Source-to-sink analysis of rift basin tectonics and sedimentation*. PhD thesis, Manchester University, UK, p. 193.

Spencer, J. E., & Pearthree, P. A. (2001). Headward erosion versus closed-basin spillover as alternative causes of Neogene capture of the ancestral Colorado River by the Gulf of California. *The Colorado River: Origin and Evolution: Grand Canyon, Arizona, Grand Canyon Association Monograph*, *12*, 215–219.

Stokes, M., Mather, A. E., & Harvey, A. M. (2002). Quantification of river-capture-induced base-level changes and landscape development, Sorbas Basin, SE Spain. *Geological Society, London, Special Publications*, *191*, 23–35.

Tarquini, S., Isola, I., Favalli, M., Mazzarini, F., Bisson, M., Pareschi, M. T., & Boschi, E. (2007). TINITALY/01: a new triangular irregular network of Italy. *Annals of Geophysics*, *50*, 407–425.

Toda, S., Stein, R. S., Richards-Dinger, K., & Bozkurt, S. (2005). Forecasting the evolution of seismicity in southern California: Animations built on earthquake stress transfer. *Journal of Geophysical Research*, *110*, B05S16.

Twidale, C. R. (2004). River patterns and their meaning. *Earth Science Reviews*, *67*, 159–218.

van der Beek, P., & Bishop, P. (2003). Cenozoic river profile devel-
opment in the Upper Lachlan catchment (SE Australia) as a test
of quantitative fluvial incision models. *Journal of Geophysical
Research*, *108*, 2309.

Wedmore, L. N. J., Faure Walker, J. P., Roberts, G. P., Sammonds, P. R., McCaffrey, K. J. W., & Cowie, P. A. (2017). A 667 year record of coseismic and interseismic Coulomb stress changes in central Italy reveals the role of fault interaction in controlling irregular earthquake recurrence intervals. *Journal of Geophysical Research: Solid Earth*, *122*, 1–21.

Wegmann, K. W., & Pazzaglia, F. J. (2009). Late Quaternary fluvial terraces of the Romagna and Marche Apennines, Italy: Climatic,

lithologic, and tectonic controls on terrace genesis in an active
orogeny. *Quaternary Science Reviews*, *28*, 137–165.

Whittaker, A. C. (2012). How do landscapes record tectonics and climate? *Lithosphere*, *4*(2), 160–164.

Whittaker, A. C., Attal, M., & Allen, P. A. (2010). Characterising the origin, nature and fate of sediment exported from catchments perturbed by active tectonics. *Basin Research*, *22*, 809–828.

Whittaker, A. C., Attal, M., Cowie, P. A., Tucker, G. E., & Roberts, G. P. (2008). Decoding temporal and spatial patterns of fault uplift using transient river long-profiles. *Geomorphology*, *100*, 506–526.

Zanchetta, G., Borghini, A., Fallick, A. E., Bonadonna, F. P., & Leone, G. (2007). Late Quaternary palaeohydrology of Lake Perugusa (Sicily, southern Italy) as inferred by stable isotopes of lacustrine carbonates. *Journal of Paleolimnology*, 38, 227–239.

SUPPORTING INFORMATION

Additional supporting information may be found online in
the Supporting Information section at the end of the article.

**How to cite this article:** Geurts AH, Cowie PA, Duclaux G, et al. Drainage integration and sediment dispersal in active continental rifts: A numerical modelling study of the central Italian Apennines. *Basin Res.* 2018;30:965–989. https://doi.org/10.1111/bre.12289