2.1. Study site and ice core drilling

The ice core site is located on the SGP (73° 51.10′S; 163° 41.22′E; 1,623 m above sea level, a.s.l.), a 150 km² area situated along the western Ross Sea margin of the Transantarctic Mountains (TAM), ∼60 km inland from the Ross Sea coast, and ∼85 km away from the Korean Jang Bogo Antarctic Research Station in northern Victoria Land, East Antarctica (Figure 1). The SGP lies within a glaciated saddle area flanked by rugged TAM peaks, where ice flows southeastward through a valley outlet toward the Ross Sea. This region is characterized by prevailing southerly/southwesterly katabatic winds descending from the elevated TAM slopes (Udisti and others, Reference Udisti, Becagli, Castellano, Traversi, Vermigli and Piccardi1999), with air masses originating from the Ross Sea (Sinclair and others, Reference Sinclair, Bertler and Trompetter2010; Markle and others, Reference Markle, Bertler, Sinclair and Sneed2012; Ro and others, Reference Ro2022). Despite the dominance of katabatic winds in northern Victoria Land, the snow layers of the SGP are well-preserved due to the absence of summer melt (Udisti and others, Reference Udisti, Traversi, Becagli and Piccardi1998, Reference Udisti, Becagli, Castellano, Traversi, Vermigli and Piccardi1999; Nyamgerel and others, Reference Nyamgerel, Hong, Han, Kim, Lee and Hur2021) and the reduced influence of katabatic winds because the SGP is sheltered by surrounding topography (Stenni and others, Reference Stenni2000). The annual SARs are relatively high, ranging from 130 to 226 mm in water equivalent (w.e.) per year (Udisti, Reference Udisti1996; Stenni and others, Reference Stenni2000; Han and others, Reference Han2015; Kwak and others, Reference Kwak2015; Nyamgerel and others, Reference Nyamgerel, Han, Kim, Hong, Lee and Hur2020, Reference Nyamgerel, Hong, Han, Kim, Lee and Hur2021, Reference Nyamgerel2024; Ro and others, Reference Ro2022), making the SGP an ideal location for preserving high-resolution glaciochemical records. The mean annual temperature is −32.5°C, with a horizontal ice flow velocity of 0.9–2.3 m yr⁻¹ and an estimated ice thickness of 550 m (Han and others, Reference Han2015; Yang and others, Reference Yang2018; Kim and others, Reference Kim, Prior, Han, Qi, Han and Ju2020).

Figure 1. (A) Map of the Antarctic continent showing the locations of the Styx-M and Roosevelt Island Climate Evolution (RICE) ice cores from the Styx Glacier plateau (SGP) (163° 41.22′E, 73° 51.10′S, 1,623 m a.s.l.) and Roosevelt Island (161° 42.36′W, 79° 21.84′S, 550 m a.s.l.), respectively. (b) Enlarged scale of the western Ross Sea region showing the SGP, the Korean Jang Bogo Station (164° 13.7′E, 74° 37.4′S, 36 m a.s.l.), and various ice core sites: Hercules Névé (HN; 165° 24′E, 73° 6′S, 2,960 m a.s.l.), Talos Dome (TD; 158° 45′E, 72° 22′S, 2,316 m a.s.l.), and GV7 (158° 51.14′E, 70° 41.05′S, 1,947 m a.s.l.). Maps were modified from images generated in software QGIS 3.24.2 (https://www.qgis.org) using a visualization platform Quantarctica (Matsuoka and others, Reference Matsuoka2021).

The drilling operation was conducted from 12 December 2014 to 1 January 2015 as part of ice core drilling program of Korea Polar Research Institute (KOPRI) in Antarctica. A 210.5 m-long ice core, hereafter referred to as the Styx-M core, was recovered in 300 runs using a shallow ice coring system (Geotech Co. Ltd., Japan), with section lengths ranging from 19 cm to 100 cm (70 cm on average). Each section was weighed in the field to determine its bulk density, providing the density profile of the ice core as a function of depth. The ice core samples were sealed in polyethylene bags, transported frozen to KOPRI, and stored at −15°C in a cold room until further analysis. Additional details on the study site and drilling operations are reported elsewhere (Han and others, Reference Han2015; Yang and others, Reference Yang2018; Jang and others, Reference Jang2019; Kim and others, Reference Kim, Prior, Han, Qi, Han and Ju2020; Nyamgerel and others, Reference Nyamgerel2024).

2.2. Analytical method

This study focused on the top 98.5 m of the Styx-M core, for which high-resolution glaciochemical and isotopic analyses has been completed. Measurements included liquid conductivity and insoluble particle concentrations using a continuous ice core melting system, major ions (Na⁺, Ca²⁺, and SO₄^2–) using a two-channel ion chromatography (IC) system, and stable oxygen isotope ratios (δ¹⁸O) of meltwater using cavity ring-down spectroscopy (CRDS). All samples were decontaminated following the procedures described below, except for δ¹⁸O subsamples, cut parallel to the core axis without additional decontamination. Sampling resolution varied depending on the analytical procedures, as detailed in the following sections.

2.2.3. Ion chromatography system

The subsamples obtained from the top 12.2 m and the 50.9 m to 98.5 m depth interval were melted at room temperature in a class 10 laminar flow clean bench located in a cleanroom of class 1000. These samples were aliquoted and analyzed for major ionic species (Na⁺, Ca²⁺, and SO₄^2–) using Dionex (Thermo Fisher Scientific, Sunnyvale, CA, USA) ICS-2100 and ICS-5000 with an IonPac CS12A (2 × 250 mm) analytical column, methanesulfonic acid eluent, and a CSRS-300 suppressor for Na⁺ and Ca²⁺ and Dionex ICS-2000 and ICS-5000 equipped with an IonPac AS15 (2 × 250 mm) analytical column, potassium hydroxide eluent, and an ASRS-300 suppressor for SO₄^2–. The method detection limits, defined as three times the standard deviation (3σ) of 10 measurements of standard solutions (0.4, 1, and 5 μg L^–1), were 1.2, 0.46, and 0.58 μg L^–1 for Na⁺, Ca²⁺, and SO₄^2–, respectively.

The major ions in meltwater from a 12.2 m to 50.9 m depth were analyzed using the online multi-IC system (see section 2.2.2), comprising six Dionex IC sets, ICS-5000 and two ICS-1100s installed with an IonPac CS12A (3 × 150 mm) analytical column, methanesulfonic acid eluent, and a CERS-500 suppressor for Na⁺ and Ca²⁺ and ICS-5000, ICS-2100, and ICS-2000 installed with an IonPac AS15 (3 × 150 mm) analytical column, potassium hydroxide eluent and an AERS-500 suppressor for SO₄^2–. The method detection limits (3σ of 10 measurements of 0.2, 0.5 and 25 μg L^–1 standard solution) were 0.08, 0.06 and 3.0 μg L^–1 for Na⁺, Ca², and SO₄^2–, respectively, with a mean sample resolution of 1.8 cm.

The ion concentrations were calibrated using diluted standard solutions from Dionex (Cat No. P/N 046070 for Na⁺ and Ca²⁺ and P/N 057590 for SO₄²⁻). Data validation was carried out using the calibration parameters and ionic balances to ensure the accuracy and quality of data. Data suspected of being contaminated were excluded after the measurement was completed. The analytical precision was 4.9% for Na⁺, 6.3% for Ca²⁺ and 5.6% for SO₄²⁻. Further details on these analytical methods are provided elsewhere (Hong and others, Reference Hong2012, Reference Hong2015; Ro and others, Reference Ro2020, Reference Ro2022). Table S1 of the supplementary material lists all ion data.

In this study, non-sea salt (nss–) fractions of Ca²⁺ and SO₄²⁻ were calculated using the following equations (Delmas, Reference Delmas1992):

(1)\begin{equation}{\lbrack\mathrm{nssCa}^{2+}\rbrack}_{\mathrm{sample}}=\;{\lbrack\mathrm{Ca}{}^{2+}\rbrack}_{\mathrm{sample}}\; - \;{\lbrack\mathrm{Na}^+\rbrack}_{\mathrm{sample}}\;\times\;{\lbrack\mathrm{Ca}^{2+}/\mathrm{Na}^+\rbrack}_{\mathrm{seawater}}\;\end{equation}

(2)\begin{align}\left[{\text{nssSO}^{2-}_{4}}\right] _{\text{sample}} & = \left[{\text{SO}^{2-}_{4}}\right]_{\text{sample}} - \left[\text{Na}^{+}\right]_{\text{sample}}\nonumber\\ &\quad\, \times \left[{\text{SO}^{2-}_{4}}/ \text{Na}^{+}\right]_{\text{seawater}}\end{align}

where [Ca²⁺]_sample, [Na⁺]_sample and [SO₄^2–]_sample are the measured concentrations of total Ca²⁺, Na⁺ and SO₄^2–, respectively, in the samples. [Ca²⁺/Na⁺]_seawater and [SO₄^2–/Na⁺]_seawater are the ratios of Ca²⁺/Na⁺ (0.038 w/w) and SO₄^2–/Na⁺ (0.25 w/w) in bulk seawater (Kester and others, Reference Kester, Euedall, Connors and Pytkowicz1967). Note that Na⁺ can originate from both sea salt and continental dust, and therefore, some studies first estimate the sea-salt-derived Na⁺ budget before calculating the nss-fractions of other ions (e.g., Bigler and others, Reference Bigler, Rothlisberger, Lambert, Stocker and Wagenbach2006). However, it is widely accepted that Na⁺ in Antarctic aerosols, snow, and ice is predominantly of marine origin, with negligible contribution from dust-derived leachable Na⁺, particularly during interglacial periods. For example, the nss-fraction of total Na⁺ at Dome C during the Holocene was estimated to be only around 2% (Röthlisberger and others, Reference Röthlisberger2002). Therefore, Na⁺ concentrations can reasonably be interpreted as a proxy for sea salt input in this study.

2.3. Calculation of snow accumulation rate

Based on the established ice chronology, the annual SAR at the SGP was calculated in millimeters of water equivalent per year (mm w.e. yr^–1) using the density of the core sections evaluated by weighing each section directly in the field (see section 2.1) and a correction for the layer thickness. This correction was accomplished using thinning functions representing the reduced thickness (λ) to the initial thickness (λ_H) ratio of an annual layer. The SGP was characterized by low horizontal ice flow velocities (0.9–2.3 m yr⁻¹) (Yang and others, Reference Yang2018; Kim and others, Reference Kim, Prior, Han, Qi, Han and Ju2020). Kim and others (Reference Kim, Prior, Han, Qi, Han and Ju2020) also reported that the depth interval analyzed in this study was affected by a simple vertical strain, whereas the strain regime became complex below ∼ 110 m due to the bedrock topography near a nunatak located ∼ 4 km southeast of the Styx-M core site. Accordingly, the effects of the vertical strain-induced layer thinning were accounted for using linear least squares fitting models developed by Nye (Reference Nye1963) (Nye model) and Dansgarrd and Johnsen (Reference Dansgarrd and Johnsen1969) (D–J model). The Nye model assumes a constant vertical strain rate throughout the entire depth of the ice sheet, while the D–J model assumes a constant vertical strain rate from the surface to a certain depth (kink height, h) and a linear decrease in vertical strain below this depth to zero at the base of the ice sheet (H). The thinning function used for the SAR calculation was derived by averaging six thinning scenarios, incorporating the Nye model and D–J model with five different h (100, 200, 300, 400, and 500 m), consistent with the parameters used by Yang and others (Reference Yang2018). The uncertainty of the thinning function was assumed to increase linearly with depth at a rate of 0.06% m^–1 (1σ of the thinning functions calculated from the six scenarios divided by the corresponding depth). Additional uncertainty associated with the assumption of simple vertical strain at the SGP cannot be excluded. Despite these uncertainties, the thinning functions for depths between 68.1 m and 98.2 m showed good agreement with those derived from gas-based ice ages and their corresponding depths, which were also influenced by the thinning effects (Figure S3), suggesting that the layer thinning corrections had been properly performed.

3.1. Ice core chronology (1259–2014 CE)

Manual layer counting was extended to a depth of 98.5 m using the previous approach building upon the previously established chronology for the top 9.89 m (1979–2014 CE), determined through annual layer counting based on the seasonally varying profiles of δ¹⁸O, nssSO₄^2–, and nssSO₄^2–/Na⁺ with summer maxima and Na⁺ and liquid conductivity with winter/spring maxima (Nyamgerel and others, Reference Nyamgerel, Han, Kim, Hong, Lee and Hur2020, Reference Nyamgerel2024; Ro and others, Reference Ro2020, Reference Ro2022). A key methodological difference in this study compared to previous studies is the exclusion of δ¹⁸O as a seasonal marker, due to the significant signal smoothing below 9.89 m likely caused by post-depositional alteration and molecular diffusion in the firn (Johnsen and others, Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander and Creyts2000; see Figures 2 and S4). Figures 2 and S4 show the manual identification of annual layers based on the nssSO₄^2–, nssSO₄^2–/Na⁺, Na⁺, and liquid conductivity profiles. These impurity records exhibited apparent seasonal variations, as shown in Figure 3. Although the sampling resolution was sufficiently high in the upper 50.9 m (∼8 samples per year) to capture seasonal signals in the ion records, the interval below this depth had a lower sampling resolution (∼4 samples per year). Considering the effects of ice-thinning at greater depths, the ionic profiles may exhibit seasonal bias, particularly toward seasons with higher snow accumulation. Despite this, the impurity records exhibited well-defined seasonal patterns, highlighting their reliability as seasonality indicators across the core section (see Figures 2 and 3).

Figure 2. Profiles of δ¹⁸O, nssSO₄^2–, nssSO₄^2–/Na⁺, Na⁺, liquid conductivity, dust number concentrations (uncalibrated), and nssCa²⁺ at the depth intervals of (a) 2–7 m (from Nyamgerel and others, Reference Nyamgerel2024), (b) 17–22 m, and (c) 54–58 m of the Styx-M core. For the full depth profiles, see Figure S4. Horizontal dotted grey lines represent the mid-summer (1 January) of each year, which was assigned based on the ion data points. Horizontal dashed blue line in (b) represents an uncertain layer having potentials but not assigned as an annual layer. Horizontal dashed red lines represent the depth intervals showing the anomalous peaks in nssSO₄^2– coinciding with high concentrations of Na⁺ (see section 3.1).

Figure 3. Seasonal (summer, December–February; autumn, March–May; winter, June–August; and spring, September–November) mean values of the chemical impurities used for counting annual layers in the Styx-M core records. Annual signals in the impurity records are used to calculate a three-month seasonal mean after interpolating values at monthly intervals over a one-year period. The ion records from the depth intervals exhibiting anomalous, simultaneous peaks in Na⁺ and nssSO₄^2– concentrations were excluded

During layer counting, challenges were encountered at specific depth intervals, where multiple nssSO₄²⁻ peaks appeared within a single year (e.g., around 17.4, 18.2, and 21.8 m; see Figure 2b). This may result from sporadic biogenic sulfur inclusions from the adjacent oceans (Goursaud and others, Reference Goursaud2017; Nardin and others, Reference Nardin2021; Emanuelsson and others, Reference Emanuelsson, Thomas, Tetzner, Humby and Vladimirova2022; Hoffmann and others, Reference Hoffmann2022). In such cases, summer layers were determined based on troughs in winter markers or subjectively decided when seasonal signals were ambiguous to optimize consistency with absolute age markers (see below). In addition, nssSO₄²⁻ concentrations typically exhibited depleted to negative values during colder seasons (see Figure 2), likely due to sulfate depletion in sea salt aerosols from the sea ice surface, caused by mirabilite precipitation (Na₂SO₄ · 10H₂O) (Wagenbach and others, Reference Wagenbach1998; Legrand and others, Reference Legrand, Preunkert, Wolff, Weller, Jourdain and Wagenbach2017). On the other hand, sporadic nssSO₄²⁻ peaks were observed within one to two data points, accompanied by abnormally high Na⁺ concentrations (often exceeding 1 mg L^–1), which are likely linked to storm events originating from the Ross Sea (Ro and others, Reference Ro2022). These anomalies did not correspond to known volcanic events and were found alongside the troughs in the δ¹⁸O and nssSO₄^2–/Na⁺ profiles, which are indicative of cold seasons (e.g., at a depth of 5.24 m; see Figure 2a). Such spikes have also been recorded in snowpit/firn cores from the Victoria Land glaciers (Gragnani and others, Reference Gragnani, Smiraglia, Stenni and Torcini1998; Stenni and others, Reference Stenni2000; Kwak and others, Reference Kwak2015; Ro and others, Reference Ro2022). Although the intermittent influences of unidentified volcanism in northern Victoria Land cannot be excluded, we propose an alternative explanation for these phenomena based on the typical concurrent spikes in insoluble particles and nssCa²⁺: an abnormal influx of sulfate-containing particles such as gypsum (CaSO₄ · 2H₂O), calcium carbonate (CaCO₃) reacted with atmospheric H₂SO₄, or H₂SO₄ hydrates. Gypsum and carbonate minerals can be supplied from local rock outcrops in the TAM (Iizuka and others, Reference Iizuka2013) and deposited anomalously at the SGP under strong wind conditions. H₂SO₄ hydrates formed in polar stratospheric clouds during the colder seasons could be transported downslope via katabatic flows (Carslaw and others, Reference Carslaw1998; Watanabe and others, Reference Watanabe, Sato and Takahashi2006), consistent with the geographical features of the SGP located downhill from the TAM (see section 2.1). However, a detailed investigation of these processes is beyond the scope of this study. Therefore, nssSO₄²⁻ peaks coinciding with Na⁺, nssCa²⁺, and insoluble particle spikes were attributed to the colder seasons (highlighted by dashed red lines in Figures 2 and S4 and blue stars in Figure 4).

Figure 4. Profiles of nssSO₄^2– concentration (black) and ²³⁹Pu concentration (green; Shin and others, Reference Shin, Lee, Han, Hwang, Ro and Hur2025) at the top 98.5 m of the Styx-M core. Temporal horizons include known (red numbers and letters) and unknown (red stars) major volcanic signals and Rittmann tephra layer (brown; Lee and others, Reference Lee, Kyle, Iverson, Lee and Han2019) (see Table 1). Red numbers (#1–5) indicate known volcanic events used as age-depth tie points (see text), and their dates are listed in the box below the figure. Red letters denote known eruptions whose timing is inferred exclusively from chronological considerations and thus are not utilized as tie points (a: 1955 CE Carrán-Los Venados, b: 1886 CE Tarawera, c: 1846 CE Armagura, d: 1835 CE Cosiguina, e: 1822 CE Gallungung, f: 1673 CE Gamkonora, g: 1640 CE Parker Peak). RMN and SD are the running mean and standard deviation, respectively, of nssSO₄^2– biogenic background levels after removing values greater than the 95% percentile in the total nssSO₄^2– concentrations. Blue stars represent the depth intervals showing the anomalous peaks in nssSO₄^2– coinciding with high concentrations of Na⁺ (see section 3.1). The signals of the Pinatubo (N1) and El Chichón (N2) eruptions were already identified in the Styx-M core by Nyamgerel and others (Reference Nyamgerel2024).

Table 1. Age-depth tie points assigned to constrain the styx-m core chronology

^a Numbers of each volcanic event marked in Figure 4.

^b V: Volcanic eruption; N.T.: Nuclear test.

^c Annual-layer dating using the StratiCounter layer-counting algorithm (Winstrup and others, Reference Winstrup, Svensson, Rasmussen, Winther, Steig and Axelrod2012).

^d From Nyamgerel and others (Reference Nyamgerel2024).

^e From Shin and others (Reference Shin, Lee, Han, Hwang, Ro and Hur2025).

^f From Lee and others (Reference Lee, Kyle, Iverson, Lee and Han2019).

The layer counting result was further refined using age-depth tie-points such as atmospheric nuclear tests indicated by the ²³⁹Pu record (Shin and others, Reference Shin, Lee, Han, Hwang, Ro and Hur2025) and major volcanic eruptions identified in the nssSO₄²⁻ record based on the following criteria (Nardin and others, Reference Nardin2020, Reference Nardin2021): (1) at least one nssSO₄²⁻ sample exceeding the running mean (RMN) plus three times the standard deviation (RMN + 3σ) of the biogenic nssSO₄²⁻ background data (produced by excluding the top 5% of total nssSO₄²⁻ concentrations) and (2) at least two consecutive data points above the RMN + 2σ (Figure 4). Previous studies identified two ²³⁹Pu peaks at 12.7 m and 14.3 m, corresponding to the nuclear tests conducted by the USSR in 1961–62 CE and the USA in 1952–54 CE, respectively (Hwang and others, Reference Hwang, Hur, Lee, Han, Hong and Motoyama2019; Shin and others, Reference Shin, Lee, Han, Hwang, Ro and Hur2025). These peaks were assigned to 1965 CE and 1955 CE, respectively (Table 1), considering the delayed fallout deposition caused by long-range transport and the assigned years of volcanic eruptions recorded near these ²³⁹Pu peaks (see below). In addition, 12 volcanic eruptions were identified in this study, along with the previously reported Pinatubo and El Chichón eruptions at depths of 5.71–6.59 m and 8.20 m, respectively (Nyamgerel and others, Reference Nyamgerel2024) (Figure 4). Although the nssSO₄²⁻ signal for the Agung eruption did not meet the volcanic identification criteria, it matched the annual layer counting and the ²³⁹Pu-derived date at 12.7 m, supporting the rigor of these criteria. We note here that although all volcanic events identified in this study have been previously reported, several volcanic signals are difficult to definitely link to specific events. As shown in Table S2, Antarctic ice cores commonly contain records of the well-known 1963 CE Agung, 1815 CE Tambora, 1809 CE unknown, 1600 CE Huaynaputina, 1458 CE unknown (likely from Kuwae), and 1257 CE Samalas eruptions, whereas several eruptions (e.g., 1955 CE Carrán-Los Venados, 1846 CE Armagura, and 1822 CE Gallunggung eruptions) have been detected in a limited number of ice cores. Interestingly, the 1809 and 1458 CE eruption signals were not detected in our records. The absence of signals is most likely due to high biogenic sulfur backgrounds at coastal sites or the uncertainty of the expected depth intervals (39.9‒40.2 m and 80.0‒80.7 m) where no data are available (see Table S1 and Figure S4). To ensure robust age-depth tie points, we only used the 1963 CE Agung, 1815 CE Tambora, 1600 CE Huaynaputina, and 1257 CE Samalas eruption signals (Table 1 and Figure 4). Based on these absolute age tie points, the Styx-M core chronology was estimated to cover the last 755 years (1259–2014 CE).

The constrained manual layer counting result was compared with ice chronologies produced by the automated StratiCounter algorithm (Winstrup and others, Reference Winstrup, Svensson, Rasmussen, Winther, Steig and Axelrod2012) and a fifth-degree polynomial age-dating model to verify the robustness of the established chronology. The polynomial model was extrapolated based on fixed points, including the surface age (2015 CE), identified age-depth tie-points, Rittmann tephra layer, and gas chronology of the Styx-M core. The Rittmann tephra, a widely distributed 1252 ± 2 CE tephra layer in Antarctic ice, was identified at a 99.2 m depth (Lee and others, Reference Lee, Kyle, Iverson, Lee and Han2019) and was assigned to 1252 CE in this study (Table 1). Gas ages from 177 points between 68.1 m and 169.8 m, determined by Yang and others (Reference Yang2018), were converted to ice ages by subtracting the ice–gas age difference (Δage) of 320 years, estimated at the depth of the Rittmann tephra layer. The Δage was assumed to be constant with time, with an uncertainty of 20 years (Yang and others, Reference Yang2018). The fifth-degree polynomial was selected because it best matched the known age of the fixed points. Consequently, the manual layer counting-derived chronology agreed well with the age-dating results from the automated and model-based approaches (Figure 5), reinforcing the robustness of the established ice chronology.

Figure 5. The Styx-M core chronologies based on annual layer counting in the depth interval between 9.89 and 98.5 m established in this study (red line) and for the top 9.89 m by Nyamgerel and others (Reference Nyamgerel2024) (light blue line), as well as the StratiCounter program (blue line) and the fifth-degree polynomial age model (black line) (see text). More details on the statistic information applied in the StratiCounter program are given elsewhere (Winstrup and others, Reference Winstrup, Svensson, Rasmussen, Winther, Steig and Axelrod2012). Also shown are the depth intervals of ²³⁹Pu peaks (green squares; Shin and others., Reference Shin, Lee, Han, Hwang, Ro and Hur2025), well-defined signals from major volcanic eruptions (dark blown circles), Rittmann tephra layer (gold diamond; Lee and others, Reference Lee, Kyle, Iverson, Lee and Han2019), and methane tie-points (purple triangles; Yang and others, Reference Yang2018) (see Table 1 and Figure 4). the methane tie-points below 98.5 m are not shown.

The uncertainty of the annual layer counting was estimated by summing individual layer uncertainties, each assigned as 0.5 ± 0.5 years (Rasmussen and others, Reference Rasmussen2006; Nardin and others, Reference Nardin2021), between consecutive tie-points. Layer uncertainties were based on ambiguous layers (i.e., layers with seasonal signals but not confidently assigned as annual) and assigned layers within depth intervals of missing records of impurities, corresponding to the scraped sections of the core during decontamination or samples for which the measurement results could not be obtained because of analytical issues (see section 2.2 and Figure S4). The number of years assigned within these missing sections was estimated using the average ratio between the number of years and the associated depth intervals before and after each gap (Nardin and others, Reference Nardin2021). Table 2 lists calculated age uncertainties, ranging from 1.5 to 16 years (2.2 to 4.4%) for the respective depth intervals. Furthermore, comparisons of the assigned years of the age–depth tie-points (except CH₄ tie-points) with those from the automated layer counting and the polynomial model revealed differences of ± 1 year and − 8 to 4 years, respectively (Table 1). The CH₄ tie-points were excluded from the comparison because Δage, which was assumed to be constant in this study, actually fluctuates with time (Blunier and others, Reference Blunier, Spahni, Barnola, Chappallaz, Loulergue and Schwander2007; Jang and others, Reference Jang2019), likely introducing potential uncertainties in the gas-derived ice ages, as shown by the slight distortions in age between 87.4 m and 97.9 m (see Figure 5). Considering all factors, the maximum dating uncertainty was estimated to be ± 8 years.

Table 2. Uncertainties in annual layer counting for the depth intervals between two consecutive tie-points from the styx-m core (also see Table 1). uncertainties have not been quantified when the impurity records show clear annual cycles between consecutive tie points. No uncertainties exist for dating the top 9.89 m (Nyamgerel and others, Reference Nyamgerel2024)

3.2. Temporal snow accumulation variability

Based on the dating achieved, the mean annual SAR between 1259 CE and 2014 CE was calculated using the methodology described in section 2.3. Table S3 lists the calculated annual SARs. The annual SAR during this period was averaged at 113 ± 39 (1σ) mm w.e. yr^–1 (Table S4), which was slightly lower than a firn densification model-based estimate of 130 mm w.e. yr^–1 (Han and others, Reference Han2015). However, this model-based SAR estimation may include uncertainty, primarily due to its assumption of a constant SAR throughout the simulation period (Han and others, Reference Han2015).

The robustness of the observed results was evaluated by comparing the present data-driven annual SAR record with the annual mean precipitation–minus–evaporation (P–E) budget at the grid point (163.75°E and 73.75°S) nearest to the SGP, calculated using the ERA5 reanalysis from the European Centre for Medium Range Weather Forecasts (ECMWF) for the period between 1940 and 2014 CE (Hersbach and others, Reference Hersbach2020) (Figure S5). The mean SARs derived from the Styx-M core and ERA5 reanalysis were comparable during the corresponding period, with mean values of 117 ± 41 and 109 ± 26 mm w.e. yr^–1, respectively (Table S4). Compared to the less variable temporal fluctuations in the ERA5 P–E budget, however, the Styx-M core SAR record exhibited enhanced values from the mid-1970s to the mid-1980s and during the 1990s (Figure S5). This suggests that the ERA5 model may not properly capture the orographic effects on small-scale SAR patterns in mountainous regions (Daloz and others, Reference Daloz2020; Blau and others, Reference Blau, Kad, Turton and Ha2024). Nevertheless, no significant differences in the overall temporal trends between the different datasets reflect the confidence in the reliability of these results, validating its utility for reconstructing the temporal patterns of the centennial-scale SAR changes at the SGP.

Figure 6a shows the full ice core record of annual SAR and its variability over the past 755 years, from 1259 CE to 2014 CE. Considering the maximum absolute uncertainty (± 8 years) in the depth-age scale, an eight-year running-averaged SAR record was used to examine the temporal SAR variability at the SGP. The temporal variability of the record highlights the dominance of negative SAR anomalies, defined as values lower than the average over the entire period of record, lasting from the mid-13th century to the end of the 16th century. Since ∼ 1600 CE, above-average snow accumulation values were observed, with intermittent negative anomalies, particularly during the 18th century. Despite the existence of inter-decadal and multi-decadal fluctuations, this record revealed a statistically significant increase (∼30%; p < 0.001) in the mean SARs throughout the entire period of 755 years (Figure 6a), with an increase in the average SAR values from 86 ± 29 mm w.e. yr⁻¹ at the beginning of the record (1259‒1300 CE) to 122 ± 32 mm w.e. yr⁻¹ in recent years (2000‒2014 CE) (Table S4). The long-term increasing trend was partly consistent with the ice flow thinning model-based SAR record reconstructed from the Styx-M core, showing a continuous increase of the SAR at the depth from 98.5 (1259 CE) to ∼ 40 m (∼1813 CE) (Yang and others, Reference Yang2018). On the other hand, in contrast to a clear positive trend of the annual SAR over the past two centuries in this record (Figure 6a), the model-based SAR time series revealed a significant decrease in the SAR during the same period (Yang and others, Reference Yang2018). This apparent discrepancy between the two SAR records is probably because Yang and others (Reference Yang2018) calculated the annual SAR with the ice age-depth model based on a linear interpolation of the gas-derived chronology that was determined by matching the CH₄ age tie-points at the 85.7 to 167.1 m depth interval with the gas chronology of the WAIS Divide ice core (WD2014 scale) (Buizert and others, Reference Buizert2015). Such interpolation could introduce significant errors in reconstructing the SAR record from the upper part of the core because the gas-derived SAR estimation does not resolve the high-frequency variations (Winstrup and others, Reference Winstrup2019). This emphasizes the importance of achieving well-constrained chronologies of ice cores using a multiple-proxy approach for reconstructing the long-term variations in snow accumulation at a given site.

Figure 6. Time series of the annual (thin line) and 8-yr running averaged (thick line) snow accumulation rates from (a) the Styx-M (1259‒2014 CE; this study), (b) GV7 (1179‒2009 CE; Nardin and others, Reference Nardin2021), (c) Talos Dome (1217‒1996 CE; Stenni and others, Reference Stenni2002), (d) Hercules Névé (1770‒1992 CE; Stenni and others, Reference Stenni1999), and (e) Roosevelt Island Climate Evolution (RICE) (700 BCE‒2012 CE; Winstrup and others, Reference Winstrup2019) ice cores. Color shadings represent above- (red) and below-average (blue) accumulation rates, respectively. Also shown are the linear regression lines over the entire period (dotted green in a–e) and the period between 1725 and 2014 CE (dotted yellow in a–b).

The overall increasing trend in the Styx-M SAR record for the last 755 years was generally in agreement with the 800-year SAR variability record from the GV7 ice core, recovered from the coast of Oates Land (Figure 1), which showed a consistent increase in the accumulation rate during the corresponding period (Figure 6b) (Nardin and others, Reference Nardin2021). Nevertheless, a different spatial pattern was observed in the trend of the SAR variability in Victoria Land for the Talos Dome (TD) ice core record (Figures 1 and 6c), providing the long-term SAR variability back to ∼ 800 years (Stenni and others, Reference Stenni2002). A long-term accumulation record from this core showed no consistent increase in the SAR over time (Figure 6c). On decadal to centennial timescales over the last 300 years, our SAR variability record exhibited decadal-scale fluctuations with the dominance of above-average accumulation values after the early 1800s and a long-lasting centennial scale increase in the period between 1725 to 2014 CE (Figure 6a). This feature was observed in the SAR variability record from the GV7 core (Nardin and others, Reference Nardin2021). For comparison, the SAR variability record from the Hercules Névé firn core, northern Victoria (Figure 1), spanning at least 200 years, showed mean accumulation rates greater than the average after the mid-19th century (Figure 6d) (Stenni and others, Reference Stenni1999). In contrast, the TD core record was characterized by SAR values larger than the long-term average value during the 20th century (Figure 6c) (Stenni and others, Reference Stenni2002). Compared to the SAR variability elsewhere in the western Ross Sea region, the accumulation rate record of the Roosevelt Island Climate Evolution (RICE) ice core from the eastern Ross Sea (Figure 1), covering the past 2700 years, showed a clear long-term decreasing trend from the mid-13th century to the present with the dominance of below-average accumulation values after the early 1700s (Figure 6e) (Winstrup and others, Reference Winstrup2019). This contradictory feature exhibited a distinct dipole pattern of opposed accumulation variability between the eastern and western Ross Sea region, referred to as the Ross Sea Dipole, which can be attributed to the complex interaction effects of the oceanic and atmospheric impact factors (Scarchilli and others, Reference Scarchilli, Frezzotti and Ruti2011; Bertler and others, Reference Bertler2018; Emanuelsson and others, Reference Emanuelsson2018; Winstrup and others, Reference Winstrup2019).

Overall, the observed differences between the spatial and temporal variations in the SARs core records accentuate the importance of reconstructing long-term snow accumulation patterns from the coast to the interior across the Ross Sea region. Snow accumulation in the Ross Sea region may be related to large-scale atmospheric dynamics such as those induced by the El Nino-Southern Oscillation (ENSO) and the Southern Annular Mode (SAM) as well as synoptic-scale cyclones that can have significant impacts on the extent of sea ice (Markle and others, Reference Markle, Bertler, Sinclair and Sneed2012; Bertler and others, Reference Bertler2018; Emanuelsson and others, Reference Emanuelsson2018; Winstrup and others, Reference Winstrup2019). Future studies will be needed to investigate the dominant factors controlling the decadal- to centennial-variability of snow accumulation and unravel the detailed links of the SAR variability with those for the various components of the ocean, atmosphere, and climate preserved in the Styx-M core.