Reconstructing atmospheric circulation and sea-ice extent in the West Antarctic over the past 200 years using data assimilation

The West Antarctic climate has witnessed large changes during the second half of the twentieth century including a strong and widespread continental warming, important regional changes in sea-ice extent and snow accumulation, as well as a major mass loss from the melting of some ice shelves. However, the potential links between those observed changes are still unclear and instrumental data do not allow determination of whether they are part of a long-term evolution or specific to the recent decades. In this study, we analyze the climate variability of the past two centuries in the West Antarctic sector by reconstructing the key atmospheric variables (atmospheric circulation, near-surface air temperature and snow accumulation) as well as the sea-ice extent at the annual timescale using a data assimilation approach. To this end, information from Antarctic ice core records (snow accumulation and δ18O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta ^{{18}}\hbox{O}$$\end{document}) and tree-ring width records situated in the mid-latitudes of the Southern Hemisphere are combined with the physics of climate models using a data assimilation method. This ultimately provides a complete spatial reconstruction over the West Antarctic region. Our reconstruction reproduces well the main characteristics of the observed changes over the instrumental period. We show that the observed sea-ice reduction in the Bellingshausen-Amundsen Sea sector over the satellite era is part of a long-term trend, starting at around 1850 CE, while the sea-ice expansion in the Ross Sea sector has only started around 1950 CE. Furthermore, according to our reconstruction, the Amundsen Sea Low pressure (ASL) displays no significant linear trend in its strength or position over 1850–1950 CE but becomes stronger and shifts eastward afterwards. The year-to-year sea-ice variations in the Ross Sea sector are strongly related to the ASL variability over the past two centuries, including the recent trends. By contrast, the link between ASL and sea-ice in the Bellingshausen-Amundsen Sea sector changes with time, being stronger in recent decades than before. Our reconstruction also suggests that the continental response to the variability of the ASL may not be stationary over time, being significantly affected by modification of the mean atmospheric circulation. Finally, we show that the widespread warming since 1958 CE in West Antarctica is unusual in the context of past 200 years and is explained by both the deeper ASL and the positive phase of the Southern Annular Mode.


Introduction
The Antarctic Ice Sheet is the biggest reservoir of fresh water on Earth that would potentially raise the global sealevel by 58m if the entire ice sheet melted (Shepherd et al. 2018). Despite the remote location, the Antarctic and the Southern Ocean have experienced major climate changes over the past decades (e.g., Jones et al. 2019;Bromwich et al. 2013;Medley and Thomas 2019;Pritchard et al. 2012), demonstrating large variability and their vulnerability to the global climate change.
The largest changes have been found over the Antarctic Peninsula and the West Antarctic Ice Sheet (hereafter 1 3 WAIS) (e.g., Jones et al. 2019), together forming West Antarctica. A major warming has been observed since the International Geophysical Year (i.e., 1958 CE), far exceeding the global warming during the same period (Vaughan et al. 2003;Turner et al. 2005a;Nicolas and Bromwich 2014). Although the warming is particularly pronounced over the Peninsula, the warming over the central WAIS since the 1950s is one of the fastest recorded on Earth (Steig et al. 2009;Bromwich et al. 2013). Moreover, in the last decades, glaciers from West Antarctica (especially the Thwaites and Pine Island glaciers) have been losing mass at the ocean-ice shelf (i.e., the floating part of the glacier) interface at an accelerated rate (Pritchard et al. 2012;Rignot et al. 2019). This ice shelf melting positively contributes to the global sea-level rise, by directly enhancing the ice flow to the Southern Ocean. In contrast with these spatially homogeneous changes, snow accumulation integrated over West Antarctica has increased at an increased pace over the twentieth century, but with strong regional differences (Medley and Thomas 2019;Wang et al. 2019). While the Antarctic Peninsula and the Eastern WAIS have gained mass at the surface through snow accumulation, the Western WAIS has displayed a snow accumulation decrease (Medley and Thomas 2019;Wang et al. 2019). Additionally, a large reduction of the sea-ice extent has been noticed in the Bellingshausen/Amundsen Sea sector that contrasts with an expansion in the Ross Sea sector over the satellite era (starting from 1979 CE) (Parkinson 2019).
The observed changes in West Antarctica are not independent of each other, and have been widely attributed to changes in atmospheric circulation (e.g., Bromwich et al. 2013; Thomas et al. 2013;Steig et al. 2013), characterized by an increase in the intensity and a poleward shift of the westerly winds (Westerlies) (Marshall 2003). This has been observed as a deepening in the Amundsen Sea Low (ASL), a quasi-stationary low-pressure system located off the Amundsen coast (e.g., Raphael et al. 2016). Due to its location, variations in both the strength and position of the ASL strongly modulate the West Antarctic climate (e.g., Hosking et al. 2013). For instance, a stronger ASL enhances the northerly flow (onshore winds) to the continent, which warms the Antarctic Peninsula and the eastern WAIS (e.g., Hosking et al. 2013) and increases snow accumulation (Turner et al. 2005b; Thomas et al. 2008). At the same time a stronger ASL enhances the southerly flow (offshore winds) of cold dry air masses over the western WAIS leading to lower temperatures and reduced snow accumulation. In addition to continental changes, a stronger ASL is also associated with reduced sea-ice cover in the Bellingshausen/Amundsen Sea sector and a sea-ice cover expansion in the Ross Sea sector (e.g., Raphael and Hobbs 2014). A deeper ASL also drives oceanic changes, enhancing warm ocean water upwelling towards ice shelves, inducing basal melting (e.g., Mankoff et al. 2012;Dotto et al. 2020).
The ASL is thus a major feature of the atmospheric circulation in the Southern Hemisphere, which displays one of the greatest year-to-year variations on Earth (e.g., Turner et al. 2019). The ASL is strongly modulated by the atmospheric circulation of the high latitudes of the Southern Hemisphere. Particularly, the Southern Annular Mode (hereafter SAM) describes the latitudinal movement of the Westerlies as well as their strength (e.g., Fogt and Marshall 2020). When the SAM is in its positive phase, the Westerlies are stronger and located further south than in the negative phase, which in turn deepens the ASL (e.g., Raphael et al. 2016). The ASL is thus strongly related to the SAM. Wind changes in the West Antarctic are also closely linked to the climate variability in tropical oceanic basins. Several studies (e.g., Ding et al. 2011;Steig et al. 2013;Thomas et al. 2013;Meehl et al. 2016;Holland et al. 2019) have suggested the important tropical teleconnections that exist between the mid-latitudes/ tropics and the West Antarctic. More specifically, positive anomalies in sea surface temperature in the tropical Pacific Ocean induce a convective heat event that further propagates to the West Antarctic region via the formation of Rossby wave train.
Despite the large changes observed recently, the impact of the ASL variability on West Antarctic climate changes has not been fully studied because in-situ observations are lacking (Turner et al. 2004). Over the satellite era (from 1979 CE), atmospheric reanalyses seem to be reliable for studying the ASL variability and its impacts on the West Antarctic climate (Bracegirdle 2013). However, because of the high internal variability in the Amundsen region (e.g., Raphael et al. 2016), longer time-series are needed to understand the processes controlling the variability in the region.
In this study, we aim to better understand the West Antarctic climate variability on decadal to centennial timescales by analyzing the relationships between the main atmospheric variables with surface conditions and sea-ice cover. Furthermore, by making use of those potential links between variables, we provide historical estimations over the past two centuries that could not be directly observed. Special attention is given to the ASL, because of its major importance in the ongoing observed changes in West Antarctica. We investigate to which extent the ASL variability explains the West Antarctic climate change over the past two centuries, and if the role of the ASL has changed over time. The sea-ice extent (in particular the Bellingshausen/Amundsen Sea and Ross Sea sector) is also examined due to its key role in modulating West Antarctic climate (e.g., Lefebvre and Goosse 2008;Thomas et al. 2013;Turner et al. 2017). Up to now, no spatial estimate of wind fields or sea-ice extent has been specifically validated for this region before the instrumental era and thus we aim to fill this gap. In addition to improving the general long-term changes in the West Antarctic climate, we assess the representativeness of the climate changes occurring during the second half of the twentieth century on a longer time perspective. Finally, a complementary goal is to highlight the mechanisms that explain the warming over West Antarctica as a whole and the strong regional asymmetry in sea-ice extent trends over the past decades.
To this end, we provide a spatially complete and multivariate reconstruction of the West Antarctic climate changes over the past two centuries (the period for which most paleoclimate proxies are available; see below) by optimally combining paleoclimate records and a climate model using a data assimilation approach. Over last years, data assimilation has been increasingly applied in paleoclimatology for estimating past climate changes (e.g., Goosse et al. 2009Goosse et al. , 2010Hakim et al. 2016;Steiger et al. 2017). Compared to simpler statistical methods (e.g., Stenni et al. 2017), data assimilation does not assume a stationary relationship between the proxy and the climate. Data assimilation ensures that the resulting climate reconstruction is dynamically coherent. Such a method is also particularly relevant for regions like the West Antarctic where a strong coupling exists between the observed changes from various variables Fan et al. 2014;Goosse et al. 2009). Here, the model is constrained by 18 O and snow accumulation records from Antarctic ice cores. Over the Antarctic, ice core is the most widespread climate archive for assessing historical climate changes at high temporal resolution (i.e., annual scale). Compared to the instrumental network, the ice core network is well developed over West Antarctica thanks to a number of international drilling efforts. In addition to the ice cores, we also make use of the tree-ring width records from the mid-latitudes that capture changes in past climate variability (e.g., Emile-Geay et al. 2017). The inclusion of non-specific-Antarctic proxies guarantees a large-scale coherence. This avoids forcing a local agreement with the Antarctic data that would not be based on processes compatible with largerscale changes.
This study is organized as follows. Section 2 includes a description of the data assimilation method, the paleoclimate proxies and the model used for the reconstruction. We discuss the performance of our reconstruction in Sect. 3 by comparing it with instrumental observations and other climate reconstructions, and discuss our main results, before concluding.

Data assimilation method
Data assimilation (DA) optimally combines the information from observations and the climate physics as included in climate models. In practice, the data assimilation process updates the initial state of the climate given by the model (i.e., the prior) according to the available observations to provide the best estimate of the climate state (i.e., the posterior), while also considering the errors associated with both the data and model. Paleo DA methods thus spread the local and temporal information from proxy records in space but also into other variables by relying on the modelled covariance in space and among variables, respectively. For instance, the 500-hPa geopotential height can be therefore reconstructed only by assimilating near-surface air temperature. In that situation, the reconstructed 500-hPa geopotential height relies on the covariance relationships between the 500-hPa geopotential height and near-surface air temperature in the model. If there is no covariance between the variable of interest and the assimilated variable, the reconstruction skill of the variable of interest will be null. Consequently, the resulting multi-variate climate reconstruction guarantees a dynamical consistency in space and between variables, which is provided in a natural way with the climate model as it directly relies on the physics included in the climate model.
The DA method employed in this study is an offline approach (also called a no-cycling method). This method uses existing climate model simulations to draw the prior. In contrast to the standard online data assimilation approach, the prior is estimated by selecting different years of a climate model simulation. Therefore, no information is propagated in time. This approach is totally appropriate when the predictability time-scale of the system is much smaller than the data assimilation time step (here, one year). This is particularly true when assimilating atmospheric variables because of the small correlation time in the atmosphere. In such cases, no additional relevant information is gained when using an online approach, as shown by Matsikaris et al. (2015).

Particle filter
The DA method used here is based on a particle filter following the implementation of Dubinkina et al. (2011). This method has been recently applied to reconstruct successfully large-scale near-surface air temperature and snow accumulation over Antarctica during the past millennia (Klein et al. 2019;Dalaiden et al. 2020a). We thus describe briefly the method here.
The particle filter (van Leeuwen 2009) aims at updating an initial estimate of the state of the system, referred to as the prior, using additional information provided by available observations, following the classical approach of Bayesian reconstructions. Specifically, the prior for each variable is represented by a probability density function (pdf), described in a discrete way using a set of independent model states, which are called particles. These model states are given by all the years (annual mean) of an ensemble of three simulations performed with the isotope-enabled Community Earth System Model version 1 (iCESM1; Brady et al. 2019;Stevenson et al. 2019) spanning the 850-1850 CE period (see Sect. 2.5 for the model evaluation). The prior thus includes 3003 particles and is constant throughout the data assimilation process. At each assimilation step (here every year), the particles forming the prior are evaluated against the available observations while taking into account the error associated with the observations as well as the inconsistency between the observations and model related to unresolved processes in the model (i.e., the observation error; see Sect. 2.6). As a result of this comparison with observations, each particle of the prior receives a weight that is proportional to its likelihood knowing the observations. The weights are thus higher for the particles closer to the observations than for the particles further to the observations. In other words, the particle filter redistributes the weights of all the particles at each time step to obtain posterior distribution in better agreement with observations. Since only the weights of the particles can change during the data assimilation process, the climate dynamics as represented in the model is fully respected. For each time step, the climate reconstruction is defined as the weighted mean of all the particles. We define the reconstruction uncertainty as the weighted standard deviation of all the particles. For more details on the implementation of the particle filter, see Dubinkina et al. (2011).

Evaluation metrics
In order to quantify the performance of our reconstruction, we use three different metrics that compare the results with observations or other reconstructions. The first one is the Pearson correlation coefficient (r), which determines the strength of the linear relationship between two time-series: where n is the number of samples, y stands for the predicted vector while x is the true vector. Overbars indicate the mean over the time. y and x are the standard deviations of y and x, respectively. Although a high value of r indicates that the predicted values follow the direction of the true values, it does not ensure that the variations of the amplitude are well reproduced. Therefore, the second metric is the coefficient of efficiency (Nash and Sutcliffe 1970), which depends on the amplitude of the signal: A high value for CE thus indicates that both the timing and the amplitude are right. Finally, we need to verify if the ensemble reconstruction is similar to the reconstruction error relative to observations. For this purpose, we compute the ensemble calibration ratio (ECR) defined here as: where i varies from 0 to n and 2 x,i is the variance of the reconstruction ensemble. An ECR of 1 means that the reconstruction is well calibrated (i.e., the reconstruction error equals the uncertainty of the reconstruction). When ECR<1, the estimation of the uncertainty of the reconstruction is too large compared to the error over the analyzed period. On the contrary, when ECR>1, the estimate of the uncertainty is likely too low compared to the real reconstruction error.

Proxy data
From the 79 ice core snow accumulation records included in the database of Thomas et al. (2017), we only keep the annually-resolved records (48 out of 79) as in Medley and Thomas (2019). Medley and Thomas (2019) provide a spatially complete annual Antarctic snow accumulation over the past two centuries using this ice core snow accumulation database (Fig. 1). This excludes a majority of snow accumulation records in the continental Dronning Maud Land region. As in Medley and Thomas (2019), we also add the B40 ice core record of snow accumulation , which was published after the compilation of Thomas et al. (2017). Additionally, we also employ the information from the ratio of stable isotopes oxygen in the ice core ( 18 O ). Stenni et al. (2017) provide a compilation of the Antarctic precipitation-weighted 18 O records (n=112) covering the last two millennia. We only keep the 29 annuallyresolved 18 O records. It is worth noting that the majority of the ice cores (both snow accumulation and 18 O records) are located in West Antarctica (WAIS; Fig. 1).
In addition to the Antarctic proxies, we also utilize the continental annually-resolved proxies situated in the Southern Hemisphere from the global database of the PAGES2k working group that includes all the proxy records suitable for reconstructing the global temperature over the last two millennia (PAGES2k Consortium et al. 2013;Emile-Geay et al. 2017). Among all the continental proxies included in the global database for the Southern Hemisphere, only the tree-ring width (hereafter TRW) are annually-resolved (in addition to some Antarctic ice cores already included in our study). A total of 18 sites presenting TRW time-series over past centuries are located in the Southern Hemisphere : three in New Zealand, three in Tasmania and 12 in South America. Among the 12 South American sites, only six are x,i positively correlated with climate (either near-surface air temperature or precipitation) according to our dataset (see Sect. 2.4) over the 1941-1990 CE period, and thus included. The geographical distribution of the selected sites is presented in Fig. 1.
Most proxy records cover the last two centuries, but their number starts to decline from around 1990 CE (Fig.  S1). Therefore, the period analyzed throughout the study is 1800-2000 CE. This provides a good compromise in terms of the number of records and the overlap with the satellite period (starting in 1979 CE) that is used for the evaluation of the reconstruction (see Sect. 3.1). This period was also used in the study of Medley and Thomas (2019) who reconstructed snow accumulation in Antarctica using the same ice core snow accumulation records as here.
In order to reduce the non-climatic noise and provide equal weighting in each region (i.e not bias a particular region with multiple records), a 500 km grid (square cells of 500 km side) was established and records from the same grid were averaged together. For snow accumulation records, they are all normalized over the 1941-1990 CE period (mean zero and unit standard deviation) and the records situated in the same cell are averaged. The mean of normalized snow accumulation records for each grid cell produces the normalized composites. For each composite, the variance is then scaled to the variance of the spatial snow accumulation Fig. 1 Annually-resolved climate/proxy records used in this study: snow accumulation records (blue dots), 18 O records (purple squares), Automatic Weather Stations (AWSs; green stars) and tree-ring width records (green dots). The map background represents the surface elevation (m) reconstruction from Medley and Thomas (2019), linearly interpolated on the 500 km grid over the 1941-1990 CE period. Compared to snow accumulation, the 18 O records display a weaker spatial variability as they are less dependent on the topography and do not need normalization. Instead, 18 O composites for each grid cell are obtained by averaging the anomalies of 18 O records over the  CE. Finally, we apply the same methodology as for the snow accumulation records to TRW time-series, but without correcting the variance (z-scored composites) since TRW data is originally normalized.

Proxy System Models
Snow accumulation records are directly compared to the precipitation minus sublimation/evaporation ( P − E ) from the climate model simulation. Snow accumulation recorded in ice cores is the sum of both precipitation and post depositional changes (wind erosion, sublimation and melt). Several studies have shown that precipitation dominates and that snow accumulation can be directly compared with the precipitation minus sublimation/evaporation ( P − E ) from the climate model simulation Souverijns et al. 2018;van Wessem et al. 2018;Agosta et al. 2019). This is particularly true when working at a large spatial scale -as in our study (500 km grid)-and over the grounded ice sheet (i.e., the Antarctic Ice Sheet without ice shelves) (Agosta et al. 2019). Additionally, since iCESM1 explicitly simulates the ratios of water isotopes in all the climate components, the 18 O records are also directly compared to the 18 O simulated by the model. A Proxy System Model (PSM) specifically designed for the TRW is required in order to compare the model results to the proxy. To this end, a TRW PSM is built for each TRW time-series. Instead of a mechanistic model, which simulates tree-ring growth by explicitly including the biological processes that drives the relationship between climate and tree growth (e.g., Guiot et al. 2014;Misson 2004;Dufrêne et al. 2005;Rezsöhazy et al. 2020), we use a simple statistical model. It has the advantage of being easily implemented but the relationship between climate and tree growth is estimated empirically, without taking into account any biological processes. This method is similar to the one used in the Last Millennium Reanalysis (Tardif et al. 2019).
In this study, we consider that trees are sensitive to the temperature or precipitation, or to both. Therefore, we introduce two types of models: where y i corresponds to the observed z-scored TRW timeseries (i.e., the dependent variable) for site i; 1i and 2i are the slopes associated with the X 1i and X 2i of the i timeseries, which are the explanatory variables; i are the errors, assumed to be normally (Gaussian) distributed with a zero mean and a unitary variance (0, 2 ). Parameters of the Eqs.
(4) and (5) are estimated using the ordinary least squares method.
Similarly to the Last Millennium Reanalysis (Tardif et al. 2019), we use the near-surface air temperature and precipitation as explanatory variables ( X 1 and X 2 ). In addition to the annual mean of near-surface air temperature and precipitation, the climate variables are also averaged over different months during the calendar year (January to December), in order to take into account the seasonal response of trees: JJA and OND. All the possible combinations are tested. The calibration of the PSM is performed over the 1941-1990 CE period using the Global Meteorological Forcing Dataset for land surface modeling (v2; http:// hydro logy. princ eton. edu/ data. php, last access: 22 June 2020; hereafter GMF) (Sheffield et al. 2006) interpolated on the 500 km grid excluding the oceanic grid cells. The sensitivity to the climate data has been assessed by also performing the calibration with the surface temperature from the NASA Goddard Institute for Space Studies Surface Temperature Analysis (Hansen et al. 2010) and precipitation from the Global Precipitation Climatology Centre . Results display no major difference (not shown).
To select the best PSM for each TRW composite, we rely on the Bayesian information criterion (BIC) (Schwarz 1978) value defined as : where LL is the natural logarithm of the likelihood for the model -estimated as the mean squared error of the linear regression model (Watkins and Mardia 1992)-and k corresponds to the number of parameters in the regression model (i.e., two and three for the uni and bi-variate models, respectively). The BIC is a particularly relevant metric to select the best model among uni and bi-variate models as the most complex models are penalized. Accordingly, the PSM displaying the lowest BIC is selected.

Climate model simulation
The initial estimate of the state of the climate system (i.e., the prior) is derived from an ensemble of three simulations performed with the isotope-enabled iCESM1 spanning the 850-1850 CE time period Stevenson et al. 2019). iCESM1 is a coupled atmosphere-ocean model including a sea-ice component. The atmosphere is resolved at approximately 2° and the ocean at 1°. The iCESM1 (6) BIC = −2 * LL + log(n) * k simulations include the orbital changes, the solar variability, the volcanic forcing through changes in stratospheric aerosols and greenhouse gases, the land use changes and finally the human-induced greenhouse gases .
Although numerous studies have shown that CESM is well-suitable for studying the Antarctic climate over the past, present and future (Lenaerts et al. 2016Fyke et al. 2017;England et al. 2016;Dalaiden et al. 2020b), this version has not yet been evaluated over the Antarctic continent. As the information from proxy records is spread spatially and into other variables using the climate model, it is important to assess the performance of the model used in the data assimilation in simulating the near-surface climate over the Antarctic. Therefore, we briefly evaluate the Antarctic climate as simulated by iCESM1 by comparing it to the latest ECMWF's atmospheric reanalysis, ERA5 (Hersbach et al. 2020). It is considered as one of the best reanalyses in simulating the Antarctic climate over the satellite-era (Gossart et al. 2019;Tetzner et al. 2019). Since the paleoclimate proxies used in this study are annually-resolved, the evaluation is carried out on annual averages.
Over the 1979-2005 CE period, the geopotential height at 500-hPa is well simulated in iCESM1 ( Fig. S2; below south of 45°S, R 2 = 97%; relative bias = -0.2% computed as the mean relative difference between the 500-hPa geopotential height from iCESM1 and ERA5 south of 45°S) and includes the three main low-pressure systems in Amundsen Sea, Dronning Maud Land and Wilkes Land. Regarding the nearsurface air temperature, iCESM1 simulates the temperature gradient with the highest temperatures along the coasts and the lowest temperatures in the interior of the ice sheet (over the Antarctic continent, R 2 = 76%; relative bias = -2.0%). Like temperature, snow accumulation is highly dependent on the topography. iCESM1 reproduces the Antarctic snow accumulation pattern well from ERA5 ( R 2 = 85% -based on log values as the snow accumulation distribution is log-normal (Agosta et al. 2019)-; relative bias = -6.9%). However, iCESM1 is not able to reproduce the high spatial variability of snow accumulation at local scale, because of its coarse horizontal spatial resolution. Besides, iCESM1 displays a total Antarctic sea-ice extent of 11.7 10 6 km 2 against 11.9 10 6 km 2 for observations from the National Snow and Ice Data Center (NSIDC; data available here: https:// nsidc. org/ data/ NSIDC-0051/ versi ons/1, last access: 5 September 2018). Compared with the NSIDC observations, iCESM1 overestimates by 5% the mean sea-ice extent in the West Antarctic sector (160°W-60°E) over 1979-2005 CE, with a similar bias in both the Bellingshausen/Amundsen Sea (130°W-60°E) and Ross Sea (160°E-130°W) sectors (5.4% and 4.8%, respectively; Tab. S1). This suggests that iCESM1 is able to simulate relatively well the mean state of Antarctic sea-ice extent at present. Finally, although a quantitative evaluation of the simulated precipitation-weighted 18 O is not possible because of the insufficient amount of Antarctic observations (three records in the Global Network of Isotopes in Precipitation), the precipitation-weighted 18 O pattern over Antarctica captures the gradient between the coasts (with the highest values, less negative) and the Plateau (with the lowest values, more negative) due to isotopic fractionations (Fig. S3). Based on the skill of iCESM1 in simulating present-day Antarctic climate, we deem this model suitable for building the prior of the data assimilation experiment. As proxy data are averaged over a 500 km regular grid (see Sect. 2.3), the prior has also been linearly interpolated onto this 500 km regular grid.

Observation error
Data assimilation requires an error of the assimilated data. This observation error plays a crucial role in data assimilation because it determines the extent to which each assimilated record is influencing the reconstruction. The records with a low observation error will thus have a larger weight in the data assimilation than the records with a high observation error.
Three types of observation errors are usually mentioned (e.g., Badgeley et al. 2020). The first type of observation error directly comes from the accuracy of the measurement, i.e., the instrumental error. The second type of observation error is related to the processes in observations that are unresolved at the spatial scale of the model because of its coarse resolution, the so-called representativeness observation error (e.g., Oke and Sakov 2008;Janjić et al. 2018). Finally, the last type of observation error arises from the performance of the PSM in simulating the relationship between the proxy and the climate variable(s). Overall, the instrumental error is much smaller than the two other observation errors (e.g., Oke and Sakov 2008), especially in paleoclimatology (Tardif et al. 2019;Steiger et al. 2018;Badgeley et al. 2020). This observation error is thus ignored in our study. Finally, regarding the snow accumulation and 18 O records, we assume that the observation error related to the PSM is non-existent, since the model simulates both variables. The representativeness error is thus considered as the largest contributor to the observation error. Therefore, for ice core snow accumulation and 18 O records, the observation error is taken equal to the representativeness error. However, for tree-ring width records, the observation error combines the representativeness error and the PSM-related observation error as the climate model used in the data assimilation does not simulate tree growth. In contrast, the climate model simulates the snow accumulation and 18 O . As a consequence, a PSM simulating those variables is not required.
To estimate the representativeness error of the snow accumulation composite associated with the ice core records, we use an Antarctic simulation performed with the latest version of the Regional Atmospheric Climate MOdel (RACMO2.3p2., hereafter RACMO) at 27 km horizontal resolution (van Wessem et al. 2018). RACMO is a polaroriented regional climate model that specifically resolves near-surface processes over Antarctica. For each snow accumulation composite, we proceed as follows. We retrieve the annual time-series of snow accumulation from RACMO at the location of all the ice core snow accumulation records included in the composite over the 1979-2016 CE period (period over which RACMO results are available). In order to remove the dependence on the elevation, the calculation is performed on time-series anomalies. All RACMOtime series are then averaged in time. Next, we compute the difference between this averaged snow accumulation time-series and the snow accumulation time-series from RACMO interpolated on the 500 km grid (i.e., the resolution of snow accumulation composites; see Sect. 2.3). Eventually, the representativeness error is defined as the standard deviation of the time-series difference. Note that the representativeness error is calculated for each year of the composite (1800-2000 CE), based on the snow accumulation RACMO data (1979. As the number of ice cores in the composite decreases when going further back in time, the representativeness error tends to increase. Finally, since RACMO does not simulate the 18 O , we use the iCESM1 linearly interpolated on the RACMO grid over the 1950-2005 CE period to estimate the 18 O representativeness error with the same methodology. A 50-year period minimizes the potential biases induced by the variability of the 18 O when using a too short time period and thus ensures the computation of a robust estimate. We have selected the final part of the series to have the largest overlap with the period over which the snow accumulation error is computed. This in part guarantees consistency in the calculation of the representativeness error for the snow accumulation and 18 O records. However, several studies (van Wessem et al. 2018;Lenaerts et al. 2016;Cavitte et al. 2020) have shown that the spatial variability of snow accumulation and 18 O is underestimated in models (i.e., RACMO and iCESM1) compared to ice core observations. Therefore, using these models to estimate the observation error in our study probably also underestimates the observation errors. In order to provide an observation error closer to the reality, we take into account the missing processes occurring between the scale represented in the models and the local scale of ice cores. This is achieved by performing a side-by-side comparison between adjacent snow accumulation records as in Fisher et al. (1985). More specifically, we compute the standard deviation of the difference between each pair of annual time-series of snow accumulation records (in anomalies) within the same grid cell of the 500 km grid over the 1941-1990 CE (which is the most recent 50-year period for which all the ice core records are available; Fig. S1). The calculation is made for all the possible combinations within the grid cell. According to the grid cell containing more than five ice cores (n = 2; both situated in West Antarctica), the side-byside comparison suggests that the error at the ice core-scale is higher by a three-factor compared to the one at the scale of RACMO. The same exercise for 18 O has been carried out and shows similar results (not shown).
Regarding the TRW composites, the error is estimated by taking the standard deviation of the residuals resulting from the linear regression (see Eqs. (4) and (5)) as in the Last Millennium Reanalysis (LMR) (Hakim et al. 2016;Tardif et al. 2019). The error thus encompasses the representativeness error and the error related to the PSM since the PSM has been calibrated using the climate data interpolated on the 500 km grid and not using the local data. After performing several sensitivity tests, we have found that our results present a small sensitivity to reasonable modification of the estimation of the observation error (not shown).

Instrumental data assimilation-based reconstruction
Before using the paleo records to reconstruct the Antarctic climate over the past centuries (hereafter referred to as the paleo-reconstruction), we first reconstruct the Antarctic climate over the 1958-2000 CE period using data assimilation with near-surface air temperature records and snow accumulation records from ice cores. Applying a methodology that is as close as possible to the one selected for the longer timescales allows us to validate the approach, in particular the data assimilation method and the estimations of the errors using observations with lower uncertainties than the paleo records. This reconstruction based on instrumental and snow accumulation records thus provides an upper bound for the skill we could expect with the paleoclimate network.
To this end, we assimilate the near-surface air temperatures from the Automatic Weather Stations (AWSs) over the Antarctic continent (Turner et al. 2004) (n = 14) starting from 1958 CE (see Fig. 1 for their geographical localization). The model is also constrained by the snow accumulation records from ice cores, which are the best estimate of long-term snow accumulation, since AWSs do not record this variable. Additionally, TRW records are replaced by the near-surface air temperatures from the gridded dataset used for the calibration of the TRW PSM (i.e., GMF; only where TRW sites are available). This reconstruction is referred to as the instrumental reconstruction. The observation error on the near-surface air temperature from AWSs is calculated similarly to the snow accumulation error (see Sect. 2.6) but using near-surface air temperatures from RACMO. Like for the snow accumulation and 18 O records, we assume that the observation error for the Antarctic near-surface air temperature is mainly due to the representativeness error. This error should take into account the missing processes occurring between the local scale (i.e., the AWS scale) and the scale of RACMO (i.e., 27-km resolution). Nevertheless, the network of AWS starting from 1958 CE is much less dense than the paleo records (Fig. 1). Therefore, a side-by-side comparison of adjacent near-surface air temperature time-series as we did for paleo records is not possible (see Sect. 2.6). We are thus not able to provide an accurate estimate of the factor to make the link between local observation error and the error at the scale on the climate model grid. However, several studies (e.g., Agosta et al. 2019;van Wessem et al. 2018) showed a weaker spatial variability of near-surface air temperature compared to snow accumulation. Therefore, we assume that the factor applied to take into account the unresolved processes at the RACMO scale is smaller than for snow accumulation (see Sect. 2.6) and multiply the local observation error of near-surface air temperature for the Antarctic AWSs by a two-factor. Finally, the situation is very different for the assimilated near-surface air temperature of the gridded dataset GMF associated with the TRW sites as it is a gridded dataset at a resolution similar to the one of the model that is thus expected to represent processes at a similar scale. Therefore, a constant error of 0.5 °C is chosen, which is a typical value when assimilating gridded near-surface air temperature (e.g., Brennan et al. 2020;Dubinkina et al. 2011).

Validation of the instrumental paleo reconstructions over the last decades
Before analyzing the long-term changes using the paleoreconstruction over the past centuries, we first need to ensure that our data assimilation method works. We establish this by testing the skill of the instrumental and paleo reconstructions at reproducing interannual climate variability and trends over the last decades, in particular the atmospheric circulation, near-surface continental climate (temperature and snow accumulation) and sea-ice cover. To this end, Fig. 2 presents the spatial correlation coefficients between our instrumental and paleo reconstructions and the latest ECMWF atmospheric reanalysis ERA5 (Hersbach et al. 2020) over the 1979-2000 CE period for the 500-hPa geopotential height, near-surface air temperature and snow accumulation, as well as sea-ice concentration using the dataset from the National Snow and Ice Data Center (NSIDC) (Parkinson 2019). Although the period is short (22 years), this evaluation gives a first indication of the performance of our reconstruction when compared to independent data.

Near-surface air temperature
For continental near-surface air temperature, our spatial instrumental reconstruction is highly correlated with ERA5 over West Antarctica (Fig. 2). Our instrumental reconstruction is also close to the reconstruction of Nicolas and Bromwich (2014) Fig. 3 a, b). The high CE and r values mean that both the interannual variability and the amplitude of changes are similar between the two reconstructions. ECR values are largely less than one, which suggests that our reconstruction ensemble has too much variance compared to the mean squared error. However, this result is not surprising since both our instrumental reconstruction and Nicolas and Bromwich (2014) used the same near-surface air temperature dataset for the assimilation. Although the skill metrics are lower for the AP compared to the WAIS, our reconstruction correctly displays the largest warming in AP (0.38 °C per decade) as in Nicolas and Bromwich (2014) (0.46 °C per decade). Comparing our paleo-reconstruction to the instrumental reconstruction shows that using paleo records instead of actual observations degrades the skill of the overall reconstruction, which is expected (Figs. 2 and 3a, b). However, near-surface air temperatures in the paleo-reconstruction still remain highly correlated with ERA5. This is in line with the good agreement with the near-surface air temperature reconstruction from Nicolas and Bromwich (2014) for the AP (r = 0.62 (p-value<0.05), CE = 0.21, ECR = 0.95) and WAIS (r = 0.53 (p-value<0.05), CE = 0.28, ECR = 0.37). For the linear trends over the 1958-2000 CE period, similarly to Nicolas and Bromwich (2014), our paleo-reconstruction displays warming of 0.42 °C per decade and 0.15 °C per decade for AP and WAIS respectively, which is closer than for our instrumental reconstruction (Fig. 3a, b). Our paleo-reconstruction therefore captures the trend better than the instrumental reconstruction, which could be related to the denser paleo network, but not the spatial variability, nor the interannual variability.

Snow accumulation
The snow accumulation from our instrumental reconstruction is less correlated with ERA5 than the near-surface air temperature, and displays a reduced area with statistically significant correlation coefficients (Fig. 2). This could be partially explained by the much higher spatial variability of snow accumulation over the Antarctic continent, compared to near-surface air temperature (van Wessem et al. 2018;Agosta et al. 2019). Additionally, as all atmospheric reanalyses, ERA5 suffers from larger biases in snow accumulation than in near-surface air temperature (Gossart et al. 2019). At first order, the spatial variability of both near-surface air temperature and snow accumulation follows the topography (high temperature and snow accumulation on coastal areas, slowly decreasing towards the interior of the ice sheet), but snow accumulation is more affected by local processes, for example related to the wind (snow erosion/deposition and blowing snow) (Agosta et al. 2019). As a consequence, some snow accumulation records may be poorly correlated with snow accumulation derived from atmospheric reanalyses (Medley and Thomas 2019). This is particularly prevalent Since these snow accumulation records are used in our data assimilation, it is not surprising to obtain lower correlation coefficients for those areas.
When evaluating our snow accumulation reconstruction compared to the snow accumulation reconstruction from Medley and Thomas (2019), the skill metrics show that our instrumental reconstruction displays good and similar performance for the WAIS (r = 0.64 (p-value<0.05), CE = 0.40 and ECR = 0.21) and AP (r = 0.54 (p-value<0.05), CE = 0.24 and ECR = 2.37; Fig. 3c, d). Moreover, both our instrumental reconstruction and the reconstruction of Medley and Thomas (2019) show a positive linear trend in West Antarctic-wide snow accumulation (i.e., the sum of the WAIS and AP regions) over the 1958-2000 CE period: 18.48 Gt year −1 per decade and 11.46 Gt year −1 per decade, respectively (both statistically significant; p-value<0.05). However, it is worth noting that our instrumental reconstruction for the AP displays a smaller linear trend than in the reconstruction of Medley and Thomas (2019) (5.06 Gt year −1 per decade against 16.88 Gt year −1 per decade, respectively). Additionally, the ECR for the AP reaches 2.37, which means that the uncertainty in our reconstruction is underestimated.
Finally, the performance of our paleo-reconstruction for snow accumulation is close to the instrumental reconstruction (Fig. 3c, d). The degradation of the performance between the instrumental and paleo reconstructions is very limited for snow accumulation compared to near-surface air temperature, which is expected as we use the same observational dataset in the two reconstructions of snow accumulation. Moreover, our paleo-reconstruction displays a 18.

Atmospheric circulation
Unlike near-surface air temperature and snow accumulation, no pressure data is used in the data assimilation process for both the instrumental and paleo reconstructions. The skill of the reconstruction of the non-assimilated variables is thus expected to be lower than for the assimilated variables. Spatial correlation coefficients between our instrumental reconstruction and the ERA5 reanalysis for the 500-hPa geopotential height during the 1979-2000 CE period show that interannual changes in the atmospheric circulation are well reconstructed over the West Antarctic and in the midhigh latitudes of the Pacific Sector of the Southern Hemisphere (Fig. 2). Furthermore, for the paleo-reconstruction, the stronger correlations over the ocean, compared to the continent, reflect the fact that changes in this region are driving the variability recorded in ice cores. This confirms that marine intrusions predominantly govern the near-surface climate variability in West Antarctica (Nicolas and Bromwich 2011). This also gives confidence in the use of ice core records for reconstructing the atmospheric circulation in that sector.
We also analyze the Amundsen Sea Low pressure (ASL), which can be characterized by several indexes (e.g., Hosking et al. 2013), for example related to the change in the ASL position (latitudinal and longitudinal positions) or to the intensity of the minimum pressure in the ASL region. Here, the ASL index is defined as the annual average of the 500-hPa geopotential height from 75°S to 60°S and 170°E to 70°W (Fogt et al. 2012;Turner et al. 2013;Hosking et al. 2013). Note that geographical definitions can differ between studies but the resulting indexes are relatively similar. This index is then normalized (zero mean and unity standard deviation) over the 1941-1990 CE period. Following this definition, a negative (positive) ASL index corresponds to a deeper/stronger (weaker) ASL. Because no reconstruction of the ASL index exists, we use the ASL index derived from the ERA5 reanalysis as the reference ASL index for the evaluation.
Over the 1979-2000 CE period, the ASL index from both the instrumental and paleo reconstructions displays a high correlation coefficient with ERA5 (0.87 (p-value<0.05) and 0.65 (p-value<0.05), respectively), as well as CE values (0.74 and 0.30, respectively; Fig. 4a). Although no estimation of the ASL index change exists over the last decades prior the satellite era, (2019) and our instrumental reconstruction are highly sensitive to the analyzed period, reflecting the high internal variability prevailing in this region (e.g., Raphael et al. 2016).
In addition to the ASL, we also analyze the main atmospheric mode of variability in the Southern Hemisphere, i.e., the Southern Annular Mode (SAM). SAM represents the intensity and position of the westerly winds. The SAM index (defined as the normalized difference in mean sea-level pressure between the 40° and 65° South longitude bands) derived from our instrumental reconstruction has a correlation coefficient r of 0.52 (p-value<0.05) with the Marshall index (Marshall 2003) that is only based on atmospheric pressure observations (Fig. 4b). However, our instrumental reconstruction does not display the observed linear trend over 1958-2000 CE as in observations. On the contrary, our paleo-reconstruction is able to reconstruct the linear trend (albeit statistically insignificant) but the correlation coefficient drops to 0.31 (p-value<0.05). The SAM being an annular mode, a high number of geographically evenly distributed observations is required to correctly reproduce this mode. This partly explains the limited skill of our SAM index since the observational network is unequally distributed geographically. Consequently, we mainly focus our analysis on the ASL and not on the SAM.
The best performance of our reconstructions for simulating the ASL compared to the SAM could be explained by two factors. The ASL is a more local atmospheric feature in comparison with the SAM. The observational network (instrumental and paleo observations) used in our reconstructions is limited to the Antarctic continent, with few records in East Antarctica, plus Tasmania, New Zealand and South America, which makes difficult to reconstruct the SAM. In contrast, our observational network includes a large number of records in West Antarctica, which is the region where the ASL has the largest influence. Additionally, the correlation coefficient between the SAM and ASL index in our paleo-reconstruction over the 1800-2000 CE period reaches -0.97 (p-value<0.001; non-sensitive to the analyzed period), meaning our reconstructed SAM variability mainly captures the ASL-related variability. Although our reconstruction seems to overestimate the ASL-SAM relationship, several studies have demonstrated that the ASL is highly correlated with the SAM in observations Hosking et al. 2013;Raphael et al. 2016;Fogt and Marshall 2020).

Sea-ice cover
As for the atmospheric circulation, the sea-ice cover reconstructions rely totally on the relationships between the assimilated variables and sea-ice cover as represented by the model since this variable is not assimilated.  1979-2000CE and 1958-2000 CE periods for the ASL and SAM, respectively. c, d Comparison between our instrumental (red) and paleo (blue) seaice extent (10 6 km 2 ) reconstructions (in anomalies) and the observed sea-ice extent from the National Snow and Ice Data Center (NSIDC; in black; only available from 1979 CE) for the Bellingshausen/ Amundsen Sea (130°W-60°W) and Ross Sea (160°E-130°W) sec-tors over the 1958-2000 CE period. Regional definitions are identical as in Parkinson (2019). Anomalies are computed over the 1979-2000 CE period. The linear sea-ice extent trends are displayed (expressed in 10 6 km 2 decade −1 ) and the standard deviation (std) of the timeseries (expressed in 10 6 km 2 ) over the 1979-2000 CE period for each region. Additionally, the correlation coefficient (r), coefficient of efficiency (CE) and ensemble calibration ratio (ECR) are computed over the respective periods depending on the observations. Error bands correspond to the reconstruction uncertainty as defined in Sect. 2.1. Our SAM index is computed as the normalized difference in mean sea-level pressure between the 40°S and 65°S longitude bands. Note that by construction, the standard deviation of all the SAM and ASL time-series are equal to 1 This is especially true in the offshore regions where the interannual variability of the sea-ice concentration is the highest. Indeed, the ocean along the West Antarctic coast is nearly always covered by sea-ice (not shown). It is especially the case in the Ross Sea sector and to a lesser extent in the Amundsen Sea sector. This means that the interannual variability of the sea-ice concentration in this region is weak. Consequently, the sea-ice concentration in that region cannot explain the variability of the signal recorded in ice cores (e.g., snow accumulation and 18 O ). Therefore, given that the relationship between the ice core signal and sea-ice concentration in that region is weak, the sea-ice concentration in that region cannot be successfully reconstructed.
In order to quantify the sea-ice extent (i.e., oceanic surface covered by at least 15% of ice) at the regional scale, we divide the West Antarctic sector into two sectors following the geographical definitions of Parkinson (2019): the Bellingshausen/Amundsen Sea sector (130°W-60°E) and Ross Sea sector (160°E-130°W). It is important to stress that the overlapped period of 22 years between the NSIDC dataset and our reconstructions is short to evaluate in detail our sea-ice extent reconstruction. Therefore, this limits our ability to assess the quality of our reconstruction in terms of trends over the past decades.
Our instrumental reconstruction displays a correlation coefficient r of 0.43 (p-value<0.05) and 0.51 (p-value<0.05) with the NSIDC dataset over the 1979-2000 CE period for the Bellingshausen/Amundsen Sea and Ross Sea sectors, respectively (Fig. 4c, d). Regarding our paleo-reconstruction, we observe an overall reduction of the local skill compared to our instrumental reconstruction (Fig. 2). However, our paleoreconstruction still displays high correlation coefficients for both the Ross Sea and Bellingshausen/Amundsen Sea sectors when compared to the NSIDC dataset (r = 0.45 (p-value<0.05) and r = 0.35 (p-value>0.05), respectively; Fig. 4c, d). Additionally, the sea-ice extent for these two regions are highly correlated between the two reconstructions over the 1958-2000 CE period (r = 0.67 (p-value<0.05) and r = 0.56 (p-value<0.05), respectively; Fig. S4). Finally, both the instrumental and paleo reconstructions show a decrease in sea-ice extent for the Bellingshausen/Amundsen Sea sector over the 1958-2000 CE period (both -0.057 10 6 km 2 per decade; p-value<0.05), and an increase for the Ross Sea sector (0.018 10 6 km 2 per decade (p-value>0.05) and 0.026 10 6 km 2 per decade (p-value<0.05), respectively for the instrumental and paleo-reconstruction). These changes over 1958-2000 CE are consistent with observations over the 1979-2000 CE period (Parkinson 2019). However, although our instrumental and paleo reconstructions of the sea-ice extent are in agreement with observations regarding the direction of the change, our reconstructions underestimate the magnitude.

Reconstruction of the ASL and sea-ice in the West Antarctic sector over the last 200 years
Now that we have evaluated the skill of our reconstructions during the recent period, we can extend the reconstruction back in time. According to our annual paleo-reconstruction of the ASL index (Fig. 5), the ASL displays three main phases over the last two centuries. From 1800 CE to around 1840 CE, we observe a positive linear trend (trend = 0.34 std per decade; p-value<0.001) following by a period of about one century (until around 1940 CE) with virtually no trend (trend = 0.05 std per decade; p-value>0.05). From around 1940 CE, the ASL displays a strong negative linear trend (trend = -0.30 std per decade; p-value<0.001). This strong negative linear trend is consistent with the observed positive linear trend in the Southern Annular Mode, which deepens the ASL (e.g., Raphael et al. 2016). In addition to being stronger since 1940 CE, our results show that the ASL has also moved eastwards at the same time ( Fig. S5; in our paleo-reconstruction, the longitudinal position of the ASL is anti-correlated with the ASL index over the 1800-2000 CE period; at r = -0.88; p-value<0.001) but no change in latitudinal position is noticed (not shown). Along with the ASL changes over the past two centuries, we analyze annual sea-ice extent changes in the West Antarctic sector, in particular in the Bellingshausen/ Amundsen Sea sector (130°W-60°E) and in the Ross Sea sector (160°E-130°W; Fig. 6a, b). These two regions are anti-correlated at the annual scale in observations over 1979-2000 CE (r = -0.26 (p-value>0.1). In our paleo-reconstruction, this correlation coefficient reaches -0.50 (p-value<0.05) over the same period. However, our paleo-reconstruction indicates that this correlation coefficient varies over the past two centuries. In the Bellingshausen/Amundsen Sea sector, the sea-ice has grown from 1800 CE to around 1860 CE (trend = 0.037 10 6 km 2 per decade; p-value<0.001) before starting a longterm reduction until present (trend = -0.021 10 6 km 2 per decade; p-value<0.001). On the contrary, the Ross Sea sector displays a different picture during the past two centuries. Over the first forty years of the nineteenth century, our reconstruction indicates a reduction of the sea-ice extent in the Ross Sea sector (trend = -0.017 10 6 km 2 per decade; p-value>0.1), which is followed afterwards by a slight long-term reduction of sea-ice extent until around 1940 CE (trend = -0.009 10 6 km 2 per decade; p-value<0.01). Finally, a positive linear trend is observed over the last sixty years (trend = -0.021 10 6 km 2 per decade; p-value<0.05). Therefore, over the last two centuries, the only significant 50-year trend of sea-ice extent in the Ross Sea sector is the sea-ice expansion which has started from about 1940 CE.

Comparison with independent sea-ice reconstructions
Our paleo sea-ice reconstructions can be compared to two independent reconstructions, which are based on phytoplankton productivity emissions from marine algae at the sea-ice edge (Methane Sulphonic Acid; hereafter MSA; Thomas et al. 2019). First, Thomas and Abram (2016) have reconstructed the position of the winter sea-ice edge at 146°W in the Amundsen-Ross sector since 1702 CE using an ice core MSA record situated in West Antarctica. Over the 1800-2000 CE, the sea-ice edge at 146°W derived from our paleo-reconstruction presents a correlation coefficient r of 0.22 (p-value<0.05) with the MSA-based reconstruction from Thomas and Abram (2016)  over the 1900-2000 CE period (Fig. 6d). Over the twentieth century, our reconstruction displays a sea-ice extent northward expansion of 0.13° per decade, which is very similar to Thomas and Abram (2016) ( 0.12 ± 0.02 • per decade). Therefore, although the interannual correlation is relatively low between both reconstructions, the long-term changes are consistent with each other. We also notice that our paleoreconstruction is not able to reconstruct the exceptional seaice expansion for the 1997-1999 CE period (Hanna 2001) in contrast with the MSA-based reconstruction from Thomas and Abram (2016).
The weaker agreement between two reconstructions at the interannual timescale is something to expect when comparing two very different paleo proxies. The MSA emissions arise from the algal blooms that occur in the spring, during the period of sea-ice break-up. Thus, this record is indicative of seasonal changes, closely linked to winter sea-ice extent, rather than annual sea-ice changes. The differences between the reconstructions can also partly be explained by the fact that the MSA record is influenced by conditions at the source as well as the transport pathways and potentially includes a non-climatic signal. In contrast, our data assimilation approach is unable to capture these small-scale local processes and, by construction, our method tends to minimize the non-climatic signal.
Second, using a similar methodology as Thomas and Abram (2016), Abram et al. (2010) have estimated a −0.08 ± 0.02 • per decade southward retreat in the Bellingshausen/Amundsen Sea sector (70°W-100°W) over the 1901-1990 CE period, while our reconstruction suggests a slightly greater retreat (-0.10° per decade; p-value<0.05). Over the twentieth century, the correlation coefficient between the two reconstructions is 0.29 (p-value<0.05). Thus, despite the differences on the seasonal and interannual scale, our reconstruction independently confirms the long-term changes in sea-ice extent over the twentieth century as suggested by Thomas and Abram (2016), including the persistent dipole between the Bellingshausen/Amundsen Sea and Ross Sea sectors (r = -0.70 (p-value<0.001) in our reconstruction over 1900-1999. However, our reconstruction indicates that the sea-ice extent retreat in the Bellingshausen/Amundsen Sea sector started around one century before the sea-ice extent expansion in the Ross Sea sector (Fig. 6a, b).
The evolution of the sea-ice extent in the Ross Sea sector over the past two centuries coincides with the evolution of the ASL, with a correlation coefficient r of -0.82 (p-value<0.001) in our paleo-reconstruction. This confirms the dominant role of the ASL in modulating annual seaice variations in this sector, as highlighted by several studies (e.g., Hosking et al. 2013;Raphael and Hobbs 2014;Raphael et al. 2019). The ASL has deepened (negative ASL index) over the past few decades, associated with enhanced cold southerly winds coming from the Ross Ice Shelf, resulting in sea-ice expansion in the Ross Sea sector (e.g., Lefebvre and Goosse 2008). Conversely, when the ASL is deeper, this atmospheric situation tends to bring more warm air into the Bellingshausen/Amundsen Sea sector in concert with the horizontal movement of the sea-ice to the coast (i.e., ice compaction), which leads to a reduction of sea-ice cover.
However, the sea-ice variability in the Bellingshausen/ Amundsen Sea sector is less correlated with the ASL index than for the Ross Sea sector in our reconstruction over the same period (r = 0.37; p-value<0.001). Additionally, unlike the Ross Sea sector, the correlation between the sea-ice extent and the ASL depends on the analyzed time period. For instance, this correlation coefficient drops to -0.16 (p-value>0.1) over 1851-1950CE (against 0.37 over 1800-2000). During this period, the sea-ice extent for the Bellingshausen/Amundsen Sea sector is decreasing while the ASL index does not present a linear trend. This suggests that the ASL is not the unique driver controlling the variability of sea-ice extent in the Bellingshausen/Amundsen Sea sector. Raphael and Hobbs (2014) pointed out that the Zonal Wave Three and El-Nino Southern Ocean partly govern the variability of the sea-ice extent in the Bellingshausen/Amundsen Sea sector. Furthermore, the analysis of the spatial correlation coefficients between the sea-ice extent in the Bellingshausen/Amundsen Sea sector and the 500-hPa geopotential height based on our paleo-reconstruction (Fig. 7a) indicates that the circulation changes along the Drake Passage (between Cape Horn and the South Shetland Islands) and in the Weddell Sea are correlated with the sea-ice extent in this region. This is confirmed by performing the same analysis with observations (NSIDC and ERA5) over 1979-2000 CE that also gives statistically significant correlation coefficients, but with smaller values (Fig. 7b). Nevertheless, the correlation coefficient between the sea-ice extent in the Bellingshausen/Amundsen Sea sector and the ASL over the 1951-2000against -0.16 over 1851against -0.16 over -1950.1)), pointing out the strong role of the ASL in the reduction of the sea-ice extent in this region over the last decades.

Is the current widespread warming over West Antarctica representative of changes in the past 200 years?
Numerous studies have pointed out a strong west-east asymmetry in the surface warming since 1958 CE as West Antarctica has undergone a large widespread warming, unlike East Antarctica (e.g., Steig et al. 2009;Nicolas and Bromwich 2014;Schneider et al. 2012;Jun et al. 2020). As a result, West Antarctica is often seen as a region that behaves homogeneously in terms of climate change and variability.
However, it is still unknown if this observed widespread warming in West Antarctica over the past decades is representative of the behavior of the region on longer timescales. Our paleo-reconstruction allows examination of this representativeness by analyzing the near-surface air temperature changes over the past 200 years along with changes in other key variables such as snow accumulation, sea-ice extent and atmospheric circulation. Figure 8 displays the linear trends for the 500-hPa geopotential height, 700-hPa air temperature, near-surface air temperature, snow accumulation and sea-ice concentration over 50-year periods between 1801 and 2000 CE. Over 1801-1850 CE, we observe a clear continental dipole in the trends of near-surface air temperature and snow accumulation between the Peninsula and a part of eastern WAIS compared to the western WAIS. The 500-hPa geopotential height field shows that the ASL weakened during 1801-1850 CE (trend = 0.30 std per decade; p-value<0.01), which is consistent with the time-series of the ASL index over the last two centuries (Fig. 5). These atmospheric conditions are associated with weaker northerly winds in the Peninsula that result in a cooling, while the opposite occurs in the western WAIS (i.e., the west of the WAIS Divide where the ice flow is null, which separates the WAIS into two regions). It has been shown that a strong positive correlation exists between temperature and snow accumulation over the Antarctic Ice Sheet (Frieler et al. 2015;Dalaiden et al. 2020a;Cavitte et al. 2020). Higher temperature induces an increase of the saturation vapor pressure that potentially leads to more continental precipitations (i.e., the Clausius-Clapeyron relationship). This is consistent with the cooling in the Peninsula (which corresponds to a snow accumulation decrease) and the warming in the western WAIS (which corresponds to a snow accumulation increase). In addition to the atmospheric circulation changes, continental changes are highly influenced by oceanic conditions (Kittel et al. 2018;Krinner et al. 2008). Particularly, when the ocean is ice-free, the air potentially contains more humidity because of the enhanced evaporation above the ocean. This eventually leads to a warming and an increase in snow accumulation over the continent when the airflow hits the coastline. Accordingly, we observe an increase in sea-ice concentration in the Bellingshausen/Amundsen Sea sector in contrast with a decrease in the Ross Sea sector.
During the second interval (i.e., 1851-1900 CE period), the continental and atmospheric circulation patterns are relatively similar to the one over the 1801-1850 CE period. However, the amplitude of the surface changes is weaker and the ASL seems slightly shifted eastwards (Fig. 8). This induces an eastward shift of the warming pattern over the continent compared to the period 1801-1850 CE. Moreover, the typical dipole in sea-ice is no longer observed since a general reduction of the sea-ice extent is noticed in the West Antarctic sector. The 1901-1950 CE period is characterized by small and insignificant changes, except that a widespread sea-ice cover reduction is noticed as during the 1851-1900 CE period. Over the last fifty years of the twentieth century, the atmospheric circulation corresponds almost to the opposite 1801-1850 CE pattern (Fig. 8). With a deepening ASL (trend = -0.22 std per decade; p-value<0.01), warmer and moister air is brought in the Peninsula resulting in a decrease of sea-ice cover in the Bellingshausen/Amundsen Sea sector, as well as a warming associated with an increase in snow accumulation in the Peninsula. As expected from the deeper ASL, the Ross Sea sector gained sea-ice during the same period while the western WAIS displays a decrease in snow accumulation. However, in contrast to the 1801-1850 CE period, the temperature change has the same sign over the whole West Antarctic region, in agreement with previous studies (e.g., Steig et al. 2009;Nicolas and Bromwich 2014). The boundary between regions showing a warming and those displaying a cooling at 700-hPa is shifted offshore over the ocean while it was located close to the WAIS Divide for the period 1801-1850 (Fig. 8). The dipole is still present but the entire continent is experiencing a warming. The trends of 700-hPa air temperature over 1951-2000 CE on both sides of the Amundsen coasts are thus opposite.
Since the change in the atmospheric circulation during this period roughly corresponds to the opposite pattern of the 1801-1850 CE period, the homogeneous warming over the continent could either indicate a change in the role of the ASL on the West Antarctic climate or a relationship between the ASL and the West Antarctic climate that may have changed over time. Figure 9 displays spatial correlation coefficients between the ASL index and near-surface air temperature, snow accumulation and 700-hPa air temperature over the 1801-1850CE and 1951-2000 Over the 1801-1850 CE period, the spatial patterns of correlation coefficients between the ASL index and the three variables are very similar to the spatial pattern of trends for the same variables and over the same period. This indicates that over that period, the ASL is the major contributor to those changes.
In contrast to the ASL-snow accumulation relationship, the relationship between the ASL and near-surface air temperature has changed between the two periods ( Fig. 9). Correlation coefficients between the ASL index and near-surface air temperature over 1951-2000 CE still exhibit a dipole Fig. 9 Annual ASL index linearly correlated with nearsurface air temperature, snow accumulation and 700-hPa air temperature in our paleo-reconstruction for the 1801-1850CE and 1951-2000 Stippling indicates statistically significant correlations (95% confidence level) Fig. 10 Change in the 500-hPa geopotential height (m) between the 1951-2000 CE and 1801-1850 CE periods in our paleo-reconstruction like for the 1801-1850 CE period, but this dipole is shifted eastward. Additionally, the strength of the correlation over the western WAIS is weaker compared to the 1801-1850 CE period. This change in the link between the ASL and near-surface air temperature could be related to a change in the mean state of the atmospheric circulation. According to our reconstruction and in agreement with previous studies, the atmospheric conditions prevailing in the Antarctic over 1951-2000 CE correspond to a positive phase of the SAM in regards to the 1801-1850 CE period (Fig. 10). Based on the SAM index derived from our reconstruction and normalized over the 1800-2000 CE period, the SAM index averaged over 1801-1850 CE is -0.11 against 0.87 over 1951-2000 CE.
To examine the role of the SAM in the relationship between the ASL and the near-surface air temperature and snow accumulation independently of our reconstruction, we analyze the relationships between the ASL and these two variables in the ERA5 reanalysis over 1979-2019 CE. More specifically, we perform the correlation of climate variables with the ASL index over negative SAM (SAM -) years compared with positive SAM (SAM+) years (Fig. 11). The analysis confirms that depending on the large-scale atmospheric background (i.e., the sign of the SAM index), the impact of the ASL on the near-surface air temperature and snow accumulation is different. During the SAM-years, the fingerprint of the ASL on the continent displays a clear dipole in both near-surface air temperature and snow accumulation. During the SAM+ years, the dipole in near-surface air temperature over the continent is much less clear with weaker correlation coefficients between the ASL and the near-surface air temperature in the western WAIS. This potentially indicates that snow accumulation is a better proxy for reconstructing the ASL than the near-surface air temperature as shown in Thomas and Bracegirdle (2015).
As a consequence of the general warming in West Antarctica over the last 50 years, the western part of the region warms but has lower snow accumulation. This is opposite to the classical link between the two variables due to Clausius-Clapeyron relationship. Cavitte et al. (2020) have shown that some areas over the Antarctic continent in fact display a negative relationship between snow accumulation and near-surface air temperature because of winds-topography interactions. This is for instance observed over the Peninsula (Elvidge and Renfrew 2016;Cavitte et al. 2020). We argue that the same phenomenon at larger spatial scale has occurred in WAIS during the 1951-2000 CE period because of the presence of the WAIS Divide. When warm moist air driven by the ASL is brought in the eastern WAIS towards the WAIS Divide, the air moisture is released as snowfall during the uplift towards the Divide. At that point, the air masses start to descend towards the Ross Ice Shelf on the leeside. This results in an adiabatic warming, which is not associated with an increase in snow accumulation. However, this feature is mainly observed over the 1951-2000 CE period due to the position of the winds at that time (Figs. 8 and S6).

Conclusions
In this paper, we have investigated the West Antarctic climate variability over the past two centuries by providing a complete reconstruction for the main atmospheric variables, as well as the sea-ice concentration. This is achieved by using a data assimilation approach that combines optimally the information brought by paleoclimate proxies and the climate physics from climate models. The reconstruction is based on all the annually-resolved snow accumulation and 18 O records from the Antarctic ice cores, as well as tree-ring width records situated in the mid-latitudes from the Southern Hemisphere. The use of the isotopeenable climate model iCESM1 in the data assimilation allows to avoid strong assumptions on the relationship between the 18 O in ice cores and climate variables.
We have evaluated our reconstruction using various independent data over the instrumental era and shown that our reconstruction is able to reproduce well the climate variability over West Antarctica (near-surface air temperature and snow accumulation), the atmospheric circulation in the West Antarctic sector (in particular, the ASL), and the seaice cover. Snow accumulation and 18 O records from ice cores, as well as tree-ring width, are thus providing valuable information to reconstruct the West Antarctic climate over the past 200 years at least. This is further confirmed by comparing our sea-ice extent reconstructions with two independent sea-ice reconstructions based on MSA records from ice cores, which leads to similar results as our reconstruction at the decadal scale.
According to our reconstruction, the overall observed reduction of sea-ice extent in the Bellingshausen/Amundsen Sea sector over the instrumental era (e.g., Parkinson 2019) is part of a long-term linear trend starting at the beginning of the industrial period. In contrast, the increase of sea-ice extent in the Ross Sea sector has started one century later, around 1950 CE and is preceded by a small sea-ice reduction between 1850 and 1950 CE. Our results show that the deeper ASL over the second half of the twentieth century generally explains these long-term changes by modulating winds in the West Antarctic sector. However, we observe that the strength of the relationship between the sea-ice extent in the Bellingshausen-Amundsen Sea sector and the ASL strongly varies over the past two centuries.
We have also shown that the observed general warming since 1958 CE in West Antarctica is not representative of changes over the past 200 years, which usually present a dipole on the continent due to the ASL. Our reconstruction suggests that this widespread warming is explained by both the stronger ASL and the positive phase of the SAM. Therefore, this could not be simply related to the fingerprint of the ASL impact on the continent or ocean surface, but also the SAM. The continental response to a change in the atmospheric circulation may not be stationary. Care must be thus taken when inferring changes at the West Antarcticascale from local proxies. Care should also be taken when relying on the short instrumental period to infer changes in Antarctic climate variability and its drivers. Future climate scenarios are often reliant on our understanding of climate variability during the instrumental period. But our study demonstrates that this time period may not be the most representative of the climate variability over multi-decadal to centennial timescales.
Additionally, those time-dependent relationships, as well as the potential complex response of the climate in West Antarctica to the atmospheric circulation, raise concerns about simple statistical methods that assume, by construction, constant relationships between the proxy and the reconstructed climate variable. Using a more sophisticated method like data assimilation is a clear advantage because it directly relies on the covariance between the reconstructed variables and the observations as described by the model. Data assimilation-based reconstruction thus ensures that all reconstructed variables are dynamically consistent.
This present work is focused on the ongoing West Antarctic climate by analyzing local changes. Our study demonstrates that the recent period appears unusual in the longer context but further work is needed to understand the drivers and the underlying mechanisms of those changes. Atmospheric circulation changes in the high latitudes of the Southern Hemisphere over the past decades have been attributed to various factors (Thompson et al. 2011), most notably the atmospheric ozone depleting and the increase of the atmospheric greenhouse gases (e.g., Fogt et al. 2009;Polvani et al. 2011). Additionally, an increasing number of studies (e.g., Ding et al. 2011;Steig et al. 2013;Clem et al. 2016;Holland et al. 2019) have shown that changes in the West Antarctic climate are closely linked to the climate variability in the tropical Pacific and other basins through tropical-Antarctic teleconnections. These teleconnections influence the West Antarctic climate by modifying the atmospheric circulation around Antarctica. Therefore, quantifying the contribution of these teleconnections on the West Antarctic climate is a next step to improve the knowledge of the climate changes in West Antarctica over the past centuries.
Assimilating in addition tropical proxies (in particular corals; e.g., Tierney et al. 2015), would likely allow improving the reconstruction of the variability of the Antarctic climate by relying on the tropical-Antarctic teleconnections, particularly in the West Antarctic sector where they are known to be strong (e.g., Thomas et al. 2013. For instance, the tropical proxies could help to reproduce the exceptional sea-ice expansion in the Ross sector observed between 1997 and 1999 CE (e.g., Hanna 2001) since this event has been largely attributed to tropical variability in the western Pacific Ocean. In addition, this new data assimilation-based reconstruction would provide a good framework to study the impact of tropical variability on the Antarctic climate over the past centuries.