Constraining Basin Parameters Using a Known Subsidence History

: Temperature history is one of the most important factors driving subsidence and the overall tectono ‐ stratigraphic evolution of a sedimentary basin. The McKenzie model has been widely applied for subsidence modelling and stretching factor estimation for sedimentary basins formed in an extensional tectonic environment. Subsidence modelling requires values of physical parameters (e.g., crustal thickness, lithospheric thickness, stretching factor) that may not always be available. With a given subsidence history of a basin estimated using a stratigraphic backstripping method, these parameters can be estimated by quantitatively comparing the known subsidence curve with modelled subsidence curves. In this contribution, a method to compare known and modelled subsidence curves is presented, aiming to constrain valid combinations of the stretching factor, crustal thickness, and lithospheric thickness of a basin. Furthermore, a numerical model is presented that takes into account the effect of sedimentary cover on thermal history and subsidence modelling of a basin. The parameter fitting method presented here is first applied to synthetically generated subsidence curves. Next, a case study using a known subsidence curve from the Campos Basin, offshore Brazil, is considered. The range of stretching factors estimated for the Campos basin from this study is in accordance with previous work, with an additional estimate of corresponding lithospheric thickness. This study provides insight into the dependence of thermal history and subsidence modelling methods on assumptions regarding model input parameters. This methodology also allows for the estimation of valid combinations of physical lithospheric parameters, where the subsidence history is known.


Introduction
The thermal history of a sedimentary basin is inherently associated with the basinʹs tectonostratigraphic evolution history. Changes in temperature in a sedimentary basin are governed by its geodynamic evolution, and vice versa. The geodynamic evolution history of a sedimentary basin has significant implications for its hydrocarbon bearing potential. McKenzie [1] provided the first physical concept for the formation of sedimentary basins in an extensional setting. According to McKenzie's model, basin development undergoes two phases: (1) Uniform instantaneous stretching of the lithosphere by pure shear. The thinning of lithosphere results in passive upwelling of the asthenosphere and a rise in the thermal gradient of the stretched lithosphere. (2) Gradual thermal cooling of the lithosphere. The latter phase results in increasing the density of the lithosphere as a result of cooling. Both phases are associated with subsidence owing to isostasy, and thus result in basin formation.
Following this seminal work, a number of modifications to the McKenzie model were suggested. For the lithospheric stretching, simple shear was considered by some researchers [2,3], whereas another model used a combination of simple shear in the crust and pure shear in the lithosphere [4]. Other variations in the mode of the lithospheric stretching include non-uniform (depth-dependent) stretching [5,6] and stretching over a finite (non-instantaneous) time [7]. The following factors primarily affecting the geothermal gradient were considered in some models: the blanketing effect of sediments [8,9], the effect of radiogenic heat production [10], and the effect of magmatic activity [11]. Finally, some models incorporated flexural response of the lithosphere to sediment loading [12,13] and depth of lithospheric necking along with rift shoulder erosion [14,15].
McKenzie's model presents an analytical solution for basin subsidence [16]. Numerical modelling allows relaxation of certain assumptions that were necessary to derive an analytical solution [10,17]. For example, the effect of sedimentary cover on thermal cooling is considered to be negligent in the McKenzie model. A numerical model in 1-dimension is presented here, which considers the effect of sediment cover during thermal subsidence. Subsidence models, whether analytical or numerical, are as good as the assumptions they use. All of these models depend upon the values taken for physical parameters such as density and thickness of crust and lithosphere. The aim of this paper is to compare subsidence profiles estimated through the backstripping method [18,19] with forward subsidence models (i.e., McKenzieʹs analytical model and the numerical model introduced here) and to explore optimum estimation of values for stretching factor, crustal, and lithospheric thickness. Some earlier studies have attempted the curve fitting approach; however, their aim was limited to the evaluation of only the stretching factor [20][21][22]. The multi-parameter estimates from this study may feed into more comprehensive and sophisticated basin models (e.g., [23,24]).
In this paper, a numerical approach that incorporates the effect of sedimentary cover is first presented. Next, a simple quantitative method for comparing measured subsidence curves with numerical and/or analytical outputs is proposed over the sample space of stretching factor, crustal, and lithospheric thickness. Finally, the methodology is tested using synthetic data and is then applied to a natural subsidence dataset obtained from the Campos basin, offshore Brazil [25], to estimate a range of stretching factors and lithospheric thicknesses suitable for basin modelling.

Numerical Model
The subsidence of a sedimentary basin in an extensional setting is modelled by two separate stages: (i) instantaneous stretching of the lithosphere followed by (ii) gradual thermal cooling. It is assumed that the basin is nourished with sediments and subsidence is closely accompanied by sedimentation. Furthermore, the following is assumed: 1) Uniform instantaneous stretching for the initial phase is by pure shear; 2) A constant asthenospheric temperature is maintained at a fixed depth equal to initial position of the base of the lithosphere; 3) Radiogenic heat production is ignored. The assumptions listed above are inherited from the classical McKenzie's model [1]. They help to (a) simplify the thermo-mechanical basin model and (b) present a preliminary generalized version to which further complexities can be added, for future studies, according to any basin specific requirements.

Initial Stretching
First, the instantaneous stretching phase is considered (see Figure 1).

Figure 1.
Thermal gradient before and after uniform instantaneous stretching by a factor β. The integrals of Equation 3 are annotated in front of the respective layers (sediment layer, crust, lithospheric mantle, and asthenospheric mantle are represented by yellow, green, orange, and purple colour, respectively). The crustal thickness is taken as 32 km, the lithospheric thickness before stretching is 120 km, and β is 2. The initial subsidence from Equation 7 is estimated to be 5.078 km.
The sedimentary layer deposited during the instantaneous lithospheric extension is assumed to remain at the surface temperature during the stretching phase, whereas the temperature of the asthenosphere remains constant at . The density of crust and mantle at any temperature is given by the first order approximation: , 1 (2) Where is the thermal expansion factor and , and , are the density of crust and mantle, respectively, at 0 ℃. The density of sediments is assumed to remain constant. Owing to the rise of lithospheric temperature after stretching (see Figure 1), the density of the crust and mantle will decrease. However, the thinning of the lithosphere leads to an initial subsidence in order to maintain isostatic equilibrium. The isostasy equation for the pressure at the depth is given by the following: are thickness of crust and lithosphere, respectively, prior to stretching. The stretching factor is defined to be the pre-stretched lithosphere thickness divided by the post initial stretching lithospheric thickness. The left hand-side of Equation 3 represents the pressure at the boundary of lithosphere and asthenosphere at depth before stretching. The temperature dependence of crust and mantle is based on thermal profile , which is the thermal gradient at undisturbed state prior to stretching Equation 5. The right hand-side of Equation 3 denotes the pressure after uniform instantaneous stretching. The thermal gradient at this state is for density estimation of the crust and mantle. The integrals terms of Equation 3 are shown in Figure 1 for visualisation. Replacing the density terms in Equation 3 from Equation 1 and Equation 2 gives the following: The thermal gradients and (see Figure 1) as depth-dependent functions are given by the following: Replacing the terms and in Equation 4 by Equation 5 and Equation 6 gives the initial subsidence as follows:

Gradual Subsidence
After instantaneous stretching, there is a gradual conductive cooling of the lithosphere, which leads to an increase in crustal and mantle density (Equations 1 and 2). In order to maintain an isostatic balance, the basin gradually subsides. The model developed here maintains two boundary conditions: (1) at a depth and below, the asthenospheric mantle is at a constant temperature throughout the stretching and cooling phase; and (2) the surface temperature remains constant at . The 1-dimensional equation for conductive cooling is given by the following: where is the thermal diffusivity, a material property, which is taken to be the same for all layers. This thermal diffusion equation is solved using an implicit Euler's method [26].
Ongoing basin development owing to thermal subsidence leads to the deposition of a thickening sedimentary layer ( at any given time ) on top of the initially deposited sediments, (see Figure  2). During each small time step ∆ , the basin subsides by a small depth ∆ and is instantaneously infilled by sediment. Therefore, is the sum of each ∆ for each ∆ up to . The calculation of ∆ , at each time step, depends upon the isostatic balance of pressure at the depth . The isostasy equations at time and ∆ are given by Equations 9 and 10, respectively. The equations compare pressure at depth after initial stretching at the start of thermal cooling (on the left hand-side) and pressure at the same depth at time (Equation 9) and ∆ (Equation 10) respectively. At any time , the isostasy equation is as follows: , , The integral terms are annotated in Figure 2 to their corresponding layers present in the system. They include the term , , the thermal transient profile at any time . The isostasy equation at the time ∆ is given by the following: Replacing density terms by Equation 1 and Equation 2 in Equation 10 gives the following: The thermal profile of the isostatic system considered during the cooling stage is affected by the sediments deposited after initial stretching and during the thermal subsidence phase. This can be noted from the inter-dependency of terms and , in Equation 9 and , . , In other words, the temperature profile in the lithosphere and asthenosphere throughout the time step, that is, from to ∆ , is approximated by the profile at time . Using the approximations in Equation 12 and Equation 13, the thermal subsidence for a given time step is calculated as follows: where the terms and are constants given by the following: .
The thickness of the sedimentary layer deposited during thermal subsidence varies from 0 at the start of cooling to at → ∞. The depth in the mantle is assumed to remain at a constant temperature according to the boundary condition. Thus, at → ∞ with the thermal equilibrium of the cooling system, the final thermal profile for this model will be / . Please note that the final thermal profile is same as the thermal profile of the system before initial stretching . Another numerical model was considered (see Supplementary File: Numerical Model 2.docx; Numerical Model 2 code.ipynb) in which the uppermost point in the asthenosphere, which remains at a constant temperature , subsides with the rate of thermal subsidence. However, it was found that the subsidence results do not vary significantly from the model presented here.

Curve Fitting Method
Subsidence modelling, either analytical (McKenzie's model) or numerical, requires constant numerical values of model parameters for calculation. Some of the typical parameters used in the models are as follows: stretching factor , lithospheric thickness , crustal thickness , density of sedimentary layer , crust and mantle , constants for thermal expansion , and thermal diffusivity . The density of the different layers as well as thermal expansibility and diffusivity constant can be provided with known and estimated values of the parameters [16]. In the case of crustal thickness prior to stretching, estimates from deep seismic imaging and gravity profiling of adjoining unstretched crust can be used as proxies [27,28]. In the absence of such geophysical measurements, it is difficult to estimate either the crust or the stretching factor. The initial lithospheric depth, however, is difficult to estimate and there appears to be no immediate guidance on how to find location-specific values.
The fitting method described below uses a known subsidence history that can be estimated using the stratigraphic backstripping technique. Values of lithospheric thickness, stretching factor, or crustal thickness can be estimated using the fitting method described here. For a fixed value of crustal thickness, lithospheric thickness varies from 100 km to 200 km with a 5 km step. In the case of the stretching factor, the range is 1.25 to 4.00 with an increment of 0.05 at each step. The subsidence models are run for 1176 (21 × 56) combinations of lithospheric depth and stretching factor and for a fixed value of crustal depth.
For each parameter combination, the fit between the calculated subsidence curve and the known subsidence curve is estimated. The error between the two curves is estimated by the following: (17) (18) Where and are the known subsidence and modelled subsidence at time step , respectively. The modelled subsidence values may be from either the McKenzie model or the numerical model or some other model variation (not considered here). The variable goes from 0 to , which is the total number of time steps in the model. In the domain space for lithospheric thickness and stretching factor values, the minimum region represented by a zero-line on a contour plot corresponds to the region where the sum of the difference between subsidence and modelled subsidence curves is zero. It is not necessarily the best fit curve. Meanwhile, the minima of in the parameter domain representing the best fit curve can be quite an extensive region. Our results indicate that the minimum curve typically lies within the region of the minimum best fit curve (the minimum region) and that the zero-line of represents a good approximation for the lithospheric thickness and the corresponding stretching factor value for the best fit case. However, it must be checked that the zero-line of lies within the minima region of contour plot. The above steps can be repeated for different values of crustal thickness.
Synthetic subsidence profiles are generated and used to estimate the suitable best fit lithospheric thickness and stretching factors for a set of crustal thickness values (32,40, and 48 km). The subsidence profiles were created by randomly generating intervals of subsidence rates for the time period 130 My to present. The first 20 My is assumed to be for initial subsidence by stretching and the remaining 110 My for thermal subsidence. These values are chosen because they are of the same order of magnitude of typical extensional basin systems [29,30] and match the primary example system investigated in this paper.
The curve fitting method is applied to a synthetic subsidence profile shown in Figure 3(A). The contour-plots in Figure 3 It is to be noted that, with an increasing crustal thickness, a lower stretching factor is required for the same lithospheric thickness. On the other hand, a higher value of lithospheric thickness is required in the case of a constant stretching factor for increasing crustal thickness. Furthermore, four synthetic subsidence curves are generated (see Figure 4(A), 4(D), 4(G), and 4(J)) and their subsidence history is compared for the models. The summary of best fit zero curve of values for both models for each synthetic curve is shown in Figure 4. It is worth noting that, with the same lithospheric thickness, a higher value of stretching factor is required for fitting the subsidence curve. This is because of the effect of sediments considered during the thermal cooling stage and the different crustal and mantle density considered in the numerical model compared with the McKenzie model.

Case Study
In this section, subsidence history from the Campos Basin, offshore Brazil, is used to apply the curve fitting method described in the previous section. Thanks to its hydrocarbon prospects, the geodynamic evolution and stratigraphy are well studied [25,[31][32][33]. This basin experienced a rapid initial rifting phase [25], thus suggesting that the assumption of instantaneous stretching is valid [7].
The formation of Campos Basin is related to the Mesozoic break up of Gondwana and opening of the South Atlantic. The rifting within Gondwana initiated from the southern South Atlantic during Triassic-Early Jurassic (220-200 Ma) and propagated northwards along the North Argentinian margin during the middle Jurassic (180-160 Ma). The intra continental rifting reached the southeastern Brazilian margin during late Jurassic-Early Cretaceous (140-132 Ma). The geodynamic evolution of the sedimentary basins along the southern Brazilian Atlantic margin is described in five tectonic phases [34,35].
The first is the pre-rift extensional phase that led to the onset of the separation between South America and Africa during Late Jurassic-Berriasian. The second phase represents the rapid lithospheric stretching and passive upwelling of asthenosphere during the Berriasian to late Barremian. This phase is associated with magmatic activity resulting in tholeiitic basalts flows along the newly formed continental margin. The third phase from Late Barremian to Late Aptian characterises the syn rift sag period and marks the end of extension in the region. Transitional to marine sediments were deposited during this phase. The fourth is the post rift stage from early to mid-Albian, which marks the development of mid Atlantic ridge and spreading of the sea floor. The thermal subsidence in the basins was typical of passive continental margins. Carbonates deposited during this phase signify deepening of the basin. The final phase stretches from late Albian to Holocene, characterising an increase in the bathymetry with time and resulting in deposition of deep water carbonates.
The subsidence rate for the depo-centre of the Campos basin was recently estimated by a study using backstripping and an inverse numerical modelling method [25]. This subsidence rate is used in this study to compare with the thermo-tectonic subsidence generated by the McKenzie model and the numerical model introduced in this study (see Figure 5(A)). The initial crustal thickness used for forward basin modelling is 32 km based on the result of deep seismic survey in the onshore region [28]. From the same study, values of sediment, crustal, and mantle density were used as 2300 kg/m 3 , 2800 kg/m 3 , and 3300 kg/m 3 , respectively. The constants for thermal expansion and thermal diffusivity are used as 3.28 × 10 −5 /°C and 10 −6 m 2 /sec, respectively. The result of comparing the subsidence trends provided from backstripping and the forward subsidence modelling for varying lithospheric thicknesses and stretching factors is shown in Figure  5 (B). For 32 km of crustal thickness, the zero curve of values vary within the range of 1.25 to 1.8. The corresponding lithospheric thickness increases with an increase in stretching factor from 100 km for 1.25 β value up to 200 for 1.8 β value. For a stretching factor in the range of 1.5 to 1.8, the suitable range of lithospheric thickness lies within 180 to 200 km.
In the past, using a visual comparison between the backstripped subsidence and McKenzie subsidence curve, the authors of [32] estimated the stretching factor to be 1.7 for the Campos basin. This lies within the range given in this study. However, an important distinction between the two comparisons is that thermal subsidence, in addition to initial subsidence, is calculated starting from the rift phase by [32], whereas it is estimated from the post rift phase in this study.
The numerical model, in comparison with the McKenzie model, is a more realistic model by incorporating the effect of sedimentary cover. The resulting difference can be seen in Figure 3, Figure  4, and more prominently in Figure 5. It can be noted from Figure 5(B) that, the higher the value of lithospheric thickness used, the greater the difference between the two models. The numerical model requires a lower stretching factor for the same lithospheric thickness for subsidence. The sedimentary cover considered in the numerical model accounts for the change in geothermal transients from instantaneous stretching until the basin reaches the state of geothermal equilibrium. This results in a comparatively lower thermal subsidence in the numerical model as compared with the McKenzie model for the same cooling period. Consequently, a higher amount of thermal subsidence is achieved through the McKenzie model for the same stretching factor as compared with the numerical model. This explains how a lower lithospheric thickness is required to achieve the same subsidence level in the numerical model with respect to the McKenzie model for a constant stretching factor.
Another major outcome from this study is the determination of constraints for acceptable lithospheric thickness together with the stretching factor. It can be noted from the case of the 1D model presented here that a variation in the value of lithospheric thickness can change the model output substantially. Thus, the estimation of lithospheric thickness is important for using the value as model input for further complex and sophisticated 2D and 3D thermo-tectonic numerical models. The curve fitting method presented here can be further generalised and applied for any number of physical parameters. This is a major deviation from earlier attempts at subsidence curve fitting where only the stretching factor was estimated (for example, [20][21][22]). In a comparative study [30], it has been shown that the subsidence patterns of sedimentary basins follow a pattern for certain tectonic settings. Thus, subsidence history can be used for tectonic interpretation based on subsidence curves.

Conclusions
A numerical approach to subsidence modelling in an extensional tectonic regime has been presented. The numerical model incorporates the effect of sedimentary cover in basin formation, which is absent in the McKenzie model. This difference leads to a comparatively lower thermal subsidence in the numerical model for (a) the same time period and (b) the same stretching factor, as compared with the McKenzie model. This discrepancy can have implications on basin subsidence histories and heat flow analysis.
Further, synthetically generated subsidence curves are compared with curves generated by the mathematical models using a least squares type approach. This allows estimation of valid combinations of stretching factor, crustal, and lithospheric thickness in basin models. If one of these parameters is known, that is, crustal thickness from seismic reflection analysis, then the other two parameters can be tightly constrained. This approach is used to fit model parameters to a known subsidence curve calculated from data gathered in the Campos Basin. The results indicate possible ranges of stretching factor values between 1.25 to 1.8 and corresponding lithospheric thickness between 100 and 200 km for the basin.
The work presented here can be applied to any other extensional sedimentary basin. Further complexities can be added in the model to include some modifications; for example, radiogenic heat production and finite stretching time. Further extension of this model, based on basin-specific requirements, in 2D and 3D could include the effect of spatial variations in rheological strength and lateral heat flow.