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.

Figure 1Workflow to build the stratigraphic architecture in pyBadlands. We present three different ways of stratigraphic interpretation. The new post-processing workflows designed to automatically interpret stratigraphic sequences are shown in the blue boxes.

Download

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 $\mathit{\delta }A/\mathit{\delta }S>\mathrm{1}$, the progradation to aggradation stacking (PA) for $\mathit{\delta }A/\mathit{\delta }S<\mathrm{1}$ and increasing, and the aggradation to progradation (even to degradation) stacking (AP or APD) for $\mathit{\delta }A/\mathit{\delta }S<\mathrm{1}$ 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).

Figure 2Stratigraphic sequence interpretation approaches used in this study. In the trajectory analysis method, four shoreline trajectory classes and three shelf-edge trajectory classes are delineated based on the lateral and vertical migrations of shoreline and shelf edge. The shoreline trajectory classes include descending regressive trajectory class (DRTC), ascending regressive trajectory class (ARTC), transgressive trajectory class (TTC) and stationary trajectory class. The shelf-edge trajectory classes include descending (DTC), ascending (ATC), transgressive (T) and stationary shelf-edge trajectory classes (Helland-Hansen and Gjelberg, 1994; Helland-Hansen and Martinsen, 1996; Helland-Hansen and Hampson, 2009). In the accommodation succession method, three types of stacking patterns and their bounding surfaces are defined based on observations of stratal terminations (e.g. onlap, offlap) and stratal geometries, including aggradation to progradation to degradation (APD) stacking, progradation to aggradation (PA) stacking and retrogradation (R) stacking. The bounding surfaces are sequence boundary (SB), maximum transgressive surface (MTS) and maximum regressive surface (MRS). Each stacking pattern reflects changes in the rate of accommodation creation (δA) with respect to the rate of sediment supply (δS) at the shoreline. APD stacking corresponds to $\mathit{\delta }A/\mathit{\delta }S<\mathrm{1}$, with a decreasing trend that can be negative; PA stacking represents $\mathit{\delta }A/\mathit{\delta }S<\mathrm{1}$ and increasing; finally, R stacking occurs for $\mathit{\delta }A/\mathit{\delta }S>\mathrm{1}$ (Neal and Abreu, 2009; Neal et al., 2016).

Download

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⁻¹. 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.

Figure 3Schematic diagram showing the sedimentation from time 1 (T1) to time 2 (T2) under sea level variation and basement subsidence. The depth, relative elevation and thickness of stratigraphic layer are recorded at every time step (dT= T2 − T1). Post-processing tools extract the shoreline and shelf-edge positions and calculate the rate of accommodation creation (δA) and the rate of sedimentation (δS) through time. The shoreline position is recorded by tracing the topographic contour that corresponds to sea level. The shelf edge is calculated by assuming a critical slope 0.025^∘. δA is calculated as changes in relative sea level at shoreline position over (T2 − T1): (sea level change +  subsidence) ∕ (T2 − T1); δS is calculated as deposited strata thickness at shoreline position over (T2 − T1): (strata thickness) ∕ (T2 − T1).

Download

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.

Figure 4Model setup. (a) The initial surface elevation of the mountain region ranges from 200 to 1000 m over a width of 100 km, while the elevation of the alluvial plain ranges from 0 to 200 m over a width of 50 km. The sink area includes a continental margin in which the elevation of the continental shelf ranges from −250 to 0 m over a width of 80 km, the elevation of the continental slope ranges from −1000 to −250 m over a width of 40 km and a flat oceanic basin whose depth is −1000 m over a width of 30 km. Black lines are isobaths with an interval of 100 m. The initial elevations of shoreline and shelf edge are 0 and −200 m. (b) Eustatic sea level and its rate of change over 30 Myr. The eustatic sea level scenario modelled using a sinusoidal curve consisting of three 10 Myr cycles of 50 m peak-to-peak amplitude. (c) Distance-dependent stretching factor β and the resulting thermal subsidence at 10, 20 and 30 Myr across the continental margin.

Download

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⁻¹ 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 $E_{\mathrm{0}}=\mathrm{4}a{\mathit{\rho }}_{\mathrm{0}}\mathit{\alpha }{T}_{\mathrm{1}}/{\mathit{\pi }}^{\mathrm{2}}({\mathit{\rho }}_{\mathrm{0}}-{\mathit{\rho }}_{\mathrm{w}})$, $r=\frac{\mathit{\beta }}{\mathit{\pi }}\sin \left(\frac{\mathit{\pi }}{\mathit{\beta }}\right)$, a=125 km is the thickness of lithosphere, ρ₀=3300 kg m⁻³ the density of the mantle at 0 ^∘C, ρ_w=1000 kg m⁻³ the density of seawater, $\mathit{\alpha }=\mathrm{3.28}\times {\mathrm{10}}^{-\mathrm{5}}$ K⁻¹ the thermal expansion coefficient for both the mantle and the crust, T₁=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.

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.

Figure 9Stratal stacking pattern extracted at 10 Myr (a) and 20 Myr (b). Between 3.5 and 5 Myr, the shelf-edge trajectory changes from rising at 10 Myr to falling at 20 Myr, as a result of basement subsidence. The descending trajectory class thus extends to 5 Myr.

Download

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.

Figure 10Stratal stacking patterns on five dip-oriented cross sections (CS1 to CS5) and two along-strike cross sections. Key stratigraphic surfaces on CS1 to CS5 are identified and coloured accordingly.

Download

Table 1Timing of stratigraphic surfaces on CS1 to CS5 (Myr).

Download Print Version | Download XLSX

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 (m³ yr⁻¹). 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.

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.

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.