1. Introduction
Knowledge of the internal structure of polythermal glaciers is fundamental to the study of their dynamics. Ground-penetrating radar (GPR) is a suitable tool for determining such structure. This paper focuses on the analysis of the internal structure of two polythermal glaciers located in the Antarctic Peninsula (AP) region. The climate in this region has experienced significant changes over recent decades. During the second half of the 20th century, the AP experienced one of the largest warming trends on the planet, with a rate of $\sim 0.5\, ^\circ$C per decade at Faraday-Vernadsky station (Turner and others, Reference Turner, Lachlan-Cope, Colwell and Marshall2005). Most glaciers in the area retreated over that period (Rau and otherts, Reference Rau2004; Cook and others, Reference Cook2016). This warming period was followed, particularly in the northern AP and the South Shetland Islands (SSI), by a brief but sustained cooling period lasting from the late 20th century to the mid-2010s (Turner and others, Reference Turner2016; Oliva and others, Reference Oliva2017), in turn followed by a return to warming conditions (Carrasco and others, Reference Carrasco, Bozkurt and Cordero2021). Climate variations have a direct effect on the hydrothermal structure of glaciers and hence on their dynamic response (Rabus and Echelmeyer, Reference Rabus and Echelmeyer2002; Pettersson and others, Reference Pettersson, Jansson, Huwald and Blatter2007; Gusmeroli and others, Reference Gusmeroli, Jansson, Pettersson and Murray2012). Regional climate also drives surface mass balance (SMB), through both accumulation and ablation processes. Snow accumulation in this region has shown an increasing trend since the 1930s (Medley and Thomas, Reference Medley and Thomas2019). Surface melting and runoff, in turn, are fundamentally governed by temperature changes and have therefore followed the above-mentioned temperature trends (Costi and others, Reference Costi2018).
Hurd and Johnsons are two polythermal glaciers (Navarro and others, Reference Navarro2009; Sugiyama and others, Reference Sugiyama2019) located on Hurd Peninsula, Livingston Island, which is the second-largest island of the SSI archipelago (Fig. 1). Polythermal glaciers have also been identified at other locations within the SSI, such as King George Island (Wen and others, Reference Wen1998; Breuer and others, Reference Breuer, Lange and Blindow2006; Blindow and others, Reference Blindow2010). A wealth of glaciological studies are available for these glaciers due to their proximity to the Spanish Antarctic station Juan Carlos I (JCI). These span from mass balance observations (Molina and others, Reference Molina, Navarro, Calvet, García-Sellés and Lapazaran2007; Navarro and others, Reference Navarro, Jonsell, Corcuera and Martín-Español2013) and modeling (Jonsell and others, Reference Jonsell, Navarro, Bañón, Lapazaran and Otero2012) to glacier dynamics observations (Machío and others, Reference Machío, Rodríguez-Cielos, Navarro, Lapazaran and Otero2017; Sugiyama and others, Reference Sugiyama2019) and modeling (e.g. Otero and others, Reference Otero, Navarro, Martin, Cuadrado and Corcuera2010), through other subjects such as geomorphology (Ximenis and others, Reference Ximenis, Calvet, Garcia, Casas and Sb́bat2000; Molina, Reference Molina2014) or front position changes (Rodríguez and Navarro, Reference Rodríguez and Navarro2015). Fundamental to our purposes are the local studies focused on radar and seismic applications (Benjumea and others, Reference Benjumea, Macheret, Navarro and Teixidó2003; Navarro and others, Reference Navarro, Macheret and Benjumea2005, Reference Navarro2009). A further interest in Hurd and Johnsons glaciers is that their SMB records are part of the World Glacier Monitoring Service (WGMS) database since the hydrological year 2002.
GPR has been widely used for the investigation of glaciers, including their ice-thickness distribution, internal structure (layering, deformation structures), thermal regime, ice properties (density, water content) and ice-bedrock interface characteristics (presence of water, bed roughness), among others (Schroeder and others, Reference Schroeder2020). In the case of polythermal glaciers, knowledge of the internal distribution of cold and temperate ice layers is critical for thermo-mechanical modeling because of the distinct rheological properties of both types of ice (Greve and Blatter, Reference Greve and Blatter2009). Such distribution is defined by the so-called cold-temperate transition surface (CTS). Multiple studies have demonstrated the potential for detecting thermal regimes within glaciers through radar surveys (e.g. Eisen and others, Reference Eisen, Bauder, Lüthi, Riesen and Funk2009; Gusmeroli, Reference Gusmeroli2010; Schannwell and others, Reference Schannwell2014). Although several studies on the numerical modeling of Johnsons and Hurd glaciers have been carried out (Martín and others, Reference Martín, Navarro, Otero, Cuadrado and Corcuera2004; Otero and others, Reference Otero, Navarro, Martin, Cuadrado and Corcuera2010; Molina, Reference Molina2014), all of them consider the glacier as an isothermal ice mass. Having available the current spatial distribution of cold and temperate ice layers would allow a full thermo-mechanical modeling of the evolution of these glaciers using models such as those by Greve (Reference Greve1997); Aschwanden and others (Reference Aschwanden, Bueler, Khroulev and Blatter2012).
GPR is also a useful tool for determining the thickness of the snow layer. Such data, when referred to the end-of-winter layer, are essential for accurately estimating the accumulation component of the SMB (e.g. Sold and others, Reference Sold2013). GPR allows snow thickness measurements to be made much faster than conventional methods (snow probing, snow pits, shallow cores), which are usually very time-consuming (e.g. Harper and Bradford, Reference Harper and Bradford2003; Grabiec and others, Reference Grabiec, Puczko, Budzik and Gajek2011; McGrath and others, Reference McGrath2018; Laska and others, Reference Laska, Grabiec, Ignatiuk, Uszczyk and Kuhn2019).
The above considerations motivate the dual focus of this paper, consisting of determining: (1) the thickness distribution of the cold ice layer of Hurd Glacier; and (2) the spatial distribution of the end-of-winter snow depth of Hurd and Johnsons glaciers. This will be achieved through the use of very high and ultra high frequency (VHF and UHF, 200 and 750 MHz respectively, in our case) GPR data collected during the 2003, 2008 and 2016 field campaigns.
2. Geographical setting
Our study area is Hurd Peninsula (62° 39–42’S, 60° 19–25’W; Fig. 1), Livingston Island, SSI, Antarctica. It is covered by an ice cap with an area of ~13.5 km2 and spans an altitude range from sea level to ~370 m a.s.l. This ice cap consists of two main basins, Hurd and Johnsons glaciers, plus two additional, smaller basins draining to the southeast that have seldom been studied because they are heavily-crevassed icefalls. Hurd Glacier (4.03 km2) is a land-terminating glacier, with three main lobes: Argentina, flowing northwestwards, Las Palmas, flowing westwards and Sally Rocks, flowing southwestwards. The later has an average surface slope of about 3°, while Argentina and Las Palmas tongues are much steeper, around 13°. Hurd Glacier is located between 250 and 330 m a.s.l. Hurd Peninsula ice cap has three additional tongues flowing eastwards. These tongues are highly crevassed, preventing in situ field measurements. Johnsons (5.36 km2) is a tidewater glacier, flowing north-westwards, that terminates at a calving front 50 m in height extending 570 m along the coast. Just a few meters are submerged. Johnsons Glacier typically has surface slopes between 6° and 10°.
Both Hurd and Johnsons glaciers have a polythermal structure showing an upper layer of cold ice, several tens of meters thick in the ablation zone. In the various lobes of Hurd Glacier, the ice thickness gradually decreases to zero at the front, and the cold ice layer reaches the bedrock (Navarro and others, Reference Navarro2009).
The average ice thickness in 1999–2001 of the Hurd-Johnsons ensemble, determined from GPR data, was 93.6 ± 2.5 m (Navarro and others, Reference Navarro2009). Maximum thickness values of Hurd Glacier are ~200 m in the accumulation zone, while those of Johnsons Glacier reach ~160 m. The ice surface velocity of Johnsons Glacier increases towards the ice front, attaining values up to 65 m a−1 (Otero and others, Reference Otero, Navarro, Martin, Cuadrado and Corcuera2010; Machío and others, Reference Machío, Rodríguez-Cielos, Navarro, Lapazaran and Otero2017). The highest ice velocities for Hurd Glacier are just about 5 m a−1 at its central basin, decreasing towards the frozen-to-bed terminal zones (Machío and others, Reference Machío, Rodríguez-Cielos, Navarro, Lapazaran and Otero2017).
The ice cap on the Hurd Peninsula is exposed to the maritime weather of the western AP. The warmest temperatures in Antarctica have been recorded in the AP (Francelino and others, Reference Francelino2021). The annual average temperature at JCI station during the period 1994–2014 was −1.2 $^\circ$C, with average summer and winter temperatures of 1.9 $^\circ$C and −4.7 $^\circ$C (Bañón and Vasallo, Reference Bañón and Vasallo2016). On Hurd and Johnsons glaciers, the mass gain is mostly by snowfall and wind-driven snow redistribution (Navarro and others, Reference Navarro, Jonsell, Corcuera and Martín-Español2013). The predominant winds on Hurd Peninsula, measured at Johnsons Glacier automatic weather station, are from the north-northeast and the southwest (Fig. 2). Since 2002, the seasonal and annual SMB of both glaciers have been measured using the glaciological method, with average values over the period 2002–2021 very close to zero (Zemp and others, Reference Zemp2021).
3. Field data and processing
For this work we used data from the austral summer campaigns 2003/04, 2008/09 and 2016/17, collected using different radar systems as detailed in Table 1. The higher vertical resolution of the 750 MHz radar is suitable for the analysis of the shallowest structures, namely the snow cover. The VHF 200 MHz radar, in turn, provides greater penetration, at the expense of a lower resolution, and allows analysis of deeper structures, such as the CTS. Common offset profiling was carried out in all cases. In the 2003 and 2016 campaigns, the radar equipment was installed on a sledge towed by a snowmobile. In 2008 campaign it was towed by a skier. Radar trace positioning was achieved by means of GNSS measurements. In the 2003 and 2008 campaigns, single-frequency GNSS was used, implying a horizontal positioning accuracy of ~3 m, but occasionally up to ~5 m. Differential GNSS measurements (centimetric accuracy) at each profile endpoints served as an additional control of the positioning (GPR profiles are carried out at nearly-uniform speed). In these campaigns, the positioning by single-frequency GNSS often provided non-reliable coordinates for the profiles (e.g. radar profiles in a zigzag), so we opted for assigning to the traces coordinates based on the accurately-positioned endpoints and the assumption of a regular traction speed. In the 2016 campaign, differential GNSS was used, involving centimetric accuracy for the positioning. The GNSS-measured coordinates were directly inputted to the GPR recording unit.
In all three campaigns, fieldwork also included snow probing (i.e. measuring the snow depth using an avalanche probe) in the network of ca. 50 mass balance stakes deployed on both glaciers plus 50 other locations, making a total of about 100 measurements. These data were used primarily for the snow depth map but are also useful to determine the thickness of each layer in the cases of the snow-firn column and snow-cold ice column, for the accumulation and ablation zones, respectively. These measurements were complemented by snow pits at 5 locations (Fig. 1). The snow pits were measured twice per campaign, close to the beginning and to the end of the summer season. This allows estimation of the changes in density of the snowpack (including internal ice layers and lenses) during the melting season.
Radar data processing was similar in all cases. The time-zero was adjusted according to the timing of the first arrivals. The DC offset was corrected using a bandpass filter, with cutoff frequencies 0.5 and 1.5 times that of the central frequency. To reduce the disturbing effects from the air and ground (snow/ice surface) direct waves, as well as the antenna ringing effects, we used background removal. Amplitude scaling was applied to remove attenuation effects with depth. Migration was used to move the position of the reflectors to their true positions and to collapse the diffractors, increasing the spatial resolution.
4. Methods
The time-to-depth conversion was done using
where v is the radar wave velocity (RWV) and τ is the two-way travel time from the transmitter to the reflector and back to the antennas, which is the field-recorded quantity. From the thickness measurements, corresponding maps were created using Surfer software and ordinary kriging interpolation.
4.1 Seasonal snow depth
The 2016 end-of-winter snow cover of Hurd and Johnsons glaciers was obtained by combining GPR-derived snow thickness and snow probing data. Note that, in snow studies, snow depth denotes the total height of the snowpack measured in the vertical direction, while the thickness of the snowpack is measured perpendicular to the surface slope (Fierz and others, Reference Fierz2009). If the surface slope is small, both quantities are close to each other. When using GPR to measure the snowpack things become more complex. The wave path of the signal received upon reflection at the bed is not dictated by the perpendicular to the surface but by the perpendicular to the bed at the reflecting point. Therefore, if snow surface and snow bed are not parallel, this introduces distortions in the radar signals that could be corrected through procedures such as migration and topographic correction. Fortunately, in most cases snow surface and base are nearly parallel and, under such conditions, what we get from the GPR measurements is the thickness (not depth) of the snowpack. In our case study, the surface slopes are gentle, and therefore snow depth and snow thickness are very similar, which justifies combining measurements from snow probing and GPR data.
Each trace of the GX750 GPR spans 274 ns, discretized into 2636 samples taken every 0.1 ns. This time window is sufficient to sample the whole column of the seasonal snow cover.
The depth-averaged RWV in snow was determined from the two-way travel time at 48 points in the radar profiles close to the locations of the snow-thickness measurements by snow probing. This RWV was used for the conversion from two-way travel time to snow thickness at all radar profiles.
A fundamental problem for the determination of the end-of-winter snow depth, no matter whether from GPR or snow probing data, is the presence of multiple internal ice layers and lenses within the end-of-winter snowpack. A sample radargram illustrating these difficulties is shown in Fig. 3. The availability of ca. 50 accumulation/ablation stakes across Hurd and Johnsons glaciers allows us to determine a lower bound of the end-of-winter snow depth. This bound is simply obtained as the difference between the stake readings at the end-of-winter and at the end of the previous summer. Actually, this difference is in itself the definition of the winter accumulation layer thickness. The reason for considering it as a lower bound is that Juan Carlos I station is operated only during the austral summer. The last measurements in the campaign are therefore done before station closing, which usually takes places a couple of weeks before the actual end of the summer melting. Snow probing is also made next to these 50 stakes, allowing us to validate the additional 50 snow probe measurements done at locations different from those of the stakes. Further validation data is provided by the 5 snow pits mentioned in the Field data and processing section. Summarising, snow pits, stake readings and snow probing data help in identifying the snow-firn/ice interface from the GPR data.
4.2 Cold ice layer thickness
The cold ice layer of Hurd Glacier was analyzed using profiles taken in the 2003/04 and 2008/09 campaigns. All of the profiles collected in the 2008/09 campaign were used, whereas, for the 2003/04 campaign, only those profiles for the ablation zone of Hurd Glacier were employed (14 out of the 88 shown in Table 1). Since its thickness was measured from GPR profiling at the glacier surface, it actually includes the thickness of the thin overlying seasonal snow cover (~1.2 m on average for the 2008/09 campaign). In certain areas, the time window used for the 2003/04 profiles (~430 ns) was insufficient to detect the CTS. This problem is absent in the 2008/09 data, in which a time window of ~1100–1500 ns with a sampling period of 1 ns was used. For the RWV of cold ice a standard value of 168 m μs−1 was used (Macheret and others, Reference Macheret, Moskalevsky and Vasilenko1993). Note that the cold ice layer is only present in the ablation zone. Therefore, we are just neglecting the higher RWV velocity of a thin (1–2 m) snow layer. Consequently, only the zones with a very thin cold ice layer will be affected by our assumption of a constant RWV. The CTS was manually determined as the transition between the clean, upper part of the radargram devoid of internal reflections/diffractions (cold ice), and the lower part with plenty of diffractions caused by the liquid water in temperate ice (Fig. 4).
4.3 Error estimation
The uncertainty in depth of the GPR measurements was derived from equation (1) using error propagation:
e denotes the error in the variable indicated by the corresponding subscript.
The error in RWV for the seasonal snow cover was taken as the standard deviation mentioned earlier. This error was found to be within the vertical resolution of the instrument, which is theoretically λ/4, with λ the wavelength, but can in practice reach up to λ/2 (Navarro and Eisen, Reference Navarro and Eisen2009). To construct the snow-depth map (Fig. 5), point values were spatially interpolated using kriging. Due to the abundance of quite homogeneously distributed measurement points, we took the interpolation error as the standard deviation of the leave-one-out cross-validation (a single value for the whole glacier).
In the case of cold ice, the error in RWV was taken from Navarro and others (Reference Navarro2009). Measurement errors for the cold ice thickness higher than the vertical resolution were obtained in zones of deep ice, reaching values up to 1 m. The determination of the CTS involves interpretation errors due to the diffuse boundary between the cold and temperate ice layers. We estimated this interpretation error following the technique described in Gusmeroli (Reference Gusmeroli2010), whereby 25 sets of pickings of the CTS were performed by three different researchers. The CTS point value was taken as the average of the 25 picks (shown as red line in Fig. 4), and its standard deviation (blue-shadowed band in Fig. 4) was assigned as error $e_{H_{p}}$. Measurement and interpretation errors were combined as
The construction of the cold ice thickness map was also based on a kriging interpolation. However, as the source data comes from GPR profiling, there is a high data density along the radar profiles, while there exist large data gaps between the profiles. Therefore, we followed the procedure for computing the interpolation error, and combining it with $e_{H_{d}}$, described by Lapazaran and others (Reference Lapazaran, Otero, Martín-Español and Navarro2016).
5. Results and discussion
5.1 Seasonal snow depth
The mean of the depth-averaged RWV in snow determined at 48 points (see Methods section) was v = 210.0 m μs−1, with a standard deviation σ = 5.3 m μs−1, which is within the expected range of values for such media (e.g. Moore and others, Reference Moore1999). This velocity value is typical of dry snow. This is consistent with the fact that our measurements were done before the onset of summer melting. The presence of thin internal ice layers and lenses within the snowpack reveals that some surface melting with subsequent meltwater percolation and refreezing happened, likely in high-temperature peaks during the late spring. However, the winter snowpack was mostly dry when the end-of-winter snow pits were measured, which justifies the mentioned velocity typical of dry snow. Furthermore, Recio-Blitz and others (Reference Recio-Blitz, Navarro, Otero, Lapazaran and Gonzalez2018) reported a mean end-of-winter density of 510 ± 14 kg m−3 at Hurd and Johnsons glaciers over the period 2004–2016. If we use this value with Robin's equation $\varepsilon = ( 1 + 0.851 \rho ) ^2$ (Robin and others, Reference Robin, Evans, Bailey and Bullard1969), relating density ρ and permittivity $\varepsilon$ of dry snow, and then apply the relationship $v = c/\sqrt {\varepsilon }$ (Dowdeswell and Evans, Reference Dowdeswell and Evans2004), with c the radar wave velocities in free space, we get a radar wave velocity in snow of 209.2 ± 1.8 m μs−1.
The 2016 seasonal snow depth maps, based on manual snow probing and snow probing combined with GPR are shown in Figures 5(a,b) respectively. The average snow depth from snow probing measurements was 1.44 m ± 0.09 m (the quoted error is the standard deviation). The largest values were found in the higher-elevation areas of the glaciers. The maximum snow depth, of 2.45 ± 0.21 m, was found in the northern part of Johnsons Glacier, where snow redistribution by wind accumulates a large amount of snow (see Fig. 2) . The interpolation error associated with the snow depth map is 0.42 m.
The distribution of snow accumulation on both glaciers depends on both meteorological and orographic conditions. The altitude range of these glaciers is in the order of 300 m, so there is a moderate difference between the amount of snow initially deposited at high and low altitudes. The linear fit of snow depth vs elevation shown in Fig. 6 spans a range of ~1 m. However, the glacier hypsometry implies that most of the glacier area is concentrated in a relatively narrow altitude range (150–300 m), and therefore most data points are clustered in the mentioned elevation range, for which the snow depths vary between 1.0 and 1.6 m. Consequently, the variation of snow depth with elevation is not a critical issue. On the other hand, snow redistribution by the prevailing north-north-easterly and south-westerly winds (Fig. 2) causes slightly higher accumulations on Johnsons Glacier, and in general in the concave parts of both glaciers. In particular, the south-westerly winds remove by erosion part of the snow initially deposited on the upper reaches of Hurd glacier and deposit it on the more concave Johnsons Glacier. This results in systematically more positive winter (and annual) surface mass balances in Johnsons Glacier as compared with Hurd Glacier (Navarro and others, Reference Navarro, Jonsell, Corcuera and Martín-Español2013). On the other hand, higher temperatures at lower elevations imply a higher degree of compaction of the snow, and hence a thinner snow layer. Compaction alone, however, is insufficient to explain the difference in snow depth between the upper and lower elevations. Some additional factors have to be taken into account. Rain should in principle be discarded, as it normally occurs much later in the season. However, some surface melting and runoff at the lowermost part of the ablation zone should not be discarded; another intervening process could be the erosion of snow away from the glacier at the lowermost elevations. Unfortunately, for most of these processes, we lack observational evidence, so we cannot quantify their influence.
The data recorded with the GX750 in 2016 present a large amount of noise and stratification, mostly due to refreezing of percolating surface meltwater. Although the field measurements were carried out in the early austral summer, before the onset of strong surface melting, some melting and subsequent percolation and refreezing of meltwater happen in the glaciers in this region during the late spring (Recio-Blitz and others, Reference Recio-Blitz, Navarro, Otero, Lapazaran and Gonzalez2018). The wavelength (λ) of our 750 MHz radar in snow (RWV = 210 m μs −1) is 0.28 m, so its vertical resolution is between 0.07 m (λ/4) and 0.14 m (λ/2). Ice layers and lenses observed in the snow pits are in general thinner than this (typically ca. 1–3 cm). However, the radar wave actually integrates the energy reflected from all discontinuities in physical properties within a spatial range equal to its resolution. At the glaciers under study, quite often adjacent thin ice layers are separated by distances lower than the radar resolution and can therefore be detected, though not resolved individually. Also thin individual layers will be detectable when they have strongly altered electrical properties. However, this is not the usual case in glacier environments. Because of such noise, the interpretation of the seasonal snow layer would not have been possible without the additional information provided by the data analysis of stake readings, snow probing, and snow pit measurements. In fact, in our case study it was often impossible to track the GPR internal layers (in particular, the snow–ice and snow–firn transition surfaces) between manual point measurements (snow pits, stake readings, snow probing). This is illustrated in Fig. 7, where we can see that, although close to the snow pits, the snow–firn interface can be identified with the help of the snow depth determined at the pit, but when we move away from the pit such an interface cannot be longer tracked. The central panel of Fig. 7 shows a section of a radargram at an intermediate distance between two snow pits, and we are unable to identify which is the layer corresponding to the snow–firn interface. In fact, the distribution of layers and lenses of refrozen ice has a high spatial variability even at short-scale distances. Note that there is often no correspondence between the ice layer location determined in the snow pits (shown as blue lines in the two lateral panels of Fig. 7) and the occurrence of internal reflections in radargrams taken just 5–10 m from the snow pits. Consequently, the snow-thickness product that includes the GPR data does not improve significantly on previous observations from stake measurements such as that by Navarro and others (Reference Navarro, Jonsell, Corcuera and Martín-Español2013). This contrasts with other studies in which high-frequency radar profiling added important information (Harper and Bradford, Reference Harper and Bradford2003), or even constituted the main source of data, with snow pits/probing providing calibration and validation data (Grabiec and others, Reference Grabiec, Puczko, Budzik and Gajek2011).
Calibrating the snow layer thickness identified in the radargrams with snow depth from the probing measurements could result in small differences between the snow depth maps shown in Figures 5a,b. The map incorporating the GPR data (Fig. 5b) provides additional detail due to greater data coverage.
5.2 Cold ice layer thickness
The polythermal structure of Hurd Glacier, first suggested by Navarro and others (Reference Navarro2009), can be inferred from Fig. 8a, which shows the thickness of the cold ice layer, while Fig. 8b shows the error in thickness. The cold ice layer is only present at lower elevations, basically corresponding to the ablation zone. Its average thickness is 29.1 ± 1.5 m (this error is the mean of the individual measurement errors), and the highest values are found along the central flowline in the lower ablation zone of Sally Rocks tongue, where the thickness of the cold layer reaches its maximum value of 80.8 ± 2.5 m. The ice below the firn layer in the accumulation zone is temperate, and the thickness of the corresponding layer gradually decreases towards the glacier snouts in all lobes of Hurd Glacier. The ice in the frontal zones is cold, meaning that the glacier is frozen to bed. The purple lines in Fig. 8a indicate where the CTS reaches zero value, meaning that the glacier is frozen to bed from those lines to the glacier margins. This distribution of the cold and temperate ice layers is characteristic of so-called Scandinavian type glaciers, typical of the Scandinavian and Svalbard regions (Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012), but seldom reported in sub-Antarctic islands such as the South Shetland Islands. In this archipelago, in addition to our studied glaciers in Livingston Island, polythermal glaciers have been reported or suggested in King George Island, either from borehole temperature data (Wen and others, Reference Wen1998), from modeling studies (Breuer and others, Reference Breuer, Lange and Blindow2006), or from GPR data (Blindow and others, Reference Blindow2010). Other glaciers and ice caps in this archipelago with a similar morphology, particularly those ending on land, could be expected to be polythermmal. Having the ice frozen to bed in the terminal zone of a land-terminating glacier has strong implications for glacier dynamics. The largest velocities found in the central part of the glacier sharply decrease when approaching the terminus, resulting in strong compression, often accompanied by reverse faulting (Molina and others, Reference Molina, Navarro, Calvet, García-Sellés and Lapazaran2007).
The errors in cold ice layer thickness shown in Fig. 8b, calculated combining the method by Lapazaran and others (Reference Lapazaran, Otero, Martín-Español and Navarro2016) and the interpretation error, are lower along the locations of the radar profiles, as a consequence of the interpolation error, except in the middle part of the central flow line, where the maximum CTS interpretation error occurred. The dashed line in the figure indicates the zero depth of the CTS, which approximately corresponds to the location of the equilibrium line. As the ice below the firn is temperate, the cold ice layer thickness is zero to the east of the dashed line. Consequently, the errors in cold ice thickness in that zone are meaningless and we have blanked them.
The largest values of the cold layer thickness of Hurd Glacier have been confirmed by thermistor chain measurements in boreholes. In Fig. 8 we have indicated the location of a borehole drilled in January 2015, where a string of 32 thermistors, spaced 2 m apart, was installed. They were first measured one year later, to ensure that the thermistor temperatures were in equilibrium with the neighboring ice. All temperatures were below 0 °C. Unfortunately, the CTS depth was larger than expected for this zone, so the lowermost thermistor, at a depth of 60.55 m below the ice surface, did not reach melting temperature; it measured ~−0.15 °C. Taking into account the readings of the above thermistors, the CTS could be located at a depth between 60 and 70 m. This CTS depth is larger than that measured in a borehole drilled in Johnsons Glacier in January 2015, where a CTS depth lower than ~50 m was found (at another borehole, closer to the glacier calving front, only temperate ice was identified) (Sugiyama and others, Reference Sugiyama2019).
The temperature of the cold ice measured at 15 m depth at Hurd borehole (located at 116 m a.s.l.) was of −1.1 °C. Corresponding values for ice cores at King George Island recovered in the summer of 1991/92 were −1.9 °C (site at 100 m a.s.l.) and −1.4 °C (site at 150 m a.s.l.) (Wen and others, Reference Wen1998). For comparison with other regions, these values are similar or slightly warmer than those measured at comparable elevations on Scandinavian-type glaciers in Svalbard, where e.g. 15-m depth temperatures at East Grönfjordbreen, at sites between 160 and 225 m a.s.l., provided values of ~−2°C (Chernov and others, Reference Chernov, Vasilyeva and Kudikov2015); at Hansbreen, values were between −1.6 and −1.8 °C at a site at ~160 m a.s.l., and between −2.4 and −3.1 °C at another site at ~240 m a.s.l. (Jania and others, Reference Jania, Mochnacki and Gadek1996). The reason for the slightly lower borehole temperatures in Svalbard is probably that the mean annual air temperatures in Svalbard (e.g. −5.2 °C at Ny-Alesund and −4.6 °C at Svalbard Airport over the period 1981–2010; Førland and others, Reference Førland, Benestad, Hanssen-Bauer, Haugen and Skaugen2011) are colder than those at the South Shetland Islands (e.g. −2.2 °C at Bellingshausen Station, King George Island, over the period 1986–2015; Oliva and others, Reference Oliva2017). Nevertheless, in both cases the warmer temperatures at 15-m depth as compared with the mean annual temperature probably reflect the importance of meltwater refreezing on the release of latent heat (Jania and others, Reference Jania, Mochnacki and Gadek1996).
6. Conclusions
From the above discussion, the following conclusions can be drawn:
1. Although GPR has been shown elsewhere to be an efficient tool for determining the end-of-winter snow cover thickness with wide coverage (Harper and Bradford, Reference Harper and Bradford2003; Grabiec and others, Reference Grabiec, Puczko, Budzik and Gajek2011; McGrath and others, Reference McGrath2018; Laska and others, Reference Laska, Grabiec, Ignatiuk, Uszczyk and Kuhn2019), in our case study it hardly provided any further information to that given by manual snow probing. The reason was the abundance of internal reflections by multiple snow layers and lenses formed from refreezing of percolating meltwater, which prevented the identification of the snow–ice/firn interface without the help of the snow-probed depths. The set of around one hundred snow probing points alone provided a satisfactory snow depth map.
2. The available data suggest that the distribution of the end-of-winter snow cover of Hurd and Johnsons Glaciers (with average and maximum thicknesses of 1.44 ± 0.09 and 2.45 ± 0.21 m, respectively, in 2016/17) is governed by a combination of snow redistribution by the prevailing winds, the greater compaction rates at lower elevations, the erosion of snow away from the glacier and possibly some surface melting and runoff, the two latter processes at the lowermost part of the ablation zone.
3. The spatial distribution of the cold ice layer thickness of Hurd Glacier (with average and maximum measured values of 29.1 ± 1.5 and 80.8 ± 2.5 m, respectively) reveals a Scandinavian-type polythermal structure. Such a structure has often been reported at other locations (e.g. Scandinavia and Svalbard), but seldom in Antarctica. This polythermal structure has important dynamical and geomorphological implications in terms of compressional regime at the glacier snout and associated reverse faulting. The polythermal structure could be a more common feature of South Shetland Island glaciers and ice caps than suggested by the scarce literature on the matter, in particular for land-terminating ice masses.
Due to the climate evolution of this region (Turner and others, Reference Turner2016; Oliva and others, Reference Oliva2017), it would be of interest to repeat some of the GPR soundings reported in this paper, in particular those regarding the polythermal structure. This would allow analysis of the evolution of the CTS in response to regional climate changes.
Data
All raw and processed data are directly available from the authors upon request.
Acknowledgements
This research was funded by grant PID2020-113051RB-C31 from Agencia Estatal de Investigación. UL was supported by grant PRE2021-100049 financed by MCIN/AEI/10.13039/501100011033 and by FSE+. We thank L. Copland (scientific editor) and two anonymous reviewers for their many suggestions, which greatly improved our manuscript.