Contrasting Geomorphic and Stratigraphic Responses to Normal Fault Development During Single and Multi-Phase Rifting

Understanding the impact of tectonics on surface processes and the resultant stratigraphic evolution in multi-phase rifts is challenging, as patterns of erosion and deposition related to older phases of extension are overprinted by the subsequent extensional phases. In this study, we use a one-way coupled numerical modelling approach between a tectonic and a surface processes model to investigate topographic evolution, erosion and basin stratigraphy during single and multi-phase rifting. We compare the results from the single and the multi-phase rift experiments for a 5 Myr period during which they experience equal amounts of extension, but with the multi-phase experiment experiencing fault topography inherited from a previous phase of extension. Our results demonstrate a very dynamic evolution of the drainage network that occurs in response to fault growth and linkage and to depocentre overfilling and overspilling. We observe profound differences between topographic and depocenter development during single and multi-phase rifting with implications for sedimentary facies architecture. Our quantitative approach, enables us to better understand the impact of changing extension direction on the distribution of sediment source areas and the syn-rift stratigraphic development through time and space.


INTRODUCTION
Unravelling the long-term interactions between surface processes and tectonics is key to understanding basin evolution; yet it remains challenging, especially for rift basins that are characterized by complex multiple phases of extensional histories. Insights into how normal fault arrays behave during multi-phase extension have been derived from numerous studies of seismic reflection, well and outcrop data (e.g., North Sea Rift: Whipp et al., 2014, Bell et al., 2014East African Rift System: Korme et al., 2004;Barmer Basin rift, India: Bladon et al., 2015; Gulf of Thailand: Morley, 2007;North West Shelf, Australia: Frankowicz and McClay, 2010). Also, scaled physical model experiments demonstrate the impact of pre-existing faults on the development of newly formed faults (e.g., Bellahsen and Daniel, 2005;Henza et al., 2010Henza et al., , 2011Chattopadhyay and Chakra, 2013;Henstra et al., 2015;Zwaan and Schreurs, 2017;Molnar et al., 2019;Maestrelli et al., 2020;Wang et al., 2021). Other papers however, indicate that fault reactivation may not depend only on multi-phase extension but on mantle and crustal weaknesses among other factors (e.g., Zwaan et al., 2021;Samsu et al., 2021). Nonetheless, resolving surface processes evolution in time and space and their response to fault network growth during multi-phase rifting remains unclear. Our lack of understanding of how erosionaldepositional patterns evolve under multi-phase extension is largely due to the fact that the topography and the stratigraphic patterns associated with the older phase of extension are overprinted by the subsequent extensional phase. For example, in the active Mygdonia Rift in northern Greece ( Figure 1B) understanding the interplay between tectonics and surface processes is ambiguous since patterns of erosiondeposition related to the older NE-SW extensional phase (Psilovikos, 1977;Dinter and Royden, 1993) are overprinted by the subsequent N-S extensional phase (Chatzipetros and Pavlides, 1998).
Conceptual models for the large-scale rift evolution are based on observations from simple rifts, i.e., rift basins that have evolved in response to a single-phase of extension. For marine/coastal environments, these evolutionary models suggest a long-term transition from overfilled basins during the early stages of rifting to underfilled basins during the later stages of rifting (Leeder and Gawthorpe, 1987;Prosser, 1993;Ravnas and Steel, 1998;Gawthorpe and Leeder, 2000). This pattern has been associated with an increase in fault displacement and localization of deformation from the "rift initiation stage" to the "rift climax" stage, as faults grow and gradually link (Cowie, 1998;Gupta et al., 1998). Nevertheless, the role of surface processes and their interaction with normal fault growth has not been considered in the above evolutionary pattern. Cowie et al. (2006) show how normal fault growth, interaction and linkage affects drainage network evolution and sediment dispersal in rift basins using a tectonic model that is coupled to a landscape evolution model. They describe the development of a major linked structure from a diffuse array of active fault segments that results in a simple rift-related topography. Here we elaborate on Cowie et al. (2006) by investigating the geomorphic and stratigraphic evolution of rift basins subjected to either simple and or multi-phase extensional histories, following a similar one-way coupled numerical modelling approach. We use rift-related topographies derived from the tectonic model of Finch and Gawthorpe (2017) that simulates the development of both single and multi-phase rift basins in three dimensions to drive the surface process model pyBadlands (Salles et al., 2018). We evaluate and compare landscape evolution and basin stratigraphy between single-phase rifts and multi-phase rifts that form in response to two phases of extension, with an angle between the two extension directions of 60°( Figure 1A), similar to what observed in several natural rifts such as in the Mygdonia Rift, Greece ( Figure 1B).
Our modelling study develops understanding of surface processes and tectonic interactions during multi-phase rifting by quantifying the magnitude and distribution of erosional and depositional processes, and by making a direct comparison to single-phase rifting. We address specific questions that concern: 1) drainage network evolution and the location of drainage divides across the developing rift, 2) spatial and temporal variations in sediment supply and the implications for syn-rift stratigraphy, and 3) shifts between erosional and depositional processes during fault array development in both single and multi-phase rifts. Our quantitative approach allows us to investigate the impact of normal fault growth on the distribution of sediment source areas and the stratigraphic development through time and space during single and multiphase rifting.

MATERIALS AND METHODS
We combine surface deformation maps produced in a selfconsistent tectonic code with a surface processes modelling code in order to explore the geomorphic and stratigraphic evolution of single and multi-phase rifts. Rift-related topographies resulting from sequential steps of vertical displacements are generated using a 3D discrete element model that is based on the modelling method of Finch and Gawthorpe (2017). In the tectonic model dynamic normal fault development occurs in response to either a single-phase of extension (single-phase experiment) or two non-colinear phases of extension (multi-phase experiment) ( Figure 2B, Supplementary Figures S1, S2). These vertical displacements, adjusted to maintain the tectonic models' surface at a constant mean elevation of 0 m, are read as input files into the surface process model pyBadlands (Salles et al., 2018).

Tectonic Model
The 3D discrete element model of Finch and Gawthorpe (2017) simulates the nucleation, propagation and linkage of normal faults in response to imposed far-field extension ( Figure 2A). Faults develop spontaneously in the crust and are defined as pairs of broken bonds between juxtaposed discrete elements (Finch et al., 2004). Once a bond breaks it is not allowed to heal. The continental crust consists of spherical elements and acts as an elastic-brittle-plastic plate floating hydrostatically on a dense and inviscid mantle held in equilibrium around a specified depth (cf. King et al., 1988;Finch and Gawthorpe, 2017). Elements within the upper crust interact through linear elastic repulsive-attractive forces whereas elements of the lower crust interact through linear viscous (Newtonian fluid) forces ( Figure 2A and Supplementary Figure S1). The model domain extends 60 km NS and 40 km EW and has a spatial resolution of about 50 m in both directions. The focus of this study is an evaluation of the geomorphic and stratigraphic response to topographic changes in an evolving rift, therefore an outline of the tectonic model methodology and the parameters used are presented in the Supplementary Materials (see "Tectonic model" and Supplementary Table S1), while a full description and the constitutive equations can be found in Finch and Gawthorpe (2017).
In the single-phase experiment, fault development occurs during a 5 Myr period in response to N-S extension ( Figure 2B). In the multi-phase experiment, there is first a 2 Myr period of N060 extension and then a second 5 Myr phase of N-S extension ( Figure 2B and Supplementary Figure S2). At the end of the first phase of extension there is 6.7% N060 extension that corresponds to an initial 3.35% N-S extension (see inset in Figure 2B). However, and in order to provide a framework for comparing between the models we have intentionally set the reference frame for extension parallel to the single rift case, i.e., N-S extension. This means we are comparing the geomorphic-stratigraphic evolution of two rift Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 748276 systems after equal amounts of N-S extension, however, in which the evolution of the multi-phase rift system is affected by tectonic structures inherited from an earlier phase of extension. The extension rate during both phases is 3.0 mm/ yr and the angle between the two extensional phases is 60°( Supplementary Figure S1). Rift topographies are produced by the tectonic model at 0.5 Myr intervals for both experiments (Supplementary Figure  S2). The tectonic model area is increasing as rifting progresses (Supplementary Figure S2), hence generating surfaces from increasing areas with time makes it difficult to understand the drainage network evolution and sedimentation. Using elements that exist within an identical area in the two models at the end of the experiment and mapping their elevations through time, allows to address surface processes response to fault growth (see also Supplementary Materials). This approach preserves changes in the throw on faults within the region and replicates interpretations that are based on using current topographies to infer past rifting, erosion and sedimentation. The generated topographies, adjusted to a mean elevation providing a sea level at 0 m at each time interval, are used to construct maps of incremental surface subsidence and uplift (following Cowie et al., 2006), which are then passed as input parameters to the surface processes model.

Surface Processes Model
We use the surface processes model pyBadlands (Salles and Hardiman, 2016;Salles et al., 2018) to calculate erosion, sediment dispersal and deposition during single and multiphase rift development. The amount of N-S extension is the same for both experiments ( Figure 2B), which allows us to compare the surface processes model results. pyBadlands uses a mass balance approach to allow sediment transport and deposition under varying tectonic and climatic forcing. It integrates hillslope diffusion and river incision by means of the stream power law given by: ė k d P l (PA) m S n where ė, is the erosion rate, A is the upstream drainage area, S, is the slope, k d is the erodibility coefficient, P is net precipitation, and m, n and l are positive empirical coefficients. Sediment deposition occurs offshore or in topographic depressions perched above a user defined base level following Planchon and Darboux (2002).
Climatic variability is not considered in our model and we assume a constant precipitation rate of 1 m/yr. We also assume a spatially uniform bedrock erodibility. Bedrock erodibility coefficient, k d , is set to 2 × 10e −6 which has been previously used to simulate landscape development in rift settings composed of a mixture of carbonate and clastic sedimentary rocks (e.g., Pechlivanidou et al., 2018;Pechlivanidou et al., 2019). Moreover, we have not allowed ultimate base level to fluctuate over time and kept it constant at 0 m elevation. We acknowledge the fact that temporal and spatial variations in precipitation and bedrock erodibility would impact the overall erosion rates and sediment volumes produced in the models. However, having uniform precipitation rate and bedrock erodibility and also keeping base level constant allows us to focus on the impact of tectonics on surface processes during single and multi-phase rifting. Also, by performing a series of sensitivity tests using reduced and increased bedrock erodibility (k d 1 × 10e −6 and k d 4 × 10e −6 , respectively; see Supplementary Materials), we show that such changes do not impact our overall results concerning surface processes behavior in single versus multi-phase rifts. The model set-up represents a one-way coupling such that sediment erosion, transport and deposition do not affect the tectonic evolution of the rift systems (Cowie et al., 2006). Although feedbacks between surface mass transfer and fault evolution have been demonstrated in extensional settings (e.g., Maniatis et al., 2009;Olive et al., 2014), the impact of surface processes on normal faulting in our models is considered to be minimal since there is no sufficient space created to accommodate large sediment loads (sensu Zwaan et al., 2018).
To simulate rift evolution during the single-phase rifting we use an initial surface topography with a random roughness (Supplementary Figure S3A) that is produced by the tectonic model. For the multi-phase rift experiment, we use the topography and the drainage network configuration that develops from this initial stage at the end of the first phase of extension as the initial surface topography (Supplementary Figure S3B). A summary of the input parameters that we used in pyBadlands are shown in Supplementary Table S2.

Structural and Topographic Evolution
At the start of the single-phase rift experiment, a large number of isolated E-W trending fault segments nucleate across the model domain ( Figure 2B). These initial segments grow by means of tip propagation allowing many of them to have link already by ca. 1 Myr ( Figure 2B). However, extension remains distributed among the faults until ca. 2 Myr, and as a result most of the model domain remains submerged below ultimate base level during this first phase of the experiment (first panel; Figure 3A). From ∼2 Myr onwards, fault interaction and linkage progressively localize deformation onto a smaller number of linked fault systems, the total lengths of which varying between 10 and 40 km, limiting small grabens or halfgrabens ( Figure 2B). These major fault systems rapidly mature and accumulate topographic relief as observed from 2 to 2.5 Myr where most of the hanging wall basin area rises above ultimate base level ( Figure 3A). After 5 Myr, maximum elevation in the single-phase rift experiment is approximately 800 m ( Figure 4A).
Overall, topographic relief development drives higher erosion, and sedimentation rates that increase by an order of magnitude from 2 to 5 Myr (averaged over 0.  Figure 5 and capital letter "L" marks the depocentre shown in inset in Figure 10B. Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 748276 5 these large depocentres continue to evolve above ultimate base level, two major depocentres remain submerged until the end of the experiment ( Figure 3A). We have calculated sediment supply, Q s , for one of the major depocentres that remains submerged (see depocentre marked with letter "L" in Figure 3A) and compare it against fault-driven accommodation space, Q acc , and upstream drainage area (see inset in Figure 10B). Our results show that Q s almost equals Q acc after 5 Myr, allowing this depocentre to become nearly filled with sediment (see also the Discussion).

Drainage Network Evolution
In the single-phase rift experiment, drainage network develops from a state in which the subsiding hanging wall areas were submerged below ultimate base level (see first panel in Figure 3A). The geometry of the river network that arises around 2-2.5 Myr is strongly controlled by the fault pattern as the main rivers preferentially flow parallel to fault strike and the divides of the numerous drainage basins largely follow the crests of the major footwall blocks ( Figure 3A).
Although, the drainage divides are mostly fault-controlled and therefore more or less fixed, we observe a dynamic evolution of the drainage network as the hydrological connectivity between adjacent drainage basins changes over time. Some initially isolated drainage basins become hydrologically connected with their neighbors (red stars in Figure 3A), which is explained by two different mechanisms: fault linkage and basin overspilling. Linkage of adjacent fault segments and their associated depocentres, results in drainage integration events as illustrated in the example shown in Figure 5. Due to subsidence of an initial segment boundary topographic high, two submerged and initially isolated depocentres are integrated around 2.5-3 Myr ( Figures 5A,B). In other cases, the infilling and subsequent overspilling of a basin drives drainage integration. The most upstream basins shown in  Drainage network isolation and the development of internally-drained (endorheic) basins also occurs during the single-phase rift evolution (yellow stars in Figure 3A). Fluvial connections become interrupted by fault interactions, as for example shown in Figure 5, where endorheic conditions are re-established in two depocentres around 4.5-5 Myr. Even though both drainage integration and isolation events occur during the evolution of the single-phase rift experiment, there is an overall long-term trend of increasing drainage network connectivity.

Sediment Dispersal and Basin Stratigraphy
In order to understand stratigraphic evolution during the singlephase rifting, we extracted three N-S orientated stratigraphic profiles that cross the full model domain and the major depocentres ( Figure 6). Whereas some depocentres are controlled by two fault systems with opposing dip, the majority of the depocentres are asymmetric and characterized by a halfgraben geometry. These cross-sections also show that some faults become dominant at the expense of others that results in significant spatial variability in fault-related relief and the depth and elevation of the individual depocenters ( Figure 6).
Due to submergence of most of the model domain below ultimate base level during the first ca. 2 Myr, mean erosion and sedimentation Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 748276 8 rates (averaged over 0.5 Myr time intervals) remain very low during this initial stage of the single-phase rift experiment (<0.005 mm/yr; Figure 4B). From 2 Myr onwards, however, emergence of most parts of the model domain produces a marked increase in erosion and sedimentation rates ( Figure 4B). In turn, this results in a transition from marine/lacustrine to fluvial deposition. Stratigraphy varies significantly between the different hangingwall depocentres during the single-phase rifting ( Figure 6). There is a large variability in total sediment thickness among the depocentres, fluctuating between ∼50 m and ∼500 m, and as the majority of deposition occurs during the last 3 Myr of the experiment this implies long-term average sedimentation rates in the order of ∼0.015-0.15 mm/yr. Time plots of cumulative sediment thickness show a pronounced temporal variability in sedimentation rates within the deepest parts of each of the developing depocentres (see panels on the right-hand side in Figure 6). Although most of the basins show a gradual increase in sedimentation rates as the accumulation of fault-related relief drives higher erosion rates along the basin margins, there are some depocentres where sediment accumulation rates become constant (e.g., Depo. 2 and 5; Figure 6A and Depo. 1; Figure 6B). The observed temporal variability within individual depocentres also relates to drainage reorganization events that cause marked changes in the upstream drainage area and therefore, have an impact on the amount of sediment supply. For example, the cumulative sediment thickness in depocentre 3 in cross-section CC' (Figures 3A, 6C) increases rapidly after 3 Myr due to increased drainage area driven by multiple drainage reorganization events.

Structural and Topographic Evolution
In the multi-phase rift experiment, the second phase of extension commences from a state with fault structures (left most panel in Figure 2B) and associated topographic relief (mean elevation ∼50 m, Supplementary Figure S3) inherited from the first 2 Myr phase of extension. From the beginning of the second phase of extension, strain becomes distributed on both the pre-existing faults oriented oblique to the second phase of extension direction as well as a growing number of newly developing faults striking perpendicular to the extension direction ( Figure 2B). These two groups of faults link and develop characteristic zig-zag planform pattern of faulting, and both groups of faults remain active during the full experiment (see blue arrows in Figure 2B). So, despite the change in extension direction, the inherited structures continue to accommodate part of the deformation during the second phase of rifting and strongly influence the geometry of resulting basins.
Even though there is inheritance of topography in the second rift phase, maximum faulting-induced elevations remain low, < 250 m, during the first ∼2.5 Myr of the second rift phase ( Figure 4D). From ∼2.5 Myr onwards, however, there is a rapid increase in maximum elevations up to ca. 1,600 m ( Figures 4D,  7A). Mean erosion and sedimentation rates follow similar trends, with low values, < 0.01 mm/yr, during the first ∼2.5 Myr, followed by a rapid increase up to 0.04 mm/yr and 0.02 mm/yr respectively, after 5 Myr (Figures 4E, 7B). A key observation during the second rift phase is the development of a small number of very large depocentres limited by the interacting pre-existing and neo-formed faults from ∼3.5 Myr onwards, which become submerged below ultimate base level ( Figure 7A). The total surface area of these what we call "megadepocentres" with depocenter volumes >40 km 3 , reaches up to ca. 70% of the total depocenter surface area after 5 Myr ( Figure 4F). These mega-depocentres remain underfilled at the end of simulation, despite the pronounced increase in the average sedimentation rates after 2.5 Myr ( Figure 4E). For example, mega-depocentre marked with "L" in Figure 7A remains underfilled after 5 Myr as sediment supply, Q s , cannot outpace fault-controlled accommodation space, Q acc (see inset in Figure 10D in the Discussion).

Drainage Network Evolution
In the multi-phase rift experiment, the drainage network configuration is inherited from the first phase of extension (Supplementary Figure S3B). During the second phase of extension, this inherited drainage network progressively adjusts to the new structures developed as a result of the evolving fault network. The resultant drainage network is relative complicated with streams in some parts oriented parallel to the first rift phase faults, whereas in other parts streams are parallel to the newly formed second phase faults, resulting in characteristic hook-shaped rivers (e.g., blue major rivers at 5 Myr; Figure 7A).
Dynamic drainage reorganization in the multi-phase rift experiment is expressed by both drainage integration (red stars) and isolation (yellow stars) events between adjacent basins. For example, the depocentres shown in Figure 8 become fluvially integrated with one another between 3 and 4 Myr due to a combination of basin overfilling and structural linkage of depocentres. Both the drainage integration and isolation events modify the hydrological connectivity of the drainage network and produce large shifts in the dimensions of the main drainage basins. However, over the long-term there is clearly a progressive increase in the hydrological connectivity of the drainage network, leading to the formation of a smaller number of large catchments ( Figure 7A).

Sediment Dispersal and Basin Stratigraphy
While there is a large number of relative small depocentres (<20 km 3 ) around 2 Myr ( Figures 4F, 7B), only some of them develop into very large depocentres that are located in the hangingwalls of the most active faults ( Figure 7B). The formation of the very large mega-depocentres (>40 km 3 ) results from some depocentres merging with their adjacent ones located across-strike. The stratigraphic profiles in Figure 9 for example show that after 5 Myr, some of the largest depocentres of the multi-phase rift experiment, e.g., depocentre 2 in cross-section AA′, depocentre 3 in crosssection BB′, and depocentre 1 in cross-section CC' (Figure 9) formed by gradual across-strike linkage of initially individual depocentres. We also observe that some depocentres are  Figure 7B). For example, the depocentres shown in the x-y profile in Figure 8 develop above ultimate base level until ca 3.5 Myr in the immediate hanging wall of first phase faults, however, they are transformed into areas of erosion after ca. 4 Myr in response to fault linkage and rapid surface uplift. These shifts from deposition to erosion are easily observed in the stratigraphic record, where sediment accumulation terminates and is replaced by sediment bypass and erosion (Figure 9). The inverse may also occur where an area changes from erosional to depositional processes. In this experiment, an early area of erosion is transformed into an area of deposition once it subsides onto the hanging wall of a reactivated first phase fault that is linked to a second phase fault (see area marked with green dashed line at 2-2.5 Myr in Figure 7B). Furthermore, we observe cases where the locus of sediment accumulation changes over time, for example, relict deposits within depocentre 4 ( Figure 9A) indicate a shift in the location of sediment accumulation after ca 2.5 Myr.
There is significant variability between depocentres in terms of total sediment thickness (varying between 40 and 550 m for the depocentres shown in Figure 9) and, therefore, in long-term average sedimentation rates. A key observation is the low sedimentation rates for most depocentres until approximately 2-2.5 Myr (see plots of cumulative sediment thickness with time in Figure 9). For the major depocentres, these relatively low sedimentation rates are followed by a pronounced increase associated with the growth of topographic relief and the increase in average erosion rates. The increase in sedimentation rates is also Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 748276 associated with long-term increase in fluvial connectivity of the drainage network. For instance, depocentre 1i in cross section CC' shows a marked increase in sedimentation rates from 3.5 Myr in response to an upstream drainage area increase (of the order of 200 km 2 ) as a consequence of drainage integration ( Figures 7A,  9C). For many of the smaller depocentres, on the other hand, sedimentation slows down over time (e.g., depocentres 1 and 3; Figure 9A and depocentre 2; Figure 9B).

DISCUSSION
Most studies on multi-phase rifts focus on structural inheritance and controls on the structural style and evolution and not on the effects of multiple phases of extension on erosion-deposition and the resultant syn-rift stratigraphic evolution. The motivation of this study is to investigate the differences in topographic evolution, erosion and basin stratigraphy between single and multi-phase rifting, following a numerical modelling approach. We compare the results from two experiments for a 5 Myr period during which they experience equal amounts of N-S extension, but with one experiment (i.e., the multi-phase experiment) characterized by fault structures inherited from a previous, 2 Myr phase of extension ( Figure 2B and Supplementary Figure S2). In the multi-phase rift experiment the overall fault growth mimics the multi-phase analogue experiments of Henza et al. (2010), Henza et al. (2011) and Wang et al. (2021) where two phases of noncolinear extension are involved. In their studies, the secondphase normal faults link with reactivated first-phase faults leading to the development of zig-zag fault patterns, similar to the results from the fault growth model that we use here. Moreover, in our fault growth model the second phase faults are significantly shorter than the first phase faults (see pink arrows in Figure 2B), which also matches multi-phase analogue model results (e.g., Henza et al., 2011). In the tectonic model that we use in this study, a readjustment period early in the second phase of extension is observed. First phase faults still have a component of growth oriented to the first extensional phase that is stored within them which does not instantaneously switch off when rift orientation changes. Hence, "hybrid" fault growth (cf. Rotevatn et al., 2018Rotevatn et al., , 2019 exists early in the second extensional phase with first phase faults propagating perpendicular to the first extensional phase and at the same time are growing through tip propagation in the second extensional phase. This readjustment phase can last ∼1-2 Myr in our tectonic model and it may not be observed in analogue models due to the fact that occurs relatively quickly by comparison. A schematic summary of the topographic, drainage network and stratigraphy evolution at early and later stages of rifting for the single and the multi-phase rift experiments are shown in Figure 10.

Topographic Evolution
During the early stages of rift evolution (i.e., 2 Myr), topographic relief remains very low in both experiments ( Figures 4A,D), however, maximum elevations are slightly higher during multi-phase rifting ( Figure 4D). This is due to inherited structures from the first phase of extension contributing to the second phase of extension (i.e., from ∼0.5 Myr, see Figure 2B) forming longer faults segments with larger throw compared to the single-phase rift. Even though this difference is subtle, it results in more than three times higher average erosion rates and, therefore, higher sediment flux in the multi-phase rift compared to the singlephase rift (i.e., 0.002 mm/yr and 0.007 mm/yr after 2 Myr, respectively; see dashed lines in Figures 4B,E). Overall, average sedimentation rates are three times higher during the early stages of multi-phase rifting compared to the single-phase rifting (i.e., 0.001 mm/yr and 0.004 mm/yr, respectively; see solid lines in Figures 4B,E). Furthermore, inheritance of structures from the first extensional phase allows the accumulation of fault offsets and associated relief to accelerate faster during the later stages of multi-phase rifting (i.e., ∼1,250 m and ∼750 m after 5 Myr, respectively; Figures 4A,D), and consequently, results in higher mean erosion and sedimentation rates compared to the single-phase rifting (i.e., mean erosion rates 0.025 and 0.04 mm/ yr, mean sedimentation rates 0.013 and 0.02 mm/yr at 5 Myr, respectively; Figures 4B,E).

Drainage Network Reorganization
A key characteristic of the single and multi-phase rifts is the strong structural control on drainage network. During the early stages of the single-phase extension (i.e., ∼2 Myr), drainage network development is limited to small, low relief catchments draining transversely fault controlled footwalls ( Figure 3A). In contrast, the multi-phase drainage network over the same time interval is characterized by larger catchments that drain relatively high pre-existing relief ( Figure 7A). The pattern of this well-developed drainage network shows a general direction that is parallel to the first phase faults. This is partially due to the small number of newly second phase faults that have formed by this time (see Figure 2B). However, given that the duration of the first extensional phase in our experiment is 2 Myr, this also suggests that more time is needed for the drainage network to readjust to the new state once the sediment routes have been established.
After 5 Myr, large axial river systems dominate the singlephase rift topography with relatively small transverse catchments draining the uplifting footwalls ( Figures 3A, 10B). Gawthorpe and Leeder (2000) show that during the linkage and through-going fault stage the drainage network is characterized by major axial drainages as well as by transverse catchments that mark breached fault segments, an observation that agrees with our model results for the single-phase rift. In the multi-phase rift experiment, on the other hand, fault interactions between first phase and second phase faults lead to the development of more complex drainage network patterns, such as characteristic hook-shaped patterns ( Figures 7A,  10D). Excellent examples of such drainage network patterns are found in the multi-phase Mygdonia Rift, where rivers are in some parts oriented parallel to the first phase, NW-SE faults, whereas in other parts are parallel to the second phase, E-W, faults ( Figure 1B). Although such patterns are a common feature of active, single-phase rifts that result from across-strike fault interactions (e.g., Eliet and Gawthorpe, 1995;Cowie et al., 2006), our modelling results show that the development of these drainage network patterns also emerges from fault interactions between first and second phase faults during multiphase rifting.
An important outcome of this study is the very dynamic evolution of the drainage network ( Figures 3A, 7A). In both rift experiments, fluvial connections between adjacent basins can, first of all, develop as a consequence of faults interaction that causes the basins to integrate. Faults propagation and linkage lead topographic ridges that initially separate adjacent basins to subside, and, therefore, fluvial connections are being established ( Figures 5B, 8B). A second mechanism leading to drainage integration between adjacent basins is the overfilling of the upstream basin with sediment, allowing it to overspill and establish a fluvial connection with its downstream neighbor ( Figures 5B, 8B). The importance of basin overfill for drainage integration has been inferred for natural extensional systems, for example the central Italian Apennines (Geurts et al., 2018(Geurts et al., , 2020 and various valley systems in the Basin and Range (e.g., Hilgendorf et al., 2020). The opposite trend towards isolation of basins also occurs in both rift topographies. In all the cases observed, drainage isolation is caused by fault interaction and linkage and uplift of new footwall topography across an initially ongoing river system. Interestingly, both drainage integration and isolation can occur at the same locality during the ongoing development of the fault network (e.g., Figure 5B). Overall, there is a gradual increase in the hydrologic connectivity of the drainage network in both experiments as faults become progressively linked. From 4 Myr onwards, however, the multi-phase drainage network remains fairly stable. Cowie et al. (2006) show that drainage stabilization marks the phase when fault arrays become fully linked and slip rates are higher and more uniform. Our results suggest that this phase is reached ∼1 Myr earlier during the multi-phase rift evolution, as pre-existing structures facilitate fault linkage.
Drainage reorganization has a direct impact on sediment flux and accumulation rates. Temporal variability in sediment accumulation rates reflects changes in the upstream area that occur in response to drainage reorganization events (Figures 6,  9). However, a profound difference between the single and multiphase rift experiments resulting from drainage reorganization is the formation of larger drainage catchments, and thus higher sediment supply into the basins formed during multi-phase rifting ( Figures 3A, 7A, 10B,D). show time evolution plots from 2.5 to 5 Myr of sediment supply, Q s (in km 3 ), fault-controlled accommodation creation, Q acc (in km 3 ), and upstream drainage area (in km 2 ) for two basins formed during the single and multi-phase rifting, respectively (for location see Figures 3A, 7A). Note that the upstream drainage area remains fairly constant and variations in sediment production are due to topographic relief development that leads to increase in catchment average erosion rates over time.

Depocentre Evolution and Implications for Sedimentary Facies Development
An important difference between single and multi-phase rift experiments is the development of significantly larger syn-rift depocentres during multi-phase rifting ( Figures 4C,F, 10B,D). The surface area of depocentres larger than >40 km 3 increases significantly in the multi-phase rift experiment from ca. 3 Myr onwards, reaching up to 70% of the total surface area ( Figure 4F). These "mega-depocentres" are not developing during the singlephase rift experiment, where there is a relatively equal contribution of depocentres with volumes <40 km 3 after 5 Myr ( Figure 4C).
Syn-rift stratigraphy also shows significant differences between single and multi-phase rift experiments. Our results show that sediment accumulates within syn-rift depocentres from the early stages of the second phase of multi-phase rifting, whereas sediment accumulation is limited during the single-phase rifting (Figures 6,9,10). Sediment that is deposited in the multi-phase rift during this stage is likely to be mature, as it is being transported over relative long distances from the well-developed drainage system inherited from the first stage of rifting (see 2 Myr in Figure 7A). During the later stages of rifting, sedimentary successions with variable thicknesses accumulate in both rift settings (i.e., 2-5 Myr, Figures 6, 9, 10B,D). However, in the single-phase rift experiment, the majority of depocentres that develop below ultimate base level become filled to spill-point or overfilled by 5 Myr (Figure 3A, 10B). In contrast, despite higher average sedimentation rates and significant increase in sediment accumulation over time in the multi-phase rift experiment, depocentres that develop below ultimate base level remain underfilled (Figures 7, 10D). This difference between the two rift settings highlights competing roles of sediment supply and fault-controlled subsidence in controlling syn-rift stratigraphic evolution. Our analysis shows that rift basins do not necessarily grow into underfilled basins as generic models suggest (c.f., Gawthorpe and Leeder, 2000), if sediment supply keeps pace with the formation of faultcontrolled accommodation during a simple, single-phase rift evolution. Inset in Figure 10B shows an example of temporal evolution of sediment supply, Q s (in km 3 ), and fault-driven accommodation space, Q acc (in km 3 ), for one representative basin from the single-phase rift experiment (for location see Figure 3A). Q s increases as higher topographic relief drives higher erosion rates from ca. 3 Myr onwards. Fault-controlled accommodation space Q acc also increases over time, however, it remains approximately constant from 4.5 Myr onwards. This allows Q s to almost equal Q acc after 5 Myr, in turn allowing the basin to become nearly filled with sediment. In more complex rift settings, such as the multi-phase rift experiment in this study, basins likely develop under sediment-starved conditions. For example, Q acc for the basin shown in inset in Figure 10D (for location see Figure 7A) increases significantly from 3 Myr as first phase faults link to second phase faults. The result is the development of a large "megadepocentre" (Q acc > 40 km 3 ), where sediment supply cannot outpace accommodation space and thus, the basin remains underfilled.
The profound differences in depocentre development between the single and multi-phase rift experiments have implications for syn-rift stratigraphic evolution. In the single-phase rift experiment, stratigraphic patterns are characterized by aggradational stacking and most depocentres are filled with sediment ( Figure 6). In depocentres which initially develop below ultimate base level and emerge later, sedimentary facies shift from marine/lacustrine to fluvial. These stratigraphic patterns have also been observed in natural extensional settings, such as the central Apennines, where the transition from lacustrine to fluvial sedimentation is commonly observed in basin stratigraphy (e.g., Geurts et al., 2020). In contrast, the development of mainly large depocentres (>40 km 3 ) during multi-phase rifting that gradually subside below ultimate base level, imply a shift in sedimentary facies from alluvial to lacustrine/marine. Sedimentary infill within these basins likely consists of reworked material as pre-existing topography is eroded during the subsequent extensional phase and areas of deposition change into areas of erosion ( Figure 10D). For example, sedimentary fill in the multi-phase Mygdonia Rift shows a transition from fluvial sediments (e.g., conglomerates, sandstones and red-beds) that were deposited at the hanging wall of first phase faults to deltaic/lacustrine sediments mainly deposited at the hanging wall of second phase faults (Psilovikos, 1977). At the southern margins of this rift, deposits of the first extensional phase are incising during the second phase of extension leading to the formation of large fan deltas in Lake Volvi ( Figure 1B).

CONCLUSION
This numerical modelling study investigates the geomorphic and stratigraphic evolution of rift basins that develop in response to single (single-phase rift experiment) and two phases of extension (multi-phase rift experiment). We compare the results from the single and multi-phase rift experiments which experience similar amounts of extension during a 5 Myr period, with the multiphase rift experiment characterized by structures inherited from a previous, 2 Myr phase of extension. We conclude that: 1) Dynamic drainage network characterizes single and multiphase rift evolution. Drainage integration events occur when adjacent depocentres combine in response to fault growth and linkage, or when depocentres become overfilled with sediment and overspill. Fluvial isolation and the formation of endorheic basins also occurs, however, in the long-term there is a progressive increase in hydrologic connectivity as faults link. 2) Inherited structures from the preceding extensional phase accelerate the accumulation of topographic relief and promote the development of large depocentres that become underfilled with sediment as fault-controlled accommodation outpaces sediment supply during multi-phase rifting. Conversely, lower relief and small to medium-sized depocentres that grow into overfilled basins dominate the single-rift topography.
3) Temporal variability in sedimentation rates in single and multi-phase rifts reflects changes in upstream drainage area that occur in response to drainage reorganization events. However, during multi-phase rifting areas can experience shifts from erosion to deposition and vice-versa, which results in incomplete stratigraphic records and the reworking of sediments. 4) Syn-rift stratigraphic development show reverse trends during single and multi-phase rifting, with sedimentary facies changing from marine/lacustrine to fluvial in the singlephase rift and from fluvial to marine/lacustrine in the multi-phase rift.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
SP designed the study, performed the experiments and wrote the article. AG designed the study and wrote the article GD helped in performing the experiments RG wrote the article CP designed the study and wrote the article EF provided the data to perform the experiments.

FUNDING
This study received funding from the MultiRift project founded through the PETROMAKS 2 program of the Research Council of Norway (Project number: 215591).