Amundsen Sea Embayment accumulation variability measured with global navigation satellite system interferometric reflectometry

Understanding how water vapor is transported to Antarctica and leads to snowfall is crucial for improving ice-sheet mass balance projections and interpreting ice-core records. Global navigation satellite system (GNSS) receivers distributed across Antarctica to monitor ice velocity and solid Earth motion can be used to understand accumulation, ablation, and snow redistribution at the ice-sheet surface on sub-daily timescales. Here, we present a forward model for reflector height change between the GNSS antenna phase center and the snow surface and an inverse framework to determine accumulation rate and near-surface firn densification from the reflector height time series. We use this model to determine accumulation at three long-term GNSS sites located in the Amundsen Sea Embayment (ASE) and at an array of GNSS receivers deployed in the 2007–2008 and 2008–2009 austral summers. From the GNSS interferometric reflectometry (IR) accumulation time series, we find that extreme precipitation (i.e., the 90th percentile of accumulation events) dominates total precipitation and that extreme event frequency varies seasonally. We use our GNSS-IR accumulation time series together with reanalysis products to characterize the atmospheric conditions that promote extreme snowfall in the ASE. The blocking pressure systems that promote extreme accumulation on Thwaites Glacier are facilitated by tropical teleconnections, specifically convection that promotes Rossby wave trains from the western Pacific, Indian, and Atlantic oceans to the Amundsen and Bellingshausen seas.

1 Introduction

The mass balance of ice sheets and glaciers is defined as the difference between ice discharge, surface mass balance (the sum of surface melt, sublimation, and snow accumulation), and grounded ice basal melt or accretion. Currently, the Antarctic Ice Sheet is losing mass at an accelerating pace due to changes primarily in West Antarctic ice discharge driven by oceanic warming and submarine ice-shelf melt (Velicogna et al., 2020; Smith et al., 2020). These changes in ice discharge have partially been offset by an increase in surface mass balance over much of the Antarctic Ice Sheet since the early 19th century due to a net increase in vapor moisture transport to the poles (Kaspari et al., 2004; Frezzotti et al., 2013; Thomas et al., 2015; Medley and Thomas, 2019). The positive trend in surface mass balance is estimated to have mitigated sea-level rise by up to ∼10 mm in the 20th century (Medley and Thomas, 2019).

Positive trends in surface accumulation may be particularly important in the Amundsen Sea Embayment (ASE), where enhanced accumulation may have the potential to mitigate ASE glacier mass loss due to ocean melt (Davison et al., 2023). Observations from shallow ice cores (e.g., Steig et al., 2005) and climate reanalysis (e.g., Davison et al., 2023) suggest that accumulation has increased in the Amundsen Sea Embayment over the last century and that accumulation is strongly linked to spatial patterns of geopotential height in the Amundsen and Bellingshausen Sea sectors (Thomas et al., 2015). Regional variations in geopotential height have been linked to interactions between strong westerly circumpolar flow and ASE topography (Baines and Fraedrich, 1989; Walsh et al., 2000; Lachlan‐Cope et al., 2001) and global tropical teleconnections (Connolley, 1997; Renwick, 2002; Lachlan-Cope and Connolley, 2006; Ding et al., 2012; Turner et al., 2013; Li et al., 2014, 2015, 2021). These interactions combine to form a persistent low-pressure system known as the Amundsen Sea Low (ASL; Raphael, 2004; Turner et al., 2013). The ASL is one of three low-pressure systems distributed along the circumpolar trough (Connolley, 1997; Lachlan‐Cope et al., 2001; Raphael, 2004) and exhibits a seasonal cycle in location and depth moving between the Amundsen and Bellingshausen Sea sectors between austral summer and austral winter (Fogt et al., 2012; Hosking et al., 2013; Turner et al., 2013). Changes in the position of this low-pressure system drive seasonal changes in moisture transport to the Antarctic interior in climate reanalysis (Hosking et al., 2013; Maclennan and Lenaerts, 2021). On average, precipitation in West Antarctica increases in the austral winter when the ASL depth is greatest and the center of the low shifts southeastward closer to the ASE. Precipitation decreases in the summer when the low is weakest and shifts northwestward to the Amundsen–Ross sector (Fogt et al., 2012; Hosking et al., 2013; Turner et al., 2013).

The Amundsen and Bellingshausen Sea sectors are also characterized by large interannual geopotential height variability (Connolley, 1997; Renwick, 2002; Lachlan-Cope and Connolley, 2006; Ding et al., 2012; Turner et al., 2013). The same mechanisms that combine to form the ASL also promote long-lived high-pressure anticyclonic blocking systems (Sinclair, 1996; Renwick, 1998, 2005; Emanuelsson et al., 2018). Interannual blocking variability has been linked to marine air intrusion and accumulation variability recorded in the Roosevelt Island ice core, West Antarctica (Emanuelsson et al., 2018); however, these processes are less well studied in the Amundsen Sea Embayment, where coastal weather station observations of accumulation are limited (Lazzara et al., 2012). Deep and shallow ice cores in the ASE region, though well suited to characterize West Antarctic accumulation variability on interannual and longer timescales (Kaspari et al., 2004), cannot resolve individual snowfall events. Single snowfall events do not directly impact the long-term behavior of polar ice sheets, because each event represents only a small fraction of the integrated decadal to millennial surface mass balance forcing required to convey surface mass balance anomalies to the coast and drive marine outlet glacier advance or retreat (Christian et al., 2020; Robel et al., 2018). Single snowfall events can be used to understand and link snowfall to moisture transport mechanisms that represent the characteristic state of the atmosphere recorded in the snow deposited at the surface. By diagnosing composites of the atmosphere during precipitation events, we aim to increase understanding of the atmospheric dynamics that promote accumulation – particularly extreme accumulation (90th percentile of events) that represent the majority of accumulation in the ASE (Turner et al., 2019; Maclennan and Lenaerts, 2021; Chittella et al., 2022).

GNSS interferometric reflectometry (GNSS-IR) offers an exciting method to determine accumulation rates with high temporal resolution in remote and undersampled regions of the West Antarctic interior. GNSS positioning uses coded signals transmitted by GNSS satellites at different L-band frequencies to determine three-dimensional position estimates of GNSS receiver antenna locations. The timing of received GNSS signals provides range information important for establishing antenna positions, but it is distorted by electromagnetic wave interactions with the neutral atmosphere, the ionosphere, and reflections off nearby surfaces before the wave reaches the antenna (i.e., multipath reflections).

Figure 1 (a) Overview map of long-term GNSS-IR sites (LTHW in orange, KHLR in purple, and UTHW in blue), integrated vapor transport calculated between 300 and 900 hPa pressure levels, the 500 hPa geopotential height (contours), and an atmospheric river event detected on 13 May 2016. (b) Subpanels show the receiver position as a function of time. (c) Photo of the LTHW receiver and power system (credit: POLENET team). (d) Reflector height change in 2016 for all three long-term GNSS sites. Station-matched colors demark reflector heights observed at each station, and black points are the daily median used as input to the accumulation inverse algorithm. Gray bars indicate the timing of atmospheric river events that made landfall in the ASE identified using methods described in the Appendix corroborated by the catalog developed by Wille et al. (2021).

The distorted travel times associated with multipath reflections have conventionally been treated as noise and informed early GNSS antenna gain pattern designs, which minimized multipath signals at high satellite elevation angles. Larson et al. (2009) established that multipath signals embedded in the signal-to-noise ratio (SNR) of GNSS positioning data contain useful information about the reflection surface geometry (i.e., the height of the antenna above the surface) and the dielectric properties of the surface. Many subsequent studies have shown that the distance between the reflection interface and the antenna phase center can be calculated and monitored over time from the SNR recorded at the GNSS antenna in order to understand local surface evolution (Larson et al., 2015; Siegfried et al., 2017; Shean et al., 2017; Larson et al., 2020; Pinat et al., 2021).

Despite the proven effectiveness of GNSS-IR methods to measure surface height and reflector amplitude changes associated with snow deposition, only two studies (Larson et al., 2015; Siegfried et al., 2017) have used GNSS-IR to link surface height changes to multi-year precipitation and climate variability on ice sheets. This is largely due to the logistical challenges associated with supporting multi-year on-ice GNSS receiver networks, where stations must be raised regularly (every 1–3 years) to prevent burial due to snow accumulation. Over the past 15 years, GNSS receivers located in the ASE (Figs. 1 and 2) have been maintained to support several different projects investigating ice flow (e.g., Hoffman et al., 2020) and the interaction between glacial isostatic adjustment and solid Earth deformation (e.g., Barletta et al., 2018). Through the application of GNSS-IR, more than 10 years of nearly continuous GNSS operation also offers a unique opportunity to explore atmospheric dynamics that promote and control precipitation and accumulation in the ASE.

We use GNSS-IR SNR data from three long-term sites deployed and maintained in the Amundsen Sea Embayment, as well as GNSS campaign arrays that were deployed across Thwaites Glacier as part of geophysical surveys in the austral summers of 2007–2008 and 2008–2009 to link extreme precipitation events and synoptic-scale atmospheric rivers observed in reanalysis to highly localized accumulation records at individual GNSS sites (Figs. 1, 2). Using our unique GNSS-IR records, we quantify the seasonal cycle of precipitation and elevation dependence of accumulation in the ASE. We also verify that extreme precipitation events supply the majority of the total accumulation in the ASE. We then use the GNSS-IR records of extreme accumulation events to composite near-surface temperature, atmospheric pressure, and integrated vapor moisture transport simulated in reanalysis in order to better understand the drivers of extreme precipitation in the ASE. We use these composites to understand the synoptic conditions and global teleconnections that promote extreme snowfall in the ASE. Our results suggest that extreme snowfall events observed in the ASE are coincident with blocking conditions in the Bellingshausen Sea.

2 GNSS interferometric reflectometry

Here we report on multi-year deployments for three on-ice GNSS sites in West Antarctica's Amundsen Sea Embayment: two along the trunk of Thwaites Glacier (upper Thwaites, UTHW, and lower Thwaites, LTHW) and one on the nearby Kohler Glacier (KHLR). The UTHW site was deployed in the 2008–2009 austral summer. The LTHW and KHLR sites were deployed in the 2009–2010 austral summer. Although there are intermittent data gaps due to power loss, all three sites have been regularly maintained and remain active today. In addition to the long-term sites, we also apply GNSS-IR to GNSS sites that were deployed on Thwaites Glacier during the 2007–2008 and 2009–2009 austral summers. During the 2007–2008 and 2008–2009 austral summers, 17 sites were deployed along two flowlines, which aimed to measure short-timescale (minutes to days) ice-velocity fluctuations from the coast towards the West Antarctic Ice Sheet central ice divide and Bentley Subglacial Trough. When possible, GNSS monuments installed in the 2007–2008 austral summer were reoccupied in the 2008–2009 austral summer; if a monument could not be found, the GNSS site was deployed at the position of the 2007–2008 monument as calculated from satellite ice-velocity data (Joughin et al., 2009). The 2007–2008 array was operational in January 2008. The 2008–2009 array was operational in December 2008 and January 2009. All GNSS site locations are plotted in Fig. 2a.

Reflector height time series were calculated using the workflow described by Roesler and Larson (2018), which is based on the earlier work of Larson et al. (2015). Daily solutions for receiver position were calculated kinematically relative to the sites of nearby base stations (TF01 and WAIS for the campaign sites; permanent GNSS sites at Backer Island and Howard Island for the long-term sites) using differential carrier-phase positioning estimates from the TRACK software (Chen, 1998; King, 2004; Hoffman et al., 2020). Daily averaged position solutions for the base station sites were determined in GAMIT/GLOBK following standard GNSS processing methods for glaciological applications (King, 2004). The reflector height was calculated from the signal-to-noise ratio (SNR) data assuming a planar reflector surface geometry. After direct signal effects are removed, the SNR data for a single satellite and receiver can be modeled as

$$\begin{array}{}\text{(1)}& \text{SNR}\left(e\right)=A\left(e\right)\sin \left({\displaystyle \frac{\mathrm{4}\mathit{\pi }{H}_{\mathrm{R}}}{\mathit{\lambda }}}\sin \left(e\right)+\mathit{\varphi }\right),\end{array}$$

where A is the amplitude, which depends on the transmitted GNSS signal power, the satellite elevation angle above the horizon, the surface dielectric constant, and its roughness; H_R is the reflector height; λ is the wavelength of the transmitted GNSS signal; and ϕ is a phase constant. This representation means a periodogram can be used to estimate H_R from the dominant frequency:

$$\begin{array}{}\text{(2)}& f={\displaystyle \frac{\mathrm{2}{H}_{\mathrm{R}}}{\mathit{\lambda }}}.\end{array}$$

To resolve this frequency, the antenna must be at least two wavelengths above the surface (0.4–0.5 m, for GPS L1 and L2 frequencies). The elevation and azimuth angles were calculated from the broadcast GPS ephemeris. The dominant SNR frequency was determined using a Lomb–Scargle periodogram. On average, this results in ∼70 estimates of SNR frequency every day, which can be related to reflector height using Eq. (2).

We eliminated weak signals for each time increment that are 0.3 m from the median reflector height for that time increment. Longer sample intervals for calculating reflector heights provide a more robust result due to the larger number of samples in the increment. For the long-term sites, we used daily averaged reflector heights to determine daily accumulation rates. See Table S1 in the Supplement for a chart detailing the days when we have snowfall observations from GNSS, with previous observations from acoustic depth gauges available from the Automatic Weather Station project (Lazzara et al., 2012). We used 3 h estimates of reflector height to snow accumulation moving across the campaign GNSS array because station spacing was sufficient to map storms that moved up Thwaites Glacier (Video supplement; Hoffman, 2023). The mean and standard deviation of the filtered signals were used to derive accumulation records using the inverse method described in the following section.

Figure 2 (a) Locations of all the GNSS sites (long-term sites indicated with stars and campaign sites indicated with circles) with reflector height observations from the (b) 2007–2008 array, (c) 2008–2009 array, and (d) long-term GNSS sites. For the long-term records, all reflection heights are shown with daily median reflectors in bold. Atmospheric river events are marked with dashed red vertical bands (see Appendix A2 for more details on event detection methods).

3 Methods

3.1 Accumulation and densification forward model

To determine the accumulation record from observations of reflector height, we use an inverse method (Metropolis–Hastings algorithm) that requires a forward model that translates the parameters of interest (in our case, an accumulation and densification time series) into observables (reflector height time series). Our forward reflector height model for on-ice GNSS sites represents one-dimensional multi-layer snow compaction and grain settling from the surface to the depth of the GNSS antenna monument base. Each snow layer in the model is described by its thickness and density. Layers evolve following the compaction constitutive relation:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial {\mathit{ϵ}}_i}{\partial t}}=-{\displaystyle \frac{{\mathit{\sigma }}_i}{\mathit{\eta }}},\end{array}$$

where ϵ_i is the layer strain, η is the snow viscosity, and σ_i is the overburden stress of the layers above (Fig. 3). The vertical stress of the overlying layers is computed for each layer, i, as

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_i=\sum _{\mathrm{1}}^{i-\mathrm{1}}g{\mathit{\rho }}_i{h}_i\cos \left(\mathit{\theta }\right),\end{array}$$

where θ is the local surface slope; g is the terrestrial gravitational constant; and ρ_i and h_i are the layer density and thickness, respectively. The snow viscosity, η, which describes how efficiently stress is transmitted through the column, has been linked to snow temperature, liquid water content, and snow grain geometry. Without in situ observations of snow properties, we fix the viscosity for the duration of each reflector height experiment.

Figure 3 The model schematic and convergence pattern of the Metropolis–Hastings algorithm used to produce an accumulation time series. In the main plot, we show the convergence of the modeled RH (colored lines) and the observed reflector height time series (black points) as a function of model iterate. The objective function we use to score solutions is shown as an inset plot. The misfit of the model compared to observations decreases and converges as a function of iteration in the inverse procedure.

3.2 Reflector height inverse method

The number of unknown parameters in our model scales with the length of the reflector height record. We reset the inference of accumulation, viscosity, and initial snow density each time the antenna mast is raised (quasi-annually). We use Bayes' theorem to compute the distributions of accumulation at each time step and specify a likelihood model, assuming that observations of reflector height are independent in time and normally distributed around the true value at each data point with variance ${\mathit{\sigma }}_{R_{h_{{\mathrm{obs}}_k}}}$. The objective functional defined by the least-squares difference between observed, $R_{h_{{\mathrm{out}}_k}}$, and modeled, $R_{h_{{\mathrm{obs}}_k}}$, reflector heights for each increment, k, in the time series when observations are available is described by

$$\begin{array}{}\text{(5)}& {\mathit{\varphi }}_k\left(m_k\right)=\left(R_{h_{{\mathrm{out}}_k}}\right(m_k)-R_{h_{{\mathrm{obs}}_k}}{)}^{\mathrm{2}}.\end{array}$$

We use this objective functional to define the likelihood of the modeled reflector heights, given the observations:

$$\begin{array}{}\text{(6)}& p\left(R_{h_{\mathrm{obs}}}\right|m)\propto \prod _{k\in n}\exp \left({\displaystyle \frac{-B\mathit{\varphi }\left(m_k\right)}{{\mathit{\sigma }}_{R_{h_{{\mathrm{obs}}_k}}}^{\mathrm{2}}}}\right).\end{array}$$

The priors for the accumulation record are modeled with Gaussian distributions with parameters listed in Table S3. The Gaussian distributions for the snow accumulation are centered on an initial estimate of mean snow accumulation that takes the maximum of the difference of snow surface height divided by the mean snow density (300 kg m⁻³) and a prior of no snowfall ($\mathrm{300}\cdot max(R_{h_{{\mathrm{obs}}_{k+\mathrm{1}}}}-R_{h_{{\mathrm{obs}}_k}},\mathrm{0})$). With a likelihood model and a prior estimate of the model parameters, we use Monte Carlo methods to sample and thereby estimate the posterior distribution of the accumulation, density, and snow viscosity, and we evaluate the statistical properties of the samples as a proxy for the joint posterior distribution of the model parameters. We use the Metropolis–Hastings algorithm to build estimates of the posterior distribution, map model degeneracy, and understand parameter uncertainty. The Metropolis–Hastings algorithm samples the prior estimate of each model parameter and independently updates each parameter at each iteration by randomly drawing a sample from the proposed distribution. If the sampled posterior probability is greater for the proposed distribution than the current distribution, the new parameter value is accepted. Otherwise, it is accepted with probability (P_accept) proportional to the ratio of the current probability (p(m_i)) and that of the proposed value (p(m_i+1)).

$$\begin{array}{}\text{(7)}& {\displaystyle }{P}_{\mathrm{accept}}=min\left(\mathrm{1},{\displaystyle \frac{p\left(m_{i+\mathrm{1}}\right)}{p\left(m_i\right)}}\right)\text{(8)}& {\displaystyle }\begin{array}{rl}{P}_{\text{accept}}& =min\left(\mathrm{1},{\displaystyle \frac{A\exp (-B\mathit{\varphi }(m_{i+\mathrm{1}}\left)\right)}{A\exp (-B\mathit{\varphi }(m_k\left)\right)}}\right)\\ & =min\left(\mathrm{1},\exp \left(-B\left(\mathit{\varphi }\right(m_{k+\mathrm{1}})-\mathit{\varphi }(m_k\left)\right)\right)\right)\end{array}\end{array}$$

Here B is a tunable regularization parameter to facilitate convergence. We use these time series to explore possible drivers of precipitation by connecting the events we observed using GNSS-IR to extreme precipitation events in reanalysis.

3.3 Reanalysis composites

Following many previous studies, we use MERRA-2, RACMO2, and ERA5 reanalysis products to contextualize our observations of local surface height change from GNSS with large-scale circulation (Gelaro et al., 2017; van Wessem et al., 2018; Hersbach et al., 2020). We chose to focus on ERA5 because past studies have shown that observed snow accumulation from GNSS-IR agreed better with ERA5 snow accumulation records than with MERRA-2 and RACMO2 snow accumulation records (Fig. S1; Table S3 in the Supplement; Maclennan et al., 2022). ERA5 is the newest global atmospheric reanalysis product from the European Centre for Medium-Range Weather Forecasts (ECMWF) and has 0.25° latitude (28 km) by 0.25° longitude (∼7 km²) resolution (Hersbach et al., 2020). We focus our study on the southern Pacific but include reanalysis data from tropical latitudes (20° N to 90° S) to evaluate the possibility of tropical teleconnections that have been suggested to affect precipitation variability in accumulation records preserved in firn and ice cores.

Using the GNSS-IR time series, extreme precipitation was classified using a seasonally corrected 90th percentile threshold of daily accumulation rates for the lower Thwaites and upper Thwaites site accumulation events. The union of these two indices was used to create a Thwaites Glacier catchment extreme accumulation event dataset. The results do not change substantially when the Kohler Glacier site is included in this index; however, because of the observed offset in the accumulation seasonality, we exclude Kohler from the merged index of the sites on Thwaites Glacier. The GNSS-IR catalog of extreme accumulation events was used to calculate composite maps of geopotential height and wind anomalies at the 500 hPa pressure level, blocking anomalies calculated from anomalies of geopotential height at the same pressure level, surface temperature anomalies (2 m), near-surface wind anomalies (10 m), integrated vapor moisture transport (IVT) anomalies, and Rossby source anomalies using ERA5 daily averages. Rossby source anomalies were calculated from 200 hPa winds following the methods described by Trenberth et al. (1998) and Holton (2004), as implemented by Dawson (2016). Daily average outgoing longwave radiation (OLR) anomaly composites and sea surface temperature anomalies were calculated from NOAA SST data and OLR satellite data (Liebmann and Smith, 1996). Composites were calculated for the 10 d preceding the onset of extreme precipitation and corrected for anomalies calculated relative to seasonal mean conditions. The interval of 10 d was chosen to evaluate the importance of both proximal atmospheric dynamics near the Amundsen Sea and tropical and extratropical teleconnections that precede changes in the ASE by up to 10 d (Li et al., 2015, 2021). As part of this analysis, we also built catalogs of blocking events and the atmospheric rivers that made landfall in the ASE, both of which have been connected to extreme precipitation and changes in the surface elevation of the West Antarctic Ice Sheet (Adusumilli et al., 2021). More information on the methods used to create these catalogs can be found in Appendix A1.

4 Results

4.1 GNSS-IR accumulation time series

Using GNSS-IR, we created the first ∼10-year time series of snow accumulation at three sites in one of the most remote areas of West Antarctica: the Amundsen Sea Embayment (Fig. 4). Currently, all precipitation on Thwaites Glacier at the GNSS sites falls as snow with no discernible rain events recorded in reanalysis or GNSS-IR surface reflection amplitude. The seasonal accumulation climatology is significant for all sites (p value <0.001; see Supplement). On Thwaites, precipitation is most frequent and intense in the austral winter (June, July, and August; JJA), with peak accumulation in August that tapers into the austral summer to a minimum in January (Fig. 5). At LTHW, precipitation is ∼50 % higher in austral winter than in austral summer. The seasonal signal is more muted at UTHW, with only a ∼30 % difference between seasons. At the KHLR site, we observe peak accumulation in March and minimum accumulation in November. The interannual variability and the seasonal cycle are driven primarily by changes in the frequency of extreme precipitation events. Extreme events account for over half the total accumulation at all three long-term GNSS sites and explain 80 % of the snow accumulation variability in reanalysis. These events also change seasonally (Fig. 5). Over the observation time period (2009–2020), extreme precipitation events occur on 13 % of wintertime days (96 d), 6.2 % of summertime days (45 d), 10.4 % of springtime days (71 d), and 10.2 % of fall days (82 d).

Figure 4 Observed cumulative accumulation at the LTHW, UTHW, and KHLR receiver sites with the timing of atmospheric river events that made landfall in the ASE (dashed red vertical lines) and the extreme accumulation events measured with GNSS-IR that occur within 2 d of the atmospheric river events (gray vertical lines).

Figure 5 (a) The seasonal accumulation cycle at LTHW, UTHW, and KHLR GNSS sites. (b) Distribution of the difference between seasonally high accumulation event frequency and low accumulation event frequency at the upper and lower Thwaites Glacier (DJF–JJA) and Kohler Glacier (SON–MAM).

Using the short-duration GNSS arrays that were deployed for 1 month in the 2007–2008 austral summer and 2 months in the 2008–2009 austral summers, we also document the distributed patterns of snowfall across 17 stations in the ASE (Figs. 2, 6). During these short-duration campaigns, we observed seven extreme precipitation events during these two summers (Fig. 6). Across the short-term arrays in both summers, we find that it snows 2–4 times more at lower Thwaites (800 m above sea level) than in the higher-elevation regions of the catchment (1200 m above sea level) during extreme events (Fig. 6). Elevation gradients in accumulation are less strong on days when we do not observe extreme events.

Figure 6 (a) The guiding 500 hPa geopotential height (gray contours) and the integrated vapor transport for the 2007–2008 array atmospheric river event with the cumulative accumulation at the site of each receiver deployed as part of the (b) 2007–2008 and (c) 2008–2009 GNSS arrays. Red lines indicate the timing of atmospheric rivers that make landfall in the ASE. Gray shading indicates the timing and duration of extreme precipitation events determined from all array GNSS sites. Average daily snowfall as a function of elevation for the (d) 2007–2008 and (e) 2008–2009 austral summers.

4.2 Composite investigation of reanalysis

Summer and winter time blocking in the Amundsen and Bellingshausen seas

The mean blocking conditions from the days preceding the onset of observed extreme accumulation at the GNSS sites for summer (DJF) and winter (JJA) events are illustrated in Fig. 7. We found that between 2 d prior to (day −2) and the onset of snowfall (day 0) blocking is present in the Bellingshausen and Amundsen seas (defined by bounds 55° S, 75° S, 120° W, and 45° W) during 42 % of winter extreme precipitation events and 37% of summer extreme precipitation events. Snowfall events with atmospheric blocking also have statistically higher mean precipitation rates than those without blocking.

The blocking frequency anomalies associated with both summer and wintertime extreme precipitation events reveal significant blocking of the westerly flow in the Bellingshausen and Amundsen seas before and during extreme snowfall events (Fig. 7). The composites of 500 hPa winds are also consistent with blocking over large areas of the Bellingshausen Sea during extreme snowfall events. During the summer extreme precipitation events, blocking frequency anomalies are higher (note the difference in color scale) compared to wintertime blocking anomalies associated with extreme wintertime precipitation (Fig. 7). Austral summer blocking anomalies also cover a significantly larger area than the austral winter anomalies and extend westward to much of the Amundsen Sea. In the austral winter, anomalously negative blocking frequency is observed closer to the Thwaites Glacier catchment. This low-pressure austral winter anomaly drives anonymous northerly mid-level winds that impinge on the coast of West Antarctica and appear to strengthen along the eastern flank of the anticyclone.

Figure 7 (a) Blocking frequency anomaly, (b) 500 hPa geopotential height shown with vectors of the 500 hPa winds, and (c) 500 hPa geopotential height anomaly and wind anomaly for extreme precipitation events observed during the austral summer (DJF). (d) Blocking frequency anomaly, (e) 500 hPa geopotential height and winds, and (f) 500 hPa geopotential height and wind anomaly observed during the austral winter (JJA). The day indices are relative to the day of precipitation onset at the GNSS receivers. Stippling of anomalies indicates significance at 95 % confidence.

Figure 8 (a) Seasonal blocking frequency cycle and the (b) blocking frequency anomalies associated with extreme precipitation on Thwaites Glacier.

The mean blocking frequency (2009–2020, Fig. 8) across all seasons in the Southern Ocean is largest in the Amundsen and Bellingshausen seas. Blocking frequency in this region increases and shifts zonally from the Amundsen Sea in the summer to the Bellingshausen Sea in the winter. The amplitude of blocking frequency change also increases with a maximum in winter and a minimum in the preceding fall. The blocking frequency anomalies associated with extreme accumulation events over the period of observations are highest during summer and spring when the total frequency of blocking is low. This clear pattern in blocking frequency anomaly suggests that blocking of moist mid-latitude westerly flow north of the Amundsen and Bellingshausen seas in summer and spring drives significant summer and spring precipitation. In the next section, we use composites of near-surface temperature and integrated vapor moisture transport to explore these relationships further.

Figure 9 (a) Near-surface (2 m) temperature anomaly, (b) mean integrated vapor moisture transport, and (c) integrated moisture transport anomaly for extreme precipitation events during the austral summer (DJF). (d) Near-surface temperature anomaly, (e) mean integrated vapor moisture transport, and (f) integrated vapor moisture transport anomaly for extreme precipitation events during the austral winter (JJA). The day indices are relative to the day of precipitation onset at the GNSS receivers. Stippling of anomalies indicates significance at 95 % confidence.

Figure 10 The detected atmospheric rivers (cyan-green contour) with the guiding 500 hPa geopotential height (gray contours) and the integrated vapor transport for the events (color map) that occur within 3 d of the 95th percentile of height changes observed with GNSS-IR. Graphs to the right of each figure show the accumulation determined using GNSS-IR reflector height change during each event. Red vertical lines indicate the timing of atmospheric river landfall.

4.3 Summer and winter time vapor transport in the Amundsen Sea

The mean surface temperature and vapor transport conditions on the days preceding the onset of extreme precipitation events we observe with GNSS-IR are illustrated in Fig. 9. Positive near-surface temperature anomalies are present across the Amundsen Sea and the southern Pacific 2 d before extreme precipitation is observed on Thwaites Glacier. These temperatures steadily increase and peak on the day of the extreme event with positive temperature anomalies exceeding 1.5 °C (Fig. 9). During austral summer (DJF), these near-surface temperature and wind anomalies appear to come from the northwest. During the austral winter (JJA), anomalously high near-surface temperatures increase in the Amundsen Sea from the northeast (Fig. 9a, d).

The mean integrated vapor transport (IVT) associated with these events also shows enhanced values in the Amundsen and Ross seas. The circum-Antarctic jet that transports the majority of high-latitude moisture in the Southern Hemisphere is diverted southward expressed as a negative easterly IVT anomaly where IVT is highest (Fig. 9b, e, c, f). This negative moisture anomaly is coherent with a positive westerly and northerly IVT anomaly at higher latitudes that extends into the West Antarctic interior (Fig. 9c, f). In austral summer, the positive IVT anomalies are more uniform and diverted to higher latitudes by similarly uniform negative westward IVT anomalies. In austral winter, IVT anomalies integrated over the area of the entire Pacific are stronger than in austral summer; however, the spatial structure of the anomalies is more heterogeneous and insignificant at 95 % confidence.

Atmospheric river event frequency also varies seasonally, with larger events making landfall in austral winter months. These events make up only a small fraction of the extreme events that appear to control precipitation on Thwaites. Similar to the conclusions of other studies (i.e., Maclennan and Lenaerts, 2021), we find that atmospheric rivers account for 8 % of the total accumulation measured at the sites of the long-term GNSS receivers over the 10 years of available observations. In total, we identify 12 atmospheric river events, defined as continuous periods (multiple hours to multiple days), when rivers make landfall within 5 d of an extreme accumulation event observed with GNSS-IR (90th percentile). The 12 distinct events are shown in Fig. 10 with the accumulation that we observe at the inland GNSS-IR sites.

4.4 Summer and winter teleconnections

Composite sea surface temperature and 500 hPa geopotential height anomalies were analyzed across the Southern Hemisphere to test the significance of tropical teleconnections that may drive extreme winter and summertime precipitation in the Amundsen Sea. In the austral winter, we find significant positive tropical Atlantic, southeastern Atlantic, and southwestern Indian Ocean sea surface temperature anomalies preceding extreme snowfall on Thwaites Glacier (Fig. 11, Fig. S3). We also identify a clear train of alternating high- and low-pressure anomalies consistent with Rossby waves that migrate eastward, building up to extreme precipitation events (Figs. 12, S3, S4, S5, S6, S7). In the austral winter, these Rossby wave trains, as expressed in 500 hPa pressure anomalies, extend from large anomalies off the coast of Africa in the eastern Atlantic 4–5 d preceding the extreme precipitation event (Fig. 12). These waves propagate across the Southern Ocean to the Amundsen and Bellingshausen seas. In austral summer, the anomalies associated with extreme precipitation events are stronger and more significant in the Bellingshausen Sea and appear to be driven by significant summer sea surface temperature anomalies in the central Pacific, in addition to significant anomalies in the Atlantic and western Indian oceans (Figs. S2, S3).

Figure 11 Austral winter sea surface temperature composite anomalies. The numbering of the days indicates the days before the extreme precipitation event observed at the GNSS-IR stations located on Thwaites Glacier, which are plotted as orange and blue stars.

Figure 12 Austral winter geopotential height composite anomalies at the 500 hPa pressure level. The numbering of the days indicates the days before the extreme precipitation event observed at the GNSS-IR stations located on Thwaites Glacier, which are plotted as orange and blue stars.

5 Discussion

The character of the Amundsen Sea Low, specifically the magnitude of the low-pressure anomaly and the location of the anomaly, has historically been the climatological index most well linked to accumulation variability in the Amundsen Sea Embayment and has been used to suggest connections between synoptic atmospheric variability and West Antarctic precipitation. These indices represent a mean anomaly integrating pressure variations over many days and many different synoptic conditions, which complicates the interpretation and attribution of signals that might promote extreme events. Chittella et al. (2022) showed that the second principal component of the geopotential height at 850 hPa, which is correlated to the 500 hPa geopotential height analyzed in our study, is highly correlated (r²=0.9) with the zonal position of the Amundsen Sea Low and explains 43 % of the extreme precipitation events observed over this period. Our data, in conjunction with interpretations of regional reanalysis and other studies that use similar reanalysis products (Maclennan and Lenaerts, 2021; Chittella et al., 2022), suggest blocking highs in the Bellingshausen Sea drive extreme vapor transport into the West Antarctic interior (Figs. 10, 7, and 8). This pattern has previously been shown to be responsible for high-snowfall events over Thwaites Glacier (Maclennan and Lenaerts, 2021; Chittella et al., 2022), and these reanalysis studies agree with accumulation observations from GNSS interferometric reflectometry presented here, though, in general, reanalysis underestimates extreme event frequency.

The GNSS-IR time series suggests that enhanced snowfall events observed in the ASE are associated with the blocking of the westerly polar jet stream and resulting perturbations in the composite integrated water vapor transport. Patterns in vapor transport anomalies associated with extreme precipitation on Thwaites Glacier suggest that the polar storm track is blocked and that storms are rerouted poleward with distinct southward transport toward the Amundsen Sea. This is similar to the behavior observed over Greenland in the Northern Hemisphere, where mid-latitude high-pressure systems over Europe and the eastern Atlantic have been observed to divert storm tracks that then promote snowfall over Greenland (Pettersen et al., 2022).

Blocking in the Amundsen and Bellingshausen seas exhibits a seasonal cycle that may explain the frequency of extreme precipitation, which controls over half of the total accumulation we observe with GNSS-IR in the Amundsen Sea Embayment. Indeed, 40 % of extreme precipitation events were accompanied by blocking in the Amundsen and Bellingshausen sectors, and blocking frequency anomalies associated with extreme precipitation become more significant in summer and spring when blocking is generally less frequent than in winter. Increases in winter extreme precipitation event frequency compared to summer appear to be the result of Rossby wave breaking from convection northwest of the Amundsen and Bellingshausen seas in the Atlantic, Indian, and western Pacific Ocean basins. These Rossby waves form a low-pressure anomaly to the west of Thwaites Glacier and a high-pressure anomaly to the east of Thwaites Glacier, funneling moisture from the tropics between the cyclone and anticyclone. In the austral summer, extratropical Rossby wave propagation still appears to play a role in routing moisture to the Amundsen Sea; however, the anomaly associated with the blocking high-pressure system in the Bellingshausen and Amundsen seas extends further westward than the austral winter high-pressure system. This is consistent with Rossby wave propagation from the central Pacific and the significant positive sea surface temperature anomalies observed there preceding extreme summer precipitation.

Previous numerical modeling experiments investigating tropical teleconnections to the Southern Hemisphere have suggested that Rossby wave propagation to the Amundsen Sea is driven primarily from the Indian and Atlantic basins northwest of the Amundsen Sea (Li et al., 2014; Simpkins et al., 2014; Li et al., 2015). Li et al. (2015) used idealized CAM4 simulations to show that anomalous convection events in the extratropical ocean basins in the Southern Hemisphere promote Rossby wave propagation that sums linearly in the southern Pacific Ocean. This complicates connecting the pressure anomalies and mechanisms that promote blocking in the Amundsen Sea to Rossby wave source mechanisms from specific basins. There are likely contributions from many different basins. This may explain the more complex pressure anomaly pattern observed in the summer when extreme precipitation is preceded by positive sea surface temperature anomalies in the central and western Pacific Ocean in addition to the Indian Ocean. During our study period, the pressure anomalies in winter, when it snows the most in the Amundsen Sea, resemble Rossby wave trajectories initiated from the subtropical Atlantic. Li et al. (2014) found the trajectory of the Rossby wave ray path to be similar in all seasons, first propagating poleward from the subtropical South Atlantic until waves were reflected by weak vorticity gradients in the polar region. In the austral summer, Li et al. (2015) found that rays continue propagating poleward until they dissipate at high latitudes; in austral winter, however, the rays interact with a second reflecting surface created by anticyclonic shear associated with the poleward flank of the subtropical jet and the equatorial flank of the mid-latitude jet stream. This reflecting surface is sufficiently strong in idealized model representations of Rossby wave propagation to steer the ray path of the waves to the Amundsen/Bellingshausen region. Significant sea surface temperature composite anomalies observed in the Atlantic and western Indian Ocean preceding extreme precipitation are consistent with the propagation paths of eastward propagating Rossby waves modeled in the Southern Ocean (Li et al., 2015).

The advection of moist air masses by blocking anticyclones leads to precipitation in West Antarctica and represents the dominant positive contribution to the Amundsen Sea Embayment surface mass balance. The anomalously warm 2 m air temperatures and water vapor transport along the coast of the continent also suggest that these pressure systems promote the transport of warm air and clouds that can insulate the continent when the snow is deposited. This has implications for the interpretation of ice cores in West Antarctica and the assimilation of ice-core data into paleoclimate reanalyses because the fractionation of isotopes is typically related to the local temperature of the precipitating cloud. The seasonality in accumulation we observe in the GNSS-IR time series and reanalysis suggests that austral winter conditions and associated inferences of atmospheric temperature are sampled more than austral summer conditions in firn and ice cores. Our results also suggest that archives from West Antarctica for all seasons are likely biased toward the conditions associated with extreme precipitation.

The blocking pressure field that has been identified here and in previous studies as promoting extreme precipitation has also been connected to melting in West Antarctica and may increase the likelihood of rainfall instead of snowfall at the surface as this sector warms (Niwano et al., 2021). At present, the processes contributing water to the near-surface cannot sustain a perennial firn aquifer, and melt refreezes due to the cold content of the firn. Meltwater refreezing raises the temperature of firn, which can promote additional surface melt in successive events. Meltwater production and refreezing also deplete the firn air content, which has preceded hydrofracture-induced ice-shelf breakup on the Antarctic Peninsula. Changes in the firn air content associated with melting and strain-induced thinning have poorly constrained impacts on the future integrity and susceptibility of the near-surface to hydrofracture and deserve further study.

Another important implication of this work is the key role of processes that are poorly captured in general circulation models, which do not reproduce the recent increasing trend in geopotential heights in the Amundsen and Bellingshausen seas and tend to underestimate extreme event frequency and blocking frequency in the broader southern Pacific region. This contributes to inaccuracies in the simulation of West Antarctic snowfall in general circulation models. These models can also have biases associated with meridional integrated water vapor transport that underestimate snow accumulation in polar regions (Espinoza et al., 2018; Guan and Waliser, 2017). In agreement with Emanuelsson et al. (2018), our study suggests that in order to accurately model snowfall over West Antarctica and, in particular, the surface mass balance of Thwaites Glacier, atmospheric blocking, moisture transport, and precipitation processes must be properly simulated.

6 Conclusions

Spatially sparse in situ meteorological observations of snowfall in Antarctica are an impediment to assessing mechanisms responsible for extreme snowfall in the ASE and the quality of snow accumulation simulated in climate reanalysis products. GNSS-IR inverse methods quantify surface mass balance and should be applied more broadly to GNSS networks in Greenland and Antarctica. These new techniques are advantageous compared to conventional automatic weather stations (AWSs) in several ways: (1) they require less equipment and power than AWSs; (2) the observations represent a larger surface area and are less likely to be biased by the station installation itself; (3) these measurements can be made on GNSS sites that are deployed for other reasons, such as ice dynamic and solid Earth studies; and (4) there exist a wealth of historical GNSS data that could be reprocessed to produce accumulation time series from GNSS-IR. In this study, we use GNSS-IR to derive accumulation records for the long-term on-ice GNSS sites located in the ASE along with two austral summers of campaign GNSS data. The GNSS-IR accumulation records reveal that precipitation on Thwaites Glacier is most frequent and intense in the austral winter (30 % to 50 % higher), with a peak accumulation in August that tapers into the austral summer to a minimum in January. The seasonal cycle is driven primarily by changes in the frequency of extreme precipitation events, which make up over half of total precipitation in the Amundsen Sea.

We use the GNSS-IR snow accumulation time series to investigate the synoptic atmospheric conditions that promote vapor transport and extreme precipitation in the Amundsen Sea Embayment. We find that blocking plays an important role in driving extreme precipitation. The connection between seasonal blocking in the southeastern Pacific and precipitation events in West Antarctica has important implications for the mass balance of vulnerable outlet glaciers located in the Amundsen Sea Embayment. We find that accumulation in the Amundsen Sea varies seasonally and interannually and that much of this change can be explained by the presence of blocking high-pressure systems that promote extratropical cyclonic activity and vapor transport to the Amundsen Sea. This suggests that atmospheric blocking, moisture transport, and precipitation processes must be properly simulated in order to accurately model snowfall and surface mass balance over West Antarctica.

Appendix A: Atmospheric catalogs

A1 Blocking catalog

Atmospheric blocking is associated with stable shifts in the position of persistent high-pressure systems, which can impact the track of cyclonic storm systems. Blocking describes many different states of the atmosphere, and thus many different methods have been developed to identify and track blocking events from reanalysis data (Schwierz et al., 2004; Barriopedro et al., 2010). We use the absolute geopotential height (AGP) reversal-based blocking index to identify blocking highs and track their evolution between 40 and 90° S, with a spatial overlap of a least 5 d based on the reversal of the geopotential height gradient at 500 hPa. The index was first implemented by Steinfeld and Pfahl (2019) and is a two-dimensional extension of the classical blocking index developed by Tibaldi and Molteni (1990), which follows from the description of Rex (1950).

A2 Atmospheric river catalog

We also catalog the atmospheric rivers that made landfall in the ASE between 1979 and 2021. Atmospheric rivers were identified using vertical integrals of specific humidity and horizontal wind from ERA5 to calculate total integrated vapor transport (IVT):

$$\begin{array}{}\text{(A1)}& \text{IVT}={\displaystyle \frac{\mathrm{1}}{g}}\underset{\mathrm{surface}}{\overset{\mathrm{trop}}{\int }}qvdp+{\displaystyle \frac{\mathrm{1}}{g}}\underset{\mathrm{surface}}{\overset{\mathrm{trop}}{\int }}qudp,\end{array}$$

where q is specific humidity, u is the zonal horizontal wind component, v is the meridional horizontal wind component, p is pressure, and g is the terrestrial gravitational constant. A catalog of events was created using a threshold value defined monthly for each grid cell. This detection scheme first identified the 98th percentile of total IVT over a 5-month period centered on the month to be cataloged for each cell. The landfall location was then identified based on the intersection of the IVT anomaly with the Antarctic coastline. The axis of the IVT anomaly was calculated following Brands et al. (2017) and then was used to determine the length and aspect ratio of the anomaly along with the mean IVT direction and the mean IVT magnitude. Events were then filtered based on their geometry and the magnitude of the meridional component of IVT relative to a poleward threshold value. IVT Anomalies are considered atmospheric rivers when the meridional to zonal aspect ratio of the IVT anomaly is greater than 2:1 and the meridional component of the IVT anomaly is larger than 100 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$. The detection algorithm was applied to hourly fields of vertically integrated vapor transport from 1000 to 300 hPa pressure levels to create an atmospheric river time series from 2008 to 2022 when continuous GNSS data on Thwaites Glacier were also available. Detected atmospheric river events presented here agree well with other catalogs made using similar methods (Wille et al., 2021).