3 D-Basin Modeling of the Changling Depression , NE China : Exploring Petroleum Evolution in Deep Tight Sandstone Reservoirs

The Changling Depression is the largest and most resource-abundant reservoir in the South Songliao Basin, NE China. The petroleum evolution rules in the Lower Cretaceous deep tight sandstone reservoir are unclear. In this study, 3D basin modeling is performed to analyze the large-scale petroleum stereoscopic migration and accumulation history. The Changling Depression has a complex fault system and multiple tectonic movements. The model is calibrated by the present formation temperatures and observed maturity (vitrinite reflectance). We consider (1) three main erosion episodes during the burial history, one during the Early Cretaceous and two during the Late Cretaceous; (2) the regional heat flow distribution throughout geological history, which was calibrated by abundant measurement data; and (3) a tight sandstone porosity model, which is calibrated by experimental petrophysical parameters. The maturity levels of the Lower Cretaceous source rocks are reconstructed and showed good gas-generation potential. The highest maturity regions are in the southwestern sag and northern sag. The peak hydrocarbon generation period contributed little to the reservoir because of a lack of seal rocks. Homogenization temperature analysis of inclusions indicated two sets of critical moments of gas accumulation. The hydrocarbon filling in the Haerjin and Shuangtuozi structures occurred between 80 Ma and 66 Ma, while the Dalaoyefu and Fulongquan structures experienced long-term hydrocarbon accumulation from 100 Ma to 67 Ma. The homogenization temperatures of the fluid inclusions may indicate a certain stage of reservoir formation and, in combination with the hydrocarbon-accumulation simulation, can distinguish leakage and recharging events.


Introduction
The Changling fault depression is the largest petroliferous block in the southern Songliao Basin (Figure 1), which covers a total area of 1.3 × 10 4 km 2 [1].Hydrocarbon accumulation has mostly occurred in Cretaceous strata at a burial depth of 2000~4000 m.Upper Cretaceous hydrocarbon deposits were discovered and explored years ago and the petroleum evolution processes have been discussed [2][3][4].However, Lower Cretaceous reservoirs did not receive sufficient attention because of their deep burial depth until unconventional reservoirs and tight gas became the research hotspots [5][6][7][8][9][10][11].The Denglouku Formation (K 1 d) is a typical tight sandstone reservoir.According to petrophysical analysis, the porosity of K 1 d is between 0.3% and 10.4% and the average porosity is 5.4%, The large-scale research of hydrocarbon generation, migration and accumulation has significance for studying the evolution of a tight gas reservoir.Former research focused on the Upper Cretaceous discovered the self-generation and self-storage oil accumulation rule [3,14].Few works have focused on the Lower Cretaceous hydrocarbon systems.With burial depth over 3000 m, only deep drill wells could discover the Lower Cretaceous layers [8].Moreover, massive compact sandstone has made defining Lower Cretaceous deep tight gas reservoirs difficult.The vitrinite reflectance (VR) of the Lower Cretaceous source rocks reaches 1.3 %Ro, which indicates high maturity gas generation phase [15].
Basin modelling is an effective method to recover the geological evolution history including burial, thermal, maturity, petroleum saturation etc. [16,17].A 1D model is based on the single well [18][19][20][21][22][23] and a 2D model is based on the section line [24][25][26].1D and 2D basin modelling are convenient and rapid approaches to study the petroleum generation processes.Dong et al. used BasinMod 1D to simulate the Upper Cretaceous source rock maturation and address the potential accumulation area in the Changling Depression using the BasinFlow subsystem of the BasinMod software [3].However, 1D and 2D modelling could not show the stereoscopic migration and accumulation history as 3D basin modelling do [27][28][29].
Specifically, a 3D basin model can (1) accurately establish the strata structure of the basin and describe the important petroleum system elements (PSE), including the source rock, migration paths, reservoirs and seals; (2) consider the distribution of erosion thickness and recover the reliable heat flow history through calibration with a series of measured data; and (3) research the petroleum system as a 3D process, meaning that both vertical and horizontal hydrocarbon migration, entrapment and leakage are visually represented in the model [30].
In this study, we constructed a 3D model to address quantitatively basin-scale understanding of the Lower Cretaceous hydrocarbon evolution.This model was a dynamic expression of the sedimentary structure and geological evolution [31].Moreover, this model could reveal the organic maturation and petroleum-generation rules throughout various geological periods [28,[32][33][34].After studying scholars' research on 3D basin modeling, we believe that the model accuracy can be improved by comprehensive calibration methods [31,35,36].We investigated the hydrocarbon charging events on the basis of a fluid inclusion study.The fluid inclusion homogenization The large-scale research of hydrocarbon generation, migration and accumulation has significance for studying the evolution of a tight gas reservoir.Former research focused on the Upper Cretaceous discovered the self-generation and self-storage oil accumulation rule [3,14].Few works have focused on the Lower Cretaceous hydrocarbon systems.With burial depth over 3000 m, only deep drill wells could discover the Lower Cretaceous layers [8].Moreover, massive compact sandstone has made defining Lower Cretaceous deep tight gas reservoirs difficult.The vitrinite reflectance (VR) of the Lower Cretaceous source rocks reaches 1.3 %Ro, which indicates high maturity gas generation phase [15].
Basin modelling is an effective method to recover the geological evolution history including burial, thermal, maturity, petroleum saturation etc. [16,17].A 1D model is based on the single well [18][19][20][21][22][23] and a 2D model is based on the section line [24][25][26].1D and 2D basin modelling are convenient and rapid approaches to study the petroleum generation processes.Dong et al. used BasinMod 1D to simulate the Upper Cretaceous source rock maturation and address the potential accumulation area in the Changling Depression using the BasinFlow subsystem of the BasinMod software [3].However, 1D and 2D modelling could not show the stereoscopic migration and accumulation history as 3D basin modelling do [27][28][29].
Specifically, a 3D basin model can (1) accurately establish the strata structure of the basin and describe the important petroleum system elements (PSE), including the source rock, migration paths, reservoirs and seals; (2) consider the distribution of erosion thickness and recover the reliable heat flow history through calibration with a series of measured data; and (3) research the petroleum system as a 3D process, meaning that both vertical and horizontal hydrocarbon migration, entrapment and leakage are visually represented in the model [30].
In this study, we constructed a 3D model to address quantitatively basin-scale understanding of the Lower Cretaceous hydrocarbon evolution.This model was a dynamic expression of the sedimentary structure and geological evolution [31].Moreover, this model could reveal the organic maturation and petroleum-generation rules throughout various geological periods [28,[32][33][34].After studying scholars' research on 3D basin modeling, we believe that the model accuracy can be improved by comprehensive calibration methods [31,35,36].We investigated the hydrocarbon charging events on Energies 2019, 12, 1043 3 of 27 the basis of a fluid inclusion study.The fluid inclusion homogenization temperature in conjunction with burial and thermal history can determine the petroleum charge history [37][38][39].By using the 3D basin modelling approach, we tried to research the unconventional hydrocarbon migration and accumulation rules in a tight sandstone reservoir.Therefore, we expected to observe the formation rules of tight gas reservoirs in the Changling Depression.

Geological Evolution
The Changling fault depression is located in the southern Songliao Basin, Northeast China.This depression lies in the central fault depression zone of Songliao Basin to the south of Nen River and Songhua River (Figure 1).Its overall distribution is in the NNE direction.This area can be divided into five sub-structural units: a western slope zone, eastern slope zone, central depression, northern fault sag, and southern fault sag [12,40].The geological evolution of the Changling fault depression has been elaborated in previous studies [41,42].Figure 2 depicts the generalized lithostratigraphy and the primary tectonic events of the Changling Depression.The structural evolution history of the Changling Depression can be divided into four stages as detailed below.

1.
Basement thermal uplift stage: the paleo-continent split up at the end of the Triassic and began to disassemble inside the continent during the Jurassic.During this geological period, the Eurasian Plate drifted southeastward while the Pacific Plate began to expand and resulted in plate subduction.This sequence triggered mantle upwelling, and the upper mantle and crust extended to form the rift.Mantle materials erupted through the fissure and formed the widespread pyroclastic reservoir [43,44].

2.
Rifting stage: the previously accumulated heat energy suddenly released and created a crack in the crust.Then, the faulted period began and the Songliao Basin gained the embryonic form of a graben or half-graben morphology.In the meantime, volcanic eruptions, magmatic intrusions and regional metamorphic events occurred through deep fractures.The Upper Jurassic Huoshiling Formation (J 3 h) and the Lower Cretaceous Shahezi Formation (K 1 sh) and Yingcheng Formation (K 1 y) exhibited sedimentary systems such as alluvial fans, flooding plains, pluvial deltas, and fan deltas (Figure 2).A series of NNE faults formed during the deposition of the J 3 h stratum and caused massive volcanic eruption with andesite and basalt extrusions.During the deposition of K 1 sh, construction of the Changling depression was controlled by boundary faults.A set of dark mudstone with coal seams formed the lacustrine source rocks of the K 1 y and K 1 sh strata.The rising rift, intense fault activity and rapid tectonic subsidence resulted in the deepening and expansion of the basin.The isolated basins began to connect during the transgressive period.The main stratum deposited was K 1 d, and the sedimentary environment was similar to that during the early rifting stage [45][46][47].

3.
Depression stage: the basin entered a depression stage when the crust achieved differential settlement.Volcanic activity dissipated and the depression achieved a uniform sedimentary setting.Widespread lakes appeared in the depression alongside a range of depositional systems.
A full combination of spatial-deposition systems developed from the land to the lacustrine area, including alluvial fans, flooding plains, shallow lacustrine deltas, and deep water turbidite.The strata included the Quantou Formation (K 1 q), Qingshankou Formation (K 2 qn), Yaojia Formation (K 2 y) and Nenjiang Formation (K 2 n).During this stage, the paleo heat flow frequently changed because of fault activity [41].

4.
Tectonic inversion stage: during the last phase of the basin evolution, the Japan Sea expanded and the basin continued to exhibit extrusion.Tectonic inversion occurred when the strata were lifted and folded.The remaining lacustrine area shrank but the overall shape of the depression remained similar.The strata included the Upper Cretaceous Sifangtai Formation (K 2 s) and Mingshui Formation (K 2 m), with the depositional systems including shallow lacustrine, alluvial fan and flooding plain facies.After the Cretaceous, the lakes almost disappeared and the strata uplifted and descended in slow motion.Alluvial and diluvial deposits have been discovered in this sequence [13,41,42].

Construction of the Basin Model
There are two sets of petroleum systems in the Changling Depression and we named them upper group and lower group (Figure 2).The upper group includes K1q, K2qn, K2y and K2n which act as both source rock and reservoir rock [3] and produces oil and condensate gas [2,4].The lower group consists of J3h, K1sh, K1y and K1d which follows the regular pattern: gas migrated vertically upward and was stored in the K1d reservoir.In this study, we focus on the lower set of the petroleum system.Until now, four Lower Cretaceous tight sandstone gas fields are discovered and we selected four representative deep wells to indicate the locations: The Changling 1 gas field (Well CS1), which is located in the central depression; the Shuangtuozi gas field (Well TS6), which lies in the southern fault sag; and the Dalaoyefu gas field (Well LS1) and Fulongquan gas field (Well F10), which are both situated in the eastern slope zone (Figure 1b).Over 200 wells have been drilled in the The paleo water depth curve was modified after Zhang and Ren [49] (PSE = petroleum system elements, SR = source rock, RR = reservoir rock, PWD = paleo water depth).

Construction of the Basin Model
There are two sets of petroleum systems in the Changling Depression and we named them upper group and lower group (Figure 2).The upper group includes K 1 q, K 2 qn, K 2 y and K 2 n which act as both source rock and reservoir rock [3] and produces oil and condensate gas [2,4].The lower group consists of J 3 h, K 1 sh, K 1 y and K 1 d which follows the regular pattern: gas migrated vertically upward and was stored in the K 1 d reservoir.In this study, we focus on the lower set of the petroleum system.Until now, four Lower Cretaceous tight sandstone gas fields are discovered and we selected four representative deep wells to indicate the locations: The Changling 1 gas field (Well CS1), which is located in the central depression; the Shuangtuozi gas field (Well TS6), which lies in the southern fault sag; and the Dalaoyefu gas field (Well LS1) and Fulongquan gas field (Well F10), which are both situated in the eastern slope zone (Figure 1b).Over 200 wells have been drilled in the study area, and we have selected approximately 50 representative wells for investigation and discussion in this paper.The model was created with the PetroMod v.2015.1 ® software (Schlumberger, Oslo, Norway).The major processes were strata deposition, compaction, erosion, source rock maturation, hydrocarbon generation, expulsion, migration and accumulation.The basic basin-modeling rules and analytical equations that were used in the simulation processes are detailed in the introduction to basin modeling [30].Among the variety of petroleum-migration algorithms, we chose the hybrid model, which can solve Darcy flow equations in low-permeability areas and apply a flow-path analysis in high-permeability areas [30,50].The more detailed and precise the parameters that we imported into the model, the closer to reality the results would become.So we decided to import as much basic data as possible within the limits of the model.This model was set to be 700 grids in width and 1000 grids in length.The actual size of the model was 135 km in width and 195 km in length.The grid cell size was thus approximately 193 × 195 m (Figure 3).
study area, and we have selected approximately 50 representative wells for investigation and discussion in this paper.
The model was created with the PetroMod v.2015.1 ® software (Schlumberger, Oslo, Norway).The major processes were strata deposition, compaction, erosion, source rock maturation, hydrocarbon generation, expulsion, migration and accumulation.The basic basin-modeling rules and analytical equations that were used in the simulation processes are detailed in the introduction to basin modeling [30].Among the variety of petroleum-migration algorithms, we chose the hybrid model, which can solve Darcy flow equations in low-permeability areas and apply a flow-path analysis in high-permeability areas [30,50].The more detailed and precise the parameters that we imported into the model, the closer to reality the results would become.So we decided to import as much basic data as possible within the limits of the model.This model was set to be 700 grids in width and 1000 grids in length.The actual size of the model was 135 km in width and 195 km in length.The grid cell size was thus approximately 193 × 195 m (Figure 3).

Input Data
The structure model was established based on upper horizon structure maps of each formation which were provided in previous research [51].After digitizing and importing the depth maps to PetroMod software, we used the layer-division results from well-logging data and drilled cores for 133 individual wells to refine the layer model.We selected several through-well surveys in the model to check the conformity of the model with the well tops.After investigating the stratigraphic chronology, we found that scholars provided different classification schemes [3,4,52].We analyzed their different opinions, combined with zircon dating results [53,54], and ultimately defined the ages of each layer by using Ma et al.'s study [52] (Figure 2).

Lithology and Petroleum System Elements
Studies have been performed on the stratigraphic characteristics of the Changling Depression deep rift strata by means of sequence stratigraphy, logging geology, sedimentology and seismic stratigraphy [42,45].The definition of the lithological properties of the Early Cretaceous reservoir was based on Jiang's research [55].The J3h can be divided into three lithological types, including a combination of a large section of andesite, carbon-mudstone strata, conglomerate and clastic rocks.The K1sh comprises thin sand with interbedded mud and coal seams.The K1y comprises three sections: acidic volcanic rocks, clastic sediments, and mafic-intermediate volcanic rocks.The K1d, which is the studied reservoir rock, contains a set of sandstone and sandy conglomerates.The Lower K1q is the seal-rock layer and consists of thick shale, and the Upper K1q comprises conglomerate and sandstone (Figure 3).The facies settings for the model are shown in Table 1.

Input Data
The structure model was established based on upper horizon structure maps of each formation which were provided in previous research [51].After digitizing and importing the depth maps to PetroMod software, we used the layer-division results from well-logging data and drilled cores for 133 individual wells to refine the layer model.We selected several through-well surveys in the model to check the conformity of the model with the well tops.After investigating the stratigraphic chronology, we found that scholars provided different classification schemes [3,4,52].We analyzed their different opinions, combined with zircon dating results [53,54], and ultimately defined the ages of each layer by using Ma et al.'s study [52] (Figure 2).

Lithology and Petroleum System Elements
Studies have been performed on the stratigraphic characteristics of the Changling Depression deep rift strata by means of sequence stratigraphy, logging geology, sedimentology and seismic stratigraphy [42,45].The definition of the lithological properties of the Early Cretaceous reservoir was based on Jiang's research [55].The J 3 h can be divided into three lithological types, including a combination of a large section of andesite, carbon-mudstone strata, conglomerate and clastic rocks.The K 1 sh comprises thin sand with interbedded mud and coal seams.The K 1 y comprises three sections: acidic volcanic rocks, clastic sediments, and mafic-intermediate volcanic rocks.The K 1 d, which is the studied reservoir rock, contains a set of sandstone and sandy conglomerates.The Lower K 1 q is the seal-rock layer and consists of thick shale, and the Upper K 1 q comprises conglomerate and sandstone (Figure 3).The facies settings for the model are shown in Table 1.To determine the source rock properties from the input model, we processed 178 Rock-Eval pyrolysis analyses from 19 wells in the Changling Depression.All these measured data were acquired from the PetroChina Jilin Oilfield Administration Bureau (PCJOAB).The lack of measured data in the J 3 h meant that the source rock evaluation was based on seismic inversion and well-logging analysis [56], alongside selected geochemistry tests [57,58].In the Rock-Eval pyrolysis, S 1 (free hydrocarbon components), S 2 (hydrocarbons generated by the pyrolysis) and S 3 (carbon dioxide content in the pyrolysis) are measured in milligrams of hydrocarbon per gram of rock [59,60].In Figure 4a, K 1 sh had a total organic carbon content (TOC) that mostly surpassed 1% and an average S 1 + S 2 of 3.3 mg/g.Abundant samples from K 1 y varied from 0.03% to 9.01%, with a mean value of 1.1%.We determined K 1 sh to be a good to excellent source rock.The S 1 + S 2 value of K 1 y varied from 0.02 mg/g to 11.5 mg/g, with a mean value of 0.97 mg/g.Defining the organic-matter abundance in the source rock of K 1 y was difficult because the sample distribution was too dispersed from poor to excellent areas.
The van Krevelen diagram is often used to determine the kerogen types and genetic potential by comparing the atomic-ratio parameters [61,62].As shown in Figure 4b, most of the points were concentrated in the lower-left area of the chart, while the VR of the samples was over 1.0 %Ro.The diagram could not distinguish among the kerogen types while the points grew close under metagenesis conditions.On the other hand, several points were dispersed around the lower maturation area, and most of the points were in the type II and III areas.In this chart, we could not accurately distinguish the source rock types from different formations.
We established a graph of the hydrogen index (HI) versus the pyrolysis peak temperature (Tmax) of different source rock layers to distinguish the source rock type of each layer [63].In Figure 4c, the scatter map shows that the source rock samples from K 1 y were dispersed and followed the type II and III source rock distribution curve.When distinguishing between different areas, source rock samples from the Fulongquan and Yaoyingtai regions followed the trend of type III source rock and exhibited high maturity.Samples from the Dalaoyefu, Haerjin and Shuangtuozi regions mainly fit the trend of type III source rock.In Figure 4d, samples from K 1 sh in the Yaoyingtai region followed the trend of type II 2 , while those from the Haerjin region fell within the type III source rock area.Samples from K 1 sh appeared to have low maturity.We optimized the decoupled source rock definition to precisely define the source rocks' properties.We created thickness maps, TOC maps, and HI maps to refine the model.We collected dark-shale thickness data from the drill cores and analyzed the well-logging curves to form the isopach maps referring to former research [15,64,65].Then, we created source rock isopach maps of K1y, K1sh and J3h (Figure 5).The TOC maps and HI maps were created by using the average measured data from wells, alongside the well-logging analysis prediction results and data from prior studies [56,65,66].The figures show that hydrocarbons were mainly distributed in the southern fault sag and northern fault sag, displaying the largest thickness in these regions.The K1sh had the most abundant source rock, with a mean TOC value of 1.38%, while the maximum TOC reached 4.35%.The dark mudstones in the K1y and K1sh, which formed in a lacustrine environment, were the main source rocks of the study area [67].We optimized the decoupled source rock definition to precisely define the source rocks' properties.We created thickness maps, TOC maps, and HI maps to refine the model.We collected dark-shale thickness data from the drill cores and analyzed the well-logging curves to form the isopach maps referring to former research [15,64,65].Then, we created source rock isopach maps of K 1 y, K 1 sh and J 3 h (Figure 5).The TOC maps and HI maps were created by using the average measured data from wells, alongside the well-logging analysis prediction results and data from prior studies [56,65,66].The figures show that hydrocarbons were mainly distributed in the southern fault sag and northern fault sag, displaying the largest thickness in these regions.The K 1 sh had the most abundant source rock, with a mean TOC value of 1.38%, while the maximum TOC reached 4.35%.The dark mudstones in the K 1 y and K 1 sh, which formed in a lacustrine environment, were the main source rocks of the study area [67].Energies 2019, 12, 1043 10 of 27

Kinetic Models
The generation rates, components and timing of hydrocarbons are controlled by the kinetic properties of the kerogen [68].In the study area, the Early Cretaceous source rocks had entered the overly mature stage, and the transformation ratio (TR) of the samples was close to 100%.Thus, we did not perform a pyrolysis experiment on the Early Cretaceous source rocks; kinetics research was reasonable only for the Late Cretaceous source rocks.We constructed two set of kinetic models that fit the Changling deep source rock layers by referring to previous studies and following the basic rules.
After analyzing the hydrocarbon-component data, we discovered that the source rock properties and the hydrocarbon-expulsion rules were a good fit for the Burnham TIII kinetic model [69].
The pyrolysis experiments on the northern Songliao Basin oil shale provided a suite of kerogen kinetics data [70], which we used as a kerogen oil-production kinetic model.After examining the gas-source correlation studies, we found that the gas was mainly converted from oil cracking and kerogen pyrolysis in the Changling Depression [15,71].Thus, we also considered secondary cracking.Previous studies contributed to the creation of a gas-generation kinetic model, as shown in Figure 6a [70,[72][73][74].The kinetic model of a coal seam in K 1 sh (Figure 6b) was established from a model that was based on samples, which included Mesozoic and younger coal that formed from plants [75], combined with reports on coal-seam pyrolysis experiments in the Songliao Basin [76].

Kinetic Models
The generation rates, components and timing of hydrocarbons are controlled by the kinetic properties of the kerogen [68].In the study area, the Early Cretaceous source rocks had entered the overly mature stage, and the transformation ratio (TR) of the samples was close to 100%.Thus, we did not perform a pyrolysis experiment on the Early Cretaceous source rocks; kinetics research was reasonable only for the Late Cretaceous source rocks.We constructed two set of kinetic models that fit the Changling deep source rock layers by referring to previous studies and following the basic rules.
After analyzing the hydrocarbon-component data, we discovered that the source rock properties and the hydrocarbon-expulsion rules were a good fit for the Burnham TIII kinetic model [69].The pyrolysis experiments on the northern Songliao Basin oil shale provided a suite of kerogen kinetics data [70], which we used as a kerogen oil-production kinetic model.After examining the gas-source correlation studies, we found that the gas was mainly converted from oil cracking and kerogen pyrolysis in the Changling Depression [15,71].Thus, we also considered secondary cracking.Previous studies contributed to the creation of a gas-generation kinetic model, as shown in Figure 6a [70, [72][73][74].The kinetic model of a coal seam in K1sh (Figure 6b) was established from a model that was based on samples, which included Mesozoic and younger coal that formed from plants [75], combined with reports on coal-seam pyrolysis experiments in the Songliao Basin [76].

Fault Properties
Multiple faults with varying and complex scales and characteristics are distributed in the Cretaceous strata in the Changling Depression.The major faults mostly formed during the chasmic period of magmatic activity, so tensional faults are clearly dominant in the Changling Depression [77].As shown in Figure 7a, 23 representative and large faults were selected as study subjects to examine the processing capacity of the model [78].The fault intersections on each horizon were modified to create the model fault (Figure 7b).The fault strikes were mostly in the NE and NW

Fault Properties
Multiple faults with varying and complex scales and characteristics are distributed in the Cretaceous strata in the Changling Depression.The major faults mostly formed during the chasmic period of magmatic activity, so tensional faults are clearly dominant in the Changling Depression [77].As shown in Figure 7a, 23 representative and large faults were selected as study subjects to examine the processing capacity of the model [78].The fault intersections on each horizon were modified to create the model fault (Figure 7b).The fault strikes were mostly in the NE and NW directions, and the faults formed during the deposition of the J 3 h until the deposition of K 1 d.We imported the interpreted fault lines (intersections on each surface) into the model and created model faults.The border of the Changling Sag mostly consisted of boundary fractures, and we assumed that the basin sides were closed.

Erosion Events
Erosion events and the denudation thickness are important parameters in the recovery of burial and thermal histories.After reviewing the literature, we used the acoustic time difference and vitrinite-reflectance method to relate the main uplift and erosion of the basin to three main regional events: 1. Tectonic inversion occurred at the end of K1y deposition (128.8~124Ma).The basin was regionally uplifted and experienced differing degrees of erosion.The average erosion thickness was 417.5 m [46].2. During the end of K2n deposition (74.6~73Ma), the strata uplifted because of intense extrusion stress and underwent continuous denudation to a thickness of approximately 205 m. 3.During the late Mingshui period (66.2~23.03Ma), the Changling Depression underwent basin shrinkage, tectonic inversion and regional uplift [66].The average denudation thickness was 342.5 m.The main uplift and erosion center is in the Fulongquan region around Well F10.
Table 2 lists the denudation thicknesses of 8 wells, which were calculated from the sedimentation rate [41] to control the horizontal distribution of erosion.Apatite fission track analysis (AFTA) indicated four episodes in the denudation history [79].Furthermore, a Monte Carlo random simulation of the AFTA data revealed four turning points of the thermal evolution at 65 Ma, 43.5 Ma, 28 Ma and 15 Ma.This research suggested that the denudation period occurred at the end of the deposition of the Mingshui Formation [79].

Erosion Events
Erosion events and the denudation thickness are important parameters in the recovery of burial and thermal histories.After reviewing the literature, we used the acoustic time difference and vitrinite-reflectance method to relate the main uplift and erosion of the basin to three main regional events: 1.
Tectonic inversion occurred at the end of K 1 y deposition (128.8~124Ma).The basin was regionally uplifted and experienced differing degrees of erosion.The average erosion thickness was 417.5 m [46].

2.
During the end of K 2 n deposition (74.6~73Ma), the strata uplifted because of intense extrusion stress and underwent continuous denudation to a thickness of approximately 205 m.

3.
During the late Mingshui period (66.2~23.03Ma), the Changling Depression underwent basin shrinkage, tectonic inversion and regional uplift [66].The average denudation thickness was 342.5 m.The main uplift and erosion center is in the Fulongquan region around Well F10.
Table 2 lists the denudation thicknesses of 8 wells, which were calculated from the sedimentation rate [41] to control the horizontal distribution of erosion.Apatite fission track analysis (AFTA) indicated four episodes in the denudation history [79].Furthermore, a Monte Carlo random simulation of the AFTA data revealed four turning points of the thermal evolution at 65 Ma, 43.5 Ma, 28 Ma and 15 Ma.This research suggested that the denudation period occurred at the end of the deposition of the Mingshui Formation [79].The paleo water depth (PWD) curve (Figure 2) was modified based on previous studies [49,80].The water-level changes were revealed by obtaining the distribution characteristics of sediments, sedimentary structures, paleontological sets and ecology, authigenic minerals and wave-base surfaces [49].
The sediment-water interface temperature (SWIT) is the upper-boundary condition of heat transfer in a sedimentary basin [81].The PetroMod software contains a module to calculate the SWIT throughout geological time.It calculates the variations in average surface paleo temperatures from the sediment-water interface depth versus local latitude and considers PWD changes [82].
The heat flow (HF) history of the Songliao Basin revealed the paleo-basement and surface HF values from 19 different ages [83].The HF reached a peak value of 99 mW/m 2 from 120 Ma to 100 Ma.The HF evolution curve represents the total HF trend of the Songliao Basin (Figure 8).According to paleo-temperature studies, the geothermal gradient during the Early Cretaceous was 50~70 • C/km [84].The geothermal gradient was 42.6~48 • C/km during the Late Cretaceous, with a HF of 95~107 mW/m 2 .From the Paleogene to the present, the geothermal gradient stabilized at 37 The paleo water depth (PWD) curve (Figure 2) was modified based on previous studies [49,80].The water-level changes were revealed by obtaining the distribution characteristics of sediments, sedimentary structures, paleontological sets and ecology, authigenic minerals and wave-base surfaces [49].
The sediment-water interface temperature (SWIT) is the upper-boundary condition of heat transfer in a sedimentary basin [81].The PetroMod software contains a module to calculate the SWIT throughout geological time.It calculates the variations in average surface paleo temperatures from the sediment-water interface depth versus local latitude and considers PWD changes [82].
The heat flow (HF) history of the Songliao Basin revealed the paleo-basement and surface HF values from 19 different ages [83].The HF reached a peak value of 99 mW/m 2 from 120 Ma to 100 Ma.The HF evolution curve represents the total HF trend of the Songliao Basin (Figure 8).According to paleo-temperature studies, the geothermal gradient during the Early Cretaceous was 50~70 °C/km [84].The geothermal gradient was 42.6~48 °C/km during the Late Cretaceous, with a HF of 95~107 mW/m 2 .From the Paleogene to the present, the geothermal gradient stabilized at 37 °C/km with a HF of 69 mW/m 2 [85].

Petrophysical Settings
After studying previous research on basin modeling, we found that the reservoir was defined as the default sandstone when following the common workflows [24,28,55].We attempted to define the properties of the sandstone within the scope of the Changling Depression and distinguish among the different formations.If we study the conventional reservoirs, the lithological properties would make a slight difference to the accumulation results.However, we had to consider the porosity and permeability parameters when modeling tight gas reservoirs.If the porosity was below the cutoff value, the sandstone would lose its ability to store hydrocarbons.The K1d tight sandstone reservoir provided the accumulation conditions for hydrocarbons, but the physical properties of the reservoir should be defined manually.Additionally, the porosity of the deeper layers may influence the migration path significantly.
Therefore, Athy's Law (Equation ( 1)) was used to calibrate the porosity model of K1d, K1y and K1sh.Athy's Law manifests as a porosity versus depth curve.This law is a practical trend that assumes hydrostatic pressures and deposition with the same lithology [86]:

Petrophysical Settings
After studying previous research on basin modeling, we found that the reservoir was defined as the default sandstone when following the common workflows [24,28,55].We attempted to define the properties of the sandstone within the scope of the Changling Depression and distinguish among the different formations.If we study the conventional reservoirs, the lithological properties would make a slight difference to the accumulation results.However, we had to consider the porosity and permeability parameters when modeling tight gas reservoirs.If the porosity was below the cutoff value, the sandstone would lose its ability to store hydrocarbons.The K 1 d tight sandstone reservoir provided the accumulation conditions for hydrocarbons, but the physical properties of the reservoir should be defined manually.Additionally, the porosity of the deeper layers may influence the migration path significantly.
Therefore, Athy's Law (Equation ( 1)) was used to calibrate the porosity model of K 1 d, K 1 y and K 1 sh.Athy's Law manifests as a porosity versus depth curve.This law is a practical trend that assumes hydrostatic pressures and deposition with the same lithology [86]: Energies 2019, 12, 1043 13 of 27 where ϕ(z) refers to the porosity at depth z; ϕ 1 is the minimum porosity in the deposition; ϕ 0 is the initial porosity of the rock before deposition; and k is Athy's factor, a constant that depends on the rock type and physical properties.
The initial porosity was estimated from pore-recovery research.The sorting features and grain sizes of rocks indicate the general interval of the original porosity.We used the average rock-density test values and the minimum porosity data in Table 3 to shape the initial curve.Next, we calibrated Athy's porosity versus depth curve by using 851 points of measured porosity data from 40 wells in the Changling Depression (Figure 9).where   refers to the porosity at depth z;  is the minimum porosity in the deposition;  is the initial porosity of the rock before deposition; and k is Athy's factor, a constant that depends on the rock type and physical properties.The initial porosity was estimated from pore-recovery research.The sorting features and grain sizes of rocks indicate the general interval of the original porosity.We used the average rock-density test values and the minimum porosity data in Table 3 to shape the initial curve.Next, we calibrated Athy's porosity versus depth curve by using 851 points of measured porosity data from 40 wells in the Changling Depression (Figure 9).We also established a correlation between the porosity and permeability, Equation ( 2) is the formula for K1d (Figure 10).
where K is the permeability, log(mD);  is the porosity, %.  = 0.77.We also established a correlation between the porosity and permeability, Equation ( 2) is the formula for K 1 d (Figure 10).
where K is the permeability, log(mD); ϕ is the porosity, %.R 2 = 0.77.where   refers to the porosity at depth z;  is the minimum porosity in the deposition;  is the initial porosity of the rock before deposition; and k is Athy's factor, a constant that depends on the rock type and physical properties.The initial porosity was estimated from pore-recovery research.The sorting features and grain sizes of rocks indicate the general interval of the original porosity.We used the average rock-density test values and the minimum porosity data in Table 3 to shape the initial curve.Next, we calibrated Athy's porosity versus depth curve by using 851 points of measured porosity data from 40 wells in the Changling Depression (Figure 9).We also established a correlation between the porosity and permeability, Equation ( 2) is the formula for K1d (Figure 10).

Calibration of the Model
We established a reference HF trend by comparing the entire basin trend and depression dispersal data (discussed in Section 3.1.6).We may call the reference trend the "plausible HF model" because the trend may be close to the actual conditions, but the actual HF value should fit the formation temperature and maturity of each well [87][88][89].The HF values were homogeneously distributed in the initial HF distribution maps.The HF calibration was based on measured data from 55 wells that covered the main hydrocarbon-bearing structures.The temperature curve that we obtained from the simulation represents only the current HF.The HF trend from 157 Ma to 0 Ma was mainly calibrated by the VR data.The HF trends for individual wells could be precisely restored from the VR and temperature-test data.
In Figure 11a,b, the green curve indicates the temperature and VR of well TS6 as calculated by the input data.We conducted reliable reconstruction of the HF history by using the Easy%Ro algorithm [89,90].The red crosses are the measured well data and the black curves represent the results of the best-fit simulation.The orange line indicates the HF trend that was simulated based on the original input trend.The black line is the suggested trend, which fits the measured data well, as shown in Figure 11c.In the simulation of well TS6, the original HF decreased by 12.18 mW/m 2 to fit the measured data.The blue line represents the smoothed HF trend for drawing the HF map, which slightly deviated from the best fit trend.The smooth process is necessary for creating a practical contour map by decreasing the HF difference between wells.A random simulation was conducted to fit the test data and show the corrected HF value.After interpolation and smoothing, the calibrated HF value distribution maps were created (Figure 12).

Calibration of the Model
We established a reference HF trend by comparing the entire basin trend and depression dispersal data (discussed in Section 3.1.6).We may call the reference trend the "plausible HF model" because the trend may be close to the actual conditions, but the actual HF value should fit the formation temperature and maturity of each well [87][88][89].The HF values were homogeneously distributed in the initial HF distribution maps.The HF calibration was based on measured data from 55 wells that covered the main hydrocarbon-bearing structures.The temperature curve that we obtained from the simulation represents only the current HF.The HF trend from 157 Ma to 0 Ma was mainly calibrated by the VR data.The HF trends for individual wells could be precisely restored from the VR and temperature-test data.
In Figure 11a,b, the green curve indicates the temperature and VR of well TS6 as calculated by the input data.We conducted reliable reconstruction of the HF history by using the Easy%Ro algorithm [89,90].The red crosses are the measured well data and the black curves represent the results of the best-fit simulation.The orange line indicates the HF trend that was simulated based on the original input trend.The black line is the suggested trend, which fits the measured data well, as shown in Figure 11c.In the simulation of well TS6, the original HF decreased by 12.18 mW/m 2 to fit the measured data.The blue line represents the smoothed HF trend for drawing the HF map, which slightly deviated from the best fit trend.The smooth process is necessary for creating a practical contour map by decreasing the HF difference between wells.A random simulation was conducted to fit the test data and show the corrected HF value.After interpolation and smoothing, the calibrated HF value distribution maps were created (Figure 12).After defining the HF based on different scales, namely, horizontal divergence and temporal variations, we ran the entire 3D model simulation.We extracted several individual wells from the results to check the accuracy of the model.The simulated temperature and VR fit the measured data well, as shown in Figure 13.

Source Rock Maturation, Petroleum Generation and Expulsion
We set a pseudo well from the southern fault sag, where the three series of potential source rocks are distributed in depth.Because of the rapid burial speed, organic matter transformed into hydrocarbons shortly after the maturation exceeded the threshold of petroleum generation (Ro = 0.5%) [61].The petroleum generation in the J3h began at 142 Ma, that in K1sh began at 132 Ma, and that in K1y began at 125 Ma.The source rocks entered the oil window shortly after the deposition of each layer, demonstrating rapid sedimentation, and the burial depth exceeded 2000 m before the formations finished deposition (Figure 14a).On the other hand, the high paleo-formation temperature gradient also contributed to the swift maturation.An obvious mitigation phase occurred, during which the TR and VR increased from 129 Ma to 112 Ma, corresponding to the first After defining the HF based on different scales, namely, horizontal divergence and temporal variations, we ran the entire 3D model simulation.We extracted several individual wells from the results to check the accuracy of the model.The simulated temperature and VR fit the measured data well, as shown in Figure 13.After defining the HF based on different scales, namely, horizontal divergence and temporal variations, we ran the entire 3D model simulation.We extracted several individual wells from the results to check the accuracy of the model.The simulated temperature and VR fit the measured data well, as shown in Figure 13.

Source Rock Maturation, Petroleum Generation and Expulsion
We set a pseudo well from the southern fault sag, where the three series of potential source rocks are distributed in depth.Because of the rapid burial speed, organic matter transformed into hydrocarbons shortly after the maturation exceeded the threshold of petroleum generation (Ro = 0.5%) [61].The petroleum generation in the J3h began at 142 Ma, that in K1sh began at 132 Ma, and that in K1y began at 125 Ma.The source rocks entered the oil window shortly after the deposition of each layer, demonstrating rapid sedimentation, and the burial depth exceeded 2000 m before the formations finished deposition (Figure 14a).On the other hand, the high paleo-formation temperature gradient also contributed to the swift maturation.An obvious mitigation phase occurred, during which the TR and VR increased from 129 Ma to 112 Ma, corresponding to the first

Source Rock Maturation, Petroleum Generation and Expulsion
We set a pseudo well from the southern fault sag, where the three series of potential source rocks are distributed in depth.Because of the rapid burial speed, organic matter transformed into hydrocarbons shortly after the maturation exceeded the threshold of petroleum generation (Ro = 0.5%) [61].The petroleum generation in the J 3 h began at 142 Ma, that in K 1 sh began at 132 Ma, and that in K 1 y began at 125 Ma.The source rocks entered the oil window shortly after the deposition of each layer, demonstrating rapid sedimentation, and the burial depth exceeded 2000 m before the formations finished deposition (Figure 14a).On the other hand, the high paleo-formation temperature gradient also contributed to the swift maturation.An obvious mitigation phase occurred, during which the TR and VR increased from 129 Ma to 112 Ma, corresponding to the first denudation of the entire region at the end of the deposition of K 1 y.At 128 Ma, the Huoshiling source rock reached a TR value of 85% and VR value of 2.5% (Figure 14b), which means that the organic matter in the J 3 h transformed into gas before the trap had formed.The increase in generated mass was low after 128 Ma, and hydrocarbon activity tended to stop gradually.Thus, the hydrocarbon expulsion from the J 3 h likely contributed little to the gas accumulation in the Changling Depression.As shown in Figure 14c, the generated mass of K 1 sh experienced two increasing peak stages from 136 Ma to 128 Ma and from 112 Ma to 100 Ma.The TR of the first stage reached 25% and the VR reached 0.9%.The lower seam generated gas, while the mudstone generated oil; however, the early hydrocarbon-expulsion event did not have proper preservation conditions.The TR of the K 1 y and K 1 sh rapidly increased from 112 Ma to 100 Ma, as did the generated mass' growth rate (Figure 14b,c).The K 1 d reservoir and K 1 q seal rock were deposited during this period, and the conditions were suitable for the accumulation of hydrocarbons.
Energies 2019, 12, x FOR PEER REVIEW 16 of 27 denudation of the entire region at the end of the deposition of K1y.At 128 Ma, the Huoshiling source rock reached a TR value of 85% and VR value of 2.5% (Figure 14b), which means that the organic matter in the J3h transformed into gas before the trap had formed.The increase in generated mass was low after 128 Ma, and hydrocarbon activity tended to stop gradually.Thus, the hydrocarbon expulsion from the J3h likely contributed little to the gas accumulation in the Changling Depression.As shown in Figure 14c, the generated mass of K1sh experienced two increasing peak stages from 136 Ma to 128 Ma and from 112 Ma to 100 Ma.The TR of the first stage reached 25% and the VR reached 0.9%.The lower seam generated gas, while the mudstone generated oil; however, the early hydrocarbon-expulsion event did not have proper preservation conditions.The TR of the K1y and K1sh rapidly increased from 112 Ma to 100 Ma, as did the generated mass' growth rate (Figure 14b,c).The K1d reservoir and K1q seal rock were deposited during this period, and the conditions were suitable for the accumulation of hydrocarbons.Two main sources were present: the southwestern sag and the northern sag (Figure 15a-c).The maturity of the northern sag is currently higher than that of the southwestern sag.All three source rocks had entered the oil window, and the maturity in the sag reached 2.0 %Ro.Alongside the dark-mudstone distribution areas (Figure 5), the degree of thermal evolution had reached a high level.Hydrocarbons formed shortly after the strata were deposited; however, the reservoir and seal rocks had not yet formed, and hydrocarbons may have escaped.Any gas that was generated before 100 Ma was not preserved and the rapid increase period may have contributed little to the final accumulation.The slow growth of generated mass after 100 Ma may have been related to the current oil and gas resources in the basin (Figure 14d-f).The magnitude of hydrocarbon generation from the K1y and K1sh source rocks was significantly larger than that from the J3h.The gas production from K1sh was larger because of the deposition of a coal seam.Two main sources were present: the southwestern sag and the northern sag (Figure 15a-c).The maturity of the northern sag is currently higher than that of the southwestern sag.All three source rocks had entered the oil window, and the maturity in the sag reached 2.0 %Ro.Alongside the dark-mudstone distribution areas (Figure 5), the degree of thermal evolution had reached a high level.Hydrocarbons formed shortly after the strata were deposited; however, the reservoir and seal rocks had not yet formed, and hydrocarbons may have escaped.Any gas that was generated before 100 Ma was not preserved and the rapid increase period may have contributed little to the final accumulation.The slow growth of generated mass after 100 Ma may have been related to the current oil and gas resources in the basin (Figure 14d-f).The magnitude of hydrocarbon generation from the K 1 y and K 1 sh source rocks was significantly larger than that from the J 3 h.The gas production from K 1 sh was larger because of the deposition of a coal seam.

Hydrocabon Accumulation in the Main Reservoir: Denglouku Formation
The modeled oil-and gas-migration paths and accumulation in the four main gas fields, namely, Haerjin, Shuangtuozi, Dalaoyefu and Fulongquan, are shown below (Figure 16).The green lines are the migration flow paths of the oil, and the red lines indicate the gas-migration paths.The migration rules of the entire basin flowed from the depression to the margin because of differences in the fluid potential.The hydrocarbon streams vertically migrated to K1d through opening faults and then horizontally migrated to structural highs.

Hydrocabon Accumulation in the Main Reservoir: Denglouku Formation
The modeled oil-and gas-migration paths and accumulation in the four main gas fields, namely, Haerjin, Shuangtuozi, Dalaoyefu and Fulongquan, are shown below (Figure 16).The green lines are the migration flow paths of the oil, and the red lines indicate the gas-migration paths.The migration rules of the entire basin flowed from the depression to the margin because of differences in the fluid potential.The hydrocarbon streams vertically migrated to K 1 d through opening faults and then horizontally migrated to structural highs.

Hydrocabon Accumulation in the Main Reservoir: Denglouku Formation
The modeled oil-and gas-migration paths and accumulation in the four main gas fields, namely, Haerjin, Shuangtuozi, Dalaoyefu and Fulongquan, are shown below (Figure 16).The green lines are the migration flow paths of the oil, and the red lines indicate the gas-migration paths.The migration rules of the entire basin flowed from the depression to the margin because of differences in the fluid potential.The hydrocarbon streams vertically migrated to K1d through opening faults and then horizontally migrated to structural highs.The petroleum-migration rules revealed four main migration-accumulation units, as shown in Figure 17: a central uplift zone, a southern fault sag zone, a western slope zone, and a northeastern slope zone.The central uplift zone included a southeastern slope zone.The Haerjin structure is located in the central uplift zone, and the other three gas fields are distributed in the southeastern slope zone.We found few migration paths that could have formed the Haerjin gas field.The faults acted as a channel for vertical hydrocarbon migration and a sealing border to prevent the lateral leakage of gas.The source rocks that contributed to the gas zone that is marked by well CS1 should be located in the K 1 y and K 1 sh of the southern sag.The eastern portion of the K 1 sh source rock may have also generated hydrocarbons and moved to the Haerjin structure (Figure 5b,e,h).The Dalaoyefu and Shuangtuozi structures are located above the eastern source rock, and hydrocarbons that formed in this area migrated to nearby traps.The source rocks reached high maturity; the VR was approximately 2.0%, reaching the gas-generation peak (Figure 15).Additionally, hydrocarbons in surrounding areas would have amassed in structural highs.The Fulongquan structure is located near the margin of the Changling Depression, and the underlying source rock has the potential for abundant gas generation.The thickness of the source rocks was approximately 400 m, and the HI and TOC values were suitable for abundant gas generation, although this area is not large (Figure 5).The different tectonic belts had diverse characteristics in terms of hydrocarbon-charging rules.We will discuss these main gas-bearing structures in Section 5 below.The petroleum-migration rules revealed four main migration-accumulation units, as shown in Figure 17: a central uplift zone, a southern fault sag zone, a western slope zone, and a northeastern slope zone.The central uplift zone included a southeastern slope zone.The Haerjin structure is located in the central uplift zone, and the other three gas fields are distributed in the southeastern slope zone.We found few migration paths that could have formed the Haerjin gas field.The faults acted as a channel for vertical hydrocarbon migration and a sealing border to prevent the lateral leakage of gas.The source rocks that contributed to the gas zone that is marked by well CS1 should be located in the K1y and K1sh of the southern sag.The eastern portion of the K1sh source rock may have also generated hydrocarbons and moved to the Haerjin structure (Figure 5b,e,h).The Dalaoyefu and Shuangtuozi structures are located above the eastern source rock, and hydrocarbons that formed in this area migrated to nearby traps.The source rocks reached high maturity; the VR was approximately 2.0%, reaching the gas-generation peak (Figure 15).Additionally, hydrocarbons in surrounding areas would have amassed in structural highs.The Fulongquan structure is located near the margin of the Changling Depression, and the underlying source rock has the potential for abundant gas generation.The thickness of the source rocks was approximately 400 m, and the HI and TOC values were suitable for abundant gas generation, although this area is not large (Figure 5).The different tectonic belts had diverse characteristics in terms of hydrocarbon-charging rules.We will discuss these main gas-bearing structures in Section 5 below.Oil and gas migrated through the enormous regional faults in the southern fault sag zone.Leakage from the source rock was a significant problem and acted as a restraining factor for gas accumulation.Several sets of fault traps could have been present, but the preservation conditions were unstable because of frequent tectonic activity (Figure 17).Hydrocarbons moved horizontally to the NW in K1d in the western slope zone.Hydrocarbons that formed in the source rock in the central depression moved to the western slope and may have accumulated in the border of the basin and fault traps.The generated mass was large, but the reservoir layer was nearly flat in that area and could not adequately hold gas.The model was set as a closed-border basin because the Changling Depression is surrounded by border faults, which cut through all the Cretaceous layers.Thus, hydrocarbons that migrated to the edge of the depression would be blocked by these faults.Enormous oil and gas accumulation occurred in the NE slope zone, where the migration rules were similar to those in the western slope zone.A structural high was present in the NE slope zone, and hydrocarbon generation in the northern source rock filled up the trap in this location.Oil and gas migrated through the enormous regional faults in the southern fault sag zone.Leakage from the source rock was a significant problem and acted as a restraining factor for gas accumulation.Several sets of fault traps could have been present, but the preservation conditions were unstable because of frequent tectonic activity (Figure 17).Hydrocarbons moved horizontally to the NW in K 1 d in the western slope zone.Hydrocarbons that formed in the source rock in the central depression moved to the western slope and may have accumulated in the border of the basin and fault traps.The generated mass was large, but the reservoir layer was nearly flat in that area and could not adequately hold gas.The model was set as a closed-border basin because the Changling Depression is surrounded by border faults, which cut through all the Cretaceous layers.Thus, hydrocarbons that migrated to the edge of the depression would be blocked by these faults.Enormous oil and gas accumulation occurred in the NE slope zone, where the migration rules were similar to those in the Energies 2019, 12, 1043 19 of 27 western slope zone.A structural high was present in the NE slope zone, and hydrocarbon generation in the northern source rock filled up the trap in this location.

Definition of Key Hydrocarbon-Accumulation Moments
Fluid inclusions, or microscopic vesicles of liquid and gas, are trapped in growing crystals.Little addition of irrelevant matter or reduction of the original substances occurs after the formation of a fluid inclusion.Hydrocarbon and liquid inclusions can be trapped concurrently during hydrocarbon migration [91].Therefore, these inclusions offer relatively accurate information regarding the temperature of mineral growth at the time of the fluid migration and entrapment, which has been widely used in petroleum geology [1,67,[92][93][94].
Prior studies have measured the homogenization temperatures of fluid inclusions from 19 reservoir samples from 10 wells in K 1 d in the Changling Depression [12,95].The frequency histogram of the homogenization temperatures of inclusions shows different distributions among the four gas-bearing structures, as discussed in Section 4.3 (Figure 18A-b, B-b, C-b and D-b).The trapping temperatures of different structures varied.The different burial depths and temperature gradients may have created differences in the inclusions' homogenization temperatures, as did the unique filling timing.We examined individual wells near the highest point of each petroliferous trap and extracted their 1D models.When we projected the dominant iso-temperature lines and strata horizons, these features intersected at several points.The crossing point in K 1 d, the main reservoir layer, marks when and where hydrocarbons moved through or in the sandstone.Therefore, we could reveal the critical period of the charging event, which is an important component of reconstructing petroleum system elements.

Definition of Key Hydrocarbon-Accumulation Moments
Fluid inclusions, or microscopic vesicles of liquid and gas, are trapped in growing crystals.Little addition of irrelevant matter or reduction of the original substances occurs after the formation of a fluid inclusion.Hydrocarbon and liquid inclusions can be trapped concurrently during hydrocarbon migration [91].Therefore, these inclusions offer relatively accurate information regarding the temperature of mineral growth at the time of the fluid migration and entrapment, which has been widely used in petroleum geology [1,67,[92][93][94].
Prior studies have measured the homogenization temperatures of fluid inclusions from 19 reservoir samples from 10 wells in K1d in the Changling Depression [12,95].The frequency histogram of the homogenization temperatures of inclusions shows different distributions among the four gas-bearing structures, as discussed in Section 4.3 (Figure 18A-b, B-b, C-b and D-b).The trapping temperatures of different structures varied.The different burial depths and temperature gradients may have created differences in the inclusions' homogenization temperatures, as did the unique filling timing.We examined individual wells near the highest point of each petroliferous trap and extracted their 1D models.When we projected the dominant iso-temperature lines and strata horizons, these features intersected at several points.The crossing point in K1d, the main reservoir layer, marks when and where hydrocarbons moved through or in the sandstone.Therefore, we could reveal the critical period of the charging event, which is an important component of reconstructing petroleum system elements.According to Figure 18A-b, the homogenization temperatures of the fluid inclusions were mainly distributed from 110 °C to 140 °C.In Figure 18A-a, we plotted iso-temperature lines at 110 °C and 140 °C, and the iso-temperature lines and K1d stratum are encircled in the red region, indicating the accumulation period and depth.Two stars mark the beginning and end of the main reservoir-formation stage from 80 Ma to 66 Ma, corresponding to the Middle Nengjiang Formation to Early Mingshui Formation deposition.Additionally, we could verify the result in Figure 18A-c.The red line represents the formation temperature, while the blue line represents the petroleum-saturation change in K1d.The petroleum saturation was calculated from the migration and accumulation simulations.The red area was limited by certain extreme temperature values and matched the width of the blue area, which was defined by the petroleum-saturation increase period.Thus, the moment of the petroleum-charging event conformed to the time that aqueous inclusions were captured.Long-term hydrocarbon accumulation occurred in the Haerjin structure from 80 Ma to 66 Ma.An original charging event also occurred from 100 Ma to 88 Ma, when the hydrocarbon saturation rapidly rose but then subsequently dropped because the trap had not yet formed (Figure 18A-c).18A-a, we plotted iso-temperature lines at 110 • C and 140 • C, and the iso-temperature lines and K 1 d stratum are encircled in the red region, indicating the accumulation period and depth.Two stars mark the beginning and end of the main reservoir-formation stage from 80 Ma to 66 Ma, corresponding to the Middle Nengjiang Formation to Early Mingshui Formation deposition.Additionally, we could verify the result in Figure 18A-c.The red line represents the formation temperature, while the blue line represents the petroleum-saturation change in K 1 d.The petroleum saturation was calculated from the migration and accumulation simulations.The red area was limited by certain extreme temperature values and matched the width of the blue area, which was defined by the petroleum-saturation increase period.Thus, the moment of the petroleum-charging event conformed to the time that aqueous inclusions were captured.Long-term hydrocarbon accumulation occurred in the Haerjin structure from 80 Ma to 66 Ma.An original charging event also occurred from 100 Ma to 88 Ma, when the hydrocarbon saturation rapidly rose but then subsequently dropped because the trap had not yet formed (Figure 18A-c 18B-a shows the two hydrocarbon-charging periods.The first charging event occurred between 106 Ma and 101 Ma, which is marked in the left-hand triangle in Figure 18B-c.The petroleum saturation rose at almost the same period until 100 Ma, when the saturation dropped to 0% over the following 6 Ma, which marks petroleum leakage.After 78 Ma, the petroleum saturation rose significantly, corresponding to the temperature range.The hydrocarbon charged over the following 10 Ma and dropped during the uplifting and erosion period.In fact, the petroleum saturation was relatively low in the Shuangtuozi structure, and the present-day hydrocarbon aggregation was mainly derived from the lateral filling from 23 Ma to 0 Ma.

Dalaoyefu Structure: Well L14
The samples from the Dalaoyefu structure recorded a dominant homogenization temperature from 90 • C to 120 • C (Figure 18C-b).A long-term hydrocarbon-charging event occurred from 101 Ma to 67 Ma (Figure 18C-a).The K 1 d is not particularly thick but did record high petroleum saturation because of continuous gas filling and lower leakage during the structural-inversion period.The maximum saturation was approximately 90% at 66 Ma (Figure 18C-c).The petroleum-charging period lasted from 112 Ma to 66 Ma according to the saturation change, and rapid growth began at 100 Ma, which matched the temperature interval.

Fulongquan Structure: Well F3
The model covered a portion of the Fulongquan structure, and the homogenization-temperature analysis results were mainly distributed from 70 • C to 110 • C (Figure 18D-b).The depth of K 1 d in the Fulongquan structure is approximately 1000 m, but the HF in this area was relatively higher (Figure 12).The charging period in K 1 d lasted from 107 Ma to 67 Ma (Figure 18D-a).The evolution of the petroleum saturation continuously changed from 112 Ma to 66 Ma (Figure 18D-c).First, the saturation increased from 112 Ma to 84 Ma during hydrocarbon filling.Then, the saturation decreased from 84 Ma to 74 Ma because of leakage and then rose again between 74 Ma and 66 Ma.During the long-term denudation from 66 Ma to 23 Ma, the petroleum saturation continuously dropped until the start of Paleogene deposition.At present, the saturation has slightly increased but is relatively lower than that of the Dalaoyefu structure.

Migration and Accumulation Rules in the Changling Depression
We extracted a 2D section model by using a section line in the Changling Depression (Figure 17) to summarize the overall rules of the migration and accumulation of hydrocarbons in tight sandstone in the Changling Depression (Figure 19).In the Haerjin and Dalaoyefu structures, humic-type pyrolysis gas that formed from source rock in the K 1 y and K 1 sh migrated through deep faults and accumulated in the K 1 d reservoir.The gas accumulation in Well CS1, Well CS2 and Well L14 mainly followed the vertical-migration rule.The gas pools in the Shuangtuozi and Fulongquan structures originated from hydrocarbons that formed in K 1 y and migrated horizontally in K 1 d.Gas accumulation in the slope areas, including Well TS1 and Well F3, primarily abided by the horizontal-migration rule.

Conclusions
The gas accumulation zones in the Changling Depression obtained using the 3D basin modelling approach matched the discovered gas fields.We could pay more attention to the potential petroleum accumulation areas in further exploration.The hydrocarbon-generation center is in the southwestern sag and the northern sag of the Changling Depression, and K1sh and K1y are the major source rocks, which have high gas-generation intensity, but the peak petroleum-generation period (112-100 Ma) was earlier than the formation of the trap.The storage configuration is the critical condition of the hydrocarbon accumulation.
Homogenization-temperature analysis of fluid inclusions indicated two sets of critical moments of hydrocarbon accumulation.The Dalaoyefu and Fulongquan structures experienced long-term hydrocarbon accumulation from 100 Ma to 67 Ma.The Haerjin and Shuangtuozi structures followed the lateral reservoir-formation stage, and hydrocarbon filling occurred between 80 Ma and 66 Ma.The Shuangtuozi structure had two phases of accumulation, but the first period did not contribute to the present gas accumulation.
A combination of inclusion analysis and hydrocarbon-accumulation simulations indicated different stages of reservoir formation, gas leakage and gas-pool reformation.The major gas-bearing structures initially formed during the deposition of the Early K1q (110 Ma) and then underwent varying degrees of leakage.By the middle of the deposition of the Nenjiang Formation (80 Ma), the structures had developed into the primary gas-accumulation reservoir.Long-term uplifting and regional erosion during the later deposition of the Nenjiang and Mingshui formations caused petroleum losses.Gas accumulation steadily increased after denudation (23 Ma) and formed sets of secondary gas reservoirs.The remaining kerogen in the K1y and K1sh supported the gas intensity during the Neogene.
The Changling Depression could be divided into four migration and accumulation units, and the distribution of the source rocks controlled the accumulation locations of hydrocarbons.Faults acted as major channels for hydrocarbon migration, and hydrocarbons could migrate horizontally into the reservoirs over short distances.
Author Contributions: Conceptualization by J.Z., J.G., model construction by J.G., source rock analysis by Y.L., data processing by Z.S.; J.G. and J.Z. contributes to the writing of the manuscript.

Conclusions
The gas accumulation zones in the Changling Depression obtained using the 3D basin modelling approach matched the discovered gas fields.We could pay more attention to the potential petroleum accumulation areas in further exploration.The hydrocarbon-generation center is in the southwestern sag and the northern sag of the Changling Depression, and K 1 sh and K 1 y are the major source rocks, which have high gas-generation intensity, but the peak petroleum-generation period (112-100 Ma) was earlier than the formation of the trap.The storage configuration is the critical condition of the hydrocarbon accumulation.
Homogenization-temperature analysis of fluid inclusions indicated two sets of critical moments of hydrocarbon accumulation.The Dalaoyefu and Fulongquan structures experienced long-term hydrocarbon accumulation from 100 Ma to 67 Ma.The Haerjin and Shuangtuozi structures followed the lateral reservoir-formation stage, and hydrocarbon filling occurred between 80 Ma and 66 Ma.The Shuangtuozi structure had two phases of accumulation, but the first period did not contribute to the present gas accumulation.
A combination of inclusion analysis and hydrocarbon-accumulation simulations indicated different stages of reservoir formation, gas leakage and gas-pool reformation.The major gas-bearing structures initially formed during the deposition of the Early K 1 q (110 Ma) and then underwent varying degrees of leakage.By the middle of the deposition of the Nenjiang Formation (80 Ma), the structures had developed into the primary gas-accumulation reservoir.Long-term uplifting and regional erosion during the later deposition of the Nenjiang and Mingshui formations caused petroleum losses.Gas accumulation steadily increased after denudation (23 Ma) and formed sets of secondary gas reservoirs.The remaining kerogen in the K 1 y and K 1 sh supported the gas intensity during the Neogene.
The Changling Depression could be divided into four migration and accumulation units, and the distribution of the source rocks controlled the accumulation locations of hydrocarbons.Faults acted as major channels for hydrocarbon migration, and hydrocarbons could migrate horizontally into the reservoirs over short distances.

Figure 1 .
Figure 1.Location and geographic setting of the study area.(a) General location map, with the borders and tectonic divisions of the Songliao Basin and Changling Depression modified after Chen et al. [13].(b) Major structures and well locations in the Changling Depression.

Figure 1 .
Figure 1.Location and geographic setting of the study area.(a) General location map, with the borders and tectonic divisions of the Songliao Basin and Changling Depression modified after Chen et al. [13].(b) Major structures and well locations in the Changling Depression.

Figure 3 .
Figure 3. 3D model view and present-day dimensions of the Changling Depression.The table describes the general stratigraphy, alongside the PSE and the lithology types that were assigned for the model.Fm. = formation, OR = overburden rock, UR = underburden rock."Igneous rock mixture" represents a mix of rhyolite, tuff and granite.

Figure 3 .
Figure 3. 3D model view and present-day dimensions of the Changling Depression.The table describes the general stratigraphy, alongside the PSE and the lithology types that were assigned for the model.Fm. = formation, OR = overburden rock, UR = underburden rock."Igneous rock mixture" represents a mix of rhyolite, tuff and granite.

Figure 4 .
Figure 4. Source rock evaluation plots of the main source rocks in the deep Changling Depression.(a) Variation in Rock-Eval S1 + S2 with the total organic carbon content (TOC) .(b) Kerogen type and maturation definition with a van Krevelen diagram.(c) Hydrogen index (HI) with pyrolysis peak temperature (Tmax) for samples from K1y.(d) HI with Tmax for samples from K1sh.

Figure 4 .
Figure 4. Source rock evaluation plots of the main source rocks in the deep Changling Depression.(a) Variation in Rock-Eval S 1 + S 2 with the total organic carbon content (TOC).(b) Kerogen type and maturation definition with a van Krevelen diagram.(c) Hydrogen index (HI) with pyrolysis peak temperature (Tmax) for samples from K 1 y.(d) HI with Tmax for samples from K 1 sh.

Figure 6 .
Figure 6.Kinetic models of the source rocks of the Changling Depression in this study.(a) Kinetic model of the shale based on the Burnham TIII model [69]; (b) kinetic model of the coal seam based on the Pepper and Corvi TIIIH(DE) model [75].

Figure 6 .
Figure 6.Kinetic models of the source rocks of the Changling Depression in this study.(a) Kinetic model of the shale based on the Burnham TIII model [69]; (b) kinetic model of the coal seam based on the Pepper and Corvi TIIIH(DE) model [75].
Energies 2019, 12, x FOR PEER REVIEW 11 of 27 directions, and the faults formed during the deposition of the J3h until the deposition of K1d.We imported the interpreted fault lines (intersections on each surface) into the model and created model faults.The border of the Changling Sag mostly consisted of boundary fractures, and we assumed that the basin sides were closed.

Figure 7 .
Figure 7. Fault distribution in the Changling Depression.(a) Fault intersections on each layer surface.(b) Distribution of the gridded faults of the Changling Depression in the model.

Figure 7 .
Figure 7. Fault distribution in the Changling Depression.(a) Fault intersections on each layer surface.(b) Distribution of the gridded faults of the Changling Depression in the model.

Figure 8 .
Figure 8.General heat flow (HF) history that was calculated from the surface HF and basement HF of the Songliao Basin based on the research of Li, 1995 [83].

Figure 8 .
Figure 8.General heat flow (HF) history that was calculated from the surface HF and basement HF of the Songliao Basin based on the research of Li, 1995 [83].

Figure 10 .
Figure 10.Correlation between the porosity and permeability in K1d reservoir.

Figure 9 .
Figure 9. Porosity versus depth curve (red curve) as calibrated by measured data (black dots).(a) Porosity calibration curve of K 1 d.(b) Calibration curve of K 1 y.(c) Calibration curve of K 1 sh.

Figure 10 .
Figure 10.Correlation between the porosity and permeability in K1d reservoir.Figure 10.Correlation between the porosity and permeability in K 1 d reservoir.

Figure 10 .
Figure 10.Correlation between the porosity and permeability in K1d reservoir.Figure 10.Correlation between the porosity and permeability in K 1 d reservoir.

Figure 11 .
Figure 11.Results of HF calibration (Well TS6).(a) Single well temperature simulation result using original HF (green curve) and the best-fit HF (black curve), calibrated by bore-hole temperatures (red crosses).(b) Single-well vitrinite reflectance (VR) simulation result using original HF and best-fit HF, calibrated by measured Ro (red crosses).(c) Optimized HF among the original HF, best-fit trend and calibrated trend after implementing the map-smoothing method.

Figure 11 .
Figure 11.Results of HF calibration (Well TS6).(a) Single well temperature simulation result using original HF (green curve) and the best-fit HF (black curve), calibrated by bore-hole temperatures (red crosses).(b) Single-well vitrinite reflectance (VR) simulation result using original HF and best-fit HF, calibrated by measured Ro (red crosses).(c) Optimized HF among the original HF, best-fit trend and calibrated trend after implementing the map-smoothing method.

Figure 12 .
Figure 12.Heat-flow distribution definition.(a) Calibrated HF distribution map for the present day (0 Ma).(b) Calibrated HF distribution map for 112 Ma, after K1d deposition.

Figure 13 .
Figure 13.Validation of the simulation results.The temperature and vitrinite-reflectance measurement data (red cross) and simulation results of wells CS1 (a,b) and SN301 (c,d) showed good fit.

Figure 12 .
Figure 12.Heat-flow distribution definition.(a) Calibrated HF distribution map for the present day (0 Ma).(b) Calibrated HF distribution map for 112 Ma, after K 1 d deposition.

Figure 13 .
Figure 13.Validation of the simulation results.The temperature and vitrinite-reflectance measurement data (red cross) and simulation results of wells CS1 (a,b) and SN301 (c,d) showed good fit.

Figure 13 .
Figure 13.Validation of the simulation results.The temperature and vitrinite-reflectance measurement data (red cross) and simulation results of wells CS1 (a,b) and SN301 (c,d) showed good fit.

Figure 14 .
Figure 14.Burial history, maturity evolution and hydrocarbon-generation mass of the main source rocks in the Changling Depression.(a) Burial history of the three main source rocks: Yingcheng (red), Shahezi (blue) and Huoshiling (green) Formations.(b) Evolution of the VR and TR throughout geological time (pseudo well location: X: 21568.264km, Y: 4927.123km).(c) Amount (in mass) of hydrocarbons that were generated from the three main source rocks over time (from the entire study area).

Figure 14 .
Figure 14.Burial history, maturity evolution and hydrocarbon-generation mass of the main source rocks in the Changling Depression.(a) Burial history of the three main source rocks: Yingcheng (red), Shahezi (blue) and Huoshiling (green) Formations.(b) Evolution of the VR and TR throughout geological time (pseudo well location: X: 21568.264km, Y: 4927.123km).(c) Amount (in mass) of hydrocarbons that were generated from the three main source rocks over time (from the entire study area).

Figure 15 .
Figure 15.Maturity distribution and hydrocarbon generation from the main source rocks in the Changling Depression.(a) Present day maturity maps for K1y.(b) Present day maturity maps for K1sh.(c) Present day maturity maps for J3h.(d) Generation mass of gas (red curve) and oil (green line) from K1y.(e) Generation mass from K1sh.(f) Generation mass from J3h.

Figure 15 .
Figure 15.Maturity distribution and hydrocarbon generation from the main source rocks in the Changling Depression.(a) Present day maturity maps for K 1 y.(b) Present day maturity maps for K 1 sh.(c) Present day maturity maps for J 3 h.(d) Generation mass of gas (red curve) and oil (green line) from K 1 y.(e) Generation mass from K 1 sh.(f) Generation mass from J 3 h.

Energies 2019 , 27 Figure 15 .
Figure 15.Maturity distribution and hydrocarbon generation from the main source rocks in the Changling Depression.(a) Present day maturity maps for K1y.(b) Present day maturity maps for K1sh.(c) Present day maturity maps for J3h.(d) Generation mass of gas (red curve) and oil (green line) from K1y.(e) Generation mass from K1sh.(f) Generation mass from J3h.

Figure 16 .
Figure 16.Current hydrocarbon-migration paths and accumulation zones in K1d.Figure 16.Current hydrocarbon-migration paths and accumulation zones in K 1 d.

Figure 16 .
Figure 16.Current hydrocarbon-migration paths and accumulation zones in K1d.Figure 16.Current hydrocarbon-migration paths and accumulation zones in K 1 d.

Figure 18 .
Figure 18.Petroleum-charging moment analysis for four key wells that represent the gas fields: Well CS1 (A), Well TS1 (B), Well L4 (C) and Well F3 (D).(a) Hydrocarbon-charging period according to the homogenization temperatures of fluid inclusions and the burial history curve.(b) Histogram of the homogenization temperatures of the fluid inclusions from Meng, 2009 [95].(c) Petroleum saturation and formation-temperature trend of K1d, which were used to distinguish the charging event.5.1.1.Haerjin Structure: Well CS1

Figure 18 .
Figure 18.Petroleum-charging moment analysis for four key wells that represent the gas fields: Well CS1 (A), Well TS1 (B), Well L4 (C) and Well F3 (D).(a) Hydrocarbon-charging period according to the homogenization temperatures of fluid inclusions and the burial history curve.(b) Histogram of the homogenization temperatures of the fluid inclusions from Meng, 2009 [95].(c) Petroleum saturation and formation-temperature trend of K 1 d, which were used to distinguish the charging event.5.1.1.Haerjin Structure: Well CS1 According to Figure 18A-b, the homogenization temperatures of the fluid inclusions were mainly distributed from 110 • C to 140 • C. In Figure18A-a, we plotted iso-temperature lines at 110 • C and 140 • C, and the iso-temperature lines and K 1 d stratum are encircled in the red region, indicating the accumulation period and depth.Two stars mark the beginning and end of the main reservoir-formation stage from 80 Ma to 66 Ma, corresponding to the Middle Nengjiang Formation to Early Mingshui Formation deposition.Additionally, we could verify the result in Figure18A-c.The red line represents the formation temperature, while the blue line represents the petroleum-saturation change in K 1 d.The petroleum saturation was calculated from the migration and accumulation simulations.The red area was limited by certain extreme temperature values and matched the width of the blue area, which was defined by the petroleum-saturation increase period.Thus, the moment of the petroleum-charging event conformed to the time that aqueous inclusions were captured.Long-term hydrocarbon accumulation occurred in the Haerjin structure from 80 Ma to 66 Ma.An original charging event also occurred from 100 Ma to 88 Ma, when the hydrocarbon saturation rapidly rose but then subsequently dropped because the trap had not yet formed (Figure18A-c).

Figure 19 .
Figure 19.The cross section through wells, which shows the generalized hydrocarbon migration and accumulation model in K1d in the Changling Depression (section location shown in Figure 17).

Funding:
This paper was supported by the Major National R&D Projects of China (No.2016ZX05027-001-006) and Fundamental Research Funds for the Central Universities (No. 2015KJJCB11).

Figure 19 .
Figure 19.The cross section through wells, which shows the generalized hydrocarbon migration and accumulation model in K 1 d in the Changling Depression (section location shown in Figure 17).

Table 1 .
Facies settings for PetroMod ® software, including ages, formations, PSE, lithology settings with thermal parameters and mechanical properties.The Athy's factors are discussed in Section 3.1.7.

Table 2 .
Erosion thicknesses of individual wells in the Changling Depression (unit: m).

Table 2 .
Erosion thicknesses of individual wells in the Changling Depression (unit: m).

Table 3 .
Parameters that were used in Athy's law to calibrate the porosity versus depth curve.

Table 3 .
Parameters that were used in Athy's law to calibrate the porosity versus depth curve.

Table 3 .
Parameters that were used in Athy's law to calibrate the porosity versus depth curve.