Hierarchical Milankovitch and sub-Milankovitch cycles in the environmental magnetism of the lower Shahezi Formation, Lower Cretaceous, Songliao Basin, northeastern China

SK-2 borehole in Songliao Basin provides unprecedented geological materials for investigating the Early Cretaceous continental paleoenvironment and paleoclimate in northeastern China. The lacustrine successions of the lower Shahezi (K1sh) Formation at the depth from 4,542 to 5,695 m was systematically studied using environmental magnetism and cyclostratigraphy in this study. Magnetic analysis reveals an inverse correlation between magnetic susceptibility (MS) and lithological ranks in fine clastic sediments, with the highest values in mudstones and the lowest in sandstones. The main magnetism carriers in the lower K1sh are pseudo-single-domain (PSD) and/or multi-domain (MD) magnetite with minor presence of hematite. MS was used to further explore the genesis of the environmental and climatic variations through cyclostratigraphic analysis. Sedimentary cycles of 113 m, 34 m, 13 m and 6 m can be identified in the power spectrum, which were interpreted as long and short eccentricity, obliquity, and precession cycles, demonstrating the impact of astronomical cyclicity on sedimentary rhythmicity. Floating astronomical time scale (FATS) of 4,090 kyr and 4,148 kyr were established by tuning the inferred long and short eccentricity cycles to the artificial 405-kyr and 105-kyr orbital eccentricity curves respectively. The estimated sediment accumulation rate around 28 cm/kyr confirms the rapid deposition process within the faulted lacustrine basin. Based on this study, the lake level oscillations in Songliao Basin are assumed to be shaped by long and short eccentricity, precession and semi-precession cycles during Early Cretaceous. This study also indicates that the sand-mudstone alternations deposition in K1sh is most likely driven by the seasonal discrepancies of summer insolation during semi-precession periods.


Introduction
The Early Cretaceous is a complex period characterized by a series of widespread oceanic anoxic events (OAEs), cold intervals interrupting the warm greenhouse climates, and several marine biotic crises, which are largely influenced by the global paleoceanographic and paleoclimatic changes (Pucéat et al., 2003;McAnena et al., 2013;Yang et al., 2013;Herrle et al., 2015;Matsumoto et al., 2020;Xu et al., 2022). As an essential component of the Earth's climate system, continental deposits play critical roles in documenting global climate changes and interactions with ocean systems. Songliao Basin (SB) in northeastern China is well-known as a typical continental lacustrine basin deposit during the Cretaceous. This region has attracted tremendous attention in multidisciplinary investigations related to the Cretaceous Earth evolution. However, most paleoenvironmental and paleoclimatic studies conducted in SB are generally focused on the Late Cretaceous or the Cenozoic Eras in the past decade (Wu et al., 2013a;Chamberlain et al., 2013;Wang et al., 2013;Wu et al., 2014;Gao et al., 2015;Gao et al., 2021a;Gao et al., 2021b;Wu et al., 2022), and studies focused on the Early Cretaceous remain insufficient. Fortunately, the successful completion of Songke 2 well (SK-2) provides an unprecedented opportunity to systematically recover the global evolution from the Early to Mid-Cretaceous from the terrestrial perspective (Feng et al., 2013;Gao et al., 2019) with coring from Permian basement to the Yingcheng Formation (Zhu et al., 2018). Marked by rapidly deposited lacustrine strata and thick lacustrine deposits and interbedded with surrounding volcanic rocks, Shahezi Formation (K 1 sh) is considered to be a valuable sequence that systematically documented the terrestrial climate evolution from Early to Mid-Cretaceous. In contrast to the slow depositional processes of marine sedimentation or turbulent fluvial/alluvial depositional environments, lacustrine sediments contain stable, continuous and high-resolution archives of the continental paleoenvironment (Liu et al., 2012). Compared with the complex clay mineralogy and (radio-)isotopic stratigraphy testing procedures, rock magnetism and the derived environmental magnetic interpretations are favored by paleoenvironmentalists and paleoclimatologists for their simple, rapid, cost-effective and informative measurements (King and Channell, 1991;Dekkers, 1997;Liu et al., 2012;Zhang et al., 2020). Most importantly, environmental evolution can be tracked through time from the properties of magnetic minerals, concentrations and domain states of the magnetic substance contained in rocks or sediments (Thompson and Oldfield, 1986;Verosub and Roberts, 1995).
Analyzing the magnetic parameters using cyclostratigraphy is one of the most effective methods to reconstruct the paleoenvironment and explore the astronomical origin which potentially controls the paleoclimate changes. In the study of Cretaceous strata in SB, Li et al. (2013) found that the main magnetization-carrying minerals in the Qingshankou Formation of SK-1s Well is clastic magnetite, and the lake level was affected by both tectonics and Milankovitch forcing. From the same well, Wu et al. (2013a) established an astronomical age-depth model for the Late Cretaceous lacustrine sequences using MS data, which provides a precise chronological basis for important geological events. These results indicate that magnetic parameters can be used as a reliable paleoclimatic proxy of continental successions for further interdisciplinary research. However, the paucity of related work in the Early Cretaceous has created a knowledge gap in our understanding of the continental paleoenvironment, preventing comprehensive coverage of the coupling between land and ocean climatic processes during the Early Cretaceous.
In this study, detailed rock-magnetic measurements were conducted on the samples from the lower part of the K 1 sh, SK-2, and then lithology and cyclostratigraphy were applied to preliminarily analyze the distinct sedimentary cycles. The main scientific objectives of this study are to 1) describe the magnetic characteristics of the sediments in the K 1 sh, and compare these characteristics among different lithologies, 2) provide preliminary temporal constraints on the studied strata, which can be used as a reference for future correlation studies between marine and terrestrial records, and 3) probe the response of sedimentation to the orbital and sub-orbital forcing mechanisms recorded in the K 1 sh.
2 Geological setting and SK-2 well 2.1 Geological setting With a N-NE trending axis, Songliao Basin (SB) is located on the Jiameng Block across three provinces (Heilongjiang, Jilin and Liaoning) of northeast China at the convergence of the Siberian Plate, North China Plate and the western margin of the paleo-Pacific tectonic belt (Zhou et al., 2009) (Figures 1A-C). It is one of the longest-evolved and largest Cretaceous continental basins in the world (Hou et al., 2018). SB is 760 km long, 330-370 km wide and covers an area of 260,000 km 2 . Tectonically, the formation of the basin in the Cretaceous period was closely related to the subduction of the paleo-Pacific plate beneath the paleo-Asian Plate and active magmatism due to both plume activity and subduction (Okada, 2000;Zhou et al., 2009). Following the mantle upwelling stage during the Middle and Late Jurassic, the Cretaceous SB underwent three episodes: syn-rifting stage (~150-110 Ma), post-rift thermal sagging stage (~110-79.1 Ma), and structural inversion stage (~79.1-64 Ma) (Wang et al., 2007;2016;Zhou et al., 2009;Feng et al., 2010).
The Cretaceous period in the SB had a humid to sub-humid subtropical paleoclimate with four cooling events, three warming events, and two semi-arid events . Based on the literature references referring to geology (Wang et al., 1995; Frontiers in Earth Science frontiersin.org 2011; Zhang et al., 2018), paleontology (Gao et al., 1992;Gao et al., 1999b) and oxygen isotopic materials (Chamberlain et al., 2013), these climatic events are in agreement with the climatic pattern documented by synchronous marine formations. The consistency of the cold Early Cretaceous climate patterns reconstructed from both continents, especially western Liaoning of northeastern China (Amiot et al., 2011;Wang et al., 2017), and marine enviroment (Pucéat et al., 2003;Gutiérrez-Puente et al., 2021), may indicate that the atmosphere-ocean system had a global impact on continental climate and environment.   , the black line represents interval in this study (E) The structural section across the central SB (A-A′ in (C)) and sites of the boreholes which had been drilled in Continental Deep Scientific Drilling of SB up to 2018. The yellow bars are intervals with continuous coring, and white bars are uncored intervals (modified from Wang et al., 2013).
Frontiers in Earth Science frontiersin.org 03 2.2 The SK-2 well and the K 1 sh Cretaceous International Continental Scientific Drilling Project of SK-2 started at April 2014, and was officially completed on 26 May 2018, with a total drilling depth of 7,018 m. It obtained 4,134.8 m length of sediments with a recovery ratio of 96.61% (Hou et al., 2018;Zhu et al., 2018). SK-2 well is located at the Xujiaweizi faulted depression in Songzhan area, where the lacustrine mudstones dominate the section (Wang et al., 2017). The project collected continuous sediments from the bottom of the K 1 d down to the basement of SB, with K 1 yc, K 1 sh, and K 1 h as the main target strata, which were deposited during the Early Cretaceous syn-rifting stage. Unlike the volcaniclastic sequences in K 1 yc and K 1 h, the intercalated K 1 sh is characterized by lacustrine clastic deposits containing abundant hydrocarbon rocks and coal seams (Wang et al., 2015;Fu et al., 2019).
The Shahezi (K 1 sh) Formation (5,695.1-3336 m) of SK-2 is generally composed of dark rocks interbedded with coals, and is interpreted as river, lacustrine, and delta deposits formed in reducing environment during a greenhouse condition (Wang et al., 2017). K 1 sh unconformably overlies the K 1 h, and it can be divided into two members. The Member 1 (5,695.1-4,542 m) consists of variegated conglomerates intercalated with gray coarse sandstones, grayish-black mudstones of unequal thickness at the bottom and gray silty sandstones, dark gray mudstones in the upper part, suggesting that sedimentary facies shifted from fan delta to deep lake in ascending order. The Member 2 (4,542-3,336 m) is composed of a large amount of fine-grained sediments, such as grayish-black silty sandstone and mudstone with occasional variegated conglomerate-bearing sandstone that were deposited in lacustrine environment. Li. (2018) divided the K 1 sh of the SK-2 into six sedimentary cycles according to detailed core description, with each of them displaying fining-upward trend that reveals multiperiodic climate cycles from relatively arid to humid. It is a succession with normal grading in the grain size deposition, indicating a weakening hydrodynamic condition and gradually wetter paleoclimate.

Geochronological constraints of the K 1 sh
Previous investigations have extensively estimated the chronostratigraphy for K 1 sh strata distributed in northeastern China basins by using diverse dating techniques. In the view of biostratigraphically, K 1 sh has been commonly considered belonging to the lower Cretaceous (Gao et al., 1999a). Wang et al. (2002) obtained ages from volcanic rocks of 135.96 ± 0.28 Ma at the bottom of the K 1 yc and 147.77 ± 0.28 Ma at the top of the K 1 h using 40 Ar/ 39 Ar dating method, providing a~10 Myr duration for K 1 sh.  constructed a 11.14 Myr floating astronomical time scale for the K 1 sh of the Songshen-4 well located near SK-2.
Afterwards, further precise zircon U-Pb ages of 111-115 Ma and 105-112 Ma in K 1 yc performed by Shu et al. (2007) and Huang et al. (2011) respectively, demonstrated the base of K 1 sh is supposedly formed during Aptian period. Recently, more new zircon U-Pb data derived from the SK-2 borehole restrict K 1 sh into~118-111 Ma Yu et al., 2020;Wang et al., 2022) ( Figure 1D) and correlate this non-marine succession to late Aptain-early Albian stage. This~118-111 Ma age scale is currently the mostly acceptable chronological constraint for the K 1 sh in present research.

Sampling and data
A total of 1,097 samples were collected from the lower K 1 sh (5,690-4,542 m) of the SK-2 except for the 4,900-4,918 m interval with an average interval of 1 m. Each sample was broken into pieces and placed into a 8 cm 3 non-magnetic, cubic plastic box for detailed magnetic measurements. Low-field mass-specific magnetic susceptibility (χ) is the magnetization induced by the applied magnetic field versus the value of the applied field per unit mass. The concentration of all magnetic minerals (i.e., ferromagnetic, paramagnetic, and diamagnetic minerals) affects the χ in sediments making it challenging to interpret the climate or environment fluctuations. Anhysteretic remanent magnetic susceptibility (χ ARM ) is obtained by applying a direct current (DC) magnetic field in the presence of an alternating magnetic field that is decreased from some peak value to 0. It is a measure of the concentration of only the fine low coercivity ferrimagnetic minerals in a rock sample which usually has a detrital origin in lacustrine sediments. Saturation isothermal remanent magnetization (SIRM) measures the concentration of all ferrimagnetic minerals in a sample by using a stable DC field at room temperature to apply the maximum remanence magnetization from the specimen. Numerous studies have shown that the magnetic parameters can serve as a sensitive indicator of the climate cycles (Ellwood et al., 2000;Latta et al., 2006;Boulila et al., 2010;Kodama et al., 2010;Wu et al., 2012). We used magnetic data on these three parameters to identify the most suitable climate proxy and better understand the paleoclimate evolution.

Magnetic testing
χ was measured in a low magnetic field (approximately 0.256 mT) using the AGICO MFK1-FA Kappabridge with a low frequency of 974 Hz. Anhysteretic remanent magnetization (ARM) was obtained by the ASC D-2000 AF demagnetizer with a 100 mT alternating field superimposed with a 0.05 mT DC field and its value would be divided by 0.05 mT to get χ ARM . Specimens were remagnetized in a positive 1 T field with the ASC IM-10-30 impulse magnetizer to approach the saturation isothermal remanent magnetization (SIRM). All remanences were measured with a JR-6A spinner magnetometer. Additionally, it is essential to weigh rock samples before magnetic measurements to ensure mass normalization can be conducted. Thermomagnetic κ-T curves for 10 representative samples were obtained using AGICO KLY-4S Kappabridge in argon atmosphere with a CS-3 heating component between room temperature and 700°C. These experiments were performed in the Paleomagnetism and Environmental Magnetism Laboratory at China University of Geosciences, Beijing.

Frontiers in Earth Science
frontiersin.org 04 IRM acquisitions, hysteresis loops and first order reversal curves (FORCs) were performed to the 10 samples using MicroMag Model 3,900 vibrating sample magnetometer (VSM, Princeton MicroMag 3,900) in the Institute of Geophysics, China Earthquake Administration. The maximum applied fields for IRM acquisitions are 1 T. Hysteresis loops were imparted to maximum fields of 200 mT for nine samples, 500 mT for one sample with 500 m averaging time. FORCs with 160 segments were diagramed using the FORCinel software (Harrison and Feinberg, 2008).

Cyclostratigraphy
Cyclostratigraphic analysis concerning the environmental magnetic records provides an effective method to investigate the mechanism of the climate variation imprinted in the basinal strata (Wu et al., 2013a;Wu et al., 2013b). Before time series analysis, the data series were interpolated onto an evenly sampled grid of~1.05 m and detrended to remove linear trend. Multi-taper method (MTM) was applied to identify the presence of spectral peaks in the data (Thomson, 1982). Classical red noise estimation was conducted to interpret the significance of spectral peaks (Mudelsee, 2014). Then FIGURE 2 (A) Probability density distribution of mass-specific magnetic susceptibility (χ) values (left) and boxplot (right) from the SK-2 Well for the ten lithological ranks. Mudstone (Rank = 1), silty mudstone (Rank = 2), muddy siltstone (Rank = 3) and siltstone (Rank = 4), fine sandstone (Rank = 5), coarse sandstone (Rank = 6), gravel-bearing fine sandstone (Rank = 7) and gravel-bearing coarse sandstone (Rank = 8), sandy conglomerate (Rank = 9), and conglomerate (Rank = 10) (B) Relationship between χ values and lithological ranks of the lower K 1 sh Formation in the SK-2 Well. Rank values are the same as (A).

Frontiers in Earth Science
frontiersin.org cycle ratio method was used to estimate the possible existence of the Milankovitch cycles (Weedon, 2003). The evolutive harmonic analysis (EHA) was used to evaluate the persistence and/or transience of orbital frequencies and to track the possible changes in sediment accumulation rates in depth domain (Meyers, 2014). The Gaussian bandpass filters were designed to extract the signals of the interpreted Milankovitch and sub-Milankovitch cycles (Bloch et al., 1969;Meyers, 2014). Stratigraphic series were transformed from depth to age scale using specified anchor points with tune function. Hilbert transform was performed to examine the amplitude modulation (AM) of the studied signals. These spectral analyses and astronomical calibration processes were implemented using the Astrochron package in R (Meyers, 2014). The correlation coefficient (COCO) analysis was carried out to quantitatively identify the astronomical forcing on stratigraphic cycles and accurately assess the sediment accumulation rate (SAR) with ACycle v2.1 . Evolutionary correlation coefficient (eCOCO) analysis was used to track the changes in SAR along the studied stratigraphic series .

Lithological statistics
The ten lithological ranks, from mudstone to conglomerate, are numbered from 1 to 10 in order with increasing grain size (see Figure 2A for the classification). These lithological ranks were used to fit a lake level curve, with lower ranks denoting high lake level and higher ones denoting low lake level ( Figure 3A). Histogram. m, fitdist. m, and boxplot. m functions to individually and collectively describe the χ distributions associated to each rank were applied in MATLAB.

FIGURE 3
Magnetic parameter profiles for sediments of the lower K 1 sh (4,542-5,694 m) from the SK-2 borehole (A) Lake level fluctuations (the blue shading) are modeled from the variations of lithological ranks (B) SIRM (C) χ ARM (D) χ (E) Stratigraphy, depositional facies and lithology of the studied strata, with red arrowheads pointing to the corresponding depths of the selected representative samples for rock magnetic experiments. FD: Fan-delta. L: Lake. Depositional facies are modified from . Dashed blue lines divide the succession into Subset 1, Subset 2, and Subset 3.

Frontiers in Earth Science
frontiersin.org 4 Results

Magnetic series of the lower K 1 sh formation
According to the trends of magnetic data combined with lithological ranks, the studied sequence can be divided into three subsets, including Subset 1 extending from 5,695 m to 5,400 m, Subset 2 spanning from 5,400 m to 4,900 m and Subset 3 corresponding to the interval from 4,900 m to 4,542 m ( Figure 3).
Subset 1 (5,695 m-5,400 m): The magnetic parameters had the highest values for the entire section. The average values of the three concentration-dependent parameters are 8.24×10 -8 m 3 /kg for χ ( Figure 3D), 17.64×10 -8 m 3 /kg for χ ARM (Figure 3C), and 13.94×10 -5 Am 2 /kg for SIRM ( Figure 3B) respectively. The values for the upper conglomerates are larger than the bottom mudstones, indicating the magnetic contents in the conglomerate are significantly higher than those in other clastic sedimentary rocks.

Rock magnetic characteristics 4.2.1 κ-T thermomagnetic curves
The critical temperatures at the peaks and/or inflections exhibited in the κ-T thermomagnetic curves can be used to determine the categories of magnetic minerals contained in the samples (Thompson and Oldfield, 1986;Bloemendal et al., 1988). The magnetic susceptibility (κ) of all samples increases considerably in the cooling curves from 600°C to 400°C and is higher apparently than the peak value in the heating curves ( Figures 4A, F, K), implying that the iron-bearing clay minerals (and/or silicate) were converted into magnetite (titanomagnetite) in the high temperature ranging from 600°C to 700°C. The increase of the magnetic susceptibility (κ) of the sample S-150-10 during 200°C-300°C in the heating process can be attributed to the formation of maghemite from some iron hydroxide ( Figure 4A). The value of κ descends significantly at a temperature range of 300°C-400°C, which is possibly related to the transformation from maghemite to hematite, and is consistent with the hematite signal above 600°C. The magnetic susceptibility of the samples like S-159-30 decreases slowly with the increase of temperature before 300°C ( Figure 4F), which is indicative of the presence of paramagnetic minerals in the samples (Li and Zhang, 2005). The magnetic susceptibility gradually increases from 300°C, reaches the peak value at about 500°C, then drops to the base value at around 580°C. These variations show the Curie temperature of magnetite at 580°C, which may result from the decomposition of iron-bearing paramagnetic minerals transformation into magnetite (Zhu et al., 2001) or the original existence of magnetite. The slight attenuation between 600°C and 700°C suggests the presence of hematite. The magnetic susceptibility values of samples like S-168-24 are one order of magnitude larger compared with other samples and show clear Curie temperature of magnetite at 560°C-600°C ( Figure 4K), which indicates the high concentration of magnetite in these samples.

IRM acquisition and reverse demagnetization curves
The IRM acquisition curves and reverse demagnetization curves are usually utilized to estimate the remanence coercivities, which are then used preliminarily to determine the main ferromagnetic minerals of the natural materials (Wu et al., 2013b). The IRM of the majority of the samples, achieves 50% of the SIRM at~44 mT IRM acquisition curve decomposition is a mathematically statistical method using the linear additivity of the cumulative logarithmic Gaussian distribution of each component in the curve (Robertson and France, 1994) to efficiently distinguish magnetic mineral components (Kruiver et al., 2001;Heslop et al., 2002). We used the Excel workbook file provided by Kruiver (Kruiver et al., 2001) and the Irmunmix program written by Heslop (Heslop et al., 2002) to fit and optimize the first derivative of the IRM acquisition curve. The different B 1/2 value (i.e., the abscissa value of the peak), coercivity distribution and SIRM of each of the curves are shown as the criteria to discriminate individual magnetic components contained in the samples. Three components were recognized (Figures 4E, J, O) including component 1: maghemite or coarse magnetite particle (<20 mT mean coercivity), component 2: fine magnetite particle (~40-60 mT mean coercivity with wide distribution), and component 3: hematite (B 1/2 > 100 mT), which is consistent with the κ-T results. Component 3 accounts for only 11% and 4% of the total IRM of the samples S-150-10 and S-168-24 respectively, indicating the predominance of low-coercivity magnetic minerals and few highcoercivity magnetic minerals in the sediments and that highcoercivity magnetic minerals are few. In sample S-159-30, the IRM contribution of component 3 is about 29%. Therefore, lowcoercivity minerals coupled with considerable amounts of highcoercivity components make up the magnetic minerals in sediments from this depth.

Hysteresis loops and first order reversal curves (FORCs)
The hysteresis loops, corrected for paramagnetics (green lines in Figures 4C, H, M), provide comprehensive information about the types of magnetic minerals and the states of magnetic domains (particle sizes) (Roberts et al., 2000), which is mainly expressed by the shape of the loops and key data such as saturated magnetization (M s ), saturated remanence magnetization (M rs ), coercivity (B c ), remanence coercivity (B cr ), M rs /M s , and B c /B cr (Dunlop and Frontiers in Earth Science frontiersin.org Özdemir, 1997). All loops are thin and long in shape, making it clear that the magnetic particles included are relatively coarse with low coercivity ( Figures 4C, H, M). In addition, the hysteresis loops of some samples did not close at high field ( Figure 4C, M), which may be the consequence of the presence of paramagnetic minerals and diamagnetic components.
As an extension of hysteresis measurements, FORC diagrams are used to reveal the coercivity and domain state of the magnetic particles. Most samples have coercivity distribution between 0-100 mT with peak H c at 30 mT, and vertical distribution H u within the range of ±50 mT ( Figures 4D, I, N). The common shapes are triangles truncated by y-axis, which suggests that the main magnetic particles in these samples are PSD particles.

Stratigraphic power spectral analysis
Here we present the integrated spectral analyses of mass-specific magnetic susceptibility data, because χ data are closely related to lithologies and sedimentary facies and only this proxy provides the best cycle interpretation in the lower K 1 sh (detailed information see section 5.1). In the MTM power spectrum of χ data sequence, statistically significant peaks at 113 m, 34 m, 13 m and 6 m sedimentary cycles above 95% confidence level can be distinguished ( Figure 5A). The ratio of these cycles 18.8:5.7:2.2: 1 is consistent with the spectral peak ratio (18.4:5.7:2.1:1) for the astronomical periods of 405 kyr, 125 kyr, 47 kyr, and 22 kyr presented in the La2004 solution ( Figure 8A). This matching preliminarily provides an estimation of SAR of~27.20-27.90 cm/ kyr. The COCO tests of the χ data suggest the optimal sediment accumulation rates (OSARs) are 21.39 and 27.72 cm/kyr, with null H 0 (hypothesis) significance level <0.001 (Figures 5C-E). Thẽ 28 cm/kyr OSAR makes those cycles identified in depth domain can be interpreted to represent the long eccentricity, short eccentricity, obliquity, and precession cycles, respectively. The frequency shifts in the EHA spectrogram ( Figure 5B) are linked to the minor peaks centered around the major peaks found in the MTM spectrogram ( Figure 5A), which may be caused by the SAR changes over the stratigraphic domain. This was quantitatively estimated by the eCOCO method, which suggests that the SARs from the lower K 1 sh Formation most likely varied between 20 and 40 cm/kyr and centered on~28 cm/kyr (Figures 5F, G).
In addition, MTM spectral analyses were also performed on the magnetic indexes χ ARM and SIRM ( Figures 6B, C) with a larger maximum frequency of 0.4 cycles/m. According to the above OSAR,  (Kruiver et al., 2001;Heslop et al., 2002): red lines are the fitted results of the scattered blue diamonds which denote the gradients of IRM acquisition curves. Magenta, green, and grey shadings represent the maghemite (and/or coarse magnetite), fine magnetite, and hematite, respectively. The B 1/2 and DP values of each component have been labeled.

Frontiers in Earth Science
frontiersin.org the results of the power spectral analyses of the three magnetic parameters are shown below.
The power spectrogram of χ shows long eccentricity period at 113 m, short eccentricity at 34 m, obliquities at 13.2 m, precession at 6.7 and 6.0 m, and semi-precession at 2.67-3.21 m ( Figure 6A). The MTM spectrum of χ ARM are likely to reveal 514 kyr at 144 m, short eccentricity at 34 m, obliquity at 10.2 m, precession at 5.4-6.1 m, and semi-precession at 2.62-3.16 m ( Figure 6B). The prominent peaks discerned in the SIRM power spectrum at 104 m, 35 m, 9.3 m, 4.7-6.1 m, and 2.51-2.89 m correspond to the 371 kyr, short eccentricity, obliquity, precession, and semi-precession cycles, respectively ( Figure 6C). The interpreted 514-kyr cycle observed on the MTM power spectra of χ ARM series, which is different than the 405-kyr eccentricity cycle, may be related to the variations for the SAR or the frequency decomposition of the eccentricity of the earth (Laskar et al., 2004).
The SARs constrained by 405-kyr eccentricity cycles change between 26.7 cm/kyr and 30.9 cm/kyr, with an average of 28.1 cm/ kyr ( Figure 7B). The SARs analyzed on a scale of 105-kyr eccentricity cycles range from 18 to 40 cm/kyr, with an average of 27.8 cm/kyr ( Figure 7F). The high SAR in this study compared with the 6.55 cm/ kyr SAR obtained from well SS4 for the K 1 sh Formation  may be attributed to the frequent faulting activities in basin center during the Early Cretaceous. The SARs of 28.1 cm/kyr and 27.8 cm/kyr obtained by astronomical tuning (Figures 7B, F) are basically consistent with those obtained by COCO and eCOCO

FIGURE 5
Power spectra, COCO and eCOCO analysis for the χ series of the lower K 1 sh Formation (A) 2π MTM power spectrum of the χ series. The solid red line, dashed green, purple, and blue lines represent the 50%, 90%, 95%, and 99% confidence levels (B) Evolutive harmonic analysis (EHA) results with a sliding window of 120 m and a step rate of 1 m. The red to blue colors represent high to low power normalized to 1. The dashed white lines marked with E, e, O and P show the orbital long eccentricity, short eccentricity, obliquity, and precession cycles, respectively (C-E) COCO analysis of the untuned χ series. The results are shown as the correlation coefficient spectrum (ρ) (C), H 0 (null) hypothesis test (D), and the number of contributing astronomical frequencies (#) (E). Tested SARs range from 1 to 100 cm/kyr (F, G) eCOCO analysis of the untuned χ series. The evolutionary ρ and evolutionary H 0 (null) hypothesis test are overlapped with SAR curves (black lines) based on 405-kyr tuning and 105-kyr tuning. Tested SARs range from 1 to 40 cm/kyr. The sliding window size is 150 m with a step of 1.0513 m. For both the COCO and eCOCO, the running step of the SAR is 0.28778 cm/kyr and the number of Monte Carlo simulations is 2000.

Frontiers in Earth Science
frontiersin.org analysis ( Figure 5C) of~28 cm/kyr, and also falls into the range calculated by radiometric ages (12.14-46.46 cm/kyr) ( Figure 1D). These results provide conclusive evidence for the cyclostratigraphic interpretation based on the χ data.

Magnetic susceptibility as paleoenvironmental proxy
The statistical analysis aims to discover the correlations between lithology and magnetic susceptibility in our specimens. The lithological ranks were distinguished by the χ values. Except for the sandy conglomerate (rank = 9) showing a normal distribution, all the other lithological ranks show skewed distributions. The χ distributions of mudstone (rank = 1), silty mudstone (rank = 2), fine sandstone (rank = 5), coarse sandstone (rank = 6), gravel-bearing fine sandstone and gravel-bearing coarse sandstone (rank = 7 and 8) are positively skewed, while the χ distributions of muddy siltstone and siltstone (rank = 3 and 4) and conglomerate (rank = 10) are negatively skewed (Figure 2A). This indicates that χ for all lithological ranks is fairly concentrated. Except for the outliers of lithological ranks of 7-10, fine mudstones has the highest χ values, coarse sandstones has the lowest ( Figure 2B). This inverse correlation within the range of 1-6 proves that the χ data for the lacustrine faces of the lower K 1 sh Formation is a reliable and continuous proxy for the lithological ranks, and consequently for paleoenvironmental change.
The sedimentary facies of the investigated K 1 sh section can be determined as three major development phases (Figure 3): 1) The Subset 1 contains thick conglomerates with higher magnetic parameters over the dark grey mudstones in the bottom, indicating a fan-delta phase evolved from the lake phase due to the strong rifting in SB. 2) The Subset 2 is dominated by lake deposits, which exhibits a humid climate and a weak hydrodynamic environment, with a propensity for shallowing upward in the bottom portion and deepening upward in the upper part. 3) The Subset 3 begins with a minor tectonic activity, which is followed by a lake environment that is shallowing upward. Coal seam intercalations suggest a lakeswamp plain facies. Conglomerates and coarse sands are evidence of advancing fan-delta phases with a decline in base level and lake recession due to increased aridity. In contrast, an expanding lakeshore with a rise in base level and higher humidity produces fine sands, silts, muds, and coal. The rock magnetic characteristics of the representative samples demonstrate the main magnetic carriers of the lower K 1 sh Formation are lower coercivity PSD ferrimagnetic magnetite particles and a small amount of high coercivity hematite (Figure 4). The fact that the χ ARM and SIRM series and the lake level are coherently fluctuating with χ series (Figures 3B-E) suggests that χ is dominated by fine-grained ferrimagnetic magnetite minerals and can be used for the paleoenvironmental changes reconstruction. Since only χ exhibits all astronomical periods on its power spectrum ( Figure 6A), χ is ultimately selected as the proxy for the subsequent temporal calibration.

Astronomically forced paleoclimate and paleoenvironment changes in Early Cretaceous Songliao Basin
Lacustrine sediments provide ideal geologic materials for studying climate cycles caused by orbital forcing. It has been evidenced that Milankovitch cycles, particularly long and short eccentricity cycles can be widely detected in global lacustrine sedimentary records (Kent and Clemmensen, 1996;Kent and Olsen, 1999;Olsen et al., 2003;Wu et al., 2009;2013a;Deenen et al., 2010;Li et al., 2017;Xu et al., 2019;Zhang et al., 2019;Du et al., 2020). Olsen et al. (2003) proposed that deep lake facies containing dark shales correspond to the maximum eccentricity in the Newark Basin, North America, while subaerial deposition of red mudstones was developed during the low eccentricity period. This agrees with the sedimentary model that the lacustrine/alluvial facies in the Green River Basin coincided with the   (Laskar et al., 2004), over the periods of 110-120 Ma (B) the 405-kyr tuned χ series, and (C) the 105-kyr tuned χ series. E1: the longer eccentricity (~2 Myr); E2: long eccentricity; e: short eccentricity; O1: the longer obliquity (~1.2 Myr); O2: obliquity; P: precession; P/2: semi-precession.

Frontiers in Earth Science
frontiersin.org maximum/minimum eccentricity periods, which was validated by Smith et al. (2010) and Aswasereelert et al. (2013). Among the multi-proxies used in this study, χ is the best indicator, because it can detect the Milankovitch signals completely, continuously and accurately ( Figures 5A, 6A). Two sets of Gaussian bandpass filters for long eccentricity (113 m) and short eccentricity (34 m) were applied to filter χ throughout the study interval ( Figure 9A). Coal seams are typically matched with the maximum of short eccentricity cycle, while coarse-grained sediments are matched with the lowest value of long eccentricity cycles ( Figure 9A). This proves that the presence of coarse clastic sediments such as FIGURE 9 (A) Depositional facies, lithological ranks, and χ data in depth domain of the lower K 1 sh, with the interpreted 405-kyr and 100-kyr cyclic signals. The Gaussian bandpass filters set for the interpreted 405-kyr (red) and 105-kyr (blue) cycles were 0.009 ± 0.002 cycles/m and 0.035 ± 0.007 cycles/m, respectively. Red "E 405 " represents 405-kyr long eccentricity, blue "e" represents~105-kyr short eccentricity (B) Upper part: 2π MTM power spectrum and the precession signals in 4,542-4,740 m interval with its amplitude modulation, the interpreted precession cycles (purple) were extracted with a passband of 0.15 ± 0.007 cycles/m; lower part: 2π MTM power spectrum and the precession signals in 4,918-5,380 m interval with its amplitude modulation, precession cycles (purple) were extracted with a passband of 0.166 ± 0.003 cycles/m (C) The precession and semi-precession cycles of χ data of 4,650-4,670 m interval. Passband width of precession cycles (solid purple) was the same as (B), passband of 0.319 ± 0.003 cycles/m was applied to obtain the semi-precession cycles (dashed green). The yellow shades represent 4,542-4,740 m and 4,918-5,380 m intervals, the sand shades represent 4,650-4,670 m interval. Black dashed lines associate the coal seams with the 105-kyr filtering maximum.

Frontiers in Earth Science
frontiersin.org sandstones is likely due to low lake level with low rainfall, which generated in the low eccentricity and low insolation periods. The coal deposition occurred during a period of high 105-kyr eccentricity, indicating that the sediment was developed under warm, humid, and vegetated environments with high precipitation and high lake level. Precession signals with wavelength of 6.68 m and 6.01 m are highly visible between~4,740-4,542 m and~5,380-4,918 m intervals, respectively ( Figure 9B). The AM analyses indicate that~4-6 precession cycles were modulated by one envelope cycle, which is consistent with the 5:1 ratio of short eccentricity cycles to precession cycles ( Figure 9B). This is in line with the depositional model in the Newark Basin, and suggest that the lake level is forced by precession cycles that is modulated by the short eccentricity cycles (Olsen et al., 2003;Kent et al., 2017), confirming the astronomic origin of the high frequency cycles found in the study interval.
High-resolution laminated interbedded sandstone-mudstone alterations are regularly cycled between 4,650 and 4,670 m in lake deposition facies (Figures 9B, C). This suggests that the seasonality controlled by precession influences summer precipitation and is supposed to be the main driving mechanism for lake level oscillation and paleoenvironment variations (Morrill et al., 2001). When the Earth moves to perihelion at the summer solstice for the Northern Hemisphere, the strong seasonality (i.e., hot summer and cool winter) increases rainfall in summer, which leads to a high lake level that favors mudstone deposition. Alternatively, the perihelion occurs at boreal winter solstice (i.e., aphelion occurs at boreal summer solstice), the weak seasonality favors the deposits of coarse clastic sediments in cool summer.

Possible origins of the semi-precession cycles
Litho-phase alternations combined with the precession and semi-precession signals ( Figure 9C), untuned χ spectrum ( Figures  6A, 9B), and tuned χ spectra ( Figures 8B, C) all reveal the presence of semi-precession cycles in the interval studied. Semi-precession forcing of lacustrine deposition has already been demonstrated from the Late Triassic-Early Jurassic Newark Basin (Olsen and Kent, 1996) and the early Eocene Green River Basin, North America (Pietras et al., 2003). The semi-precession peaks observed in the 405-kyr and 105-kyr tuned χ spectra (Figures 8B, C) are not observed in the power spectrum of the La2004 solution ( Figure 8A). It may indicate that the~10-kyr signals are non-linear climatic responses to the Milankovitch insolation forcing during this period. AM analysis was applied to the semi-precession signals derived from the 405-kyr tuned χ series to investigate the nonlinear link between the interpreted semi-precession cycle and Milankovitch cycles ( Figure 10C). The MTM power spectrum of the AM envelope curve reveals peaks at 23 kyr, 33 kyr, 47 kyr, 55 kyr, 68 kyr, 84 kyr, and 110 kyr ( Figure 10E), indicating that precession, obliquity, and eccentricity cycles impacted the semi-precession cycles indirectly. Furthermore, we also performed AM envelope to the precession ( Figure 10B) which is regarded as the modulator of semi-precession cycles. The MTM power spectrum of precession AM analysis curve demonstrates that the semi-precession cycles were additionally influenced by~400 kyr and~1 Myr long eccentricity cycles ( Figure 10D).
It is well-known that semi-precession is a rather sensitive signal in tropical regions, which is modulated by precession, as the Sun passes the equator twice annually (Verschuren et al., 2009;FIGURE 10 Amplitude modulations of the (C) semi-precession and (B) precession signals from the 405-kyr tuned χ series (A), and the corresponding 2π MTM spectra of the individual amplitude modulation curve (D, E).

Frontiers in Earth Science
frontiersin.org Anderson, 2011;Ulfers et al., 2022). Monsoonal amplifications caused by the maximum summer insolation gradient in both hemispheres appear to be the origin of semi-precession forcing in low latitudes (Verschuren et al., 2009;De Vleeschouwer et al., 2012;Jian et al., 2020). However, SB located at the mid-latitudes of the Northern Hemisphere in the Early Cretaceous also captured the semi-precession signals. This supports the theory mentioned by previous research that the heat and energy advective processes play a considerable role in transferring this sub-Milankovitch signal from low latitudes to high latitudes (McIntyre et al., 1989;Short et al., 1991). This causes the tropical pattern to occur in the northern temperate paleoenvironment as well. The flux of magnetic material not only reaches its maximum during precession minima but also intensifies during precession maxima (De Vleeschouwer et al., 2012). Low summer insolation received for northern hemisphere when the Earth moves near perihelion at equinoxes weakens the corresponding rainy season (Berger et al., 2006;Verschuren et al., 2009). We speculate that teleconnections such as atmospheric circulation and ocean current dynamics can explain how this sub-Milankovitch cycle is transmitted from low to high latitudes (Berger and Loutre, 1991;Short et al., 1991;Hagelberg et al., 1994;Arz et al., 1998;Ulfers et al., 2022). Furthermore, the detailed mechanism of climate or whether the semi-precession force is the direct high-latitude reaction still needs more extensive geologic data and model simulation at various sites located at different latitudes.

Conclusion
(1) Low coercivity ferrimagnetic magnetite is the dominant magnetic mineral with the minor presence of hematite in sediments of the lower K 1 sh. Given the prominent relationship between lithological ranks and magnetic susceptibility, the mudstones have higher χ values than sandstones, the differences in magnetic concentrations in the sediments are most likely caused by variations in hydrodynamic conditions (e.g., changes in lake level) caused by climate change. Therefore, magnetic susceptibility can be used as an alternative and reliable proxy for cyclostratigraphic studies to reflect the changes in paleoenvironment and paleoclimate in the Early Cretaceous SB.
(2) Cyclostratigraphic study indicates that the sedimentation environment of the lower K1sh was forced by multiscale astronomical orbital parameters. It suggests that long and short eccentricity controlled the litho-facies, with coarse detrital sediments deposited during the minima of long eccentricity, and coal seams developing during the maximum of short eccentricity. Meanwhile, precession and semiprecession cycles were the potential factors affecting the lake level. In particular, the mudstone-sandstone alternations are most likely controlled by the semi-precession cycles. (3) The duration of the 4542-5695 m interval of lower K1sh is estimated as 4,090 kyr-4,148 kyr, which is consistent with the radiometric range of 118.2-113.9 Ma, providing a basis for the future correlation of stratigraphic evolution and geological events between marine and land systems. The estimated sedimentary accumulation rate of~28 cm/kyr, may support an abundant supply of sediments to the lake basin during the K 1 sh deposition period.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.