Seasonal variability in Antarctic ice shelf velocities forced by sea surface height variations

Antarctica's ice shelves resist the flow of grounded ice towards the ocean through “buttressing” arising from their contact with ice rises, rumples, and lateral margins. Ice shelf thinning and retreat reduce buttressing, leading to increased delivery of mass to the ocean that adds to global sea level. Ice shelf response to large annual cycles in atmospheric and oceanic processes provides opportunities to study the dynamics of both ice shelves and the buttressed grounded ice. Here, we explore whether seasonal variability of sea surface height (SSH) can explain observed seasonal variability of ice velocity. We investigate this hypothesis using several time series of ice velocity from the Ross Ice Shelf (RIS), satellite-based estimates of SSH seaward of the RIS front, ocean models of SSH under and near RIS, and a viscous ice sheet model. The observed annual changes in RIS velocity are of the order of 1–10 m a⁻¹ (roughly 1 % of mean flow). The ice sheet model, forced by the observed and modelled range of SSH of about 10 cm, reproduces the observed velocity changes when sufficiently large basal drag changes near the grounding line are parameterised. The model response is dominated by grounding line migration but with a significant contribution from SSH-induced tilt of the ice shelf. We expect that climate-driven changes in the seasonal cycles of winds and upper-ocean summer warming will modify the seasonal response of ice shelves to SSH and that nonlinear responses of the ice sheet will affect the longer trend in ice sheet response and its potential sea-level rise contribution.

1 Introduction

The Antarctic Ice Sheet discharges mass via outlet glaciers and ice streams flowing into the ocean across the grounding lines, forming ice shelves several hundreds of metres thick, surrounding about half of the Antarctic coastline (Allison et al., 2011; Fretwell et al., 2013). Ice shelves play critical roles in ice sheet dynamics by providing back stresses that impede the gravity-forced flow of grounded ice towards the grounding line (Thomas, 1979). Ice shelf extent, thickness, and mass can vary over time (e.g. Cook and Vaughan, 2010; Paolo et al., 2015; Adusumilli et al., 2020), leading to changes in ice velocity for both grounded and floating ice (e.g. Scambos et al., 2004; Fürst et al., 2016; Reese et al., 2018; Gudmundsson et al., 2019). Persistent ice shelf thinning or retreat over years or decades can lead to a significant increase in the rate of mass loss of grounded ice (e.g. Velicogna et al., 2014; Joughin et al., 2014; Gudmundsson et al., 2019; Smith et al., 2020) and an associated increase in the rate of Antarctica's contribution to global sea level.

Time series of ice velocity (u_ice) from Global Navigation Satellite System (GNSS) receivers mounted on grounded and floating ice are, typically, of fairly short duration, limited to ∼1–3 months over austral summer. These short records reveal a strong tidal-band signal (e.g. Makinson et al., 2012) but cannot resolve annual cycles. However, a few longer GNSS records (e.g. Klein et al., 2020) and satellite-based estimates of u_ice (Greene et al., 2018, 2020) show variability on intra-annual (monthly to seasonal) timescales. Given that the seasonal cycle dominates variability in atmospheric and oceanic forcing of ice shelves, understanding how this forcing cycle affects ice shelf flow may provide important insights into the processes affecting the ice shelves and ice sheets and how they might respond to the weaker but more persistent forcing at longer timescales, from interannual variability (e.g. Dutrieux et al., 2014; Paolo et al., 2018) to multi-decadal trends (Jenkins et al., 2018).

Two mechanisms have been proposed to explain seasonal variability in ice shelf flow, linked to seasonal variability in (i) basal melt rates and (ii) sea ice. Klein et al. (2020) investigated the hypothesis that a seasonal cycle of spatially varying basal melt rates on the Ross Ice Shelf (Tinto et al., 2019; Stewart et al., 2019) might result in seasonality of u_ice; however, their modelled variability of u_ice was much smaller than GNSS measurements indicated. Greene et al. (2018) proposed that changes in buttressing from sea ice could explain the satellite-derived seasonal cycle of Totten Glacier's ice shelf; however, their uncertainties in satellite-derived intra-annual u_ice estimates were large, and the mechanism of ice shelf buttressing by sea ice is poorly understood.

In this paper, we investigate an alternative hypothesis: seasonal variability of sea surface height (SSH) modifies ice velocity through a combination of sea surface tilt and changing basal stresses at the grounding zone. This hypothesis is motivated by an extension of the role of tides in ice shelves and grounded-ice motion (Gudmundsson, 2007, 2013; Brunt and MacAyeal, 2014; Rosier and Gudmundsson, 2020), evidence from open ocean satellite altimetry that SSH around Antarctica has a pronounced seasonal cycle (Armitage et al., 2018; Rye et al., 2014) and the recent development of ocean models from which estimates of seasonal variability of SSH under ice shelves can be extracted. We explore our hypothesis by running a viscous model of the ice sheet and ice shelf in the Ross Sea sector with forcing from the modelled seasonal cycle of SSH under the Ross Ice Shelf and comparing the model output with GNSS time series of ice shelf velocity. We selected the Ross Ice Shelf because variability in ice shelf mass balance at longer timescales is known to be small (Das et al., 2020; Adusumilli et al., 2020), and there are several GNSS records exceeding 1 year in length that reveal intra-annual variability (e.g. Siegfried et al., 2014a; Bromirski et al., 2017; Blewitt et al., 2018). We show that the ice sheet model reproduces the observed annual cycle of the GNSS records if a sufficiently large cycle of SSH-induced basal shear stress change near the grounding line is parameterised in our viscous model. SSH-induced tilt of the ice shelf provides a small but significant additional contribution to velocity changes.

2 Data and models

We explore our hypothesis using a combination of in situ and satellite-derived observations, as well as ocean and ice sheet modelling. We take advantage of several existing GNSS records from the Ross Ice Shelf (RIS) collected during various field campaigns (Sect. 2.1), focusing on the ones that are sufficiently long to identify intra-annual velocity variations. We combine these records with estimates of intra-annual variations in SSH fields for the open ocean in front of the ice shelves from an existing satellite altimetry data set and from ocean models that include ice shelves (Sect. 2.2). We then compare the GNSS records to an ice flow model forced with the varying SSH from the ocean models (Sect. 2.3).

2.1 GNSS data

We use several long (5–19 months) time series of ice shelf motion from GNSS deployments on RIS (Fig. 1). These records were collected during different time intervals between 2014 and 2019 (Table 1). GNSS data from all stations were processed with a precise point positioning (PPP) approach (Zumberge et al., 1997; Geng et al., 2012, 2019).

DRRIS 2015–2016. An array of 13 GNSS stations was deployed on RIS from November 2015 to December 2016 as part of the Dynamic Response of the Ross Ice Shelf to Wave-induced Vibrations (DRRIS) project (Bromirski et al., 2017; Klein et al., 2020). Three stations were deployed along the ice front and nine along a flowline from the central ice front station to about 400 km upstream. One station (RS03) was located 100 km to the west of the along-flowline array and another (RS08) was on grounded ice on the western margin of Roosevelt Island. Only one station (DR10) recorded position data for a full year; however, the intra-annual signals in positions and velocities at the other DRRIS stations on floating ice were highly correlated with DR10 observations (Klein et al., 2020, their Fig. 6).

WISSARD 2014–2016. An array of GNSS stations was deployed as part of the Whillans Ice Stream Subglacial Access Research Drilling (WISSARD; Siegfried et al., 2014a; Tulaczyk et al., 2014) project. We used the record from station GZ19 located about 3 km offshore of the Whillans Ice Stream grounding line, which acquired data between November 2014 and November 2016 (Begeman et al., 2020).

Antarctica PI continuous network 2017–2019. Two GNSS stations (BATG and LORG) acquired data in the northwestern RIS. We obtained the time series for these sites from the GNSS database processed by the Nevada Geodetic Laboratory (NGL; Blewitt et al., 2018). Station BATG was located about 100 km east of Minna Bluff and acquired data from February 2017 to August 2018. Station LORG was located about 100 km east of Ross Island and about 90 km from BATG; the station recorded from November 2018 to November 2019 with a few interruptions, for a total of 289 d. The vertical components of tidal variability at these stations were reported by Ray et al. (2021).

Table 1 Station latitudes and longitudes at the time of deployment, mean speed, project and/or database which collected the data, duration (number of days of available data), and periods of deployment for GNSS stations. The primary stations used in this study are indicated in bold.

Figure 1 Map of the Ross Ice Shelf and its surrounding principal outlet glaciers and ice streams. The locations of GNSS stations used in this study and their names are indicated; see Table 1 for more details. Our focus is on long time series from DR10, BATG, LORG, and GZ19 (yellow stars). BYRD (orange square) is not a GNSS site but identifies the area analysed in Fig. 9e. The background image shows time-averaged surface velocities measured by satellites (Rignot et al., 2017). The grounding line and the ice front, from Depoorter et al. (2013), are plotted with black lines. The 1500 m isobath, separating regions defined as the open continental shelf (OCS) and the deep Ross Sea (DRS), is plotted in dark grey.

2.2 SSH measurements and model estimates

SSH can be estimated using satellite radar altimetry, and monthly SSH estimates are available for the period 2011–2016 for regions north of the Antarctic coastline and ice shelves using measurements from the European Space Agency's CryoSat-2 radar altimeter (Armitage et al., 2018). These SSH estimates cover fully open water (free of ice shelves) and leads in the ice pack but do not extend under the ice shelves. Measuring SSH variations in the ocean cavities under ice shelves is challenging because they are small compared with other contributors to height changes, such as uncertainties in seasonal cycles of basal mass balance (e.g. Stewart et al., 2019; Tinto et al., 2019), snow and firn density changes (e.g. Zwally and Jun, 2002; Arthern and Wingham, 1998), and penetration of radar signals into the surface snow and firn layers (Ridley and Partington, 1988; Davis and Moore, 1993). Therefore, it is not currently possible to accurately estimate SSH variability under ice shelves.

Instead, we investigated the representation of intra-annual variability in SSH from five existing ocean models with thermodynamically active ice shelves (Mathiot et al., 2017; Tinto et al., 2019; Naughten et al., 2018; Dinniman et al., 2020; Richter et al., 2022). We used their SSH output relative to the Armitage et al. (2018) open water data set to determine the most realistic model for analyses and to assess the likely variability in SSH under ice shelves. More information on these models, and assessment of their performance, is provided in the Supplement. From these analyses we determined that the Ross Sea regional model described by Tinto et al. (2019) provides a seasonal cycle that is most consistent with the Armitage et al. (2018) satellite-based results for the Ross Sea continental shelf north of RIS, suggesting that it is also the best model for SSH variability under RIS.

2.3 Ice sheet/ice shelf model

2.3.1 Model summary and initialisation

We used the open-source ice sheet and ice flow model Elmer/Ice (Gagliardini et al., 2013), the glaciological extension of the Elmer finite-element software developed at the Center for Science in Finland (CSC-IT). The modelling framework is similar to that described by Klein et al. (2020). We added variability in SSH in both time and space, relative to the initial static sea level, focusing on SSH output from the Tinto et al. (2019) ocean model as justified in the Supplement. Our ice model uses the vertically integrated shallow-shelf approximation (SSA; MacAyeal, 1989), a simplification of the Stokes equations (usually used for resolving viscous flow problems) in which the ice velocity is considered constant throughout the ice thickness. This approximation is well suited to ice shelves and ice streams where vertical shear stresses are negligible relative to other stresses acting on the ice. The ice rheology is based on a nonlinear constitutive relationship between strain rates and deviatoric stresses, classically used in ice flow modelling and known as Glen's flow law (Glen, 1958). The shear stress at the ice–bed interface, τ_b, is modelled with a Weertman friction law (Weertman, 1957) at the ice–bed interface:

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\mathrm{b}}=C{\left|u_{\mathrm{b}}\right|}^{\frac{\mathrm{1}}{m}-\mathrm{1}}{u}_{\mathrm{b}},\end{array}$$

with C being the friction coefficient, u_b the sliding velocity, and exponent $m\in [\mathrm{1}-\mathrm{\infty }]$ where increasing values of m are characteristic of a more plastic bed. We use a value m=1 in this study and discuss this choice in Sect. 4.1.

Following the same procedure used by Klein et al. (2020), we initialised our model by inferring the basal shear stress (on grounded ice) and the ice viscosity, using an inverse model that optimises the two parameters by minimising the difference between the model and observed surface ice velocities as well as the difference between ice flux divergence and observed mass balance. Details on the model setup can be found in Appendix A1.

There are two main effects of SSH variability on the ice shelf velocities: (i) changes in driving stress and (ii) changes in basal stress through grounding line migration.

Figure 2 Conceptual model of the SSH effect on the ice shelf slope and grounding line position: combination of (a) a positive ice shelf tilt and a negative ΔSSH close to the grounding line and (b) a negative ice shelf tilt and a positive ΔSSH close to the grounding line. The average annual state of the ice shelf is shown by dashed lines, while the perturbed state is shown by plain lines. The combinations shown are for seasonal velocity changes from grounding line migration and tilt being roughly in phase as suggested by SSH models for RIS (see Fig. 4).

• i.

  Driving stress change.

  Changes in gradients of SSH locally impact the driving stress, σ_g (in MPa), acting on the ice flow. This stress is a direct function of the surface gradient, ∇z_s, with z_s being the ice shelf surface height (assuming solid ice from the surface to base) relative to the background unperturbed sea surface, following, for example, Morland (1987), MacAyeal (1989), and Gudmundsson (2013):

  $$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_{\mathrm{g}}={\mathit{\rho }}_{\mathrm{ice}}gh\mathrm{\nabla }(z_{\mathrm{s}}+\mathrm{\Delta }\text{SSH}).\end{array}$$

  In Eq. (2), ρ_ice is the density of ice (917 kg m⁻³, assumed constant over the ice thickness), g is the gravitational acceleration (9.81 m s⁻²), h(xyt) (m) is the ice shelf thickness, and ΔSSH(xyt) is the SSH perturbation. A decrease in the ice shelf seaward gradient leads to a decrease in driving stress and a deceleration in the ice flow (Fig. 2a). An increase in the ice shelf seaward gradient leads to an increase in driving stress and an acceleration in the ice flow (Fig. 2b).

• ii.

  Change in basal stress through grounding line migration.

  SSH variations lead to changes in bed stresses in the grounding zone, as they raise and lower the ice shelf and influence the subglacial hydrology near the grounding line. A negative ΔSSH at the grounding line causes a downstream migration of the grounding line, increasing the grounded-ice area and potentially slowing down ice movement through an increase in basal drag (Fig. 2a). Conversely, a positive ΔSSH at the grounding line leads to an upstream migration of the grounding line, decreasing the area affected by basal stresses and accelerating the ice flow (Fig. 2b). The grounding line migration distance (ΔL) upstream and downstream is influenced by viscoelastic deformation of the ice shelf. The mechanism has been studied in the context of tidal deformation by treating it as an elastic and hydrostatic beam problem (e.g. Sayag and Worster, 2011, 2013; Walker et al., 2013). This analytical solution agrees reasonably well with grounding line migration calculated by solving the contact problem in a viscoelastic, tide-forced model (Rosier et al., 2014).

  In a purely hydrostatic framework, the grounding line migration (ΔL) depends on both the surface and bed slopes (Eqs. B7 and B8; Appendix B) (Sayag and Worster, 2014): as surface and bed slopes decrease, ΔL increases. This inverse relationship directly affects the magnitude of the change in friction in the grounding zone and also the ice flow response. The implications of uncertainties in our knowledge of bed slope and ΔL and the mechanical processes involved are discussed further in Sect. 4.2.

  Grounding line migration ΔL has also been treated as an elastic fracture problem, accounting for water pressure variations at the ice base as the grounding line migrates. Using this framework, Tsai and Gudmundsson (2015) showed that the magnitude of upstream ΔL is larger than in the hydrostatic or purely elastic case and depends nonlinearly on parameters such as the ice thickness and ΔSSH. For thick ice (e.g. in the grounding zone of Byrd Glacier), ΔL can be more than twice the value obtained using the hydrostatic framework, and for small ΔSSH (typically, a few centimetres), ΔL can be as much as 1 order of magnitude higher than in the hydrostatic framework.

2.3.2 Model runs

We ran 100 inversions of both the basal friction and the ice viscosity, constraining the fit to velocity and thickness rate of change observations, as well as the degree of smoothness of the solution. The set of inversions explores the effect of each constraint by varying their respective weight. From this ensemble of initial states, we selected an optimal (in terms of velocity and ice flux divergence fit) sub-ensemble of 15 members (Ω₁₅). The details of the initialisation procedure and the selection of Ω₁₅ are discussed in Appendix A1.

Using the sub-ensemble Ω₁₅ as a reference, we applied monthly averaged SSH anomalies (ΔSSH) from five different ocean models (see Supplement) as a steady-state perturbation, raising or lowering the ice surface, and ice base and computing the flow change with respect to the reference (see Appendix A2). For each run, we kept the ice shelf thickness $h(x,y,t)$ constant and assumed that the ice shelf and the grounding line location adjust instantaneously to the ΔSSH.

To assess the importance of methods for representing grounding line migration in our viscous ice sheet model, we ran three different parameterisations of ΔL (for a total of 225 simulations =15 members × 5 SSH models × 3 grounding line parameterisations), as follows:

1. ΔL_B2 is based on the hydrostatic equilibrium of the grounding line and Bedmap2 (a gridded product describing surface elevation, ice thickness, and the basal topography of the Antarctic; Fretwell et al., 2013) bed slopes at the grounding line.

2. ΔL_C (constant bed slope) is a significantly larger migration that corresponds to values used by Rosier and Gudmundsson (2020) for their study of the Filchner–Ronne Ice Shelf when treating the grounding line migration with elastic fracture mechanics introduced by Tsai and Gudmundsson (2015).

3. ΔL_B2L is also a larger value of grounding line migration but accounting for the Bedmap2 surface and bed slope variations along the grounding line.

To account for subgrid-scale migration of the grounding line, our model implementations parameterise ΔL as a change in friction, rather than as a change in floatation state at specific grid nodes (Appendix B). The implications of the two larger migration parameterisations are discussed in more detail in Sect. 4.2.

3 Results and discussion

We first review the intra-annual variability in ice flow recorded by the GNSS receivers on RIS (Sect. 3.1.1) and the measured (Armitage et al., 2018) and modelled (Tinto et al., 2019) seasonal cycles of SSH for the Ross Sea including under RIS (Sect. 3.1.2). We then compare the variability in driving stresses due to SSH anomalies and grounding line migration (Sect. 3.2) and the effect of both processes on the ice speed flow (Sect. 3.3).

3.1 Intra-annual signals in GNSS displacement and SSH records

3.1.1 GNSS displacement

All long-duration GNSS stations on RIS (Sect. 2.1) show variability in horizontal displacement on various timescales including diurnal (∼1 d period), fortnightly (∼2-week period), and intra-annual (Fig. 3). As reported by Klein et al. (2020), data from the DRRIS stations show evidence of an annual cycle with a displacement anomaly amplitude of about 1 m, alternating between a negative trend during December–May and a positive trend during June–November. GZ19 shows no apparent annual cycle, but its displacement shows a similar range of variability (±1 m) to DR10 during the 2-year record. The time series at BATG, which is not concurrent with the DRRIS stations and GZ19, shows a smaller amplitude range (about 0.2–0.3 m) that appears to have a periodicity of about 6 months. The LORG time series in 2019 shows a similar pattern to BATG in 2018.

Figure 3 GNSS horizontal displacement anomalies in the north direction (approximately parallel to the time-averaged flow) for GNSS stations used in this study. The time interval for each panel is 2 years; however, years differ between panels. (a) DR02, DR04, DR10, DR16, and RS16 (for legibility, other DRRIS sites are not shown here but exhibit a similar trend; the complete array can be found in Klein et al., 2020); (b) GZ19; (c) BATG; and (d) LORG. Note that (a) and (b) are plotted on the same timescale, while (c) and (d) have 2-year and 3-year shifts with respect to the two upper panels. The black lines are smooth versions of displacement anomalies with a 1 d Gaussian RMS width.

The diurnal lateral displacement signal is caused by the fundamental tides of the region, which are almost entirely diurnal (e.g. Padman et al., 2003; Ray et al., 2021). We attribute the fortnightly signal in displacement at all GNSS sites (and, possibly, also the ∼6-month periodicity at BATG and LORG) to nonlinear response of the ice sheet and ice shelf to variability in the tidal range, leading to viscoelastic flexural adjustments of the ice sheet at the grounding zone, as the range of the diurnal tide varies through the fortnightly spring–neap modulation (e.g. Rosier and Gudmundsson, 2020). We removed the fortnightly tide-forced variability by filtering to monthly and longer timescales (by using a sliding Gaussian filter with a 2-week standard deviation); however, any ∼6-month tidal signal remains as a source of noise in our interpretation of intra-annual ice shelf flow changes driven by non-tidal SSH variability.

Figure 4 Seasonal sea surface height deviation from the annual mean (ΔSSH): (top row) satellite observations averaged over the period 2011–2016 (ΔSSH_CS2, Armitage et al., 2018) and (bottom row) modelled for the period 2002 (SSH₂₀₀₂, Tinto et al., 2019). The ice front and grounding line are represented by black lines. The outer edge of the open continental shelf (OCS) is along the 1500 m isobath, shown with a grey line. Ice speeds are shown in shades of grey, with darker shades being faster.

3.1.2 Satellite-derived and modelled SSH

The seasonal cycle of ΔSSH in satellite-derived SSH fields around Antarctica, for the period 2011–2016, shows a typical range of about 5 cm on the open continental shelf (OCS; see Fig. 1) of the Ross Sea and comparable changes offshore in the Deep Ross Sea (DRS); see Fig. 4, top row, and Fig. 5. For the OCS, a positive SSH anomaly occurs in winter (April–September). The Tinto et al. (2019) model, based on annually repeating forcing for 2002, shows similar phasing of the ΔSSH cycle (Fig. 4, bottom row; Fig. 5a) but with larger amplitude than for the satellite-derived fields. The qualitative and quantitative (see Table S2 of Pearson correlation coefficients of the different models) agreement between the model and the observations offshore of RIS provides support for the use of this ocean model for predicting SSH variability under RIS, even though the ocean model does not overlap in time with either the observed SSH fields or the GNSS observations.

Figure 5 (a) Annual cycle of monthly mean ΔSSH over the open continental shelf (OCS – plain lines) and beneath the ice shelf (RIS – dotted lines) for SSH₂₀₀₂ (blue) and for CryoSat-2 measurements (SSH_CS2, red) averaged over 2011–2016, for the open continental shelf (OCS) only. (b) Mean ΔSSH for the deep Ross Sea. The grey shade shows the winter period. See Fig. S3 for similar comparisons that include all available ocean models of SSH.

3.2 Comparing driving stress change and grounding line migration

RIS thickness decreases from ∼800 m close to the grounding line to ∼300–400 m at the ice front, over a distance of ∼800 km (see Tinto et al., 2019, their Fig. S2a). This results in mean thickness and surface gradients of about $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$, respectively. Since we are interested primarily in along-flow variations in ice velocity, we calculate the along-flow Lie derivatives (Yano, 2020) of the ice shelf surface height (${\mathrm{\nabla }}_{\widehat{u}}{z}_{\mathrm{s}}$) and ΔSSH (${\mathrm{\nabla }}_{\widehat{u}}\text{SSH}$). Values of ${\mathrm{\nabla }}_{\widehat{u}}{z}_{\mathrm{s}}$ range from 10⁻⁵ to 10⁻² over most of the ice shelf (Figs. 6a and S4a). Gradients of ΔSSH in Tinto et al. (2019; SSH₂₀₀₂) can reach 10⁻⁶ to 10⁻⁵ in February (Figs. 6b and S4b). This means that local tilting of the ice shelf by ${\mathrm{\nabla }}_{\widehat{u}}\text{SSH}$ can modify the local driving stress of the ice shelf (Eq. 2) typically by 0.1 %–1 % and sometimes up to several percent, with substantial spatial variability (Fig. 6c). ${\mathrm{\nabla }}_{\widehat{u}}\text{SSH}$ also varies by month (not shown). For example, in February, about 30 % and 6 % of the ice shelf experiences a fractional change of driving stress exceeding 0.1 % and 1 %, respectively (Fig. 6d). The largest fractional change in driving stress occurs away from the grounding line where the ice surface height gradients are smaller than closer to the grounding line and where the SSH gradients are the larger.

The complex spatial variability of the along-flow derivatives of ΔSSH (Fig. 6b) arises from changes in orientation and magnitude of the sub-ice-shelf circulation relative to ice flow. This circulation is itself complex: see, for example, Supplementary Video 1 in Tinto et al. (2019).

For most months there is a strong along-flow gradient in SSH close to the ice front (Figs. 4b and 6b), which directly impacts driving stress (Eq. 2). These variations in driving stress lead to ice velocity changes, which we present as anomalies with respect to the annual average velocity field. In general, months with a regionally averaged (i.e. over the ice shelf) negative ΔSSH (e.g. January–March period in Figs. 4 and S2) that slows ice flow as the grounding line migrates seaward also experience a relative uplift of the surface close to the ice front, leading to an additional slowdown (Fig. 2a). Conversely, the months experiencing a regionally averaged positive ΔSSH generally show a relative surface drop close to the ice front and an upstream migration of the grounding line, both contributing to an acceleration of the ice shelf (Fig. 2b).

Figure 6 Comparison of ice shelf surface gradients and SSH gradients from SSH₂₀₀₂ (Tinto et al., 2019), both calculated in the direction of ice flow ($\widehat{u}$) in February. (a) Ice shelf surface gradient (${\mathrm{\nabla }}_{\widehat{u}}{z}_{\mathrm{s}}$), (b) SSH gradient (${\mathrm{\nabla }}_{\widehat{u}}\text{SSH}$), and (c) their ratio. Gradients are filtered with a 5 km standard deviation Gaussian smoothing. (d) Gradient values, for each 1×1 km cell, plotted as a function of each other. The colour map represents the ice flow speed. A total of 6 % and 30 % of the model nodes over the ice shelf experience a driving stress variation of more than 1 % (left of the plain line) and 0.1 % (left of the dashed line), respectively.

Figure 7 February anomaly in velocity ΔU, averaged over an ensemble of 15 initial states (Ω₁₅), formed as the difference between the annual average for ΔSSH₂₀₀₂ and the three different parameterisations of the grounding line migration: (a) ΔL_B2, (b) ΔL_C, and (c) ΔL_B2L (see Sect. 2.3.2). The locations of DR10, GZ19, BATG, and LORG (identified in Fig. 1) are indicated by white stars. The grounding line and the ice front are shown by black lines. The background annual average flow velocity for grounded ice is plotted in shaded grey, with darker grey being faster.

Our modelling indicates that the amplitude of the grounding line migration, ΔL, is the primary control on the amplitude of the seasonal velocity signal. In February, for example, the model ensemble using ΔL_B2predicts the smallest amplitude of velocity deviation of the three cases, with $\mathrm{\Delta }{U}_{\mathrm{B}\mathrm{2}}\sim -\mathrm{1}$ m a⁻¹ over most of the ice shelf (Fig. 7a). Larger values of ΔL (parameterisations ΔL_C and ΔL_B2L) allow the grounding line to move farther downstream during summer, leading to deviations $\mathrm{\Delta }U\sim -\mathrm{3}$ m a⁻¹ in the centre of the ice shelf (Fig. 7b and c). The largest differences between the effects of ΔL_C and ΔL_B2L are generally found close to the grounding line in the deep and narrow fjords such as the floating extension of Byrd Glacier where ΔL_C leads to a slowdown ΔU_C<5 m a⁻¹ compared with ΔU_B2L∼3 m a⁻¹ (Fig. 7b and c). These are regions where true bed slopes are steeper than the average around the RIS perimeter and which are also more sensitive to the initial state as the ensembles show a larger standard deviation in these areas with respect to the rest of the domain (Fig. 8, bottom row).

We regard the ΔL_B2 parameterisation, which yields small grounding line migration, as an approximation of ice shelf response to SSH gradients alone.

3.3 Seasonal cycle in ice flow

All ensembles forced with ΔSSH₂₀₀₂ exhibit a maximal seasonal negative flow speed anomaly during summer and maximal positive anomaly during winter (Fig. 8); however, ΔL_B2 simulations tend to switch to positive anomalies later than simulations using ΔL_C and ΔL_B2L. Simulations using ΔL_B2 produce maximal amplitudes of speed anomaly at the ice front that progressively decrease farther upstream, while ΔL_C and ΔL_B2L produce maximal speed anomaly amplitudes in the deep fjords along the base of the Transantarctic Mountains. The amplitudes of speed anomalies of ΔL_B2 are about 2–4 times smaller than for ΔL_C and ΔL_B2L simulations, depending on location.

To validate the results of the three grounding zone parameterisations, we extracted the modelled ice velocity anomalies at the GNSS locations and compared these to velocity variations (Fig. 9) estimated from the time derivative of measured displacement anomalies (Fig. 3).

At DR10, the range of the observed velocity anomaly (ΔU) was about 10 m a⁻¹ with a minimum in February–March and a maximum in July (Fig. 9a). The other DRRIS GNSS stations located in the centre of the ice shelf did not record during austral winter (see Fig. 5a), preventing us from properly identifying the timing of maximum velocity for these stations. The ΔL_C and ΔL_B2L ensembles both give similar ΔU estimates that are qualitatively similar to observations, with velocity variations of about 50 % to 70 % of the observed amplitude and minima and maxima in summer and winter, respectively. The ΔL_B2 grounding zone parameterisation has a much lower amplitude than observations and gives a maximum velocity in October, about 2 months later than the other ensembles and 4 months later than the observations. However, the timing of the summer ΔU minimum is close to the observations and the other grounding zone parameterisations. Expanding our analysis to the entire GNSS array of DRRIS, similar seasonal phasing occurred at each GNSS station located approximately along the central flowline of the ice shelf. ΔU amplitude generally decreases with increasing distance from the ice front (Fig. 10), although with some variability that may result from the proximity of the DRRIS array to the Byrd Glacier flow and its impact on RIS flow.

At GZ19, close to the grounding line of Whillans Ice Stream, there is no seasonal cycle visible in the GNSS observations of displacement anomaly (Fig. 3b). The measured velocity anomaly (Fig. 9b) shows an overall slowdown, consistent with previous observations of slowdowns of Whillans and Mercer ice streams and the adjacent region of RIS over the last few decades (e.g. Joughin et al., 2005; Thomas et al., 2013), and shorter periods of deceleration and acceleration that could be due to the inherent variability in the two ice streams (e.g. Winberry et al., 2009). This trend was not captured by our ice flow models, which do not account for varying forcing other than the annual cycle of SSH. The modelled anomalies at GZ19 are weak, with $\mathrm{\Delta }{U}_{\mathrm{B}\mathrm{2}}\sim \pm \mathrm{0.5}$ m a⁻¹ and $\mathrm{\Delta }{U}_{\mathrm{C}}\sim \mathrm{\Delta }{U}_{\mathrm{B}\mathrm{2}\mathrm{L}}\sim \pm \mathrm{1}$ m a⁻¹ over the year.

At station BATG, about 100 km east of Minna Bluff, the velocity time series shows an approximately 6-month periodicity, with a ΔU range of about 2.5 to 3 m a⁻¹ (Fig. 9c). ΔL_B2 provides a poor fit to these observations, in both ΔU amplitude and phase, with the amplitude better reproduced by ΔL_C and ΔL_B2L. However, the pattern of observed velocity anomaly changes between the first and second year of the record. In the first year, the 6-month cycle shows a large velocity drop in July–August (reaching a minimum in September), corresponding to the second minimum of the year. In the second year, the observed velocity reached a maximum in May and remained relatively high until the end of August, fitting the modelled velocities. While the record terminated at the end of August, this marked a particularly long plateau of high velocities (from May to August), suggesting that the record includes a seasonal signal that is added to the 6-month cycle that we tentatively attribute to semiannual changes in tidal range (see Sect. 1).

Figure 8 Ensemble mean seasonal (3-month average) ice flow anomaly for ΔSSH₂₀₀₂ and three parameterisations of the grounding line migration: (first column) Bedmap2 (ΔL_B2), (second column) a constant bed slope (ΔL_C), and (third column) a flatter version of Bedmap2 (ΔL_B2L). The seasonal anomalies are computed from monthly model outputs. The standard deviation over the ensembles (bottom row) shows variability in space and time over the year. The locations of DR10, GZ19, BATG, and LORG are indicated by white stars. The ice front and the grounding line are indicated by the black line. Ice surface velocities over the grounded ice are plotted with a grey scale, from white (slow flow) to dark grey (fast flow).

Figure 9 Comparison between GNSS and model velocity anomaly when applying ΔSSH values from SSH₂₀₀₂ for (a) DR10, (b) GZ19, (c) BATG, and (d) LORG and (e) at the Byrd Glacier outlet (see locations in Fig. 1). The annual model cycle is repeated over 2 years. The average model velocity anomalies (over Ω₁₅ ensembles) – ΔU_B2 (grey), ΔU_C (red), and ΔU_B2L (blue) – are displayed with 1 and 2 standard deviations of the 15 estimates in each ensemble (dark and light shades, respectively). If not visible, the standard deviation is statistically insignificant. ΔU_B2, ΔU_C, and ΔU_B2L are plotted from monthly time series with a continuous line drawn between each snapshot. In (b), ΔU_C (red) and ΔU_B2L (blue) are so similar that we cannot distinguish them. The observed velocities (green) are obtained as the time derivative of the measured displacement anomaly (the period of observation is given in green in each panel) from GNSS, with a Gaussian filter with a 2-week standard deviation. See Fig. S6 for similar comparisons that include all available ocean models of SSH.

The time series of ΔU at LORG (Fig. 9d) for the period November 2018 to November 2019 is highly correlated (p=0.95) with the time series at BATG over the second year (from November 2017 to October 2018). The predicted velocity anomalies for ΔL_C and ΔL_B2L at these two stations agree especially well with the observations over the entire LORG times series and the second year of the BATG time series. More specifically, the model is able to reproduce the month-to-month accelerations and decelerations and the overall longer span of positive anomalies visible in LORG observations.

To examine the relative effect of the variations in driving stress and basal friction through grounding line migration, we consider a key region of RIS, the floating extension of Byrd Glacier near its grounding line. Byrd Glacier is the fastest and the deepest outlet glacier feeding RIS and is the region of RIS where the outputs from the three ensembles (ΔL_B2, ΔL_C, and ΔL_B2L) deviate the most. Observations show that, over a time span of a few years, flow upstream of the grounding line can increase by 10 %, coinciding with the discharge of subglacial lakes lubricating the bed (Stearns et al., 2008; Yuan et al., 2023). At seasonal timescales, variations in Byrd Glacier remain poorly constrained due to the lack of year-round GNSS measurements; however, Greene et al. (2020) used feature tracking in satellite imagery to estimate ice velocities and characterise the magnitude and timing of seasonal ice dynamic variability. For a region close to the grounding line of Byrd Glacier, they estimated seasonal variability of ΔU with a range of roughly 45 m a⁻¹. Our ensemble using the ΔL_B2L representation (Fig. 9e) shows a phase that is consistent with Greene et al. (2020); however, our modelled range in ΔU is always less than 10 m a⁻¹. This difference could be explained by the substantial uncertainty due to irregular, seasonally biased sampling of the satellite data (see Fig. 4 of Greene et al., 2020). Our modelling may also underestimate either the ΔSSH or the basal condition changes that the ΔSSH changes trigger. There may also be other processes at play that we do not account for in our modelling; for example, the change in seasonal melt (explored in Klein et al., 2020) which, while expected to be small, could slightly increase the ΔU signal.

4 Sources of uncertainty in SSH and ice flow response

Our ice sheet modelling results suggest that seasonal variations in SSH beneath RIS are sufficient to drive ice velocity variations of several metres per year over a large portion of the ice shelf when using the ΔL_C and ΔL_B2L parameterisations to represent basal stress change in the migrating grounding zone. The modelled velocity variability generally decreases with increasing distance from the ice front, although large variability is also associated with several major outlet glaciers flowing through the Transantarctic Mountains. However, the correlation between model and GNSS observations depends on the model initialisation, friction law, grounding line parameterisation, and the source of the SSH forcing. In this section, we discuss the sensitivity of the model and the uncertainty of each of these parameters.

4.1 Model initialisation and friction law

The inverse model used to generate the initial steady-state solution is under-constrained. Because we infer two parameters with multiple constraints during the inversion, an initial state with a minimal velocity misfit will not necessarily lead to a minimal ice thickness rate of change. Different combinations of friction and viscosity parameters can lead to similar misfits. Using the ensemble Ω₁₅, consisting of the 15 optimal initial states (see Sect. 3.3 and Appendix A1), to estimate the effect of the initialisation on the forward model helps quantify this effect. For the ensemble of simulations using the ΔL_B2 parameterisation, the impact of the initial state on the velocity intra-annual cycle is minimal over the ice shelf, with an average relative standard deviation under ∼5 % over most of the ice shelf (Fig. 8 (bottom row) and Figs. 9 and 10). The ensemble responses for the ΔL_C and ΔL_B2L parameterisations, while providing more realistic estimates of intra-annual velocity changes, show more sensitivity to the initial state with year average relative standard deviations of ∼15 %–20 % (±0.1–0.15 m a⁻¹) at DR10 and ∼25 % (±0.4 m a⁻¹) at Byrd Glacier. We attribute the relatively high variance of the ensemble in these regions to the sensitivity of the model to the initial basal friction, while the relatively low variance of the ensemble over most of the ice shelf indicates low sensitivity of the model to the initial viscosity parameter.

Figure 10 Peak-to-peak seasonal range of velocity anomaly (ΔU) when forcing the ice flow model with ocean model ΔSSH₂₀₀₂. The error bars (in shades) correspond to 1 and 2 standard deviations in each ensemble. Observed peak-to-peak range is also plotted for GNSS stations with data records longer than 1 year (i.e. DR10, GZ19, BATG, and LORG).

The friction law used in the model will also influence ice flow response, even for the same value of ΔL. Friction laws of different complexity have been proposed in the literature (Weertman, 1957; Budd et al., 1979; Schoof, 2005; Tsai et al., 2015) and have been shown to have different impacts on grounding line dynamics (e.g. Brondex et al., 2019). In our study, we only used the most common friction law (Weertman, 1957).

The results described in Sect. 3.2 were obtained with a linear version (m=1) of Eq. (1); i.e. stress is proportional to velocity. We also tested the value m=3 (e.g. Brondex et al., 2019; Gudmundsson et al., 2019), which only changes modelled velocity anomalies by a few percent. More complex friction laws (e.g. Schoof, 2005; Tsai et al., 2015; Joughin et al., 2019) that include the impact of water pressure change at the ice base as the grounding line migrates could increase the amplitude of our intra-annual velocity variations. However, such friction laws introduce additional poorly constrained parameters (Gillet-Chaulet et al., 2016) and, therefore, are not considered in this study.

4.2 Grounding line migration and basal stress change

The parameterisation of ΔL directly controls the amplitude of the grounding line migration which, in turn, controls the change in the friction coefficient we apply at the grounding line (see Sect. 3.3 and Appendix B). ΔL_B2 leads to small migration of the grounding line (typically a few metres) so that most of the impact of SSH variability on the ice flow comes from changes in ΔSSH gradients. While driving stress variations from these SSH gradients and small grounding line migration (ΔL_B2) due to ΔSSH can slow down or accelerate the ice flow by about ±1 m a⁻¹ (Figs. 6 and 7), these modelled variations are only ∼20 % of the observed ΔU at DR10. Incorporating a larger grounding line migration in the model (ΔL_C and ΔL_B2L) gives values of ΔU that are consistent with our GNSS observations. Such grounding line migrations with respect to the hydrostatic case (ΔL_B2) are, arguably, too strong, but are in line with observations by Brunt et al. (2011) and the values used by Rosier and Gudmundsson (2020) on the Filchner–Ronne Ice Shelf. The surface and bed slope are key parameters of the grounding line migration parameterisation (Tsai and Gudmundsson, 2015; Appendix B). The bed slopes around the RIS perimeter, estimated by Brunt et al. (2011) by applying the hydrostatic assumption to observed migration of the inner margin of tidal ice flexure in repeat-track satellite altimetry, are likely to be biased low, based on the modelling of Tsai and Gudmundsson (2015). While the ΔL values given by ΔL_C and ΔL_B2L are in the upper range, they remain consistent with previous studies of tidal migration of the grounding line.

Another explanation for our need for a large migration of the grounding line is the relatively low value of the basal friction coefficients we inferred at the grounding line during the model initialisation. Our initialisation scheme relies on the optimisation of the friction coefficient C in Eq. (1). This friction law does not include a direct dependency on the effective pressure as a Coulomb law would (e.g. Brondex et al., 2019; Urruty et al., 2022). However, as C is determined by inversion, it includes a dependency on the effective pressure close to the grounding line and reduces the value of C at the grounding line to match observations of velocity and thickness rate of change (e.g. Urruty et al., 2022). The inferred value represents an average annual value of the friction coefficient. The distribution of the seasonal variation around this annual average cannot be exactly determined without proper knowledge of the subglacial hydraulic system, but one can assume that the variation could be larger than the variations we estimate through our hydrostatic parameterisation (i.e. a change in C directly proportional to ΔL). Seawater intrusion at the ice–bed interface and in sediments has been shown to have a high impact on the ice flow response (e.g. Robel et al., 2022). Subglacial models depending on subglacial water pressure decrease effective pressure significantly near the grounding line, leading to an increased sensitivity for a given power in the sliding law (e.g. Kazmierczak et al., 2022). Seawater intrusion could also be enhanced by a highly retrograde slope (e.g. Byrd Glacier; see Fig. S6). Retrograde bed slope will enhance both the migration of the grounding line and the intrusion of seawater in the subglacial hydrologic system. The consequences of ΔL on this effective pressure are difficult to determine, but incorporating such a mechanism in our modelling could lead to a larger impact of ΔSSH on the ice flow, even in the purely hydrostatic case.

A final explanation relates to the potential effect of ΔL on the subglacial water system. If the grounding line retreats elastically over a short period, before relaxing to a position closer to the hydrostatic equilibrium, then this short retreat could modify the subglacial water system over a long distance upstream the grounding line. The proper modelling of the water system is largely out of the scope of this study but could help validate our theories in the future.

4.3 Estimating the sea surface height anomalies

The SSH anomalies (ΔSSH) computed in the different ocean models (see Sect. 2.2 and Supplement) result from temporal variability in ocean currents driven by wind stress and lateral density gradients. However, these models do not account for the steric changes due to thermal and haline expansion and contraction or the ocean's response to atmospheric pressure variations. Both the ROMS and NEMO modelling frameworks (see Supplement) use the Boussinesq approximation based on the Navier–Stokes equations: the models conserve volume rather than mass and therefore do not properly account for steric changes. At the same time, variations in atmospheric pressure also lead to isostatic adjustments of the ocean (e.g. Goring and Pyne, 2003), while ice shelves have been shown to respond similarly (Padman et al., 2003). This effect, known as the “inverse barometer effect” (IBE), is not considered in the simulations used in this study. Combining the effect of Boussinesq SSH variations (ΔSSH_boussinesq), the steric effect (ΔSSH_steric), and the IBE (ΔSSH_IBE), we obtain the total ΔSSH monthly deviation:

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }\text{SSH}=\mathrm{\Delta }{\text{SSH}}_{\mathrm{boussinesq}}+\mathrm{\Delta }{\text{SSH}}_{\mathrm{steric}}+\mathrm{\Delta }{\text{SSH}}_{\mathrm{IBE}}.\end{array}$$

Some efforts were made in the 1990s to evaluate the effect of steric sea level due to thermal expansion, concluding that a globally uniform, time-dependent correction of sea level can correct a non-Boussinesq solution (e.g. Greatbatch, 1994). Mellor and Ezer (1995) showed that the seasonal variation in this term is about 1 cm in the Atlantic Ocean, which represents about 10 % of our modelled amplitude of SSH variation over the ice shelf. At the spatial scale of RIS, this correction is roughly spatially uniform and, therefore, would not modify the driving stresses over the ice shelf but could affect the grounding line migration.

Seasonal changes in surface air pressure take place over the Antarctic continent, resulting in a decrease in surface pressure (loss of atmospheric mass) from January to April and an increase in surface pressure (gain of atmospheric mass) from September to December (Parish and Bromwich, 1997). Since most of the ocean models we presented use ERA-Interim reanalysis (Dee et al., 2011) as an atmospheric forcing, we therefore use ERA-Interim surface pressure over RIS to estimate the IBE effect contribution to ΔSSH and its potential effect on the ice flow. ERA-Interim is an older product than the currently recommended ERA5-Land surface air pressure (Hersbach et al., 2020), but both give similar surface pressures over RIS for the period we study, which limits the uncertainty of the IBE effect.

We simulate the effect of IBE on ΔSSH following Eq. (2) and apply the full ΔSSH as a forcing to the ice flow model. Due to the smaller isostatic adjustment of ice shelves to ΔSSH_IBE close to the grounding line, we do not include its effect in the grounding line migration. The relative effect of the IBE on the seasonal ice flow is maximal at DR10 and BATG due to their relative proximity to the ocean. In contrast, GZ19 and the region of Byrd Glacier are less affected, since the IBE does not impact grounding line migration (Fig. 10). Overall, accounting for the IBE modifies the peak-to-peak amplitude of ice flow variations by up to ∼1.5 m a⁻¹ (Fig. 10) without significantly impacting the seasonal pattern and phase of the ice flow velocity change. We note that if the IBE was to have a significant impact on the grounding line migration on average, it would most likely increase the amplitude of the grounding line migration with a similar phase to the one we observe without IBE. On a 38-year record of IBE (Fig. S5) the negative inverse barometer signal observed from December to June would lead to downstream migration of the grounding line and a deceleration in the ice shelf flow. Conversely, the positive signal observed from July to November would lead to an upstream migration of the grounding line and an acceleration in the ice shelf flow.

4.4 Ice rheology and timescales

Our ice flow model uses the shallow-shelf approximation (SSA), a viscous rheology, which is well suited for studying long-timescale mechanisms involving ice creep (more than a few days). At the same time, two of our parameterisations of the grounding line migration are based on an elastic rheology, which is more appropriate for short-timescale mechanisms such as tidal effects (less than a few days). In reality, both rheologies are at play but either can sometimes be disregarded with respect to the other, depending on the Maxwell time:

$$\begin{array}{}\text{(4)}& t_{\mathrm{M}}={\displaystyle \frac{\mathit{\eta }}{E}},\end{array}$$

with E being Young's modulus and η the characteristic viscosity of ice. Using a value η∼40 MPa a⁻¹ (obtained from the inferred viscosity parameter and strain rates averaged over the ensemble Ω₁₅ for RIS) and $E={\mathrm{10}}^{\mathrm{3}}-{\mathrm{10}}^{\mathrm{4}}$ MPa (Cuffey and Patterson, 2010) gives t_M∼2 d to ∼2 weeks. Given the seasonal timescale of the variability under consideration in this paper, our viscous ice flow model should adequately represent the real viscoelastic rheology of ice. The elastic migration of the grounding line is, therefore, less representative of the actual viscoelastic rheology for the timescale changes we are observing (SSH anomalies remain relatively stable in periods shorter than a month). However, the elastic parameterisation has previously been successfully applied in a viscoelastic ice flow model studying ice flow response to fortnightly tidal forcing (Rosier and Gudmundsson, 2020). As mentioned in Sect. 4.2., the elastic parameterisation is also a proxy to simulate unrepresented mechanisms that might trigger SSH-induced basal stress change in the grounding zone. Moreover, although the use of an elastic rheology to study a viscous problem usually requires decreasing the effective Young's modulus of ice (which could decrease ΔL), Tsai and Gudmundsson (2015) suggest that their parameterisation of the grounding line migration may also apply to a purely viscous case. This could also explain why grounding line positions in Stokes models, which are not constrained to the hydrostatic approximation, are generally more sensitive than in SSA models such as the one used in this study (e.g. Pattyn et al., 2013).

5 Conclusions

We have used an ice sheet model to investigate our hypothesis that sea surface height (SSH) variations can explain observed seasonal variability in ice velocity measured with four GNSS records of roughly 1- to 2-year duration on the Ross Ice Shelf (RIS). The model was forced with monthly SSH fields obtained from ocean models that include thermodynamically active ice shelves. Varying SSH fields can affect ice flow through two processes: changing the driving stress by locally tilting the ice shelf and changing the basal condition in the grounding zone. In ocean models that include the Ross Ice Shelf, the two sources of ice shelf acceleration – surface SSH sloping downwards towards the ice front and positive SSH anomalies along the grounding zone (Fig. 2b) – are roughly in phase. We found that the ice sheet model is able to reproduce the approximate phasing and magnitude of measured seasonal changes in ice velocity, given appropriate parameterisation of induced changes of basal stresses in the grounding zone.

In our model, the changes in bed stress due to grounding line migration as SSH changes are based on a parameterisation of viscoelastic processes, but these may also be interpreted as poorly understood effects on the subglacial hydraulic system just upstream from the grounding line. When this parameterised migration and/or basal shear stress change is sufficiently large, the combination of varying driving stress and grounding zone friction produces seasonal responses that are consistent with the data records at the GNSS station locations (Fig. 9). Station DR10 in the central northern RIS experienced the largest annual cycle, about 1 % of the annual mean flow, while station GZ19, located close to the grounding line of Whillans Ice Stream, does not include a substantial seasonal cycle. Modelled intra-annual ice flow changes at two stations in the northwestern RIS, BATG and LORG, are smaller than at DR10 but still significant. There is some evidence in the data from these sites to confirm the predicted annual cycles (Fig. 9c and d); however, these data records also include substantial variability at ∼6-month periodicity that is not apparent in the modelled signal. We tentatively attribute the ∼6-month signal to the astronomically forced, semiannual variability in daily tidal height range that results in time-averaged changes in grounded ice flow-through viscoelastic processes (Rosier and Gudmundsson, 2020). However, in the absence of concurrent measurements of SSH variability near the grounding line, we cannot rule out the presence of an SSH forcing signal with ∼6-month periodicity that is not represented in the SSH forcing models. We note that ocean models with annually repeating forcing, from which SSH forcing can be obtained, vary widely in their estimates of seasonal variations (Fig. S2), while multi-year simulations with realistic forcing that varies on interannual timescales produce large year-to-year changes in SSH (Fig. S1).

The largest modelled seasonal cycle in RIS ice flow occurs in the inlet close to the Byrd Glacier grounding line (Figs. 8, 9e). There are no long-term GNSS records from this region to confirm the modelled values; however, a previous study using satellite-derived variations in ice flow for Byrd Glacier confirms that this region experiences large seasonal flow variability (Greene et al., 2020). The high amplitude of the modelled velocity anomaly in this region is determined by the bed geometry, the associated amplitude of the grounding line migration, and basal shear stress variations.

Our finding that seasonal signals in ice flow velocity may be linked to SSH implies that an improved understanding of ocean-driven ice shelf velocity variations at intra-annual timescales can provide valuable insights into the most efficient and accurate methods for modelling the likely future dynamic response of ice shelves and grounded ice sheets as climate and sea level changes. Similar to modelling the integrated effect of tidal loading in longer simulations, integrating the SSH effect would allow us to estimate the change in seasonal ice response associated with changes in seasonality of SSH. This may be important in the future if, for example, summer acceleration coincides and interacts nonlinearly with other seasonal forcings such as the near-ice-front basal melting (Klein et al., 2020). The small seasonal SSH changes that we observe and model here are actually similar in amplitude to the annual rates of sea-level rise that this ice shelf will experience in the future. Our results are directly relevant to other studies showing that a sea-level rate of ∼10 cm a⁻¹ could affect the grounding line migration by about 40 % with respect to models that do not include such processes (Larour et al., 2019).

Progress is needed in four areas: (1) seasonally resolved measurements of open ocean SSH; (2) ocean modelling, including all components (mass, steric height change, and inverse barometer) that contribute to SSH changes under ice shelves; (3) improved multi-year records of seasonally resolved ice velocity changes through either long-term continuous GNSS records or satellite-based methods; and (4) improved representation of grounding zone processes including subglacial hydrology, basal friction, and grounding line migration. Current satellite altimetry missions such as NASA's ICESat-2 can provide the SSH data close to ice fronts for validating and improving ocean models of SSH including under ice shelves, while concurrent GNSS measurements and reliable, data-constrained model estimates of sub-ice-shelf SSH can be used to identify optimal configurations for viscous models and for tuning grounding line parameterisations used in longer time integrations of ice shelf response to SSH changes.

Appendix A: Inverse and direct ice flow model

A1 Ice flow model and initialisation

Following Klein et al. (2020), all our simulations were conducted at the scale of the RIS basin, which encompasses the ice shelf and the grounded ice catchments that drain into RIS (Rignot et al., 2011; Fig. S7). We used a triangular finite element mesh with a spatial resolution that varies from 0.5 km at the grounding line to 20 km in regions of slow flow. The model spatial resolution on the ice shelf is typically ∼2 km.

The SSA model uses a vertically averaged effective ice viscosity with a nonlinear dependence on strain rate, and assuming isotropic material properties

$$\begin{array}{}\text{(A1)}& \mathit{\eta }={\mathit{\eta }}_{\mathrm{0}}{\mathit{\epsilon }}_e^{(\mathrm{1}-n)/n},\end{array}$$

where ε_e is the second invariant of the strain-rate tensor, η₀ is a vertically integrated apparent viscosity parameter, and n=3 is the value most consistent with field data and most commonly used in ice sheet modelling (Cuffey and Paterson, 2010). Bedrock elevation and ice thickness were taken from Bedmap2 (Fretwell et al., 2013), with a surface elevation correction applied to the floating ice to ensure floatation for an ice density of ρ_i=917 kg m⁻³ and a water density of ρ_w=1028 kg m⁻³. A Neumann condition, resulting from the hydrostatic water pressure exerted by the ocean on the ice, was applied at the calving front (Gagliardini et al., 2013), and a Dirichlet condition forced the normal velocities to zero on the inland boundary of the basins adjacent to RIS.

Our model inversion optimises both the basal friction coefficient (C) and the effective viscosity of the ice (η₀) by minimising multiple cost functions:

$$\begin{array}{}\text{(A2)}& J_{\mathrm{total}}=J_u+{\mathit{\lambda }}_{\mathrm{d}h/\mathrm{d}t}{J}_{\mathrm{d}h/\mathrm{d}t}+{\mathit{\lambda }}_{\mathrm{C}}{J}_C+{\mathit{\lambda }}_{{\mathit{\eta }}_{\mathrm{0}}}{J}_{{\mathit{\eta }}_{\mathrm{0}}},\end{array}$$

where J_u measures the difference between observed and modelled velocities, and $J_{\mathrm{d}h/\mathrm{d}t}$ measures the misfit between modelled and observed thickness rates of change, computed as the difference between flux divergence and mass balance (e.g. Brondex et al., 2019; Mosbeux et al., 2016). J_C and $J_{{\mathit{\eta }}_{\mathrm{0}}}$ are two regularisation functions added as constraints on the smoothness of the solution by penalising the first spatial derivatives of C and η₀. Three of the four cost functions are weighted by a regularisation parameter λ to allow us to give more or less weight to a function.

We ran an ensemble of 100 inversions, varying the different regularisation parameters $({\mathit{\lambda }}_{\mathrm{d}h/\mathrm{d}t}$, λ_C, ${\mathit{\lambda }}_{{\mathit{\eta }}_{\mathrm{0}}})$ as follows:

$$\begin{array}{c}{\displaystyle }{\mathit{\lambda }}_{\mathit{\beta }}=\mathit{\{}{\mathrm{10}}^{\mathrm{4}},\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}},{\mathrm{10}}^{\mathrm{5}},\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}},{\mathrm{10}}^{\mathrm{6}}\mathit{\}},\\ {\displaystyle }{\mathit{\lambda }}_{\mathrm{C}}=\mathit{\{}{\mathrm{10}}^{\mathrm{4}},\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}},{\mathrm{10}}^{\mathrm{5}},\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}},{\mathrm{10}}^{\mathrm{6}}\mathit{\}},\\ {\displaystyle }{\mathit{\lambda }}_{\mathrm{d}h/\mathrm{d}t}=\mathit{\{}{\mathrm{10}}^{-\mathrm{4}},\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}},{\mathrm{10}}^{-\mathrm{3}},\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}\mathit{\}},\end{array}$$

which leads to $N_{\mathrm{simulations}}=N_{{\mathit{\lambda }}_{\mathrm{C}}}\times N_{{\mathit{\lambda }}_{{\mathit{\eta }}_{\mathrm{0}}}}\times N_{{\mathit{\lambda }}_{\mathrm{C}}}=\mathrm{100}$.

The best members of the ensemble exhibit an ice flow pattern very close to observations, with an RMS velocity misfit (RMS(u)) as low as ∼10.1 m a⁻¹ and an RMS misfit on the ice thickness rate of change (RMS($\mathrm{d}h/\mathrm{d}t$)) as low as ∼0.7 m a⁻¹ over the grounded ice and the ice shelf combined (Fig. A1). From this ensemble, we obtained a sub-ensemble of 15 members (Ω₁₅) with misfit values below 15 m a⁻¹ on velocities and 1 m a⁻¹ on ice thickness rate of change (Fig. A1). Although this threshold on velocity is slightly higher than the data uncertainty reported by Rignot et al. (2011, 2017), both thresholds are close to the RMS misfits in other studies based on similar techniques (e.g. Gudmundsson et al., 2019; Brondex et al., 2019; Reese et al., 2018; Fürst et al., 2015). This ensemble of initial states, Ω₁₅, is then used for each of our simulations of grounding line migration (i.e. ΔL_B2, ΔL_C, and ΔL_B2L) for each model of SSH variability.

Figure A1 Ensemble of inversions (100 members, grey and blue points) in $\text{RMS}\left(u\right)-\text{RMS}(\mathrm{d}h/\mathrm{d}t)$ space. The vertical and horizontal grey boxes represent the sub-spaces RMS(u)<15 m a⁻¹ and $\text{RMS}(\mathrm{d}h/\mathrm{d}t)<\mathrm{1}$ m a⁻¹. The intersection of the two boxes represents the optimal sub-space (Ω₁₅), which contains 15 members (blue points).

A2 On the use of a diagnostic ice flow model

Klein et al. (2020) reported that the initial state obtained after inversion is not perfectly stable because of remaining uncertainties in other ice sheet parameters (see also e.g. Seroussi et al., 2011), which leads to locally large and unphysical ice thickness rates of change when running transient simulations (e.g. Brondex et al., 2019; Gillet-Chaulet et al., 2012; Klein et al., 2020). This problem is usually overcome by running a relaxation experiment, where the model is allowed to evolve under a constant forcing until a more stable state is reached and before applying the desired perturbation (e.g. Brondex et al., 2019; Gillet-Chaulet et al., 2012). However, this procedure sometimes incurs a significant cost in terms of the differences between observations and the modelled ice thickness and velocities. Although our initial states are similar to those in Klein et al. (2020), our experiment differs by the nature of the perturbation we apply. The basal melting investigated by Klein et al. (2020) directly affects the ice thickness, leading to a modification of the ice flow. The SSH deviations used here do not directly modify the ice thickness but rather modify the driving stress and grounding line position, which leads to a modification of the ice flow, eventually leading to a dynamical change in ice thickness. These changes in ice thickness are fairly small and can be neglected compared with changes in driving stress and grounding line position. Therefore, our model does not actually vary in time; instead, we apply the monthly averaged ΔSSH as a perturbation to the shallow-shelf model and calculate the difference in the velocity field between the perturbed model and the reference model. Monthly velocity change can therefore be determined and compared with the GNSS velocity variations.

Appendix B: Parameterisation of the grounding line migration

B1 Theory and equations

The grounding line migration under tidal variation is usually treated as a purely elastic and hydrostatic problem (Tsai and Gudmundsson, 2015). In this context, at the grounding line, the ice is lifted due to floatation, and the upward buoyancy force in the water column is compensated by the downward gravitational force in the ice column:

$$\begin{array}{}\text{(B1)}& F_{\mathrm{i}}=F_{\mathrm{w}}\iff {\mathit{\rho }}_{\mathrm{i}}gH={\mathit{\rho }}_{\mathrm{w}}g(z_{\mathrm{sl}}-z_{\mathrm{b}}),\end{array}$$

where z_sl is the sea level, z_b is the ice bed elevation, ρ_i is the ice density, and ρ_w is the water density.

Adapting Tsai and Gudmundson (2015), upstream from the grounding line, we can approximate the bed elevation at the point of migration by

$$\begin{array}{}\text{(B2)}& z_{\mathrm{b},\mathrm{\Delta }L}=z_{\mathrm{b},\mathrm{GL}}+\mathit{\beta }\mathrm{\Delta }L,\end{array}$$

with β being the bed slope (equal to the ice base slope if located upstream the grounding line) and ΔL the grounding line migration we try to estimate. Similarly, the ice thickness upstream the grounding line can be estimated as

$$\begin{array}{}\text{(B3)}& H_{\mathrm{\Delta }L}=H_{\mathrm{GL}}+(\mathit{\alpha }-\mathit{\beta })\mathrm{\Delta }L,\end{array}$$

with α being the surface slope. From there, we can rewrite

$$\begin{array}{}\text{(B4)}& {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}(H_{\mathrm{GL}}+(\mathit{\alpha }-\mathit{\beta }\left)\mathrm{\Delta }L\right)=\mathrm{\Delta }\text{SSH}-Z_{\mathrm{b},\mathrm{GL}}-\mathit{\beta }\mathrm{\Delta }L\end{array}$$

and estimate

$$\begin{array}{}\text{(B5)}& \mathrm{\Delta }{L}^{+}={\displaystyle \frac{\mathrm{\Delta }\text{SSH}}{{\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}}(\mathit{\alpha }-\mathit{\beta })+\mathit{\beta }}}.\end{array}$$

For the downstream migration, our assumption leads to a reduction of the ice base slope by a factor $\mathrm{1}/(\mathrm{1}-{\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}})\sim \mathrm{9}$ and therefore a potential grounding line migration:

$$\begin{array}{}\text{(B6)}& \mathrm{\Delta }{L}^{-}=\mathrm{\Delta }{L}^{+}\times (\mathrm{1}-{\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}}).\end{array}$$

Combining Eqs. (B5) and (B6), we obtain the following parameterisation:

$$\begin{array}{}\text{(B7)}& \mathrm{\Delta }{L}^{\pm }={\displaystyle \frac{\mathrm{\Delta }{S}^{\pm }}{{\mathit{\gamma }}^{\pm }}},\end{array}$$

where ΔS^± is the SSH perturbation in the grounding zone and

$$\begin{array}{}\text{(B8)}& {\mathit{\gamma }}^{+}=\mathit{\beta }+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}(\mathit{\alpha }-\mathit{\beta });{\mathit{\gamma }}^{-}={\displaystyle \frac{{\mathit{\gamma }}^{+}}{\mathrm{1}-{\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}}}}.\end{array}$$

This parameterisation assumes the surface and bed slope to be constant in the grounding zone, while average surface and bed slopes are potentially different immediately upstream versus immediately downstream of the grounding line. However, these differences are unlikely to ever be ∼10 times different. This is especially true for small migrations such as the ones of our hydrostatic model ΔL_B2, i.e. about a few tens of metres except in some areas of the Siple Coast and some Transantarctic Mountain glaciers (Fig. S8a); these migration distances are less than our model resolution and below the scales at which we expect large surface and bed slope variations. We also note that the effect of the larger ΔL_B2 over the Siple Coast (Fig. S8a) is mitigated by a relatively low basal shear stress, limiting the effect of the migration on the ice flow (Fig. S8b).

Figure B1 Schematic representation of the grounding line migration with sea surface height change (ΔS). (a) Flowline view with Δx being the element edge size at the grounding line and ΔL⁺ the upstream migration of the grounding line. (b) Two-dimensional plan view of the virtual migration (dotted blue line) of the grounding line (blue line) to an upstream (ΔL⁺) and downstream location (ΔL⁻); the evolution of the friction coefficient (C) is proportional to ΔL^±. Black and grey elements are initially grounded and floating.

B2 Parameterisations applied to the ice sheet model

The three parameterisations used in our study and presented in Sect. 2.3.2 are further detailed here:

1. ΔL_B2. We calculated ΔL_B2 by applying γ_B2 values corresponding to Bedmap2 bed slopes (e.g. ${\mathit{\beta }}_{\mathrm{B}\mathrm{2}}\sim [\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}-\mathrm{5}\times {\mathrm{10}}^{-\mathrm{2}}])$ and surface slopes (e.g. ${\mathit{\alpha }}_{\mathrm{B}\mathrm{2}}\sim {\mathit{\beta }}_{\mathrm{B}\mathrm{2}}/\mathrm{10}$ on the ice shelf and at the grounding line and ${\mathit{\alpha }}_{\mathrm{B}\mathrm{2}}\sim {\mathit{\beta }}_{\mathrm{B}\mathrm{2}}/\mathrm{40}$ when averaged over the entire basin) in Eq. (B8), where γ controls the length of the grounding line migration for a given ΔSSH. In the hydrostatic case, γ_B2 is calculated as a function of α and β.

2. ΔL_C. Following Rosier and Gudmundsson (2020), we calculated ΔL_C by applying constants for positive ${\mathit{\gamma }}_{\mathrm{C}}^{+}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ and negative ${\mathit{\gamma }}_{\mathrm{C}}^{-}={\mathit{\gamma }}_{\mathrm{C}}^{+}/\mathrm{9}$ bed slopes in Eq. (B8).

3. ΔL_B2L. We calculated ΔL_B2L by applying a coefficient ${\mathit{\gamma }}_{\mathrm{B}\mathrm{2}\mathrm{L}}={\mathit{\gamma }}_{\mathrm{B}\mathrm{2}}/\mathrm{20}$, with γ_B2L capped to ${\mathit{\gamma }}_{\mathrm{B}\mathrm{2}\mathrm{L}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ to limit extremely large grounding line migration in regions with very small γ_B2L values. This scaling factor was chosen so that the mean migration distance around the RIS perimeter was similar to ΔL_C.

B3 Subgrid-scale parameterisation

For $\mathrm{\Delta }{S}^{\pm }=\mathrm{10}$ cm (roughly the maximal modelled ΔSSH for RIS), $\mathit{\alpha }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$, and $\mathit{\beta }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$, Eqs. (B7) and (B8) lead to a $\mathrm{\Delta }{L}^{+}\sim \mathrm{100}$ m upstream and $\mathrm{\Delta }{L}^{-}\sim \mathrm{10}$ m downstream migration of the grounding line. These values are much smaller than the Δx∼500 m spacing of our model grid nodes in the vicinity of the grounding line.

We overcome this problem by parameterising the grounding line migration as a variation of the friction coefficient at the grounding line (Fig. B1). We define the initial basal shear force (“i” subscript – before migration) over the element edges surrounding grounding line nodes as

$$\begin{array}{}\text{(B9)}& F_{\mathrm{i}}={\mathit{\tau }}_{\mathrm{i}}\mathrm{\Delta }x,\end{array}$$

with τ_i=C_i|u_i|^m, where C_i is the reference friction coefficient, and u_i is the velocity on an element edge; we can write the shear force over a fraction Δx−ΔL of the last grounded element edge as

$$\begin{array}{}\text{(B10)}& F_{\mathrm{f}}={\mathit{\tau }}_{\mathrm{i}},\left(\mathrm{\Delta }x-\mathrm{\Delta }L\right).\end{array}$$

Equation (B4) can also be written as a function of a final shear stress (“f” subscript – after migration) integrated over the entire element:

$$\begin{array}{}\text{(B11)}& F_{\mathrm{f}}={\mathit{\tau }}_{\mathrm{f}}\mathrm{\Delta }x,\end{array}$$

with τ_f=C_f|u_f|^m where C_f is the friction coefficient at the grounded line node after migration of the grounding line.

Assuming $\left|\frac{u_{\mathrm{f}}}{u_{\mathrm{i}}}\right|\sim \mathrm{1}$, we can rewrite

$$\begin{array}{}\text{(B12)}& C_{\mathrm{f}}={\displaystyle \frac{\left(\mathrm{\Delta }x-\mathrm{\Delta }L\right)}{\mathrm{\Delta }x}}{C}_{\mathrm{i}},\end{array}$$

with C_f<C_i for ΔL>0 and C_f>C_i for ΔL<0.