**Firn cores**

Three shallow firn cores were drilled in the vicinity (in a distance of 1-3 km) of Vostok Station in 2016-2019 with a use of a light mechanical auger. The cores’ lengths were 70, 55 and 65 m. In the glaciological laboratory of Vostok we first defined the density of the firn by thorough measurements of the cores’ volume and mass. The procedure and the results of the density measurements are published in (Ekaykin et al., 2022)17.

Then along the core’s main axis samples for the stable water isotope analysis were cut (see below) with the resolution of 10 cm.

The electrical conductivity (ECM) was measured continuously along each core to define the position of the layers containing the products of the volcanic eruptions (see below). The rest of one of the cores was shipped to Limnological Institute of Russian Academy of Sciences (Irkutsk, Russia) for chemical analyses. The methods and results are published in (Veres et al., 2023)18.

**Surface mass balance**

For the central parts of the East Antarctic plateau the surface mass balance, SMB (a term “net snow accumulation rate” is used as well as a synonymous) is equal to the difference between bulk snow accumulation (precipitation) and ablation (sublimation)19. Based on the firn and ice core data a reconstruction of only the SMB is possible. Usually it is assumed that the SMB is a reliable proxy of the precipitation in central Antarctica, since the total annual sublimation is relatively small: in the case of Vostok, in the present-day conditions sublimation removes roughly 15-20% of the total annual precipitation20.

In the high-accumulation zones (like in central Greenland or in the coastal regions of Antarctica) where annual firn layers can be observed (based on seasonal cycles in the concentration of heavy water isotopes or certain chemical compounds), one can reconstruct annual values of the SMB. In the case of low-accumulation sites like Vostok, only average SMB values between the absolute age markers can be calculated.

Here we use the volcanic age markers identified in all the three cores based on the data on the electrical conductivity and non-marine sulfate concentrations18. In total, 68 volcanic peaks were discovered, of which 22 were attributed to well-dated eruptions. Note that the dates of the firn layers containing the products of the volcanic explosions are 1-2 years younger than the eruptions themselves.

The mean SMB between two adjacent volcanic peaks are calculated based on the available firn density profile17. In particular, in order to define the mean SMB between, e.g., layers of Tambora and Agung, one should divide the relative mass (in g cm-2) of this core section (integrated density profile) by the number of years comprised between these two markers (in this example it equals to 1964-1816 = 148 years).

Two corrections are usually applied to the SMB time series derived from ice core data: 1) for ice advection and 2) for the layer thinning. In this case both of them are negligibly small. With the mean ice velocity of 2 m per year, the oldest firn layers originate about 4 km upstream from the drilling site, where the accumulation rate is the same3. The layer thinning is proportional to *h*/*H*, where *h* is the depth of a layer, and *H* is the total thickness of the glacier, both are in ice equivalent. For the deepest firn layers the correction is about 1 %.

The time-series of SMB is shown in Figure 2. The blue shading depicts the uncertainty bars (±2 standard errors of mean) related to the spatial variability in SMB values. Such spatial variability (‘depositional noise’) is typical for the snow accumulation in central Antarctica3. As an example, the mean SMB between Agung and Pinatubo layers (1964-1992) is 1.89 g cm-2 yr-1 in one core, 1.99 g cm-2 yr-1 in the second one and 2.085 g cm-2 yr-1 in the third core. The mean of the three cores is 1.99 g cm-2 yr-1 with the standard deviation (STD) of 0.10 g cm-2 yr-1 and the standard error of mean (SEM, defined as STD divided by the square root of the number of observations) is 0.06 g cm-2 yr-1.

Note that this ‘depositional noise’ is typical not only for the SMB time-series, but also for the other cores’ properties including the stable water isotopic content (see below).

**Air temperature reconstruction from stable water isotope data**

The measurements of the concentration of heavy water isotopes (d18O and dD) were performed in Climate and Environmental Research Laboratory of Arctic and Antarctic Research Institute with the use of Picarro L2130-*i* and L2140-*i* analyzers. The analytical precision was estimated by re-measurement of 10 % of randomly chosen samples and was equal to 0.05 and 0.5 ‰ for d18O and dD, correspondingly.

As a result of the measurements we obtained the vertical profiles of the d18O, dD and deuterium excess (*dxs* = dD - 8*d18O) for each core with the resolution of 10 cm. These vertical profiles were then transformed into the time-series of the same parameters using the depth-age scale developed in (Veres et al., 2023)18.

Then based on these time-series we calculated the mean isotopic values for the same time intervals for which the SMB data are available (Supplementary Figure 1).

Note a rather wide uncertainty envelopes around the mean values of d18O and *dxs*. This uncertainty is not related to the instrumental error of the isotopic measurements (which is 1-2 orders of magnitude smaller), but with the spatial variability of the isotopic values between different cores due to a high fraction of the ‘depositional noise’ in the total variance of the isotopic values in individual cores.

It is argued that logarithmic formulation of deuterium excess rather than classical linear formulation is preferable when studying the past temperature based on the ice core stable water isotopic data15. However, when the climate variability is relatively small (as for the last 2000 years in central Antarctica) and the range of isotopic values is small as well, then the linear *dxs* produces the same result as the logarithmic one15. Thus, in this study we use the linear *dxs* for simplicity.

The shading around the mean isotopic values depict the uncertainty bars (±2 SEM) related to the spatial variability of the firn isotopic values. Note that the typical values of SEM are 1-2 orders of magnitude larger than the instrumental error of the determination of the d18O and dD values.

To define the past local air temperature from the isotopic data we applied a simplified approach developed by K. Cuffey and F. Vimeux21, 22. The method is based on the assumption that both isotopic content (either d18O or dD) and deuterium excess (*dxs* = dD - 8*d18O) depend on the local condensation temperature (*Tc*) and the near-surface air temperature (*Ts*) in the moisture source:

ΔδD = *a* Δ*T**c* + *b* Δ*T**s*

and

Δ*dxs* = *c* Δ*T**c* + *d* Δ*T**s*,

where Δ denotes the deviation from the present-day values (i.e., from the mean over the 1988–2018 period).

Another parameter affecting δD and *dxs* is the relative humidity in the moisture source, but it is considered to be as a function of *T**s*. It is not possible to reconstruct the relative humidity independently from *T**s* unless a third independent isotopic characteristic (‘17O-excess’) is involved23.

When reconstructing the temperature during cold glacial epochs, the change in the sea water isotopic composition must be taken into account, as well, but for the last 2000 years this factor is negligible24.

From this system of equations it is easy to derive *T**c* and *T**s* as a function of δD and *dxs*:

$$\Delta {T}_{c}=\frac{d\Delta \delta D-b\Delta dxs}{ad-cb}$$

and

$$\Delta {T}_{s}=\frac{a\Delta dxs-c\Delta \delta D}{ad-cb}$$

The values of coefficients *a*, *b*, *c* and *d* are taken from a simple isotope model25: *a* = 10.2‰ °C− 1, *b* = -3.2‰ °C− 1, *c* = -1.55‰ °C− 1 and *d* = 1.6‰ °C− 1. The present-day values (the means for 1988–2018 period) of δ18O, δD and *dxs* are − 57.09‰, -440.05‰ and 16.6‰, correspondingly.

It should be noted that these coefficients are not constant, they are dependent on the temperature15, but for the late Holocene, when the temperature variability was relatively small, they can be taken as constants with sufficient accuracy.

The time-series of the Δ*T**c* is then reduced to the time-series of the near surface air temperature Δ*T**s* applying a scaling factor of 0.67: Δ*T**s* = Δ*T**c* / 0.6715.

One can note that the uncertainty of the temperature reconstructions (Figs. 2b and 2c) are relatively larger than the uncertainty of the isotopic time-series (Supplementary Fig. 1), which is explained by the summation of the uncertainties of δD and *dxs*.

**Likelihood of similar to modern SMB and temperature values in the past**

According to the instrumental measurements3, the mean SMB value at Vostok over the past 52 years (1970-2021) is 2.25±0.13 g cm-2 yr-1. This value is considerably higher than those observed before 1816: the maximum SMB in the pre-industrial epoch (2.1±0.02 g cm-2 yr-1) was observed in 1286-1345. However, one could argue that a similar to present SMB value could take place during, e.g., the 426-year interval between 682 CE and 1108 CE. In order to estimate a possibility of this, we first clustered the available pre-1816 SMB values (Supplementary Table 1) into 4 groups according to the length of the time intervals: 1) 20-40 years (4 values, average length is 32 years), 2) 40-60 years (3 values, average length is 51 years), 3) 90-170 years (8 values, average length is 120 years) and 4) 300-500 years (2 values, average length is 382 years). Then for each group we calculated standard deviation of the SMB values. As expected, the variability of the mean snow accumulation rates decreases when the period of observation increases (Supplementary Figure 2).

The climatic variability of the SMB values averaged over a 50-year time interval is equal to 0.19 g cm− 2 yr− 1 (1 STD). The mean SMB value in 682–1108 was 1.83 g cm− 2 yr− 1. A maximum 50-year SMB value within this 426-year time interval can be estimated as the mean + 2 STD and equals to 2.21 g cm− 2 yr− 1, which is less than 2.25 g cm− 2 yr− 1, observed in the last 50 years. The likelihood that a similar to present-day SMB value happened in 682–1108 is 1.3%. There is even less likelihood that mean 50-yr SMB value reached 2.25 g cm− 2 yr− 1 during the 337-year time interval between 168 BCE and 169 CE.

Similarly, it can be shown that the mean SMB value in 1816–1964 (2.06 g cm− 2 yr− 1) is unprecedented for the previous 2200 years. The likelihood than in any 150-yr interval before 1816 the mean SMB reached 2.06 g cm− 2 yr− 1 is < 0.1%.

As for the air temperature, although the present-day values are significantly higher than in the pre-1816 period, we cannot conclude that they are unprecedented for the past 2200 years because of large uncertainties (Fig. 2). For example, in 1641–1695 the mean Δ*T**g* was − 0.66°C, but the upper limit of the error bars is + 0.13°C, i.e., warmer than the present-day value.

**Sensitivity of the SMB to air temperature**

It is assumed that the sensitivity of the Antarctic surface mass balance (SMB) to the air temperature in central Antarctica is proportional to the Clausius-Clapeyron relationship, which describes the saturation water vapor pressure *e*s as a function of temperature *T* (Nicola et al., 2023)2:

$$\frac{dln{e}_{s}}{dT}=\frac{L}{{R}_{v}{T}^{2}}\equiv \alpha \left(T\right)$$

,

where *L* is the latent heat of vaporization (2.5 × 106 J kg− 1), *R**v* is the specific gas constant for water vapor (461 J K− 1 kg− 1) and *T* is in K.

Here *T* is the effective condensation temperature *T**c*, which is in the case of the central Antarctic is usually approximated by the air temperature at the top of the inversion layer *T**inv* 14. At Vostok the mean value of *T**inv* for the period of 1963–1991 (when the balloon-sounding data is available) is -38.1°C (= 235.05 K)16, and the corresponding SMB-temperature sensitivity *α* is 9.8% per 1°C. However, it is argued that so-called ‘diamond dust’ (a form of precipitation typical for central Antarctica) is formed not on the top of the inversion, but rather in the whole inversion layer, as soon as air is cooling during its descent towards the surface. In this case the effective condensation temperature would be much colder, of the order of -49°C (224.15 K), and the corresponding sensitivity is 10.8% per 1°C.

Then, if we want to calculate the sensitivity of the SMB to the near-surface air temperature, we need to consider different amplitude of temperature variability in the inversion layer and in the near-surface layer. It is assumed that the slope between temporal variability of *T**g* and *T**c* is 0.67 14, 15, i.e., the amplitude of condensation temperature is 0.67 times that of the near-surface temperature. Applying this coefficient, we obtain the slope between the SMB and *T**g* equal to about 16% per 1°C.

Thus, a relatively strong SMB-temperature sensitivity obtained for Vostok can be easily explained by a very cold temperatures typical for this region.

When considering the relationship between the SMB and temperature, one should remember that the Clausius-Clapeyron equation describes the sensitivity of precipitation to temperature, rather than of the SMB. As mentioned above (see ‘Surface mass balance’ in the Methods) the SMB in central Antarctica is the difference of precipitation and sublimation. The sublimation rate is also related to temperature20: the warmer is air, the stronger is sublimation. Thus the total sensitivity of the SMB to temperature may be slightly weaker than that predicted by the Clausius-Clapeyron relationship.