1. Introduction
Ice cores are of crucial importance in the study of Earth’s climate and are fundamental to understanding recent climate changes attributed to anthropogenic forcing, in the context of a long history of natural climatic variability. Ice domes are particularly well-suited for drilling ice cores to study past climates because both horizontal velocity and shear are assumed to be minimal at a stationary ice dome. As a result, the stratigraphy of cores extracted from ice dome sites is not expected to be significantly distorted by ice dynamics, allowing for a monotonically increasing ice age with depth.
The Earth’s oldest ice core was retrieved as part of the European Project for Ice Core in Antarctica (EPICA) at the French/Italian Concordia Station (EPICA community members, 2004) in Dome C (DC), East Antarctica. The DC is the fourth highest dome in Antarctica, with an elevation of 3233 m, and is located ∼1200 km from the Southern Ocean and the Ross Sea (see Figures 1 and 2). This site represents the eastern culmination of the East Antarctic Ice Sheet (EAIS). The EPICA Dome C ice core (EDC) extends to a depth of 3193 m, yielding dated records corresponding to ages of 808 000 years (Bazin and others, Reference Bazin2013), making it the oldest continuously retrieved ice core to date. The European Beyond EPICA consortium project aims to obtain a continuous ice-core record of the past 1.5 million years from Little Dome C (LDC), a secondary dome located about 35 km southeast of Concordia Station (Lilien and others, Reference Lilien2021). In the domain of ice core research, age-depth models play a central role in constraining the timing, sequence and duration of past climate events and stages (e.g., Kawamura and others, Reference Kawamura2007; Parrenin and others, Reference Parrenin, Remy, Ritz, Siegert and Jouzel2004; Reference Parrenin and Hindmarsh2007; Bazin and others, Reference Bazin2013; Chung and others, Reference Chung2024). The ideal approach is to establish a chronology of ice core timelines through the counting of annual layers based on chemical and physical characteristics. However, the cumulative effects of low accumulation, wind redistribution, vertical compression and ice flow distortion prevent the application of these techniques on a universal scale. The counting of layers is frequently unfeasible for deep Antarctic ice cores from areas of low accumulation, and absolute time markers are generally not available.
Map showing the Dome C area. The map illustrates subglacial bedrock elevation above mean sea level (in metres), surface elevation contours (in increments of 100 m), drainage systems, subglacial water flux, and the grounding line position during the Last Glacial Maximum. This map was created using the Quantarctica GIS package (Matsuoka and others, Reference Matsuoka2021), developed by the Norwegian Polar Institute and published under the Creative Commons Attribution 4.0 International License.
Map of Antarctica and the Southern Ocean. The map illustrates the horizontal velocity of bedrock movements attributed to plate tectonics, represented by red arrows (in mm a−1). The black arrow indicates the absolute movement measured at Concordia Station (DCRU) due to ice dynamics and plate tectonics at the bedrock. Additionally, the blue arrow represents the estimated movement of the ice summit at DCRU in relation to the bedrock.
The chronology of ice core records is contingent on the precision of ice flow models that correlate core depth with age. A comparison between the EDC record, other ice cores (Parrenin and others, Reference Parrenin2012; Bazin and others, Reference Bazin2013; Fujita and others, Reference Fujita, Parrenin, Severi, Motoyama and Wolff2015) and the LR04 marine stack records (Lisiecki and Raymo, Reference Lisiecki and Raymo2005) reveals discrepancies with regard to the timing and duration of events along the core (e.g., Bazin and others, Reference Bazin2013; Dreyfus and others, Reference Dreyfus2007; Durand and others, Reference Durand2007; Tison and others, Reference Tison2015).
Ice sheet flow is a crucial parameter in numerical modelling, as it can validate a model’s ability to reproduce present-day geometry. The gravitational driving stress that generates this flow is proportional to the surface slope and the input boundary conditions (Bamber and others, Reference Bamber, Vaughan and Joughin2000). Models of depth–age relationships for deep ice cores are sensitive to the migration of the dome position (Anandakrishnan and others, Reference Anandakrishnan, Alley and Waddington1994). The behaviour of an ice divide is influenced by its accumulation rate history, as well as the spatial patterns and conditions at the ice sheet boundaries (e.g., Hindmarsh and others, Reference Hindmarsh, King, Mulvaney, Corr, Hiess and Gillet-Chaulet2011; Nereson and others, Reference Nereson, Raymond, Waddington and Jacobel1998; Frezzotti and others, Reference Frezzotti2004; Urbini and others, Reference Urbini2008).
It is worth noting that surface elevations at DC, Vostok and Dome Fuji have exhibited fluctuations of up to 100–150 m between glacial and interglacial periods, driven by changes in accumulation (Ritz and others, Reference Ritz, Rommelaere and Dumas2001). Cavitte and others (Reference Cavitte2018) reconstructed paleo-accumulation rates for the DC region over the past 73 000 years using a one-dimensional pseudo-steady ice flow model (Parrenin and others, Reference Parrenin2017), complemented by isochronal constraints derived from Radio Echo Sounding (RES) surveys. The authors highlight that the large-scale surface accumulation gradient (spanning hundreds of kilometres) has remained spatially stable over the past 73 000 years, reflecting the current atmospheric reanalysis modelled by the European Centre for Medium-Range Weather Forecasts and the observed snow accumulation gradients in the region (Urbini and others, Reference Urbini2008; Le Meur and others, Reference Le Meur, Magand, Arnaud, Fily, Frezzotti, Cavitte, Mulvaney and Urbini2018). However, Urbini and others (Reference Urbini2008) noted a shift in the spatial accumulation pattern at a local scale (tens of kilometres) during the last four centuries, as observed by their Ground Penetration Radar (GPR) surveys.
Determining the summit point of the dome and its migration over time presents a significant challenge due to the extremely low slope (less than 0.1 m km−1) and the surface micro-morphology (up to 0.2 m) at the metre scale (e.g., dunes, sastrugi). The extraction of information regarding the dynamics of the dome summit is only possible through the repeated measurements of ice velocity (Vittuari and others, Reference Vittuari2004; Yang and others, Reference Yang2014) at specific locations using Global Navigation Satellite System (GNSS). This is necessary due to the extremely slow velocities at the summit, which approach zero and are comparable to continental movements (on the order of mm a−1). Additionally, the reference geodetic markers located on stable rock outcrops are at considerable distance (over 1000 km).
As part of the EPICA-EDC, Concordia Station, and ITASE (International Trans-Antarctic Scientific Expedition) projects, a strain network composed of 37 poles was established and measured three times (in 1996, 1999 and 2012/14). The first repetition of the strain network observations (from 1996 to 1999) was documented by Vittuari and others (Reference Vittuari2004), who reported that poles closest to the DC summit moved at rates of up to a few mm a−1, while those located 25 km from the summit moved at rates of 211 ± 7 mm a−1.
In addition, a pole equipped with a “Submergence Velocity - Coffee Can” system (Hulbe and Whillans, Reference Hulbe and Whillans1994; Hamilton and Whillans, Reference Hamilton and Whillans2002) was installed in 1999 and has since been measured five times (Dome C Coffee Can, DCCC) to monitor long-term thickening or thinning of the ice sheet. This method involves measuring the submergence velocity of subsurface markers and comparing that velocity with both present and long-term local snow accumulation. Since 2005, a permanent GNSS station (DCRU) has been installed on the roof of the Concordia Station noisy tower, better known as “Rumoroso” building.
To improve understanding of the present dynamics of DC and LDC, it is essential to establish how the dome responds to both current and past environmental conditions. The objective of this paper is twofold: firstly, to analyse GNSS data at both continental and local scales; secondly, to compare the results with the current glaciological conditions of the DC area and its catchments. Additionally, this study aims to provide insights into the behaviour of ice dynamics at DC based on ice velocity measurements collected over the past two decades. A numerical model analysis of ice dynamics falls outside the scope of this paper and will be addressed in future research.
2. Geographical and glaciological setting of dome C
Legrésy and others (Reference Legresy, Rignot and Tabacco2000) were the first to provide a detailed description of ice flow around DC using synthetic aperture radar (SAR) interferometry. Since the 1990s, various authors (Tabacco and others, Reference Tabacco, Passerini, Corbelli and Gorman1998; Hodgkins and others, Reference Hodgkins, Siegert and Dowdeswell2000; Rémy and Tabacco, Reference Rémy and Tabacco2000; Siegert and others, Reference Siegert, Eyers and Tabacco2001; Forieri and others, Reference Forieri, Zuccoli, Bini, Zirizzotti, Remy and Tabacco2004; Zirizzotti and others, Reference Zirizzotti, Cafarella and Urbini2012; Passalacqua and others, Reference Passalacqua, Ritz, Parrenin, Urbini and Frezzotti2017; Urbini and others, Reference Urbini, Cafarella, Tabacco, Baskaradas, Serafini and Zirizzotti2015; Cavitte and others, Reference Cavitte2016; Reference Cavitte2021; Young and others, Reference Young2017) have contributed valuable information regarding the three-dimensional structure of the ice sheet in the DC area using satellite and airborne RES data. The bedrock within 12.5 km of the DC summit is relatively flat, with an ellipsoidal height nearing zero relative to the World Geodetic System 1984 (WGS84). Additionally, the ice thickness at the EDC site has been measured at 3309 ± 22 m (The EPICA Dome C 2001-02 Science and Drilling Teams, 2002 report).
The most prominent characteristic of DC’s surface morphology is its positioning along an ice divide that extends from Vostok Lake (3488 m) to the Adélie Coast (see Figure 1). An ice saddle, situated at an elevation of 3200 m, 115 km south-southeast of the dome summit, separates DC from the East Antarctic ice divide originating from Vostok Lake. The elevation difference between the summit and the saddle is a mere 33 m. This saddle is influenced by the deep Aurora Trench, which contains one of the thickest ice covers in Antarctica, measuring 4,755±16 m, along with a significant presence of subglacial lakes and basal melting (Caffarella and others, Reference Cafarella, Urbini, Bianchi, Zirizzotti, Tabacco and Forieri2006).
DC serves as the culmination of different drainage basins: the northwest-west side drains into the Aurora subglacial basin and the Southern Ocean, with the Totem Glacier representing the main glacier in this area. On the contrary, the southeast-east side drains into the Wilkes Subglacial Basin, leading to the Ross Sea, with the Byrd Glacier constituting the primary ice outlet glacier.
Another significant characteristic of DC’s surface orography is the elliptical shape (see Figures 1 and 3), where the minor axis (NW–SE) is about 70% shorter than the major axis (SW–NE) (Rémy and Tabacco, Reference Rémy and Tabacco2000). This morphology can be explained by Nye’s (Reference Nye1991) theory on the topology of ice-sheet centres. The axes of the ellipse are asymmetrical; the northeastern part of the minor axis is about 20% wider than its southwestern counterpart, while the southeastern flank of the major axis is steeper than the northwestern flank by 10% at 25 km from the divide. The direction of the major axis deviates by 20° from the prevailing SSW–NNE (340°) wind direction, as reported by (Frezzotti and others, Reference Frezzotti2005).
Surface ice velocity and topographic map. The map displays surface elevation contours at 1.0 m intervals, derived from Howat and others (Reference Howat2022). Ice velocity contours and errors from the 1996 to 2012/14 period are shown at 25 mm a−1 intervals. Vectors indicate the velocity and flow direction for the same period (see Table 1).
Rémy and Tabacco (Reference Rémy and Tabacco2000) observed that the strain rate in the SE–NW direction is three times lower than in the perpendicular direction. This hypothesis is further substantiated by GNSS velocity data (Vittuari and others, Reference Vittuari2004) along the minor (SE–NW) and major (SW–NE) axes of the dome. The average velocities along each axis are directly proportional to the distance from the centre (R 2 > 0.9), thereby confirming that the velocity in the SW–NE direction is three times lower than in the NW–SE direction (Vittuari and others, Reference Vittuari2004).
The morphology of DC resembles a ‘regular’ ellipse within the first 10 km from the summit, beyond which this regularity is interrupted by the presence of two distinct ‘valleys’. The first of these extends in the direction of the GNSS poles E13, E18 and E19 (NNW–SSE), while the second corresponds to the Concordia Subglacial Trench, aligned with the GNSS poles B11 and C13 (S-N) (see Figures 1 and 3). The surfaces of these valleys are influenced by the underlying bedrock morphology, the presence of humid conditions at the ice/rock interface and the existence of subglacial lakes (Zirizzotti and others, Reference Zirizzotti, Cafarella and Urbini2012; Urbini and others, Reference Urbini, Cafarella, Tabacco, Baskaradas, Serafini and Zirizzotti2015; Young and others, Reference Young2017).
The bedrock morphologies in the DC area are likely the result of the region’s tectonic setting (Cianfarra and others, Reference Cianfarra, Forieri, Salvini, Tabacco and Zirizzotti2009). In the southwestern part, along the ice divide originating from Vostok, the summit area is bordered by a complex chain of mountains oriented NW–SE, reaching a maximum altitude of ∼600 m above sea level, with the ice thickness reduced to 2700 m (Rémy and Tabacco, Reference Rémy and Tabacco2000; Young and others, Reference Young2017; Cavitte and others, Reference Cavitte2021). The surface expression of this mountain range is the secondary dome, referred to as ‘LDC’, positioned ∼40 km southwest of the DC summit.
3. Materials and methods
The GNSS network consists of 38 aluminium poles (each with a diameter of 12 cm and a length of 3 m) arranged in four concentric rings at increasing distances from the summit. The configuration includes six poles in each of the first two rings at 3 km and 6 km from the summit, 12 poles at 12.5 km and 25 km. Additionally, a reference service pole was installed during the first campaign, along with a pole for the ‘coffee can’ submergence velocity measurement site at the dome summit (see Table 1 and Figure 3).
The GNSS pole network at Dome C
This table shows the locations of the poles, and their annual movements derived from measurements taken between 1996 and 2014. Surface Mass Balance (SMB) is indicated in water equivalent (we) per year. The measurements from 2014 are highlighted in bold. Detailed information about the errors associated with the estimates can be found in Supplementary Tables 1 and 2.
The DC network was established and surveyed for the first time in late December 1995 (hereafter referred to as 1996) and re-measured in early January 1999. The third observation period started in January 2012 and completed between late December 2013 and early January 2014 (hereafter referred to as 2014). It is important to note that all measurements were conducted during the same annual period, with the objective of significantly mitigating potential seasonal variation. The GNSS poles were installed at a minimum depth of 1 m in the snow and equipped with devices for forced centring of the GNSS antennas, which acted as three-dimensional reference points.
The DC network was designed using detailed surface topography data derived from satellite radar altimetry (Brisset and Rémy, Reference Brisset and Rémy1996), supplemented by a continuous kinematic GNSS surface survey conducted during the 1993/94 expedition (Cefalo and others, Reference Cefalo, Tabacco, Manzoni and Bologna1996). The summit position, which was recorded at approximately 75°06'06'' S, 123°23'43'' E, was confirmed by subsequent detailed continuous kinematic GNSS surface elevation surveys carried out in 1995/1996 (Capra and others, Reference Capra, Cefalo, Gandolfi, Manzoni, Tabacco and Vittuari2000). Due to logistical constraints, the Concordia Station base summer camp, which housed the EPICA-EDC drilling site, was establishe during the same period at a location ∼1400 m west of the orographic summit dome (Vittuari and others, Reference Vittuari2004).
In consideration of the temporal context in which these measurements were undertaken, the receivers employed were predominantly equipped to receive signals from the Global Positioning System (GPS) constellation. The 1996 survey was performed using five dual-frequency Trimble 4000 series receivers, with a receiver set to continuous acquisition throughout the measurement period at pole 1000, located near the French traverse vehicle caravan, where a power generator was available. At the time of the study, the construction of the Concordia station had not yet been built. The data acquired at pole 1000 facilitated the connection to GNSS permanent stations belonging to the International GNSS Service (IGS). The duration of measurements at the remaining poles ranged from a minimum of 1 hour for proximate poles (3 and 6 km rings) to a maximum of 2 and a half hours for more distant poles (12.5 and 25 km rings), thereby ensuring consistent precision across the network.
During the 1999 campaign, seven dual-frequency Trimble 4000 series receivers were utilized, thereby enabling the execution of longer measurement sessions, with acquisition periods ranging from 2 to 24 hours. This was made possible by the presence on site of housing modules and tracked vehicles from the Italian leg of the ITASE traverse. As in the case of the previous campaign, a GPS station installed at the Concordia summer camp was maintained in continuous acquisition throughout the measurement period.
In the following 2012/14 campaign, the measurement points were occupied through multiple 24-hour daily sessions, enabled by advancements in instrumentation that featured much larger internal memory and the capability to recharge batteries on-site using solar panels. Additionally, a permanent GNSS station was activated on the roof of Concordia Station (DCRU). The instruments employed in this campaign included four Trimble 5700 receivers, two Trimble R7 GNSS receivers, and one Topcon GB100 receiver, all equipped with Choke Ring antennas. Prior to the 2012/14 measurements, the original poles installed in 1996, which exhibited minimal protrusion from the ice, were extended by 1.5 m. These extensions were executed with great care, and the slant to the vertical was measured in two orthogonal directions (N–S and E–W). These extensions are expected to keep the network operational for the next 20 years (see Supplementary Table 2).
The DCRU station, which was installed at a height of more than 17 m above the snow surface, has continuously recorded data at a rate of 15 seconds since January 2005. Three significant gaps in data acquisition occurred: from April 25, 2007, to January 18, 2008; from June 12, 2009, to December 5, 2009; and from January 11, 2014, to December 7, 2014 (see Figure 5). Levelling measurements of the tower, conducted since 2007, have not revealed any significant tilting, indicating a negligible effect on the GNSS permanent station measurements (personal communication from C. Le Calvez, 2018).
To accurately determine the ice velocity in the DC drainage area, a series of poles were installed during ITASE traverses conducted between 1998 and 2002. These poles were subsequently measured using GNSS between 2003 and 2006 (see Supplementary Table 3; Vittuari and others, Reference Vittuari2004; Frezzotti and others, Reference Frezzotti, Capra and Vittuari1998; Reference Frezzotti, Tabacco and Zirizzotti2000; Reference Frezzotti2004; Reference Frezzotti, Urbini, Proposito, Scarchilli and Gandolfi2007; this paper).
The “coffee can” submergence velocity measurement system was installed in the DCCC borehole and equipped with GNSS instrumentation to evaluate the long-term thickening and thinning of the dome area (Hulbe and Whillans, Reference Hulbe and Whillans1994; Hamilton and Whillans, Reference Hamilton and Whillans2002). In January 1999, a 12-cm diameter and 3-m-long pole was sunk 1.5 m into the snow and secured to the surface (see Figure 4). Subsequently, the pole was measured by GNSS in January 1999, December 2005, December 2007, January 2012, December 2013, January 2014 and January 2020.
Diagram of a ‘coffee can’ or submergence velocity measurement system (adapted from Hamilton and Whillans, Reference Hamilton and Whillans2000) and the DCCC density-depth profile measured on the core recovered from the 43 m bore hole of 12 cm diameter; each core recovered diameter and mass were measured.
The long-term vertical velocity component was derived from GNSS measurements of a marker that was anchored in January 1999 at a depth of 42.3 m in firn with a density of 635 kg m−3. This marker, which was frozen at the bottom, is a piece of steel rod secured to the surface by a 5 mm diameter non-stretchable steel wire cable. This piece of rod is held in place by a plastic bag filled with water, which is dropped into the borehole and allowed to freeze at the bottom at −55°C (see Figure 4).
The ‘long-term’ change in ice sheet thickness was calculated using the formula provided by Hamilton and others (Reference Hamilton, Spikes and Stearns2005) (see Figure 4), which uses snow accumulation rates derived from ice core stratigraphy, including atomic bomb markers (1955–1998), the Tambora volcanic eruption (1815 AD) and the Salamas eruption (1259 AD). This calculation applied to the DCCC site indicates the long-term behaviour at the measurement locations on the ice sheet, as both velocity and accumulation measurements represent extended timescales (102–103 years) (Hamilton and others, Reference Hamilton, Spikes and Stearns2005).
Typically, the most significant contributors to the error budget are density and accumulation rate uncertainties (Hamilton and Whillans, Reference Hamilton and Whillans2000). Long-term and current accumulation patterns are confirmed by several measurements. Furthermore, the density profile demonstrates significant variability at the surface, rapidly decreasing to a relative minimum of variability at a mean density of 600–650 kg m−3 (Frezzotti and others, Reference Frezzotti2005; Hörhold and others, Reference Hörhold, Kipfstuhl, Wilhelms, Freitag and Frenzel2011; Gautier and others, Reference Gautier, Savarino, Erbland, Lanciki and Possenti2016).
Since 1996, the GNSS network data has been processed using the Bernese GNSS software (BSW52) with a classical double difference strategy. The utilization of GPS signals exclusively ensured consistency between the various campaigns and receivers employed. The geodetic reference frame has been established using a classical regional network within the International Terrestrial Reference Frame 2014 (ITRF2014), which includes three permanent IGS GNSS stations located on the Antarctic Plate and in operation since 1996: these are Casey (CAS1), Mawson (MAW1), and McMurdo (MCM4). This approach was used to estimate the relative ice movements at DC with respect to the bedrock (see Figures 2 and 3). In subsequent campaigns, loosely constrained daily solutions were computed using the 1996 a-priori coordinates of the reference stations. This methodology differs from the initial calculations described by Vittuari and others (Reference Vittuari2004), where the comparison of pole positions between 1996 and 1999 was based on the presence of two points positioned using the DORIS technique.
Following the installation of the GNSS permanent station DCRU in 2005, the study of its movement in the global reference frame was refined, and its movement relative to the bedrock was estimated, considering the plate rotation models. The DCRU station was utilized as a local reference point to integrate the 2012–2014 measurements into the ITRF2014, in a manner analogous to the utilization of the 1000 and DOMEC stations as reference points in 1996 and 1999, respectively. The availability of a continuous time series also facilitated an examination of the evolution of the 3D velocity of DCRU over 14 years (see Figure 5).
Absolute horizontal and vertical displacement of the DCRU Concordia permanent GNSS station between 2005 and 2019. (a): North component; (b) East component; (c): detrended horizontal North–East component (along the ice flow direction); (d) vertical component; (e): detrended vertical component. The grey line represents the estimated seasonal model, while the red lines on either side indicate the estimated 1σ uncertainty (weighted mean square error, WMSE) of the noise within the model, calculated using the realistic sigma option implemented within the TSVIEW package developed for the GAMIT-GLOBK (MIT) GNSS analysis software (Herring, Reference Herring2003). This approach employs an autocorrelated noise model for the time series, rather than assuming white noise, thus avoiding the assumption of temporally independent errors. The green horizontal lines represent the bounds of three times the WRMS scatter of the detrended residuals.
The estimated maximum uncertainty between two GNSS horizontal positions was found to be up to 0.8 mm a−1 for the 1996–2012/14 measurements (see Figure 3; Supplementary Table 1), up to 5.2 mm a−1 for the 1996–1999 measurements, and up to 0.9 mm a−1 for the 1999–2012/14 measurements. Increasing the time interval between measurements reduces the uncertainties in velocity estimates, as the positional uncertainties from each survey are divided by the number of years elapsed between the measurements.
The uncertainty in the horizontal velocity changes was estimated using the variance propagation law considering the number of years between measurements, yielding values between 4.3 and 5.3 mm a−1 for the measurements taken from 1996 to 1999 and 1999 to 2012/14 (see Figure 5; Supplementary Table 1).
4. Results
4.1 Concordia Station motion: absolute versus relative velocities
The surface ice velocities observed at the top of DC are extremely small, measuring only a few centimetres per year. Consequently, it is essential to carefully analyse the measured velocities to study the movements of the ice with respect to the bedrock. To understand the influence of plate rotation on the relative movement of the ice relative to the bedrock, a detailed study of the time series obtained from the permanent station DCRU was conducted. The movement of the permanent GNSS station DCRU relative to a terrestrial reference frame was assessed by estimating the expected plate rotation motion through the calculation of a Eulerian pole. The velocities of several GNSS stations located in Northern Victoria Land (VLNDEF, Victoria Land Network for Deformation Control project) and from various Antarctic GNSS permanent stations situated between 1200 and 2200 km from Concordia, which are part of the IGS network (see Figure 2), were utilized to derive the data.
The 14-year time series of the DCRU station (from 2005 to 2019), incorporating seasonal signals and temporal correlations in the calculations, permitted the estimation of an absolute velocity of 10.4±0.4 mm a−1 at an azimuth of 276°. The North component was measured at 1.01 ± 0.25 mm a−1, while the East component was -10.36 ± 0.27 mm a−1 (see Figures 2 and 5).
The average vertical velocity of the DCRU station (Figure 5) between 2005 and 2019 was found to be -114.3 ± 2.6 mm a−1. During this period, several velocity rates were observed: -144.4±1.1 mm a−1 from 2005 to 2007, -118.2±0.6 mm a−1 from 2008 to 2014, and -102.3±0.4 mm a−1 from 2014 to 2019.
To verify the robustness of the results obtained using the classical double-difference approach, a recalculation of the entire DCRU time series was performed using a recent undifferentiated Precise Point Positioning Ambiguity Resolution (PPP-AR) method. The velocities obtained with the latter method were found to be within uncertainty to those described above: North component of 1.20 ± 0.28 mm a−1, East component of −10.15 ± 0.3 mm a−1, and Up component of −114.55 ± 2.38 mm a−1.
As illustrated in Figure 2, the rock basement beneath DC exhibits a homogeneous horizontal displacement rate induced by the rigid rotation of the Antarctic plate (Aitken and others, Reference Aitken2014; Golynsky and others, Reference Golynsky2018). The absolute motion of the DCRU station can be attributed to two primary components: plate tectonics and ice dynamics. According to the most recent model (Zanutta and others, Reference Zanutta2018), the bedrock movement at DCRU is estimated to be 11.9 ± 1.0 mm a−1 at an azimuth of 173°, relative to ITRF2014. Conversely, the flow velocity resulting from ice dynamics is estimated to be 17.4 ± 0.7 mm a−1 at an azimuth of 317°, calculated as the vector difference between the GNSS absolute movement at DCRU and the estimated bedrock motion due to plate tectonics around the same point.
In addition, the seven-parameter transformation (three translations, three rotations and the scale factor) was applied, epoch by epoch, to the results of the PPP-AR analysis of ∼20 years of data collected by six permanent GNSS stations in Antarctica (see Supplementary Figure 1) in order to validate the hypothesis of linear movement of the continental region of East Antarctica.
4.2 Surface horizontal velocity at Dome C
The horizontal pole displacements over the 18-year interval (1996–2012/14) are presented in Table 1 and Figure 3. GNSS measurements indicate that at 25 km from the summit area, poles along the major axis (NE–SW) are moving at 63–65 mm a−1, while those along the minor axis (NW–SE) are moving at a faster rate of 210–220 mm a−1. On the minor axis, the pole displaying the fastest movement is E19, whose rate is 221 ± 0.7 mm a−1 at an altitude of 3222.54 m and a slope of about 0.45 m km−1. This pole is located on the NW slope of the dome, which exhibits higher surface mass balance (SMB) (see Figure 6). Another pole, designated C18, located on the opposite slope characterized by a steeper slope and a decrease in SMB, has a velocity of 214 ± 0.7 mm a−1 at an altitude of 3221.25 m.
Present SMB in mm we a-1 from GNSS poles (black dots; 1996–2012/14) and SMB from the GLACIOCLIM SAMBA stake network (blue dots; 2004/6-2020). The map also features a snow paleo-accumulation representation based on snow radar data (grey scale ranging from 1602 to 1739) and dash contour lines indicating changes from the present back to 1739 (Urbini and others, Reference Urbini2008). Surface elevation change (in cm) is shown from GNSS measurements taken between 1996 and 2012/14.
Along the ice divide (SW–NE), at a distance of 25 km from the summit, where the slope is gentle (about 0.2 m km−1), the velocity measurements are approximately one-third of those observed on the steeper slopes on the NW and SE sides (NE A18: 63 ± 0.7 mm a−1 at 3228.1 m; SW E13: 65 ± 0.7 mm a−1 at 3229.8 m). A transect oriented NW to SE (from E19–F11 to C13–C18) shows that the NW slope exhibits velocities that are >10% greater than those observed on the “steeper” SE slope. This is highlighted by an elevation differential of ∼1 m at a distance of 25 km from the summit, between the two slopes of the dome.
A comparison of ice velocities during the periods 1996–1999 and 1999–2012/14 indicates changes in both velocity and azimuth (see Figure 7).
Surface ice velocity contours (in mm a−1) for the periods 1995–1999 (green) and 1999–2012/2014 (magenta). The point sizes indicate the changes in horizontal velocity between the two periods, with increasing velocities represented in red and decreasing velocities shown in blue.
The poles located on the SE slope (C19, C14, C17, D11) demonstrate an increase in velocity of about 11–14 mm a−1, while those located along the W–NW–N slope (E15, E18, A10) exhibit a comparable decrease in velocity. On the contrary, the poles in the NE area exhibit negligible change in velocity. The poles located over the summit area, which show low velocities, demonstrate significant azimuth variations, with some exhibiting counterclockwise rotation (e.g., 1000, E12, E10) and others displaying clockwise rotation (e.g., A14, C14, A15, A16).
The observed velocity changes are outside the margin of error due to measurement uncertainty (∼±5.0 mm a−1), and the geographic distribution of all increasing and decreasing horizontal velocity measurements is consistent at each site, with their values exceeding the level of uncertainty (see Figure 7).
4.3 Surface elevation and vertical velocity at Dome C
The surface elevation of the poles (ellipsoidal, GNSS-derived) and the relative distance from the snow surface between 1996 and 2014 allow estimation of changes in snow surface elevation over the 18 year period and the sinking rate of the poles.
The average vertical displacement of the poles from 1996 to 1999 was approximately −282 mm, corresponding to a vertical velocity of −94 mm a−1 (with a standard deviation of 3 mm a−1). The maximum value recorded during this interval was −101 mm a−1 (A18), while the minimum was −89 mm a−1 (A11). From 1999 to 2012, the average displacement increased to −1084 mm. Considering the time interval from 1999 to 2014, the average displacement increases to −1265 mm, resulting in an overall average vertical velocity of −84 mm a−1, with a standard deviation of 4 mm a−1 for the 1999–2012/14 period. The maximum absolute vertical velocity recorded in this interval was −91 mm a−1 (A13) and the minimum was −75 mm a−1 (A16).
Over the entire 18-year period (1996–2012/14), the average vertical velocity was −85 mm a−1, with a standard deviation of 2.3 mm a−1.
The velocity of submersion (measured at both surface and borehole markers) at the DCCC site was recorded seven times between 1999 and 2020 (see Table 2). The calculated rates of thickness change are relatively small, and the formal uncertainties are close to the estimated thickness changes. It has been observed by several authors (Hulbe and Whillans, Reference Hulbe and Whillans1994; Hamilton and others, Reference Hamilton, Whillans and Morgan1998; Hamilton and Whillans, Reference Hamilton and Whillans2002) that errors in measurements of submergence ice velocity are influenced primarily by the determination of temporal variability in snow accumulation rates.
DCCC ‘coffee can’ measurements, with a horizontal velocity of 9.7 ± 1.0 mm a−1 (1999–2014) and a slope of 0.0001 rad
* Source: Gautier and others (Reference Gautier, Savarino, Erbland, Lanciki and Possenti2016).
The rates of snow accumulation in the DC area have been calculated using several methods and time span intervals. These include a 15-year period (2004–2020) from the GLACIOCLIM SAMBA stake farm; a 40-year span based on atomic bomb tritium marker levels (1955–1998); an 180-year period informed by the Tambora eruption (1815 AD); and 750-year period based on the Salamas eruption (1259 AD). These studies provide comparable average snow accumulation rates, with a slight increase of ∼10% observed in recent decades. For periods from 1816 to 1998 and from 1259 to 2011, the accumulation rate was 25 ± 1.3 mm water equivalent (we) per year (Frezzotti and others, Reference Frezzotti2005; Gautier and others, Reference Gautier, Savarino, Erbland, Lanciki and Possenti2016). For the period from 1955 to 1998, the rate was recorded at 26 ± 1.3 mm we a−1 (Frezzotti and others, Reference Frezzotti2005). However, recent findings (Genthon and others, Reference Genthon, Six, Scarchilli, Ciardini and Frezzotti2016; GLACIOCLIM SAMBA) have revealed that the snow accumulation rate at the AWS-US stake farm measurement was 30 ± 1.8 mm we a−1 for the period from 2004 to 2020 using a density of 340 kg m3.
5. Discussion
5.1 Surface mass balance of Dome C
The temporal variability at a secular scale and the spatial distribution of SMB have been assessed using firn cores, GPR transects (Frezzotti and others, Reference Frezzotti2005; Urbini and others, Reference Urbini2008; Verfaillie and others, Reference Verfaillie, Fily, Le Meur, Magand, Jourdain, Arnaud and Favier2012; Le Meur and others, Reference Le Meur, Magand, Arnaud, Fily, Frezzotti, Cavitte, Mulvaney and Urbini2018), and, since 2006, by employing a stake network and Automated Weather Station (AWS) observations integrated with atmospheric reanalysis data dating back to 1979 (Scarchilli and others, Reference Scarchilli, Frezzotti and Ruti2011; Genthon and others, Reference Genthon, Six, Scarchilli, Ciardini and Frezzotti2016). The present spatial variability of SMB, as detectable through the network of buried poles (see Figure 6), is generally consistent with the snow accumulation gradients identified by GPR (Urbini and others, Reference Urbini2008) and with results from three stake farm measurements (Genthon and others, Reference Genthon, Six, Scarchilli, Ciardini and Frezzotti2016).
The variability of snow accumulation and density at small scales, resulting from relief-related spatial variability (Eisen and others, Reference Eisen2008; Ekaykin and others, Reference Ekaykin, Tebenkova, Lipenkov, Tchikhatchev, Veres and Richter2020), is a primary factor contributing to the uncertainty in stake-based SMB determinations. Utilizing a density of 340 ± 20 kg m3 (Frezzotti and others, Reference Frezzotti2005), the GLACIOCLIM SAMBA stake on the North (A13 GNSS Pole) recorded an SMB of 30 ± 1.8 mm we a−1, while the South stake farm (C19 GNSS Pole) reported a value of 26.6 ± 1.6 mm we a−1. Concurrently, the stake farm in proximity to the AWS-US and DCCC exhibited an SMB of 30 ± 1.8 mm we a−1 over the period from 2004 to 2020. The southern stake farm of GLACIOCLIM SAMBA confirms a significant 12% gradient in SMB compared to the other two stake farms located near the ice divide (AWS-US) and 25 km to the North, as measured by GPR.
Furthermore, the GNSS poles established in 1996 in the West, NW and NE sectors were found to be completely or nearly buried during the 2012/2014 survey, whereas the poles from the NE to the South were still protruding above the snow surface by 20-40 cm. Utilizing an average snow density of 340 ± 20 kg m3, an average SMB value of 31 ± 1.8 mm we a−1 was calculated, with a standard deviation of 5.9 mm we a-1. The minimum recorded value was 27.6 ± 1.6 mm we a−1 (C19, South), while the maximum was 35.6 ± 2.1 mm we a−1 (E16, NW) during the period from 1996 to 2012/14.
Individual GNSS poles generally exhibited slightly higher SMB values compared to the average values recorded by the nearby GLACIOCLIM SAMBA stake farm, which consists of 50 stakes. This slightly higher SMB observed at the GNSS pole (12 cm in diameter) compared to the adjacent stake farm poles (2 cm in diameter) could be attributed to the influence of the larger pole diameter on snow-drift driven accumulation.
The interval between 1996 and 2012/14 was subject to analysis (Supplementary Table 2), revealing an average increase in surface elevation of +70 mm over an 18-year period, with a standard deviation of 110 mm, as determined by GNSS measurements of pole top positions and their distances from the surface. The maximum value observed change in surface elevation was +330 mm (E16), while the minimum value was −210 mm (A14).
The comparison between the changes in surface elevation and the SMB measurements derived from GPR indicates a consistent correlation (R 2 = 0.4). The SMB surveys reveal an asymmetric distribution of snow accumulation (see Figure 6), reflecting the elliptical orography of the dome. The SMB gradient is oriented at an angle of approximately 45° relative to the elliptical shape of the dome. Although a slight slope-related divergence in katabatic flow can be detected, the area is generally flat, with an elevation difference of about 10 m or less within 25 km from the summit, resulting in a homogeneous snow surface.
Urbini and others (Reference Urbini2008) observed that DC is a critical region, signifying the transition in accumulation distribution along the East Antarctic ice divide, with a decrease in the accumulation rate from NE to SW directions. Le Meur and others (Reference Le Meur, Magand, Arnaud, Fily, Frezzotti, Cavitte, Mulvaney and Urbini2018) further confirmed that the most striking characteristic is the accumulation pattern, which remains consistent over time, exhibiting persistent gradients; for example, there is a marked decrease from 26 mm we a−1 at DC to 20 mm we a−1 at the 300 km SW end of the profile over the last 234 years on average.
The ‘coffee can’ system has also yielded consistent results, with the signal at a depth of 42.3 m aligning with the long-term (millennial) accumulation rate of 25 mm we a−1. Conversely, the DCCC pole/surface data aligns with the contemporary (decadal) increase in snow accumulation, recorded at 30 mm we a−1, and signifies a modest elevation increase of 70 mm over the span of 20 years (∼3.5 mm a−1). The observed escalation in surface elevation at DCCC is consistent with the data collected from the pole network.
The effect of glacial isostatic adjustment (GIA) of the crust due to glacial loading/ unloading on vertical velocity is negligible (<1 mm a−1). Recent compilations of GIA models indicate that the uplift rate in the DC area is negative (Shepherd and others, Reference Shepherd2018).
The ice divide at DC is characterized by a pronounced gradient in SMB values along the SE-East to NW-West direction (Frezzotti and others, Reference Frezzotti2005; Urbini and others, Reference Urbini2008; Genthon and others, Reference Genthon, Six, Scarchilli, Ciardini and Frezzotti2016). The SE-East region is considered one of the coldest and driest areas on Earth. The southern part of the ice divide is characterized by lower moisture levels and is depleted in water isotopes, as well as a deficiency in the chemistry and dust content of the snow (Proposito and others, Reference Proposito2002; Magand and others, Reference Magand, Frezzotti, Pourchet, Stenni, Genoni and Fily2004), primarily due to the shadow effect on the SE side of the DC ice divide.
This SE-East catchment area (Byrd/Mullock/David) is characterized by very low SMB, averaging between 22 and 37 mm we a−1 (see Table 3), and has the most extensive megadune and wind crust ablation areas (Scambos and others, Reference Scambos2012; Frezzotti and others, Reference Frezzotti, Gandolfi and Urbini2002; Traversa and others, Reference Traversa, Fugazza and Frezzotti2023). On the contrary, the NW-West region exhibits SMB values ranging from 116 to 228 mm we a−1, which is an order of magnitude higher than those in the SE-East region.
Mass balance of the principal glacier systems draining the Dome C area based on the catchment area defined by Rignot and others (Reference Rignot, Mouginot, Scheuchl, van den Broeke, van Wessem Mj and Morlighem2019)
The table includes the following: basin name, total SMB for the basin area, average SMB (1979–2008), and ice discharge (2009–2017) as reported by (*) Rignot and others (Reference Rignot, Mouginot, Scheuchl, van den Broeke, van Wessem Mj and Morlighem2019). The average subglacial melting has been calculated using data from (°) Willis and others (Reference Willis, Pope, Gwendolyn, Arnold and Long2016). (^) The difference in the position of the grounding line between the LGM and the Present is sourced from Bentley and others (Reference Bentley2014).
A similar asymmetry in SMB and gradient distribution has also been observed along the East Antarctic ice divide between Dome Fuji and the EPICA DML (Fujita and others, Reference Fujita2011). In situ measurements of snow accumulation, combined with atmospheric reanalysis, indicate that most snow accumulation occurs during a limited number of significant precipitation events, with prevailing winds generally blowing from the N-NW quadrant to the S-SE. These air masses, which originate from coastal regions, carry and deposit more moisture along their trajectory from N-NW to S-SE (Scarchilli and others, Reference Scarchilli, Frezzotti and Ruti2011; Genthon and others, Reference Genthon, Six, Scarchilli, Ciardini and Frezzotti2016).
5.2 Recent Dome C ice divide dynamics behaviour
A comparison of ice velocities between the periods 1996–1999 and 1999–2014 (see Figure 7) shows a reduction in velocity in the NW-West sector, alongside an increase in velocity in the SE-East sector up to ∼14 mm a−1. In addition, a rotation of the velocity vectors was observed in the summit area, particularly within 6 km of the summit and at 12.5 km (E12 and A16), along the main axes of the dome (see Table 1). On the contrary, the poles located in the northern part of the dome did not exhibit significant changes in either velocity or direction.
The detrended time series of the DCRU station (see Figure 5) over the period 2005 to 2014 shows a similar trend to that of the nearby GNSS poles (E10, E14). By analysing the velocity maps for the intervals 1996–1999 and 1999–2012/14 (see Figure 7), it is possible to estimate the neutral point of the dome’s dynamics and to assess its migration towards the NW, which amounts to about 1 km over a period of 10 years.
The driving stresses, as calculated by Young and others (Reference Young2017) using the latest available digital terrain model (Helm and others, Reference Helm, Humbert and Miller2014), the ice thickness map and measured ice velocities were then compared. The results demonstrate a substantial correlation (R 2 = 0.88) for points showing a decrease in velocity (NW-West sector) and a moderate correlation (R 2 = 0.51) for points showing an increase in velocity (SE-East sector).
However, a comparison of digital elevation models (DEMs) derived from satellite data collected between 1993-1996 (Remy and others, Reference Rémy, Shaeffer and Legrésy1999) and GNSS surveys (Capra and others, Reference Capra, Cefalo, Gandolfi, Manzoni, Tabacco and Vittuari2000) with later (RAMP2; Bamber and others, Reference Bamber, Gomez-Dans and Griggs2009) and more recent (Cryosat-2, Slater and others, Reference Slater2018) models shows no detectable changes in dome shape or elevation.
A 15-year period was considered, during which the vertical velocity at the DCCC site was found to be −83 mm a−1 at the surface and −39 mm a−1 at the depth of the cable (42.3 m). In comparison, between 2005 and 2019, the vertical velocity at the Concordia Station (DCRU) was −114 mm a−1, stabilizing at −102 mm a−1 during the period from 2014 to 2019. A comparison of vertical velocity measurements between stations DCRU (continued) and DCCC indicates that the surface at station DCRU subsided at a rate ∼30–35% faster than that at station DCCC.
The average vertical velocity at DCRU (see Figure 3) exhibited higher values during the 2005–2007 period, with detrended reductions in vertical velocity ranging from ∼26 to 42 mm a−1, in comparison to the 2008–2014 and 2014–2019 periods, respectively. These variations suggest that the occurrence of increased firn compaction phenomena due to the loading of the station structure, which stabilized 2 years after the completion of the station, resulted in a continued high rate of subsidence (−114 mm a−1 vs. −83 mm a−1) relative to the undisturbed surface area at DCCC.
The Concordia Station was built on a 100 m by 50 m firn platform, which was created in 1998 by compacting the snow to a depth of 4 m and raising it about 1.5 m above the surface. Construction of the two main buildings began in 2000 and was completed in 2006.
The density of the platform firn was measured to be in the range of 500–600 kg m3 over a thickness of ∼6.5 m. Snow and firn density measurements carried out in 2008 and 2021 down to a depth of 15 m utilizing firn cores, (personal communication by F. Possenti, 2022) at the platform near the Concordia Station tower showed similar values to those of the station platform at the time of its installation, with densities of 500–600 kg m3 for the upper 6.5 m and decreasing to around 450 kg m3 below. Notably, no detectable alterations were observed in the density profile over a period of 13 years. The density profile obtained at a depth of <6.5 m was highly comparable to the density profile obtained under undisturbed conditions at the DCCC site. The analysis of the DCRU time series demonstrates seasonal variability in the horizontal and vertical velocity components. Maximum vertical displacement occurs in winter and autumn, while a decrease is observed in the summer months (see Figure 5). Conversely, the horizontal velocity shows an antithetical trend, with a minimum in summer, an increase in autumn and, again, a decrease in winter and spring. The seasonal amplitude of vertical displacement is ∼2 cm, whereas the seasonal variation in horizontal velocity is ∼1 cm a−1. This displacement value is consistent with previous estimates for the interior of East Antarctica (e.g., Veldhuijsen and others, Reference Veldhuijsen, van de Berg, Brils, Kuipers Munneke and van den Broeke2023), although it differs with respect to the timing of the maximum and minimum vertical velocity, which is expected to occur during the spring/summer season (Stevens and others, Reference Stevens, Lilien, Conway, Fudge, Koutnik and Waddington2023).
5.3 Glaciological overview of the Dome C catchments
The drainage basins of the DC are characterized by strongly asymmetric glaciological conditions between the SE-East and NW-West catchments, including differences in SMB, subglacial bedrock morphology and melting rates, subglacial hydrological systems, grounding line positions and ice thickness during the Last Glacial Maximum (LGM) (see Figure 1, Table 3).
The marine based nature of DC’s catchments within the Aurora and Wilkes subglacial basins is indicative of their sensitivity to present and past climate perturbations (Wilson and others, Reference Wilson2018; Crotti and others, Reference Crotti2022). Using DEMs from the radar Cryosat-2 (Slater and others, Reference Slater2018) and radar-lidar satellite altimetry (Bamber and others, Reference Bamber, Gomez-Dans and Griggs2009), the drainage systems of DC have been analysed with tools available in the QGIS open-source Geographic Information System (GIS) software. Due to the very shallow slopes, the catchment basins are separated by less than a few metres.
The drainage catchment areas (DB-QGIS) derived from the two DEMs (radar Cryosat-2 and radar-lidar satellite altimetry) using QGIS software tool are found to be highly congruent with the results presented by other authors (e.g., Zwally and others, Reference Zwally, Giovinetto, Beckley and Saba2012). Rignot and others (Reference Rignot, Mouginot, Scheuchl, van den Broeke, van Wessem Mj and Morlighem2019) utilized slope and flow direction data from satellite InSAR measurements (Rignot and others, Reference Rignot, Mouginot and Scheuchl2017) to derive drainage boundaries, which significantly differ for the basins draining the SE area of DC into the Ross Sea (see Figure 1).
There is a significant discrepancy in the drainage basins delineated by Rignot and others (Reference Rignot, Mouginot, Scheuchl, van den Broeke, van Wessem Mj and Morlighem2019) and Zwally and others (Reference Zwally, Giovinetto, Beckley and Saba2012), as well as the basin derived using DB-QGIS, particularly in the David Glacier, Mullock Glacier and Cook Glacier basins (see Figure 1). According to Rignot and others (Reference Rignot, Mouginot, Scheuchl, van den Broeke, van Wessem Mj and Morlighem2019), the David Glacier and Cook Glacier basins extend into the NNE part of DC, whereas the same area in Zwally and others (Reference Zwally, Giovinetto, Beckley and Saba2012) and the DB-QGIS analysis terminates ∼300 km east of DC, corresponding to the Adventure Trench, which is elongated in an SSE-NNW direction. In this region, several authors (Wingham and others, Reference Wingham, Shepherd, Muir and Marshall2006; Carter and others, Reference Carter, Blankenship, Young and Holt2009) have identified large subglacial lakes that have undergone periodic outbursts, while Le Brocq and others (Reference Le Brocq2013) modelled the subglacial water flux draining to the Mullock and Byrd outlet glaciers.
The Adventure Trench is defined by a bedrock gorge with an approximate orientation of SSE-NNW, with a depth of ∼1 km and an ice thickness that exceeds 4 km. Its width is ∼30–40 km, while the surface slope exhibits a westward direction. Rignot and others (Reference Rignot, Mouginot, Scheuchl, van den Broeke, van Wessem Mj and Morlighem2019) propose that the NE part of DC drains to the Cook Glacier, highlighting the main difference at the Adventure Trench’s orographic basin. The difference in DEM elevation between upstream and downstream at the Adventure Trench ice divide is only a few metres (Frezzotti and Flora, Reference Frezzotti and Flora2002).
The 49 GNSS measurements of ice velocity conducted along the ITASE traverse between Terra Nova Bay and DC (see Supplementary Table 3), and from DC to Talos Dome, show a strong correlation in both magnitude and direction of flow with those derived from MEaSUREs and inferred from InSAR data (Mouginot and others, Reference Mouginot, Rignot, Scheuchl and Millan2017). The correlation coefficient, R 2, was calculated to be 0.99 for all data points (see Figure 8). However, the correlation coefficient shows a substantial decrease for velocities below 10 m a−1 (R 2 = 0.72) and further deterioration for velocities below 3 m a−1 (R 2 = 0.10), aligning with the estimation errors reported by Mouginot and others (Reference Mouginot, Rignot, Scheuchl and Millan2017). These findings support the results of Yuande and others (Reference Yuande2018), who conducted a comparative analysis of 71 GNSS velocity measurements with InSAR data along the transect from Zhongshan Station to Dome Argus. The authors noted that the results align when ice surface velocities exceed 5 m a−1.
GNSS versus InSAR MeASUREs surface velocity along the ITASE traverse transect. (a) Velocity values from 0.1 to 714 m a−1. (b) Velocity values below 10 m a−1. (c) Velocity values below 3 m a−1.
The direction of ice flow along the Terra Nova Bay–DC transect corresponds with the maximum surface slope, as determined by a DEM. Higher velocities relative to this general trend were observed at the D4 site (Vittuari and others, Reference Vittuari2004). At the D4 GNSS site (see Figure 1), located ∼5 km downflow from the ice divide (as defined by Zwally and others (Reference Zwally, Giovinetto, Beckley and Saba2012) and the DB-QGIS of David Glacier), the ice flow measured by GNSS was 4.3 ± 0.04 m a−1, with a direction of 116° (Vittuari and others, Reference Vittuari2004). This value is in close agreement with the MEaSUREs value at the D4 GNSS site of 5.27 ± 0.63 m a−1 at a direction of 120°. The differences in ice velocity measured by GNSS and InSAR at nearby sites are within the confidence intervals of the respective measurements. The ice flow data from D4 supports the drainage basin delineation proposed by Rignot and others (Reference Rignot, Mouginot, Scheuchl, van den Broeke, van Wessem Mj and Morlighem2019) rather than the orographic slope basin delineated by the DB-QGIS and Zwally and others (Reference Zwally, Giovinetto, Beckley and Saba2012) at the Adventure Trench. This provides evidence that David Glacier and Cook Glacier drain the NE part of DC. Frezzotti and Flora (Reference Frezzotti and Flora2002) suggested that the ice flow at the Adventure Trench ice divide has changed direction, possibly since the mid-Holocene, due to the encroachment of the Ross Ice Shelf system. The encroachment may be attributed to alterations in the slope profile of the Ross Ice Shelf outlet glaciers, the perpendicular alignment of the subglacial topography, and the presence of water-saturated sediments.
Given the inaccuracy of InSAR in determining velocity direction and magnitude below 3 m a−1, the flow directions closer to the ice divide or dome (within about 100–150 km) are more accurately represented by the DB-QGIS than in the model presented by Rignot and others (Reference Rignot, Mouginot, Scheuchl, van den Broeke, van Wessem Mj and Morlighem2019).
Utilizing a combination of the DB-QGIS and Rignot catchment systems (see Figs. 1 and 7), the drainage pattern reveals that the DC area within 25 km from the summit flows from SW to NE, with D11 draining into the Ross Sea. The area between D11 and B11 contributes to the Byrd Glacier (1250 km along the flow) and Mullock Glacier (1390 km along the flow). The drainage from pole A19 leads to David Glacier (1090 km along the flow). The northern and western parts of DC (from poles A18 to E13) drain into the Southern Ocean. Specifically, the area at pole A18 drains into Cook Glacier, while a limited section between poles A18 and A13 contributes to Moscow Glacier (945 km along the flow, 870 km straight). From pole A13 to pole E13, the flow is directed to Totem Glacier (1160 km along the flow, 915 km straight).
Within the defined boundaries of the DC’s GNSS network area (25 km radius from the summit), most of the drainage is directed to Byrd Glacier (about 47%) and Totem Glacier (about 44%). The remaining three basins—Mullock, David and Cook—collectively account for a minimal proportion of the northern area’s drainage (about 9%). It is important to note that all poles located on the SE slope of DC drain into the Byrd Glacier Basin.
The Byrd–Mullock–David drainage system flows mainly from West to East (see Figure 1) and is generally perpendicular to the primary tectonic features of the Wilkes Subglacial Basin, Concordia Trench, Adventure Trench and Transantarctic Mountains, which trend South–North. The South–North flow directions of the Cook and Totem Glaciers are consistent with the orientation of the Wilkes and Aurora subglacial basins, respectively. The drainage systems originate mainly from bedrock situated below sea level (see Figure 1).
When a surface profile is traced along the drainage flows derived from MeASUREs ice velocity data, significant differences emerge between the Cook–Moscow–Totem basins and the Byrd–Mullock–David basins (see Figure 9; Table 3). The former shows a very regular elliptical equilibrium shape, characteristic of an ice sheet, whereas the Byrd–Mullock–David profile shows a distinct shape with a bend in the central part that occurs between 200 and 700 km from DC. The surface profile exhibits a steeper gradient in the final section closer to the grounding line, a feature that may be attributed to the dynamics associated with the outlet glacier traversing the Transantarctic Mountains. Instead, the velocity profiles for Byrd Glacier exhibit a pronounced increase in velocity at distances between 200 and 300 km from DC, which is not observed in any of the other five velocity profiles (see Figure 9).
Profiles along the flow paths of glaciers draining into the Ross Sea (Byrd, Mullock and David) and the Southern Ocean (Totem, Moscow and Cook). The profiles are derived from BEDMAP2 surface elevation data, ice velocity measurements from MEaSUREs and bedrock elevation from BEDMAP2.
The velocity profile in the Cook–Moscow–Totem region is quite regular and primarily follows the slope. The model-derived subglacial melting rates (Willis and others, Reference Willis, Pope, Gwendolyn, Arnold and Long2016; Le Brocq and others, Reference Le Brocq2013) and the distribution of subglacial lakes (Wright and Siegert, Reference Wright and Siegert2012; Wright and others, Reference Wright2012; Reference Wright, Young, Bamber, Dowdeswell, Payne, Blankenship and Siegert2014; Smith and others, Reference Smith, Fricker, Joughin and Tulaczyk2009) show a significant asymmetry between the NW-West (Cook–Moscow–Totem) and SE-East (Byrd–Mullock–David) catchment areas. The SE-East basin is characterized by an extensive presence of subglacial lakes, and its modelled subglacial melting rate is twice that of the NW-West basin (see Figure 1, Table 3). The ratio of SMB to subglacial melting rate for the Cook–Moscow–Totem catchment area is only about 1%, due to the very low accumulation in the SE-East sector. In contrast, the Byrd–Mullock–David region exhibits a much higher ratio of about 10%.
While the surface ice flow in the NW-West catchment area (Cook–Moscow–Totem) is consistent with potential subglacial meltwater drainage pathways, the subglacial meltwater drainage pathways in the SE-East catchment area are often oriented perpendicular to the surface flow.
Rignot and others (Reference Rignot, Mouginot, Scheuchl, van den Broeke, van Wessem Mj and Morlighem2019) estimated a slightly negative mass balance for the Byrd–Mullock–David system and a negative mass balance for the Cook–Moscow–Totem glaciers. Byrd Glacier velocity exhibited a change between 2005 and 2007, attributed to the drainage of subglacial lakes (Streans and others, Reference Stearns, Smith and Hamilton2008). The temporary acceleration of Byrd Glacier, caused by subglacial flooding, has been assessed by Scheuchl and others (Reference Scheuchl, Mouginot and Rignot2012) and references therein. Despite this increase in speed, the glacier has maintained a relatively stable flow regime over the last few decades.
The velocity of Cook Glacier has shown dynamic fluctuations over the past 45 years, with only minor changes observed in the last decade (Miles and others, Reference Miles, Stokes and Jamieson2018). A brief increase in velocity was observed following a subglacial flooding event, and a more pronounced long-term increase of ∼20% was documented between 1989 and 2001. In relation to its equilibrium speed, Totem Glacier accelerated from 1989 to 1996 and slowed in 2000, thus bringing its ice flow back into equilibrium with accumulation. Thereafter, a period of acceleration occurred, peaking in 2007, after which a period of stability followed. The observed changes in ice flow acceleration appear to be directly linked to changes in ocean temperatures. Specifically, an increase in ocean temperatures has been found to correlate with an increase in ice velocity, while cooler temperatures has been found to correlate with a decrease in velocity (Li and others, Reference Li, Rignot, Mouginot and Scheuchl2016).
The ice discharge of the East Antarctic glaciers that drain into the Ross Sea area has remained remarkably stable, according to observations conducted by Gardner and others (Reference Gardner2018) and Rignot and others (Reference Rignot, Mouginot, Scheuchl, van den Broeke, van Wessem Mj and Morlighem2019). In contrast, the region west of the Law Dome has shown subtle increases in ice flow velocity and discharge, although these signals are close to the detection limit, partly due to larger errors associated with previous radar mosaics for this area. However, larger features such as Totten Glacier and the tributaries of the Moscow University Ice Shelf, which drain the western sector of DC, show localized variations in ice velocity.
Schröder and others (Reference Schröder, Horwath, Dietrich, Helm, Van Den Broeke and Ligtenberg2019) observed that several coastal regions of the Wilkes Coast have undergone significant elevation changes, as determined using multi-mission satellite altimetry collected between 1978 and 2017. In particular a thinning of the Totten Glacier was observed at an average rate of 72 ± 18 cm a−1 at the grounding line and 33 ± 12 cm a−1 near the Cook Ice Shelf, while David and Byrd glaciers showed stable or positive elevation changes.
The current asymmetry of the DC drainage system also reflects the paleo-evolution of the ice sheet between the NW-West sector (Southern Ocean) and the SE-East sector (Ross Sea). During the LGM, grounding ice extended to the continental shelf break off Wilkes Land, but did not reach the western Ross Sea area (Mackintosh and others, Reference Mackintosh2014). The positions of the LGM grounding lines (see Figure 1) along the Byrd–Mullock–David glaciers have shifted northwards between 475 and 875 km in the western Ross Sea when compared to their current positions, with ice thicknesses increasing from 500 to 2000 m. In contrast, the grounding line positions along the Cook–Moscow–Totem glaciers have migrated only 80 to 200 km northwards, with only minimal thickening (Bentley and others, Reference Bentley2014; Golledge and others, Reference Golledge2013).
Finally, the deglaciation history since the LGM differs markedly between the two catchments. The present grounding line was reached relatively suddenly on the Southern Ocean coast (Cook–Moscow–Totem) around 14 ka BP, whereas the grounding line in the Ross Sea area was reached 7-9 ka later, around 7 to 5 ka BP (Golledge and others, Reference Golledge2013; Stutz and others, Reference Stutz2021). Hillebrand and others (Reference Hillebrand2021) suggest that the Byrd and/or Mullock glaciers may have captured about half of the catchment area of the Darwin and Hatherton glaciers during the last deglaciation. Furthermore, Yokoyama and others (Reference Yokoyama2016) indicated that the breakup of the Ross Ice Shelf started around 5 ka BP, with the ice shelf reaching its current configuration around 1.5 ka later, suggesting a reorganization of the catchment area in the eastern part of DC.
The present asymmetry observed in the DC drainage system may also be related with the paleo-evolution of the ice sheet between the NW-West sector (Cook–Moscow–Totem, Southern Ocean) and the SE-East sector (David–Mullock–Byrd, Ross Sea) since the deglaciation began. Additionally, the differences in the subglacial hydrological systems between the SE-East and NW-West regions may be attributed to the distinct relationships between surface ice flow/orography and subglacial water flow/tectonic structures.
In the NW-West part, the alignment of surface and subglacial flow directions facilitates a more efficient hydrological system. On the contrary, the perpendicular orientation of the surface and subglacial flow in the Byrd/David catchment area promotes the formation of subglacial lakes and contributes to the observed subglacial floods. Urbini and others (Reference Urbini2008) observed that repeated GNSS measurements at Talos Dome have revealed analogous changes in ice velocity associated with dome summit migration, which appear to correlate with changes in accumulation distribution, similar to observed DC spatial accumulation asymmetry.
6. Conclusion
The vertical velocity of the station, as estimated from the measurements collected at the permanent GNSS station at Concordia Station (DCRU) between 2005 and 2019, is −114.3 ± 2.6 mm a−1. In addition, the flow rate due to ice dynamics, evaluated with respect to a fixed part of bedrock, is measured to be 17.4 ± 0.7 mm a−1. The average vertical velocity of Concordia Station has fluctuated over time due to the loading of the station, stabilizing around −102.30 ± 0.4 mm a−1 since 2014. This value is ∼30–35% faster than that of the undisturbed surface at the DC summit (DCCC). The DCRU GNSS data also show seasonal asynchronous variability in both sinking and horizontal velocity, with variations of a few centimetres.
The SMB rate at DC has exhibited an approximate 10% variation over the preceding decade in comparison with earlier centuries. The spatial variability of snow accumulation at the kilometre scale, derived from snow radar data and measurements from poles/stake farms, is greater than its temporal variability (20–30%) observed at the century scale. A significant temporal increase in accumulation has been observed since the mid-20th century, which may correlate with a shift in snowfall patterns reflecting changes in snowfall trajectories (Frezzotti and others, Reference Frezzotti2005; Urbini and others, Reference Urbini2008; Genthon and others, Reference Genthon, Six, Scarchilli, Ciardini and Frezzotti2016). This phenomenon may be attributed to an increased frequency of blocking anticyclones, which have been observed to transport moist air to the higher elevations of the East Antarctic Plateau (Scarchilli and others, Reference Scarchilli, Frezzotti and Ruti2011; Frezzotti and others, Reference Frezzotti, Scarchilli, Becagli, Proposito and Urbini2013).
Cavitte and others (Reference Cavitte2018) reported that the large-scale (hundreds of kilometres) surface accumulation gradient has remained spatially stable over the past 73 000 years, mirroring current modelled and observed precipitation gradients in the region. The surface elevation of the dome has shown slight increases with a vertical velocity of ∼3.5 mm a−1. Concurrently, the horizontal ice velocity field of the dome has migrated north-westward at a rate of ∼100 m a−1. This migration is accompanied by an increase in velocity on the southeastern slope of the dome and a comparable reduction in the northwestern sector. This movement brings the neutral dynamic position of the dome closer to the EDC core site than it was in 1996.
The observed migration of the dome can be attributed to two key processes: the snow accumulation gradient at the dome/ice divide (short-term) and/or the ongoing dynamics of the catchment area of Byrd Glacier, which has been undergoing reorganization since the middle Holocene in the SE-East sector draining into the Ross Sea. The results do not provide sufficient information to determine whether the current migration of DC is primarily influenced by surface mass balance, dynamic adjustment, or a combination of both. Additionally, the temporal response remains uncertain, as it is unclear whether it is driven by long-term (millennial) or short-term (centennial-decade) factors. In January 2025, the Beyond EPICA team successfully reached the bedrock at a depth of about 2800 metres at LDC. Preliminary field analysis of the retrieved ice core samples appears to represent a continuous climate record of ∼1.2 million years—and possibly even older. Chung and others (Reference Chung2024) used a 2.5D inverse model to determine the age–depth profile along a flow line from DC to LDC that is assumed to be stable in time. Surface velocity measurements are used to determine the flow line and ascertain the flow tube width, allowing the model to account for lateral divergence. The results show that the deepest ice at LDC is predicted to be about 1.12 Ma old, which is in close agreement with the initial results from the preliminary analysis conducted on site. The findings of this paper suggests that even the inner regions of East Antarctica are undergoing substantial changes over decadal to millennial time scales, and highlight the limited perturbations of the climate signal on the ice core due to glaciological variability as per the dome position, with some minor discrepancies restricted in the timing and duration of climate events along the ice core compared to those revealed by geological dating. Martin and others (Reference Martín, Hindmarsh and Navarro2009) observed that the formation of the Raymond arch at DC requires stable ice divide position for a period of approximately 144 000 years. The relatively high variability observed at DC may account for the lack of detected Raymond bumps, as indicated by extensive radar surveys of the area (Cavitte and others, Reference Cavitte2021, and data therein).
These data will be valuable in future assessments of changes in mass balance and ice thickness in these areas, as well as for detecting potential effects of climate change on ice cores. Horizontal and vertical velocity measurements at Dome C provide the longest and most accurate record of ice dome flow. These measurements also offer supplementary information for comparison with satellite observations for mass balance assessments and numerical models.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.28.
Acknowledgements
This research was conducted within the framework of the Concordia Station and the EPICA-EDC projects, and it was financially supported by the PNRA and IPEV. The present study was made possible by the joint French–Italian Concordia Program, which established and operates the permanent station at Dome C. We acknowledge GLACIOCLIM SAMBA Project, (https://glacioclim.osug.fr/-Antarctique-) for sharing Dome C dataset.
Author contributions
L.V., S.G., C.R., M.F. and S.U. interpreted and analysed the radar GNSS data. L.V., L.M., S.G., S.U. and A.Z. participated in survey design and data acquisition. Data processing was conducted by L.V., L.M., S.U. and A.Z. All authors contributed to discussions on depth uncertainty. L.V. and M.F. prepared the manuscript with contributions from all co-authors.
Financial Support
This work was supported by PNRA and IPEV. In addition, M.F. acknowledges support for this research from the Department of Science, Roma Tre University (MIUR-Italy Dipartimenti di Eccellenza 2023-2027).
Competing Interests
The authors declare that they have no conflict of interest.
Open access
