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². 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).

As a superimposed rift basin, the SB overlies a basement that consists mainly of Precambrian to Paleozoic metamorphic and igneous rocks and Paleozoic to Mesozoic granites. The Mesozoic and Cenozoic continental strata as thick as 10,000 m are recovered in the central part of the basin (Figure 1E) (Wang et al., 2006; Pei et al., 2007; Zhou et al., 2009; Zhang et al., 2011; Hou et al., 2018). The Cretaceous lithologic strata are divided into ten formations, namely, Huoshiling (K₁h), Shahezi (K₁sh), Yingcheng (K₁yc), Denglouku (K₁d), Quantou (K₁q), Qingshankou (K₂qn), Yaojia (K₂y), Nenjiang (K₂n), Sifangtai (K₂s), and Mingshui (K₂m) in ascending order (Wang et al., 2007; 2013; Feng et al., 2010) (Figure 1E).

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 (Wang et al., 2013). Based on the literature references referring to geology (Wang et al., 1995; Du et al., 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.

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₁d down to the basement of SB, with K₁yc, K₁sh, and K₁h as the main target strata, which were deposited during the Early Cretaceous syn-rifting stage. Unlike the volcaniclastic sequences in K₁yc and K₁h, the intercalated K₁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₁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₁sh unconformably overlies the K₁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₁sh of the SK-2 into six sedimentary cycles according to detailed core description, with each of them displaying fining-upward trend that reveals multi-periodic 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.

Previous investigations have extensively estimated the chronostratigraphy for K₁sh strata distributed in northeastern China basins by using diverse dating techniques. In the view of biostratigraphically, K₁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₁yc and 147.77 ± 0.28 Ma at the top of the K₁h using ⁴⁰Ar/³⁹Ar dating method, providing a ∼10 Myr duration for K₁sh. Liu et al. (2020) constructed a 11.14 Myr floating astronomical time scale for the K₁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₁yc performed by Shu et al. (2007) and Huang et al. (2011) respectively, demonstrated the base of K₁sh is supposedly formed during Aptian period. Recently, more new zircon U-Pb data derived from the SK-2 borehole restrict K₁sh into ∼118–111 Ma (An, 2018; Liu et al., 2020; 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₁sh in present research.

A total of 1,097 samples were collected from the lower K₁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³ 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.

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³/kg for χ (Figure 3D), 17.64×10^-8 m³/kg for χ_ARM (Figure 3C), and 13.94×10^-5Am²/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.

Subset 2 (5,400 m–4,900 m): The mean values of χ, χ_ARM and SIRM are 6.67×10⁻⁸ m³/kg, 7.28×10⁻⁸ m³/kg and 6.21×10⁻⁵ Am²/kg respectively. They are comparable to the bottom of Subset 1 and remain steady with periodic fluctuations (Figures 3B–D). While χ and χ_ARM do not exhibit apparent trends (Figures 3C, D), SIRM has a clearly increasing trend (Figure 3B), indicating high coercivity magnetic minerals increases during this interval.

Subset 3 (4,900 m–4,542 m): Mean values of χ, χ_ARM, and SIRM are 6.87×10⁻⁸ m³/kg, 6.55×10⁻⁸ m³/kg and 6.13×10⁻⁵ Am²/kg. The value of χ is comparable to that of Subset 2 and χ_ARM and SIRM are lower than those of Subset 2 (Figures 3B–D). χ_ARM and SIRM show distinctly declining trends with shallowing depth (Figures 3B, C).

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₁sh Formation is a reliable and continuous proxy for the lithological ranks, and consequently for paleoenvironmental change.

The sedimentary facies of the investigated K₁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 lake-swamp 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₁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.

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; 2014; 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 sub-aerial 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 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 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.

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 non-linear 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; 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.