Quantitative stratigraphic analysis in a source-to-sink numerical framework
The sedimentary architecture at continental margins reflects the interplay between the rate of change of accommodation creation (δA) and the rate of change of sediment supply (δS). Stratigraphic interpretation increasingly focuses on understanding the link between deposition patterns and changes in δA∕δS, with an attempt to reconstruct the contributing factors. Here, we use the landscape modelling code pyBadlands to (1) investigate the development of stratigraphic sequences in a source-to-sink context; (2) assess the respective performance of two well-established stratigraphic interpretation techniques: the trajectory analysis method and the accommodation succession method; and (3) propose quantitative stratigraphic interpretations based on those two techniques. In contrast to most stratigraphic forward models (SFMs), pyBadlands provides self-consistent sediment supply to basin margins as it simulates erosion, sediment transport and deposition in a source-to-sink context. We present a generic case of landscape evolution that takes into account periodic sea level variations and passive margin thermal subsidence over 30 million years, under uniform rainfall. A set of post-processing tools are provided to analyse the predicted stratigraphic architecture. We first reconstruct the temporal evolution of the depositional cycles and identify key stratigraphic surfaces based on observations of stratal geometries and facies relationships, which we use for comparison to stratigraphic interpretations. We then apply both the trajectory analysis and the accommodation succession methods to manually map key stratigraphic surfaces and define sequence units on the final model output. Finally, we calculate shoreline and shelf-edge trajectories, the temporal evolution of changes in relative sea level (proxy for δA) and sedimentation rate (proxy for δS) at the shoreline, and automatically produce stratigraphic interpretations. Our results suggest that the analysis of the presented model is more robust with the accommodation succession method than with the trajectory analysis method. Stratigraphic analysis based on manually extracted shoreline and shelf-edge trajectory requires calibrations of time-dependent processes such as thermal subsidence or additional constraints from stratal terminations to obtain reliable interpretations. The 3-D stratigraphic analysis of the presented model reveals small lateral variations of sequence formations. Our work provides an efficient and flexible quantitative sequence stratigraphic framework to evaluate the main drivers (climate, sea level and tectonics) controlling sedimentary architectures and investigate their respective roles in sedimentary basin development.
Since its introduction in 1970s, sequence stratigraphy has been widely used to interpret depositional architectures in terms of variations in eustatic sea level or relative sea level (i.e. accommodation) (Vail et al., 1977a; Pitman, 1978; Posamentier et al., 1988; Posamentier and Vail, 1988; Jervey, 1988). With recognition of the role of sediment supply in affecting stratal stacking patterns, the rate of change of accommodation creation (δA) versus the rate of change of sediment supply (δS) – the δA∕δS ratio – has been widely accepted as the main control of sequence formations (Schlager, 1993; Muto and Steel, 1997; Catuneanu et al., 2009; Neal and Abreu, 2009; Neal et al., 2016). The δA∕δS concept offers several advantages compared to conventional stratigraphic models as it directly relates depositional patterns to the main contributing geological drivers (e.g. eustasy, tectonics and sediment supply). Yet, the inherent difficulties in accurately describing accommodation and reconstructing sediment supply limit the quantification of the δA∕δS ratio and, as a result, the practical application of the δA∕δS concept in stratigraphic interpretations (Muto and Steel, 2000; Burgess et al., 2016).
Experimental stratigraphy (e.g. analogue experiments, stratigraphic forward modelling, SFM) plays a significant role in exploring the development of sedimentary architecture under various forcing conditions (Martin et al., 2009; Burgess et al., 2012; Muto et al., 2016). Over the past few decades, SFM has been widely used to investigate the interplay between major sequence-controlling factors (e.g. eustasy, tectonics, flexural isostasy, sediment supply, sediment compaction, basin physiography) and their influences in the formation of stratigraphic sequences (Reynolds et al., 1991; Posamentier and Allen, 1993; Steckler et al., 1993; Carvajal and Steel, 2009; Burgess et al., 2012; Granjeon et al., 2014; Csato et al., 2014; Sylvester et al., 2015; Harris et al., 2015, 2016). SFM can therefore provide insights into the quantification of links between sequence formation and the changing δA∕δS ratio. Here we use SFM as a tool to investigate the development of stratigraphic architecture under the interplay between accommodation variations and sediment supply. We use the landscape evolution code pyBadlands that describes sediment transport from source to sink in a self-consistent manner (Salles, 2016; Salles and Hardiman, 2016; Salles et al., 2018). In pyBadlands, the erosion from upstream catchments is linked to the sedimentation on basin margins through sediment routing resulting from a combination of channelling and hillslope processes. As a consequence, sediment supply to basin margins is dynamically determined without user control; it results from the interaction of imposed tectonics, climatic and eustatic variations as well as autogenic changes in upstream catchment physiography.
We then apply the rules of sequence stratigraphy to interpret predicted depositional cycles. In sequence stratigraphy, various sequence models exist and subdivide the stratal successions into unconformity- and/or correlated-conformity-bounded stratigraphic units (Galloway, 1989; Mitchum Jr., 1977; Van Wagoner et al., 1988; Helland-Hansen and Gjelberg, 1994). These models have been proved to be useful for a number of cases. However, the interpretation of systems tracts and sequences based on these sequence models can be non-objective and non-unique (Burgess et al., 2016; Burgess and Prince, 2015). Observation-based methods to interpret stratigraphic sequences have the advantage of being objective and independent of spatial and temporal scales. Hence, over the years, it has been recognized that stratigraphic interpretations should be guided by physical observations and be independent of depositional models and associated assumptions (Abreu et al., 2014; Burgess et al., 2016; Neal et al., 2016). Here, we focus on two such methods: the trajectory analysis and the accommodation succession method. Helland-Hansen and Gjelberg (1994), Helland-Hansen and Martinsen (1996), and Helland-Hansen and Hampson (2009) proposed the trajectory analysis method that correlates depositional units with the lateral and vertical migrations of the shoreline and shelf-edge trajectories. Neal and Abreu (2009) and Neal et al. (2016) proposed a refined sequence stratigraphic framework known as the accommodation succession method in which sequence sets are interpreted based on offlap break trajectory and stratal geometries. The temporal evolution of accommodation change and sediment supply can then be derived from interpreted sequence sets and key stratigraphic bounding surfaces.
This study also attempts to evaluate the performance of these two approaches to interpret the stratal architecture predicted with pyBadlands. To illustrate the workflow, we build a synthetic source-to-sink model that includes a mountain range (sediment source), an alluvial plain (sediment transfer zone) and a passive continental margin (sink area) where relatively well understood sequence-controlling drivers such as eustasy and thermal subsidence (Watts and Steckler, 1979; Bond et al., 1989; Steckler et al., 1993) are imposed. We first present the temporal evolution of predicted stratal stacking patterns and map out key stratigraphic surfaces, which serves as a reference for comparison to stratigraphic interpretations. We then follow the trajectory analysis and the accommodation succession method to interpret the predicted stratal architecture. Finally, we design a suite of numerical tools to extract of shoreline and shelf-edge trajectories, as well as accommodation change and sedimentation evolution over space and time, with the aim to integrate the trajectory analysis method and the accommodation succession method within pyBadlands to derive quantitative interpretations. These new capabilities make it possible to quickly interpret synthetic depositional cycles in a consistent manner using either the trajectory analysis or the accommodation succession method.
The workflow to build a quantitative framework of stratigraphic analysis in pyBadlands is illustrated in Fig. 1. In this study, we focus on the post-processing of model outputs. pyBadlands records the depth, relative elevation (depth relative to sea level) and thickness of each stratigraphic layer through time. With this information, we are able to extract cross sections and to reconstruct the temporal evolution of stratal stacking patterns and 3-D Wheeler diagrams at any location. The 3-D Wheeler diagram contains information of distance along the cross section, time and deposited sediment thickness. This allows us to identify key stratigraphic surfaces based on observations of stratal geometry and facies relationships.
We then interpret the synthetic depositional cycles in two ways. First, we follow the workflow proposed in the trajectory analysis and the accommodation succession method to subdivide the stratigraphic record. The trajectory analysis technique defines different trajectory classes based on the trajectories recorded at either shoreline or shelf-edge positions (Helland-Hansen and Gjelberg, 1994; Helland-Hansen and Martinsen, 1996; Helland-Hansen and Hampson, 2009). Though both represents a break in slope, the shelf edge evolves at larger spatial (Fig. 2) and over longer temporal scales than the shoreline, making it easier to define the shelf edge on seismic data. By investigating the migration direction of the shoreline and the shelf edge, four shoreline trajectory classes and three shelf-edge trajectory classes are used to characterize the successive depositional packages. The shoreline trajectory classes include the transgressive trajectory class (TTC), the ascending regressive trajectory class (ARTC), the descending regressive trajectory class (DRTC) and the stationary (i.e. non-migratory) shoreline trajectory class. The shelf-edge trajectory classes include descending trajectory class (DTC), ascending trajectory class (ATC), transgressive trajectory class (T) and stationary or flat trajectory class (Fig. 2).
Neal and Abreu (2009) and Neal et al. (2016) presented a refined hierarchy framework, known as the accommodation succession method, in which the subdivision of depositional units is entirely based on the stratal geometry resulting from the evolution of accommodation change and sediment infill. Three stacking patterns and their bounding surfaces are defined and subsequently used to assess the changing history of accommodation and sediment supply (Fig. 2). These include the retrogradation stacking (R) for , the progradation to aggradation stacking (PA) for and increasing, and the aggradation to progradation (even to degradation) stacking (AP or APD) for and decreasing (and possibly negative). The three key surfaces bounding these stacking patterns are the sequence boundary (SB), the maximum transgressive surface (MTS) and the maximum regressive surface (MRS).
We manually mark the shelf-edge (or offlap break) trajectory and stratal terminations on the final output of stratal stacking pattern reconstructed from a simulated cross section. Key stratigraphic surfaces and stacking patterns are then defined.
Second, a series of post-processing tools are implemented in pyBadlands to numerically extract the shoreline and shelf-edge positions, as well as the temporal evolution of δA and δS through time and space (Fig. 3). The shoreline position is recorded by tracking the topographic contour that corresponds to sea level. The shelf edge is calculated by assuming a critical slope of 0.025∘ in this case. Changes in relative sea level and sedimentation rate are used as proxies for δA and δS (Poag and Sevon, 1989; Galloway and Williams, 1991; Liu and Galloway, 1997; Galloway, 2001). Therefore, δA at any given location from time T1 to T2 integrates changes in eustatic sea level and basement subsidence, and δS at any given location between time T1 and T2 is given by deposited strata thickness. In this study, both δA and δS are in m Myr−1. We then extract δA and δS at shoreline positions through time. Trajectory classes, stacking patterns and stratigraphic surfaces are defined automatically based on calculated shoreline, shelf-edge trajectories, and time-dependent δA and δS.
We provide an example to illustrate how our workflow can be used to automatically generate stratigraphic sequences and analyse them in an integrated numerical toolbox.
We create an initial synthetic landscape of dimensions 300 km by 200 km with a spatial resolution of 0.5 km. The region includes a mountain range, an alluvial plain and an adjacent continental margin consisting of a continental shelf, a continental slope and an oceanic basin. Details of the geometry are presented in Fig. 4a. The model duration is 30 Myr, generating 300 stratigraphic layers with display intervals every 0.1 Myr. Our model setting mimics the first-order, long-term landscape evolution and associated stratigraphic sequence development along a passive continental margin, with forcing conditions including climate, eustatic sea level change and thermal subsidence.
In this study, we use a single flow direction, detachment-limited stream power law, and a simple downslope creep law to describe erosion, sediment transport and deposition. Equations and model parameters are provided in the Supplement. Considering that this study focuses on long-term stratigraphic evolution due to sea level change, both climate and subsidence patterns are considered constant. Climate is assumed to be directly related to precipitation with a spatially and temporally uniform precipitation rate of 2.0 m yr−1 over 30 Myr. The sediment input varies through time, depending on the dynamic evolution of source area.
Sea level fluctuations are a major driver of changes in accommodation and thus stratigraphic sequence development. They act at different temporal scales, resulting in various stratigraphic cycles with first-order cycles of duration around 200–300 Myr, second-order cycles of duration around 10–80 Myr and third-order cycles of duration around 1–10 Myr (Vail et al., 1977b). Here we consider the effects of long-term eustatic cycles and assume that eustasy is independent of climate and local tectonics. For simplicity, eustatic sea level is modelled using a sinusoidal curve consisting of three 10 Myr cycles of 50 m amplitude (Fig. 4b), which correspond to second- to third-order eustatic fluctuations (Vail et al., 1977b).
Thermal subsidence is an important process leading to the deepening of the basin floor due to isostatic adjustment during lithosphere cooling. A simple stretching model from McKenzie (1978) is applied in this study, in which subsidence is produced by thermal relaxation following an episode of extension. The equation to derive the subsidence at time t is
where , , a=125 km is the thickness of lithosphere, ρ0=3300 kg m−3 the density of the mantle at 0 ∘C, ρw=1000 kg m−3 the density of seawater, K−1 the thermal expansion coefficient for both the mantle and the crust, T1=1333 ∘C the temperature of the asthenosphere and τ=62.8 Myr the characteristic thermal diffusion time. The stretching factor β characterizes the extension of the lithosphere. These parameters are taken from McKenzie (1978). In our experiments, thermal subsidence is imposed on the continental margin, starting from the initial shoreline (Fig. 4a), which experiences the least subsidence, to the outward edge where subsidence is maximum. We take β as distance-dependent and calculate the thermal subsidence accumulated at 10, 20 and 30 Myr, respectively (Fig. 4c). The subsidence rate is constant at each single position but increases along the x axis. In our model, relative sea level is the combined result of eustatic sea level and thermal subsidence and thus varies spatially due to basement subsidence.
4.1 Temporal evolution of stratal architecture
Our post-processing tools quickly extract simulated stratal stacking patterns as well as Wheeler diagrams in any region of the simulated domain. Figure 5 presents the development of stratal stacking pattern generated along a dip-oriented cross section through the middle of the domain. The stratal architecture is coloured based on depositional environments defined using six paleo-depth windows. We infer depositional facies, changes in depositional trends, stratal terminations and stratal geometries from the temporal evolution of stratal stacking patterns. This information is then used to identify stratigraphic surfaces and to define systems tracts and sequences (Van Wagoner et al., 1988). Sediment flux is extracted along the cross section by computing the total volume of deposited sediments per unit width in 0.5 Myr intervals (Fig. 5d). As shown in Fig. 5a, an oblique prograding clinoform develops due to the falling eustatic sea level, with strata toplapped by a subaerial unconformity (SU). This subaerial unconformity terminates downdip at the offlap-break at 4.0 Myr and it transfers to a marine correlative conformity (CC*), which together forms a sequence boundary (Fig. 5a). The shoreline steps back and strata onlap the subaerial unconformity from 4.0 Myr. The prograding packages then shift to aggradational pattern until 6.5 Myr when the shelf edge reaches its most basinward location, marked by the maximum regressive surface at 6.5 Myr. Eustatic sea level first falls and then rises at a relative slow changing rate between 4.0 and 6.5 Myr, and the clinoform formed during this period is characterized by thin topsets and thick foresets. The sediment flux constantly increases from 0.7 to 2.0 km2 Myr−1 over the first 6.5 Myr (Fig. 5d), with small variations due to the lateral shifts in the position of the river mouth. Eustatic sea level rises from 6.5 to 9.0 Myr, and strata continue onlapping the subaerial unconformity and fill incised channels with fluvial sediments to form a maximum flooding surface (MFS) at 9.0 Myr. The clinoform formed during this time period is characterized by thick topsets and an absence of foresets. From 9.0 to 10.0 Myr, the shoreline steps back and strata continue onlapping while the stacking pattern changes from retrogradation to aggradation. During the phase of sea level rise from 6.5 to 10.0 Myr, the sediment flux slightly decreases to 1.6 km2 Myr−1 (Fig. 5d). Eustatic sea level falls from 10 to 13.5 Myr and strata stacking changes to progradation, forming a second subaerial unconformity (Fig. 5b). The surface formed at 10.0 Myr that separates aggradation from progradation is defined as correlative conformity (CC). Mid-slope sediments start to accumulate within the progradational clinoform between 10.0 and 13.5 Myr. The sediment flux from 10.0 to 13.5 Myr reveals a gentle increasing trend with large variations of up to 0.8 km2 Myr−1, followed by a significant drop in sediment flux at 13.5 Myr (Fig. 5d). The stratigraphic units accumulated between 4.0 and 13.5 Myr form a complete sequence (no. 2) bounded by two composite surfaces consisting of subaerial unconformities and correlative conformities.
Following the same criteria for identification of stratigraphic surfaces, a complete sequence (no. 3) is defined from 13.5 to 23.5 Myr, with the maximum regressive surface formed at 16.5 Myr, the maximum flooding surface formed at 19.0 Myr and the correlative conformity formed at 20.0 Myr (Fig. 5b and c). The sediment flux from 13.5 to 23.5 Myr shows an increasing trend with large variations of up to 1.1 km2 Myr−1, followed by a significantly drop from 2.8 to 1.5 km2 Myr−1 at 23.5 Myr (Fig. 5d). From 23.5 to 30 Myr, an incomplete sequence (no. 4) develops with a maximum regressive surface at 26.5 Myr and a maximum flooding surface at 28.5 Myr. The sediment flux during this time period shows two anomalous peaks at 24.0 and 26.5 Myr (Fig. 5d).
We also reconstructed the final stratal stacking pattern by computing stratal thickness in 100 kyr increments (Fig. 6a), which constitutes the basis for a 3-D Wheeler diagram (Fig. 6b). The 3-D Wheeler diagram shows the horizontal distribution of depositional environments and the accumulated sediment thickness along the cross section through time. Deposition hiatuses and condensed sections, which are essential to recognize bounding surfaces such as subaerial unconformities or sequence boundaries, are also denoted on the Wheeler diagram, as well as transgressive and flooding surfaces (Payton et al., 1977; Miall, 2004). The stratal thickness pattern shows three cycles, with thicker sediment accumulation during progradational and aggradational stratal stacking accompanied with a low rate of eustatic sea level change, and thinner sediment accumulation during retrogradational stratal stacking accompanied with rising eustatic sea level. The stacking of sequences no. 1 to no. 4 shows an aggradational progradation pattern, with at least 10 km seaward progradation between successive cycles. Eustatic sea level fluctuations control sequence formations, as expected, and the effect of thermal subsidence and offshore sedimentation can be recognized when taking the stratigraphic package as a whole: the aggradational stacking reveals contributions of basement subsidence in creating accommodation, while the progradational stacking reveals the contributions of offshore sedimentation in decreasing accommodation. When correlating the formation of stratigraphic surfaces with the horizontal migration of depositional packages, we find that the timing of maximum regressive surface agrees well with the onset of decreasing stratal thickness. Sequence boundaries correspond to the shift of shoreline from forestepping to backstepping, while correlative conformities correspond to the shift of shoreline from backstepping to forestepping.
4.2 Interpretation of depositional sequences
We now focus on the interpretation of the stratal architecture using both the trajectory analysis and accommodation succession methods.
4.2.1 Trajectory analysis
On the stratal stacking pattern extracted from the final output at 30 Myr, we manually pick the break in slope of the shelf-slope-scale clinoforms as shelf-edge positions (magenta dots in Fig. 7a). Shoreline positions are difficult to pick on the cross section because shoreline clinoforms are not well developed with this model setting. We therefore focus on the analysis of shelf-edge trajectory evolution. According to lateral and vertical shifts of the shelf edge through time, descending shelf-edge trajectory classes are identified from 0 to 5, 10 to 14, and 20 to 23 Myr; ascending shelf-edge trajectory classes are recognized from 5 to 6.5, 14 to 17, and from 23 to 26.5 Myr; and transgressive shelf-edge trajectory classes are defined from 6.5 to 10, 17 to 20, and 26.5 to 30 Myr (Fig. 7a).
In addition to manually picking the shelf-edge trajectory on the final output, we developed post-processing tools to extract time-dependent shelf-edge and shoreline positions from successive outputs and interpret predicted depositional cycles accordingly. Figure 7b displays the extracted lateral and vertical migrations of the shelf-edge position and interpreted shelf-edge trajectory classes. The descending shelf-edge trajectory class is identified from 0 to 5.5, 10 to 13, and 20 to 22.5 Myr; the ascending shelf-edge trajectory class is identified from 5.5 to 6.5, 13 to 16.5, and from 22.5 to 25.5 Myr; and the transgressive shelf-edge trajectory class is identified from 6.5 to 10, 16.5 to 20, and 25.5 to 30 Myr (Fig. 7b). The shelf-edge trajectory classes identified through time generally resemble the ones identified manually from the final output, with differences in the timing of transition from one class to another by up to 1.5 Myr, which is greater than the temporal resolution (0.5 Myr) used to represent stratigraphic sequences (Fig. 7a and b).
We now use the post-processing toolbox to identify changes in shoreline trajectories that are difficult to pick manually for this case. Figure 7c displays the automatically detected lateral and vertical migrations of the shoreline position and interpreted shoreline trajectory classes accordingly. The descending regressive trajectory class is identified from 0 to 4, 10 to 13.5, and 20 to 23 Myr; the transgressive trajectory class is identified from 5 and 10, 15 and 20, and 25 to 30 Myr (Fig. 7c). DRTC and TTC are separated by an erosional surface and a depositional hiatus, whereas the transition from TTC to DRTC is related to condensed stacking sections in the distal area induced by limited sediment supply (Fig. 5c). We note that between 4 and 5, 13.5 and 15, and 23 and 25 Myr, the shoreline migrates landward even though sea level is falling – we call this trajectory type the “descending transgressive trajectory class” (DTTC) as this shoreline evolution is not described in the literature. In our models, the DTTC occurs because the basement subsidence overrides the falling sea level and thus creates positive accommodation (Fig. 7e). This phenomenon is due to the model forcing conditions, and its identification directly linked to the time-dependent analysis of shoreline trajectory carried out here.
4.2.2 Accommodation succession analysis
Next, we apply the accommodation succession method to analyse the sequence formation in terms of changes in accommodation and sedimentation. Following the workflow proposed by Neal et al. (2016), we first manually marked stratal terminations (i.e. onlap, downlap, toplap) and offlap breaks on the simulated stratal stacking pattern (Step 1 in Fig. 8a). Here, the offlap break is defined as the shelf edge rather than the shoreline as shoreline-scale clinoforms do not develop with these model settings. Three stacking patterns and their bounding surfaces are then defined (Steps 2 and 3 in Fig. 8a). For example, toplaps and downlaps are observed in the first 3.5 Myr and correspond to progradational (P) stacking. The stratal geometry associated with P stacking is characterized by either erosion or by a lack of topset deposition and relatively thick clinoform front. Though strata keep downlapping, onlap terminations start replacing toplap terminations after 3.5 Myr, which reflects successive phases of progradation and aggradation. This depositional trend is defined as progradation to aggradation stacking, which is characterized by thin topset deposits and thick clinoform fronts. The unconformity formed between P and PA stacking is interpreted as a sequence boundary. Retrogradation stacking corresponds to the onlapping of stratal deposits and landward shift of offlap break around 6.5 Myr. The thick topset deposit and condensed distal stacking characterize the stratal geometry deposited during this stage. Maximum regressive surfaces are defined at the transition between PA and R stacking. At 10 Myr, the offlap break starts migrating seaward and downward. Toplap and downlap terminations are also observed after this time. The geometry of deposited strata also changes and is characterized by the formation of thicker clinoform fronts. This stacking corresponds to the AP class (aggradation to progradation). Maximum transgressive surfaces separate R stacking from AP stacking. At 13 Myr, onlap terminations are visible and the offlap break starts migrating upward, corresponding to the transition towards PA stacking just above the SB surface. Following the same criteria, successive stacking of PA, R and AP as well as bounding surfaces (SB, MTS and MRS) are defined on the cross section (Fig. 8b). Finally, the interpreted R, AP and PA stacking patterns are used to assess the changing history of δA and δS (Step 4 in Fig. 8b).
Instead of calculating δA∕δS as a ratio (Neal et al., 2016), we compute δA−δS at the shoreline over time (Fig. 8d), because δS can be equal to zero (Fig. 3). Under this approach, corresponds to R stacking and is equivalent to ; and decreasing corresponds to APD stacking and is equivalent to and decreasing; and and increasing corresponds to PA stacking and is equivalent to and increasing. The stacking patterns can then be defined automatically (Fig. 8c) and are almost identical to the manually identified ones: differences in the timing of the change in stacking pattern are 0.5 Myr at most, which is equal to the temporal resolution (0.5 Myr) with which stratigraphic sequences are represented in Figs. 7 and 8.
We compare and contrast the interpretations resulting from observations of temporal evolving stratal architecture (Sect. 4.1) with the interpretations resulting from the manual application of both trajectory analysis and accommodation succession methods (Sect. 4.2.1) and from quantitative analysis of these two methods using our post-processing tools (Sect. 4.2.2).
The manual application of shelf-edge trajectory analysis presents reliable interpretations of transgressive stratigraphic units (T) but displays notable discrepancies in separating descending stratigraphic units (DTC) from ascending stratigraphic units (ATC) especially during in the early stages of deposition (Fig. 7a). Note that here we only discuss discrepancies beyond 0.5 Myr, which is the time interval of reconstructing stratal stacking patterns in Sect. 4.2. The imposed thermal subsidence modifies the position of the shelf-edge positions after its formation (Fig. 9): the shelf-edge trajectory appears to rise between 3.5 and 6.5 Myr on the snapshot at 10 Myr (Fig. 9a), whereas it appears to fall between 3.5 and 5 Myr on the snapshot at 20 Myr (Fig. 9b), because of ongoing thermal subsidence between 10 and 20 Myr. As a consequence, identifying strata on the final output artificially extends the duration of descending trajectory class (Figs. 7a and 9b). This should be kept in mind when identifying shelf-edge trajectories on seismic profiles, which are by nature a snapshot of the evolution of a basin. Constraining time-dependent processes such as thermal subsidence from sedimentary packages would be useful to correct shoreline and shelf-edge trajectories for these processes.
The quantitative shelf-edge trajectory analysis reveals discrepancies in defining ascending trajectory classes at 5.5, 22.5 and 25.5 Myr (Fig. 7b), as extracting the position of the shelf edge is more uncertain for steep strata. The tool we have developed can be applied to actual sections as long as seismic timelines are accessible. Again, we emphasize that additional constraints from observation of stratal geometry would improve the interpretations of stratigraphic units. In terms of quantitative shoreline trajectory analysis, the distinction between the descending transgressive trajectory class and the transgressive trajectory class is controlled by the eustatic sea level fluctuations (Fig. 7c). Since the shoreline and the shelf edge represent the break in slope of clinoforms at different scales (Patruno et al., 2015, Fig. 2), it is not appropriate to apply shoreline trajectory analysis to shelf-slope clinoforms. However, the numerical tools we provided to extract time-dependent shoreline positions based on a given sea level forcing would also be useful and applicable to short-term numerical and analogue experiments (Martin et al., 2009; Granjeon et al., 2014).
The stratigraphic interpretations from the accommodation succession method indicate that there is no significant difference between analysing the final output and time-dependent outputs (Figs. 8b and 5). Therefore, the analysis of the presented model is more robust with the accommodation succession method than with the trajectory analysis method.
The quantitative accommodation succession analysis also shows largely consistent interpretations (Figs. 8c and 5), except for a 1 Myr discrepancy in demarcating aggradation to progradation to degradation stacking and progradation to aggradation stacking at 3 Myr. We find that the δA−δS curve (Fig. 8d) presents trends similar to the rate of eustatic sea level change (Fig. 4b). This suggests that the evolution of δA−δS is a proxy for the derivative of sea level change with respect to time rather than a direct proxy for sea level change. However, a difference exists between the δA−δS curve and the rate-of-sea-level-change curve: the δA−δS curve shows asymmetrical fluctuations with larger amplitude below zero than above zero, which is attributed to the sediment supply. Discrepancies of <0.5 Myr are observed between the δA−δS curve and the rate of eustatic sea level change curve, which are likely to be related to the temporal resolution (0.5 Myr) used to compute δA−δS.
A common issue when calculating the ratio δA∕δS is the lack of unique approaches and common dimensional units to define these two metrics (Muto and Steel, 2000; Burgess, 2016). Both metrics represent changes in volume over a specific time interval and thus should be defined in cubic metres per year (m3 yr−1). However, difficulties in delineating the potential space for sediment deposition require the use of proxies to quantify accommodation. Although the sedimentation rate is often used as a proxy for δS (Poag and Sevon, 1989; Galloway and Williams, 1991; Liu and Galloway, 1997; Galloway, 2001), it only provides information about the deposition at a single location and does not reflect the spatial distribution of sedimentation (Petter et al., 2013). Furthermore, the distribution of sediment deposition is not determined solely by sediment supply but is a combined result of the distance to sediment source, basin physiography and sediment transport efficiency (Posamentier et al., 1992; Posamentier and Allen, 1993). Recently, new methods have been proposed to improve the quantification of δS. Petter et al. (2013) proposed an approach that directly reconstructs sediment paleo-fluxes from stratigraphic records. Ainsworth et al. (2018) used a technique termed “TSF analysis” in which parasequence thickness (T) is used as proxy for accommodation at the time of deposition while parasequence sandstone fraction (SF) is used as proxy for sediment supply. Our work could be used to integrate and test these new quantitative interpretations based on the evolution of accommodation and sedimentation in a stratigraphic modelling framework. The quantification of δA and δS presented here could be extended in future work to investigate the interplay between accommodation change and sediment supply in 3-D.
Our source-to-sink numerical scheme generates 3-D landscape evolution and stratigraphic development, though only 2-D stratigraphic architectures are extracted along dip-oriented cross sections. To evaluate the lateral variations in sequence formation potentially induced by sediment diffusion in offshore environment and lateral migrations of river mouth, we extract five dip-oriented cross sections (CS1 to CS5) that are parallel to each other and two along-strike cross sections (Fig. 10). Cross section CS3 is the same one as presented in the result section. The timing of key stratigraphic surfaces on CS1 to CS5 is shown in Table 1, with differences varying from 0.5 to 1.5 Myr. The timing of sequence boundaries shows the most variations, compared to other stratigraphic surfaces. The timings of correlative conformities of sequence no. 2 (CC1) and sequence no. 3 (CC2) are consistent on cross sections CS1 to CS5 and correspond to the onset of eustatic sea level fall. The stacking of depositional environments along strike shows increasing variations towards the basin. For the presented case there is overall little variation in stratigraphic sequences across strike and along strike, which is expected from the model setup. The presented tools can be used for the 3-D stratigraphic analysis of more complex cases.
We have explored stratigraphic development in a source-to-sink context in which the dynamic sediment supply to the passive continental margin depends on climatic and tectonic evolution in the source area. Therefore, the physiographic evolution of the upstream region controls the depositional patterns in the sink area together with accommodation change (Ruetenik et al., 2016; Li et al., 2018). Most stratigraphic forward models focus on simulating sediment transport and deposition in the sink area, which limits the interpretation of the effect climatic and tectonic evolution on the stratigraphic record. In our framework, sediment transport and supply to the margin is dynamically related to autogenic and allogenic processes. As an example, though forced with uniform rainfall pattern in the source region, the rate of sediment supply to the sink area fluctuates with time (Fig. S1 in the Supplement). This highlights the complex relationships between erosional signals and the preservation of a depositional record (Van Heijst et al., 2001; Kim et al., 2006; Kim, 2009). The source-to-sink numerical scheme used in pyBadlands makes it possible to explore important questions for the future of sequence stratigraphy, such as the role of sediment supply variations in the generation of stratigraphic sequences at different temporal scales (Burgess, 2016), as well as the importance of allogenic and autogenic processes in the formation and evolution of stratal record (Paola, 2000; Paola et al., 2009).
We modelled the formation of stratigraphic sequences over 30 Myr, which represents second- to third-order stratigraphic cycles (Vail et al., 1977b). Over this temporal scale, long-term eustatic sea level changes and dynamic uplift or subsidence induced by tectonics or deep-Earth processes (e.g. mantle flow-driven dynamic topography) might drive deposition (Burgess and Gurnis, 1995; Burgess et al., 1997), moderated by higher-frequency fluctuations in climate and sea level. The long-term stratigraphic record along continental passive margins thus potentially contains important constraints on the evolution of the structure of the deep Earth (Mountain et al., 2007; Braun, 2010; Flament et al., 2013; Salles et al., 2017). The framework we have introduced in this study integrates both long-term surface processes and stratigraphic modelling and can be used to quantitatively investigate the influences of long-term tectonics and deep-Earth dynamics on stratal geometries and depositional pattern evolution as well as their feedback mechanisms (Jordan and Flemings, 1991; van der Beek et al., 1995; Rouby et al., 2013). Note that the tools we have introduced here are not specific to any temporal or spatial scale and can also be used for short-term stratigraphic modelling.
We used pyBadlands to model sediment erosion, transport and deposition from source to sink, as well as to investigate the formation of stratigraphic sequences on a passive continental margin in response to long-term sea level change, thermal subsidence and dynamic sediment supply. We analysed the predicted stratigraphic architecture based on observations of shelf-edge or offlap break trajectory, stratal terminations and stratal geometries, following the workflow of two stratigraphic interpretation approaches: the trajectory analysis and the accommodation succession method. We introduced a set of post-processing tools to extract the temporal evolution of shoreline, shelf edge, rate of accommodation change (δA) and sedimentation (δS), based on which automatic interpretations can be obtained. Our results suggest that the stacking patterns defined with the accommodation succession method provide more robust reconstructions of the changing history of accommodation and sedimentation than the trajectory analysis method, because the interpretation of sequences with the trajectory analysis depends on time. In contrast, the accommodation succession method is not affected by the time dependence of processes controlling the evolution of deposition. As a consequence, seismic data should be backstripped before stratigraphic sequences are interpreted using the trajectory analysis method. Our work presents an integrated workflow that can be used to generate 2-D and 3-D stratal architectures on basin margins and to interpret stratigraphic sequences produced by large-scale and long-term numerical experiments.
pyBadlands is an open-source package distributed under the GNU GPLv3. The source code is available on GitHub (https://github.com/badlands-model, last access: 18 June 2019). The easiest way to use the code is via our Docker container (https://github.com/badlands-model/pyBadlands-Docker-Serial, last access: 18 June 2019, pyBadlands-serial), which is shipped with the complete list of dependencies, the model companion and a series of examples. The code, inputs and post-processing functions used in this study are available on GitHub (https://github.com/XuesongDing/GMD-models, last access: 18 June 2019).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd-12-2571-2019-supplement.
XD designed the experiments and analysed the outputs with all co-authors. TS developed the model code. XD prepared the manuscript with contributions from all co-authors.
The authors declare that they have no conflict of interest.
Xuesong Ding, Tristan Salles and Patrice Rey were supported by Australian Research Council grant IH130200012 (Basin GENESIS Hub), and Nicolas Flament was supported by Australian Research Council grant DE160101020. Ding acknowledges supports from the International Association for Mathematical Geosciences (IAMG) grant. This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government. We thank Zoltan Sylvester and Jack Neal for their constructive comments on this paper. We also thank Kyle Straub, Cornel Olariu and Jack Neal for their valuable comments on an earlier version of the manuscript.
This research has been supported by the Australian Research Council (grant nos. IH130200012 and DE160101020).
This paper was edited by Thomas Poulet and reviewed by Zoltan Sylvester and Jack Neal.
Abreu, V., Pederson, K., Neal, J., and Bohacs, K.: A simplified guide for sequence stratigraphy: nomenclature, definitions and method, in: William Smith Meeting, 22–23 September 2014, The Geological Society, Burlington House, London, Paper number: 4, 2014. a
Ainsworth, B. R., McArthur, J. B., Lang, S. C., and Vonk, A. J.: Quantitative sequence stratigraphy, AAPG Bull., 102, 1913–1939, 2018. a
Bond, G. C., Kominz, M. A., Steckler, M. S., and Grotzinger, J. P.: Role of thermal subsidence, flexure, and eustasy in the evolution of early Paleozoic passive-margin carbonate platforms, in: Controls on Carbonate Platform and Basin Development, SEPM, Special Publication, 44, 39–61, 1989. a
Braun, J.: The many surface expressions of mantle dynamics, Nat. Geosci., 3, 825–833, 2010. a
Burgess, P. M. and Gurnis, M.: Mechanisms for the formation of cratonic stratigraphic sequences, Earth Planet. Sc. Lett., 136, 647–663, 1995. a
Burgess, P. M. and Prince, G. D.: Non-unique stratal geometries: implications for sequence stratigraphic interpretations, Basin Res., 27, 351–365, 2015. a
Burgess, P. M., Gurnis, M., and Moresi, L.: Formation of sequences in the cratonic interior of North America by interaction between mantle, eustatic, and stratigraphic processes, Geol. Soc. Am. Bull., 109, 1515–1535, 1997. a
Burgess, P. M., Roberts, D., and Bally, A.: A brief review of developments in stratigraphic forward modelling, 2000–2009, Regional Geology and Tectonics: Principles of Geologic Analysis, Chapter 14, 379–404, 2012. a, b
Carvajal, C. and Steel, R.: Shelf-edge architecture and bypass of sand to deep water: influence of shelf-edge processes, sea level, and sediment supply, J. Sediment. Res., 79, 652–672, 2009. a
Catuneanu, O., Abreu, V., Bhattacharya, J., Blum, M., Dalrymple, R., Eriksson, P., Fielding, C. R., Fisher, W., Galloway, W., Gibling, M., Giles, K. A., Holbrook, J. M., Jordan, R., Kendall C. G. St. C., Macurda, B., Martinsen, O. J., Miall, A. D., Neal, J. E., Nummedal, D., Pomar, L., Posamentier, H. W., Pratt, B. R., Sarg, J. F., Shanley, K. W., Steel, R. J., Strasser, A., Tucker, M. E., and Winker, C.: Towards the standardization of sequence stratigraphy, Earth-Sci. Rev., 92, 1–33, 2009. a
Csato, I., Catuneanu, O., and Granjeon, D.: Millennial-scale sequence stratigraphy: numerical simulation with Dionisos, J. Sediment. Res., 84, 394–406, 2014. a
Flament, N., Gurnis, M., and Müller, R. D.: A review of observations and models of dynamic topography, Lithosphere, 5, 189–210, 2013. a
Galloway, W. E.: Genetic stratigraphic sequences in basin analysis I: architecture and genesis of flooding-surface bounded depositional units, AAPG Bull., 73, 125–142, 1989. a
Galloway, W. E. and Williams, T. A.: Sediment accumulation rates in time and space: Paleogene genetic stratigraphic sequences of the northwestern Gulf of Mexico basin, Geology, 19, 986–989, 1991. a, b
Granjeon, D., Martinius, A., Ravnas, R., Howell, J., Steel, R., and Wonham, J.: 3D forward modelling of the impact of sediment transport and base level cycles on continental margins and incised valleys, Depositional Systems to Sedimentary Successions on the Norwegian Continental Margin: International Association of Sedimentologists, Special Publication, 46, 453–472, 2014. a, b
Harris, A. D., Carvajal, C., Covault, J., Perlmutter, M., and Sun, T.: Numerical stratigraphic modeling of climatic controls on basin-scale sedimentation, in: AAPG Annual Convention and Exhibition, 31 May–3 June 2015, Denver, Colorado, USA, 2015. a
Harris, A. D., Covault, J. A., Madof, A. S., Sun, T., Sylvester, Z., and Granjeon, D.: Three-dimensional numerical modeling of eustatic control on continental-margin sand distribution, J. Sediment. Res., 86, 1434–1443, 2016. a
Jervey, M.: Quantitative geological modeling of siliciclastic rock sequences and their seismic expression, Spec. Publ. SEPM, 42, 47–69, 1988. a
Jordan, T. and Flemings, P.: Large-scale stratigraphic architecture, eustatic variation, and unsteady tectonism: A theoretical evaluation, J. Geophys. Res., 96, 6681–6699, 1991. a
Kim, W.: Decoupling allogenic and autogenic processes in experimental stratigraphy, AGU Fall Meeting Abstracts, EP53A-0609, 2009. a
Kim, W., Paola, C., Voller, V. R., and Swenson, J. B.: Experimental measurement of the relative importance of controls on shoreline migration, J. Sediment. Res., 76, 270–283, 2006. a
Li, Q., Gasparini, N. M., and Straub, K. M.: Some signals are not the same as they appear: How do erosional landscapes transform tectonic history into sediment flux records?, Geology, 46, 407–410, https://doi.org/10.1130/G40026.1, 2018. a
Martin, J., Paola, C., Abreu, V., Neal, J., and Sheets, B.: Sequence stratigraphy of experimental strata under known conditions of differential subsidence and variable base level, AAPG Bull., 93, 503–533, 2009. a, b
Miall, A. D.: Empiricism and model building in stratigraphy: the historical roots of present-day practices, Stratigraphy, 1, 3–25, 2004. a
Mitchum Jr., R. M.: Seismic stratigraphy and global changes of sea level: Part 11. Glossary of terms used in seismic stratigraphy: Section 2. Application of seismic reflection configuration to stratigraphic interpretation, Mem. Amer. Assoc. Petrol. Geol., 26, 205–212, 1977. a
Mountain, G. S., Burger, R. L., Delius, H., Fulthorpe, C. S., Austin, J. A., Goldberg, D. S., Steckler, M. S., McHugh, C. M., Miller, K. G., Monteverde, D. H., Orange, D. L., and Pratson, L. F.: The long-term stratigraphic record on continental margins, Continental margin sedimentation: From sediment transport to sequence stratigraphy: International Association of Sedimentologists Special Publication, 37, 381–458, 2007. a
Muto, T. and Steel, R. J.: Principles of regression and transgression: the nature of the interplay between accommodation and sediment supply: perspectives, J. Sediment. Res., 67, 994–1000, 1997. a
Muto, T., Steel, R. J., and Burgess, P. M.: Contributions to sequence stratigraphy from analogue and numerical experiments, J. Geol. Soc. London, 173, 837–844, 2016. a
Neal, J. E., Abreu, V., Bohacs, K. M., Feldman, H. R., and Pederson, K. H.: Accommodation succession (δA∕δS) sequence stratigraphy: observational method, utility and insights into sequence boundary formation, J. Geol. Soc. London, 173, 803–816, 2016. a, b, c, d, e, f, g, h
Paola, C.: Quantitative models of sedimentary basin filling, Sedimentology, 47, 121–178, 2000. a
Paola, C., Straub, K., Mohrig, D., and Reinhardt, L.: The “unreasonable effectiveness” of stratigraphic and geomorphic experiments, Earth-Sci. Rev., 97, 1–43, 2009. a
Patruno, S., Hampson, G. J., and Jackson, C. A.: Quantitative characterisation of deltaic and subaqueous clinoforms, Earth-Sci. Rev., 142, 79–119, 2015. a
Payton, C. E. et al.: Seismic stratigraphy: applications to hydrocarbon exploration, American Association of Petroleum Geologists Tulsa, OK, 26, 1–516, 1977. a
Petter, A. L., Steel, R. J., Mohrig, D., Kim, W., and Carvajal, C.: Estimation of the paleoflux of terrestrial-derived solids across ancient basin margins using the stratigraphic record, GSA Bulletin, 125, 578–593, 2013. a, b
Pitman, W. C.: Relationship between eustacy and stratigraphic sequences of passive margins, GSA Bulletin, 89, 1389–1403, 1978. a
Poag, C. W. and Sevon, W. D.: A record of Appalachian denudation in postrift Mesozoic and Cenozoic sedimentary deposits of the US middle Atlantic continental margin, Geomorphology, 2, 119–157, 1989. a, b
Posamentier, H. and Vail, P.: Eustatic controls on clastic deposition II – Sequence and systems tract models, in: Sea-Level Changes: an Integrated Approach, edited by: Wilgus, C. K., Hastings, B. S., Posamentier, H., Van Wagoner, J., Ross, C. A., and St. C. Kendall, C. G., Spec. Publ. Soc. Econ. Paleont. Miner., 42, 125–154, 1988. a
Posamentier, H., Jervey, M., and Vail, P.: Eustatic controls on clastic deposition I – Conceptual frameworks, in: Sea-Level Changes: an Integrated Approach, edited by: Wilgus, C. K., Hastings, B. S., St. C. Kendall, C. G., Posamentier, H. W., Ross, C. A., and Van Wagoner, J. C., Spec. Publ. Soc. Econ. Paleont. Miner., 42, 109–124, 1988. a
Posamentier, H. W., Allen, G. P., James, D. P., and Tesson, M.: Forced regressions in a sequence stratigraphic framework: concepts, examples, and exploration significance (1), AAPG bulletin, 76, 1687–1709, 1992. a
Reynolds, D. J., Steckler, M. S., and Coakley, B. J.: The role of the sediment load in sequence stratigraphy: The influence of flexural isostasy and compaction, J. Geophys. Res., 96, 6931–6949, 1991. a
Rouby, D., Braun, J., Robin, C., Dauteuil, O., and Deschamps, F.: Long-term stratigraphic evolution of Atlantic-type passive margins: A numerical approach of interactions between surface processes, flexural isostasy and 3D thermal subsidence, Tectonophysics, 604, 83–103, 2013. a
Salles, T.: Badlands: A parallel basin and landscape dynamics model, SoftwareX, 5, 195–202, 2016. a
Salles, T. and Hardiman, L.: Badlands: An open-source, flexible and parallel framework to study landscape dynamics, Comput. Geosci., 91, 77–89, 2016. a
Salles, T., Flament, N., and Müller, D.: Influence of mantle flow on the drainage of eastern Australia since the Jurassic Period, Geochem. Geophy. Geosy., 18, 280–305, 2017. a
Salles, T., Ding, X., and Brocard, G.: pyBadlands: A framework to simulate sediment transport, landscape dynamics and basin stratigraphic evolution through space and time, PloS one, 13, e0195557, https://doi.org/10.1371/journal.pone.0195557, 2018. a
Schlager, W.: Accommodation and supply – a dual control on stratigraphic sequences, Sediment. Geol., 86, 111–136, 1993. a
Steckler, M., Reynolds, D., Coakley, B., Swift, B., and Jarrard, R.: Modelling passive margin sequence stratigraphy, in: Sequence stratigraphy and facies associations, edited by: Posamentier, H. W., Summerhayes, C. P., Haq, B. U., and Allen, G. P., International Association of Sedimentologists Special Publication 18, Blackwell Sientific Publications, Oxford, UK, 19–41, 1993. a, b
Sylvester, Z., Cantelli, A., and Pirmez, C.: Stratigraphic evolution of intraslope minibasins: Insights from surface-based model, AAPG Bulletin, 99, 1099–1129, 2015. a
Vail, P. R., Mitchum Jr., R., and Thompson III, S.: Seismic stratigraphy and global changes of sea level: Part 3. Relative changes of sea level from coastal onlap: section 2. Application of seismic reflection configuration to stratigrapic interpretation, Mem. Amer. Assoc. Petrol. Geol., 26, 63–81, 1977a. a
Vail, P. R., Mitchum Jr., R., and Thompson III, S.: Seismic stratigraphy and global changes of sea level, Part 4: Global cycles of relative changes of sea level, Mem. Amer. Assoc. Petrol. Geol., 26, 83–98, 1977b. a, b, c
van der Beek, P., Andriessen, P., and Cloetingh, S.: Morphotectonic evolution of rifted continental margins: Inferences from a coupled tectonic-surface processes model and fission track thermochronology, Tectonics, 14, 406–421, 1995. a
Van Heijst, M., Postma, G., Meijer, X., Snow, J., and Anderson, J.: Quantitative analogue flume-model study of river-shelf systems: principles and verification exemplified by the Late Quaternary Colorado river-delta evolution, Basin Res., 13, 243–268, 2001. a
Van Wagoner, J., Posamentier, H., Mitchum, R., Vail, P., Sarg, J., Loutit, T., and Hardenbol, J.: An overview of the fundamentals of sequence stratigraphy and key definitions, SEPM Society for Sedimentary Geology, 42, 39-46, 1988. a, b
Watts, A. and Steckler, M.: Subsidence and eustasy at the continental margin of eastern North America, in: Deep Drilling Results in the Atlantic Ocean: Continental Margins and Paleoenvironment, edited by: Talwani, M., Hay, W., and Ryan, W. B., AGU, 3, 218–234, 1979. a