*Article* **Reconstruction of Lake-Level Changes by Sedimentary Noise Modeling (Dongying Depression, Late Eocene, East China)**

**Zhongheng Sun 1,2, Tao Jiang 1,2,\*, Hongtao Zhu 3,\*, Xinluo Feng <sup>4</sup> and Pengli Wei <sup>5</sup>**


**Abstract:** The late Eocene succession of the Dongying Depression forms a highly productive hydrocarbon source. However, due to lack of an unambiguous fine chronostratigraphic framework for the late Eocene stratigraphy, it is challenging to understand the paleolake's evolution and the driven mechanism of lake-level variation, a limitation which hinders hydrocarbon exploration. In this work, high-resolution gamma-ray logging data were analyzed to carry out the cyclostratigraphic analysis of the third member (Es3) of the Shahejie Formation in the Dongying Depression. Significant 405-kyr eccentricity cycles were recognized based on time series analysis and statistical modeling of estimated sedimentation rates. We abstracted ~57 m cycles of the GR data in the Es3 member, which were comparable with the long eccentricity cycles (~405-kyr) of the La2004 astronomical solution, yielding a 6.43 Myr long astronomical time scale (ATS) for the whole Es3 member. The calibrated astronomical age of the third/fourth member of the Shahejie Formation boundary (41.21 Ma) was adopted as an anchor point for tuning our astrochronology, which provided an absolute ATS ranging from 34.78 ± 0.42 Ma to 41.21 ± 0.42 Ma in Es3. According to the ATS, sedimentary noise modeling for the reconstruction of lake-level changes was performed through the late Eocene Es3. The lake-level changes obtained based on sedimentary noise modeling and spectrum analysis reveal significant ~1.2 Myr cycles consistent with global sea level variations which were related to astronomical forcing. Potential driven mechanisms of marine incursion and/or groundwater table modulation were linked to explain the co-variation of global sea level changes and regional lake level changes. Our results suggest global sea level fluctuations may have played an important role in driving the hydroclimate and paleolake evolution of the late Eocene Dongying Depression.

**Keywords:** cyclostratigraphy; lake level; sedimentary noise modeling; Dongying Depression

## **1. Introduction**

Cyclostratigraphy is an important branch of stratigraphy that deals with the sedimentary record of astronomically forced paleoclimate change [1–4]. The identification and deciphering of these sedimentary records improve our understanding of rates and mechanisms of biological, chemical, and physical changes, paleoclimate variations and key geological problems in Earth's distant past [4–7]. Previous cyclostratigraphy studies have mainly focused on marine environments such as deep-water facies, carbonate platforms and mixed and detrital marginal facies [8–10]. Terrestrial lake basins have recently become a hotspot for intensive petroleum exploration, with continuous successful discoveries [11–13]. Many terrestrial basins have therefore been used to conduct astronomical dating and paleoclimate reconstruction, such as Piceance Creek Basin [14], Songliao

**Citation:** Sun, Z.; Jiang, T.; Zhu, H.; Feng, X.; Wei, P. Reconstruction of Lake-Level Changes by Sedimentary Noise Modeling (Dongying Depression, Late Eocene, East China). *Energies* **2023**, *16*, 2216. https:// doi.org/10.3390/en16052216

Academic Editor: Reza Rezaee

Received: 31 December 2022 Revised: 13 February 2023 Accepted: 23 February 2023 Published: 24 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Basin [15] and Bohaibay Basin [13,16–21]. The Bohai Bay Basin (BBB), located on the eastern part of the North China Block, is one of the largest hydrocarbon-bearing basins in China, accounting for approximately one-third of the total hydrocarbon production of the country [11]. The third member of the Shahejie Formation (Es3) forms the most productive hydrocarbon source rock. Understanding the paleolake evolution and lake level fluctuation history is vital for clarifying the development mechanism of source rock and further enhancing hydrocarbon exploration [22]. Therefore, it is both scientifically and economically significant to fully understand the factors that influence the stratigraphic architecture and paleolake evolution in these basins. Reconstruction of high-resolution lake-level changes can better help us decipher these factors. Sedimentary noise modeling has been proven to be a powerful tool for reconstructing the sea- and lake-level changes in both marine and terrestrial lake basins [13,21,23,24]. Yet, lacking the accurate astrochronology constraints of the entire Es3 member (previous studies only focused on lower and/or middle of the Es3), the high-resolution lake level history has not been reconstructed [18–21].

Regarding the chronology and paleoclimate of the Paleogene Bohai Bay Basin, numerous previous studies have been employed through various methods of cyclostratigraphy and magnetostratigraphy and they have built multiple sets of chronostratigraphic frameworks [13,16–21]. The initial research postulated the boundary between the Es3 member and the fourth member of the Shahejie Formation (Es4) as 40,904 Ma [16]. Soon afterwards, the boundary age between the Es3 and Es4 was recalibrated by Liu et al. (2018) [17] according to the top age of the Chrons C18r of the GTS2012 [25], which yielded a new age of 42.47 Ma. However, updated GTS2020 determined the base age of Chron C19n as being 41.18 Ma, which resulted in the calibrated age being 42.47 Ma. This age is about 1.29 Ma older than the previously defined 41.18 Ma rising new controversy. Meanwhile, new ATS of the upper of Es4 and lower of Es3 was established based on high-resolution measurements of magnetic susceptibility of Well Fanye-1 in the same year [18]. However, lacking the accurate age anchor points was a constraint that caused these ATS to have various degrees of discrepancy. Subsequently, Shi et al. (2019) [19] updated the ATS using new acquired paleomagnetic data but lacked a statistically-based astronomical sedimentation rate estimation to better constrain the ATS. Given this, multiple statistical methodologies were employed to rebuild the ATS for Well Fanye-1 with an updated age model in GTS2020, using the short eccentricity cycles(~100-kyr) tuning [20]. Terrestrial strata are characterized by widely variable sedimentation rates and transient depositional or erosional hiatuses. The ~100 kyr cycle is thought to be unstable, which could lead to significant errors in the short-eccentricity tuning approach [21]. Synthesizing the above differences, the new ATS and magnetostratigraphy of the middle-late Eocene was further calibrated using 405 Kyr long eccentricity cycle tuning method developed by Ma et al. (2023) [21]. Compared with previous ATS and GTS2020, the model shows minor error uncertainties indicating the ATS of the Es4 and Es3 <sup>L</sup> is robust and reliable. Therefore, the ATS provided a basis for us to construct a complete astronomical time framework for the entire Es3 member.

In this paper, a cyclostratigraphy study of the late Eocene (Bartonian to Priabonian) terrestrial succession of Well S146 in the Dongying Depression is conducted using gammaray series. The overall objective of the study is to (1) construct an absolute astronomical time scale (ATS) for the late Eocene Es3 member in the Dongying Depression, Bohai Bay Basin in eastern China; (2) reconstruct lake-level changes based on sedimentary noise modeling of the tuned GR series; and (3) explore the potential driven mechanisms of million-year-scale lake level variations.

### **2. Geological Background**

The Dongying Depression is located in the southern part of the Jiyang subbasin of the Bohai Bay Basin (BBB), which is one of the most petroliferous basins in China. The Dongying Depression covers an area of approximately 5700 km<sup>2</sup> which is bounded by the Qingtuozi Uplift in the east, the Qingcheng and Binxian Uplift in the west, Chenjiazhuang Uplift in the north, and Luxi and Guangrao Uplift in the south [26] (Figure 1). It is characterized by faulting in the north and overlapping in the south, forming landscape with steep slope in the north and slow slope in the south. Additionally, several subunits can be identified, including the northern steep slope zone, the central uplift zone, the Minfeng sag, the Lijin sag, the Niuzhuang sag, the Boxing sag, and the southern gentle slope zone [26,27] (Figure 1B).

**Figure 1.** (**A**) Location map of the Bohai Bay Basin showing the major tectonic elements. (**B**) Tectonic map of the Dongying Depression showing locations of the studied Well S146.

The Depression, tectonically similar to the BBB, underwent several complex and distinct phases of Mesozoic and Cenozoic continental rifting, differential subsidence and down warping resulting in thick Paleogene synrift (65.0–24.6 Ma) and Neogene-Quaternary post-rift (24.6 Ma to the present) successions. From bottom to top, Paleogene synrift strata can be further divided into Kongdian (E*k*), Shahejie (Es) and Dongying (E*d*) formations. The Paleogene synrift experienced three stages: (1) initial rifting with lake division into sub-lake basin, (2) rift climax with lake expansion, and (3) weakened rifting with shrinkage of lake basin, forming a complete sedimentary cycle of shallow lake facies, deep lake facies, and successive fluvial and shallow lake facies [26]. The third member of the Shahejie Formation (Es3) is subdivided into the lower (Es3 L), middle (Es3 M) and upper (Es3 U) units, and the fourth member of the Shahejie Formation (Es4) is subdivided into the lower (Es4 L) and upper (Es4 U) units (Figure 2).

**Figure 2.** Generalized stratigraphic column of the Dongying Depression, modified from Feng et al. (2013) [26]. The target interval is marked in the black bold line. Fm = Formation [16–21].

The boundaries of all above members and units identified in the Shahejie Formation can be correlated and traced across the whole Depression. The Shahejie Formation was mainly deposited in lacustrine setting during the synrift period. The Es3 and Es4 members developed important Paleogene source rocks. The Es4 <sup>L</sup> is dominated by purple-red mudstone with extensively developed salt rock and gypsum. The Es4 <sup>U</sup> is dominated by gray and dark gray mudstone intercalated with carbonate, sandstone, and oil shale. The Es3 <sup>L</sup> mainly consists of dark gray mudstone and taupe oil shale and is high-quality source rock. The Es3 <sup>M</sup> comprises mudstone, argillaceous siltstone, and calcareous mudstone interbedded with thin sandstone. The Es3 <sup>U</sup> is dominated by mudstone, and contains argillaceous limestone [27,28]. Target well 146 is located in the south part of the Lijin sag, which is the sag center of the Dongying Depression. The Es4 <sup>U</sup> and Es3 <sup>M</sup> of well 146 developed massive dark mudstone, interpreted as deep and semi-deep lacustrine deposits in a stable tectonic setting. In seismic profiles, these are characterized by medium- to lowamplitude continuous parallel reflections, which are typical features of a lake depositional environment [26].

### **3. Materials and Methods**

### *3.1. Gamma Ray Logging*

In this study, cyclostratigraphic analysis was conducted on the gamma-ray (GR) logging data from the S146 well which is situated in the depocenter of the Dongying Depression where tectonic influences on sedimentation are lacking.

Gamma logging is a downhole measurement of the intensity of gamma rays emitted by the decay of naturally occurring radionuclides in sedimentary rocks. GR data record variations in the content of radioactive isotope of uranium (U), thorium (Th), and potassium (K) in the natural rocks. These three radioactive isotopes are sensitive to the variations of clay minerals and organic matter which, in abundance, are sensitive to environmental change [2]. Therefore, GR data preserve original signals related to climatic changes and can be interpreted to develop a history of the climate changes [29]. GR has proven to be one of the most efficient proxies reflecting the paleoclimatic changes in both marine and continental strata [15]. In a warm and humid environment, generally accompanied by intense chemical weathering and enhanced precipitation with higher eccentricity, more clay minerals and organic matter are input to the lake basin, corresponding to the high value of the GR logging curve. In contrast, the content of clay mineral is reduced, and the content of sandstone and the inorganic carbonate increased under cool and arid environments with lower eccentricity, which results in lower GR values. Many cyclostratigraphic studies based on the GR proxy have demonstrated the potential of GR recording the Milankovitch signals and proposed GR as a climatic proxy for identifying wet-arid climate change [15,30,31]. Meanwhile, an East Asian Lake hydrology study from four lake basins (Dongying Depression, Nanxiang Basin, Jianghan Basin, and Fushun Basin) in eastern China by Ma et al. (2023) [21] also related the hydroclimate and lake level change recorded in GR data to astronomical forcing. By reason of the foregoing, GR logging is reasonable for cyclostratigraphic analysis and lake level reconstruction in this study.

Here, we present continuous GR logging data of the Shahejie Formation in S146 well, drilling from the Lijin sag located in the middle of Dongying Depression of the BBB. The parameters for study: S146 borehole, from the base of Es3 <sup>L</sup> to the top of Es3 <sup>U</sup> with depth ranging from 3736 to 2912 m (a thickness of 824 m). The GR data of the borehole was acquired by continuous downhole measurements conducted by the Shengli Oilfield Company of SINOPEC with an average sampling rate of 0.125 m.

### *3.2. Time Series Analyses Methods*

Prior to performing spectrum analysis, a series of data preprocessing were carried out. We applied linear interpolation on the GR data to ensure that the sampling intervals were equal, because the raw sample interval of the GR data collected from the S146 borehole may be unequal owing to unknown empty sampling. The raw GR data were detrended, adopting the LOWESS method by subtracting a 30% (i.e., 247.2 m) weighted average to remove low-frequency long-term trends prior to spectral analysis [32]. To preliminarily evaluate the variations of sedimentation rates, we used a sliding-window evolutionary spectrum method (fast Fourier transform spectra) to examine the changes in orbital cycle frequencies related to varied sedimentation rates [33]. We employed the correlation coefficient (COCO) approach to objectively identify the optimal sedimentation rates in order to more properly assess the variations in sedimentation rates through the studied succession. The COCO method was developed to estimate the correlation coefficient between the power spectra of the sedimentary record in the depth domain and an astronomical solution in the time domain. Based on this, various "test" sedimentation rates were computed to convert the proxy series of sedimentary records from depth to time [34] (Li et al., 2018). The optimal sedimentation rate corresponds to the one yielding the highest correlation coefficient. Meanwhile, the evolutionary correlation coefficient (eCOCO), using a sliding-window approach were performed to track variable sedimentation rates across the Gr series and to relate power spectra of GR series to its associated astronomical solution [34]. Additionally, a null hypothesis of no astronomical force is included with significance tested using Monte Carlo simulation [34] (2000 simulations). Monte Carlo simulation is a broad class of computational algorithms that relies on repeated random sampling and statistical analysis to model the probability of various outcomes. Combined with estimated sedimentation rate, a 2π multitaper (MTM) method of spectral analysis of both depth- and time-domain GR series were carried out [35]. The robust red noise model proposed by Mann and Lees (1996) [36] was performed to statistically assess the significance of peaks obtained by MTM analyses at the mean, 90%, 95%, and 99% confidence levels. To identify which of these

spectral peaks corresponds to the orbital eccentricity, obliquity, and precession, those peaks detected above 95% confidence levels underwent additional analysis using the La2004 astronomical model. The MTM-identified long eccentricity cycles (405 kyr) and short eccentricity cycles (100 kyr) were separated using Gaussian bandpass filtering and visual interpretation. These isolated cycles were then counted and tuned to convert GR series from the depth domain to the time domain. The 405 kyr long eccentricity cycle produced by the gravitational effects of Jupiter and Venus on the Earth is the most stable astronomical orbital period since the past 250 Ma, providing a "pendulum" for geological time timing and serving as the main period for geological time calibration. Therefore, this paper uses ~405 kyr long eccentricity cycle to tune the astronomical age of the target interval, and directly compares the filters with the eccentricity curve of La2004 (Laskar et al., 2004) to establish a floating astronomical time scale (ATS).

### *3.3. Sedimentary Noise Modelling*

The sedimentary noise model was originally proposed to reconstruct past sea level change in marginal marine environments [23]. In recent years, more and more case studies have proven that the sedimentary noise model can also be used for restoring lake level changes [13,21,24]. The sedimentary noise model comprises two complementary approaches: dynamic noise after orbital tuning analysis (DYNOT) and lag-1 autocorrelation coefficient (ρ1) analysis [23]. The DYNOT model calculated the ratio of non-orbital signal variance versus the total variance via extracting noise recorded in water-depth related proxies to indicate relative sea-level changes. The ρ1 model is an independent tool for relative water-level reconstruction which is based on the lag-1 autocorrelation coefficient of both tuned and untuned dataset. Generally, high water-level (sea- or lake-level) yields a low sedimentary noise corresponding to a low DYNOT value and a high ρ1 value. In contrast, lower water-level results in high sedimentary noise level, which yields a high DYNOT value and a low ρ1value [13,21,23]. All abovementioned analyses were carried out in *Acycle v2.5* software [37].

### **4. Results**

### *4.1. Time Series Analysis*

The GR series of studied intervals of the Well S146 exhibits significant consistency with lithological variations with high GR value corresponding to mudstone whereas low GR values are dominated by fine sandstone and sandstone (Figure 3a). The MTM power spectrum analysis of the untuned GR logging data in depth domain of the Es3 interval of Well 146 in the Dongying Depression shows prominent periodic peaks, with confidence levels more than 95%, at wavelengths including 57, 28.8, 14.3, 7.0, 5.4, 4.7, 2.9, 2.7 and 2.2 m (Figure 4c).

The COCO analysis proposed by Li et al. (2018) was used to obtain the potential optimal deposition rate. A total of 2000 Monte Carlo simulations with tested sedimentation rates ranging from 0 to 25 cm/kyr were carried out to perform a null hypothesis test for non-astronomical orbital period signals. The COCO analysis of untuned GR data from the S146 borehole shows multiple peaks of sedimentation rate at 1.8, 4.0, 6.8, 9.4, 13.0, and 16.0 cm/kyr, with null hypothesis (H0) significance level no more than 0.01 (Figure 3c,e). Among them, the most remarkable cluster of peaks is from 13.0 to 16.0 cm/kyr (with the highest peak at 13.0 mc/kyr), which may indicate the potential mean sedimentation rate is between 13.0 and 16.0 cm/kyr. Notably, the sedimentation rates band ranges from 13.0 to 16.0 cm/kyr is in line with a mean sedimentation rate of 16.8 (Es3) and 10 cm/kyr (Es3 L) in the Dongying Depression as inferred by Liu et al. (2018) and Shi et al. (2018), respectively.

The eCOCO analysis of the same untuned GR series in Well S146 was performed to track variable sedimentation rates across the studied interval. We ran 2000 Monte Carlo simulations with tested sedimentation rates ranging from 0 to 25 cm/kyr and a sliding step of 3 m and sliding window of 200 m. The results of the eCOCO reveal sedimentation rates of 13.0–16.0 cm/kyr, producing a broader frequency band than those of 4.0–6.8 cm/kyr consistent with results of COCO analyses (Figure 3d). Moreover, the correlation coefficient values for sedimentation rates of 13.0–16.0 cm/kyr are greater than 0.3 and typically substantially higher than those at sedimentation rates of 4.0–6.8 cm/kyr. A sedimentation rate of 13.0–16.0 cm/kyr is relatively stable through the entire succession. To further determine which of these two groups of peaks is the optical sedimentation rate, we boldly assumed about 5.4 cm/kyr (the average value of 4.0–6.8 cm/kyr) as the mean sedimentation rate and yielded the duration of the Es3 at 15.8 Ma which deviate from our existing understanding [37]. Taken together, the estimated average sedimentation rate of 13.0 cm/kyr for Es3 member derived from COCO analysis is robust and consistent with estimations of sedimentation rate produced by the eCOCO.

**Figure 3.** COCO analysis and eCOCO sedimentation rate map of the GR series in the Well S146 (2912–3736 m): (**a**) lithological column of the Well S146; (**b**) raw GR series of the Well S146; the correlation coefficient (**c**) and evolutionary correlation coefficient (**d**) shown with sedimentation rate curve (black line) based on 405-kyr tuning of the longest (∼57 m) statistically significant cycles in the data; (**e**) null hypothesis test; and (**f**) evolutionary null hypothesis (H0) significance level.

**Figure 4.** (**a**) Detrended and raw GR series of the Well S146 in depth domain; (**b**) filters with wavelength of ~57 m and ~14.3 m; (**c**) 2π MTM power spectrum in depth domain shown with background AR(1) model, and 90%, 95%, and 99% confidence levels; and (**d**) evolutionary power spectra calculated using a 400-kyr running window.

The wavelength ratio of the clusters of 57, 14.3, 7.0-5.4-4.7, and 2.9-2.7-2.2 m (ca. 57:14:5.4:2.7) identified in the MTM power spectrum of the untuned GR series of Well S146 is near the ~20:5:2:1. Moreover, based on our objective sedimentation rate estimate by the COCO and eCOCO analysis, the estimated average sedimentation rate is 13.0 cm/kyr. This yields tentative durations of ~438 kyr for the~57 m cycle, ~110 kyr for the ~14.3 m cycle, ~54-42-36 kyr for the ~7.0-5.4-4.7 cycles, and ~22-20-17 kyr for the 2.9-2.7-2.2 m cycles. These periodicities are broadly consistent with the expected ratio of long-eccentricity (405 kyr), short-eccentricity (~100 kyr), obliquity (41 kyr) and precession

(19–23 kyr), respectively. Therefore, the sedimentary strata can be linked to the above astronomical forcing.

### *4.2. Astronomical Tuning*

According to the above analysis of MTM, evolutionary spectrum, and the COCO and eCOCO of S146 borehole, ~57 m sedimentary cycles were tied to the 405 kyr longeccentricity metronome. To extract the counterpart cycles from the GR series, we applied a bandpass filter using a ~57 m wavelength (passband of 0.0175 ± 0.004375) obtaining about 16 long-eccentricity cycles (405-kyr). Meanwhile, cycles of ~14.3 m wavelength (passband of 0.07 ± 0.0175) were band-pass filtered to isolate short eccentricity cycle (100 kyr) yielding ~64 cycles. Tuning ~57 m sedimentary cycles using the 405-kyr eccentricity tuning approach allowed us to build a floating ATS for the Es3 member with duration of ~ 6.43 Ma (Figure 5).

**Figure 5.** (**a**) Tuned GR series of the Well S146 in time domain; (**b**) filters with cycles of 405 kyr and 100 kyr; (**c**) eccentricity curve and 405/100 kyr filters of La2004 [38]; (**d**) 2π MTM power spectrum in time domain shown with background AR(1) model, and 90%, 95%, and 99% confidence levels; and (**e**) evolutionary power spectra calculated using a 400-kyr running window.

Based on this tuning method, the entire Es3 interval yields a mean sedimentation rate of ~12.8 cm/kyr which compares well with those (13.0 cm/kyr) estimated with the COCO/eCOCO analysis. The MTM power spectrum of the tuned GR series in Well S146 shows significant cycle peaks of 405, 207.4, 138.6, 100, 59.2, 49.6, 43.6, 34.3, 28.1 and 23 kyr. The evolutionary spectrum of the tuned GR series reveals obvious 405 kyr periodic signals and 100 kyr periodic signals close to a straight line. The calculated sedimentation rate with 405 kyr tuning varies from 11.1 to 15.5 cm/kyr, with an average sedimentation rate of 12.8 cm/kyr. The age of boundary between the Es4 and Es3 (BES43) is used as anchor point for calibrating our floating ATS. In addition, the Es3 <sup>L</sup> and Es3 <sup>M</sup> boundary (BES3 LM) is adopted as checked point to further verify our ATS. Recent astrochronologic and magnetostratigraphic study conducted in the Dongying Depression has provided a robust astronomical time scale and magnetostratigraphic timescale for the Es4 and Es3 <sup>L</sup> formations, with the BES43 age of 41.21 Ma and the BES3 LM age of 38.72 Ma, respectively [21]. According to the anchor point, our floating ATS of were able to transform to absolute ATS extending from 34.78 ± 0.42 Ma back to 41.21 ± 0.42 Ma from top of the Es3 <sup>U</sup> to base of Es3 <sup>L</sup> (Figure 5a). This period corresponds to the La2004 astronomical solution from 86th to 103rd 405-kyr eccentricity cycles (E86–E103) [38] (Figure 5).

### *4.3. Sedimentary Noise Modeling of Lake-Level Changes*

We reconstructed sedimentary noise curves for the Dongying Depression based on DYNOT and ρ1 modeling, with 5000 Monte Carlo simulations of the tuned GR data of Es3 using a 400 kyr sliding window. Sedimentary noise curves of DYNOT and ρ1 modeling display consistent patterns. Significant enhanced sedimentary noises were identified in five intervals and highlighted in gray bars (Figure 6).

**Figure 6.** Sedimentary noise model interpretation of lake-level changes in Well S146. (**A**) Tuned GR series of the Well S146 (black) with 405-kyr filtered output (red); (**B**,**C**) DYNOT and ρ1 models of tuned GR series; (**D**) global sea level changes modified from Miller er al., 2020 with 1.2 Myr filtered output [10]; (**E**) 1.2 Myr cycles filtered from DYNOT and ρ1 models; and (**F**) earth's obliquity solution (Black, Laskar et al., 2004) [38] and its 1.2 Myr-AM cycles (red).

According to spectral analysis of the DYNOT output series (median value), we observed five cycles in sedimentary noise series with periods of 1.2 Myr (Figure 7). Similarly, spectral detection of the ρ1 output series (median value) reveals cycles that have periods of 1.2 Myr (Figure 7).

**Figure 7.** Periodograms of the detrended median output by DYNOT (**a**) and ρ1 (**b**) models, AM of the sea level curve by Miller 2020 [10] (**c**) and AM of the obliquity (**d**) in La2004 [38] astronomical model from 34.78 to 41.21 Ma.

### **5. Discussion**

### *5.1. An Astronomical Time Scale (ATS) in the Dongying Depression*

We constructed a ~6.43 Myr long floating ATS for Es3 member in Well Shi 146 of the Dongying Depression through the late Eocene. According to our cyclostratigraphic analysis, the studied interval has a total duration of 6.43 ± 0.42 Ma, ranging from 41.21 ± 0.42 Ma to 34.78 ± 0.42 Ma. The base age of Es3 L, Es3 <sup>M</sup> and Es3 <sup>U</sup> are 41.21 ± 0.42 Ma, 38.72 ± 0.42 Ma, and 36.28 ± 0.42 Ma, respectively.

The 405 kyr eccentricity bandpass filter extracted from our constructed ATS has a significant in-phase correlation with the 405 kyr filter extracted from the La2004 eccentricity solution. The prominent eccentricity, obliquity, and precession cycles were identified based on the MTM spectra analysis of the tuned GR series. These cycles match well with astronomical models of La2004 [38]. The sedimentation curve generated from 405 kyr cycles compares well with the objective sedimentation rate estimated by COCO and eCOCO (Figure 3).

The uncertainty of the tuning and anchor point may result in additional uncertainty. Therefore, it is necessary to conduct uncertainty analysis of all potential sources. Three sources of uncertainty for our ATS and reported ages were considered including: (1) the uncertainty of BES43 age in the Dongying Depression, as reported by Ma et al. (2023) [21]; (2) the uncertainty of contact relationship of Es4 and Es3; (3) the uncertainty of the floating ATS we established. The ages of base boundaries of Chrons C18n.1n through C20n with mean errors in the FY1 borehole were reported by Ma et al. (2023) [21] and compared to those of previous six studies. The base of the Es3 roughly corresponds to the base of the C19n. Although Ma et al. (2023) [21] did not report the uncertainty of BES43 age

in the Dongying Depression, he indeed provided seven uncertainty errors spanning the boundary between the Es4 and Es3 ranging from C20n to C18n.1n. These uncertainties ranged from 0.100 Ma to 0.106 Ma, with a mean uncertainty of 0.103 Ma. Thus, we suggest 0.103 Ma as the uncertainty of the BES43 age in the Dongying Depression. Seismic reflection profiles showing the boundary between the Es4 and Es3 are characterized by stable seismic conditions even with parallel reflection and conformable across the Dongying Depression [26]. There is no report of onlap or traction seismic reflection which can indicate significant tectonic activities. Our S146 well is located in the depocenter of the Dongying Depression in a stable sedimentation setting with volcanic activities reported. Taken together, these points of evidence eliminate the second source of uncertainty. An empirical uncertainty of 0.103 Myr for the BES43 in the Dongying Depression may be obtained from superposition calculation of the first two sources of uncertainty. The identified long eccentricity cycles in our floating ATS may cause an additional uncertainty of one 405-kyr cycle for the third source. Consequently, we can obtain a total uncertainty of 0.42 Myr (calculated form <sup>√</sup> 0.1032 + 0.4052) [13] for the base of the Es3 L, Es3 <sup>M</sup> and Es3 U.

Our ATS of Es3 member appears to be compatible with previous astrochronology of the same interval performed by cyclostratigraphy studies [13,17,18,20,21]. The base ages of Es3 <sup>M</sup> reported by Shi et al. (2018) [18], Jin et al. (2022) [20] and Ma et al. (2023) [21] are 39.23 Ma, 39.556 Ma and 38.72 Ma respectively, which have discrepancies of 0–0.836 Ma. The top age of Es3 <sup>U</sup> reported by Liu et al. (2018) [17] and Wang et al. (2020) [13] leads to 1.21–1.3 Ma discrepancies. The duration of the Es3 of our ATS is 6.43 Ma, which is quite close to the 6.48 reported by Liu et al. (2018) [17]. Therefore, we argue these discrepancies could be caused by diachronous sub-member boundaries within the Es3, age differences between anchor points or erroneous identification of long eccentricity cycles in either or both those previous studies and our study [13]. Overall, we argue that we have constructed a robust and reliable absolute ATS of the Es3 of the Dongying Depression according to our uncertainty constraints discussed above.

### *5.2. Verification of Sedimentary Noise Model for Reconstructing Lake Levels*

Sedimentary noise in data series can be influenced by multiple factors including proxy related noise (e.g., dating uncertainty), unsteady sedimentation such as post-depositional reworking, local, short term tectonic, depositional, and volcanic activities in additition to lake level changes [23]. The noise associated with paleoclimate proxy indicators is affected by three factors: index sensitivity, measurement error, and dating error. Both previous studies and this study have confirmed that GR signals can well record paleoclimate changes in this area, indicating that the sensitivity of the indicator persists [20,21]. Obtaining logging data by wireline logging is a very mature technology. At the same time, the measurements are taken in a stable borehole, so the measurement errors are relatively small. Dating in this study is based on the COCO and ECOCO sedimentation rate assessment models, which are very close to the previous age anchor point, indicating our ATS is very reliable. Unstable depositional environment and short-term tectonic activities may lead to increased depositional noise. However, the sedimentary environment of the Es3 in Well Shi 146 was a semi-deep lake-deep lake facies with stable structural environment during the depositional period, and no faults were developed [39]. Therefore, unstable depositional environments and short-term tectonic activity had limited impact on sedimentary noise. A hydrological shift driven by volcanic activity may have increased additional sedimentary noise. However, lithostratigraphic observations showing no volcanic rocks were developed in Es3 of Well Shi 146, implying volcanic activities did not affect the Lijin sag and related sedimentary noise [39,40]. Post-depositional reworking mainly includes bioturbation and diagenesis, which may have a greater impact on sedimentary noise. Nevertheless, there is no report of bioturbation structure in this well. Moreover, the lithology of the Es3 is dominated by mudstone, which suffered normal compaction, being less affected by postdepositional reworking. And yet, the complexity of sedimentary noise of terrestrial strata has not yet been thoroughly understood. The similar variation patterns of the DYNOT

and ρ1 models combined with the above discussion indicates that the sedimentary noisemodeling approach is nonetheless an effective method for lake-level reconstruction in lacustrine strata.

### *5.3. Astronomical Forcing on Lake-Level Fluctuations in the Dongying Depression*

Astronomical forcing of lake-level fluctuations has been reported throughout the Phanerozoic at several basins with different latitudes across the word [13,21,24,41,42]. The robust ATS constructed in this study provides a basis for linking the astronomical cycles, lake hydraulic cycle and global sea level changes. The DYNOT modeling results shows significant anti-phase relationship with that of the ρ1 modeling. Multiple alternation of high and low sedimentary noise level was observed throughout the entire Es3 succession. The high sedimentary noise levels were interpreted as low lake levels. In opposition, the low sedimentary noise levels were interpreted to correspond to high lake levels.

Periodogram analysis of both the DYNOT and ρ1 modeling output median values showing clear long-term cycles with periods of ~1.2 Myr (Figure 7). The 1.2 Myr cycles were identified in both the DYNOT and ρ1 modeling output series, which show similar anti-phase relationships compared to their original modeling results. This may correspond to 1.2 Myr obliquity modulation cycles, indicating a potential causal link between astronomical forcing and million-year scale lake-level change. Previous studies have reported ~1.2 Myr periodicities in the Eocene [13,43–45]. The obliquity modulation cycle of ~1.2 Myr is demonstrated to affect heat and moisture transport flux from low latitudes to high latitudes and further manage the meridional insolation gradient on Earth [23,46]. During the 1.2 Myr obliquity nodes (modulation minima), heat and moisture transport flux to high latitudes would have been reduced. Consequently, lakes in mid-to-high latitudes would become cooling and drying with low water level and vice versa [13,21,23]. Our results show a clear in-phase relationship of 1.2 Myr cycles isolated form DYNOT and ρ1 with 1.2 Myr obliquity amplitude modulation cycles of La2004 solution during the late Eocene. The global sea-level curves proposed by Miller et al. (2020) [10] were used to compare with the lake level changes of the Dongying Depression. Spectral analysis of the sea-level records shows clear ~1.2 Myr cycles which is close to 1.2 Myr obliquity cycle. The 1.2 Myr cycles of lake level in Dongying Depression are generally in-phase with those of sea-level changes and obliquity cycle. These results suggest 1.2 Myr obliquity drove the changes of both lake-level and sea-level during the late Eocene. Identifying the driven mechanism of the coevolution of the global sea level and reginal lake level remains a challenge. Multitudinous hypothesizes have been made to explain this co-variation [13,21]. We believe that one at least or more of the following mechanisms may have contributed to the regulation of the coevolution in the late Eocene.

Regional transgression events during marine transgressions impose significant effects on the rise of the lake level. Four intervals of marine incursions into the Dongying Depression have been reported based a combination of interpretations with high B/Ga ratios, high S/TOC ratios and high Sr/Ba [47]. Moreover, marine fossils including foraminifera, calcareous nanofossils and trace fossil *Paleodictyon* were observed and further confirmed the marine incursion in Dongying Depression [21,47]. Our reconstructed lake-level changes appear to be well consistent with global sea level changes implying that marine incursions play an important role in regulating the lake level changes in the Dongying Depression.

Groundwater table modulation is another driven mechanism for explaining the covariation. By an integrated analysis of four lake basins (Dongying Depression, Nanxiang Basin, Jianghan Basin, and Fushun Basin) in east Asia, groundwater table modulation was used to clarify the mechanism for co-variation of the lake level and sea level during the middle-late Eocene [21]. Notably, the mechanism for driving the lake level changes of these four basins spans time periods of 48.51-38.52 not fully covering our studied interval with time periods of 41.21-34.78. Despite this, the groundwater table modulation is nonetheless a powerful potential driven mechanism. The shift in lake level reflects changes in groundwater storage and serves as a gauge for the continental groundwater

table [21,48,49]. Deep faults can serve as efficient conduits for groundwater migration connecting lakes and open sea and further regulate the lake's level [50]. Subduction of the Pacific Plate in the late Eocene has an important impact on the evolution of the Bohai Bay Basin, forming a series of deep faults [31,51]. Moreover, some of the preexisting faults are reactivated and deepen into the crust. The Tanlu fault zone is the typical deep fault in which the kinematic direction shifts from left lateral motion to right lateral motion in this instance. Taken together, these make it possible that the BBB was connected to the open sea in the late Eocene.

### **6. Conclusions**

The natural GR series of the S146 borehole through Es3 member in the Dongying Depression provides an archive of paleoclimate and environmental changes for lake level reconstruction. The following conclusions can be made:

(1) We constructed a robust astronomical time scale (ATS) of the Es3 member of the late Eocene in the Dongying Depression based on the identification of 405 kyr-long eccentricity cycles in the S146 well. The developed high-resolution and continuous ATS spans the time from 34.78 to 41.21 Ma.

(2) A quantitative assessment technique of eCOCO combined with COCO analysis generates an optimal sedimentation rate of 13.0 cm/kyr and variable sedimentation rates through the studied interval, which better supports our robust ATS.

(3) Paleolake level changes of the late Eocene were reconstructed based on sedimentary noise modeling. The lake-level changes in 1.2 Myr cycles in Dongying Depression were in phase with the same cycles of synchronous global sea-level changes. We suggest that marine incursion and/or groundwater table modulation are the main driven mechanisms of the co-variation of global sea level changes and regional lake level changes.

**Author Contributions:** Conceptualization, Z.S. and T.J.; methodology, X.F.; software, P.W.; validation, X.F. and P.W.; formal analysis, X.F.; investigation, Z.S.; resources, T.J.; data curation, H.Z.; writing—original draft preparation, Z.S.; writing—review and editing, H.Z. and T.J.; visualization, Z.S.; supervision, H.Z. and T.J.; project administration, H.Z. and T.J.; funding acquisition, Z.S. and H.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research presented in this paper was funded by China Postdoctoral Science Foundation (No. 2022M712950), the Open Funds for Hubei Key Laboratory of Marine Geological Resources, China University of Geosciences (No. MGR202213). This research was also funded by the National Natural Science Foundation of China (No. 42172127, No. 41872149).

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The Geoscience Institute of the Shengli Oilfield, SINOPEC is thanked for providing access to their geological data and permitting publication of the results. The editor and anonymous reviewers are thanked for helpful comments and constructive suggestions that greatly improved our manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
