The sill samples were collected from the Huailai section (N 40°28.7´; E 115°23.4′), Huangtugui section (N 40°34.5´; E 118°33.3′), Zhenzhuquan section (N 40°33.3´; E 116°21.3′), and Chaoyangxian section (N 40°55.7´; E 120°23.6′). Additionally, nineteen diabase samples were collected for mineralogical analysis from Huailai, Huangtugui, Zhenzhuquan and Chaoyangxian section. Twenty shale samples were collected for equivalent vitrinite reflectance testing and ten samples for TOC analysis from the Huailai section; nineteen shale samples were collected for equivalent vitrinite reflectance testing and ten samples for TOC analysis from the Huangtugui section; and eleven shale samples were collected for equivalent vitrinite reflectance testing and five samples were collected for TOC analysis from the Zhenzhuquan section. The distance of the shale samples to their surrounding rocks can be seem in the Supplementary Tables S1, S2. The weathered surfaces of all samples were cut to ensure that the samples are fresh and free of contamination, then the shale samples were ground to 100 mesh. All the experiments were conducted at the MOE Key Laboratory of Orogenic Belts and Crustal Evolution, School of Earth and Space Sciences, Peking University.

The detailed petrographic analysis of rock thin sections was conducted by the polarizing microscope (Nikon Eclipse LV100NPOL) and the scanning electron microscope (SEM, FEI Quanta 650 FEG) equipped with energy-dispersive X-ray spectroscopy (EDS, Oxford INCA X-Max 50).

To analyze the total organic carbon (TOC) content of shales, the inorganic carbon from 100 mg powdered sediments was removed by treatment with dilute HCl. The resulting residue was then repeatedly washed with deionized water and dried. Analysis of TOC was performed using a Leco CS-230 analyzer. Then, the samples were sectioned either perpendicular or parallel to the bedding and subsequently processed following the guidelines outlined in ISO 7404-2 (2009). Reflectance measurements, in accordance with ISO 7404-5 (2009) and ASTM D7708-14 (2014), were conducted using either the Leica 4500 P microscope with CRAIC microscope photometer or the Leica DM4 M microscope with the Hilgers Technisches Büero Fossil system, equipped with both white and ultraviolet light sources. The equivalent vitrinite reflectance (EqVRo) values were determined using vitrinite-like reflectance (Rv) data, applying the formula established by Luo et al. (2021).

In this study, we simulated the carbon release from sills by using the SILLi 1.0 1D model established by Iyer et al. (2018) with the help of the one-dimensional finite element method (FEM). The reliability of this method has been thoroughly validated by previous researchers (Svensen et al., 2018b; Iyer et al., 2018). In this model, each sedimentary layer, including eroded layers, is deposited sequentially in time based on the depositional age. The sedimentation rate of each layer is calculated based on the layer’s thickness and the time difference between its top age and the deposition age of the preceding layer (Svensen et al., 2018a). Deposition of an erosional layer occurs in the same way as the layers in the present-day sedimentary column. Sills are emplaced instantaneously at a given time. Based on the geothermal and energy diffusion equations, the temperature changes of each sedimentary layer in the stratigraphic column are calculated under specific boundary temperature conditions. The effective thermal capacity of the rock incorporates the latent heat of fusion associated with the crystallization process occurring within the sill (Aarnes et al., 2010; Iyer et al., 2013). The thermal diffusion equation is modified to incorporate dehydration reactions occurring in the surrounding rocks (Galushkin, 1997; Wang, 2012). The thermal maturity of each sedimentary layer is determined using the EASY%Ro algorithm, which simulates the complex process of kerogen decomposition with increasing temperature based on the Arrhenius equation (Sweeney and Burnham, 1990). Through the aforementioned process, we can obtain the thermal maturity for each stratigraphic horizon. Considering that the studied stratum belongs to the Precambrian period when higher plants had not yet appeared, it is not possible to obtain the vitrinite reflectance (Ro) of organic matter in the formation directly. The simulated vitrinite reflectance values (Ro) represent the equivalent vitrinite reflectance values (EqVRo) of organic matter in the studied strata. The equivalent vitrinite reflectance (EqVRo) values were calculated from vitrinite-like reflectance (Rv) based on the formula reported by Luo et al. (2021): EqVRo=1.07 × Rv-0.18. The release of organic CO₂ is estimated by assessing the TOC variation in the surrounding rocks based on the EASY%Ro method, with its maximum conversion value capped at 85% (Iyer et al., 2018). The release of inorganic CO₂ is estimated by using temperature and pressure curves of carbonate rocks over time, which are pre-calculated using Perplex_X (Connolly and Petrini, 2002). Note that the carbon gas release is calculated as CO₂-equivalents.

The sill sample have blocky structure, with some parts showing spheroidal weathering (Figure 4A). Under microscopic observation, the sample is identified as diabase with distinct diabasic texture and mineral particles in semi-automorphic granular structure. The main minerals include plagioclase (55–50 vol%), clinopyroxene (30–35 vol%), and amphibole (10–15 vol%). Plagioclase crystallizes earlier than clinopyroxene and exhibits a higher degree of self-shaping. The plagioclase shows slight epidotization alteration along its edges, while the clinopyroxene displays residual structures due to local amphibolization (Figures 4B,C).

At the Huailai section, within a 5.5 m range of the sill contact zone, as the distance increases, the colour of the surrounding rocks transitions from light green to dark gray-black, with an increase in chromaticity; the vitrinite reflectance values of the surrounding rocks decrease from 2.56 to 0.83, with a decrease rate of 67.58%; TOC values increase from 0.20 to 2.00, with an increase rate of 9 (Figure 5A). At the Huangtugui section, within a 3.0 m range of the sill contact zone, as the distance increases, the colour of the surrounding rocks transitions from light yellow to dark gray-black, with an increase in chromaticity; the vitrinite reflectance values of the surrounding rocks decrease from 2.86 to 0.87, with a decrease rate of 69.58%; TOC values increase from 0.03 to 2.45, with an increase rate of 80.67 (Figure 5B). At the Zhenzhuquan section, within a 2.5 m range of the sill contact zone, as the distance increases, the colour of the surrounding rocks transitions from light yellow to dark gray-blue, with an increase in chromaticity; the vitrinite reflectance values of the surrounding rocks decrease from 1.90 to 0.57, with a decrease rate of 70%; TOC values increase from 0.23 to 1.50, with an increase rate of 5.52 (Figure 5C).

Input data (Supplementary Tables S3–S11) for the numerical thermal sill model are based on nine stratigraphic columns from the Yanliao LIP (Figure 3). In order to obtain a more conservative estimation, the model utilized an emplacement scenario where individual sills were placed every 10,000 years. If the sills are placed close together, choosing to emplace them simultaneously in the model would cause the aureoles to experience a higher maximum temperature. This could potentially result in an overestimation of the CO₂ production. Pure limestone, despite containing large amounts of CO₂, has a thermal decomposition temperature that surpasses 950°C, exceeding that of impure carbonates (Ganino and Arndt, 2009; Heimdal et al., 2018; Iyer et al., 2018). Given that contact aureoles often cannot attain such high temperatures (Ganino and Arndt, 2009; Supplementary Figure S1), impure carbonates, such as marls, are more inclined to produce substantial volumes of CO₂ during sill emplacement.

Overall, the thermal simulation results indicate that sill emplacement increases the temperature (To ∼234°C in average) and vitrinite reflectance (%Ro) (From ∼0.50% to ∼2.37% in average) of the surrounding rocks, while significantly decreasing its total organic carbon content (TOC) (From ∼1.10% to ∼0.23% in average) (Supplementary Table S12). By comparing the simulation results with field measured results, at the Huailai section, a 3.0 m thick sill emplacement resulted in a decrease in vitrinite reflectance values of the surrounding rocks from 4.00 to 0.75 within a 5.5 m range, showing a consistency (R ²) of 0.89 with the field measurements; the TOC values within the 5.5 m range of the surrounding rocks increased from 0.00 to 1.11, showing a consistency (R ²) of 0.82 with the field measurements (Figure 5A). At the Huangtugui section, a 0.5 m thick sill emplacement resulted in a decrease in vitrinite reflectance values of the surrounding rocks from 4.00 to 0.70 within a 3.5 m range, showing a consistency (R ²) of 0.85 with the field measurements; the TOC values within the 3.5 m range of the surrounding rocks increased from 0.00 to 2.50, showing a consistency (R ²) of 0.91 with the field measurements (Figure 5B). At the Zhenzhuquan section, a 2.5 m thick sill emplacement resulted in a decrease in vitrinite reflectance values of the surrounding rocks from 4.00 to 0.80 within a 2.5 m range, showing a consistency (R ²) of 0.87 with the field measurements; the TOC values within the 2.5 m range of the surrounding rocks increased from 0.00 to 1.50, showing a consistency (R ²) of 0.77 with the field measurements (Figure 5C).

Besides, The peak organic CO₂ flux of this set of sill vary from 4.12 kg/m²/yr (1,321 Ma) in Chaoyangxian to 4027.78 kg/m²/yr (1,323 Ma) in Huailai, with an average value of 757.00 kg/m²/yr. The peak inorganic CO₂ flux of this set of sill vary from 63.70 kg/m²/yr (1,324 Ma) in Chaoyangxian section to 1,334.46 kg/m²/yr (1,316 Ma) in the same section, with an average value of 21.64 kg/m²/yr. The CO₂ production is calculated by subtracting the CO₂ value from background maturation from the peak CO₂ value after sill emplacement (Heimdal et al., 2018). The cumulative organic CO₂ production of this set of sill vary from 0.34 ton/m² in Zhenzhuquan section to 1,177.14 ton/m² in Chaoyangxian section, with an average value of 143.73 ton/m². The cumulative inorganic CO₂ production of this set of sill vary from 1.03 ton/m² in Lingyuan section to 58.35 ton/m² in Chaoyangxian section, with an average value of 21.64 ton/m². The cumulative total CO₂ production of this set of sill vary from 0.34 ton/m² in Zhenzhuquan section to 1,235.49 ton/m² in Chaoyangxian section, with an average value of 150.94 ton/m² (Table 1).

Next, we will analyze the thermal effects and cabon release of sill emplacement based on three stratigraphic columns located in the southwest (Huangtugui section), central (Pingquan section), and northeast (Chaoyangxian section) regions of the Yanliao LIP.

Considering the early formation age of the sill emplacement layers and their subsequent experiences of multiple tectonic transformation, it is challenging to reconstruct the outcrop’s exposure and burial history. Taking into account the measured shale in the surrounding rock, which was not affected by the emplacement, with low thermal maturity (EqVRo≤ 0.7%) (Liu et al., 2022), it indicates that the formation emplaced by sills in this area has been in a state of prolonged uplift or shallow burial, thus its burial history after Cambrian period has not been reconstructed. Figure 5 has demonstrated the reliability of the SILLi 1.0D model in simulating the thermal effects of the single-set sill. In order to validate the accuracy of our simulation results of multi-set sills emplacement in study area, we collected equivalent vitrinite reflectance data from 85 sets of sill emplacement surrounding rock samples from Huailai, Huangtugui, Zhenzhuquan and Pingquan section (Liu et al., 2011; Liu et al., 2022 and this study) and TOC data from 38 sets of sill emplacement surrounding rock samples from Huailai, Huangtugui, Zhenzhuquan and Kuancheng section (Shui et al., 2020 and this study) (Supplementary Tables S1 and S2). We analyzed the consistency between the simulated results of our model and the measured results of the samples. The results indicate a high consistency (R ²= 0.94) between the measured equivalent vitrinite reflectance and the simulated vitrinite reflectance (Figure 8A) and a high consistency (R ²= 0.97) between the measured TOC and the simulated TOC (Figure 8B). Prior researches has confirmed the precision of the SILLi 1.0 1D model (Svensen et al., 2018b; Iyer et al., 2018), and its simulated outcomes exhibit a strong consistent with the measured values of relevant indicators. This indicates that the lack of reconstructed burial history after the Cambrian period in this area does not affect the match between model predictions and actual conditions.

Based on the estimation of SILLi 1.0 1D model on nine typical stratigraphic columns (Figure 9), it was shown that the total cumulative of CO₂ production by sill emplacement into the surrounding rocks was 150.94 ton/m². Considering the distribution area of the Derim sill in Northern Australia Craton can reach up to 1.6 × 10⁵ km² (Kirscher et al., 2021; Mitchell et al., 2021), the total area of this set of sill could reach 3.0×10⁵ km². The sill emplacement resulted in an equivalent CO₂ release of approximately 4.53 × 10¹³ tons, equivalent to a carbon release of 1.24×10¹³ tons.

When using the 1.8 × 10⁶ tons/km³ CO₂ value for basaltic magma degassing from McCartney et al. (1990), the estimated gas release from volcanic eruptions in this LIP is 5.94 × 10¹¹ tons of CO₂ equivalents. The total gas release from this set of magma activity is estimated to be 4.58 × 10¹³ tons of CO₂ equivalents. This is equivalent to 1.53 times (3 × 10¹³ tons) of the CO₂ released by the Siberian Large Igneous Province (Grasby and Bond, 2023), 2.97 times (1.54 × 10¹³ tons) of the CO₂ released by the Central Atlantic Magmatic Province (Schaller et al., 2012), and 21.2 times (2.16 × 10¹² tons) of the current atmospheric concentration of CO₂ (IPCC, 2013; Jones et al., 2016). It is evident that the carbon release triggered by this set of magma activity is enormous.

In addition to the enormous release of carbon, sill emplacement into the formation can also generate toxic gases such as SO₂ and CH₃Cl through reactions with sulfur-rich, chlorine-rich, and other accessory minerals (such as pyrite) (Shui et al., 2020). The actual gas release during the magma activity period is expected to be greater than our estimated value.

In recent years, there has been considerable debate regarding the ∼1.4 Ga oxygenation event (Mesoproterozoic Oxygenation Event, MOE). Nonetheless, mounting evidence now suggests the occurrence of an oxygenation event in both the seawater and atmosphere around 1.4 Ga: (1) Widespread distribution of black shales around 1.4 Ga in North China, North Australia, and potentially in other cratons (Zhang et al., 2018; Diamond et al., 2021), suggesting significant organic matter burial and the subsequent accumulation of high oxygen in the Earth’s surface system; (2) The high Mo isotopic values (reaching up to ∼1.7‰) found in the shales of the middle Xiamaling Formation are comparable to those of modern seawater (∼2.3‰; Diamond et al., 2018), suggesting a more oxygen-rich ocean (Diamond et al., 2018; Diamond and Lyons, 2018; Zhang et al., 2019); (3) Isotope mass balance modeling of U from the northern Australia Craton during the same period suggests that over 75% of the seafloor was probably oxygenated around 1.4 Ga (Yang et al., 2017); (4) The significant negative shift in δ³⁴S_py in the upper Xiamaling Formation is understood as indicating higher seawater sulfate concentration (Wang et al., 2020; Zhang et al., 2022b), potentially attributed to increased chemical weathering on the continent due to elevated atmospheric O₂ levels (Daines et al., 2017); (5) The Fe speciation and concentrations of redox-sensitive metals in the green shales of the upper Xiamaling Formation similarly indicate oxygenated conditions in the bottom waters, suggesting a rise in marine oxygenation (Zhang et al., 2017b); (6) Elevated I/(Ca + Mg) values in the upper Xiamaling Formation (reaching up to 1.59 μmol/mol), suggesting oxygenated conditions in the surface seawaters. All this evidence suggests that while the oxygenation event ∼1.4 Ga may not have been global in scale, it likely represented at least a localized oxygenation event across the Yanliao Area (Xu et al., 2024), which lasting for 230 million years from 1.59 to 1.36 Ga (Zhang et al., 2022a; Figure 10).

The effects of Large Igneous Province (LIP) activity on Earth’s surface environment may have been largely consistent across both the Phanerozoic and Precambrian (Ernst and Youbi, 2017). The eruptions of LIPs would release substantial carbon gases into the atmosphere. The magmatism generate volcanic materials (fresh silicate material with high initial reactivity) and expand the surface area accessible for weathering (Jones et al., 2016; Ernst and Youbi, 2017), thereby providing substantial nutrients to the ocean, enhancing primary productivity (Horton, 2015; Diamond et al., 2021; Xu et al., 2024). Based on the positive Hg/TOC excursion, elevated nutrient levels (particularly P) and increased TOC in the middle Xiamaling Formation, Xu et al. (2024) suggested that the weathering of ∼1.4–1.35 Ga emplaced LIPs (with a maximum peak at ∼1.38 Ga, Zhang et al., 2022b) may have amplified the input of nutrients to the ocean, bolstering higher primary productivity, the burial of organic matter, and ultimately contributing to an oxygenation event (MOE) around 1.4 Ga.

It is worth noting that the ∼1.32 Ga sill (Yanliao LIP) and the ∼1.4–1.35 Ga LIPs are both attributed to Columbia continental rift system. The Yanliao LIP is specifically identified as the rift zone branch between North Australia and North China, which reached the drifting stage, ultimately leading to the fragmentation and final breakup of Columbia (Zhang et al., 2017a; Wang et al., 2019; Kirscher et al., 2021). Besides, the ∼1.32 Ga Yanliao LIP occurred closer in proximity to the North China Craton, which contains records of the MOE event, compared to the occurrence sites of ∼1.4–1.35 Ga LIPs (Figure 10A). Considering the calculated release of carbon gas content triggered by the ∼1.32 Ga sill emplacement above, it appears that this period of magmatic activity promote the activation of carbon in crustal sediments and may have also impacted the global environment. However, results from previous studies at that time indicate that atmospheric oxygen levels were nearly stable (Zhang et al., 2022a). Therefore, we believe that further discussion is needed on the global climate conditions during the emplacement of the ∼1.32 Ga sill. Unfortunately, due to the absence of co-depositional records in the North China Craton and the North Australian Craton during the pre-magmatic period, we lack further evidence to indicate the impact of the ∼1.32 Ga sill emplacement on the global ecosystem. Given the occurrence of Lower Riphean sedimentary formations linked to sill emplacement periods in various ancient cratons worldwide (eg. Siberian Craton), it is possible that future exploration of the geological records in these ancient cratons could reveal more information about the environmental impact of the ∼1.32 Ga sill emplacement.

The application of the SILLi 1.0 1D model to simulate the thermal consequences of sill emplacement offers additional understanding regarding the release of carbon and environmental ramifications linked to magmatic events (such as Yanliao LIP). However, the thermal simulation results still possess certain limitations: (1) The thermal simulation model is established based on one-dimensional finite element method and cannot simulate the temporal variation of the thermal effects generated by rock intrusion on a three-dimensional level; (2) The model only takes into account the influence of sill emplacement on the vitrinite reflectance (%Ro), total organic carbon content (TOC) and temperature of the surrounding rocks, and does not consider the relationship between other indicators in the formation (such as mineralogical and biomarker indicators) and sill emplacement; (3) The thermal simulation model assumes an approximate maximum conversion rate of 85% for TOC in sedimentary layers, but it does not take into account the complete loss of organic carbon in the contact zone, resulting in an underestimation of released gases. Additionally, the model does not consider the presence of inert kerogens in organic-rich layers, which reduces the amount of organic matter that can be converted into hydrocarbons (Iyer et al., 2017) and affects the transformation rate of TOC; (4) Due to the pre-magmatic activity, the co-depositional records in the North China Craton and the North Australian Craton absent during the sill emplacement. We are currently unable to estimate the amount of released CO₂ that is sequestered in reservoirs as oil and gas, or transformed into calcium carbonate cement in sedimentary rocks without being discharged into the atmosphere. How to overcome the limitations of the SILLi 1.0 1D model and obtain more precise constraints on carbon release estimation and environmental effects assessment caused by sill emplacement is an area worth exploring in future research for us.