The role of snowmelt, glacier melt and rainfall in streamflow dynamics on James Ross Island, Antarctic Peninsula

The Antarctic Peninsula is experiencing a rapid increase in air temperature, which has a major impact on the entire ecosystem, including the runoff process. Understanding the water balance dynamics in such a fragile and sensitive environment is therefore crucial. Although water availability plays an important role in polar ecosystems, runoff generation in the Antarctic Peninsula region is still poorly understood. We analysed the variability in rain, snow and glacier contributions to runoff in relation to climate variability in a small, partly glaciated catchment on James Ross Island in the north-eastern Antarctic Peninsula. We used the hydrological model HBV to simulate the runoff process from June 2010 to May 2021 at a daily resolution. The model was calibrated against both measured discharge and glacier mass balance. Model simulations showed the negative mass balance of Triangular Glacier for 9 out of 11 study years, with an average annual mass loss of 49 mm water equivalent. About 92 % of the annual runoff occurred between October and May. On average, peak runoff occurred in the second half of the summer season due to the combination of glacier melt and snowmelt. The majority (76 %) of runoff originated from snowmelt, 14 % originated from glacier melt and 10 % from rainfall. Higher total annual runoff occurred in warmer and drier years associated with high glacier melt runoff. The contribution of snowmelt to total runoff was higher in colder years with more precipitation. Our simulation showed the presence of significant runoff-generating events outside the usual high summer runoff measurement season.

1 Introduction

The Antarctic Peninsula region has experienced a rapid increase in air temperature during the second half of the 20th century (Vaughan et al., 2003; Turner et al., 2005). Although the warming of this region was interrupted in the early part of the 21st century (Turner et al., 2016; Oliva et al., 2017), the cooling period ended in the mid-2010s (Carrasco et al., 2021); since then, the region has experienced several warm events with record high temperatures (González-Herrero et al., 2022; Gorodetskaya et al., 2023). Regardless of short-term climate oscillations, future climate projections indicate that air temperature is expected to rise significantly in both the near and far future (Bozkurt et al., 2021; Zhu et al., 2022).

Increasing air temperature leads to a shift in the phase of precipitation from solid to liquid (Vignon et al., 2021) and accelerates snow cover melting (Abram et al., 2013), glacier mass loss (Cook et al., 2005; Vaughan, 2006; Chuter et al., 2022; Seehaus et al., 2023) and changes in the permafrost active layer (Hrbáček and Uxa, 2020; Kaplan Pastíriková et al., 2023). These changes have a significant impact on the runoff regime of catchments in proglacial areas (Gooseff et al., 2007; Nowak et al., 2021). Streams in proglacial environments represent a key driver for changes in the entire ecosystem, as they form the landscape through fluvial erosion and sediment transport (Rosa et al., 2014; Kavan et al., 2017). Streamflow variations and water availability have a major influence on the evolution of terrestrial ecosystems (Gooseff et al., 2017), and changes in the amount of freshwater and sediment released into coastal seas can affect marine ecosystems (Meredith et al., 2018; Braeckman et al., 2021). Available liquid water is a necessary condition for the existence of life in extremely cold and dry environments, and shifts in the hydrological cycle introduced by climate change may significantly alter local fragile ecosystems and organisms (Gutt et al., 2021). It is therefore essential to quantify water balance components and identify their relative contribution to the overall hydrological regime of the environment.

In Antarctica, runoff is mainly studied in the McMurdo Dry Valleys (Chinn and Mason, 2016; Gooseff et al., 2022). In this area, an increase in the flow season duration has been observed (Gooseff et al., 2007). However, climatic conditions in the McMurdo Dry Valleys differ significantly from those in the Antarctic Peninsula region. Annual precipitation in the McMurdo Dry Valleys is less than 100 mm (Doran et al., 2002; Fountain et al., 2010), several times lower than the modelled amount of 300–700 mm on James Ross Island in the north-eastern Antarctic Peninsula (van Wessem et al., 2016). In addition, rainfall events occur in this region, unlike the Dry Valleys, where only snowfall has been observed (Fountain et al., 2010).

In the Antarctic Peninsula region, runoff has been studied on King George Island, James Ross Island, the South Orkney Islands and Vega Island. Nonetheless, these studies covered only one summer season and focused mainly on sediment transport (Rosa et al., 2014; Hodson et al., 2017; Sziło and Bialik, 2017; Kavan et al., 2017, 2023; Stott and Convey, 2021; Kavan, 2021). A limited number of studies have focused on the runoff generation process (Moreno et al., 2012; Lyons et al., 2013; Lee et al., 2020). The most complex study in this field so far was conducted by Jung et al. (2022), who described the different contributions of old and new water during warm and cold periods on King George Island based on the chemical and stable isotopic water composition. Falk et al. (2018) analysed the variability of surface runoff in a small partly glaciated catchment on King George Island during one summer season using a simple hydrological model developed by the authors. Kavan et al. (2017) analysed the runoff variability of two streams on James Ross Island for the 2014/15 summer season. Their results indicate that the runoff variability of these streams was mainly driven by air and ground temperatures, with maximum runoff occurring at the end of January. The subsequent study (Kavan, 2021) focused primarily on fluvial sediment transport but also described runoff variability for the 2017/18 summer season.

Although changes in the runoff process affect both terrestrial and marine ecosystems, as runoff variability plays an essential role in erosion, sediment transport and water availability (Foreman et al., 2004; Kohler et al., 2015; Hodson et al., 2017), runoff generation in the Antarctic Peninsula region is still poorly understood. This is partly due to the severe limitations of direct measurements of individual components of the precipitation-runoff process, which can be made only during the short austral summer season and are often subject to large measurement errors (Tang et al., 2018; Seefeldt et al., 2021). Therefore, the objectives of our research were (1) to reconstruct streamflow in a small, partly glaciated catchment on James Ross Island, Antarctic Peninsula, and (2) to assess the inter-annual variations in rain, snow and glacier contributions to streamflow in relation to recent climate variability. The study area belongs to the largest deglaciated area in the Antarctic Peninsula, with ongoing changes in water source contributions to streamflow affecting the seasonal runoff distribution, which further impacts terrestrial ecosystems. To our knowledge, studies on the runoff process in the Antarctic Peninsula region have been conducted only for individual summer seasons. Therefore, the inter-annual variability of the runoff process has not been investigated. To fill this gap, we used a conceptual hydrological model coupled with a glacier mass balance model to perform runoff simulations of individual water balance components from June 2010 to May 2021 based on available climate, glacier mass balance, snow depth and streamflow data.

2 Data and methods

2.1 Study area

The analyses were performed in the Triangular Glacier catchment, which is situated in the Ulu Peninsula in the northern part of James Ross Island (Fig. 1). The catchment area covers 2.8 km², and the elevation ranges from 70 to 588 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ The catchment is situated on the Lachman Crags plateau, its south-western slopes and an adjacent low-elevated plain forming part of the Abernethy Flats to the west. The highest part of the catchment spreads over the Lachman Crags plateau (∼ 400 to 600 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), which is separated from the rest of the catchment by up to 300 m high cliffs (Engel et al., 2023). The lower part of the catchment is covered by glaciers, ice-cored moraines and braid plains (Jennings et al., 2021). A small, shallow ice-marginal lake is located near the terminus of Triangular Glacier.

Figure 1 Location of the study catchment in the northern part of James Ross Island (Matsuoka et al., 2018; Czech Geological Survey, 2009).

The catchment area was delineated from the digital terrain model produced by GEODIS (Czech Geological Survey, 2009) for the northern part of James Ross Island at an 8 m spatial resolution. Almost one-third of the catchment (28 %) is glaciated, and 19 % of the catchment is covered by Triangular Glacier (GLIMS Glacier ID: G302151E63856S), which represents the smallest glacier on James Ross Island (Engel et al., 2023). The south-eastern part of the catchment (9 %) is covered by the northernmost tip of the Lachman Crags ice cap (GLIMS Glacier ID: G302158E63874S).

The mean annual air temperature measured on the north coast of James Ross Island near the Johann Gregor Mendel Station (Fig. 1b) was −6.7 °C for the period 2011–2021 (Ambrozova et al., 2019; Kaplan Pastíriková et al., 2023). Estimated annual precipitation ranges from 300 to 700 mm of water equivalent (van Wessem et al., 2016). Data from three stations located on the northern tip of the Antarctic Peninsula showed an increase in precipitation from 1970 to the early 1990s, followed by a decrease from 1991 to 1999 (Carrasco and Cordero, 2020). From the mid-1990s to the mid-2010s, snowfall increased and rainfall decreased due to the cooling period (Carrasco and Cordero, 2020). Snow depth in the lowland plains of James Ross Island is typically less than 30 cm (Hrbáček et al. 2016). In contrast, seasonal changes in snow depth varied from 30 to 90 cm in Triangular Glacier (Engel et al., 2023). This uneven accumulation of snow is largely due to the prevailing strong south-westerly winds, which have a significant impact on both snowfall distribution and snow redistribution by wind (Kňažková et al., 2020; Kavan et al., 2020). According to Engel et al. (2023), snow drift accounts for a 43 % reduction in snow depth in Triangular Glacier.

The evolution of Triangular Glacier followed the changes in air temperature. From 1979 to 2006, the glacier lost almost one-third of its volume. After that, the rate of retreat decreased for almost a decade. Since 2015, the retreat rate has increased again (Engel et al., 2023).

2.2 Data

The semi-distributed bucket-type hydrological model HBV – in its software implementation HBV-light version 4.0.0.25 (Seibert and Vis, 2012; Seibert and Bergström, 2022) – was used to simulate the runoff process for the period 1 June 2010 to 31 May 2021. The date of 1 June was estimated as the start of a water year. Mean daily air temperature and total precipitation were used as input to the model. Daily runoff and seasonal glacier mass-balance series were used for model calibration (see below).

Air temperature was measured at the Johann Gregor Mendel Station (Fig. 1b) using the Minikin TH data logger and EMS33 temperature sensor (EMS, Brno, Czechia), with an accuracy of ± 0.15 °C. A linear regression was used to adjust the time series based on air temperature measurements performed in the middle part of Triangular Glacier during the 2017/18 summer season.

Direct measurement of precipitation, mostly occurring as snowfall, is very difficult and inaccurate in the study area due to strong winds, which limits the usage of standard approaches, e.g. a heated rain gauge. Short-term precipitation measurements were taken in January and February 2022 at the Johann Gregor Mendel Station using a Thies laser precipitation monitor (disdrometer) and a manual rain gauge.

Due to incomplete and inaccurate measured data, daily precipitation was simulated by the Weather Research and Forecasting model 4.3 (WRF; Skamarock et al., 2019). The model was run in a two-nested-domain setup with a horizontal resolution in the inner domain of 5 km. In the vertical coordinate, 65 eta-levels were used for the simulation of atmospheric dynamics. The WRF model was forced by ERA5 reanalysis (Hersbach et al., 2020). The WRF model has been successfully applied in previous studies on James Ross Island (Matějka et al., 2021; Matějka and Láska, 2022). Following these validation studies, the model used the 3DTKE boundary layer scheme (Zhang et al., 2018), Thompson microphysics scheme (Thompson et al., 2008) and NoahMP (Niu et al., 2011) as a land surface scheme dealing with interactions between the surface and the atmosphere. Both shortwave and longwave radiation were parameterised with the RRTMG scheme (Iacono et al., 2008). The ability of WRF to simulate precipitation reliably has been confirmed by validation against measured data (Matějka et al., 2022).

The time series of runoff for the period 8 February–15 March 2018 was calculated from automatic water level measurements (in 10 min intervals; station indicated with a yellow dot, Fig. 1b and c) using a hydrostatic pressure sensor (DipperLog F100/M30, Heron Instruments, CA) and the rating curve derived from manual velocity measurements (Flowtracker Handheld Acoustic Doppler Velocimeter, SonTek, USA).

The annual surface mass balance of Triangular Glacier for the years 2014/15–2019/20 was determined by the glaciological method (Kaser et al., 2003) reported by Engel et al. (2023). The initial volume of Triangular Glacier at the start of the model simulation (2009/10) and calibration (2014/15) periods was calculated based on the ice surface digital elevation model (February 2006) and the ground-penetrating radar survey of 25 January 2017 (Engel et al., 2023). To account for the change in glacier extent between 2005/06 and 2009/10, we adopted the value of −0.3 % a⁻¹ reported by Engel et al. (2023) as the mean annual retreat rate between 2005/06 and 2014/15.

Snow depth was measured using a sonic distance sensor (Judd Communication, USA) fixed at the automatic weather station located in the central part of Triangular Glacier during the period 6 February 2017 to 23 January 2020. The zero point level was chosen arbitrarily (Engel et al., 2023). The snow depth was used only to calculate the snow water equivalent (SWE) as the observed snow depth multiplied by a constant density of 389 kg m⁻³, which represents the mean measured density in the summer season 2022. The estimated SWE time series was used only for model evaluation.

Global solar radiation was measured using the CMP11 pyranometer (Kipp and Zonen, the Netherlands) near the Johann Gregor Mendel Station from 1 January 2015 to 10 March 2018. The global radiation time series was used only for the correlation analysis between climate and runoff data.

2.3 Hydrological model HBV coupled with glacier routine

For runoff simulations, the HBV model was used (Seibert and Vis, 2012). The study catchment was divided into sub-zones based on elevation intervals of 100 m for glacier-free areas and 50 m for glaciated areas. The model inputs are formed by air temperature, precipitation and potential evapotranspiration (PET). The temperature-based method defined by Oudin et al. (2005) was used for the calculation of daily PET using the observed air temperature measured at the Johann Gregor Mendel Station. Available measured data of glacier mass balance and runoff were used for model calibration.

The HBV model consists of five routines: (1) a degree-day method-based snow routine, which controls snow accumulation and melt, (2) a degree-day method-based glacier routine, which simulates glacier mass balance, (3) a soil moisture routine, which calculates daily groundwater recharge and actual evaporation simulated as functions of daily PET and the actual soil moisture, (4) a response routine in which runoff from two groundwater boxes is calculated and (5) a routing routine that simulates the propagation of runoff through the catchment by applying a triangular weighting function. A more detailed model description is given by Seibert and Vis (2012).

The HBV model partitions snowfall and rainfall using the concept of a single threshold temperature (T_T). T_T is one of the calibrated model parameters. Because the catchment is divided into elevation zones, the precipitation phase is calculated separately for the specific elevation zone using a calibrated temperature lapse rate. The snowmelt is calculated using a degree-day method, which additionally accounts for snow water holding capacity and potential refreezing of meltwater. The most important parameter of this method is the melt (degree-day) factor, which was calibrated by the model. The snow routine does not account for either sublimation or evaporation from snow cover.

The ice water equivalent (WE) was calculated using the Δh parameterisation method that links glacier area and thickness and allows for the representation of the glacier advance. This method is implemented in the model and is described in detail by Seibert et al. (2018). The ice melt is simulated only in the case where there is no snow on the glacier using the degree-day approach with a calibrated melt factor. The melt water is added to the liquid water content of the glacier. The outflow is then computed for each elevation zone (Seibert et al., 2018).

The snow redistribution function was set up in the model, which redistributes snow above a certain amount from elevations above a defined threshold to lower elevations. Specifically for the study catchment, snow accumulations above 1000 mm that occurred above 525 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ were redistributed evenly over the lower zones. A small section of the Lachman Crags ice cap located in the highest part of the catchment was not considered in the model. This part of the catchment was considered as glacier-free in the simulation.

The tracking algorithm implemented in the model and described by Stahl et al. (2017) was used to determine the contribution of rain, snowmelt, and glacier meltwater to runoff. The separation method is based on the complete mixing of the three components in a “virtual mixing tank” that runs parallel to the system of individual storages and fluxes of the HBV model. Unlike many other separation methods, the method applied in this study considers rain falling on snow or glaciers as the rain component (in case it does not refreeze). Similarly, melting snow on glaciers is considered as the snow contribution to runoff (Weiler et al., 2018; Jenicek and Ledvinka, 2020).

The model was calibrated against both runoff (available only for the 2017/18 summer season) and glacier mass balance using a genetic algorithm procedure (Seibert, 2000) for the period 2017/18–2019/20 with the warming-up period 2014/15–2017/18. Before the calibration, the upper and lower limits of all calibrated parameters were applied to ensure that the parameters are physically relevant. Two objective functions were used for calibration: (1) Kling–Gupta efficiency KGE (Gupta et al. 2009; 70 % weight), which combines the evaluation of the bias of mean runoff, flow variability and flow dynamics, and (2) mean absolute error of glacier water equivalent (30 % weight). To reduce the issue of parameter equifinality, a median simulation calculated from 100 calibration runs was used for further analysis. The model performance was additionally evaluated using an approximately 3 year long time series of estimated SWE. During this evaluation, the mean simulated SWE in the elevation zone, in which the weather station is located, was compared to the estimated SWE.

2.4 Runoff and climate characteristics and data analysis

Several runoff and climate characteristics were used in the analysis. Mean daily or annual air temperature (T), precipitation (P), SWE, annual sum of positive air temperature (T_positive) and average daily global solar radiation (GR) were used to describe climate variability (Table 1). All of the above characteristics represent mean catchment values. Additionally, glacier water equivalent (WE) was analysed separately for each elevation zone. Runoff variability was represented by daily, monthly and annual runoff (Q) and by rainfall, snowmelt and glacier melt runoff components (Q_rain, Q_snow and Q_glacier, respectively). Spearman's correlation coefficient (r_s) was calculated to analyse the relationship between runoff and climate characteristics.

Table 1 Runoff and climate characteristics used in the analyses.

3 Results

3.1 Evaluation of the model performance

The values of the objective functions resulting from the model calibration showed acceptable performance for runoff and good performance for glacier mass balance. The median values of 100 calibration runs were 0.33 and −11.9 mm for the Kling–Gupta efficiency and glacier absolute mean error, respectively (Fig. 2). For selected objective functions not used for calibration, the median values were 0.32 for the Nash–Sutcliffe efficiency, 0.64 for the runoff volume error, 0.83 for the runoff Spearman rank correlation coefficient and 0.98 for the glacier water equivalent (absolute mean relative error).

Figure 2 The objective function values for 100 calibration runs. Functions used for calibration of the model are shaded in grey.

Figure 3 (a, b) Simulated versus observed runoff (Q) in the Triangular Glacier catchment for the period 8 February–15 March 2018. (c) Simulated and observed mass balance. Observed mass balance with estimated uncertainty was taken from Engel et al. (2023). (d) Mean simulated SWE in the elevation zone, in which the weather station is located, and the estimated SWE, which was calculated as the observed snow depth on Triangular Glacier (with a zero point level that was chosen arbitrarily (Engel et al., 2023)) multiplied by a density of 389 kg m⁻³, which was the mean measured density in the summer season of 2022.

A comparison of the observed and simulated runoff showed that the model simulated the daily variability of the runoff relatively well, although peak runoffs were generally underestimated (Fig. 3a). The simulated glacier mass balance matched well with the observed glacier mass balance for years with available measurements (cumulative difference 7.6 mm of WE for the period 2014/15–2019/20; Fig. 3c). For all years but one (2015/16), the simulated mass balance was within the estimated uncertainty of the measured mass balance (Fig. 3c).

The model was further evaluated by a comparison of the simulated SWE and the estimated SWE. There was a strong correlation between observed and simulated values, although the absolute values differ (Pearson correlation coefficient equal to 0.82; Fig. 3d). The fact that the estimated SWE was not used for model calibration suggests the model's ability to simulate the snow storage sufficiently well.

Table 2 Mean monthly snowmelt (% of annual snowmelt).

3.2 Variability in snow cover and glacier mass balance

The model results showed that the average annual mass loss of Triangular Glacier was −49.7 mm (± 1 mm considering the variability of 50 % of the model simulations given by different calibrated parameter sets). The negative mass balance during the first two water years of the study period was followed by the highest positive mass balance (19 mm) in 2012/13, followed by two years of near-equilibrium mass balance in 2013/14 and 2014/15 (Fig. 4a). From 2015/16 until the end of the study period, the mass balance remained negative. The highest mass loss occurred in the last year of the study period in 2020/21 (−150 mm). The second-highest mass loss was found in 2016/17 (−130 mm).

Figure 4 (a) Cumulative daily changes in glacier water equivalent (WE) from the start of the simulation (cyan line), and changes in WE in individual water years (bars). (b) Daily snow water equivalent (SWE; grey line), and 10 d moving average of the air temperature (negative values in blue, positive values in red).

In terms of the variability of the glacier mass balance with elevation, the lower-middle part of the glacier, at elevations between 150 and 200 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, which is the second-largest part of the glacier by area, experienced the largest absolute decrease in glacier water equivalent, with a reduction of 217 mm. In contrast, the highest parts of the glacier did not lose any mass.

Mean catchment daily SWE ranged from 59 mm on 18 March to 216 mm on 29 September. The highest SWE (382 mm) occurred on 12 January 2013 (Fig. 4b). It was caused by an unusually cold beginning of summer, resulting in almost no snowmelt until late January (see Fig. 7a). The lowest SWE was 9 mm, which occurred on 4 March 2021, associated with very high air temperatures during the warm part of the year (Fig. 4b). Mean monthly snowmelt was highest in November (21.2 % of the total snowmelt; Table 2), followed by January and February (16.1 % and 16.4 % of the total snowmelt, respectively).

3.3 Runoff regime

The simulated mean annual runoff for the study period was 415 mm (± 45 mm). The majority of runoff (92 %) occurred from October to May. The mean annual evapotranspiration was 7 mm, and the mean annual precipitation was 369 mm (Fig. 5). About 76 % (315 ± 39 mm) of the runoff originated from snowmelt, 14 % (58 ± 5 mm) from glacier melt and 10 % (42 ± 7 mm) from rainfall. The flow regime is characterised by two peaks, one in November and the second in February. During the winter, the air temperature was well below the freezing point, and thus almost no runoff occurred. In general, as the air temperature rises at the start of the warm season and causes snow to melt, the contribution of Q_snow to the total runoff increases. In contrast, Q_glacier starts to contribute during summer, as parts of the glacier are often snow-free. Q_rain contributes to the total runoff throughout the warm season, primarily from October to March (Fig. 5). In Figs. 4 and 5, we can see that, throughout the year, Q_snow was the primary contributor to runoff. From December to March, the mean monthly Q_glacier exceeded the mean monthly Q_rain.

Figure 5 Mean daily precipitation (P; blue bars), 3 d moving average of air temperature (T; red line), mean daily rainfall runoff (Q_rain), snowmelt runoff (Q_snow), glacier melt runoff (Q_glacier; stacked areas, left axis), cumulative Q, Q_rain, Q_snow and Q_glacier (lines, right axis) in the Triangular Glacier catchment for the period 2010/11–2020/21.

Figure 6 Mean (bars), minimal and maximal (error bars) monthly runoff (Q), rainfall runoff (Q_rain), snowmelt runoff (Q_snow) and glacier melt runoff (Q_glac) in the Triangular Glacier catchment for the period 2010/11–2020/21.

In order to further analyse the runoff regime, we calculated the monthly values of runoff and its components. Several distinct peaks can be seen in the mean monthly values: Q_snow in November and February, Q_rain in November and January and Q_glacier in February (Fig. 6). The mean monthly runoff was highest in February, accounting for 18.3 % of the annual runoff. The high runoff during this month was primarily caused by the combination of the highest monthly Q_glacier and the second-highest monthly Q_snow. In November, which had the highest mean monthly Q_snow, the second-highest mean monthly runoff (17 % of the annual runoff) was recorded. Additionally, November had the highest mean monthly Q_rain. In contrast, the lowest monthly runoff occurred in August, accounting for only 0.3 % of the annual runoff. The runoff was also very low in June and July (0.9 % and 0.7 % of the annual runoff, respectively). Large variability of the monthly runoff was found between individual years (Fig. 6). The highest variability of monthly Q was observed in December (from 0.5 mm in 2013 to 173 mm in 2012), as a result of the highest variability of monthly Q_snow and Q_glacier occurring. These variations are mainly determined by the timing of seasonal snowmelt, which can occur throughout the warm season based on variations in the air temperature (Fig. 4b).

Table 3 Annual Q, Q_rain, Q_snow and Q_glacier, and the contribution of Q_rain, Q_snow and Q_glacier to runoff. Minimum and maximum values are shown in bold.

3.4 Inter-annual runoff variability

Annual runoff varied between 282 and 501 mm during the period of 2010/11–2020/21 (Table 3). In 2013/14, the lowest annual runoff occurred due to a cold summer (as can be seen in Fig. 3), resulting in the lowest annual Q_rain (17 mm), low annual Q_snow (260 mm), and almost no annual Q_glacier (4 mm). The highest annual runoff was observed in 2019/20 as a result of the second-highest annual Q_snow (375 mm) combined with above-average Q_rain (45 mm) and above-average Q_glacier (79 mm).

Figure 7 Contribution of rainfall runoff (Q_rain), snowmelt runoff (Q_snow) and glacier melt runoff (Q_glacier) to total runoff (Q) during two selected years: (a) the year 2012/13 with the lowest annual glacier melt contribution and (b) the year 2020/21 with the highest glacier melt contribution.

The contributions of rainfall, snowmelt and glacier melt to runoff varied significantly from year to year. The lowest annual Q_glacier occurred in the year 2012/13 (3 mm; Table 3), and the contribution to runoff was 1 % (Fig. 7a). Over this year, the highest annual Q_snow (382 mm) was observed. This was due to the cold beginning and middle of the summer, resulting in the largest snow accumulation during the observed period, which did not melt until the end of the summer (Fig. 4b). In both years 2012/13 and 2013/14, Q_snow contributed the most to the total runoff, accounting for 92 % of it.

The highest annual Q_glacier (151 mm) occurred in 2020/21 (Table 3, Fig. 7b). This was due to the large number of days with high air temperature from November to May (72 d with positive daily air temperature from November to February), which caused the snow cover to melt-out earlier than in other years (Fig. 4b). In 2020/21, Q_glacier accounted for 31 % of the annual runoff. Despite below-average annual Q_rain and Q_snow, the annual runoff was the second-highest (489 mm) in 2020/21. In 2016/17, annual Q_glacier also contributed 31 % to total runoff. This was due to the low snow accumulation during the winter (Fig. 4b), which caused the lowest annual Q_snow (237 mm).

In 2010/11, the highest annual Q_rain occurred (100 mm; Table 3), which accounted for 23 % of the total annual runoff. This is an unusually high amount, as the second-highest annual Q_rain was only about half of the previous year's amount (53 mm in 2011/12).

3.5 Seasonality in peak flows

In addition to examining the inter-annual variability of runoff, we analysed individual extreme runoff events and their causes. Most of the annual peak flows occurred from November to the first half of February. Only in 2013/14 did the annual peak occur at the start of October (Fig. 8). These annual peak flows were caused by different water sources contributing to the total runoff at different times. Annual peak flows were consistently associated with high Q_snow and occurred during the same runoff episode as the annual peak Q_snow, except in 2015/16. In that year, the annual peak runoff in January was due to a combination of the annual peak Q_glacier and high Q_snow, while the annual peak Q_snow occurred in October.

Figure 8 Daily rainfall runoff (Q_rain), snowmelt runoff (Q_snow), glacier melt runoff (Q_glacier) (lines) and maximal values of Q, Q_rain, Q_snow and Q_glacier (triangles) for individual water years.

The annual peak Q_rain typically occurs during October and November in most years. However, in 2012/13, 2014/15 and 2018/19, the annual peak Q_rain occurred in January, February and December, respectively (Fig. 8). The annual Q_glacier peak mostly occurred in February, except for January in 2011/12 and 2015/16 and April and March in 2012/13 and 2013/14, respectively (Fig. 8). However, in these two years (2012/13 and 2013/14), almost no Q_glacier was observed (see Table 1).

3.6 Impact of climate on runoff variability

To identify the climatic factors that show the greatest impact on runoff variability, we conducted a correlation analysis. This analysis of climate and runoff characteristics showed strong positive correlations between annual T_positive and annual Q_glacier (r_s = 0.82, p-value < 0.05; Fig. 9a). A weaker correlation was found between annual T_positive and annual Q (r_s = 0.68, p-value < 0.05) and between annual Q_glacier and annual T (r_s = 0.63; p-values < 0.05). Additionally, annual Q_glacier showed a negative correlation with annual P (r_s = −0.62, p-value < 0.05), and annual Q was similarly correlated with P but with lower statistical significance (r_s = −0.6, p-value < 0.1). In contrast, neither annual Q_rain nor annual Q_snow showed a significant correlation with annual T or annual P. Among the three runoff components, only annual Q_glacier showed a significant correlation with total annual Q (r_s = 0.59, p-value < 0.1), indicating that annual Q_glacier strongly influences the variability of total annual runoff.

Figure 9 Spearman's rank correlation coefficients for the selected meteorological and hydrological characteristics at the (a) annual and (b) daily level. The colours indicate the value of the Spearman's rank correlation coefficient. Only statistically significant correlations are shown; correlations with p-value < 0.1 are shown with a coloured rectangle without a border, and correlations with p-value < 0.05 have a black border. Correlations obtained using daily data (b) were calculated for the period from 1 January 2015 to 10 March 2018, which represents the intersection of the data available for the selected variables.

The correlation analysis for daily values was conducted on a subset of the data (from 1 January 2015 to 10 March 2018) for which measured global radiation was available. This analysis showed that daily Q, Q_rain, Q_snow and Q_glacier were more strongly correlated with daily air temperature (r_s = 0.78, 0.72, 0.73, 0.39, respectively; p-values < 0.05) than with daily global radiation (r_s = 0.44, 0.37, 0.39, 0.12, respectively; p-values < 0.05; Fig. 9b). A significant, although weak, negative correlation was found between daily runoff characteristics and daily precipitation (r_s > −0.15; p-value < 0.05). All daily runoff components were strongly correlated with daily runoff. However, correlations involving Q_glacier were much weaker than those associated with the other runoff components.

Figure 10 Relationship between absolute anomaly of T_positive, relative anomaly of P and relative anomaly of (a) mean annual runoff Q, (b) rainfall runoff Q_rain, (c) snowmelt runoff Q_snow and (d) glacier melt runoff Q_glacier for the years 2010/11–2014/15 (squares) and 2015/16–2020/21 (circles). Note that colour map ranges differ across panels.

To further shape the interpretation of the results of the correlation analysis, we investigated the relationships between the absolute anomaly of T_positive and the relative anomalies of P and Q and the individual runoff components (Fig. 10).

Higher annual runoff occurred during warmer (positive anomaly of T_positive) and drier (negative anomaly of precipitation) years (Fig. 10a). Such years rather occurred in the second half of the study period (circles in Fig. 10). The high annual runoff in these years was associated with the higher Q_glacier (Fig. 10d). However, Q_snow and Q_rain did not show the same behaviour as Q and Q_glacier (Fig. 10b an c). All four years with above-average Q were among the five years with above-average Q_glacier. However, only two (one) years with above-average Q_snow (Q_rain) were also years with above-average Q.

Of the three runoff components, annual Q_glacier had the highest inter-annual variability, with a coefficient of variation (CV) of 0.79. The second-highest CV was for annual Q_rain (0.53). The CVs for Q and Q_snow were much lower (0.15 and 0.16, respectively). Inter-annual variations in total runoff are mostly caused by variations in annual Q_glacier rather than variations in annual Q_snow, even though Q_snow is a dominant component of the total runoff.

4 Discussion

4.1 Mass balance of Triangular Glacier

In their study, Engel et al. (2023) linked the changes in the Triangular Glacier mass balance to air temperature variability. The simulated mass balance of Triangular Glacier is generally consistent with the direct mass balance measurements by Engel et al. (2023), who reported a low melt rate between 2006 and 2014, followed by a period of large mass loss. The simulated mass balance of Triangular Glacier was negative in all water years except for 2012/12 and 2013/14. In water year 2014/15, the simulated balance was slightly negative (−1 mm; Fig. 3). However, if the year would be defined by dates of mass balance measurement, the simulated mass balance would be positive (20.25 mm; Fig. 2c), which is in agreement with Engel et al. (2023). The difference between the simulated and measured mass balance may be due to the location of the glacier on the leeward side of Lachman Crags, which has a significant effect on snow accumulation (Kavan et al., 2020).

The simulated mass balance variability from 2010/11 to 2013/14 corresponds well with the mass balance pattern of the nearby Whisky and Davis Dome glaciers, as described by Engel et al. (2018). Both of these glaciers had a positive balance in these years, except for 2011/12, which was the year with the highest simulated mass loss for Triangular Glacier from 2010/11 to 2013/14. Similar to Davis Dome, the simulated mass balance of the Triangular Glacier was also slightly negative in 2014/15, while Whisky Glacier had a positive balance in that year. However, in contrast to the positive mass balance for Davis Dome and Whisky Glacier, our simulation shows a negative mass balance in 2010/11.

4.2 Runoff regime

According to Kavan et al. (2023), seven studies have been devoted to measuring runoff in the Antarctic Peninsula area between 1991 and 2020 (Inbar, 1995; Kavan et al., 2017; Sziło and Bialik, 2017; Falk et al., 2018; Stott and Convey, 2021; Kavan, 2021; Kavan et al., 2023). For most of the studies, the duration of the measurements did not exceed 2 months. Only one study (Stott and Convey, 2021) measured for 78 d. The measurements were mostly made during January and February and sometimes also in March. Only one study (Stott and Convey, 2021) covered most of December. However, no study took measurements during September, October, November or May. Therefore, our results, which include 11 years of simulated daily runoff, provide valuable information, although they are subject to the uncertainties of modelling based on limited amounts of data.

The results showed that the flow season in the study catchment lasts from September to May, with snowmelt being the main contributor to runoff throughout the year. Based on the mean monthly runoff, the seasonal hydrograph may be characterised by two peaks (Fig. 6). The first peak, occurring in November, was mainly caused by snowmelt. The second and higher peak, occurring in February, was a result of a combination of snowmelt and glacier melt contributing to the runoff. In most years, the highest daily runoff was observed in January or February, and it was associated with peak Q_snow.

The occurrence of maximum runoff at the turn of January and February is in agreement with conclusions by Kavan et al. (2017), who measured streamflow in two proglacial streams near the Johann Gregor Mendel Station during the 2014/15 summer season (Fig. 1b). However, these results are limited to measurements of less than 2 months. Gooseff et al. (2007), based on long-term measurements in the McMurdo Dry Valleys, reported that peak flows occur at the beginning of the flow season. The reason for this finding may be that the total precipitation amount is lower than that for our study area, resulting in lower snow accumulation on glaciers and thus its earlier melt-out, causing the earlier initiation of glacier melt, which is the primary source of runoff. It is important to note that the flow season in the McMurdo Dry Valleys is shorter than that observed in the study area, with an average length of 70 d (Gooseff et al., 2007).

On King George Island, runoff variability is primarily driven by air temperature associated with snowmelt at the beginning of the season (Falk et al., 2018). This is because the climate is warmer, and the region experiences higher rainfall during the summer season compared to our study area. Towards the end of the summer season, runoff variability is driven by rain events. Of the three components, rain contributed the least to total runoff in our study catchment; however, an increase in rainfall frequency and intensity is expected in the future (Vignon et al., 2021).

The results showed high inter-annual variability of runoff mainly associated with glacier melt. The persistent mass loss of small land-terminating glaciers around the Antarctic Peninsula may lead to the depletion of water storage in catchments, and a larger dependency of runoff generation on seasonal snow cover will very likely lead to increased variability of runoff (Huss and Hock, 2018). This, in turn, can affect freshwater ecosystems, as runoff variability has a significant impact on both stream (Kohler et al., 2015; Gooseff et al., 2017) and lake (Obryk et al., 2016; Foreman et al., 2004) habitats by controlling nutrient transport and modulating water levels. In polar arid areas, long-term changes to the ecosystem can be caused by even short-term extreme events associated with high runoff (Fountain et al., 2016).

Glacier melt and associated high runoff also support sediment delivery to the ocean, providing the marine ecosystems with important minerals and nutrients (Hodson et al., 2017). The marine ecosystem is highly sensitive to changes in sediment loads (Sahade et al., 2015). Sediment transport is expected to increase in the near future due to rising air temperatures (Stott and Convey, 2021). There are several studies from the Antarctic Peninsula region analysing the amount of transported sediment, but none of them analysed a period longer than 3 months (Kavan et al., 2023). Our results, showing the variability of the runoff over 11 years, may improve the estimates of the amount of transported sediment.

4.3 Impact of climate on runoff variability

Higher total annual runoff was simulated in warmer and drier years with above-average Q_glacier. Despite being the largest contributor to total runoff, Q_snow showed little inter-annual variability as the seasonal snow cover eventually melts each year.

Inter-annual variation in runoff was primarily influenced by variations in annual Q_glacier, which showed considerably higher variability compared to annual Q_snow and annual Q_rain. This increased variability in annual Q_glacier resulted from its moderate dependence on annual temperature and annual precipitation and strong dependence on annual T_positive. However, it should be noted that the analysis is based on a relatively short period of 11 years, which limited the variability of climate characteristics in individual years, such as sufficient cases of wet and warm years in the study period (see Fig. 10).

The relative contribution of snowmelt to total runoff was higher in years with relatively higher precipitation and lower air temperatures. This outcome was due to the fact that, in these years, more snow accumulated and therefore melted later in the season, leaving the glaciers snow-covered for most of the summer season and resulting in low glacier melt.

The correlation between daily runoff and daily Q_glacier was weaker compared to other runoff components. The reason is that glacier melt in each elevation zone is simulated only in the case of the snow-free glacier. In contrast to snow and rain contributions, this resulted in a reduced number of days in the study period where Q_glacier contributed to total runoff.

Additionally, Q_glacier often follows after Q_snow; therefore, total runoff on a given day includes both Q_glacier and Q_snow from that day (Q_snow comes from higher elevation zones where the snow is still present), as well as Q_snow from previous days, taking into account the response time simulated by the HBV response function. These factors cause the total runoff during glacier melt to be mixed from different water sources, resulting in a lower correlation coefficient between Q_glacier and Q.

Mean daily runoff showed a stronger correlation with air temperature than with global radiation, which is consistent with the findings of Kavan et al. (2017) and Falk et al. (2018). Additionally, Kavan et al. (2017) reported a similarly strong relationship between ground temperature at a 5 cm depth and discharge in glacier-free catchments, which contrasts with the weak relationship between ground temperature and discharge in glaciated catchments reported by Falk et al. (2018).

4.4 Input data and hydrological modelling uncertainty

The modelling of the precipitation-runoff process is subject to uncertainties, which must be considered in the interpretation of the results. Uncertainties are mostly associated with the model structure and with the amount and accuracy of input and calibration data. Our study catchment is located in a very remote, poorly accessible region, where limited data are available for model calibration. For example, streamflow data are typically available only during the austral summer, when the polar station is in operation (Kavan et al., 2017; Kavan, 2021).

It must be noted that the model calibration and simulation were carried out on different data sets, which, except for air temperature, did not cover the whole study period, and in the case of precipitation, it was a fully modelled time series. The streamflow measurements were taken during only one summer season (8 February–15 March 2018) and may not fully represent the variability of streamflow during the whole study period. The glacier mass balance measurements were carried out in 2014–2020, which, in combination with the glacier volume measurements in 2015/16, gave us sufficient information for the model calibration period. However, the volume of Triangular Glacier at the start of the simulation in 2009/10 was estimated from a digital elevation model of the ice surface of Triangular Glacier in February 2006, assuming a constant change in glacier extent of −0.3 % per year, as stated in Engel et al. (2023) as the mean annual retreat rate between 2005/06 and 2014/15. For this reason, the uncertainty of the initial state of the glacier was introduced into the simulation.

Additionally, the streamflow measurements themselves are more uncertain than in more accessible regions because it is not possible to construct a standard gauge profile. Therefore, the measured profile (cross-section) may change during the measurement period, which additionally causes lower accuracy of the respective rating curve. Due to this constraint, the HBV model was chosen for the water balance components simulations, as it has been found, despite several limitations, to be suitable for simulating the runoff process in polar environments (Wawrzyniak et al., 2017; Osuch et al., 2019, 2022). Additionally, it has been repeatedly tested in catchments where only very short or episodic streamflow time series are available or even single measurements exist (e.g. Seibert and McDonnell 2015, Pool et al. 2017). The HBV model structure is rather simple compared to physically based models, but it has been shown to achieve better results than more complex models (Seibert and Bergström, 2022; Girons Lopez et al., 2020). However, a limitation of the HBV for modelling runoff in polar regions is that it does not include soil thermal processes (Bui et al., 2020).

Another possible source of uncertainty is the WRF-modelled precipitation. However, when compared to observations in summer 2021/2022, this model performed very well (Matějka et al., 2022). The model was validated using four nested domains in resolutions 8100, 2700, 900 and 300 m. For all domains, the bias for cumulative precipitation from 2 January 2022 to 26 February 2022 reached −16 % to −1 % with the Spearman's rank correlation coefficient for daily amounts of 0.86 to 0.87. This result suggests a relatively low sensitivity of simulated precipitation accuracy to model resolution in contrast to, e.g. wind speed in complex terrain (Matějka and Láska, 2022). The impact of the WRF model resolution on simulated precipitation in the Antarctic Peninsula region was also discussed by Pishniak and Beznoshchenko (2020).

The above data limitations may be partially overcome by calibrating the model against multiple components of the runoff process, leading to better overall performance (Konz and Seibert, 2010; Finger et al., 2015; Nedelcev and Jenicek, 2021). In the case of our study catchment, the HBV model was calibrated against both runoff and glacier mass balance. In addition, simulated SWE correlated well with estimated SWE (Fig. 2d). All of the above procedures reduced model uncertainty and increased overall model performance and reliability (van Tiel et al., 2020).

The degree-day factor for snow is one of the model parameters; thus, it was calibrated. Before the calibration, the upper and lower limits were applied to ensure that the parameter is physically relevant. After the automatic calibration, a fine-tuning of the degree-day factor was performed to prevent an unrealistic increase in snow storage over the study period (so-called “snow towers”). The final median values resulting from 100 calibration runs were 6.01 $\mathrm{mm}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ (non-glaciated areas) and 2.23 $\mathrm{mm}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ (snow over the glaciers); the median value of TT was −0.21 °C.

Another source of uncertainty in our study is not considering the Lachman Crags ice cap in the hydrological modelling, although a part of the glacier is located in the highest part of the study catchment at elevations above 500 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ Unfortunately, there are no measurements of the mass balance of this glacier. Because we did not consider this glacier in the hydrological model, only snow accumulation and snowmelt were simulated from this part of the catchment. This could lead to an underestimation of Q_glacier. However, we assume the overall influence of the glacier on the total catchment water balance as negligible because our model did not simulate any glacier melt above an elevation of 350 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (results not shown).

5 Conclusions

We simulated water balance and daily runoff in a small partly glaciated catchment on James Ross Island over the period 2010/11–2020/21. Specifically, we analysed the inter-annual variability of glacier, snow and rain contributions to total runoff. The main findings are given in the following points.

The mean annual runoff for the study period from June 2010 to May 2021 simulated by the HBV model was 415 mm, the mean annual evapotranspiration was 7 mm and the mean annual precipitation was 369 mm. The water budget is largely influenced by the negative mass balance of Triangular Glacier (for 9 of 11 years), with an average annual mass loss of −50 mm WE. About 76 % of runoff originates from snow cover, 14 % from glaciers and only 10 % from rainfall.

The majority (92 %) of annual runoff occurred between October and May, with the highest mean monthly runoff occurring in the second half of summer due to a combination of glacier melt and snowmelt. Additionally, high runoff was found in November due to the rapid melting of the seasonal snow cover. In all months, Q_snow was the dominant contributor to runoff.

Of the three runoff components, only annual Q_glacier was significantly correlated with annual T and P. Annual Q_glacier also had the highest variability and was the only component significantly correlated with total annual Q, indicating its importance for inter-annual runoff variability.

The contribution of snowmelt to total runoff was higher in colder years with more precipitation (mostly snowfall), while glacier melt contributed more during warmer years with overall less precipitation.

Overall, the results showed that runoff on James Ross Island is highly variable from year to year. The main driver of this variability was the variability in air temperature, which affected the variation in annual glacier runoff more than annual snowmelt runoff. Our simulation showed the presence of runoff-generating events outside the usual summer runoff measurement season (January-March). These warm events contribute significantly to the total annual runoff and are usually omitted in conventional observational studies that cover only the peak melt season.