Basins are generally classified according to their plate tectonic position (Helwig, 1985; Allen and Allen, 1990; Poelchau et al., 1997). The South Yellow Sea Basin is a back-arc extensional sedimentary basin formed during the Meso-Cenozoic (Yi et al., 2003; Lee et al., 2006; Pang et al., 2016b). The tectono-thermal evolution modeling is one of the foundations of the modeling of extensional basins. This study used the forward modeling method and employed PetroMod basin-modeling software to combine the analysis of the basin dynamics with the modeling of the petroleum system (Baur et al., 2010). With the current geothermal data as constraining conditions, the tectono-thermal evolution history of the study area was recovered both on the lithosphere and basin scale. This method was developed from the limited extension rate model (Jarvis and McKenzie, 1980), which is based on the model of instantaneous uniform pure-shear tension proposed by McKenzie (1978). At the extension stage, not only tectonic subsidence caused by thinning of the lithosphere but also subsidence caused by a thermal relaxation effect in the lithosphere were taken into consideration; hence, the method can be well matched with the 1-D back-stripping technology. Moreover, heat loss and thermal subsidence in the post-rifting period were considered during the modeling process. The contribution of the radioactive generation of heat in the crust was also superimposed.

The process is summarized as follows:

Step 1. Based on the present-day geometry, the age assignment and defined erosions by stepwise decompaction of the sediments, the Airy theory of isostatic compensation, and the back-stripping method proposed by Steckler and Watts (1978) were employed for the back-stripping equation to calculate the burial history and theoretical trend in tectonic subsidence.

Step 2. The stretching factors are determined by applying a stretching model. The corresponding geophysical parameters were set and the stretch factors were constantly adjusted to calculate the best-fit between the theoretical and calculated tectonic subsidence curves. Accordingly, the stretch factors are used to develop the heat flow history by applying the stretching model again. Then the calibration of the stretch factors should be performed. First, on the basis of the definition of stretch factor, the initial thickness of the crust and the stretch factor were used to determine the current calculated thickness of the crust, which was confirmed by comparison with the actual thickness of the crust. Next, the measured thermometric data were used to calibrate the thermal history trend.

Step 3. The thermal history obtained by the 1-D modeling and the related parameters were used as 3-D input to simulate the paleo-geothermal field and the petroleum system.

Modeling of the tectono-thermal evolution was performed for several representative well locations to reconstruct the history of their tectono-thermal evolution. On the basis of the BHH and published data (Yi et al., 2003; Gao and Zhou, 2014; Zhang et al., 2014b), the stratigraphy, lithostratigraphy, and calibration data were input into the 1-D module. The fine strata were depicted according to the generalized geogram of a single well, in which the thickness from the well bottom to the bottom boundary of the downfaulted sedimentary basin was calculated by the interpretation of the multichannel seismic data. Lithological, geophysical, and geochemical parameters were set by the combination of data from the PetroMod standard-based rock sample library (Table 1), where the values of average thermal conductivity refer to the measured samples from the South Yellow Sea Basin and the data from adjacent regions (Wang et al., 1995; Yang et al., 2003a, b). The settings for other related constraining conditions and parameters such as the crustal thickness, mantle thickness, and corresponding density and temperature on the bottom of the lithosphere are listed in Table 2. Some data refer to the research results of the Subei Basin on land (Bao et al., 2013a, b).

The PWD was mainly deduced by the analysis of sedimentary facies and biostratigraphy (Yi et al., 2003; Cao and Ye, 2008; Cai et al., 2014; Zhang et al., 2014b). The values for the lacustrine dark mudstone group, light gray mudstone group, and sandstone group were set at 20, 10, and 3 m, respectively. The SWIT was calculated from the global mean surface temperature (Wygrala, 1989), and the paleo-latitude was set at 35°N.

In the modeling process, the trend in subsidence and the theoretical trend in tectonic subsidence were first calculated by fitting to each other to obtain the corresponding stretch factor. Taking the ZC1-2-1 well as an example, the stretch factor was determined via fitting calculations of the trend in subsidence and the theoretical trend in tectonic subsidence [β(crust)=1.21, γ(mantle)=1.49] (Fig. 3). The initial crust-mantle thickness has a significant impact on the calculation of the tectonic subsidence curve (Liu and He, 2015). After the crust-mantle stretch factor was obtained, the current calculation-based crust-mantle thickness was compared with the measured value to test the validity of the initial crust-mantle thickness and the corresponding stretch factor. The reference values of the measured crust thickness were 28.5–30.0 km in the western part of the Yellow Sea and around 29 km in the eastern part of the Yellow Sea (Xu et al., 2009). On the basis of the balanced section method, previous studies have conducted similar research on the Subei Basin in the lower Yangtze region (Bao et al., 2013a) and concluded that the stretch factor for the Subei Basin since the Cretaceous period was around 1.1–1.4, whereas the numerical modeling-based stretch factor for the North Subbasin of the South Yellow Sea was 1.15–1.5. Both factors were roughly the same, which confirmed that this method could be applied for the study of the tectono-thermal evolution of such basins.

The measured geothermal data were used for the calibration of the results of modeling on a basin scale. Considering the ZC1-2-1 and H9 wells, which provided relatively complete data, in the study area as examples, the trend in HF calculated by numerical modeling served as the initial value. The initial trend in HF that was assigned in the model was constrained by the results of the study by Qiu (2004), which were calibrated by comparing measured and calculated temperature-sensitive parameters, e.g., vitrinite reflectance. According to research into the geological conditions and HF in the lower Yangtze region and adjacent regions (Wang et al., 1995; Yang et al., 2003a, b ; Qiu, 2004), the current value of HF is about 65 mW/m², which is slightly higher than the average value of 61 mW/m² for the global land area. Therefore, the fluctuation interval of the current value of HF was set at 62–66 mW/m². The EASYRo(%) model was used for modeling calculations to obtain the corresponding Ro depth curve (Sweeney and Burnham, 1990), and then the measured Ro data were used for the calibration of the trend in HF to obtain the calibrated ideal trend line for HF (Fig. 4). The measured Ro data from the ZC1-2-1 and H9 wells were mainly derived from the lower and central parts of the Funing formation and the Sanduo and Dainan formations. The combined constraints can be used for the effective calibration of the geothermal history on a basin scale since the post-rifting period in the Eocene (Fig. 4).

The main uncertainty of the modeling results of tectono-thermal evolution lies in the input values of geophysical parameters. Constrained by previous study results (Hao et al., 2010; He, 2002; Liao et al., 2009; Xu et al., 2008, 2009), we first calibrated the stretch factors and current calculated thickness of the crust and compared them with other studies (Bao et al., 2013b; Xu et al., 2008, 2009). Then the key geophysical parameters, including the thickness of initial crust and lithospheric mantle, were determined for subsequent modeling of heat flow trend. Also the trends in HF calculated by numerical modeling were used to calculate Ro evolutionary curves to check out the validity of input values by comparing them with the measured data (Fig. 4b). So this modeling method of integrating geodynamics and geothermometers can effectively enhance the reliability of the modeling results.

The Meso-Cenozoic thermal evolution processes in the North Subbasin can be divided into three phases (Fig. 4). The initial rifting event occurred during the Middle-Late Jurassic. The thick Jurassic formations that were encountered by the RC20-2-1 well represent strong and rapid subsidence. However, current drilling wells mainly reveal strata since the Cretaceous (Fig. 2). Therefore, by modeling of the tectono-thermal evolution of a single well location it is relatively difficult to reconstruct the history of the tectono-thermal evolution before the Cretaceous. The interpretation of seismic results indicates that Jurassic strata commonly developed in the deep depression area, and their distribution range is significantly smaller than that of the widespread overlying strata (Fig. 2). The results of modeling show that the Cretaceous is a key period when the North Subbasin, in general, entered into a stage of strong rifting and underwent rapid subsidence until the Late Cretaceous, followed by a stage of moderate rifting during the Paleogene. The South Yellow Sea Basin entered into a stage of depression development in the Early Miocene, when rifting mostly ceased.

There is a close relation between syn-rifting tectonic subsidence and the history of thermal evolution (Royden and Keen, 1980; Tong et al., 2005). The results of modeling the history of thermal evolution can be checked against the measured temperature-sensitive parameters. Therefore, an important application of 1D modeling is to calibrate boundary conditions, in particular, to reconstruct the thermal history on the basis of the rapid comparison of results with measured data and then to provide constraints on 3D modeling. The final trend in HF was determined after iterative fitting-based simulations. This process was conducted for well locations (H9 and ZC1-2-1) to obtain the final optimal trend. The final trend in HF matches the measured values of vitrinite reflectance well (Fig. 4). It exhibited a slow increase from an average level of 61 mW/m² during the Middle-Late Jurassic, which represents the initial rifting phase of the North Subbasin. Then, the value of HF underwent a strong increase and reached an average value of about 80 mW/m² from circa 145 Ma to circa 74 Ma. The rapid thermal subsidence indicates the syn-rifting phase of evolution of the basin. The peak value of HF was followed by a gradual decline to an approximate value of 65 mW/m² until the end of the Oligocene, which represents the post-rifting phase. Subsequently, in the depression phase, thermal subsidence nearly ceased and the HF maintained a constant value of 65 mW/m². By contrast, the HF trend of north subbasin of the South Yellow Sea is the same as the south subbasin, however, the overall value of the north subbasin is relatively lower than the south subbasin (Yang et al., 2003a, b).

In PetroBuilder 3-D, a numerical model is established according to the stratigraphic framework, lithologies, facies, and source rock properties. On the basis of data from the BHH, a well-defined conceptual model of the framework of basin filling can be established. The conceptual model is an integrated quantified framework with all parameters or inputs defined (Welte and Yalçin, 1987; Wygrala, 1988). This step requires proper inputs of the stratigraphic framework and related chronostratigraphy, geological events of deposition and erosion, and the creation and basic attributes of faults, as well as definitions of the petroleum system (e.g., the source rock, carrier layer, seal layer, and reservoir rock) and related parameters. Several key input elements that have a significant influence on the results of modeling are discussed as follows.

Abundant data and maps covering the maturation of organic matter and generation, expulsion, migration, and accumulation of hydrocarbons can be obtained from the results of 3-D modeling. The geodynamic processes in the petroleum system can be reconstructed. The results of modeling of vitrinite reflectance show that the source rocks in and below the E_F are mature (Fig. 8), which is in accordance with the measured data and the results of 1-D modeling (Cao and Ye, 2008; Cao et al., 2009) (Fig. 4). The maximum value (Ro(%)) in the E_F can reach 0.8% at the depocenter, but most values are generally near the threshold of 0.5% (Figs 4 and 8). The Paleogene Sanduo and Dainan formations (E_s+d), however, hardly reach the threshold for hydrocarbon generation. The results of modeling also indicate that, at present, the majority of Jurassic and Cretaceous source rocks in the North Subbasin have reached or exceeded peak oil generation and most source rocks have completed the generation and expulsion of oil. The burial depth of Jurassic source rocks can reach 8 000 m at the depocenter (Segments A and C; Figs 1 and 6), and the results of modeling of vitrinite reflectance (Ro(%)) can reach 3.5% (Fig. 8). However, it should be noted that the distribution of Jurassic source rocks is still uncertain.

The results reveal that the main generation and expulsion in Jurassic source rocks took place during the syn-rifting and post-rifting phases, whereas peak generation and expulsion in Cretaceous and Paleogene source rocks took place during the post-rifting phase (Figs 3 and 9). Generation was dominated by thermal subsidence. When controlled by abundant active faults, petroleum can be expelled immediately after generation, showing a good relation between generation and expulsion on the time axis. The area yields for generation and expulsion (Mt/km²) are correspondingly controlled by depocenters during the initial rifting and syn-rifting phases. In general, the efficiency of hydrocarbon expulsion is relatively high, particularly in Jurassic source rocks, and the average value is about 90%. Statistical histograms of the generation and expulsion rates of different strata intervals are shown in Fig. 9.

Flow path simulations demonstrate that most oil and gas accumulations in the North Subbasin surround depocenters and are situated close to the modeled migration pathways (Fig. 10). Flow paths originate from deep depocenters, extend along structural ridges and faults, and then pass through or terminate in pools. The basin boundary is open in the simulation process, which leads to some losses of hydrocarbons by pathways that connect to the basin margin. The results of accumulations and possible flow paths at different times show that oil and gas generated during the Jurassic and the Early Cretaceous were almost lost from the basin margin and some deep active faults. However, there are still a considerable number of residual accumulations present in the structurally uplifted areas around the depocenter. According to the density of flow paths and accumulation sites, dominant migration pathways can be predicted.

Because of the complexity of the petroleum system, corresponding geodynamic processes, and limited amount of data available from the basin, which has been explored to a limited degree, some fundamental conditions that have a major influence on the results of modeling are still uncertain. First, the interpretation of tectonic data, such as maps of layers and fault lines, almost all of which are derived from multichannel seismic data, is unclear. There are 56 sublayers that have been subdivided in the final model of the stratigraphic framework; however, only seven surface maps and a few wells are available to provide data to constrain the subdivisions (Fig. 2). In particular, the deep Jurassic strata, which have been drilled through up to a thickness of over 2 000 m by the RC20-2-1 well, appears to be the most important petroleum-generating layer according to the results of modeling, contributing to nearly half of the gross volume of the basin. Regrettably, the bottom surface of the Jurassic (Bot_J) still lacks sufficient supporting data. Second, the development and activity periods of faults control the expulsion and migration of petroleum to a great extent and have a significant influence on the vertical distribution of petroleum accumulations. Figure 10 shows that the flow paths and accumulations have closely linked relations to the gridded fault lines. The fault spatial models and related gridding processes can be reasonably constrained for the requirements of simulations (Fig. 7). However, the properties of the faults, as discussed in Section 5.1.2, still lack credible constraints. Once the physical stratigraphic framework and assignment of facies are completed, the reservoirs where oil and gas finally accumulate are mainly dominated by the definition of the properties of faults according to the results of multiple simulations. The multiple simulations based on a single-factor analysis are an effective approach to ranking uncertainty factors. However, it is essential to merge the gridding cells to decrease the simulation time, which will correspondingly influence the accuracy of the simulation results. Third, the maturity of the E_F and E_S+D is generally near the threshold (Fig. 8). They are therefore sensitive to the history of thermal evolution and related parameters. The simulation results of thermal evolution history have been checked against calibration parameters and, if necessary, the related inputs have also been modified to obtain a better match between the simulated and measured data. We compared the modeling results of the petroleum system with the results of previous studies and observed that hydrocarbon potential and expulsion have the same order of magnitude as that obtained in the previous resource evaluation results (Xiao, 2002). However, a difference in the hydrocarbon accumulation was observed. The accumulation obviously changes each time the key parameters, particularly the faults properties, are adjusted. Although the South Yellow Sea Basin is still in the initial stage of exploration, basin modeling is an efficient approach for the evaluation and exploration of basin resources.