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

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.

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.  [10,14]). The lithosphere is composed of a brittle deformed upper layer and a ductile deformed lower layer. (b) A schematic diagram of the lithospheric extension (after [40]). The lithosphere (grey line) is extended to the thinned lithosphere (black line). : lithosphere, : asthenospheric mantle, : initial thickness of lithosphere, : stretching factor, 0 : temperature of lithosphere top, 1 : temperature of asthenospheric mantle.

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 postprocessing. Since the program supports a graphical user interface (GUI) and provides example  [10,14]). The lithosphere is composed of a brittle deformed upper layer and a ductile deformed lower layer. (b) A schematic diagram of the lithospheric extension (after [40]). The lithosphere (grey line) is extended to the thinned lithosphere (black line). L: lithosphere, A: asthenospheric mantle, y L : initial thickness of lithosphere, β: stretching factor, T 0 : temperature of lithosphere top, T 1 : temperature of asthenospheric mantle.

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, Geosciences 2020, 10, 451 3 of 14 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.
Geosciences 2020, 10, x FOR PEER REVIEW 3 of 14 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 ℃ and 1333 ℃, 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 timestepping 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 ℃ 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. 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  κ

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 (y L ) of the lithosphere are changed to y L × β for width and y L /β for thickness (Figure 1b). Lithospheric thinning raises the top of the asthenospheric mantle (A) with temperature (T 1 ) to a depth of y L /β, which changes the temperature profile and geothermal gradient of the lithosphere. The geotherm is evaluated by the following equations: where T ∆t is the geotherm at rifting duration, ∆t; y is the depth. The temperature is calculated by solving the heat equation, described as: Geosciences 2020, 10, 451 where t is the time; κ is the thermal diffusivity defined as κ = k ρC p consisting of thermal conductivity (k), density (ρ) and heat capacity (C P ) ( Table 1). [41] improved the geotherm evolution, which is expressed by the horizontal and vertical velocities: where u(x) is the horizontal velocity; v(y) is the vertical velocity ( Figure 2a). The steady-state geotherm before rifting (∆t = 0) is given by: The temperature during the extension is calculated using the energy equation, given by: 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 (P a ), sediment (P s ), upper lithosphere (P l1 ), lower lithosphere (P l2 ), and asthenospheric mantle (P m ). The basal pressure is calculated by: where ρ, T m , and h are the density at 0 • C, mean temperature, and thickness of the unit, respectively; g is the gravity; α is the thermal expansion coefficient. Subsidence (Z S ) during the syn-rift phase is calculated by: where P b is the total basal pressure; P S0 is the total basal pressure at the start of syn-rift (rifting commencement); ρ m is the mean density of the asthenospheric mantle; ρ a is the density of air (Table 1). Subsidence (Z P ) during the post-rift phase is calculated by: where P P0 is the total basal pressure at the start of post-rift (rifting cessation).

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 Geosciences 2020, 10, 451 6 of 14 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).
Geosciences 2020, 10, x FOR PEER REVIEW 6 of 14 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 (∆ = 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, ∆ , 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 ∆ 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 ∆ due to greater thermal contraction during the post-rift phase (Figure 3f).

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 (Figures 4 and 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.

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 ℃ 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 (Figures 4 and 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  (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. 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  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.

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.

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.