Base level changes based on Basin Filling Modelling: a case study from the Paleocene Lishui Sag, East China Sea Basin

Estimation of base level changes in geological records is an important topic for petroleum geologists. Taking the Paleocene Upper Lingfeng Member of Lishui Sag as an example, this paper conducted a base level reconstruction based on Basin Filling Modelling (BFM). The reconstruction was processed on the ground of a previously interpreted seismic stratigraphic framework with several assumptions and simplification. The BFM is implemented with a nonlinear diffusion equation solver written in R coding that excels in shallow marine stratigraphic simulation. The modeled results fit the original stratigraphy very well. The BFM is a powerful tool for reconstructing the base level, and is an effective way to check the reasonableness of previous interpretations. Although simulation solutions may not be unique, the BFM still provides us a chance to gain some insights into the mechanism and dynamic details of the stratigraphy of sedimentary basins.


Introduction
Petroleum exploration is committed to identifying the spatial and temporal distribution of petroleum system elements (such as source, reservoir and transportation, etc.) to identify more oil and gas resources (Hu et al. 2016Li et al. 2020;Zhu et al. 2020). However, the diversity of the underground world and the nature of its difficulty to observe directly surely increase the complexity of this work. Following sedimentary facies model, the advent of sequence stratigraphy has brought another revolution to this industry and offered a powerful methodology to analyze and characterize the underground space (Yang et al. 2018;Gao et al. 2018). Although many schools still exist debating on some concepts of sequence stratigraphy, base level is a universally accepted term that denotes a driving force ensemble in forming stratigraphic sequences. However, ambiguity may exist when talking about the definition of base level. To avoid misunderstanding, here we consider the base level definition below. The base level approximates the sea level and meets the fluvial profile equilibrium at the shoreline (Catuneanu 2006;Miall 2016). It is different from the "stratigraphic base level" of Cross and Lessenger (1998) that has a continental part. Comprehensively reflecting the effects of eustasy, tectonics and climate, base level is usually regarded as a representative term of allogenic controls. Since allogenic controls provide basin-scale comparable information for facies prediction, it makes base level analysis extremely important in both the academy and industry Catuneanu and Zecchin 2013;Qayyum et al. 2015;van der Meer et al. 2017).
In early ages, base level was deemed almost equivalent to the global sea level. Many efforts have been made to derive the global sea level curve in geological time scale since 1970s (e.g., Haq et al. 1987;Miller et al. 2005;Xu et al. 2006;Müller et al. 2008). With the deepening of Edited by Jie Hao understanding, tectonic movements were included into the analysis of base level in the field of sequence stratigraphy. The global sea level has gradually been replaced by relative sea level (RSL) (e.g., Hunt and Tucker 1992;Christie-Blick and Driscoll 1995), as it defines the joint effect of global sea level and tectonism. The RSL is an efficient tool to encompass the uncertainties in distinguishing the eustacy and the tectonism from stratigraphic records, and it could also approximate to the base level with some assumptions (the climate or environmental energy flux factors are neglected) (Miall 2016).
The base level has long been represented by a single curve indicating that the globle sea level value is the same everywhere at a given moment. While it is not true for the RSL or the base level, because their values can vary locally due to differential subsidence Qayyum et al. 2017). This implies that the RSL and the base level cannot be expressed by just one curve as the base level change is sensitive to the location of the observation point. Hence, reconstructing the complete base level is both a static and geomorphic problem, which poses challenges to the traditional methods (Jerolmack and Paola 2010;Li et al. 2016).
Basin Filling Modelling (BFM, sometimes also called forward stratigraphic modelling) gives a new light to this issue. It can forwardly simulate the stratigraphy with parameters, and it can also inversely estimate parameters through iterated simulations. Although we cannot completely eliminate the non-uniqueness of the solutions, we can still narrow the range of the possible scenarios through the BFM. The models for siliciclastic stratigraphy can generally be classified into 3 types, i.e., geometric models, hydraulic models and diffusion-based models (Paola 2000;Burgess et al. 2006;Sacchi et al. 2015). The first type is relatively simple and can directly show the geometric patterns of the objects, while it lacks the ability to reveal the physical processes (Strobel et al. 1989;Kendall et al. 1991;Li et al. 2018). Hydraulic and Diffusion-based models are also called dynamic models, for they are process-based and need to solve governing equations (Kaufman et al. 1991;Rivenaes 1992;Griffiiths et al. 2001;Hutton and Syvitski 2008). Hydraulic models solve degraded Navier-Stokes equations and are good at modeling hydraulic-scale (at relatively smaller scales, e.g., river avulsion) processes. Diffusion-based models solve diffusionrelated equations and are good at modeling geomorphologyscale (at relative larger scales, e.g., clinoform formation) processes. Among these three models, the diffusion-based models suit this study best, since base level reconstruction is a process-based problem while it does not need to deal with the details below a clinoform scale.
The study area Lishui Sag is one of the subordinate structural units of East China Sea Basin (ECSB), offshore east China. The study interval Upper Lingfeng Member (E 1 lU) is an important oil-bearing succession in Lishui Sag that has accumulated very rich exploration data. Despite some commercial successes, little consensus has been made regarding the base level changes during the deposition of the Paleocene successions. In this paper, taking the Paleocene E 1 lU of Lishui Sag as an example, we conduct a base level reconstruction based on a diffusion-based BFM model. It aims at not only revealing some base level details to the study area, but also gaining some insights into the mechanism and dynamic details of the stratigraphic formation.

Geological settings
The East China Sea Shelf Basin (ECSB), located in the southeast margin of the Eurasian Plate and above the South China Plate, is a part of the western Pacific tectonic system, ranging from 120° 50′ E to 129° 00′ E in longitude and from 25° 22′ N to 33° 38′ N in latitude ( Fig. 1a, b). The basin has an NNE direction and could be divided into East Depression belt, Central Uplift belt, and West Depression belt, covering an area of 26.7 × 104 km 2 . Inside the West Depression belt are several sags, including Lishui Sag (Fig. 1b, c). Lishui Sag is located in the southwest part of the entire ECSB, covering an area of 14,600 km 2 . Before the deposition of the Paleocene Mingyuefeng Formation (E 1 m), the Lingfeng Uplift divided the sag into two sub-sags: Lishui West Subsag and Lishui East Subsag, both showing half-graben characteristics (Fig. 1d).

Tectonism
Three events have influenced the ECSB most since Late Triassic in terms of tectonic evolution. They are (1) the forming of ancient Asia Continent, (2) The Pacific Plate subducting under Eurasian Plate, and (3) Indian Plate subducting under Eurasian Plate. These events could be rearranged into 4 specific tectonic phases (Fig. 2). From the Late Triassic to Present, they are: (1) Syn-rift phase (Late Cretaceous to Paleocene), during which the basin had undergone intensive movements that created many fault-depression structures; (2) Post-rift phase (Early Eocene to Late Eocene), during which extensional movements reduced; (3) Uplift phase (Late Eocene to Oligocene), during which the basin had been subjected to uplift and following erosion; and (4) Regional subsidence phase (Miocene-present), during which the basin have subsided as a whole Li et al. 2019).

Stratigraphy
The Cenozoic strata in the ECSB basin are well developed, in contrast with the poorly developed Mesozoic strata. According to a generally accepted dividing plan Sun et al. 2019), the Paleogene to Neogene strata could be divided into 9 formations. From old to new, they are Paleocene Yueguifeng Formation (E 1 y), Lingfeng Formation (E 1 l), Mingyuefeng Formation (E 1 m), Eocene Oujiang Formation, Wenzhou Formation, Pinghu Formation, Oligocene Huagang Formation, Miocene Longjing Formation, Yuquan Formation, Liulang Formation and Pliocene Santan Formation (Fig. 2).
The interval of interest is the Paleocene strata, deposited during the syn-rift stage. The lowest Yueguifeng Formation, overlying on Cretaceous Shimentan Formation (K 2 s) via an angular unconformity, has a thickness ranging from 0 to 349 m (Fig. 2). The E 1 y is of mainly lacustrine deposition, characterized by light to dark gray sandstones, siltstones, mudstones as well as their intercalations, which differs a lot from the underlying brown to red-brown deposits. Above the E 1 y is the Lingfeng Formation that could be further divided into 2 members: Lower Lingfeng Member (E 1 lL) and Upper Lingfeng Member (E 1 lU). The E 1 lU member has very good trap conditions and is believed to have a petroleum reservoir potential. The depositional environment of the E 1 lL and the E 1 lU is interpreted as shallow marine environment according to paleontology, paleogeochemistry and paleogeography evidence Li et al. 2019;Zhang et al. 2019). The Paleocene uppermost formation is the E1m formation, which demonstrates particularly prominent marine-continent transition characteristics. The lower member of the E 1 m (E 1 mL) formation comprises mainly dark gray mudstone and sandstone depositions while the upper member (E 1 mU) consists of mainly light gray to brown depositions containing coal seams.
Typical sequence stratigraphic characteristics are easily recognizable in the Paleocene successions, especially in the E 1 lU member and the E 1 mL member. Sequence stratigraphic surfaces, such as subaerial unconformity, correlative conformity, regressive surface of marine erosion, basal surface of forced regression, maximum regressive surface and maximum flooding surface were identified in our previous work (see Zhang et al. 2015) based on seismic characteristics, such as onlap, downlap, offlap, truncation, etc. Falling stage systems tracts (FSST) were identified in both the two members. Thus, a quartile systems tract framework following Hunt and Tucker (1992) was also established, respectively (Fig. 3). Lishui Sag

Assumptions and simplifications
Sequence stratigraphy is still a developing science. Although some debates may exist, many consensuses have already been made, including that allogenic factors (eustacy, tectonics, climate, etc.) are key controls on the development of paleogeomorphology and stratigraphic sequences (Catuneanu 2006). Among the allogenic factors, eustacy and tectonics in most cases are the two most important ones and are jointly known as RSL. While the term base level is slightly different and generally includes one more factor, the climate. For simplicity of the problem, several assumptions and simplifications were made. The Paleocene Lishui Sag is of shallow marine environment located in a back-arc rift basin, and the climate factors contributed much less than RSL (Zhang et al. 2012). Thus in this paper, climate factors were assumed neglectable. Under this assumption, the base level and the RSL are almost the same that these two terms are used as synonyms in the following text except for specific discussions.
Besides allogenic controls, sediment supply (both rate and content) is another important factor determining the depositional patterns, which is primarily a function of climate (Rubin and McCulloch 1980;Leeder et al. 1998;Catuneanu 2006;Fuller et al. 2009;Liu et al. 2017a). Since the climate factor was neglected, the sediment supply was assumed as a constant in this study. Although this may be far from rigorousness, it is still a commonly used mean for both experimental and numerical studies (Muto and Steel 1997;Nelson and Morgan 2018;Zhang et al. 2018).
Within a single systems tract, we assumed that the stratigraphy was continuously developed and no differential subsidence occurred. And we assumed the shoreline trajectory was continuous between two adjacent systems tracts. While due to the incompleteness nature of stratigraphic records, the shoreline trajectory was not necessarily continuous between two systems tracts. This will be discussed in Part 5.

General workflow
The general workflow of this study is shown in Fig. 4 and it can be summarized as (1) Decompaction the target section, (2) Running of the BFM model, and (3) Comparison and modification of the BFM results.
Firstly, the time-depth conversion of the target seismic section CC' (see Fig. 3b for location) was made according to the VSP data of Well1 adjacent to this section. Then decompaction of each systems tract from the section CC' was conducted through backstripping (details see "Appendix 1" section). The decompacted stratigraphic thickness was assumed as the total accommodation space during the whole deposition process. Then, time step and mesh system of the BFM model were tried and defined according to the geological background and the information of the target section CC' (see Fig. 1). Then, the BFM model was run with predefined initial topography and sediment supply and RSL parameters (for details of the BFM model please see the following Part 3.3 and "Appendix 2" section). The predefined parameters did not always fit the expected results, so a negative feedback strategy was used to modify these parameters.

The BFM model
The BFM simulations were implemented by a program named Sedapp, which was written in R codes by the authors. This diffusion-based forward stratigraphy modelling program was developed based on many similar models, such as Syvitski et al. (1988), Flemings and Jordan (1989), Kaufman et al. (1991), Rivenaes (1992Rivenaes ( , 1997, Paola (2000) and Harris et al. (2016). While it has some differences compared with these above models, especially in the nonlinear diffusion coefficient settings. For mathematical details please see "Appendix 2" section.

Results
After iterated trials, we finally got the results below, which fit the seismic data very well (Fig. 5). The present underground framework, a reference framework after decompaction, and a simulated stratigraphic framework are summarized in Fig. 5. The Wheeler diagram and 3D base level diagram is shown in Figs. 6 and 7.  Four systems tract, from bottom to top, are called "FSST", "LST", "TST" and "HST" (Fig. 5a), corresponding to previously interpreted falling stage systems tract, lowstand systems tract, transgressive systems tract and highstand systems tract of the Paleocene E 1 lU Member (see Zhang et al. 2015). Also, the simulation time periods are called ST.1, ST.2, ST.3, and ST.4, respectively. Detailed internal architectures and characteristics of each systems tract are described in the following sections.

E 1 lU "FSST"
The simulated lowest systems tract, "FSST', covers almost the whole area except for the northwest end, and has a median thickness of about 251 m. The clinoforms are well developed and are of a relatively steep gradient. Offlap and downlap features are apparent (Fig. 8a). The shorelines points move forward into the basin while their altitudes gradually fall. Shallow marine facies is the prevailing facies, while the subaerial part is seldom developed. Overall, the "FSST" demonstrates a coarsening-upward pattern. The base

E 1 lU "LST"
The simulated "LST" covers a larger area than the underlying "FSST" and has a mean thickness of about 401 m. The clinoforms are also well developed, while their slopes are relatively gentler than that of the "FSST". Both onlap and downlap features are apparent (Fig. 8c). The shorelines points continue to move forward into the basin while their altitudes gradually rise. Shallow marine facies is still well developed, while the subaerial part is also well developed in the northwest part. Overall, the "LST" also demonstrates a coarsening-upward pattern. Base level slowly rises while it is location-sensitive that the central part has relatively high values due to differential subsidence (Fig. 8d).

E 1 lU "TST"
The simulated "TST" covers a similar area with the underlying "LST" and has a mean thickness of about 296 m. The clinoforms are of a relatively gentle gradient. Onlaps are apparent, while downlaps are not (Fig. 8e). The shorelines points turn back toward the continent and their altitudes continue to rise. Shallow marine facies gradually prevails again, while the subaerial part shrinks. Overall, the "TST" demonstrates a fining-upward pattern. The base level continues to rise and is also location-sensitive due to additional differential subsidence (Fig. 8f).

E 1 lU "HST"
The uppermost "HST" covers a smaller area but occupies the northwest end for the first time. Its maximum thickness can reach 607 m while its average is at 186 m (Fig. 8g). Unlike the first three systems tracts, the simulated "HST" is not a typical HST. The shoreline points continue to move landward (even somewhere off the chart) and their altitudes also continue to rise. Shallow marine facies is the prevailing facies, while the subaerial part only shows up at the lowerleft corner. The base level is still location-sensitive while rises at a higher speed (Fig. 8h). A new interpretation of this "HST" will be discussed in Part 5.

Discussions
In most cases, it is very difficult to uniquely reconstruct a geological unit through limited stratigraphic records. Fortunately, our model is an open system, the results can be updated when new information is added. For example, if we had some new evidence that the shoreline trajectory was not as assumed in Sect. 3.1 but some other situations, then we might get somewhat different results (see Figs. 6 and 9). More accurate calibration may depend on more comprehensive shoreline evidence, especially from new drilled cores in the future. This is also true for the base level curve shapes. Higher frequency base level details, although have little impact on the overall sequence geometry, can determine the exact shoreline locations. However, due to the lack of accurate shoreline details, straight lines are the best choices for a rough approximation. Things are a little different for the ST.2 ("LST"). In order to Sequence stratigraphy frameworks of the study interval. a The present underground framework interpreted from seismic data; b A reference framework derived from a after decompaction; c Simulated sequence stratigraphy avoid an autoretreat (see Muto and Steel 2002) within a normal regression unit, the straight line is replaced by a convex up curve (Fig. 10). Similarly, the simulation time t and the sediment transport coefficient α are of a close relationship (See Eq. 9 in "Appendix 2" section). If the other parameters were kept fixed, these two can determine the final volume of the sediments. We used the same α throughout the four systems tracts, so we get different ts (Table 1). Also, t here is a relative time measure and its unit is "step". Considering the overall time span of the study interval (about 0.8 Ma, from 60 to 59.2 Ma according to Zhang et al. 2019) and taking no account of erosion or sedimentation hiatus, each step corresponds to 297 years. If more detailed aging information was added, more accurate time-stratigraphy could also be supplemented.
We indeed cannot completely eliminate the non-uniqueness of the solutions under present conditions, but we can still check the reasonableness of some previous interpretations. Unlike the lower three systems tracts, the previously interpreted "HST" has very few typical seismic reflections for identification. Onlap/downlap features are not available, neither is the direct drilled evidence. Its modelling process makes it even distinctive. The "HST" could not be reconstructed through typical HST parameters (e.g., relative slow or medium base level rise, normal regression. If we really want to do this, the parameters can be strange). Conversely, rapid base level rise and transgression could easily make it. This also coincides   with the common sense that sediments are more likely to accumulate on the landward sides during transgressions and on seaward sides during regressions (Muto and Steel 1997). Under this context, re-interpreting the "HST" as a TST (or part of a TST) seems more plausible. Global sea level has long been a controversial topic. Different people made different results through different methods (Müller et al. 2008). Among them, Haq et al. (1987) and Miller et al. (2005) are two frequently cited high-resolution curves (Fig. 11). During the forming period of the study interval E 1 lU, the two curves from different sources in Fig. 11 have little in common except their magnitudes (both are smaller than 200 m). They are both one order of magnitude smaller than that of the base level in this study (ranging about from − 500 to 1500 m). Hence, the dominant contributor to the base level change may be the basement subsidence rather than the global sea level fluctuation, considering that it was then at the intensive Syn-rift phase (Fig. 2). Differential subsidences make the accommodation creations different from place to place. This makes the base level change local-sensitive and serves a good reason that it should be plotted in a 3D graph like Fig. 7. In addition, isostasy and sediment compaction can also contribute   Fig. 10 Base level curves at different locations (The x coordinates can be referred to Fig. 6). The base level changes between two adjacent systems tracts are assumed as instantaneous in this relative time scale to the local-sensitivity of the base level change according to Airy's flexural theories (Watts 2001). Besides, like the tectonism having many phases, the base level change can also be phased. An abrupt change of the subsidence between two adjacent phases may lead to the abrupt raise of the base level (Fig. 12).

Conclusions
Base level change in an area is not necessarily a single curve. It can be a series of curves varying from location to location, which can be attributed to differential subsidences. According the BFM results, the dominant controlling factor of the base level may be tectonism. The base level change here is indeed location-sensitive, especially in the dip direction. However, besides differences, there are also some similarities that have nothing to do with the location: (1) During the period of ST.1, the base level continuously fell that made the shoreline move seaward quickly.
(2) During the period of ST.2, the base level rose relatively slowly that made the shoreline move seaward slowly.
(3) During the formation of the periods of ST.3 and ST.4, the base level rose very quickly that made the shoreline move back landward. The BFM is not only a powerful tool for reconstructing the base level, but also an effective way to check the reasonableness of previous interpretations. The previously interpreted "HST" may not be an HST, it is more like a transgressive unit.
Although simulated solutions may be of non-uniqueness due to the limitation of current data, the BFM still provides us a chance to gain some insights into the mechanism and dynamic details of the stratigraphic formation.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Appendix 1
The backstripping techniques are based on a fundamental assumption that the skeleton mass of the strata does not change with the compaction process. So we have: (1 − (z))dz  x, km z, m -1000 Between "FSST" and "LST" Between "TST" and "HST" Between "LST" and "TST"  Specifically, the porosity φ(z) could be seen as a combination of several functions of different lithologies. Since there is only a small amount of limestone in the study interval, the φ(z) function here was regarded as a combination of only sandstone and mudstone for simplicity: p s : the volume fraction of sandstones; p m : the volume fraction of mudstones; φ s (z): the porosity function of sandstones; φ m (z): the porosity function of mudstones.
The porosity functions of sandstones and mudstones were derived from well logging data of Well1. The essence of the decompaction is to solve a series of H ′ 1 and H ′ 2 along the target section. This process was implemented through a discretization of Eq. (1): here n and m are the user-defined divisors divided by the current thickness and original thickness of the study interval. The H 1 , H 2 are known numbers. Considering the original status means just at the end of the deposition, the H ′ 1 is actually 0. The left hand side (LHS) of Eq. (3) could be calculated directly, and the only unknown is the H ′ 2 on the right hand side. Through rearrangement, Eq. (3) becomes: Since both left hand side and right hand side of Eq. (4) has the unknown H ′ 2 , an iteration method was used to calculate the result. For der = 1 and n = 2, the 3D (exactly 2DH, for the h is another dimension perpendicular to x and y) scenario of Eqs. (5) and (6) can also be expressed as: here x and y are spatial coordinates. This is especially suitable for the cases dealing only with 2 classes of lithology for simplicity, where Γ 1 is the transport coefficient for sand and Γ 2 is the transport coefficient for mud.
Note that the Γ 1 or Γ 2 cannot be put outside the parentheses, because they are not constants but functions of spatial coordinates and the time. The expression of a general Γ for the 1DH scenario can be expressed as: here α is a preexponential factor (L 2 /T), which has a close relationship with the simulated time. While β is a spatial scale factor (L −1 ), that can determine the gradient of the clinoforms. The d(x,t) is a function of spatial coordinates and time, denoting the distance from the shoreline (for the marine part).
Numerical methods were used to solve the above model. The numerical solution procedure was based on the finite volume method (FVM), which has a good local conservativeness and clear physics meaning (Versteeg and Malalasekera 2007;Moukalled et al. 2016;Liu et al. 2017b;Wang et al. 2017). Cell-centered variable arrangement method was used to store the unknowns at the centroids of grid elements (Appendix Fig. 13).