Model setup

We define load as a change in the spatial distribution of mass (due to snow buildup or basal melting/accretion) or pressure at the ice-sheet surface between two points in time (before and after the change). The response to this load is evaluated in terms of the vertical displacement of two particular interfaces: the ice–water interface and the surface of the ice sheet. These interfaces can change their vertical position as a result of two mechanisms: first, a deformation, that is, the vertical displacement of firn/ice particles and, second, the accretion or erosion of firn/ice particles at the interface. In order to discriminate between both contributions in our model we introduce the particle horizon. This concept denotes a laterally continuous set of firn/ice particles within the ice sheet. This horizon participates in the deformation mechanism, through the relative vertical displacement of its constituting particles, but is not affected by the accretion/erosion contribution at the ice-sheet interfaces above and beneath. Illustrative examples for particle horizons are the surfaces of coeval ice particles, manifested as internal layers in radio-echo sounding profiles, whose constitutive particles represented the ice-sheet surface at the time of their accumulation (though a common accumulation time is not a necessary condition for a particle horizon in our case). In the following, we do not consider any (a) water volume change in Lake Vostok (i.e. no influx nor discharge, Richter and others, Reference Richter2014, and imposing a zero net basal mass balance over the lake area); (b) hydrodynamic lake-level variation; (c) change in the vertical firn/ice density profile and (d) vertical solid-earth deformation beneath the subglacial lake (e.g. due to glacial isostatic adjustment or tectonics).

Archimedes’ principle determines the vertical position of the ice-water interface at the base of ice floating in hydrostatic equilibrium on an extended water surface. It relates any load applied to the ice surface with the vertical displacement of the ice-water interface resulting from the hydrostatic adjustment to this load. In the case of Lake Vostok, the extent of the ice body exceeds that of the water surface. The hydrostatic equilibration is then confined to within the area of the subglacial lake. The hydrostatic adjustment to the load causes a bending of the ice sheet, which is assumed unfractured. This implies a finite relative vertical displacement of neighbouring ice columns, limited by the visco-elastic properties of the ice sheet. The lateral smoothness of particle horizon (e.g. ice-water interface) deformation leads to an attenuation of the hydrostatic response towards the lake shore, as well as across abrupt load gradients within the lake area. The hydrostatic ice-water interface adjustment complies, in addition, with the water volume conservation in the lake. This volume conservation leads to ice-water interface displacements also in those parts of the lake area unaffected by the load, but with inverted sign. Hence, load-induced height changes of the ice-water interface are determined by the following three constraints: (a) the hydrostatic equilibrium; (b) the flexural attenuation and (c) the water volume conservation. The surface height change is then obtained by adding the load thickness to the ice-water interface displacement. Figure 2 illustrates the effect of the three constraining processes on the considered interfaces in response to a surface load as, e.g. a snow-buildup anomaly.

Fig. 2. Schematic cross section of Lake Vostok. (a) Geometry of ice surface (black) and ice-water interface (red) in an initial situation before the application of a surface load. (b) Vertical deformation of the surface and ice-water interface (IWI) due to the hydrostatic equilibration of the load. The load is described by its height $h_{\rm l}$ and density $\rho _{\rm l}$; $h_{\rm w}$: vertical ice-water interface displacement; $\rho _{\rm w}$: water density. (c) Effect of the flexural attenuation, damping the vertical surface and ice-water interface deformation across the attenuation fringes along the lake shores and the load edge. (d) Effect of the water volume conservation in the lake on the vertical surface and ice-water interface displacements. The water volume $V$ lost due to the load-induced ice-water interface deformation (hashed area to the left) equals the volume gain (hashed area to the right). The dashed black line shows the equilibrated surface in the case of a pressure load equivalent to the mass distribution of the considered surface load.

A load acting on the surface of the floating ice sheet is described in any place $x$ by the areal load mass density $\kappa ( x)$. The areal density of a snow or ice load is determined by the thickness $h_{\rm l}$ and density $\rho _{\rm l}$ of the load:

(1)$$\kappa( x) = \rho_{\rm l} h_{\rm l}( x) .$$

In the case of atmospheric pressure loading the areal density is obtained by dividing the local pressure $P_{\rm l}( x)$ by gravity $g$:

(2)$$\kappa( x) = {P_{\rm l}( x) \over g};\; $$

i.e. the applied pressure can be represented equivalently by adding or subtracting mass atop the ice sheet.

The hydrostatic equilibration of a surface load (prior to accounting for flexural attenuation and water volume conservation) produces a vertical ice-water interface displacement $h_{\rm w}^{( 1) }( x)$ governed by the local areal density $\kappa ( x)$, expressed as the equivalent water column height by applying the density of water $\rho _{\rm w}$:

(3)$$h_{\rm w}^{( 1) }( x) = - {\kappa( x) \over \rho_{\rm w}}.$$

A positive $h_{\rm w}^{( 1) }( x)$ denotes an upward displacement of the ice-water interface. The hydrostatic response to a spatially non-uniform load produces a vertical deformation of the ice-water interface which follows the opposite of the spatial pattern of the areal density variations over the lake area. The vertical ice-water interface displacement is zero where the ice is grounded, that is, outside the subglacial lake and over islands. We keep the grounding line position fixed, neglecting the small migrations caused by load-induced ice-water interface displacements. The ice-water interface deformation differs from $h_{\rm w}^{( 1) }( x)$ due to the flexural attenuation at the transition from floating to grounded ice, as well as across steep load gradients. The attenuation is described by a factor $a( x)$ that depends on the distance to the nearest shore segment and varies from 0 (grounded, no ice-water interface displacement) to 1 (freely floating, complete hydrostatic equilibrium). Hence, the ice-water interface displacement with account for flexural attenuation becomes

(4)$$h_{\rm w}^{( 2) }( x) = a( x) h_{\rm w}^{( 1) }( x) = - {a( x) \kappa( x) \over \rho_{\rm w}}.$$

For the characteristic attenuation profile we adopt a linear shape over a length $L_{\rm a}$:

(5)$$a( x) = \left\{\matrix{ 0,\; \hfill & \rm{over\, grounded\, ice} \hfill \cr \displaystyle{d( x) \over L_{\rm a}},\; \hfill & \rm{if}\ d( x) \leq L_{\rm a} \hfill \cr 1,\; \hfill & \rm{if}\ d( x) > L_{\rm a} \hfill }\right.$$

where $d( x)$ denotes the distance to the nearest shore segment. The profile length $L_{\rm a}$ determines the width of the attenuation fringe along the shore, across which the ice-water interface deformation is damped. We adopt a profile length of 10 km. This choice is supported, on the one hand, by the surface deformation observed over Lake Vostok (Wendt and others, Reference Wendt2005) and, on the other hand, by the deviation of the static surface geometry from hydrostatic equilibrium (Ewert and others, Reference Ewert2012b; Schwabe and others, Reference Schwabe, Ewert, Scheinert and Dietrich2014); and will be elaborated on in the Discussion section.

In general, the ice-water interface deformation involves a volume excess $\Delta V$ over the lake area $S$:

(6)$$\Delta V = \iint_S h_{\rm w}^{( 2) }( x) \,\, {\rm d}S.$$

Water volume conservation is ensured through an additional vertical shift of the ice-water interface by $a( x) k$, where $k$ is constant over the lake area and $a( x)$ is, again, the flexural attenuation near the shores. Volume conservation implies that $k = - \Delta V / \tilde {S}$, where $\tilde {S} = \iint a( x) {\rm d}S$ is the effective area on which the volume excess $\Delta V$ is compensated. Hence, the final expression for the ice-water interface displacement is

(7)$$h_{\rm w}( x) = h_{\rm w}^{( 2) }( x) - a( x) {\Delta V\over \tilde{S}} = {a( x) \over \rho_{\rm w}} \left(- \kappa( x) + {\iint a( x) \kappa( x) {\rm d}S\over \iint a( x) {\rm d}S} \right).$$

Due to the subtraction of $k$, only the relative differences of the load within lake area are effective rather than the absolute value of the introduced mass or pressure load. In the case of atmospheric pressure loading, the induced surface height change $h_{\rm s}( x)$ is identical to the vertical ice-water interface displacement $h_{\rm w}( x)$. For snow loads, the surface height change is obtained by adding the load height $h_{\rm l}$ to the ice-water interface displacement:

(8)$$h_{\rm s}( x) = h_{\rm w}( x) + h_{\rm l}( x) .$$

In general, the spatial patterns and gradients of the ice-water interface and surface deformations depend on the patterns of load distribution and attenuation. Combining Eqns (1), (3) and (8) we see that, in the case of snow-buildup loads and neglecting attenuation effects, the gradients of surface vs ice-water interface displacements relate as $h_{\rm s}/h_{\rm w} = ( h_{\rm l} ( 1-\rho _{\rm l}/\rho _{\rm w}) ) / ( -h_{\rm l} \rho _{\rm l}/\rho _{\rm w}) = 1-\rho _{\rm w}/\rho _{\rm l} < 0$. The negative sign of this deformation gradient ratio denotes that the hydrostatic equilibration partitions an ice thickening between an upward growth at the surface and a downward growth at the ice-water interface. The surface snow density at Vostok ($\rho _{\rm l} =\, \sim \! 330$ kg m$^{-3}$) suggests furthermore that surface height change gradients are approximately twice as large as the simultaneous ice-water interface height change gradients.

The calculations are performed on a grid of the lake area (red box in Fig. 1) with 100 m spacing. The attenuation factors within the transition zone along the grounding line and the effective lake area $\tilde {S}$ are the same for any load scenario and are therefore computed only once. The ice-water interface displacement $h_{\rm w}( x)$ is calculated according to Eqn (7). The maximum gradient within the ice-water interface displacement grid $h_{\rm w}( x)$ is evaluated. If the gradient exceeds a threshold of 10$^{-5}$, the flexural attenuation across load edges is applied. For example, the geometrical description of load patterns can be simplified by polygons delimiting areas of constant areal density (Fig. 2). A polygon segment which separates different surface densities to either side is referred to as load edge. In this case, the load edge attenuation is implemented by smoothing the load areal density grid $\kappa ( x)$ within a 10 km wide attenuation band centred on the load edge by a linear interpolation between the $\kappa$ values on either side of the edge (consistent with the linear, 10 km long attenuation profile applied for the shore attenuation). The calculation of the $h_{\rm w}( x)$ grid according to Eqn (7) is then repeated introducing the smoothed $\tilde {\kappa }( x)$ grid. The surface height jump at the load edge shown in Figure 2c equals the height of the load in Figure 2a, $h_{\rm l}$.

Different geometrical load discretizations (patterns) are adopted in the computational experiments presented below: a uniform load within a polygon (pattern P in Table 1), a regional gradient (G) or an areal density field interpolated from the gridcell pattern of a regional climate model onto our computation grid (I).

Table 1. Summary of the diagnostic scenarios for which the hydrostatic response over Lake Vostok is modelled

^aLoad patterns denote: P, uniform load bounded by polygon; I, interpolation from RACMO cell centres onto our calculation grid; G, uniform gradient.

Mass changes at the ice base require a special treatment. In this case, the particle horizon to which the flexural attenuation is applied does not coincide with the ice-water interface. Thus, the attenuation and the volume conservation conditions are evaluated over different surfaces (attenuation: particle horizon, e.g. ice surface; volume conservation: ice-water interface). A preliminary particle horizon deformation $h_{\rm p}^{( 1) }$ is determined according to isostasy and shore attenuation:

(9)$$h_{\rm p}^{( 1) }( x) = a( x) \kappa( x) \left({1\over \rho_{\rm l}} - {1\over \rho_{\rm w}} \right).$$

The ice-water interface displacement implied by this particle horizon deformation and the basal load height $h_{\rm l}$ yields a preliminary volume excess

(10)$$\Delta V = \iint_S ( h_{\rm p}^{( 1) }( x) - h_{\rm l}( x) ) \, {\rm d}S.$$

This volume excess is corrected through an increment to the preliminary particle horizon deformation, analogous to the case of surface loads (Eqn (7)), yielding

(11)$$h_{\rm p}( x) = h_{\rm p}^{( 1) }( x) - a( x) {\Delta V\over \tilde{S}}.$$

The surface height change corresponds to the particle horizon deformation $h_{\rm s}( x) = h_{\rm p}( x)$ and the final ice-water interface displacement is obtained by subtracting the load height

(12)$$h_{\rm w}( x) = h_{\rm p}( x) - h_{\rm l}( x) .$$

In our calculations we apply a standard water density of $\rho _{\rm w} = 1000$ kg m$^{-3}$. For snow loads we adopt the mean density of the uppermost 20 cm snow layer derived from stake and snow pit observations at Vostok station of $\rho _{\rm l} = 330$ kg m$^{-3}$ (Ekaykin and others, Reference Ekaykin, Lipenkov, Petit and Masson-Delmotte2003). For basal ice loads we apply an ice density of $\rho _{\rm l} = 923$ kg m$^{-3}$ as derived from the Vostok ice core (Lipenkov and others, Reference Lipenkov, Salamatin and Duval1997).

Loading scenarios

We consider snow build-up, basal ice mass changes and air pressure changes as loads. Spatio-temporal snow-buildup variations are represented by the surface mass balance estimates of the regional climate model RACMO2.3p2 (henceforth RACMO; van Wessem and others, Reference van Wessem2018). Stake observations in the Vostok stake farm and along three profiles across the Lake Vostok area were used to validate this surface mass balance model (Richter and others, Reference Richter2021). We use RACMO's surface mass balance estimates on a grid with a spatial resolution of 27 km, over the period 1979–2019 with monthly resolution, to extract snow-buildup load fields over Lake Vostok.

The basal mass balance of the ice sheet floating on top of Lake Vostok is characterized by simultaneous melting and freezing at the ice-water interface in different parts of the lake. Models (Wüest and Carmack, Reference Wüest and Carmack2000; Mayer others, Reference Mayer, Grosfeld and Siegert2003; Thoma and others, Reference Thoma, Grosfeld and Mayer2007; Salamatin and others, Reference Salamatin, Tsyganova, Popov and Lipenkov2009) and the sparse available observations (Bell and others, Reference Bell2002; Tikku and others, Reference Tikku, Bell, Studinger and Clarke2004) agree that basal accretion prevails in the southern half of the lake area, while basal melting dominates in the northern part. Local basal mass balance estimates were derived for ten polygons distributed over the lake area from a combination of 3-D particle displacement velocities observed by GNSS with ice thickness determined by radio-echo sounding (Richter and others, Reference Richter2014). These estimates suggest that the general basal mass balance variation can be approximated by an N–S gradient of 0.2 kg m$^{-2}$ a$^{-1}$ km$^{-1}$.

Atmospheric pressure loading is effective over Lake Vostok according to the inverse barometer response to temporal changes in the air pressure field. The air pressure is continuously recorded at Vostok, but the lack of nearby additional meteorological stations impairs the retrieval of air pressure gradients representative for the lake area. According to Wendt and others (Reference Wendt2005) and Richter and others (Reference Richter2014) the N–S air pressure gradient over the lake area varies within $\pm 2$ Pa km$^{-1}$. Here, we use the air pressure fields predicted by RACMO over the period 1979–2019 with a temporal resolution of 3 h.

Our modelling scheme is applied to six diagnostic load scenarios. Table 1 summarizes the scenarios, which are briefly described in the following:

Scenario 1 simulates a snow-buildup anomaly of 24.8 mm restricted to the southern part of the subglacial lake. It is intended to visualize the principal concepts of the hydrostatic load response over the subglacial lake through a hypothetical, geometrically simple case, corresponding to Figure 2. The adopted load height represents the maximum snow buildup observed over 1 month (September 1972) within the stake farm record at Vostok station 1970–2019 (Ekaykin and others, Reference Ekaykin2020; Richter and others, Reference Richter2021).

Scenario 2 mimics the annual snow buildup over Lake Vostok according to the mean surface mass balance predicted by RACMO, averaged over the model period 1979–2019.

Scenario 3 examines the effect of the annual basal ice mass change. It applies the adapted calculation procedure in Eqns (9)–(12) and a gradient load pattern, which assumes that the basal mass balance variation over the Lake Vostok area is dominated by an N–S gradient of 0.2 kg m$^{-2}$ a$^{-1}$ km$^{-1}$, implying an ice load height ranging from $-24.4$ mm (melting, northern shore) to $+ 33.7$ mm (accretion, southern shore), while conserving the water volume in the lake.

Scenario 4 combines the annual surface (scenario 2) and basal (scenario 3) mass changes. It represents the mean annual load field accounting for both snow buildup and basal accretion/melting.

Scenario 5 represents an extreme monthly snow-buildup load. Among the 488 monthly surface mass balance predictions by RACMO over the period 1979–2019 we chose that month (August 1992) whose surface mass balance field produces the maximum response in terms of the standard deviation (SD) of the ice-water interface displacements $h_{\rm w}( x)$ over the lake.

Scenario 6 exemplifies an extreme event of atmospheric pressure load. The RACMO air pressure time series 1979–2019 were screened to select that air pressure field (22 February 2011 6:00 UTC) which produces the maximum SD of the ice-water interface displacements $h_{\rm w}( x)$ relative to the mean air pressure field over the model period.

Our model does not account for the dynamical ice mass effect of the ice flow. This is especially important for the interpretation of the effects of annual surface mass balance and basal mass balance applied in scenarios 2–4. The observational evidence presented by Richter and others (Reference Richter2014) suggests that the ice sheet in the Lake Vostok region has been very close to steady state over an observation period of 10 years. Under such conditions, the sum of surface mass balance and basal mass balance over the subglacial lake is compensated by the net ice flux into (along the western shore) and out of (eastern shore) the lake area.

Simulation of geodetic observables

We model three types of geodetic observables: (1) surface elevation changes at ICESat ground track crossover locations; (2) particle subsidence rates inferred from repeated GNSS observations and (3) vertical particle displacements derived from continuous GNSS observations at Vostok station. Observables (2) and (3) are physically similar, but differ in their spatial and temporal coverage which affects also the observational accuracy and their sensitivity to different processes. Observable (2) reflects spatially heterogeneous differences between discrete observation epochs, whereas observation set (3) yields local variations over diverse timescales.

1. (1) We apply our modelling scheme to explore the effect of spatio-temporal variations in snow buildup and air pressure at 60 crossover locations of the ICESat ground track over the lake during the 18 laser operation periods of the altimetry mission 2003–09. The adjustment of surface elevation differences in these crossover locations were used to determine laser campaign biases. Those laser campaign bias determinations assumed a constant surface elevation (Ewert and others, Reference Ewert2012b; Richter and others, Reference Richter2014) or a constant elevation change rate (Schröder and others, Reference Schröder2017) in these calibration locations throughout the mission period. The distribution of these calibration locations is shown in Figure 3. Here, we introduce the surface mass balance and air pressure variations estimated by RACMO into our modelling scheme in order to assess the impact of their neglect in the laser campaign bias estimation. For this purpose, we derive integral estimates for each laser operation period and the entire set of calibration locations, instead of computing individual corrections for all surface elevation observations introduced in the laser campaign bias adjustment. In the case of the snow load this is justified by its relatively long variation timescales. The local response to air pressure variations between the individual altimetric data acquisitions, in turn, is significant but not accounted for in our estimate. Nevertheless, the laser campaign bias are averaged in space ($\approx \!60$ calibration locations) and time (laser operation periods of $\approx \!1$ month) over multiple satellite passes across the lake area. As will be shown below, these spatio-temporal averages are dominated by the snow load, justifying our determination of integral estimates. Within each laser operation periods, $h_{\rm s}$ are calculated at the 60 calibration locations and with 3 hourly temporal resolution based on the air pressure fields simulated by RACMO. These $h_{\rm s}$ are averaged yielding, for each laser operation periods, the mean surface height change and standard deviation due to the inverse barometer effect. In order to assess the effect of snow-buildup variations, cumulative surface mass balance anomalies (csurface mass balanceA) are derived from the surface mass balance predicted by RACMO, assuming steady state for the RACMO model period 1979–2019. For each laser operation periods a mean snow-buildup field is inferred by weighting the monthly csurface mass balanceA fields proportionally to the number of days of each month contained in the laser operation periods. The $h_{\rm s}$ are computed for all laser operation periods and calibration locations, and are averaged to the mean surface height change, and its SD, due to snow buildup. The combined effect of air pressure and snow-buildup variations on each laser campaign bias results from summing the mean surface height changes derived for both processes, while the SDs are summed quadratically. These estimates can be considered as corrections to the laser campaign biases derived from crossover adjustments over Lake Vostok, improving the bias values by accounting for the loading effects. Henceforth, they are referred to as laser campaign bias corrections. In addition, the computation of the mean surface height changes among the laser operation periods is repeated based purely on the csurface mass balanceA fields, without consideration of the hydrostatic equilibration. This mimics the effect identical surface mass balance conditions would have over grounded ice, thus its comparison with the results that incorporate the hydrostatic adjustment, including the air pressure variations, reveals the impact of the presence of the subglacial lake on the laser campaign bias estimation. Moreover, the effect of the hydrostatic adjustment at individual calibration locations is exemplified by examining the $h_{\rm s}$ variations in four crossover locations. We select two of the 60 calibration locations close to the northern and southern shores, but sufficiently far away from the shoreline for an undamped hydrostatic response ($a( x)$ = 1). Both selected locations are complemented by the closest ICESat crossover locations outside the lake, over grounded ice. The two pairs of locations are shown in Figure 3.

1. (2) We assess the impact of surface load changes on particle subsidence rates inferred from repeated GNSS observations. Each of the 56 markers reported by Richter and others (Reference Richter2014, Table 1) was occupied during multiple (2–7) Antarctic field seasons between 2001 and 2013 over observation periods between 1 and 60 d. Here, we focus on those 42 markers which are located within the lake area ($a( x) > 0$; included in Fig. 3). Daily coordinate solutions for each marker and observation day were obtained from the GNSS data processing. Based on these daily solutions vertical velocities were derived which range from $-33$ to $-106$ mm a$^{-1}$ with uncertainties between $\pm 2$ and $\pm 22$ mm a$^{-1}$. The hydrostatic response to surface loads contributes to the vertical particle movement at the base of these GNSS markers and may contaminate the velocity estimates used to describe the long-term particle movement. For each day a coordinate solution is present for any of the considered markers, the hydrostatic response of the ice-water interface to snow-buildup anomalies and air pressure variations is modelled and interpreted in terms of the load-induced particle movement. The monthly csurface mass balanceA grids are interpolated in time to the GNSS observation day, whereas the pressure field for that day is obtained by averaging the 3 hourly model simulations over 24 h. In this way, coordinate corrections are obtained for each marker and observation day. We evaluate the amount by which the application of these corrections changes the reported vertical velocities.

2. (3) The vertical particle movement derived from continuous GNSS observations at Vostok station over 10 years is examined in search of observational evidence for the effects of snow-buildup anomalies and air pressure variations. For this purpose, we use the time series of daily solutions in the vertical coordinate component between early-2008 and late-2018. This represents an update of the time series presented by Richter and others (Reference Richter2014, Fig. 5) and results from a post-processing of dual-frequency GPS observations in differential mode using Bernese GNSS Software (Dach and others, Reference Dach, Lutz, Walser and Fridez2017). Non-linear signals in the vertical particle movement are extracted by detrending the time series. Monthly mean relative vertical particle positions are derived by averaging the detrended daily vertical coordinates within each month. In addition, the ice-water interface displacements our modelling scheme yields at the location of the GNSS station as response to the snow load represented by the monthly csurface mass balanceA fields derived from RACMO are detrended and compared with the observed non-linear position variations. For the detection of the inverse barometer effect on the particle movements observed at Vostok, the time series of daily vertical positions are high-pass filtered (filter width 30 d). The air pressure fields estimated by RACMO with 3 hourly resolution are averaged over 24 h to daily mean pressure load fields. These are introduced into our modelling scheme to yield a time series of the ice-water interface displacement at the GNSS station with daily resolution. This time series is filtered identically to the GNSS position time series.

Fig. 3. Map of the Lake Vostok area (model domain). Grey dots: 60 ICESat track crossover locations for which the impact of the hydrostatic load response is evaluated (Figs 10a, b); blue/red dots and diamonds: two pairs of ICESat crossover locations in the northern (blue) and southern (red) parts of the lake used to compare the surface height change between hydrostatically balanced (dots) and grounded (diamonds) ice (Fig. 10c); green triangles: location of 42 GNSS surface markers (Richter and others, Reference Richter2014) for which the impact of the hydrostatic load response on vertical particle velocities is evaluated.