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

# Autogenic versus allogenic controls on the evolution of a coupled fluvial megafan–mountainous catchment system: numerical modelling and comparison with the Lannemezan megafan system (northern Pyrenees, France)

Margaux Mouchené<fnref n="1" />,ª, Peter van der Beek<fnref n="1" />, Sébastien Carretier<fnref n="2" />,<fnref n="3" />, and Frédéric Mouthereau<fnref n="2" />

*Correspondence to:* Margaux Mouchené (mmouchene@tulane.edu)

Received: 8 August 2016 – Discussion started: 17 August 2016

Revised: 9 January 2017 – Accepted: 30 January 2017 – Published: 21 February 2017

**Abstract.** Alluvial megafans are sensitive recorders of landscape evolution, controlled by both autogenic processes and allogenic forcing, and they are influenced by the coupled dynamics of the fan with its mountainous catchment. The Lannemezan megafan in the northern Pyrenean foreland was abandoned by its mountainous feeder stream during the Quaternary and subsequently incised, leaving a flight of alluvial terraces along the stream network. We use numerical models to explore the relative roles of autogenic processes and external forcing in the building, abandonment and incision of a foreland megafan, and we compare the results with the inferred evolution of the Lannemezan megafan. Autogenic processes are sufficient to explain the building of a megafan and the long-term entrenchment of its feeding river on time and space scales that match the Lannemezan setting. Climate, through temporal variations in precipitation rate, may have played a role in the episodic pattern of incision on a shorter timescale. In contrast, base-level changes, tectonic activity in the mountain range or tilting of the foreland through flexural isostatic rebound do not appear to have played a role in the abandonment of the megafan.

## 1 Introduction

Alluvial fans and megafans are prominent geomorphic objects of remarkably conical shape, constructed by the accumulation of sediments at the outlet of mountain valleys. They occupy a key position in the sediment routing system and, as such, have been widely used as recorders of external forcing on landscape evolution in a variety of settings. Controls on the building and incision of these deposits, through alternating phases of aggradation and erosion, have been shown

to be related to climatic changes (Barnard et al., 2006; Arboleya et al., 2008; Assine et al., 2014), tectonic activity (DeCelles and Cavazza, 1999), base-level oscillations (Havvey, 2002) or a combination of those factors (Abrams and Chadwick, 1994; Dade and Verdeyen, 2007; Schlunegger and Norton, 2014). Laboratory experiments reproducing alluvial fan dynamics have helped understanding of the respective roles of these controls on fan morphology, facies changes and cyclic erosion–deposition processes (Kim and Muto, 2007; Nicholas et al., 2009; Rohais et al., 2011; Guerit et al., 2014).

<footnote n="1">¹Université Grenoble Alpes, CNRS, ISTerre, CS-40700, 38058 Grenoble, France</footnote>

<footnote n="2">²GET, Observatoire Midi Pyrénées, Université de Toulouse, CNRS, IRD, 14 avenue E. Belin, 31400 Toulouse, France</footnote>

<footnote n="3">³Department of Geology, FCFM, Universidad de Chile, Santiago, Chile</footnote>

<footnote>ªnow at: Department of Earth and Environmental Sciences, Tulane University, New Orleans, LA 70118, USA</footnote>
Both analog and numerical modelling studies have shown evidence for autogenic processes that could be of critical importance in fan evolution (Humphrey and Heller, 1995; Coulthard et al., 2002; Nicholas and Quine, 2007). Temporary sediment storage on the fan results in cyclic behaviour, with alternating phases of deposition and incision in the absence of external forcing (e.g. Coulthard et al., 2002). This behaviour is expressed in the thresholds (in run-off, slope or shear stress) and defined and implemented in the models. A critical value must be reached and exceeded for transport to be effective; after some further time steps, this parameter value decreases below the threshold and deposition occurs again (Schumm, 1979; Roering et al., 1999; Whipple and Tucker, 1999; DiBiase and Whipple, 2011).

Another level of complexity, often overlooked in previous experiments with alluvial systems, comes from the strong coupling and feedbacks between the source catchment and the basin. The specific response time and amplitude of each part of the system to a given forcing may differ and this results in a complex, oscillating erosion signal (Densmore et al., 2007; Humphrey and Heller, 1995; Babault et al., 2005; Carretier and Lucazeau, 2005). Numerical modelling by Pepin et al. (2010) suggested that autogenic processes play a key role in the evolution of such a coupled system subjected to constant external forcing. For these authors, permanent autogenic entrenchment can occur in a coupled catchment–fan system without changes in boundary conditions and external forcing when (i) the transport threshold (critical shear stress) is significant and (ii) progradation is limited by an open boundary with fixed elevation (e.g. a large river system at the foot of the fan).

In the northern foreland of the Pyrenees (France), the Lannemezan megafan was built during the Miocene by the erosional products of the mountainous Neste river catchment, and was abandoned during the Quaternary (Mouchené et al., 2017). The respective roles of climate and tectonics in this evolution remain unresolved. In this study, we seek to test hypotheses on the mechanisms at play in the abandonment and incision of the Lannemezan megafan through numerical modelling of alluvial megafan construction and abandonment. Although the complexity of this natural case might not be fully reproduced by the numerical model, we run a series of model scenarios to explore the respective effects of potential forcing factors, including autogenic dynamics, climate change, tectonic tilting, and base level change, on trends and patterns of incision (time and space scales, amplitudes) of the megafan. Disentangling the respective signals of autogenic processes and allogenic forcing requires understanding of (i) the wavelength and amplitude of each signal, (ii) the possible buffering effects of the response times of the fan and of the mountainous catchment, and (iii) the amplification–reduction factors introduced by the coupling of the system.

## 2 The Lannemezan megafan

Whereas the drainage network in the Pyrenean range is regularly spaced and mostly transverse to the structural trend, rivers of the northwestern foreland spread in a radial pattern over the convex topography of large Miocene alluvial fans (Fig. 1). The Lannemezan megafan is the most prominent geomorphic feature of the northern Pyrenean foreland, with a surface of 13 000 km² and a mean slope of 0.3°. Its characteristic semi-conical shape is outlined by the Garonne River (to the south, east and north) and by the radial river network on its surface. It was built during the middle Miocene (Biteau et al., 2006), while orogenic activity within the Pyrenean range had already waned (Vergès et al., 2002; Sinclair et al., 2005). Molasse-type deposits of middle- to late-Miocene age, with rounded pebbles and boulders in an abundant clayey and sandy matrix, make up most of the megafan volume (Paris, 1975; Azambre et al., 1989). These are capped by the Lannemezan Formation, consisting of (i) an upward-fining, stratified clay and sand sequence that contains strongly weathered gravel and pebbles in a very fine matrix, dated by a hipparion-bearing fauna at its base as late Miocene to Pliocene ("Pontico-Pliocene"; Paris, 1975; Azambre et al., 1989) and (ii) a Quaternary sheet of very similar composition.

The Neste River exits the mountain range at the apex of the megafan and thus most probably provided the material for building the megafan. Comparing sediment volumes and thermochronology-derived exhumation rates, Mouchené (2016) showed that the relatively small Neste catchment could have exported a sufficiently large sediment flux to produce all of the megafan deposits. Contrary to the neighbouring Garonne valley, the Neste valley is believed to have mostly been exempt of glaciers during the Quaternary, and no glacial features are known on the megafan itself (e.g. Delmas, 2015).

The Lannemezan megafan is currently disconnected from its source catchment and is being incised. The Neste turns sharply to the east at the megafan apex and incises the fan head ~100 m vertically before merging with the larger Garonne River at its mountainous outlet (Fig. 1); this drainage pattern suggests the capture of the Neste by the Garonne. During incision, the rivers of the northern foreland (including the Neste, Garonne and fan rivers) left a series of alluvial terraces. Mouchené et al. (2017) dated the abandonment of the fan and the onset of incision at ≥ 300 ka from ¹⁰Be and ²⁶Al cosmogenic nuclide dating of the fan surface and terraces along the Neste valley. The episodic abandonment of these river terraces during incision may be related to changing fluvial dynamics during shifts between Quaternary cold and warm phases (Mouchené et al., 2017).
<figure><img src="image_2.png" /><figcaption>Figure 1. The Lannemezan megafan and Neste catchment in the central northern Pyrenees (inset map shows location in southern France).</figcaption></figure>

## 3 Model description

We use a recent version of the CIDRE code, which models landscape evolution in a continental setting (Carretier et al., 2015). We recall here the main characteristics of the code and refer the reader to Carretier et al. (2015) and references therein for further details.

At the beginning of each time step, a specified volume of water is distributed homogeneously over the cells making up the model surface. The propagation of water and sediment is performed in cascade, from the highest to the lowest cell and following decreasing elevation, to ensure mass conservation. A multiple-flow algorithm is used to propagate the

water flux to downstream cells proportionally to the slope in each direction (Murray and Paola, 1997; Coulthard et al., 2002; Carretier et al., 2009), allowing a distributary drainage pattern to develop.

### 3.1 Mass balance

During a time step $\partial t$, the elevation $z$ of a grid cell changes as follows:

$$ \frac{\partial z}{\partial t} = -\epsilon + D + U, \qquad (1) $$
where $\epsilon$ is a local erosion (detachment or entrainment) rate, $D$ is a local deposition rate and $U$ is an uplift or subsidence rate.

The local deposition rate $D$ is defined as

$$D = \frac{q_s}{L}, \qquad (2)$$

with $q_s$ the incoming sediment flux per unit width and $L$ the transport length. The transport length $L$ determines the proportion of incoming sediment flux that is deposited in the cell: a large $L$ results in little deposition, such as a steep slope or high water discharge would favour in natural settings. The cell outflux per unit width $q_s$ is the sum of the sediment detached from a given cell plus the sediment eroded upstream that crossed this cell without being deposited; it is thus non-local (e.g. Tucker and Bradley, 2010). This approach is generalized for both hillslope and fluvial processes by specifying $\epsilon$ and $L$ in both cases.

## 3.2 Hillslope processes

The approach used by Carretier et al. (2015) is different from the non-linear diffusion model proposed by previous authors (Roering et al., 1999; Carretier et al., 2009, 2014). Instead, in this model the elevation variation results from the difference between a local detachment rate and a deposition rate using Eqs. (1) and (2), where $\epsilon$ erosion rate and transport length $L$ are defined as

$$\epsilon = \kappa S \qquad (3)$$

$$L = \frac{dx}{1 - (S/S_c)^2}, \qquad (4)$$

where $\kappa$ is an erodibility coefficient, $S$ is the steepest slope and $S_c$ is a critical slope. If the slope is steeper than $S_c$, $\epsilon$ is set such that $S = S_c$. The detachment rate is proportional to the local gradient, but the deposition rate ($q_s/L$ in Eq. 2) depends on the slope and critical slope: when $S \ll S_c$, most of the sediment entering a cell is deposited there and when $S \sim S_c$, $L$ becomes infinity and there is no deposition on the cell.

## 3.3 Fluvial processes

For fluvial processes, a detachment algorithm including a threshold is used for sediment and bedrock:

$$\epsilon = K(k_t q^m S^n - \tau_c)^p \qquad (5)$$

$$L = \xi q, \qquad (6)$$

where $K$ is an erodibility coefficient, $q$ is water discharge per unit flow width on the cell, $S$ is slope, and the exponents $m$, $n$, and $p$ are positive. $k_t$ is the shear stress parameter so that $k_t q^m S^n = \tau$ (shear stress) and Eq. (5) takes the classic form of the excess shear-stress formula (Tucker, 2004). $\tau_c$ is the critical shear stress required for clast detachment. $p$ is set

to 1 in our experiments (following Lavé and Avouac, 2001). The transport length $L$ depends on particle size and density (included in the coefficient $\xi$). This law implies that the deposition rate decreases when the water discharge per unit width $q$ increases.

For fluvial processes, the flow width $w$ can be set to the cell width $dx$ or to a river width such as

$$w = k_w Q^{0.5}, \qquad (7)$$

where $k_w$ is a coefficient depending on the lithology and $Q$ is the total water discharge at a river section. Flow-width variation is critical in the modelling of alluvial-fan evolution because it plays a role in avulsion processes, in the changing flow dynamics (a change in flux geometry may lead to overflowing and a shift to distributive flow) and in incision patterns (leaving alluvial terraces in some cases).

## 3.4 Cover effect

Erosion of sediment is different from that of bedrock (Eqs. 3 and 5), and within the bedrock different layers can be defined by their respective erodibility and detachment or slope thresholds ($\kappa$ and $S_c$ for hillslope processes and $K$ for fluvial processes). During a time step $dt$, different layers can be eroded on a given cell: the erosion of each layer consumes part of $dt$ so that less time remains to erode the underlying layer. This time reduction is taken into account by multiplying $dt$ by $(1 - \frac{\text{volume layer}}{w} dx \epsilon dt)$ between layers. In this way, the “cover effect” of a sediment layer covering the bedrock (e.g. Whipple and Tucker, 2002; Lague, 2010) can be taken into account.

## 3.5 Lateral erosion

Flowing water can erode lateral cells, which are topographically above them and placed in a lateral direction perpendicular to each downstream direction. The lateral sediment flux $Q_{sl}$ is defined as a fraction of the flux in the considered direction (e.g. Murray and Paola, 1997; Nicholas and Quine, 2007):

$$Q_{sl} = \alpha Q_s, \qquad (8)$$

where $\alpha$ is a bank erodibility coefficient; it is specified for sediment and implicitly determined for bedrock layers, proportional to their erodibility (i.e. $\alpha_{sediment}/\alpha_{bedrock} = K_{sediment}/K_{bedrock}$ with $K$ from Eqs. 5 and 7).

## 4 Model set-up

The model simulates the evolution of a 100 × 150 km region split into a foreland zone (100 × 100 km) and an up-lifting mountain zone (100 × 50 km; Fig. 2). The grid cell size is 500 × 500 m. The dimensions are chosen to allow
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The model grid consists of a mountain section, subjected to constant uplift *U*, and of a flat foreland section; both sections are subjected to precipitation. B.C. is boundary conditions; neither water nor sediment can cross a closed B.C., while an open B.C. corresponds to transverse rivers of fixed elevation capable of transporting both sediment and water fluxes out of the grid.</figcaption></figure>

megafan building on an area matching that of the Lannemezan megafan and to permit competing catchments to develop during the drainage network growth phase; they correspond to a compromise between computing time and spatial resolution. Our model has much larger dimensions than previous experiments on coupled catchment–foreland systems (Tucker, 2004; Nicholas et al., 2009; Pepin et al., 2010; Langston et al., 2015) and the foreland/mountain width ratio is much higher than in previous work (Pepin et al., 2010). The initial surface is a horizontal grid with a Gaussian elevation noise ($\sigma = 0.5$ m) so we can study the system dynamics from the start of drainage network growth.

Although the convergence stopped in the early Miocene, the Pyrenees are still a high-relief mountain range, with long-term exhumation rates and present-day uplift rates both on the order of $0.1–0.3 \text{ mm yr}^{-1}$ (e.g. Jolivet et al., 2007; Nguyen et al., 2016), possibly controlled by the ongoing isostatic response to post-orogenic erosion. Because isostasy is currently not included in the model, we impose block uplift of the mountain part of the model at a constant rate of $0.3 \text{ mm yr}^{-1}$ (except in experiments where this is explicitly modified, see below). Homogeneous precipitation is applied at a constant rate over the entire model ($P = 1 \text{ m yr}^{-1}$); this parameter is modified in some experiments (Experiment 2a, b, c).

The sides of the mountain block (southern border and southern third of the eastern and western borders) are closed; neither water nor sediment can exit the grid through these. The other boundaries are open, corresponding to transverse rivers of fixed elevation (0 m), and are able to transport both sediment and water fluxes out of the grid. This is an approximation, as these transverse rivers will obviously have a downstream slope and therefore not be at a fixed constant

elevation. However, the major rivers of the Pyrenean foreland have very low gradient (<0.5%) and lie significantly lower than the streams incising the megafan.

We conducted a series of trial runs to adjust the relevant parameters in order to reproduce the first-order morphological traits of the northern Pyrenean foreland. In particular, the values for transport length (*L*) for the erodibility of bedrock and sediments (respectively *k*<sub>br</sub> and *k*<sub>all</sub>) and the critical shear stress (τ<sub>c</sub>) need to be established. These parameters are critical in the relief evolution but are generally poorly constrained. Giachetta et al. (2015) provided a compilation of values for erodibility of a set of lithologies in the context of the Iberian peninsula. However, these were used for models where the critical shear stress is zero and should be significantly different (approximately 1 to 2 orders of magnitude larger) when τ<sub>c</sub> > 0. The erodibility coefficient also depends on the value of *m* (e.g. Carretier et al., 2009); we therefore used a different value than that of Giachetta et al. (2015). In any case, the erodibility of sediment should be larger than that of the bedrock, the ratio between the two critically influencing the landscape morphology. We thus tested this ratio and the transport length *L* in order to reproduce the first-order characteristics of the northern Pyrenean landscape. Resemblance between the model and the ASTER DEM of the Lannemezan area was evaluated based on a number of morphological parameters (which need to agree within 30% to be accepted): in the range and foreland we assessed maximum, minimum, and mean elevations, river spacing, and relief (i.e. valley-to-ridge elevation difference); in the foreland we also assessed the length, width and northward slope of the megafan. The duration of the megafan-building and -incision phases is also compared to the evolution of the Pyrenean foreland as described by previous authors (Crouzel 1957; Azambre et al., 1989; Mouchené, et al., 2017). The best-fit parameters used for the model runs are presented in Table 1. Pepin et al. (2010) suggested that the critical shear stress should be significant for permanent incision to occur. We thus fixed a positive value for τ<sub>c</sub> (15 Pa) following Lavé and Avouac (2001) and Pepin et al. (2010).

# 5 Results

## 5.1 Megafan building

We successfully reproduced the first-order morphology of a fluvial megafan constructed on a low-elevation, stable foreland from the erosional products of a slowly uplifting mountain-range-like block (Fig. 3). The drainage network initiates from the area of transition between the mountain and foreland blocks (Fig. 3, A). In the foreland, it propagates outward and fans aggrade. Evenly spaced rivers (every ~10 km) build small fans and progressively lengthen their watershed towards the hinterland through headward incision. The fans quickly merge into a bajada, on top of which the flow is distributive (Fig. 3, B). At around 7.65 Myr, the mountain range
<figure><img src="image_4.png" /></figure>

**Figure 3.** Building the megafan. Temporal evolution of the mean elevation change in the foreland and mountain and map views of the mountain and foreland landscapes (black lines are 100 m contour lines, water flux in blue shades) at the time steps marked with arrows, through the megafan building phase. (A) Drainage network initiates and propagates in the mountain block through headward incision, while sediments are deposited along the front by regularly spaced streams. (B) Deposits merge in the foreland to form a bajada fed by a decreasing number of rivers as the mountain streams enlarge their basins. (C) Mean elevation in the range stabilizes and aggradation continues in the foreland, dominated by outflux of a central, main channel as the more lateral streams are drained directly toward the borders. (D) As aggradation continues, limited incision can occur along the borders of the fan (here on the western border) but (D') those streams are quickly refilled. (E) Similarly, temporary incision can happen near the apex. (F) After 15.3 Myr the megafan is built.

becomes fully connected (i.e. all cells of the mountain block are connected to the base level through the river network) and the mountain outflux is dominated by a few large rivers (~ five). In the meantime, aggradation continues in the foreland with a markedly conical pattern. The rivers situated at the easternmost and westernmost ends of the mountain range bend sharply to follow an along-strike course and quickly reach the open model boundaries, constrained by their short distance to a base-level outlet. In the following time steps, their watershed will increase in size by retreat of the drainage divide towards the middle of the range and the mountainous outlets of these streams will migrate towards the nearest border (Fig. 3, C). At this stage, the foreland deposits

are mainly provided by a single central channel, the flow of which distributes sediments largely over the whole foreland, now clearly defining a megafan (the flow spanning 180° over the foreland).

Several episodes of temporary entrenchment (< 50 m) occur during the building phase. They either concern the lower parts of the fan being incised by headward incision (Fig. 3, D) or the apex being incised by the main stream (Fig. 3, E). In both cases, within a few hundred thousand years the main stream has brought sufficient material to the entrenched zone to refill it and to overflow and become distributive again (Fig. 3, D and D'). This cyclic pattern is expected on
<table><caption>**Table 1.** Parameters used in the experiments. (a) shows the fixed parameters for all model runs. $\tau_c$ is the critical shear stress; $K_{br}$ and $K_{all}$ are the bedrock and sediment erodibility respectively; $m$, $n$, and $p$ are coefficients for the fluvial erosion law; $L$ is the transport length; and $\alpha$ is the lateral erosion coefficient. (b) shows the model settings for the experimental runs and results.</caption><thead><tr><th>(a)<br>τ<sub>c</sub></th><th>K<sub>br</sub></th><th>K<sub>all</sub></th><th>m</th><th>n</th><th>p</th><th>L</th><th>α</th></tr></thead><tbody><tr><td>15 Pa</td><td>0.5 10<sup>-3</sup></td><td>4 10<sup>-3</sup></td><td>0.6</td><td>0.7</td><td>1</td><td>0.3</td><td>0.01</td></tr></tbody></table><table><thead><tr><th>(b)<br>Settings</th><th colspan="7">Results</th></tr><tr><th></th><th>Experiment number</th><th>Precipitation rate (m yr<sup>-1</sup>)</th><th>Precipitation occurrence (fraction of time step)</th><th>Base level (m a.s.l.)</th><th>Uplift rate (mm yr<sup>-1</sup>)</th><th>Permanent entrenchment</th><th>Time of entrenchment after end of building at 15.3 Myr (kyr)</th></tr></thead><tbody><tr><td></td><td>1-default</td><td>1</td><td>1</td><td>0</td><td>0.3</td><td>Yes</td><td>270</td></tr><tr><td rowspan="4">Climate</td><td>2a1</td><td>0.5</td><td>1</td><td>0</td><td>0.3</td><td>No</td><td>–</td></tr><tr><td>2a2</td><td>2</td><td>1</td><td>0</td><td>0.3</td><td>Yes</td><td>180</td></tr><tr><td>2b</td><td>1</td><td>sinusoidal</td><td>0</td><td>0.3</td><td>No</td><td>–</td></tr><tr><td>2c</td><td>1</td><td>0.5</td><td>0</td><td>0.3</td><td>Yes</td><td>150</td></tr><tr><td>Base level</td><td>3a</td><td>1</td><td>1</td><td>–50</td><td>0.3</td><td>Yes</td><td>250</td></tr><tr><td rowspan="3">Uplift</td><td>4a</td><td>1</td><td>1</td><td>0</td><td>0</td><td>Yes</td><td>500</td></tr><tr><td>4b</td><td>1</td><td>1</td><td>0</td><td>0.1</td><td>Yes</td><td>310</td></tr><tr><td>4c</td><td>1</td><td>1</td><td>0</td><td>1</td><td>Yes</td><td>500</td></tr><tr><td rowspan="2">Tilting</td><td>5a</td><td>1</td><td>1</td><td>0</td><td>0 to 0.68</td><td>No</td><td>–</td></tr><tr><td>5b</td><td>1</td><td>1</td><td>0</td><td>0 to 2</td><td>No</td><td>–</td></tr></tbody></table>

megafans (Leier et al., 2005) and shows that the code mimics the natural fluvial dynamics of these settings.

In the long term, the mean elevation stabilizes in both the foreland and the mountain. From about 9 Myr, mean elevation stabilizes in the range (Fig. 3, B) but the elevation change remains slightly positive, which means that the relief is eroded at a slower rate than the applied uplift rate (i.e. true topographic steady state is not reached). Aggradation continues in the foreland, although at a slow pace (0.015 to 0.02 mm yr$^{-1}$ of mean elevation change); the timescale of aggradation is thus larger than that of the relief development. This is consistent with the observations of Babault et al. (2005) from an analog model. For them, aggradation in the foreland influences erosion of the mountain range by modifying the relative uplift rate (i.e. the difference between the uplift applied to the mountain block and the aggradation rate). Erosion of the range balances the continuously varying relative uplift rate, creating a “dynamic equilibrium” (Babault et al., 2005). A steady-state equilibrium (in which erosion rate equals uplift rate in the mountain) cannot be reached in such a landscape as long as aggradation occurs in the foreland.

## 5.2 Autogenic entrenchment

If the same conditions are maintained, natural entrenchment of the main stream occurs rapidly over a timescale that is 2 orders of magnitude smaller than the fan-building timescale (Fig. 4). Contrary to the building phase, during

which episodes of temporary entrenchment occurred but were followed by refilling and overflow, the incision starting at 15.3 Myr near the apex is sufficient to constrain the main stream avulsions to the eastern half of the fan for the subsequent time steps (Fig. 4, G, H, I). A small stream that developed on the eastern foot of the fan episodically captures the main flow and is thus progressively deepened and incised through headward incision.

Still, the main flow remains highly distributive and overflows this path several times before finally being permanently captured at around 15.57 Myr (Fig. 4, J). This event triggers a rapid incision phase, reaching nearly 150 m of incision close to the apex and larger amounts further downstream. The mean elevation change in the foreland drops dramatically upon entrenchment and the fan is subsequently eroded (mean elevation change remains negative for the rest of the experiment; Fig. 4). Erosion sharply increases in the mountain, especially in the watershed of the now connected main stream (Fig. 4). Figure 5 suggests that this incision leads to an increase in the relief due to rapid incision in the riverbed and little to no increased erosion on the hillslopes and ridges.

The main stream then laterally erodes its right bank in the foreland, tending towards an along-strike flow direction without further vertical incision (Fig. 4, K, L). Incision occurs in this bank and oblique to it, eventually capturing a secondary stream of the mountain range.
<figure><img src="image_5.png" /></figure>

**Figure 4.** Autogenic entrenchment. Temporal evolution of the mean elevation change in the foreland and mountain and map views of the mountain and foreland landscapes (black lines are 100 m contour lines, water flux in blue shades) at the time steps marked with arrows, through the autogenic incision phase (note change in timescale between this figure and Fig. 3). Starting from the landscape obtained at 15.3 Myr (see F in Fig. 3), the fan is incised when the main flow reaches the position of a small stream (e.g. G, I) but continues to grow (H) when the main flux overflows and migrates again on the fan. (J) After several of these cycles, the main flow is finally captured permanently in the stream. (K) As the main river now incises laterally towards the mountain front, a secondary stream is captured. (L) At the end of the experiment a large valley is incised along the front of the range. Topographic profiles across this valley (right panel) show that about 120 m of incision occurs in the foreland at the time of capture (black profiles, inset J); subsequently the valley is mostly enlarged by lateral erosion in the foreland, deepened and enlarged near the apex, and markedly deepened in the mountain (red profiles, inset L; note that horizontal scale is different for each profile).

## 5.3 External forcing

Subsequent experiments start from the topography obtained at the end of the “building phase” at 15.3 Myr and aim at evaluating the respective roles of different external factors in the incision pattern of the megafan. We consecutively explore the influence of changing parameters related to climate (precipitation rate and frequency of precipitation events), base-level change and tectonics (uplift rate and style). These models are run for 500 kyr to evidence the effects of external factors on this specific timescale, which corresponds to the abandonment and incision timescale of both the model and the Lannemezan megafan. Parameters used for these experiments are summarized in Table 1.

### 5.3.1 Precipitation rate and style

Decreasing (Experiment 2a1) or increasing (Experiment 2a2) the precipitation rate only results in decelerating or accelerating the processes observed in the original experiment. The same evolution is observed in experiment 2a1 as using the default settings, but the evolution is slower and the model does not reach the permanent entrenchment stage after the 500 kyr simulation. In the experiment with increased precipitation rate (2a2), erosion is enhanced and results in important widening of the valleys, but scattered deposition in the lower valleys create instabilities that perturb the model results.

In experiment 2b, we set the precipitation rate to follow a sinusoidal distribution with 100 kyr cycles to simulate the Quaternary climatic cycles. In this case, the trends of mean elevation change in the mountain and in the foreland are inversely correlated (Fig. 6). The mean elevation of the fore-
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Temporal elevation change of three locations in the mountain: bed of the main feeding river (red crosses), on nearby slope (black crosses) and ridge (blue crosses, see map view in right panel for locations). Following the entrenchment (marked with vertical grey line), the river rapidly incises, increasing (temporarily) the relief, as ridge elevations are not affected by the incision episode. The hillslope response is slow and lags behind that of the riverbed.</figcaption></figure>

land slightly increases through the experiment but remains stable in periods of low run-off, as the sediment supply from the mountain is halted. The mountain is eroded in periods of maximum run-off, whereas the elevation increases (at the uplift rate) in periods of minimum run-off. There is a slight delay in the mountain response to the variations in precipitation: as the run-off starts to increase, the elevation in the mountain continues to rise at the uplift rate for another time step (10 kyr) before it starts to decrease (Fig. 6). Similarly, there is a small lag between maximum run-off and minimum mean elevation change (Fig. 6). This delay corresponds to the response time of the mountain to cyclic precipitation rate changes and is consistent with works by Carretier and Lucezaue (2005) and Braun et al. (2015), who suggested a 1 to 30 kyr offset between forcing and response to rainfall variability at orbital (Milankovitch) timescales. In our experiment, the same delay is observed in the foreland, although the signal is less clear for periods of high run-off (Fig. 6).

At the end of the 500 kyr simulation, no permanent entrenchment is observed on the megafan. The small amount of incision that occurs on the fan when precipitation decreases is more than compensated for by renewed sediment influx from the mountain as precipitations start to increase again. The incision of the riverbed in the mountain in periods of high run-off is more than compensated for by the uplift in periods of low run-off (dominated by the applied uplift; Fig. 6).

### 5.3.2 Base-level change

A 50 m drop in base level is applied at the beginning of Experiment 3a. This leads to erosion in the foreland through headward incision of a number of streams, developing mostly on the western and eastern borders and persisting until the end of the experiment. Connection between the main feeder channel and the largest incising stream on the eastern border happens earlier than in the default model (at around 250 kyr) but the subsequent landscape evolution is very similar in both cases, although more incised streams remain at the end of this experiment (Fig. 7).

### 5.3.3 Uplift rate

Experiment 4a tests a scenario where uplift stops after the megafan building phase. In this experiment, the headward-incising stream connected to the eastern border captures a secondary river (at 300 kyr) before connecting to the main central channel at the end of the experiment (500 kyr, Fig. 8). The mean erosion rate in the range decreases steadily down to $0.19 \text{ mm yr}^{-1}$ (value for last 10 kyr of the experiment).

Increasing the uplift rate to $1 \text{ mm yr}^{-1}$ (Experiment 4c) quickly and permanently increases the elevation in the range without significantly increasing the aggradation in the foreland. This may be due to erosion (detachment and/or transport) not responding rapidly enough to catch up with this increase. Permanent entrenchment occurs at the end of the experiment (500 kyr) through the same process as in the default experiment.
<figure><img src="image_7.png" /></figure>

**Figure 6.** Experiment 2b (sinusoidal precipitation). Temporal evolution of (a) run-off and (b) mean elevation change in the foreland and (c) in the mountain. The trends of mean elevation change in the mountain and in the foreland are inversely correlated, and they show a slight delay relative to the change in run-off, corresponding to the response time. (d) Temporal elevation change of the central valley floor in the mountain (location marked by red cross on inset map) during Experiment 2b. The incision (blue arrows) related to humid periods (blue vertical lines are maximum run-off) is (over-) compensated for in drier periods (grey vertical lines correspond to no run-off) by the uplift, so that the elevation generally increases through the experiment. The maximum incision is delayed from the maximum run-off (response time). Temporary deposits at the valley outlet around 230 kyr dam the valley and trigger rapid backfilling of it, responsible for the high elevation between 230 and 310 kyr; the following incision episode removes this dam.

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

**Figure 7.** Model configuration at the end ($t = 15.80$ Myr) of Experiment 3a (50 m drop in base level at the onset of the incision phase). The megafan is incised by headward incision of a number of streams on its western and eastern borders (and marginally on the northern border).

### 5.3.4 Tilting experiment

In Experiment 5, we seek to reproduce the effect of isostatic rebound on the erosional pattern of the range and its foreland. At the moment, the CIDRE model does not include flexure. We thus chose to simulate the first-order effect of the flexural response to erosional unloading of the range through simple linear tilting of the model. This corresponds to an uplift pattern that increases linearly from zero at the northern boundary to a maximum fixed value at the southern boundary.

To scale the tilting to the observed geomorphic characteristics of the northern Pyrenean foreland, we estimate the potential tilting of the Lannemezan megafan since the onset of incision. We use a scaling law between fan area and fan slope to estimate the initial depositional slope of the Lannemezan Formation that caps the Miocene deposits. We use this formation because its base is the only mapped surface effectively preserved from erosion since deposition. We use a digitized geological map and an ASTER DEM (70 m resolution) to extrapolate the basal surface of the Lannemezan Formation using ArcGIS software and we estimate its current slope at 0.5°. Figure 9a shows area and slope data for the Lannemezan megafan compared to data compilations from active and inactive alluvial fans and megafans from the Alps, the Andes and the Himalayas (Horton and DeCelles, 2002; Guzzetti et al., 1997). The discrepancy of the Lannemezan
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Experiment 4a, the incised stream first connects through headward incision to the secondary river (left, *t* = 300 kyr) before being connected to the outlet of the central river at the end of the experiment (right, *t* = 500 kyr).</figcaption></figure>

data with the scaling law suggests an estimated ~0.4° tilt, which is simulated by uplift increasing linearly from 0 at the northern boundary to 2 mm yr⁻¹ at the southern boundary of the model.

It should be noted that this scaling relationship suggests a depositional angle of ~0.1°, which is within the range of values for megafans (e.g. DeCelles and Cavazza, 1999) but is not consistent with the slope observed in the default experiment (~0.4–0.5°). We compare this result with the tilt estimated using another often-used scaling relationship, between the catchment area and the fan slope (e.g. Champagnac et al., 2008). Figure 9b shows area and slope data for the Lannemezan megafan compared to data compilations from active and inactive alluvial fans from the Alps (Guzzetti et al., 1997; Crosta and Frattini, 2004; Champagnac et al., 2008). The discrepancy of the Lannemezan data with this scaling law only suggests an estimated 0.13° tilt, which will be simulated by uplift increasing linearly from zero at the northern boundary to 0.68 mm yr⁻¹ at the southern boundary. We test both these minimum (0.13°, Experiment 5a) and maximum (0.4°, Experiment 5b) tilt scenarios.

With the linearly increasing uplift, the megafan continues to grow; the mean elevation change in the foreland is steady, positive and higher than in the default experiment (> 0.2 mm yr⁻¹, including uplift). In both experiments, connection with the headward-incising stream and entrenchment occur (at ~280 kyr in the lower tilt experiment, and at ~240 kyr in the higher tilt experiment) but only temporarily affects this trend because as tilting continues the river outflows from this path (at ~320 and ~260 kyr respectively; Fig. 10). Deposition in the main path causes the stream to overflow from this channel and resume distributive flow over

the megafan, but instabilities in the models blur the results (Fig. 10). This suggests that overall tilting of the model prevents or limits permanent entrenchment. The mean elevation change in the mountain is steady and positive (~0.25 and 1.2–1.3 mm yr⁻¹ respectively), with peaks following the transient capture.

Figure 11 compares the evolution of a north–south topographic profile across Experiment 5a (0.13° tilt) and the default model. The megafan topography on this section is rather stable over the course of the default experiment. However, the megafan slope increases significantly in experiment 5a, showing that the tilt affects the megafan slope without being fully compensated for by erosion.

# 6 Discussion

## 6.1 Megafan building

In the model, the main steps of the megafan-building phases are (i) foreland deposition starts with small fans that quickly merge into a bajada; mountain watersheds merge so only a few streams are left; (ii) the rivers situated near the boundary change their direction to reach the shortest flow path to the border, and the central stream becomes the dominant provider for foreland sedimentation; and (iii) the megafan grows in response to cyclic flow dynamics (oscillating between channelized and distributive flow) and reaches dynamic equilibrium.

The timescale of the building phase of the megafan is long (> 10 Myr) when compared to active megafans of similar volumes deposited during the Quaternary (e.g. in the Alps, Andes and Himalayas; Assine et al., 2014; Fontana et
<figure><img src="image_10.png" /></figure>

**Figure 9.** (a) Scaling relationship between fan area and fan slope for alluvial systems of the Alps, the Andes and the Himalayas; data from Guzzetti et al. (1997; crosses) and Horton and DeCelles (2002; open circles). Thick line is the best power-law fit to the combined data: $S_f = 2.42 A_f^{(-0.30)}$. Tilt (vertical arrow) for the Lannemezan megafan (red circle) is estimated as the difference between present-day slope and predicted slope from the power-law fit. (B) Scaling relationship between drainage area and fan slope for alluvial systems of the Alps; data from Guzzetti et al. (1997; crosses), Crosta and Frattini (2004; triangles) and Champagnac et al. (2008; black circle = Valensole, grey square = Chambaran). The Lannemezan megafan–Neste system (red circle) lies slightly out of the relation (fit: $S_f = 10.4 A_b^{(-0.51 \pm 0.05)}$). Quaternary tilt (vertical arrow) of the Lannemezan megafan surface is estimated at $0.13^\circ$, the difference between present-day slope and slope predicted by the power-law fit. Modified after Champagnac et al. (2008).

al., 2014; Abrahami, 2015) but compares well to older systems (e.g. Campanian–Maastrichtian Hams Fork formation in Utah; DeCelles and Cavazza, 1999) and is consistent with the Lannemezan megafan building phase encompassing the early middle Miocene to Pliocene (i.e. ~ 15 Myr).

The long foreland (foreland length / mountain length = 2) allows for a large fan to develop but requires the model parameters to be set in a way that allows transportation over such great distance; in particular the parameter *L* must be large enough. This required longer *L*, which may be interpreted as a smaller settling rate (Davy and Lague, 2009), is consistent with the downstream fining of sediment in the Lannemezan megafan (Crouzel, 1957; Azambre et al., 1989), even though sediment fining is not accounted for in CIDRE.

The boundary conditions, open in the foreland and closed in the mountain, play a key role in the development and evolution of the drainage network. In particular, open boundary conditions on all three sides of the foreland allow for (i) the central river to become dominant in the sediment flux deposited in the foreland (thus creating a megafan) as the more lateral rivers rapidly adapt their course to the shortest path reaching the base level (along-strike to reach boundaries), and (ii) the conical shape to develop (contrary to Pepin et al., 2010, where cyclic boundary condition on lateral boundaries resulted in a more bajada-like landform). Open boundary conditions in the mountain would result in strike-parallel drainage, which shows that megafan building requires a relatively large range. Thus, in natural settings transverse rivers with efficient fluvial transport (to evacuate both water and sediments) appear necessary on all sides for a river–fan sys-

tem to be singled out and grow into a megafan deposit. In the northern Pyrenean foreland, the Garonne–Ariège and Adour rivers could have played this role, which suggests that they might have existed prior to the Miocene onset of the megafan building.

In our model, the absence of subsidence in the foreland may have encouraged the development of a fan covering a large area, which imposes overfilled conditions in the foreland basin. High subsidence rate would have allowed thick accumulation close to the range and thus limited its northward extension (e.g. Allen et al., 2013). This hypothesis could be tested with the addition of an algorithm for flexure (Simpson, 2006; Naylor and Sinclair, 2008). Nevertheless, the overfill hypothesis may be justified by the deceleration of subsidence rates in the Pyrenean retro-foreland since the Eocene (Desegaulx and Brunet, 1990; Desegaulx et al., 1990). In any case, the impact of varying subsidence rate on megafan growth and abandonment remains to be evaluated (e.g. Dingle et al., 2016).

## 6.2 Autogenic incision

### 6.2.1 Time and space scales

The autogenic entrenchment happens around 15.57 Myr, which, within the framework of the Lannemezan megafan evolution, is consistent with the fan-building phase encompassing the early middle Miocene to Pliocene and incision taking place during the Quaternary (since 300 ka; Mouchené et al., 2017).
Experiment 5a

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

Experiment 5b

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Tilting experiments 5a (0.13° tilt) and 5b (0.4° tilt). Connection with the headward incising stream occurs (at 280 and 240 kyr respectively) but model instabilities in the channel, interpreted as deposition, induce overflowing (at 320 and 260 kyr respectively). Distributive flow over the megafan resumes and lasts until the end of the experiment.</figcaption></figure>

At 15.3 Myr (end of the building phase), the mean elevation in the mountain range is about 1460 m, with a maximum elevation at 3160 m near the southern border, which is consistent with the northern flank of the Pyrenees. In the foreland, a maximum elevation of 950 m is reached at the fan apex in the model, which is higher than the current elevation of the Lannemezan megafan apex (~660 m), but the mean elevation of the foreland is around 270 m, comparable to that of the Lannemezan megafan. Megafan shape and dimensions (area) agree between our model and the Lannemezan megafan. At this point, the watershed of the main feeding river is about 1100 km², which is quite large when compared to the current Neste watershed (~750 km²). The scale of the vertical

entrenchment of the river is similar in the model and in the Lannemezan case (~100–150 m near the apex).

6.2.2 Mechanism and necessary conditions

For Pepin et al. (2010), autogenic entrenchment of an alluvial fan occurs only if (i) progradation is limited by the open boundary with fixed elevation and (ii) the transport threshold (critical shear stress) is significant. Nicholas et al. (2009) also suggested that declining aggradation in the fan results from increasing fan area during progradation (building phase), and incision is triggered by the lack of accommodation space when boundary conditions are reached. In nature, some in-
<figure><img src="image_13.png" /></figure>

**Figure 11.** Evolution of a north–south topographic profile near the middle of the model during Experiment 5 (region-wide tilting). In- set shows that with default settings (Experiment 1), the slope of the megafan does not increase; regional tilting (Experiment 5) is needed to create an increasing northward slope through time.

cised fans are linked to powerful transverse rivers (Milana and Ruzycki, 1999; Dühnforth et al., 2007, 2008) but the causality is not proven and external forcing is demonstrated in some cases (e.g. Dühnforth et al., 2008).

In our model, entrenchment naturally occurs, but it hap-
pens long after the moment when sediments reach the model
boundaries (Figs. 3 and 4). This is different from Pepin et
al. (2010) since autogenic entrenchment occurred precisely
when the sediment reached the free border in their experi-
ments. We suspect that this difference comes from the lateral
erosion included in our modelling and absent in the simu-
lations of Pepin et al. (2010). Lateral erosion limits the inci-
sion by fostering lateral migration and channel widening. Al-
though our modelling seems more realistic, the comparison
of this prediction with natural settings is not straightforward
because boundary conditions are likely to change over time
(e.g. Harvey, 2002). Consistent with the findings of Pepin et
al. (2010) and Nicholas and Quine (2007), all our models
that predict autogenic entrenchment use a significant critical
shear stress (entrainment threshold). We suspect this thresh-
old to control part of the incision magnitude and the delay
between the moment when sediment reaches the free border
and the moment when incision occurs. This aspect should be
further evaluated by varying the critical shear stress in other
experiments.

Van Dijk et al. (2009) proposed that aggradation on fans
allows a critical slope to be reached, triggering the incision.
However, in the tilting experiments (Experiment 5), the fan
slope reaches greater values than in the default experiment at
the time of entrenchment, but permanent entrenchment does
not occur. Therefore, attainment of a critical slope does not
appear to constitute a threshold for entrenchment. One ex-
planation may be that tilting fosters erosion in the mountain,
with a larger incoming sediment discharge entering the fore-
land, which prevents incision from growing.

6.2.3 Incision pattern

In the model, permanent entrenchment results from (lim-
ited) incision near the apex by the feeding river and major
headward incision of a stream from the foot of the fan until
both ends meet to define a continuously entrenched pathway
(Fig. 4). In the case of the Lannemezan megafan, we cannot
provide evidence to support or disprove this mechanism, but
the drainage pattern at the apex of the Lannemezan megafan
resembles the model (Fig. 12a and b). We could therefore
envisage the following scenario:

- A preexisting river Ariège-lower Garonne River flow-
  ing through the foreland to the east of the megafan,
  while the Neste river feeds the Lannemezan megafan
  through a distributive deposition pattern (Fig. 12c).

- A tributary of the ancestral Ariège-Garonne River retreats headward toward the apex of the megafan (Fig. 12d).

- Headward incision of the tributary toward the west leads to sequential capture of (a) the upper Garonne and (b) the upper Neste, abandonment of the fan, and rapid incision and terrace formation (preferably on the left bank of the now merged Neste-Garonne River due to southward river migration; Fig. 12e).

The amount of incision is already very important in the first
time step following the connection (~100 m near the apex)
and will only be further increased by another 80 m near the
apex. Downstream, the stream erodes its right bank towards a
strike-parallel pathway but does not incise vertically (Fig. 4).
This last characteristic resembles the lateral migration of the
Neste and Garonne rivers during their incision, evidenced by
the extensive alluvial terrace staircase left almost systemat-
ically on their left banks. In our model, the sediments of
the channel bed could enhance the effect of lateral incision
and inhibit further vertical incision through their cover effect
(see also Hancock and Anderson, 2002; Brocard and van der
Beek, 2006).

However, the terraces of the northern Pyrenean foreland
prove that the ~100 m incision of the Lannemezan megafan
at its apex was episodic, which contrasts with the pulse of in-
cision predicted by the model. Cosmogenic nuclide surface-
exposure dating suggests that these incision episodes in the
northern Pyrenean foreland are linked to cold-to-warm cli-
matic transitions (Mouchené et al., 2017).

6.3 Impact of climate change

In the model with sinusoidal precipitation rates (Experiment 2b), humid periods are characterized by erosion in the mountains and deposition in the foreland (with episodic incision); both decrease in drier periods because stream power decreases and less material is transported from the mountains. The wet-to-dry transition corresponds to a decrease in
<figure><img src="image_14.png" /></figure>

**Figure 12.** Comparison between (a) DEM of the northern Pyrenees and foreland at the longitude of the Lannemezan megafan and (b) the final model (contour lines every 100 m) shown at the same scale. The proposed evolution for the Lannemezan megafan is schematized below: (c) a preexisting river Ariège-lower Garonne River flows through the foreland while the Neste transports sediments deposited in the foreland through a distributive pattern to build a megafan; (d) a tributary of the Ariège-lower Garonne retreats headward and westward towards the apex of the Lannemezan megafan while the megafan keeps growing; (e) the migration of the tributary leads to sequential capture of (1) the upper Garonne and (2) the upper Neste and thus the abandonment of the megafan; rivers incise together in the foreland deposit leaving a series of alluvial terraces (preferably on the left bank due to the direction of river migration indicated by white arrows).

sediment input but also to a decrease in fluvial efficiency as the run-off nears zero, which prevents incision.

In the northern Pyrenean foreland, incision and abandonment of alluvial terraces has been linked to cold-to-warm climatic transitions (Mouchené et al., 2017) where the rapid decrease in sediment flux and gradual transitioning of the river to a single meandering thread, with a low width/depth ratio, would encourage vertical incision (e.g. Hancock and Anderson, 2002). Warm-to-cold transitions can also be associated with incision because of the increase in run-off vari-

ability and decline in vegetation that characterizes these periods; however, in nature, they are usually more gradual than cold-to-warm transitions. During glacial (dry, cold) periods, regolith is actively produced on hillslopes by efficient frost cracking but it is mobilized only at the onset of the following interglacial (wetter) period, when rainfall increases (e.g. Carretier et al., 1998). To reproduce and further explore this effect, we would need to include a climate (temperature)-dependent law for sediment production in the model.
In nature, incision is not always related to the return of wetter conditions; Meyer et al. (1995) suggest that incision of the terraces in their study site in northwestern Yellowstone National Park happens during warmer, more drought-prone periods because of the infrequent floods scouring the channel bed. Langston et al. (2015) recently modelled a similar pattern of incision by applying more intense, longer duration precipitation events during interglacial periods, but without changing the average precipitation rate. Periglacial processes have also been suggested to be a key controlling factor for erosion (e.g. Marshall et al., 2015; Dosseto and Schaller, 2016): erosion is enhanced during cold periods in regions where they occur, whereas it is enhanced during warmer periods in regions exempt of periglacial processes. Mass wasting processes could be the main driver for erosion increase during wet periods (e.g. Bookhagen et al., 2005), although their relationship to other environmental parameters, such as vegetation cover, remains disputed (e.g. Istanbulluoglu and Bras, 2005; Carretier et al., 2013; Dosseto and Schaller, 2016). Our current model does not to take such processes into account. Aggradation and incision thus seem to be controlled by the variability in rainfall intensity and event duration but also by temperature-dependent hillslope processes, rather than by mean precipitation rate alone.

A number of studies have related terrace incision with climate changes (e.g. Barnard et al., 2006; Bridgland and Westaway, 2008). This also seems to be the case in the northern Pyrenean foreland, where terrace abandonment was related to Quaternary climatic changes, although the model does not reproduce this pattern (it does not produce terraces at all). Several experiments suggest that the longer the foreland, the more it buffers the effects of short-period variations (Métivier and Gaudemer, 1999; Babault et al., 2005; Carretier and Lucazeau, 2005). Therefore, the effect of rapid climatic changes could be dampened by the large dimensions of the foreland in our model, preventing terrace formation. The lack of temperature-dependent processes in our experiments (glacial erosion, temperature-dependent regolith production) may also prevent terrace formation. Finally, the model resolution could be insufficient to resolve alluvial terraces.

## 6.4 Uplift rate

In the experiment where uplift stops after 15.3 Myr (Experiment 4a), the mountains erode at a rate of 0.19 mm yr⁻¹, comparable to the highest values obtained through estimation of basin-averaged erosion rates using cosmogenic nuclides in river sands (0.01 to 0.16 mm yr⁻¹; Mouchené, 2016). Uplift is thought to have significantly decreased in the Pyrenees since the Miocene, with modern GPS-derived uplift rates being small (0.1 ± 0.2 mm yr⁻¹ of differential uplift of the mountain belt with respect to a regional reference frame; Nguyen et al., 2016). Our results suggest that the Lannemezan megafan could have been built in a period of reduced tectonic uplift. The evolution of the piedmont is

very similar to that of the default experiment (where uplift is maintained at 0.3 mm yr⁻¹) except for the entrenchment that is refilled in experiment 4a. Thus, it appears that tectonic activity in the mountain belt does not strongly influence incision dynamics in the foreland.

## 6.5 Flexural isostatic rebound

We attempted to simulate the effect of flexural isostatic rebound on the incision pattern through tilting of the model. In the Alps, tilting of the foreland appears related to isostatic rebound in response to accelerated glacial erosion and possibly deep-seated geodynamic processes (Champagnac et al., 2008). This pattern has not been demonstrated for the Pyrenees. Although the simplistic approach we used does not reproduce the flexural response to erosional unloading of the range in detail, the slope of the fan topographic profile increases with time through this process, as suggested for alpine fans by Champagnac et al. (2008). Quantification of this increase in slope, although complicated by poor outcrop conditions, needs to be done in the northern Pyrenean piedmont to compare with the slope angles obtained in our model. In any case, tilting prevented permanent entrenchment in the experiment so this mechanism cannot explain the abandonment of a foreland megafan.

In the model, the topographic profiles merge downstream as a consequence of tilting. The alluvial terraces along the northern Pyrenean rivers also merge downstream and this pattern is also observed in the Alpine foreland. However, this pattern does not necessarily relate to tilting of the megafan: in other settings, this characteristic has been interpreted as a climatic imprint on incision (Poisson and Avouac, 2004; Wobus et al., 2010; Pepin et al., 2013). Thus, tilting does not appear to play a major role in the abandonment of the Lannemezan megafan.

## 7 Conclusions

Numerical modelling of the evolution of a catchment–foreland system has provided (i) new insight into the building and incision of a foreland megafan and (ii) key elements to infer the driving forces in the natural evolution of the remarkable Lannemezan megafan and its mountainous catchment in the northwestern Pyrenees.

For a megafan to develop, the foreland must be large enough to provide sufficient space for the fan to expand for a long period of time; a lack of subsidence may help this process. The role of preexisting transverse rivers flowing across the foreland seems to be critical in the building and incision of the megafan. They rapidly capture the closest streams exiting the range, which allows for a central mountainous stream to be singled out and to provide for most of the foreland deposits stacked in the megafan. In the northern Pyrenean foreland, the through-flowing Adour and Garonne–Ariège rivers may have helped shaping the Lannemezan megafan:
the spacing of these preexisting major drainage axes controls
the size of the fan, limit its extension and efficiently evacuate
water and sediments out of the megafan. The megafan grows
in response to the autogenic oscillations between sheet-flow
and channelized flow. These oscillations trigger small inci-
sions that are subsequently overfilled and rapid lateral move-
ment of the flow over the whole fan surface.

Permanent entrenchment of the Lannemezan megafan
could thus be the result of autogenic processes through
(i) progressive headward incision of a stream from the foot
of the fan (not too far from the apex) and (ii) final and rapid
incision of the apex once this stream has captured the feeding
river at its mountainous outlet. No external forcing is needed
to induce long-term entrenchment on the order of magnitude
observed in the field (100 m vertical incision near the apex)
but external factors cannot be ruled out. In particular, on a
shorter timescale, incision may have been influenced by Qua-
ternary climatic variations as suggested by the abandonment
of terrace staircases along the foreland rivers, incising the
Lannemezan megafan.

Variations in precipitation rate alone do not appear to be
sufficient to produce these episodic incision and alluviation
phases, and temperature-dependent hillslope processes may
also be involved. In contrast, base-level changes, tectonic ac-
tivity in the mountain range or tilting of the foreland through
flexural isostatic rebound appear to be unimportant factors in
the abandonment of the megafan.

8 Data availability

The source code is available upon request to Sébastien Car-
retier (sebastien.carretier@get.omp.eu). The input files for
the simulations described in this paper are available in
the GIT repository https://github.com/margauxmouchene/
CIDRE_Lannemezan_input.

**Competing interests.** The authors declare that they have no con-
flict of interest.

**Acknowledgements.** This study was supported by French National Research Agency ANR (Project PYRAMID, ANR-11-BS56-0031) and forms part of Margaux Mouchené's PhD thesis funded by the French Ministry of Higher Education (MESR). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). We thank Hugh Sinclair and an anonymous referee for constructive comments that helped clarify the paper.

Edited by: S. Castellort

Reviewed by: H. Sinclair and one anonymous referee

References

Abrahami, R.: Erosion et flux sédimentaires associés dans un bassin versant soumis à un régime de mousson: l’exemple de la Tista

(Inde), du Sikkim au Mégafan, PhD Thesis, Université de Greno-
ble, 223 pp. 2015.

Abrams, M. J. and Chadwick, O. H.: Tectonic and climatic implica-
tions of alluvial fan sequences along the Batinah coast, Oman, J.
Geol. Soc. Lond., 151, 51–58, 1994.

Allen, P. A., Armitage, J. J., Carter, A., Duller, R. A., Michael, N. A., Sinclair, H. D., Whitchurch, A. L., and Whittaker, A. C.: The Qs problem: Sediment volumetric balance of proximal foreland basin systems, Sedimentology, 60, 102–130, 2013.

Arboleya, M.-L., Babault, J., Owen, L. A., Teixell, A., and Finkel, R. C.: Timing and nature of Quaternary fluvial incision in the Ouarzazate foreland basin, Morocco, J. Geol. Soc. Lond., 165, 1059–1073, 2008.

Assine, M. L., Corradini, F. A., do Nascimento Pupim, F., and McGlue, M. M.: Channel arrangements and depositional styles in the São Lourenço fluvial megafan, Brazilian Pantanal wetland, Sediment. Geol., 301, 172–184, 2014.

Azambre, B., Crouzel, F., Debroas, E. J., Soulé, J. C., and Ternet, Y.: Carte géologique de la France à 1/50000, Feuille 1053 Bagnères-de-Bigorre, Bureau des Recherches Géologiques et Minières, Orléans, 1989.

Babault, J., Bonnet, S., Crave, A., and Van Den Driessche, J.: Influence of piedmont sedimentation on erosion dynamics of an uplifting landscape: An experimental approach, Geology, 33, 301–304, 2005.

Barnard, P. L., Owen, L. A., and Finkel, R. C.: Quaternary fans and terraces in the Khumbu Himal south of Mount Everest: their characteristics, age and formation, J. Geol. Soc., 163, 383–399, 2006.

Biteau, J. J., Le Marrec, A., Le Vot, M., and Masset, J. M.: The Aquitaine Basin, Petrol. Geosci., 12, 247–273, 2006.

Bookhagen, B., Thiede, R. C., and Strecker, M. R.: Late Quaternary intensified monsoon phases control landscape evolution in the northwest Himalaya, Geology, 33, 149–152, 2005.

Braun, J., Voisin, C., Gourlan, A. T., and Chauvel, C.: Erosional response of an actively uplifting mountain belt to cyclic rainfall variations, Earth Surface Dynamics, 3, 1–14, 2015.

Bridgland, D., and Westaway, R.: Climatically controlled river ter-
race staircases: A worldwide Quaternary phenomenon, Geomor-
phology, 98, 285–315, 2008.

Brocard, G. Y. and van der Beek, P. A.: Influence of incision rate, rock strength, and bedload supply on bedrock river gradients and valley-flat widths: Field-based evidence and calibrations from western Alpine rivers (southeast France), in: Tectonics, Climate, and Landscape Evolution, edited by: Willett, S. D., Hovius, N., Brandon, M. T., and Fisher, D., Geological Society of Amer- ica, Geological Society of America Special Paper 398, 101–126, doi:10.1130/2006.2398(07), 2006.

Carretier, S. and Lucazeau, F.: How does alluvial sedimentation at range fronts modify the erosional dynamics of mountain catchments?, Basin Res., 17, 361–381, 2005.

Carretier, S., Lucazeau, F., and Ritz, J.-F.: Approche numérique des interactions entre faille active, climat et érosion: exemple de la faille de Bogd, Gobi Altay, Mongolie, C. R. Acad. Sci. II A, 326, 391–397, 1998.

Carretier, S., Poisson, B., Vassallo, R., Pepin, E., and Farias, M.: Tectonic interpretation of transient stage erosion rates at different spatial scales in an uplifting block, J. Geophys. Res., 114, F02003, doi:10.1029/2008JF001080, 2009.
Carretier, S., Regard, V., Vassallo, R., Aguilar, G., Martinod, J., Riquelme, R., Pepin, E., Charrier, R., Hérail, G., Farías, M., Guyot, J.-L., Vargas, G., and Lagane, C.: Slope and climate variability control of erosion in the Andes of central Chile, Geology, 41, 195–198, 2013.

Carretier, S., Goddéris, Y., Delannoy, T., and Rouby, D.: Mean bedrock-to-saprolite conversion and erosion rates during mountain growth and decline, Geomorphology, 209, 39–52, 2014.

Carretier, S., Martinod, P., Reich, M., and Goddéris, Y.: Modelling sediment clasts transport during landscape evolution, Earth Surface Dynamics, 3, 1221–1254, 2015.

Champagnac, J.-D., van der Beek, P., Diraison, G., and Dauphin, S.: Flexural isostatic response of the Alps to increased Quaternary erosion recorded by foreland basin remnants, SE France, Terra Nova, 20, 213–220, 2008.

Coulthard, T. J., Macklin, M. G., and Kirkby, M. J.: A cellular model of Holocene upland river basin and alluvial fan evolution, Earth Surf. Proc. Land., 27, 269–288, 2002.

Crosta, G. B. and Frattini, P.: Controls on modern alluvial fan processes in the central Alps, northern Italy, Earth Surf. Proc. Land., 29, 267–293, 2004.

Crouzel, F.: Le Miocène continental du bassin d'Aquitaine, Bulletin du Service de la Carte Géologique de France 54, 264 pp., 1957.

Dade, W. B. and Verdeyen, M. E.: Tectonic and climatic controls of alluvial-fan size and source-catchment relief, J. Geol. Soc. Lond., 164, 353–358, 2007.

Davy, P. and Lague, D.: The erosion/transport equation of landscape evolution models revisited, J. Geophys. Res., 114, F03007, doi:10.1029/2008JF001146, 2009.

DeCelles, P. G. and Cavazza, W.: A comparison of fluvial megafans in the Cordilleran (Upper Cretaceous) and modern Himalayan foreland basin systems, Geol. Soc. Am. Bull., 111, 1915–1334, 1999.

Delmas, M.: The last maximum ice extent and subsequent deglaciation of the Pyrenees: an overview of recent research, Cuadernos de Investigación Geográfica, 41, 359–387, 2015.

Densmore, A. L., Allen, P. A., and Simpson, G.: Development and response of a coupled catchment fan system under changing tectonic and climatic forcing, J. Geophys. Res., 112, F01002, doi:10.1029/2006JF000474, 2007.

Desegaulx, P. and Brunet, M. F.: Tectonic subsidence of the Aquitaine Basin since Cretaceous times, B. Soc. Geol. Fr., 8, 295–306, 1990.

Desegaulx, P., Roure, F., and Villein, A.: Structural evolution of the Pyrenees: tectonic inheritance and flexural behaviour in the continental crust, Tectonophysics, 182, 211–225, 1990.

DiBiase, R. A. and Whipple, K. X.: The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate, J. Geophys. Res., 116, F04036, doi:10.1029/2011JF002095, 2011.

Dingle, E. H., Sinclair, H. D., Attal, M., Milodowski, D. T., and Singh, V.: Subsidence control on river morphology and grain size in the Ganga Plain, Am. J. Sci., 316, 778–812, 2016.

Dosseto, A. and Schaller, M.: The erosion response to Quaternary climate change quantified using uranium isotopes and in situ-produced cosmogenic nuclides, Earth-Sci. Rev., 155, 60–81, 2016.

Dühnforth, M., Densmore, A. L., Ivy-Ochs, S., Allen, P. A., and Kubik, P. W.: Timing and patterns of debris flow deposition on

Shepherd and Symmes creek fans, Owens Valley, California, deduced from cosmogenic $^{10}$Be, J. Geophys. Res., 112, F03S15, doi:10.1029/2006JF000562, 2007.

Dühnforth, M., Densmore, A. L., Ivy-Ochs, S., and Allen, P. A.: Controls on sediment evacuation from glacially modified and unmodified catchments in the eastern Sierra Nevada, California, Earth Surf. Proc. Land., 33, 1602–1613, 2008.

Fontana, A., Mozzi, P., and Marchetti, M.: Alluvial fans and megafans along the southern side of the Alps, Sediment. Geol., 301, 150–171, 2014.

Giachetta, E., Molin, P., Scotti, V. N., and Faccenna, C.: Plio-Quaternary uplift of the Iberian Chain (central-eastern Spain) from landscape evolution experiments and river profile modeling, Geomorphology, 246, 48–67, 2015.

Guerit, L., Métivier, F., Devauchelle, O., Lajeunesse, E., and Barrier, L.: Laboratory alluvial fans in one dimension, Phys. Rev. E, 90, 022203, doi:10.1103/PhysRevE.90.022203, 2014.

Guzzetti, F., Marchetti, M., and Reichenbach, P.: Large alluvial fans in the north-central Po Plain (Northern Italy), Geomorphology, 18, 119–136, 1997.

Hancock, G. S. and Anderson, R. S.: Numerical modeling of fluvial strath-terrace formation in response to oscillating climate, Geol. Soc. Am. Bull., 114, 1131–1142, 2002.

Harvey, A. M.: The role of base-level change in the dissection of alluvial fans: case studies from southeast Spain and Nevada, Geomorphology, 45, 67–87, 2002.

Horton, B. K. and DeCelles, P. G.: Modern and ancient fluvial megafans in the foreland basin system of the central Andes, southern Bolivia: implications for drainage network evolution in fold-thrust belts, Basin Res., 13, 43–63, 2002.

Humphrey, N. F. and Heller, P. L.: Natural oscillations in coupled geomorphic systems: an alternative origin for cyclic sedimentation, Geology, 23, 499–502, 1995.

Istanbulluoglu, E. and Bras, R. L.: Vegetation-modulated landscape evolution: effects of vegetation on landscape processes, drainage density, and topography, J. Geophys. Res.-Earth, 110, F02012, doi:10.1029/2004JF000249, 2005.

Jolivet, M., Labaume, P., Monié, P., Brunel, M., Arnaud, N., and Campani, M.: Thermochronology constraints for the propagation sequence of the south Pyrenean basement thrust system (France-Spain), Tectonics, 26, TC5007, doi:10.1029/2006TC002080, 2007.

Kim, W. and Muto, T.: Autogenic response of alluvial-bedrock transition to base-level variation: Experiment and theory, J. Geophys. Res., 112, F03S14, doi:10.1029/2006JF000561, 2007.

Lague, D.: Reduction of long-term bedrock incision efficiency by short-term alluvial cover intermittency, J. Geophys. Res., 115, F02011, doi:10.1029/2008JF001210, 2010.

Langston, A. L., Tucker, G. E., and Anderson, R. S.: Interpreting climate-modulated processes of terrace development along the Colorado Front Range using a landscape evolution model, J. Geophys. Res.-Earth, 120, 2121–2138, doi:10.1002/2014JF003403, 2015.

Lavé, J. and Avouac, J.-P.: Fluvial incision and tectonic uplift across the Himalayas of central Nepal, J. Geophys. Res., 106, 26561–26591, 2001.

Leier, A. L., DeCelles, P. G., and Pelletier, J. D.: Mountains, monsoons, and megafans, Geology, 33, 289–292, 2005.
Marshall, J. A., Roering, J. J., Bartlein, P. J., Gavin, D. G., Granger, D. E., Rempel, A. W., Praskievicz, S. J., and Hales, T. C.: Frost for the trees: Did climate increase erosion in unglaciated landscapes during the late Pleistocene?, Science Advances, 1, e1500715, doi:10.1126/sciadv.1500715, 2015.

Métivier, F. and Gaudemer, Y.: Stability of output fluxes of large rivers in South and East Asia during the last 2 million years: implications on floodplain processes, Basin Res., 11, 293–303, 1999.

Meyer, G. A., Wells, S. G., and Jull, A. J. T.: Fire and alluvial chronology in Yellowstone National Park: Climatic and intrinsic controls on Holocene geomorphic processes, Geol. Soc. Am. Bull., 107, 1211–1230, 1995.

Milana, J. P. and Ruzycki, L.: Alluvial-fan slope as a function of sediment transport efficiency, J. Sediment. Res., 69, 553–562, 1999.

Mouchené, M.: Évolution post-orogénique du système couplé piémont/bassin versant: le méga-cône alluvial de Lannemezan et son bassin versant au Nord des Pyrénées, PhD thesis, Université Grenoble Alpes, 346 pp., 2016.

Mouchené, M., van der Beek, P., Mouthereau, F., and Carcaillet, J.: Controls on Quaternary incision of the Northern Pyrenean foreland: chronological and geomorphological constraints from the Lannemezan megafan, SW France, Geomorphology, 281, 78–93, 2017.

Murray, A. B. and Paola, C.: Properties of a cellular braided-stream model, Earth Surf. Proc. Land., 22, 1001–1025, 1997.

Naylor, M. and Sinclair, H. D.: Pro-vs. retro-foreland basins, Basin Res., 20, 285–303, 2008.

Nguyen, H. N., Vernant, P., Mazzotti, S., Khazaradze, G., and Asensio, E.: 3-D GPS velocity field and its implications on the present-day post-orogenic deformation of the Western Alps and Pyrenees, Solid Earth, 7, 1349–1363, doi:10.5194/se-7-1349-2016, 2016.

Nicholas, A. P. and Quine, T. A.: Modeling alluvial landform change in the absence of external environmental forcing, Geology, 35, 527–530, 2007.

Nicholas, A. P., Clarke, L., and Quine, T. A.: A numerical modelling and experimental study of flow width dynamics on alluvial fans, Earth Surf. Proc. Land., 34, 1985–1993, 2009.

Paris, J.-P.: Carte géologique de la France à 1/50000, Feuille 1054 Montréjeau, Bureau des Recherches Géologiques et Minières, Orléans, 1975.

Pepin, E., Carretier, S., and Herail, G.: Erosion dynamics modelling in a coupled catchment–fan system with constant external forcing, Geomorphology, 122, 78–90, 2010.

Pepin, E., Carretier, S., Herail, G., Regard, V., Charrier, R., Farias, M., García, V., and Giambiagi, L.: Pleistocene landscape entrenchment: a geomorphological mountain to foreland field case, the Las Tunas system, Argentina, Basin Res., 25, 613–637, 2013.

Poisson, B. and Avouac, J. P.: Holocene hydrological changes inferred from alluvial stream entrenchment in north Tian Shan (northwestern China), J. Geol., 112, 231–249, 2004.

Roering, J. J., Kirchner, J. W., and Dietrich, W. E.: Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology, Water Resour. Res., 35, 853–870, 1999.

Rohais, S., Bonnet, S., and Eschard, R.: Sedimentary record of tectonic and climatic erosional perturbations in an experimental coupled catchment-fan system, Basin Res., 24, 198–212, 2011.

Schlunegger, F. and Norton, K. P.: Climate vs. tectonics: the competing roles of Late Oligocene warming and Alpine orogenesis in constructing alluvial megafan sequences in the North Alpine foreland basin, Basin Res., 27, 230–245, 2014.

Schumm, S. A.: Geomorphic thresholds: the concept and its applications, T. I. Brit. Geogr., 4, 485–515, 1979.

Simpson, G. D. H.: Modelling interactions between fold-thrust belt deformation, foreland flexure and surface mass transport, Basin Res., 18, 125–143, 2006.

Sinclair, H. D., Gibson, M., Naylor, M., and Morris, R. G.: Asymmetric growth of the Pyrenees revealed through measurement and modeling of orogenic fluxes, Am. J. Sci., 305, 369–406, 2005.

Tucker, G. E.: Drainage basin sensitivity to tectonic and climatic forcing: implications of a stochastic model for the role of entrapment and erosion thresholds, Earth Surf. Proc. Land., 29, 185–205, 2004.

Tucker, G. E. and Bradley, D. N.: Trouble with diffusion: Reassessing hillslope erosion laws with a particle-based model, J. Geophys. Res., 115, F00A10, doi:10.1029/2009JF001264, 2010.

van Dijk, M., Postma, G., and Kleinhans, M. G.: Autocyclic behaviour of fan deltas: an analogue experimental study, Sedimentology, 56, 1569–1589, 2009.

Vergés, J., Fernández, M., and Martinez, A.: The Pyrenean orogen: pre-, syn-, and post-collisional evolution, in: Reconструкtion of the Evolution of the Alpine-Himalayan Orogen, edited by: Rosenbaum, G. and Lister, G. S., Journal of the Virtual Explorer, 8, 57–74, 2002.

Whipple, K. X. and Tucker, G. E.: Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs, J. Geophys. Res., 104, 17661–17674, 1999.

Whipple, K. X. and Tucker, G. E.: Implications of sediment-flux-dependent river incision models for landscape evolution, J. Geophys. Res., 107, 2039, doi:10.1029/2000JB000044, 2002.

Wobus, C. W., Tucker, G. E., and Anderson, R. S.: Does climate change create distinctive patterns of landscape incision?, J. Geophys. Res., 115, F04008, doi:10.1029/2009JF001562, 2010.