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)

. Alluvial megafans are sensitive recorders of landscape evolution, controlled by both autogenic processes and allogenic forcing and 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 15 and 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 at 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 at a shorter time-scale. 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.


Introduction 20
Alluvial fans 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 25 oscillations (Harvey, 2002) or to 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 the respective roles of these controls on the 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). Both analog and numerical modelling studies Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 August 2016 c Author(s) 2016. CC-BY 3.0 License. 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 in the fan results in cyclic behaviour, with alternating phases of deposition and incision in the absence of external forcing (eg. Coulthard et al., 2002). This behaviour is expressed in the thresholds (in runoff, slope or shear stress) 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 5 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 of alluvial systems, comes from the strong coupling and feedbacks between the source catchment and the basin. Specific response amplitude and time of each part of the system to a given forcing may differ and this results in a complex, oscillating erosion signal Humphrey and 10 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 submitted 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) . 15 In the northern foreland of the Pyrenees (France), the Lannemezan megafan was built since the Miocene by the erosional products of a mountainous catchment, and was abandoned during the Quaternary (Mouchené et al., submitted). 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 20 natural case might not be fully reproduced by the numerical model, we will explore trends and patterns of incision (time and space scales, amplitudes) to understand the effect of different potential forcing factors (climate, tectonics, base level change, etc.) on landscape evolution. 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 25 system.

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 30 13,000 km 2 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. Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 August 2016 c Author(s) 2016. CC-BY 3.0 License.
Molasse-type deposits of middle-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) a well-sorted, stratified clay and sand sequence, containing strongly weathered gravel and pebbles, dated by a Hipparion-bearing fauna at its base as latest Miocene to Pliocene ("Pontico-Pliocene"; Paris, 1975;Azambre et al., 1989), and (ii) a Quaternary sheet of very similar composition. 5 The Neste River exits the mountain range at the apex of the megafan and thus most probably provided the material building the megafan. However, this stream now turns sharply to the east at this point and incises the fan head ~100m vertically, before merging with the larger Garonne River at its mountainous outlet (Fig. 1). The capture of the Neste by the Garonne and related abandonment of the fan was dated by Mouchené et al. (submitted) at ~300ka, from 10 Be and 26 Al cosmogenic nuclide dating of the fan surface. 10 During incision, the rivers of the northern foreland (including the Neste, Garonne and fan rivers) left a series of alluvial terraces, the episodic abandonment of which was dated by these authors and may be related to changing fluvial dynamics during shifts between Quaternary cold and warm phases (Mouchené et al., submitted).

Model description
We use a recent version of the CIDRE code, which models landscape evolution in a continental setting (Carretier et al., 15 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 20 to downstream cells proportionally to the slope in each direction (Murray and Paola, 1997;Coulthard et al., 2002;Carretier et al., 2009). This allows for a distributary drainage pattern to develop.

Mass balance
During a time step , the elevation z of a point of a cell changes as follows: where ϵ is a local erosion (detachment or entrainment) rate, D is a local deposition rate and U an uplift or subsidence rate.
The local deposition rate D is defined as: with q s the incoming sediment flux per unit width and L the transport length.

4
The transport length L determines the proportion of incoming sediment flux that is deposited in the cell: a large L results in a little deposition, 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 and that crossed this cell without being deposited; it is thus non-local (e.g. Tucker and Bradley, 2010). This approach is generalised for both hillslope and river processes by specifying ϵ and L in both cases. 5

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;. Instead, in this model, the elevation variation results from the difference between a local detachment rate and a deposition rate using Equations 1 and 2: Where is an erodability coefficient, S is the steepest slope and S c is a critical slope. If the slope is steeper than S c , is set such that S = S c . The detachment rate is proportional to the local gradient, but the deposition rate (q s /L in Equation 2) depends on the slope and critical slope: when S << S c , most of the sediment entering a cell is deposited there and when S ~S c , L becomes infinity and there is no deposition on the cell. 15

Fluvial Processes
For fluvial processes, a detachment algorithm including a threshold is used for sediment and bedrock: where K is an erodibility coefficient, q is water discharge per unit flow width on the cell, S is slope and the exponents are positive. k t is the shear stress parameter so that ! ! ! = (shear stress) and Equation 5 takes the classic form of the excess shear-stress formula (Tucker, 2004). ! is the critical shear stress to be reached for clast detachment. p is a positive coefficient, 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 ). This law implies that the deposition rate decreases when the water discharge per 25 unit width q increases.
For fluvial processes, the flow width can be set to the cell width dx or to a river width such as: Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 August 2016 c Author(s) 2016. CC-BY 3.0 License.
where ! 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 flowing), and in incision patterns (leaving alluvial terraces in some cases). 5

Cover effect
Erosion of sediment is different from that of bedrock (Equations 3 and 5), and, within the bedrock, different layers can be defined by their respective erodibility and detachment or slope thresholds ( 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 10 ) 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.

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. 15 Murray and Paola, 1997;Nicholas and Quine, 2007): where is a bank erodibility coefficient; it is specified for sediment and implicitly determined for bedrock layers, proportionally to their erodibility (i.e. α sediment /α bedrock = K sediment /K bedrock with K from Equations 5 and 7).

Model setup 20
The model simulates the evolution of a 100-by-150 km region split into a foreland zone (100x100 km) and an uplifting mountain zone (100x50km; Fig. 2). The grid cell size is 500x500 m. The dimensions are chosen to allow 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., 25 2009;Pepin et al., 2010;Langston et al., 2015) and the foreland/mountain width ratio is much higher than in previous work Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 August 2016 c Author(s) 2016. CC-BY 3.0 License. (Pepin et al., 2010). The initial surface is a horizontal grid with a Gaussian elevation noise (σ=0.5 m) so we can study the system dynamics from the start of drainage network growth.
The mountain part of the model is uplifted as a block at a constant rate of 0.3 mm y -1 (except in experiments where this is explicitly modified, see below), which seems reasonable for the late stages of evolution of the central Pyrenees (Jolivet et al., 5 2007;Nguyen et al., in review). Homogeneous precipitation is applied at a constant rate over the entire model (P = 1 m y -1 ); this parameter is modified in some experiments (Exp. 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 (0m) and able to transport both sediment and water fluxes out of the grid. 10 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 erodibilities of bedrock and sediments (respectively k br and k all ) and for the critical shear stress (τ c ) 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 erodibilities of a set of lithologies. However, these were used for models where the critical shear stress is zero and should be 15 significantly different (approximately one to two orders of magnitude larger) when τ c >0. Also, the erodibility coefficient depends on the value of m (e.g., Carretier et al., 2009) and we used a different value than that of Giachetta et al. (2015). In any case, the sediment erodibility 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: minimum, maximum and mean elevations of the range and foreland 20 megafan, relief and drainage network pattern. The 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 τ c (15 Pa) following Lavé and Avouac (2001) and Pepin et al. (2010).

Megafan building 25
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 towards the north and fans aggrade. Evenly spaced rivers (every ~10 km) build small fans and progressively lengthen their watershed towards the south through headward incision. The fans quickly merge into a bajada, 30 on top of which the flow is distributive (Fig. 3 B). At around 7.65 Ma, the mountain range 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 Earth Surf. Dynam. Discuss., doi: 10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 August 2016 c Author(s) 2016. CC-BY 3.0 License. by a few large rivers (~5). 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, probably 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). Meanwhile, the foreland deposits 5 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 (< 50m) 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 10 overflow and become distributive again (Fig. 3 D and D'). This cyclic pattern is expected on megafans (Leier et al., 2005) and shows that the code mimics the natural fluvial dynamics of these settings.
On the long term, the mean elevation stabilizes in both the foreland and the mountain. From about 9 Ma, mean elevation stabilizes in the range (Fig. 3 B) but 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 slow 15 pace (0.015 to 0.02 mm y -1 of mean elevation change), the time scale of aggradation is thus higher than that of the relief development. This is consistent with the observation of Babault et al. (2005) from an analog model. For them, aggradation in the foreland influences the erosion of the mountain 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 20 uplift rate in the mountain) cannot be reached in this landscape as long as aggradation occurs in the foreland.

Autogenic entrenchment
If the same conditions are maintained, natural entrenchment of the main stream occurs rapidly; over a timescale that is two 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 Ma near the 25 apex is sufficient to constrain the main stream avulsions to the eastern half of the fan for the subsequent time steps ( Figure 5 suggests that this incision leads to an increase of the relief due to rapid incision in the riverbed and little to no increased erosion on the hillslopes and ridges. The main stream then erodes laterally 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, which eventually captures a secondary stream of the mountain range. 5

External forcing
Subsequent experiments start from the topography obtained at the end of the "building phase" at 15.3 Ma and aim at evaluating the respective roles of different external factors on the incision pattern of the northern Pyrenean foreland. We subsequently explore the influence of changing the parameters related to climate (precipitation rate and frequency), base level change and tectonics (uplift rate and style). These models are run for 500 ky to evidence the effects of external factors 10 at this specific timescale, which corresponds to the abandonment and incision timescale of the Lannemezan megafan.
Parameters used for these experiments are summarized in Table 1.

Precipitation rate and style
Decreasing (Exp2a1) or increasing (Exp 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 15 evolution is slower and the model does not reach the permanent entrenchment stage after the 500 ky 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.
We set the precipitation rate of experiment 2b to follow a sinusoidal distribution with 100 ky cycles, to simulate the Quaternary climatic cycles. In this case, the trends of mean elevation change in the mountain and in the foreland are 20 inversely correlated (Fig. 6). The mean elevation of the foreland slightly increases through the experiment but remains stable in periods of low runoff, as the sediment supply from the mountain is halted. The mountain is eroded in periods of maximum runoff, whereas the elevation increases (at the uplift rate) in periods of minimum runoff. There is a slight delay in the mountain response to the variations in precipitation: as the runoff starts to increase, the elevation in the mountain continues to rise at the uplift rate for another time step (10 ky) before it starts to decrease (Fig. 6). Similarly there is a small lag 25 between maximum runoff 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 Lucazeau (2005) and Braun et al. (2015), who suggested 1 to 30 ky offset between forcing and response to rainfall variability at Milankovitch periods. In our experiment, the same delay is observed in the foreland, although the signal is less clear for periods of high runoff (Fig. 6).
Overall, during more humid periods, the mountain is eroded and the sediments are transported to the foreland, which 30 provides material for fan aggradation, whereas in periods of low precipitation, the material is not eroded and/or transported from the mountain to the piedmont (Fig. 6).
At the end of the 500-ky simulation, no permanent entrenchment is observed on the megafan. The small amounts of incision that occur on the fan when precipitation decreases are more than compensated 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 runoff is more than compensated for by the uplift in periods of low runoff (dominated by the applied uplift; Fig. 6)

Base-level change 5
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 thalwegs, developing mostly on the western and eastern borders and persisting until the end of the experiment. Connection between the main stream and the largest incising thalweg on the eastern border happens earlier than in the default model (at around 250 ky) but the subsequent landscape evolution is very similar in both cases, although more incised thalwegs remain at the end of this experiment (Fig. 7). 10

Uplift rate
Experiment 4a tests a scenario where uplift stops after the megafan building phase. In this experiment, the thalweg to the east border is continuously incised and connects with the outlet of a secondary river (at 300 ky) before connecting to the main central channel at the end of the experiment (500 ky, Fig. 8). The mean erosion rate in the range decreases steadily down to 0.19 mm y -1 (value for last 10 ky of the experiment). 15 Increasing the uplift rate to 1 mm y -1 (Experiment 4c) quickly and permanently increases the elevation in the range without increasing much 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 ky) through the same process as in the default experiment.

Tilting experiment 20
In this experiment (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 0 at the north boundary to a maximum fixed value at the south boundary.
To scale the tilting to the observed geomorphic characteristics of the northern Pyrenean foreland, we estimate the potential 25 tilting of the Lannemezan megafan. 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 a base surface using ArcGIS software and 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 30 megafans from the Alps, the Andes and the Himalaya (Horton and DeCelles, 2002;Guzzetti et al., 1997). The discrepancy of Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 August 2016 c Author(s) 2016. CC-BY 3.0 License.
the Lannemezan data with the scaling law suggests an estimated ~0.4° tilt, which will be simulated by uplift increasing linearly from 0 at the north boundary to 2 mm y -1 at the southern boundary.
It should be noted that this scaling relationship suggests a depositional angle of ~0.1°, which is not surprising for a megafan (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 oft-used scaling relationship, between the catchment area and 5 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 0 at the north boundary to 0.68 mm y -1 at the southern boundary. We test both these minimum (0.13°, Experiment 5a) and maximum (0.4°, Experiment 5b) tilt estimations. 10 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 y -1 , including uplift). In both experiments, connexion with the headward incising thalweg and entrenchment occur (at ~280ky in the lower tilt experiment, and at ~240ky in the higher tilt experiment) that only temporarily affects this trend because, as the tilt continues, the river outflows from this path (at ~320ky and ~260ky respectively; Fig. 10). The mean elevation change in the mountain is steady and positive (~0.25 mm y -1 and 1.2-15 1.3 mm y -1 respectively), with peaks following the capture.
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. Figure 11 compares the evolution of a north-south topographic profile across experiment 5a (0.13° tilt) and the default run. 20 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 by erosion. 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. 5 This required longer L, which may be interpreted as a smaller settling rate (Davy and Lague, 2009) . Also, 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 system to be singled out and 15 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 condition of 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 20 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 Eocene (Desegaulx and Brunet, 1990;. In any case, the impact of varying subsidence rate on megafan growth and abandonment remains to be evaluated.

Time and space scales
The autogenic entrenchment happens around 15.57 My, which, within the framework of the Lannemezan megafan evolution,  (Milana and Ruzycki, 1999;Dühnforth et al., 2007; but the causality is not proven (and 10 external forcing is demonstrated in some cases; e.g. Dühnforth et al., 2008).
In our model, entrenchment naturally occurs, but it happens long after the moment where sediments reach the model boundaries ( Fig. 3 and 4). This is different from Pepin et al (2010). In their experiments, autogenic entrenchment occurred precisely when the sediment reached the free border. We suspect that this difference comes from the lateral erosion included in our modelling and absent in simulations of Pepin et al. (2010). Lateral erosion limits the incision by fostering lateral 15 migration and channel widening. Although our modelling seems more realistic, it is difficult to compare this prediction with natural settings because boundary conditions are likely to change over time (e.g. elevation, length;Harvey, 2002)  One explanation may be that tilting fosters erosion in the mountain, with a larger incoming sediment discharge entering the foreland, which prevents incision from growing.

Incision pattern
In the model, permanent entrenchment results from (limited) incision near the apex by the feeding river and headward incision of a thalweg from the foot of the fan until both ends meet to define a continuously entrenched pathway (Fig. 4). In 30 the case of the Lannemezan megafan, we cannot provide evidence to support or disprove this mechanism but the drainage pattern resembles the model (Fig. 12 A and B). We could therefore envisage the following scenario: Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 August 2016 c Author(s) 2016. CC-BY 3.0 License. megafan through distributive pattern (Fig. 12 C); -tributary of Ariège/Garonne river retreat headwards toward the west (toward the apex of the megafan; Fig. 12 D) -migration of the tributary toward the west leading to sequential capture of (1) the upper Garonne and (2) the upper Neste, abandonment of the fan, rapid incision and terrace formation (preferably on left bank due to river migration; 5 Fig. 12 E).
The amount of incision is already very important in the time step following the connexion (~100m 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 strikeparallel pathway but does not incise vertically. This last characteristic resembles the right-lateral migration of the Neste and Garonne rivers during their incision, evidenced by the extensive alluvial terrace staircase left almost systematically on the 10 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 incision of the Lannemezan megafan (100-m at the apex) was episodic, which contrasts with the pulse of incision 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 climatic transitions 15 (Mouchené et al., submitted).

Impact of climatic change
In the model with sinusoidal precipitation rates (Experiment 2b), more 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 20 decreases and less material is being transported from the mountains. The wet-to-dry transition corresponds to a decrease in sediment input but also to a decrease in fluvial efficiency as the runoff nears 0, 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., submitted) 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, 25 2002).
Warm-to-cold transitions can also be associated with incision because of the increase in runoff variability and decline in vegetation that characterizes these periods, but 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 30 reproduce this, we would need to include a climate (temperature)-dependant law for production of sediment.
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 Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 August 2016 c Author(s) 2016. CC-BY 3.0 License. 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. Dühnforth et al. 2010;Dosseto and Schaller, 2016): erosion is enhanced during cold periods in regions where they occur; it is enhanced during warmer periods in regions exempt of periglacial processes. Mass wasting processes could be the main 5 driver for erosion increase during wet periods (e.g. Bookhagen et al. 2005b) 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 into account such processes. Aggradation and incision thus seem to be controlled by the variability in rainfall intensity and event duration but also by temperaturedependent hillslope processes, rather than by mean precipitation rate. 10 In nature, a number of studies relate terrace incision with climatic changes (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), so the effect of rapid climatic changes could be 15 dampened by the large dimensions of the foreland in our model, preventing terrace formation. The lack of temperaturedependent 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.

Uplift rate
In the experiment where uplift stops after 15.3 My (Experiment 4a), the mountains erode at a rate of 0.19 mm y -1 , 20 comparable to the highest values obtained through estimation of basin-averaged erosion rates using cosmogenic nuclides in the northern Pyrenees (0.01 to 0.16 mm y -1 ; Mouchené, 2016). Uplift is thought to have significantly decreased in the Pyrenees since the Miocene, with modern GPS-derived uplift rates being very small (0.1 ± 0.2 mm y -1 ; Nguyen et al., in review). 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 y -1 ) except 25 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.

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 is related to isostatic rebound in response to accelerated glacial erosion (Champagnac et al., 30 2008). This pattern was not 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 Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 August 2016 c Author(s) 2016. CC-BY 3.0 License. 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, in the experiment, tilting prevented permanent entrenchment (incision was only temporary) 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 5 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 evolution of the Lannemezan megafan.

Conclusions 10
Numerical modelling of the evolution of a catchment/foreland system has provided (i) new insight in 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 and a lack of subsidence may help this process. The role of pre-existing transverse rivers flowing across the 15 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 a megafan). In the northern Pyrenean foreland, the throughgoing Adour and Garonne/Ariège rivers may have helped shaping of the Lannemezan megafan by clearing the sediment and water fluxes out of the megafan. The megafan grows thanks to the autogenic oscillations between sheet-flow and channelized flow. These oscillations trigger small 20 incisions that are subsequently overfilled and rapid lateral movement 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 thalweg from the foot of the fan (not too far from the apex) and (ii) final and rapid incision of the apex once this thalweg has captured the feeding river at its mountainous outlet. No external forcing is needed to induce longterm entrenchment on the order of magnitude observed in the field (100-m vertical incision near the apex) but external 25 factors cannot be ruled out. In particular, on a shorter time-scale, incision may have been influenced by Quaternary 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 seem 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 activity in the mountain range or tilting of the foreland through flexural isostatic rebound appear unimportant. 30 Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam.  Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-44, 2016 Manuscript under review for journal Earth Surf. Dynam.             Inset shows that with default settings (experiment 1), the slope of the megafan does not increase, regional tilting (experiment 5c) is needed to create increasing northward slope through time.