Numerical Modelling to Evaluate Sedimentation Effects on Heat Flow and Subsidence during Continental Rifting

Abstract

Sedimentation impacts thermal and subsidence evolution in continental rifting. Estimating the blanketing effect of sediments is crucial to reconstructing the heat flow during rifting. The sedimentary load affects the basin subsidence rate. Numerical investigation of these effects requires active and complex simulations of the thermal structure, lithospheric stretching, and sedimentation. In this paper, we introduce a numerical model to quantify these effects, which was developed using the COMSOL Multiphysics^® simulation software. Our numerical setting for the analytical and numerical solutions of thermal structure and subsidence is based on previous continental rifting studies. In our model, we accumulate a column of 5 m thick sediment layers with varied stretching factors and sedimentation rates, spanning the syn-rift to early post-rift phases over a period of 12 myr. Our results provide intuitive models to understand these sedimentation effects. The models show that an increase in sedimentation thickness significantly decreases surface heat flow, leading to lower geothermal temperature, and amplifies the subsidence magnitude. The findings also demonstrate that increases in the stretching factor and sedimentation rate enhance the blanketing effect and subsidence rate. Based on these results, we discuss key outcomes for geological applications and the possible limitations of our approach.

1. Introduction

Continental rifting is an important geological process described in rifted continental margins and sedimentary basins such as the Rio Grande Rift, East African Rift, and North Sea Basin (e.g., [1,2,3,4,5,6,7,8,9]). Fundamentally, it commences by stretching/thinning of the continental lithosphere and rising of the mantle asthenosphere close to the surface (e.g., pure shear model in Figure 1a). Despite debates about the relative importance of active (mantle convection-driven) versus passive (lithospheric-driven) processes and the evolution models of rifted continental margins, it is widely recognized that continental rifting eventually evolves into either seafloor spreading or failed (fossil) rifts (e.g., [10,11,12,13,14,15]). Rifting mechanisms have been investigated in many studies on geological and geophysical exploration, in addition to numerical and analogue modelling (e.g., [16,17,18,19,20,21,22]). Because rifting causes structural deformation and subsidence, accompanied by half-grabens and normal faults, extensional basins commonly develop on rifted margins [2,5,15,23]. For many years, sedimentary basins have been intensively studied in various geoscientific fields and in exploration projects for commercial interests (e.g., hydrocarbon, economic mineral deposits, geothermal energy, and CO₂ capture and sequestration) [10,24,25].

Sedimentation records in continental rifting are crucial for reconstructing the spatial–temporal evolution of the system (e.g., [26,27,28,29,30,31,32]). Rifting is affected by sedimentation, particularly during the evolution of the thermal structure and subsidence. The thermal effect of depositing cold sediments on top of the lithosphere is termed the thermal blanketing effect, which leads to a depression in the thermal gradient that extends down into the lithosphere [33,34,35,36,37]. The sedimentary load and compaction affect the rate and amount of subsidence in sedimentary basins [6,22,25,38,39]. However, their impact on the thermal gradient has often been underestimated in numerical geodynamic and basin modelling, and quantifying the effects is required. A major challenge is constructing an accurate numerical setting of continental rifting to examine and calculate the sedimentation effect. Therefore, we developed a numerical model using simulation software, COMSOL Multiphysics^®, to evaluate sedimentation effects on the evolution of continental rifting. Analytical and numerical benchmark solutions for the thermal structure and subsidence in the continental rifting system [40,41,42] are well established in our model domain. In the proposed model, constant stacking of sediment layers is applied from the syn-rift to early post-rift phases. We present principal results of surface heat flow, geothermal gradient, and subsidence models with variations in the stretching factor (β) and sedimentation rate. Previous models from stationary and non-sedimentation states are compared to our results. For future geological applications, we discuss requirements and limitations introduced in our numerical modelling.

Geosciences 10 00451 g001 550

2. Numerical Solutions 2.1. Numerical Domain Setup

The geological setting of continental rifting was reconstructed as a two-dimensional model using COMSOL Multiphysics^® (https://www.comsol.com/comsol-multiphysics). COMSOL Multiphysics^® is a numerical modelling software based on the finite element method (FEM), which provides functions ranging from designing numerical models in dimensional space to simulation and post-processing. Since the program supports a graphical user interface (GUI) and provides example models that enable users to conduct efficient numerical modelling, it has been used in a variety of scientific and engineering fields. We developed a numerical model domain consisting of five units: air, sediments, upper lithosphere (crust), lower lithosphere, and asthenospheric mantle. During the syn-rift phase, the upper and lower lithosphere were gradually thinned to the thicknesses determined by the givenβ (Figure 1b). The surface and bottom temperatures of the lithosphere, including crust, were set as 0 °C and 1333 °C, respectively. The domain consisted of rectangular elements, 2.5 m wide by 125,000 m thick (Figure 2a). To model thermal evolution occurring along the depth axis, a symmetric right-hand side of the setting was estimated to reduce numerical costs. An open boundary was applied to the right wall, and the left-hand side was described as a reflective wall. We implemented our model using a time-dependent direct solver with the generalized-alpha time-stepping method [43]. The initial domain setting was modelled by the steady-state geotherm. The lithosphere and asthenospheric mantle units were assumed as incompressible fluids (i.e., Boussinesq Approximation), and heat production by viscous dissipation and radiogenic isotope decay were neglected (e.g., [33,36]). The domain size, elements, solver, and time-stepping method were uniformly applied to the model in this study. Input parameters are shown in Table 1.

We stacked discrete sediment layers in the model domain to evaluate the sedimentation effect (Figure 2b). The sediment layers (5 m thickness; 0 °C temperature, as a deposition of cold sediments) were stacked by certain time increments to control the sedimentation rate and thickness; the time increment of 0.02 myr corresponds to the sedimentation rate of 250 m/myr (0.02 myr into 5 m). This process allowed us to establish a constant element size in the model domain, which is useful to calculate temperature and depth at each stacking stage (each time increment). The layers were stacked continuously in the predefined domain regarded as an air unit (3000 m thick accommodation). During the syn-rift phase, the velocity field was stretched into the sediment layers. The bottom and top of the stacked sediment layers correspond to the basin basement and land surface (or seafloor) in the sedimentary basin profile.

2.2. Temperature and Subsidence Calculations

We developed our numerical model of continental rifting based on the concepts outlined in [40,41]. Calculations and analytical solutions were previously examined by [42]. We accumulated sediments on top of the lithosphere in thermally steady, stretching, and cooling states. Our modelling was refined to provide points to compare previous numerical modelling studies. We note that modelling in this study is aimed at exhibiting thermal blanketing and subsidence by sedimentation in the most basic setting.

The fundamental models by [40] suggested using analytical solutions with a stretching factor (β ) to investigate post-rift thermal structures and subsidence in the extensional basin. Their initial subsidence model under the assumption of instantaneous stretching is a good approximation to the subsidence caused by finite extension rates. As a follow-up study, G.T. Jarvis and D.P. Mckenzie [41] improved the model to establish temporal evolution from stretching to thermal re-equilibration of the thinned lithosphere. They assumed an extensional basin is created by horizontal stretching of the lithosphere (L), with horizontal extension and vertical shortening determined byβ . The uniform stretching concept is defined as pure shear extension (Figure 1a). The thinned lithosphere is composed of a brittle deformed upper layer, manifested by faulting and fault-block rotation, and a ductile deformed (thinned) lower layer [10,14].

During stretching (syn-rift phase), the initial width and thickness (_yL) of the lithosphere are changed to_yL×βfor width and_yL/β for thickness (Figure 1b). Lithospheric thinning raises the top of the asthenospheric mantle (A) with temperature (_T1) to a depth of_yL/β, which changes the temperature profile and geothermal gradient of the lithosphere. The geotherm is evaluated by the following equations:

_T∆t = _T1×y×β_yL(y ≤ _yLβ)

_T∆t = _T1(y ≥ _yLβ)

where_T∆tis the geotherm at rifting duration,∆t; y is the depth. The temperature is calculated by solving the heat equation, described as:

∂T∂t=κ(^∂2T∂^x2+^∂2T∂^y2)

wheretis the time;κis the thermal diffusivity defined asκ=kρ_Cpconsisting of thermal conductivity (k), density (ρ) and heat capacity (_CP ) (Table 1). [41] improved the geotherm evolution, which is expressed by the horizontal and vertical velocities:

u(x)=lnβ∆tx (0<t<∆t), u(x)=0 (t≥∆t)

v(y)=lnβ∆ty (0<t<∆t), v(y)=0 (t≥∆t)

whereu(x)is the horizontal velocity;v(y) is the vertical velocity (Figure 2a). The steady-state geotherm before rifting (∆t= 0) is given by:

_Tsteady= _T1×y_yL

The temperature during the extension is calculated using the energy equation, given by:

∂T∂t+u∂T∂x+v∂T∂y=κ(^∂2T∂^x2+^∂2T∂^y2)

After the end of extension (post-rift phase), the energy equation is transformed to the heat equation because all the velocities are assumed as zero. Since the thinned lithosphere cools down due to heat conduction, it recovers its initial thickness and steady-state conditions (i.e., thermal re-equilibration).

We evaluate subsidence as the vertical variation of the initial lithosphere top. Tectonic subsidence by lithospheric stretching is the major mechanism of the syn-rift phase, while thermal subsidence by thermal contraction drives the post-rift phase [10,25]. Sedimentary load is added to the tectonic and thermal subsidence models, which is regarded as total subsidence (so-called basement subsidence). The basal pressures of five units are used to calculate the amount of subsidence: air (_Pa), sediment (_Ps), upper lithosphere (_Pl1), lower lithosphere (_Pl2), and asthenospheric mantle (_Pm). The basal pressure is calculated by:

P=ρ(1−α_Tm)gh

whereρ,_Tm, andhare the density at 0 °C, mean temperature, and thickness of the unit, respectively; gis the gravity;αis the thermal expansion coefficient. Subsidence (_ZS) during the syn-rift phase is calculated by:

_ZS=_Pb−_PS0(_ρm−_ρa)g

where_Pbis the total basal pressure;_PS0is the total basal pressure at the start of syn-rift (rifting commencement);_ρmis the mean density of the asthenospheric mantle;_ρa is the density of air (Table 1). Subsidence (_ZP) during the post-rift phase is calculated by:

_ZP=_Pb−_PP0(_ρm−_ρa)g

where_PP0is the total basal pressure at the start of post-rift (rifting cessation).

3. Numerical Modelling of Heat Flow and Subsidence 3.1. Benchmarking Models without Sedimentation

Using the numerical domain in Figure 2a, we reconstructed the surface heat flow and thermal subsidence models presented in [40,41]. In our experiments, the calculated thermal structure corresponded to the outcomes of previous models within a relative error of 1%. The results constitute the background set-up of our numerical modelling. This benchmarking helped us to examine the numerical domain and to refine the best setting for stacking sediment layers. The steady-state heat flow was set to 0.8 Heat Flow Units (HFU). Figure 3a presents surface heat flow models, applyingβvalues of 1, 1.25, 1.5, 2, 4, and 8 when simultaneous rifting occurs (∆t= 0). The heat flow at 0 myr was proportional toβ; 0.8 to 6.4 HFU corresponded toβ values from 1 to 8. This is associated with higher geothermal gradients (Figure 3b) becauseβ is inversely proportional to the thickness of the lithosphere. In the post-rift phase during 100 myr, all heat flows converged to near 0.8 HFU due to gradual decreases in the geothermal gradient, leading to the contraction that constitutes thermal subsidence (Figure 3e). The models ofβ= 2 were selected to test variations in the thermal structure by changing the rifting durations,∆t , to 4.3, 8.5, 21.3, 42.7, and 85.4 myr [42]. The applied durations corresponded, respectively, to 100, 50, 20, 10, and 5 of the stretching rate (G’), as defined in [41]. During the syn-rift phase, heat flow increased due to thinning of the lithosphere and rising of the asthenospheric mantle; however, models with shorter∆t showed steeper heat increases, greater values at 0 myr, and higher geothermal gradients (Figure 3c,d). This is attributed to heat loss through the surface (top of lithosphere) during lithospheric thinning. The thermal subsidence rate is faster with shorter∆t due to greater thermal contraction during the post-rift phase (Figure 3f).

3.2. Heat Flow and Subsidence Models with Sedimentation

The time frame of our numerical modelling was set as −6 to 0 myr for the syn-rift phase, 0 myr at rifting cessation, and 0 to 6 myr for the post-rift phase. A steady-state heat flow of 0.8 HFU was established at the start of rifting; and 0 °C as the temperature at the lithosphere top (full parameters in Table 1). Heat flow and subsidence were modeled usingβvalues of 1.5, 2, 4, and 8. The sedimentation rates of 0 and 250 m/myr were applied to eachβ model of heat flow, thermal gradient, and subsidence (Figure 4 and Figure 5). According to [36], when the lithosphere remains unstretched, a necessary condition for a strong blanketing effect is a sedimentation rate of 250 m/myr. Heat flow results from competing effects between the temperature increase by lithospheric extension and temperature reduction due to thermal blanketing by cold sediments. To observe temperature deviations caused by sediments, the geothermal gradients of the upper 2 km were presented in this study. The models ofβ = 2 were selected to investigate the effect of sedimentation rate on the thermal structure and subsidence, since a factor ~2 is frequently observed in well-studied basins (e.g., [6,30,32,44]). Applied sedimentation rates were 42, 83, and 250 m/myr, which correspond to 500, 1000, and 3000 m in total thickness of stacked sediment layers, deposited over 12 myr.

During the syn-rift phase, the surface heat flow increased starting from 0.8 HFU at −6 myr. The higherβ causes greater heat flow due to greater lithospheric thinning and rising of the asthenosphere top (Figure 4a). Due to heat loss during the syn-rift, however, the values in the models for sedimentation rate, 0 m/myr, were lower than those for the simultaneous rifting model, as shown in Figure 3a. During the post-rift phase, the heat flow decreased constantly over time, ultimately reaching a stationary state. When we stacked sediment layers in the models using a uniform sedimentation rate of 250 m/myr, the surface heat flow values decreased due to sedimentation (Figure 4a). This indicates that sediment layers hinder heat transfer from the previous lithosphere top (basin basement) to the sediment top (basin surface), and this effect increases withβ; at 0 myr, 0.14 HFU was lower in the model ofβ= 1.5, while 0.5 HFU was lower in the model ofβ = 8. This decreasing trend of surface heat flow was consistent with the geothermal gradient in the upper 2 km of the sediment column (Figure 4b). The surface heat flow models ofβ = 2, applying additional sedimentation rates of 42 and 83 m/myr (Figure 4c), showed that a higher sedimentation rate with thicker sediments enhances the blanketing effect. This relation is also shown in the geothermal gradient model (Figure 4d).

In the subsidence models (Figure 5), the subsided depth was relative to the initial depth of the lithosphere top, of which the major mechanisms were tectonic subsidence and thermal subsidence, respectively, for the syn-rift and post-rift phases. Increase inβenhanced the syn-rift tectonic subsidence; however, this was significant in the post-rift thermal subsidence due to the greater magnitude of lithospheric thinning and thermal contraction; at 6 myr, the subsidence depths were 110.3 m (β= 1.5) and 391.2 m (β = 8) (Figure 5a). When we added the sedimentary load with a sedimentation rate of 250 m/myr, the subsidence depths increased to 1842.2 m; at 6 myr, values were 1952.5 m (β= 1.5) and 2233.4 m (β= 8). This indicates that a constant increase in sedimentary load has no relation toβin this model, which resulted in a significant increase in subsidence rate from 9 to 163 m/myr in the model ofβ= 1.5. In the models ofβ = 2 with different sedimentation rates (Figure 5b), we can see that increasing the sedimentary load influenced variations in subsidence depth and its rate.

Since accumulated sediment loses porosity due to mechanical compaction and diagenesis, reductions in sediment thickness occur commonly as the burial depth and vertical stress increase [10,25,27,45]. However, our models applied minimal compaction conditions to show results for the maximal possible sediment thickness. If we correct the cumulative thicknesses of 500, 1000, and 3000 m for the compaction effect, applying the mean compaction trend of sandstone [46], the resulting compacted thicknesses are ~473, ~902, and ~2399 m, respectively. According to [38] (and references therein), ~25% mean compaction is estimated as sediments are altered during burial. In the case of only soft sediments, the compaction effect is higher; rates of porosity loss and thickness reduction increase to ~50%, while compaction is much lower in hard sediments. Thus, we estimate subsidence curves applying a constant thickness reduction with a maximum compaction of ~50% (dotted lines in Figure 5), which decreases the subsidence rate.

4. Geological Application, Limitations, and Requirements

Our numerical results introduce variations to the surface heat flow, geothermal gradient, and subsidence with increasing sediment thickness in the continental rifting system. Multi-parameter estimates help us to gain insight into effects of sedimentation on the thermal and subsidence evolution from syn-rift to early post-rift phases, which can be applied to improve modelling techniques and support more comprehensive and sophisticated basin studies as well as exploration projects (e.g., [47,48,49]). However, our numerical modelling did not consider some processes that occur during geodynamical evolution and sediment accumulation, which can introduce uncertainties when modelling various geological settings. Factors such as polyphase rifting, variable sedimentation rates, and sedimentary compaction need to be considered in numerical as well as field-based studies to investigate continental rifting and extensional basin evolution.

Continental rifting systems and their associated basins develop commonly in multiple syn-rift and post-rift phases (e.g., [1,5,8,15,24,26,31,50]). Since the extension time is a matter of heat flow and subsequent thermal subsidence, stepwise/variable rifting processes influence the lithospheric and basin thicknesses. When rifting does not lead to continental breakup (seafloor spreading), failed (fossil) rifts may present an alternate/anomalous syn- to post-rift evolution (e.g., [12,51]). Thus, controlling the extension rate over geologic time is an important factor during modelling. For example, the initial rift system between the Australian, Antarctic, and Indian plates (East Gondwana) was caused by a Permian extensional episode; however, it was followed by a period of tectonic quiescence in the mid-Jurassic period and minor extensional events during the Triassic and Early Jurassic. Renewed rifting during the mid-Jurassic contributed to the breakup between Australia–Antarctica and India during the Early Cretaceous [3,4,8,9,50,52]. However, continental breakup was not immediately followed by thermal subsidence along the southwest Australian rifted margin [28]. Similar cases have been reported by studies in northwestern Zealandia [53] and the Orphan Basin, Canada [54]. The lack of substantial post-rift subsidence could be induced by the effects of mantle plumes such as thermal buoyancy, mantle instability, and small-scale convection. We infer that the postponed/extended post-rift phase influences the rates of thermal contraction and subsidence. In the mid-Cretaceous period (~100 Ma), the onset of rifting between Australian and Antarctic plates was initiated by a plate motion change [3,4,8]. From the late Cretaceous, rifting transitioned diachronously to seafloor spreading from west to east along the southern Australian margin [52,55]. It was followed by normal seafloor spreading during the Eocene [56] and northward acceleration of the Australian plate, which is still ongoing [57]. The complex breakup history between Australia and Antarctica played a significant role in the evolution of the Southern Rift System and associated sedimentary basins [58,59,60]. Thus, understanding the variable rates of rifting and seafloor spreading is required to reconstruct the thermal structure, basin subsidence, as well as paleogeography, around the southern Australian rifted margins from the Cretaceous to Cenozoic.

As sediment transport commonly occurs downward from source to accommodation areas, topographic uplift and basement subsidence are crucial factors that determine sediment supply, sedimentation rate, and distribution [2,7,26,61]. The sedimentation rate usually keeps pace with accommodation and is variable by provenance, depositional system, and eustatic sea-level change (e.g., [39,62,63]). In our modelling, however, a constant sedimentation rate was applied, regardless of the topographic uplift and accommodation volume induced by tectonic and thermal subsidence as well as environmental effects. This setting is advantageous to observe self-standing variations to heat flow and subsidence resulting from a constant increase in sediment thickness during the syn-rift to early post-rift phase. However, changes in sedimentation rate and bathymetry would be crucial to capture when reconstructing the basin evolution. Since the depositional environment of the syn-rift to early post-rift phase is commonly a terrestrial to shallow marine setting (e.g., [23,28,32,64]), we assume that bathymetric increases following subsidence have little impact on this model setting. When modeling deep marine environments during rifting phases (e.g., [59]), however, significant bathymetric increases should be considered.

Another aspect to consider is the compaction of accumulated sediments. Compaction results fundamentally from the decrease in pore size (porosity) among grains, which is caused by sedimentary load due to increasing depth and vertical stress (mechanical compaction). In general, the porosity decreases to near zero at depths greater than 4 km [10,45,46], which corresponds to the thickness reduction of sedimentary layers. The porosity and thickness at the time and depth of interest are commonly calculated by applying an appropriate exponential relation between porosity and depth [6,25,32,39]. However, compaction is highly variable depending on the primary lithology, diagenesis (physicochemical compaction and cementation), as well as the quantification method [27,38,45,65]. Since the thermal diffusivity of sediments is dependent on porosity, the bulk thermal conductivity increases with depth [36]. Changes in sediment thickness can significantly influence the blanketing effect as well as subsidence rate [33]. Thus, our models show results with the maximum possible sediment thickness and minimal compaction condition, which effectively estimates an upper bound for the sedimentation effect. If the numerical model requires sophisticated controls for sediment thickness over time, the primary lithology of sediments needs to be determined, which affects the initial porosity and compaction trends [46]. Since syn-rift structures are commonly filled with terrestrial to shallow-water siliciclastics, which change to clays and carbonates as the waterdepth increased during the post-rift phase (e.g., [28,32,50,64]), a progressive transition of primary lithology, as well as associated compaction rate through time and depth, should be applied to the modelling.

5. Conclusions

In order to investigate the effects of sedimentation on the evolution of continental rifting, such as the blanketing effect and sedimentary load, we developed a numerical model that simulates heat flow and subsidence from syn-rift to early post-rift phases over 12 myr. The geological setting for continental rifting was reconstructed using simulation software, COMSOL Multiphysics^®. The surface heat flow and thermal subsidence models presented in previous studies were reconstructed, results of which constitute the background set-up for our numerical modelling. We applied a constant stacking of 5 m thick sediment layers to the numerical setting with variations ofβfrom 1 to 8 as well as sedimentation rates from 0 to 250 m/myr. Our results show that the surface heat flow values decreased as sedimentation thickness increased. This indicates that the sediment layer hinders the heat transfer from the previous lithosphere top (basin basement) to the sediment top (basin surface). The decrease in surface heat flow is consistent with changes in the geothermal gradient, which show that higherβand sedimentation rates enhance the blanketing effect. We present the subsided basin basement depth relative to the initial depth of the lithosphere top, which gives fundamental insights into tectonic and thermal subsidence, respectively, for the syn-rift and post-rift phases. When we applied sedimentation up to 3000 m thick over 12 myr, the sedimentary load had a critical effect on increasing the subsidence rate. Sedimentary compaction, thereby reducing sediment thickness, can change the subsidence rate. Although the numerical model does not consider processes such as progressive variations in rifting rate, sedimentation rate, and compaction effect, which occur during geodynamic evolution and sediment accumulation, our results present intuitive numerical models to understand the effects ofβand sedimentation rate on the thermal structure and subsidence rate.

Figure 2.(a) Numerical lithospheric extension model (modified from [41]). A schematic diagram showing boundary conditions and velocity streamlines. Model domain and geothermal gradient ofβ= 2 and extension duration∆t= 42.7 myr. Steady-state temperature profile (∆t= 0 myr) is shown (dashed line). (b) Numerical lithospheric extension model with successive stacking of 5 m thick sediment layers. A schematic diagram showing boundary condition, velocity streamlines, predefined accommodation, and sediment layers (dotted brown line). Model domain and geothermal gradient are shown with the steady-state temperature profile.L: lithosphere,A: asthenospheric mantle,S: predefined accommodation,T0: temperature of lithosphere top,T1: temperature of asthenospheric mantle,u(x): horizontal velocity,v(y): vertical velocity.

Figure 3. Numerical model settings in this study; (a) surface heat flow and (b) geothermal gradient models, applyingβ= 1, 1.25, 1.5, 2, 4, and 8. (c) Surface heat flow and (d) geothermal gradient models ofβ= 2, applying∆t= 0, 4.3, 8.5, 21.3, 42.7, and 85.4. (e,f) Post-rift thermal subsidence models, applying the above values ofβand∆t . Analytic and numerical solutions are from [40,41,42]. Negative and positive values on the time axis indicate the syn-rift and post-rift phases, respectively. Geothermal gradient models depict rifting cessation at time 0 myr; 1 Heat Flow Unit (HFU) corresponds to 41.84 mW/m2. The models whereβ= 1 (dashed lines) represents the steady-state with 0.8 HFU.

Figure 4. Surface heat flow and geothermal variation by sedimentation. Negative and positive values on the time axis indicate the syn-rift and post-rift phases, respectively; each phase duration is 6 myr in the model. The geothermal gradient model depicts the upper 2 km at time 0 myr (rifting cessation). Models of (a) surface heat flow and (b) geothermal gradient, applyingβ= 1.5, 2, 4 and 8. The dashed and solid lines are calculated with sedimentation rates of 0 and 250 m/myr, respectively. Models of (c) surface heat flow and (d) geothermal gradient ofβ= 2, applying sedimentation rates of 0, 42, 83 and 250 m/myr (total stacked sediment layers of 0, 500, 1000 and 3000 m in thickness, respectively).

Figure 5. Subsidence depth and rate variation by sedimentation. Negative and positive values on the time axis indicate the syn-rift and post-rift phases, respectively; each phase duration is 6 myr in the model. The subsidence depth is relative to the initial depth of the lithosphere top. (a) Subsidence models, applyingβ= 1.5, 2, 4, and 8. The dashed and solid lines are calculated with sedimentation rates of 0 and 250 m/myr, respectively. (b) Subsidence models ofβ= 2, applying sedimentation rates of 0, 42, 83, and 250 m/myr (total stacked sediment layers of 0, 500, 1000, and 3000 m thickness, respectively). Averaged subsidence rates during 12 myr are shown (grey numbers). Dotted lines present the subsidence curves estimated with a maximum sedimentary compaction of ~50%.

Symbol	Description	Unit	Value
yL	Thickness of lithosphere	km	125
T0	Temperature of the surface	°C	0
T1	Temperature of the asthenospheric mantle	°C	1333
ρl	Density of the lithosphere at 0 °C	kg/m3	2800
ρl1	Density of the upper lithosphere (crust) at 0 °C	kg/m3	2800
ρl2	Density of the lower lithosphere at 0 °C	kg/m3	3330
ρm	Density of the asthenospheric mantle at 0 °C	kg/m3	3330
ρs	Density of the sediments at 0 °C	kg/m3	2000
ρa	Density of the air at 15 °C	kg/m3	1.225
k	Thermal conductivity	W/(m·°C)	3.139
Cp	Heat capacity	J/kg·°C	1.402 × 103 [40] 1.172 × 103 [41]
κ	Thermal diffusivity (κ=kρl Cp )	m2/s	7.994 × 10−7 [40] 8.045 × 10−7 [41]
g	Gravity	m/s2	9.81
α	Thermal expansion coefficient	°C−1	3.28 × 10−5

Author Contributions

Conceptualization, Y.K.; Methodology, Y.K.; Software, Y.K.; Validation, Y.K. and E.Y.L.; Formal Analysis, Y.K. and E.Y.L.; Investigation, Y.K. and E.Y.L.; Resources, Y.K.; Data Curation, Y.K.; Writing-Original Draft Preparation, Y.K.; Writing-Review and Editing, Y.K., M.H., E.Y.L.; Visualization, Y.K. and E.Y.L.; Supervision, M.H.; Project Administration, M.H.; Funding Acquisition, M.H. and E.Y.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by KRF program by the Ministry of Science and ICT through the NRF of Korea (2017H1D3A1A01054745) and K-IODP by the Ministry of Oceans and Fisheries, Korea.

Acknowledgments

We acknowledge discussions with C. Lee on numerical modelling techniques. The authors would like to thank the Editor and two reviewers for their valuable comments, which significantly improved the manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the choice of research project; in the design of the study; in the collection, analyses or interpretation of data and results; in the writing of the manuscript; or in the decision to publish the results.

References

1. Corti, G. Continental rift evolution: From rift initiation to incipient break-up in the Main Ethiopian Rift, East Africa. Earth Sci. Rev. 2009, 96, 1-53.

2. Cowie, P.A.; Underhill, J.R.; Behn, M.D.; Lin, J.; Gill, C.E. Spatio-temporal evolution of strain accumulation derived from multi-scale observations of Late Jurassic rifting in the northern North Sea: A critical test of models for lithospheric extension. Earth Planet. Sci. Lett. 2005, 234, 401-419.

3. Gaina, C.; Müller, R.D.; Brown, B.; Ishihara, T.; Ivanov, S. Breakup and early seafloor spreading between India and Antarctica. Geophys. J. Int. 2007, 170, 151-169.

4. Gibbons, A.D.; Barckhausen, U.; van den Bogaard, P.; Hoernle, K.; Werner, R.; Whittaker, J.M.; Müller, R.D. Constraining the Jurassic extent of Greater India: Tectonic evolution of the West Australian margin. Geochem. Geophys. Geosyst. 2012, 13, Q05W13.

5. Morgan, P.; Seager, W.R.; Golombek, M.P. Cenozoic thermal, mechanical and tectonic evolution of the Rio Grande rift. J. Geophys. Res. Solid Earth 1986, 91, 6263-6276.

6. Sclater, J.G.; Christie, P.A. Continental stretching: An explanation of the post-mid-Cretaceous subsidence of the central North Sea basin. J. Geophys. Res. Solid Earth 1980, 85, 3711-3739.

7. Steckler, M.S. Uplift and extension at the Gulf of Suez: Indications of induced mantle convection. Nature 1985, 317, 135-139.

8. Veevers, J.J. Updated Gondwana (Permian-Cretaceous) earth history of Australia. Gondwana Res. 2006, 9, 231-260.

9. Williams, S.E.; Whittaker, J.M.; Granot, R.; Müller, D.R. Early India-Australia spreading history revealed by newly detected Mesozoic magnetic anomalies in the Perth Abyssal Plain. J. Geophys. Res. Solid Earth 2013, 118, 3275-3284.

10. Allen, P.A.; Allen, J.R. Basin Analysis: Principles and Application to Petroleum Play Assessment, 3rd ed.; John Wiley and Sons: New York, NY, USA, 2013.

11. Peron-Pinvidic, G.; Manatschal, G. The "IMAGinING RIFTING" Workshop Participants, Rifted Margins: State of the Art and Future Challenges. Front. Earth Sci. 2019, 7, 218.

12. Stein, S.; Stein, C.A.; Elling, R.; Kley, J.; Keller, G.R.; Wysession, M.; Rooney, T.; Frederiksen, A.; Moucha, R. Insights from North America's failed Midcontinent Rift into the evolution of continental rifts and passive continental margins. Tectonophysics 2018, 744, 403-421.

13. Whitmarsh, R.B.; Manatschal, G.; Minshull, T.A. Evolution of magma-poor continental margins from rifting to seafloor spreading. Nature 2001, 413, 150-154.

14. Buck, W.R.; Martinez, F.; Steckler, M.S.; Cochran, J.R. Thermal consequences of lithospheric extension: Pure and simple. Tectonics 1988, 7, 213-234.

15. Franke, D. Rifting, lithosphere breakup and volcanism: Comparison of magma-poor and volcanic rifted margins. Mar. Pet. Geol. 2013, 43, 63-87.

16. Gueydan, F.; Morency, C.; Brun, J.-P. Continental rifting as a function of lithosphere mantle strength. Tectonophysics 2008, 460, 83-93.

17. Huismans, R.S.; Beaumont, C. Symmetric and asymmetric lithospheric extension: Relative effects of frictional-plastic and viscous strain softening. J. Geophys. Res. Solid Earth 2003, 108, 2496.

18. Koptev, A.; Burov, E.; Gerya, T.; Le Pourhiet, L.; Leroy, S.; Calais, E.; Jolivet, L. Plume-induced continental rifting and break-up in ultra-slow extension context: Insights from 3D numerical modeling. Tectonophysics 2018, 746, 121-137.

19. Liao, J.; Gerya, T. From continental rifting to seafloor spreading: Insight from 3D thermo-mechanical modeling. Gondwana Res. 2015, 28, 1329-1343.

20. Storey, B.C. The role of mantle plumes in continental breakup: Case histories from Gondwanaland. Nature 1995, 377, 301-308.

21. Zeyen, H.; Volker, F.; Wehrle, V.; Fuchs, K.; Sobolev, S.V.; Altherr, R. Styles of continental rifting: Crust-mantle detachment and mantle plumes. Tectonophysics 1997, 278, 329-352.

22. Zwaan, F.; Schreurs, G.; Adam, J. Effects of sedimentation on rift segment evolution and rift interaction in orthogonal and oblique extensional settings: Insights from analogue models analysed with 4D X-ray computed tomography and digital volume correlation techniques. Glob. Planet. Chang. 2018, 171, 110-133.

23. Gawthorpe, R.L.; Leeder, M.R. Tectono-sedimentary evolution of active extensional basins. Basin Res. 2000, 12, 195-218.

24. Busby, C.; Azor, A. Tectonics of Sedimentary Basins: Recent Advances; Blackwell Publishing Ltd.: Chicester, UK, 2012.

25. Lee, E.Y.; Novotny, J.; Wagreich, M. Subsidence Analysis and Visualization for Sedimentary Basin Analysis and Modelling; SpringerBriefs in Petroleum Geoscience and Engineering; Springer: Cham, Switzerland, 2019.

26. Gac, S.; Hansford, P.A.; Faleide, J.I. Basin modelling of the SW Barents Sea. Mar. Pet. Geol. 2018, 95, 167-187.

27. Lee, E.Y.; Novotny, J.; Wagreich, M. Compaction trend estimation and applications to sedimentary basin reconstruction (BasinVis 2.0). Appl. Comput. Geosci. 2020, 5, 100015.

28. Lee, E.Y.; Wolfgring, E.; Tejada, M.L.G.; Harry, D.L.; Wainman, C.C.; Chun, S.S.; Schnetger, B.; Brumsack, H.-J.; Maritati, A.; Martinez, M.; et al. IODP Expedition 369 Science Party, Early Cretaceous subsidence of the Naturaliste Plateau defined by a new record of volcaniclastic-rich sequence at IODP Site U1513. Gondwana Res. 2020, 82, 1-11.

29. Roberts, A.M.; Kusznir, N.J.; Corfield, R.I.; Thompson, M.; Woodfine, R. Integrated tectonic basin modelling as an aid to understanding deep-water rifted continental margin structure and location. Pet. Geosci. 2013, 19, 65-88.

30. Steckler, M.S.; Watts, A.B. Subsidence of the Atlantic-type continental margin off New York. Earth Planet. Sci. Lett. 1978, 41, 1-13.

31. Steinberg, J.; Roberts, A.M.; Kusznir, N.J.; Schafer, K.; Karcz, Z. Crustal structure and post-rift evolution of the Levant Basin. Mar. Pet. Geol. 2018, 96, 522-543.

32. Zhao, Z.; Sun, Z.; Wang, Z.; Sun, Z. The mechanics of continental extension in Qiongdongnan Basin, northern South China Sea. Mar. Geophys. Res. 2014, 36, 197-210.

33. Hutchison, I. The effects of sedimentation and compaction on oceanic heat flow. Geophys. J. Int. 1985, 82, 439-459.

34. Karner, G.D. Sediment blanketing and the flexural strength of extended continental lithosphere. Basin Res. 1991, 3, 177-185.

35. Lucazeau, F.; Le Douaran, S. The blanketing effect of sediments in basins formed by extension: A numerical model. Application to the Gulf of Lion and Viking graben. Earth Planet. Sci. Lett. 1985, 74, 92-102.

36. Wangen, M. The blanketing effect in sedimentary basins. Basin Res. 1995, 7, 283-298.

37. Zhang, Y.K. The thermal blanketing effect of sediments on the rate and amount of subsidence in sedimentary basins formed by extension. Tectonophysics 1993, 218, 297-308.

38. Martín-Martín, M.; Robles-Marín, P. Alternative methods for calculating compaction in sedimentary basins. Mar. Pet. Geol. 2020, 113, 104132.

39. Gallagher, K.; Lambeck, K. Subsidence, sedimentation and sea-level changes in the Eromanga Basin, Australia. Basin Res. 1989, 2, 115-131.

40. McKenzie, D.P. Some remarks on the development of sedimentary basins. Earth Planet. Sci. Lett. 1978, 40, 25-32.

41. Jarvis, G.T.; Mckenzie, D.P. Sedimentary basin formation with finite extension rates. Earth Planet. Sci. Lett. 1980, 48, 42-52.

42. Kim, Y.; Lee, C. A benchmark for thermal structure evolution of extensional basin. J. Geol. Soc. Korea 2016, 52, 885-900.

43. Hulbert, G.M.; Chung, J. Explicit time integration algorithms for structural dynamics with optimal numerical dissipation. Comput. Methods Appl. Mech. Eng. 1996, 137, 175-188.

44. Tunwal, M.; Mulchrone, K.F.; Meere, P.A. Constraining basin parameters using a known subsidence history. Geosciences 2020, 10, 263.

45. Giles, M.R. Diagenesis: A Quantitative Perspective. Implications for Basin Modelling and Rock Property Prediction; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1997.

46. Kim, Y.; Lee, C.; Lee, E.Y. Numerical analysis of sedimentary compaction: Implications for porosity and layer thickness variation. J. Geol. Soc. Korea 2018, 54, 631-640.

47. Botor, D.; Anczkiewicz, A.A.; Mazur, S.; Siwecki, T. Post-Variscan thermal history of the Intra-Sudetic Basin (Sudetes, Bohemian Massif) based on apatite fission track analysis. Int. J. Earth Sci. 2019, 108, 2561-2576.

48. Carr, A.D.; Uguna, C.N. Some thoughts on the influence of pressure and thermal history assumptions on petroleum systems modelling. J. Pet. Geol. 2015, 38, 459-465.

49. Ceriani, A.; Di Giulio, A.; Fantoni, R.; Scotti, P. Cooling in rifting sequences during increasing burial depth due to heat flow decrease. Terra Nova 2006, 18, 365-371.

50. Harry, D.L.; Tejada, M.L.G.; Lee, E.Y.; Wolfgring, E.; Wainman, C.C.; Brumsack, H.-J.; Schnetger, B.; Kimura, J.-I.; Riquier, L.; Borissova, I.; et al. Evolution of the southwest Australian rifted continental margin during breakup of East Gondwana: Results from IODP Expedition 369. Geochem. Geophy. Geosyst. 2020. (accepted).

51. Misra, A.A.; Maitra, A.; Sinha, N.; Dey, S.; Mahapatra, S. Syn-to post-rift fault evolution in a failed rift: A reflection seismic study in central Cambay Basin (Gujarat), India. Int. J. Earth Sci. 2019, 108, 1293-1316.

52. White, L.T.; Gibson, G.M.; Lister, G.S. A reassessment of paleogeographic reconstructions of eastern Gondwana: Bringing geology back into the equation. Gondwana Res. 2013, 24, 984-998.

53. Boston, B.; Nakamura, Y.; Gallais, F.; Hackney, R.; Fujie, G.; Kodaira, S.; Miura, S.; Kaiho, Y.; Saito, S.; Shiraishi, K.; et al. Delayed subsidence after rifting and a record of breakup for Northwestern Zealandia. J. Geophys. Res. Solid Earth 2019, 124, 3057-3072.

54. Dafoe, L.T.; Keen, C.E.; Dickie, K.; Williams, G.L. Regional stratigraphy and subsidence of Orphan Basin near the time of breakup and implications for rifting processes. Basin Res. 2017, 29, 233-254.

55. Ball, P.; Eagles, G.; Ebinger, C.; McClay, K.; Totterdell, J. The spatial and temporal evolution of strain during the separation of Australia and Antarctica. Geochem. Geophys. Geosyst. 2013, 14, 2771-2799.

56. Whittaker, J.M.; Williams, S.E.; Müller, R.D. Revised tectonic evolution of the Eastern Indian Ocean. Geochem. Geophys. Geosyst. 2013, 14, 1891-1909.

57. Jones, I.; Verdel, C.; Crossingham, T.; Vasconcelos, P. Animated reconstructions of the Late Cretaceous to Cenozoic northward migration of Australia, and implications for the generation of east Australian mafic magmatism. Geosphere 2017, 13, 460-481.

58. Direen, N.G.; Stagg, H.M.J.; Symonds, P.A.; Colwell, J.B. Dominant symmetry of a conjugate southern Australian and East Antarctic magma-poor rifted margin segment. Geochem. Geophys. Geosyst. 2011, 12, Q02006.

59. MacLeod, K.G.; White, L.T.; Wainman, C.C.; Martinez, M.; Jones, M.M.; Batenburg, S.J.; Riquier, L.; Haynes, S.J.; Watkins, D.K.; Bogus, K.A.; et al. Late Cretaceous stratigraphy and paleoceanographic evolution in the Great Australian Bight Basin based on results from IODP Site U1512. Gondwana Res. 2020, 83, 80-95.

60. Sauermilch, I.; Whittaker, J.M.; Bijl, P.K.; Totterdell, J.M.; Jokat, W. Tectonic, Oceanographic, and Climatic Controls on the Cretaceous-Cenozoic Sedimentary Record of the Australian-Antarctic Basin. J. Geophys. Res. Solid Earth 2019, 124, 7699-7724.

61. Doglioni, C.; Prosser, G. Fold uplift versus regional subsidence and sedimentation rate. Mar. Pet. Geol. 1997, 14, 179-190.

62. Knierzinger, W.; Wagreich, M.; Palzer-Khomenko, M.; Gier, S.; Meszar, M.; Lee, E.Y.; Koukal, V.; Strauss, P. Provenance and palaeogeographic evolution of Lower Miocene sediments in the eastern North Alpine Foreland Basin. Swiss J.Geosci. 2019, 112, 269-286.

63. Wan, L.; Bianchi, V.; Hurter, S.; Salles, T. Evolution of a delta-canyon-fan system on a typical passive margin using stratigraphic forward modelling. Mar. Geol. 2020, 429, 106310.

64. Wainman, C.C.; Borissova, I.; Harry, D.L.; Hobbs, R.W.; Mantle, D.J.; Maritati, A.; Lee, E.Y. The Expedition 369 Scientists, Evidence for non-marine Jurassic to earliest Cretaceous sediments in the pre-breakup section of the Mentelle Basin, southwestern Australia. Aust. J. Earth Sci. 2020, 67, 89-105.

65. Perrin, C.; Marquez, X.; Flores, J.; Berthereau, G. Porosity Depth Saturation (PDS) model: Quantification of porosity preservation with burial in carbonate oil reservoirs, and application to infer oil charging time. Mar. Pet. Geol. 2020, 120, 104515.