Increasing glacier runoff in northwestern Greenland simulated from 1950 to 2023

Increased river runoff due to ice melting in Greenland contributes to sea-level rise, as well as flooding in coastal settlements, posing serious risks to local communities. To investigate fluctuations of glacier runoff in Greenland and its atmospheric drivers, long-term variations in runoff from Qaanaaq Glacier, northwestern Greenland, were reconstructed from 1950 to 2023 using a glacier energy–mass balance model and climate reanalysis dataset. Exceptionally large daily runoff (top 0.1 %) has only happened since 1990, indicative of an increasing frequency of major runoff events in recent decades. The largest (8.5 m³ s⁻¹ in 2023) and second largest (7.2 m³ s⁻¹ in 2001) runoffs resulted in the destruction of roads in the settlement of Qaanaaq, demonstrating the significant effects on the local community. The two largest runoffs have been attributed to intense rainfall contributing 79 % and 78 % of the daily runoffs, respectively, due to enhanced moisture and heat transport caused by an atmospheric river. Long-term annual glacier runoff increased at a rate of 0.7 Mt per decade, associated with the atmospheric warming with 0.24 °C per decade during 1950–2023. Rainfall exhibited 2.2-fold increase from 1951–1960 to 2010–2020. The variations in annual glacier runoff are controlled mainly by synoptic-scale atmospheric conditions represented by the Greenland Blocking Index (r=0.69). Composite analysis of the climate reanalysis dataset suggests particularly high sensitivity of air temperature in northern Greenland to anticyclonic conditions over Greenland, which lead to strengthened warm southerly winds. Accurate representation of such extreme conditions in climate models is crucial for predicting glacier runoff and flood occurrence in Greenland.

1 Introduction

The Greenland ice sheet has been losing mass over the past three decades, thereby contributing to a sea-level rise of 13.6±1.3 mm during 1992–2020 (Otosaka et al., 2023). More than half of the mass loss during 2000–2018 (55±5 %) was caused by negative surface mass balance (SMB) due to increased meltwater runoff (van den Broeke et al., 2016; Mouginot et al., 2019) associated with atmospheric warming and circulation changes over Greenland (Fettweis et al., 2013; van den Broeke et al., 2016). As a result of the mass loss, freshwater runoff from the Greenland ice sheet had been increasing at a rate of 16.9±1.8 km³ a⁻² during 1992–2010 (Bamber et al., 2012).

Local glaciers and ice caps located along the periphery of the Greenland ice sheet are susceptible to climate warming because they are located at low elevations (Mernild et al., 2011; Hanna et al., 2012). While these regions constitute a relatively small area of Greenland (4 %–7 %, depending on the connectivity to the ice sheet) (Rastner et al., 2012), their ice mass loss had a high contribution (11±2 %, when 4 % of area in Greenland was considered as local) to that across Greenland during 2003–2019 (Khan et al., 2022). Khan et al. (2022) also noted that ice mass loss across Greenland was heterogeneous, with the greatest changes occurring in north Greenland ($-\mathrm{26.1}\pm \mathrm{1.4}$ Gt a⁻¹), equivalent to 60 % of the mass loss from the local glaciers during 2018–2021. The substantial loss of ice in northern Greenland has been attributed to changes in atmospheric circulation. A recent poleward shift of the atmospheric ridge during summer has resulted in advection of a southerly warm air mass into northern Greenland (Fettweis et al., 2013; Tedesco et al., 2016). The circulation changes have also resulted in increased cloudiness and subsequent atmospheric warming due to decreased longwave heat loss, which enhanced rapid snowline retreat and an increase in meltwater runoff by increasing absorption of shortwave radiation by dark bare ice (Noël et al., 2019).

Increased runoff from local glaciers and the ice sheet poses serious problems for communities in Greenland. For example, the outlet stream from Qaanaaq Glacier in northwestern Greenland flooded in the summers of 2015 and 2016, which destroyed the roads in the settlement of Qaanaaq (Kondo et al., 2021). The increases in runoff were attributed to intense glacier ablation and rainfall. The outlet stream further flooded in the summers of 2022 and 2023. According to local residents in Qaanaaq, similar events occurred in 1997, 2001, and 2009, since the establishment of Qaanaaq Airport in 1991 (Kim Peterson, personal communication, 2018). Another flood has been reported in 2012 at Watson River, an outlet river from the ice sheet in central–western Greenland (Mikkelsen et al., 2016). The flood was attributed to an extreme melt event over the Greenland ice sheet on 11 July 2012, when surface melt covered up to 98 % of the entire ice sheet (Nghiem et al., 2012; Mikkelsen et al., 2016; van As et al., 2017). Damage on a bridge was also reported by Kobbefjord Research Station in southwestern Greenland near Nuuk after a rain event in 2022 (Kobbefjord Research Station, 2022). These events indicate an urgent need to understand the mechanisms of the rapid increase in glacier runoff to make reliable projections of future floods and mitigate their effects. These floods have been linked to extreme climatic conditions. However, no investigations have been conducted on flood occurrence in Greenland prior to 2012. Therefore, it is not clear whether the frequency and magnitude of large runoff events have been increasing owing to the atmospheric warming in the Arctic in recent decades.

Long-term variations in surface melt and runoff from the Greenland ice sheet have been related to large-scale climate variability. Previous studies have identified a significant negative correlation between the North Atlantic Oscillation (NAO) mode of atmospheric circulation and surface melting of the ice sheet (Chylek et al., 2004; Box et al., 2012; Hanna et al., 2013). The negative phase of the summer NAO index increases the prevalence of high pressure centered on Greenland, which enhances ice surface melting due to clear-sky conditions and a decrease in snowfall (Box et al., 2012; Fettweis et al., 2013). This condition also causes warm air advection to western Greenland and results in extreme ice sheet melting events (Hanna et al., 2014; Tedesco et al., 2016). Since the NAO is an indicator of the atmospheric circulation over the North Atlantic, the Greenland Blocking Index (GBI) is used as a more specific indicator of the conditions over Greenland (Hanna et al., 2013, 2016). The GBI is defined as the mean of the 500 hPa geopotential height over Greenland at 60–80° N and 20–80° W, and represents the prevalence of anticyclonic conditions over Greenland. Thus, the GBI is negatively correlated with the NAO. The correlation of the GBI with summer warmth in Greenland is stronger than that of the NAO, and thus is considered a good predictor of Greenland runoff (Hanna et al., 2013; Ahlstrøm et al., 2017). The occurrence of extreme Greenland blocking episodes has been increasing in recent decades (McLeod and Mote, 2016), and has possibly affected long-term variations in glacier runoff in the coastal regions of Greenland. However, our understanding of glacier runoff and its link to climate is limited due to the lack of river discharge measurements. Previous studies of glacier runoff have primarily focused on southern Greenland, where continuous river discharge observations are available (Rennermalm et al., 2012; Hasholt et al., 2013; Ahlstrøm et al., 2017; van As et al., 2017). Long-term records of glacier runoff using such reliable in situ measurements are scarce in northern Greenland.

In this study, we conducted glacier runoff modelling of Qaanaaq Glacier in northwestern Greenland. Long-term runoff variations were reconstructed for 1950–2023 using a glacier energy–mass balance model and a climate reanalysis dataset calibrated by field and satellite observations. We analyzed the mechanisms and atmospheric drivers of the rapid increase in runoff that potentially leads to flooding in Qaanaaq. We also investigated the effect of synoptic-scale atmospheric circulation on the annual climatic mass balance and runoff in northwestern Greenland.

Figure 1 (a) Location of the Qaanaaq region in Greenland. (b) Sentinel-2 image of Qaanaaq Ice Cap (21 August 2018), composed by SentinelFlow (v. 0.1.2). The black rectangle indicates the area shown in (c). (c) Sentinel-2 image of Qaanaaq Glacier (21 August 2018), showing the locations of the SIGMA-B weather station (diamond), runoff (cross) and surface mass balance (SMB) observations (circles). The black line shows the catchment of the outlet stream from Qaanaaq Glacier.

2 Study site

Qaanaaq Glacier is an outlet glacier from Qaanaaq Ice Cap in northwestern Greenland (77°28^′ N, 69°14^′ W; Fig. 1). The proglacial stream passes under a road connecting the settlement of Qaanaaq and Qaanaaq Airport, which is 2 km downstream from the glacier snout, before flowing into the ocean. Based on a digital elevation model (DEM) of ArcticDEM Mosaic, created by DEMs acquired within the period 2007–2022 (Porter et al., 2024), and a manually delineated glacier outline using Landsat 9 OLI-2 imagery (27 July 2023) show the catchment has a total area of 10.5 km², of which 75 % is covered with ice from elevations of 230 to 1110 m a.s.l. (Fig. S1 in the Supplement).

Various studies have been conducted on Qaanaaq Glacier since 2012. The SMB has been measured since 2012 at six stake sites along the glacier flowline at elevations of 243–968 m a.s.l. (Sugiyama et al., 2014; Tsutaki et al., 2017) (Fig. 1c). The glacier-wide mass balance in 2012–2016 was $-\mathrm{0.22}\pm \mathrm{0.3}$ $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ and the mean equilibrium line altitude was 910 m a.s.l. (Tsutaki et al., 2017). An automatic weather station (AWS) was installed at the SIGMA-B site at 944 m a.s.l. on 19 July 2012 and has since been operational (Aoki et al., 2014a). A quality-controlled dataset of observed variables at SIGMA-B (air temperature, relative humidity, air pressure, wind speed and direction, radiation components and glacier surface height) is publicly available (Nishimura et al., 2023a).

Darkening of glacial ice has been reported in the middle reaches of Qaanaaq Glacier, causing a spatially variable surface albedo and melt rate (Sugiyama et al., 2014; Takeuchi et al., 2014). The darkening of Qaanaaq Glacier is attributed to impurities that are mainly cryoconite granules (Takeuchi et al., 2014), which grow in association with mineral particles on the glacier (Uetake et al., 2016). The minerals are supplied to the glacier surface by wind-blown deposition, as well as from englacial dust (Nagatsuka et al., 2016; Matoba et al., 2020).

3 Methods

3.1 Glacier energy–mass balance model

We calculated the glacier runoff using a glacier energy–mass balance model, GLIMB (e.g., Fujita and Ageta, 2000; Fujita and Sakai, 2014), which has been previously used for calculating the mass balance and runoff of various mountain glaciers (e.g., Sakai and Fujita, 2017; Khalzan et al., 2022; Minowa et al., 2023) and the Greenland ice sheet (Fujita et al., 2021). The model calculates the climatic mass balance based on the snow/ice melt rate, snowfall rate, evaporation rate, and refreezing rate of melt and rainwater. The surface energy balance (Q; in W m⁻²) was calculated as follows:

$$\begin{array}{}\text{(1)}& Q=(\mathrm{1}-\mathit{\alpha })S+\mathit{\epsilon }L-\mathit{\epsilon }\mathit{\sigma }{\left(T_{\mathrm{s}}+\mathrm{273.15}\right)}^{\mathrm{4}}+H_{\mathrm{s}}+H_{\mathrm{l}}+G,\end{array}$$

where α is the albedo (dimensionless), S is the downward shortwave radiation (W m⁻²), L is the downward longwave radiation (W m⁻²), ε is the longwave emissivity (assumed to be 0.98; dimensionless), σ is the Stefan–Boltzmann constant ($\mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ W m⁻² K⁻⁴), and T_s is the surface temperature (°C) determined by an iterative calculation (Fujita and Ageta, 2000). The sensible (H_s) and latent (H_l) turbulent heat fluxes (W m⁻²) were estimated using a bulk aerodynamic method (Kondo, 1994). The conductive heat flux, (G; W m⁻²) was calculated from the temperature profile of the subsurface snow/ice computed within the model using the heat transfer equation (Fujita and Ageta, 2000). The energy flux toward the glacier surface was defined to be positive. The daily surface melt rate (M; m w.e. d⁻¹) was calculated when the surface energy balance was positive as follows:

$$\begin{array}{}\text{(2)}& M={\displaystyle \frac{t_{\mathrm{d}}Q}{{\mathit{\rho }}_{\mathrm{w}}{l}_{\mathrm{m}}}},\end{array}$$

where t_d represents the length of a day in seconds (86 400 s), and ρ_w and l_m represent the density of water (1000 kg m⁻³) and latent heat of fusion of ice (3.33×10⁵ J kg⁻¹), respectively. On ice-free area, M represents snowmelt when snow cover was present.

The daily runoff depth (d; m d⁻¹) was calculated as follows:

$$\begin{array}{}\text{(3)}& d=M+P_{\mathrm{r}}+max\left[E_{\mathrm{v}},\mathrm{0}\right]-R_{\mathrm{f}}-R_{\mathrm{T}},\end{array}$$

where P_r is the rainfall (m d⁻¹), E_v is the evaporation rate on the glacier surface that is taken into account only if it is positive (i.e., condensation; m w.e. d⁻¹), R_f is the refreezing rate (m w.e. d⁻¹) calculated based on the subsurface snow/ice temperature profile (Fujita and Ageta, 2000), and R_T is the retained liquid water within the snowpack given by the volumetric water content of the snow (Fujita and Ageta, 2000). The precipitation phase (snowfall or rainfall) was determined by assuming the probability of snowfall and rainfall depending on the air temperature (Fujita and Sakai, 2014). Daily river runoff (D; m³ s⁻¹) was calculated using the daily runoff depth from the glacier (d_g; m d⁻¹) and ice-free terrain (d_f; m d⁻¹) as follows:

$$\begin{array}{}\text{(4)}& D=\sum _z{\displaystyle \frac{A_{\mathrm{g},z}{d}_{\mathrm{g},z}+A_{\mathrm{f},z}{d}_{\mathrm{f},z}}{t_{\mathrm{d}}}},\end{array}$$

where A_g and A_f are the glacierized and ice-free areas within the catchment (m²) for the given elevation bin, z (Fig. S1), respectively.

The annual specific mass balance (b; m w.e. a⁻¹) was obtained from the sum of the daily mass balance over a mass balance year as follows:

$$\begin{array}{}\text{(5)}& b=\sum _t(P_{\mathrm{s},t}-M_t+E_{\mathrm{v},t}+R_{\mathrm{f},t}+R_{\mathrm{T},t}),\end{array}$$

where P_s is the snowfall rate (m w.e. d⁻¹) and t represents the days within the mass balance year. The glacier-wide annual mass balance (B; m w.e. d⁻¹) was obtained by integrating the annual mass balance at elevation bins as follows:

$$\begin{array}{}\text{(6)}& B={\displaystyle \frac{\sum _z{A}_{\mathrm{g},z}{b}_z}{\sum _z{A}_{\mathrm{g},z}}}.\end{array}$$

The model input variables were air temperature, relative humidity, downward shortwave and longwave radiation, wind speed, precipitation, and temperature lapse rate. All variables were taken from the reanalysis dataset at daily time intervals. The influence of the strong nonlinearity of the Stefan–Boltzmann law on the temporal discretization is negligible in daily calculations compared with that in hourly calculations. The calculations were conducted for every day and 50 m elevation bins. We calculated the climatic mass balance within a mass balance year defined from 25 August (in the preceding year) to 24 August, based on the annual maximum (15 June) and minimum (24 August) in the averaged cumulative mass balance. Winter and summer mass balances were computed from 25 August (in the preceding year) to 15 June and 16 June to 24 August, respectively.

3.2 Meteorological reanalysis data

We used the ERA5-Land reanalysis dataset (grid spacing of 0.1° × 0.1°) produced by the European Centre for Medium-Range Weather Forecasts (Muñoz Sabater, 2019) as the input for GLIMB. The meteorological variables were calibrated by the observed variables at SIGMA-B during 2012–2020 (Nishimura et al., 2023a). To best represent the observations at SIGMA-B, linear regression equations were applied to the ERA5-Land dataset for air temperature, relative humidity, downward shortwave radiation, and wind speed (Fig. S2). After calibration, the air temperature in each 50 m elevation bin was calculated using the temperature lapse rate obtained from the daily ERA5 pressure level data (Hersbach et al., 2023a). We calculated the lapse rate by differentiating the air temperatures at 975 and 875 hPa geopotential heights. The elevations were equivalent to 285±80 and 1117±90 m (mean ± standard deviation) over the study period, respectively, which cover 95 % of the elevation range of Qaanaaq Glacier. The lapse rate calculated with this method was consistent with the observed lapse rate at Qaanaaq Glacier (Kondo et al., 2021) within 13 % difference.

The downward longwave radiation (L) was calculated as an elevation-dependent variable. L from ERA5-Land (L_ERA) is expressed by the Stefan–Boltzmann law as

$$\begin{array}{}\text{(7)}& L_{\mathrm{ERA}}={\mathit{\epsilon }}_{\mathrm{eff}}\mathit{\sigma }(T_{\mathrm{ERA}}+\mathrm{273.15}{)}^{\mathrm{4}},\end{array}$$

where ε_eff is the effective atmospheric emissivity (van den Broeke et al., 2004; Ruckstuhl et al., 2007) and T_ERA is the air temperature from ERA5-Land. After calibrating the air temperature with SIGMA-B data (T_calib) using a linear regression equation (Fig. S2a), the corresponding downward longwave radiation from ERA5-Land (L_calib) is expressed as

$$\begin{array}{}\text{(8)}& L_{\mathrm{calib}}=L_{\mathrm{ERA}}{\left({\displaystyle \frac{T_{\mathrm{calib}}+\mathrm{273.15}}{T_{\mathrm{ERA}}+\mathrm{273.15}}}\right)}^{\mathrm{4}}.\end{array}$$

We compared L_calib with the observed L at SIGMA-B and corrected it with a linear regression equation in the same manner as the other meteorological variables (Fig. S2e). After correction, L_calib was used to calculate L in each elevation bin (L_z) with the air temperature (T_z) being

$$\begin{array}{}\text{(9)}& L_z=L_{\mathrm{calib}}{\left({\displaystyle \frac{T_z+\mathrm{273.15}}{T_{\mathrm{calib}}+\mathrm{273.15}}}\right)}^{\mathrm{4}}.\end{array}$$

A previous study noted an underestimation bias in precipitation in climate reanalysis data due to relatively coarse grid spacing that is difficult to resolve given the terrain effect on precipitation (Box et al., 2023). Precipitation is not observed at the SIGMA-B site, and thus we calibrated daily precipitation for ERA5-Land by applying a multiplier (r_p) that was tuned to minimize the discrepancy between the modelled and measured glacier surface height at SIGMA-B from 2013 to 2019 by iterative calculations. Surface height was calculated by assuming fresh snow density and snow densification parameters that were optimized using a Monte Carlo simulation in a previous study to best reproduce the observed snow surface height and albedo in the accumulation zone of the northwestern region of the Greenland ice sheet (Fujita et al., 2021). The performance was evaluated with the mean error (ME) between the modelled and observed surface height.

ERA5 reanalysis data were also used to analyze the synoptic-scale climatic conditions (grid spacing of 0.25° × 0.25°). We used the daily precipitation, surface air temperature, and vertically integrated water vapor transport from the ERA5 single-level data (Hersbach et al., 2023b), as well as the 500 hPa geopotential height and air temperature and wind speed at 850 hPa geopotential height from the ERA5 pressure-level data (Hersbach et al., 2023a). To better resolve the spatial patterns of precipitation with finer scales, we also analyzed the Copernicus Arctic Regional Reanalysis (CARRA) (Schyberg et al., 2020), which is available at a spatial resolution of 2.5 km, focusing on periods of particularly high runoff. The CARRA reanalysis was also used as a higher-resolution reference to evaluate the performance of precipitation variability represented in ERA5-Land.

3.3 Glacier surface albedo

Snow surface albedo was calculated with the scheme of Kondo and Xu (1997), which assumes an exponential reduction of snow albedo with time after a fresh snowfall. We utilized parameters optimized by a Monte Carlo simulation in the previous study in northwestern Greenland (Fujita et al., 2021). We modified the fresh snow albedo to 0.7 instead of the optimized value of 0.897, based on the snow surface albedo observed by SIGMA-B (Nishimura et al., 2023a).

The bare ice albedo of Qaanaaq Glacier is spatially variable, depending on the altitude due to ice darkening. To consider the spatial variability, we measured the bare ice albedo using satellite imagery at the six stake sites and the SIGMA-B site (Fig. 1c). We analyzed Landsat 8 OLI surface reflectance imagery using Google Earth Engine. The total shortwave albedo was calculated from the previously proposed conversion formulae using narrowband reflectance measurements (Liang, 2001). We analyzed 146 images acquired during 2013–2023 with a cloud coverage of <10 %. Pixels covered with clouds or cloud shadows were eliminated using the quality assessment band provided in the Landsat 8 dataset. Radiometric saturation was also detected and eliminated using the saturation quality assessment band. The filtered images were manually inspected to ensure the remaining pixels were not covered with clouds or cloud shadows. We used the mean albedo values after the disappearance of the snow cover to use representative value for bare ice albedo. Albedo measured by the satellite at the SIGMA-B site was compared with that observed by the SIGMA-B AWS (Fig. S3). The comparison showed agreement between the two observed values, with a root-mean-square error (RMSE) of 0.07 and ME of 0.06, which indicate there is an overestimation bias in the satellite-derived surface albedo. The bias was eliminated by applying a linear regression correction (Fig. S3g), and this albedo was subsequently used for GLIMB. After the bias correction, RMSE between the two observations was 0.02.

3.4 Validation data

In situ glaciological and hydrological data at Qaanaaq Glacier were used to validate GLIMB. We used the annual SMB data obtained during 2012–2016, which cover four mass balance years (Tsutaki et al., 2017). Runoff observations of the outlet stream from Qaanaaq Glacier were used to validate the glacier runoff calculations. The runoff was measured in the summer periods during 2017–2019 at 1.4–2.0 km from the glacier terminus (Mankoff et al., 2020; Kondo et al., 2021).

3.5 Return period of the annual maximum runoff

We analyzed the return period of the calculated annual maximum daily runoff based on probability exceedance. The annual maximum values were fitted to the two-parameter Gumbel distribution (Gumbel, 1941). The parameters of the probability distribution function were estimated with the probability-weighted moments method (Greenwood et al., 1979). The goodness of the fit was evaluated with the probability plot correlation coefficient (PPCC) (Vogel, 1986). The result confirmed that PPCC >0.97, indicative of significance at the 95 % level for the given number of samples.

4 Results

4.1 Model performance

As a result of tuning, the precipitation multiplier, r_p, was determined to be 1.41. The multiplier was consistent with that from a previous study in Greenland, which was reported to be r_p=1.4, based on rainfall measurements on the ice sheet compared with ERA5 data (Box et al., 2023). The calculated surface height at the SIGMA-B site captures the observed variations with a RMSE and ME of 0.19 and 0.12 m, respectively (Fig. S4a). The surface albedo also exhibited a good agreement with the observations (RMSE and ME of 0.07 and 0.01, respectively) (Fig. S4b). The time-series of the calibrated daily precipitation from ERA5-Land was further compared with that from the CARRA reanalysis dataset (Fig. S5). Although extreme precipitation events (>60 mm d⁻¹) were not reproduced in ERA5-Land, overall performance exhibited good agreement with RMSE of 2.1 mm. The comparison indicates that the precipitation variability is well captured in the calibrated ERA5-Land precipitation, even when compared with the higher spatial resolution dataset.

Figure 2 Observed versus calculated surface mass balance in 2012–2016. The colors indicate stake sites (Fig. 1c). The dashed line indicates a linear regression through the data. The coefficient of determination of the regression (R²) and RMSE between the calculated and observed surface mass balance are also shown.

The calculated climatic mass balance was compared with the observations made on the glacier during 2012–2016 (Fig. 2). The calculated annual SMB at each site ranged from −2.37 to 0.09 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$, which agrees with the range of observed SMBs (−2.15 to 0.16 m w.e. a⁻¹). The spatiotemporal variations in SMB were well reproduced with the RMSE between the calculated and observed SMB of 0.23 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$, indicating a sufficiently good performance for reconstructing the glacier variability. Substantially negative SMBs of below −1.5 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ were also well reproduced with a RMSE of 0.25 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ (13 % of the mean), indicating good agreement even during periods of intense ablation.

Figure 3 (a–c) Calculated (black) and observed (red) daily runoff with uncertainties (shading) for Qaanaaq Glacier in the summers of 2017–2019. (d) Observed versus calculated daily runoff. The dashed line is a linear regression of the data. The coefficient of determination of the regression (R²) and RMSE between the calculated and observed runoff are also shown.

The calculated glacier runoff agreed well with the observed proglacial discharge (Fig. 3). The Nash–Sutcliffe coefficient (Nash and Sutcliffe, 1970) of the calculated glacier runoff versus the observed discharge was 0.80, which is well above the threshold values proposed for accurate hydrological models (Moriashi et al., 2007). The RMSE between the calculated and observed runoff was 0.50 m³ s⁻¹. The model tends to overestimate smaller runoffs (Fig. 3d), resulting in an overall ME of 0.22 m³ s⁻¹. The calculated runoff agreed well with the observed seasonal variability in runoff, including the initiation and cessation of summer runoff, which was particularly well captured by the hydrograph in 2019 (Fig. 3c). The rapid increase in daily runoff on 29 July 2017 was also reproduced by the model (Fig. 3a), suggesting the calculations can be used to investigate periods of high runoff. Although we did not apply a routing model to calculate the discharge time-series, the results were reliable. This is consistent with a previous hydrological study conducted on Qaanaaq Glacier (Kondo et al., 2021) as well as Arctic glaciers (Hodgkins, 1997), which suggested a relatively short routing delay caused by the supraglacial stream rather than draining into the glacier bed. This delay is crucial in reproducing the hourly time-series of the hydrograph (Kondo et al., 2021), but is not significant for the daily runoff modelling conducted in this study.

Figure 4 Daily runoff averaged for the periods 1960–1979 (black), 1980–1999 (blue), and 2000–2019 (red). The shading represents one standard deviation for the corresponding period.

4.2 Daily and seasonal runoff

GLIMB with calibrated variables was used to calculate the daily runoff of Qaanaaq Glacier during 1950–2023. The daily values in summer (June, July, and August) averaged 0.9 m³ s⁻¹. The runoff typically increases from mid-May to mid-July, and decreases from early August, with the mean number of runoff days being 108±11 (with one standard deviation) (Fig. 4). Hydrographs averaged over 20-year intervals showed an increase in runoff from 1960 to 2020. Mean daily runoff in the summer period increased from 0.6 m³ s⁻¹ in 1960–1979 to 0.8 and 1.2 m³ s⁻¹ in 1980–1999 and 2000–2019, respectively, indicating a 2.1-fold increase over 60 years. The mean number of runoff days increased by 16 % from 1960–1979 (100±7 d) to 2000–2019 (116±11 d).

Figure 5 (a) Calculated daily runoff (blue line). Inverted triangles indicate flood events at Qaanaaq. The horizontal dashed and dotted lines are the 99th and 99.9th percentiles of non-zero daily runoff, respectively. (b) Histogram of the calculated daily runoff exceeding the 99th (black) and 99.9th (white) percentiles for the decadal periods. (c) Runoff depth on the 10 largest runoff days showing the percentage contributions of glacier melt (blue) and rainfall (orange). The colors of the circles and white numbers represent the dominant component and its percentage contribution, respectively. The numbers above the circles indicate the return period (years) of the daily runoff.

A long-term time-series of daily runoff shows multi-decadal variations in the occurrence of high runoff (Fig. 5a). A relatively high annual maximum runoff was obtained for 1950–1962, which was followed by consistently lower runoff during 1963–1979, with daily runoff never exceeding the 99th percentile (3.9 m³ s⁻¹) (Fig. 5a). Thereafter, runoff increased gradually in recent decades, resulting in a 1.7-fold greater mean annual maximum since 1990 as compared with the prior period. Exceptionally large runoff events exceeding the 99.9th percentile (5.5 m³ s⁻¹) have all occurred since 1990 (Fig. 5a and b). Furthermore, the frequency of such events increased even more after 2000, with 91 % and 88 % of days exceeding the 99th and 99.9th percentiles, respectively (Fig. 5b). These runoff peaks were triggered by both rainfall and melt, indicating a range of factors contributing to each event (Fig. 5c).

The largest daily runoff over the past 74 years was calculated for 22 August 2023 (8.5 m³ s⁻¹), and was >3σ above the mean annual maximum (Fig. 5a and c). The value had a return period of 120 years, exceeding the 100-year runoff event threshold (8.3 m³ s⁻¹) (Fig. S6). The second largest runoff was calculated for 2 August 2001 (7.2 m³ s⁻¹) and was >2σ above the mean annual maximum. These two high-runoff days coincided with the destruction of roads in Qaanaaq Village, demonstrating clear evidence of substantial societal effects. Heavy rainfall occurred on these two largest runoff days (43 and 37 mm d⁻¹ over the catchment in 2023 and 2001, respectively), which accounted for 79 % and 78 % of the runoff on these 2 d, respectively (Fig. 5c). The third largest daily runoff was equivalent to 6.4 m³ s⁻¹ on 12 July 2012 during a period of extreme surface melting over the Greenland ice sheet (Nghiem et al., 2012; Fausto et al., 2016). While this event is known to have resulted in significant surface melting and runoff from the ice sheet, approximately half of the daily runoff (52 %) was attributed to rainwater on Qaanaaq Glacier (Fig. 5c). On the same day, rainfall was observed by an in situ field campaign at the SIGMA-A site (1490 m a.s.l.) located on the Greenland ice sheet, 65 km inland from Qaanaaq Glacier (Aoki et al., 2014b). Rainfall was estimated as 100 mm during 10–14 July 2012 based on the ERA-Interim reanalysis data corrected by in situ bucket measurements (Niwano et al., 2015). The fourth, fifth, and sixth largest daily runoffs were calculated to have occurred on 13 July 2010, 1 August 1990, and 19 July 2014, with values of 5.9, 5.8, and 5.7 m³ s⁻¹, respectively. The bulk of these daily runoff events was derived from glacier melt (95 %, 100 %, and 72 %, respectively).

Figure 6 (a) Annual runoff and the contribution of snow/ice melt (blue) and rain (orange). (b) Calculated glacier-wide summer (red), winter (blue), and annual (dot) mass balances. (c) Summer mean air temperature at the SIGMA-B site taken from ERA5-Land. (d) Precipitation (black), snowfall (light blue), and rain fraction (bars) averaged over the catchment. The dashed lines in (c) and (d) are linear regressions of the data with p<0.01. The trend slope of the regression is also shown.

4.3 Annual runoff and glacier mass balance

Long-term annual runoff variations are predominantly determined by snow/ice melt, with the contribution from rainfall being consistently low (mean of 12 %) (Fig. 6a). The interannual variability of the annual runoff exhibits a standard deviation of 3.2 Mt a⁻¹, equivalent to 44 % of the average (7.1 Mt a⁻¹). The largest and second largest annual runoff were calculated for 2011 and 2019 with values of 15.0 and 14.0 Mt a⁻¹, respectively, coinciding with high summer mean air temperatures (1.4 °C at SIGMA-B for both years). The annual variations exhibited patterns similar to those of the annual maximum daily runoff, characterized by relatively high values in the 1950s, lower values in the 1960s and 1970s, and a subsequent increase in the following decades (Fig. 6a). An increasing trend of annual runoff during 1965–2023 (0.12 Mt a⁻²; p<0.001) was associated with warming air temperatures in the same period (0.33 °C per decade) (Fig. 6c). The annual variations were also characterized by substantially suppressed runoff, including the five lowest values in 1964, 1992, 1976, 1996, 1953 (Fig. 6a).

The variations in annual glacier runoff are well explained by the summer glacier mass balance, as indicated by their high correlation (r=0.96; p<0.001) (Figs. 6a, b, and S7). For example, relatively high runoff during 1957–1960 coincided with negative summer mass balance. This period corresponded with the years when the summer mean temperature remained near average (Fig. 6c), while precipitation was substantially reduced (Fig. 6d), suggesting snowfall was important for summer glacier ablation and runoff by modulating glacier surface conditions during the ablation season. An increasingly negative summer mass balance during 1965–2023 (−0.10 m w.e. a⁻¹ per decade) explains the increase in annual runoff during this period. Suppressed annual runoff, including the five lowest values (Fig. 6a), corresponded with a glacier mass budget surplus (Fig. 6b) due to the low summer air temperatures and/or greater snowfall (Fig. 6c and d).

The average annual climatic mass balance was −0.21 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ with a standard deviation of 0.33 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ for the period 1951–2023. The most negative mass balance was −0.99 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in 2019 and the most positive was 0.36 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ in 1996. The annual mass balance showed similar variations as the summer mass balance, having a significantly higher correlation (r=0.99; p<0.001) than with the winter mass balance (r=0.33; p<0.01) (Fig. S7). The annual mass balance was also correlated with the summer mean air temperature ($r=-\mathrm{0.77}$; p<0.001) and snowfall within the mass balance year (r=0.56; p<0.001). After the mass loss in 1957–1960 and mostly positive mass balance in 1961–1965, the annual mass balance decreased at a rate of −0.09 m w.e. a⁻¹ per decade by 2023 due to the increasingly negative summer mass balance. As a result of the mostly negative mass balance throughout the study period, the cumulative mass loss has been 15.7 m w.e. when averaged over the entire glacier.

The contribution of rainfall to the annual runoff was consistently low over the entire period. Overall, total precipitation has been increasing over the past 74 years at a rate of 20.8 mm per decade (Fig. 6d). This rate exceeds that of snowfall (12.4 mm per decade), and thus the rain fraction of the total precipitation has been increasing at 1.2 % per decade. As a result, rainfall exhibited a 2.2-fold increase from 1951–1960 to 2010–2020. No significant correlation was found between the annual runoff and total precipitation ($r=-\mathrm{0.07}$; p>0.05), but a weak correlation was observed with rainfall (r=0.31; p<0.01) and snowfall ($r=-\mathrm{0.25}$; p<0.05) (Fig. S7). The negative correlation with snowfall suggests suppressed/enhanced glacier runoff coincided with greater/less snowfall.

5 Discussion

5.1 Synoptic controls on high daily runoff

Our results demonstrate that the frequency and magnitude of large daily runoff have been increasing significantly since 2000 (Fig. 5a and b). These runoff events were driven by both rainfall and glacier melt, with the top three daily runoff events being caused primarily by rainfall (Fig. 5c). The two largest runoff events resulted in the destruction of roads in Qaanaaq, which highlights the importance of understanding the mechanisms of such events. We note that not all the flood events observed in Qaanaaq were related to a large runoff (Fig. 5a). For example, a moderate runoff (3.4 m³ s⁻¹) was calculated on the day of a flood in 1997 observed by the local residents, which was close to the mean annual maximum of 3.3 m³ s⁻¹. This unclear relationship is due to the complexity of flood occurrence, which is not only related to runoff magnitude, but also infrastructure vulnerability. We also note that our model may not adequately reproduce thick firn formed prior to 1950, and ice layers within the firn, which can impede downward water infiltration and thereby facilitate runoff from the accumulation area (Machguth et al., 2016; Mikkelsen et al., 2016). These complex supraglacial hydrological processes may partly explain the unclear relationship between flood occurrence and the runoff magnitudes.

Figure 7 Composites of the ERA5 reanalysis data during the three largest runoff days caused by (a–c) rain and (e–g) glacier melt. Composites of the CARRA reanalysis data for the same days were also shown (d, h). (a, e) Air temperature at 850 hPa geopotential (color) and 500 hPa geopotential height (gray contours). (b, f) Vertically integrated water vapor transport (color) and wind vectors at 850 hPa geopotential (arrows). (c, g) Daily precipitation from ERA5. (d, h) Daily precipitation from CARRA. The red circles indicate the location of Qaanaaq.

To identify the mechanisms of the particularly high peak runoff, we analyzed the atmospheric conditions during the top-ranked runoff days. We extracted top three days for which more than half of the runoff originated from rainwater (22 August 2023, 2 August 2001, and 12 July 2012) (Fig. 5c). For the extracted days, we created composites for the daily atmospheric reanalysis of the 500 hPa geopotential height, air temperature and wind speed at 850 hPa geopotential, vertically integrated water vapor transport, and precipitation (Fig. 7a–d). The results showed that an atmospheric ridge centered over Greenland strengthened southerly warm air advection onto the northern ice sheet, which delivered considerable water vapor to northwestern Greenland along Baffin Bay (Fig. 7a and b). In northwestern Greenland, the wind was directed toward the inland part of the ice sheet, as indicated by the pressure contour that crosses the Greenlandic coastline from west to east (Fig. 7a and b), where the location corresponds to the intense precipitation at Qaanaaq. The rainfall reached 23 mm d⁻¹ in the composite from ERA5 (Fig. 7c), which was twice the average annual maximum intensity. The precipitation patterns were further resolved by CARRA reanalysis dataset with a finer spatial resolution, showing intensive precipitation exceeding 60 mm d⁻¹ along the coast of the northwestern part of the ice sheet, including 28 mm of daily precipitation at Qaanaaq (Fig. 7d). From the ERA5 dataset, the air temperature at 850 hPa geopotential height (∼1400 m a.s.l.) was >2.9 °C at the nearest grid to Qaanaaq for all the analyzed days. The results suggest warm and humid air transported from the south caused heavy rainfall as a result of orographic uplift, given the geometry of the ice sheet and local ice caps located along the coast.

Exceptionally strong moisture transport to northwestern Greenland was observed on the largest runoff day of 22 August 2023. Concentrated poleward water vapor transport exceeding 600 kg m⁻¹ s⁻¹ was identified for >1500 km along Baffin Bay, due to the strong wind (>20 m s⁻¹ at an 850 hPa geopotential), enhanced by low- and high-pressure systems located over the Canadian Arctic Archipelago and Greenland, respectively (Fig. S8a and b). The moisture transport resulted in the substantial rainfall and runoff at Qaanaaq (Fig. S8c). This feature is consistent with atmospheric river formation, which is defined as a narrow band of strong water vapor transport exceeding 150 kg m⁻¹ s⁻¹ (Mattingly et al., 2018). Strong moisture transport exceeding 150 kg m⁻¹ s⁻¹ was also observed during the other rain-induced peak runoffs on 2 August 2001 and 12 July 2012, suggesting atmospheric river contributed to rainfall events in Qaanaaq. An atmospheric river resulted in light rain (3.4 mm) at Summit station on the Greenland ice sheet in the summer of 2021 (Box et al., 2022), as well as extreme melting events over the Greenland ice sheet, such as in 2012 (Neff et al., 2014; Fausto et al., 2016). Together with the recent increase in the activity of atmospheric rivers over Greenland (Mattingly et al., 2016, 2018), our results highlight the importance of such extreme climatic conditions for flood occurrence in northern Greenland.

We next investigated the conditions required for the peak runoff caused by glacier melt. Weather conditions were analyzed for the top three melt-induced peak runoff events (13 July 2010, 1 August 1990, and 19 July 2014) (Fig. 5c). The surface energy balance calculations indicate a substantial increase in the sensible heat flux during these melt events, reaching an average of 137 W m⁻². The contribution was twice greater than that from net shortwave radiation (61 W m⁻²) on these days, which differs from the typical conditions during the annual maximum runoff, when the shortwave radiation (72 W m⁻²) was 1.7 times greater than the sensible heat (43 W m⁻²). Net longwave radiation was also larger during the melt-induced peak runoff (3 W m⁻²) than the period averaged for the annual maximum (−10 W m⁻²). From the ERA5-Land reanalysis data taken from the nearest grid to SIGMA-B, strong winds (mean of 8.9 m s⁻¹) and a high mean air temperature (4.4 °C) were obtained for the glacier, suggesting warm and strong winds with cloudy conditions resulted in the elevated sensible heat flux. The weather and heat conditions are similar to the situation during the flood event in Qaanaaq during the summer of 2015, which caused intense melt of Qaanaaq Glacier (Kondo et al., 2021).

To investigate the atmospheric conditions that led to intense glacier ablation, composites of the three largest melt-induced runoff days were created for the atmospheric variables in the same way as for the rain-induced runoff (Fig. 7e–h). Similar to the periods of rain-induced peak runoff, relatively high pressure was located over Greenland. Additionally, low-pressure trough was located over the Canadian Arctic Archipelago (Fig. 7e). This pressure pattern resulted in a warm southerly air flow toward Qaanaaq along Baffin Bay, resembling the wind pattern during the rain-induced runoff events (Fig. 7f). As a result, the composite means of wind speed and air temperature at 850 hPa geopotential height reached 15.8 m s⁻¹ and 2.9 °C, respectively, at Qaanaaq. The southerly contains considerable amounts of water vapor (Fig. 7f), but not much precipitation was observed at Qaanaaq (mean of 3.1 and 8.8 mm d⁻¹ from ERA5 and CARRA, respectively) (Fig. 7g and h). This is because the wind vector was directed toward Ellesmere Island in the Canadian Arctic Archipelago, where daily precipitation reached 21 and 36 mm d⁻¹ from ERA5 and CARRA reanalysis, respectively. Therefore, dominant influences from the low-pressure system in the Canadian Arctic resulted in strong and warm winds, but less precipitation at Qaanaaq.

The foregoing analysis shows there were similar atmospheric conditions for both the rain- and melt-induced peak runoff events, characterized by the transport of warm and humid air mass from the south along the western coast of Greenland. Previous studies showed that the maximum latitude of the high-pressure ridge over Greenland shifted poleward during recent summer periods, which increased ice sheet runoff in more northern regions of Greenland (Fettweis et al., 2013; Tedesco et al., 2016). Our results demonstrate that such atmospheric conditions are important not only for the mass balance of the ice sheet, but for the occurrence of peak runoff potentially leading to flooding on the coast of northwestern Greenland. Furthermore, the northwestern region has been identified as the area experiencing the greatest increase in rainfall in Greenland, given the four-fold increase over the period 1980–2019 (Niwano et al., 2021). Although no trend has been reported for extreme rainfall intensity in the northern region, care should be taken for monitoring rainfall associated with the changes in the atmospheric circulation that favors warm air transport into northern Greenland.

5.2 Climatic controls on glacier mass loss

Our mass balance estimates for Qaanaaq Glacier were compared with the previously reported glacier-wide elevation changes measured by satellite. Elevation change rates on Qaanaaq Glacier, as obtained by laser altimetry (Khan et al., 2022), were $-\mathrm{0.49}\pm \mathrm{0.06}$ and $-\mathrm{0.39}\pm \mathrm{0.06}$ $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ (assuming an ice density of 890 kg m⁻³) during 2003–2009 and 2008–2019, respectively, which are consistent with our modelling results during these periods (−0.46 and −0.37 m w.e. a⁻¹). Another study reported a value of $-\mathrm{0.44}\pm \mathrm{0.09}$ $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ for 2010–2020 based on DEMs (Hugonnet et al., 2021), which is also consistent with our calculated value of −0.43 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$. The good agreement of the results with the geodetically derived mass loss confirms that our results quantitatively capture the key variations in Qaanaaq Glacier, although we assumed a fixed altitude–area distribution. We assessed the effect of the assumption of the fixed altitude–area distribution by calculating a runoff time-series using a DEM and the glacier extent derived from aerial photographs in 1985 (AeroDEM; Korsgaard et al., 2016). The comparison showed a 7 % increase in the mean annual runoff, which has little effect on the interannual variations of runoff as well as the increasing trend over recent decades (Fig. S9). The occurrence of large peak runoff events also remained largely unaffected by the changed hypsometry, given that the daily runoff prior to 1990 never exceeded the 99.9th percentile, which is consistent with the original calculation.

Figure 8 Summer mean air temperature (T) versus glacier-wide mass balance (MB). Symbol colors indicate the snowfall anomaly in the corresponding mass balance year. The dashed line indicates a linear regression through the data.

Glacial meltwater accounts for 88 % of the total river runoff, and thus glacier fluctuations are a key control on the runoff variations. Glacier mass balance is strongly affected by summer climatic conditions and glacier ablation. The high correlation between the annual climatic mass balance and summer air temperature or snowfall suggests these are significant controls on the runoff. Multiple regression analysis further highlighted their importance, showing these two variables explain 88 % of the variability of the annual mass balance, with a RMSE of 0.12 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ (Fig. S10; Table S1 in the Supplement). The relationships between climatic mass balance, summer air temperature, and snowfall are shown in Fig. 8. The linear regression revealed the sensitivity of the mass balance to the summer air temperature, which predicted 0.29 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ of ice mass loss with a 1 °C warming in air temperature. The plot also indicates the mass balance being positive/negative during greater/lesser snowfall as shown in the deviation from the linear regression line (Fig. 8). The amount of snowfall has a stronger effect on the mass balance by affecting glacier melt, rather than increasing the accumulation. Snowfall in the mass balance year modulates the timing of snowline retreat, which influences absorption of shortwave radiation on the dark ice surface. For example, substantially less snowfall (33 % less than average) was recorded in 1960, the year of the most negative mass balance prior to 2000, which resulted in exceptionally rapid snowline retreat. Based on the model calculations, the snowline rose to 1075 m a.s.l. by 2 July, which was 400 m above the average on this day. As a result of the rapid disappearance of snow cover, absorption of shortwave radiation in the summer of 1960 (115 W m⁻²) was 1.5-fold greater than the average (79 W m⁻²). Although increasing summer air temperature is the primary driver of the decreasing mass balance, the results suggest the snowfall amount is a secondary control on glacier fluctuations.

Figure 9 (a) Modelled annual runoff (blue line) and the Greenland Blocking Index (GBI) in summer (orange line). Circles indicate the five years with the highest (red) and lowest (blue) summer GBI. Composites of ERA5 reanalysis data for 500 hPa geopotential height (white lines) and its anomaly from the temporal trend (color) are shown for the five years with the (b) highest and (d) lowest summer GBI. Composites of ERA5 reanalysis data for the surface air temperature anomaly from the temporal trend are shown for the five years with the (c) highest and (e) lowest summer GBI. The red circles indicate the location of Qaanaaq.

To identify the synoptic-scale climatic conditions affecting annual glacier runoff, we examined climate indices. We found the most significant correlation was with the GBI during summer (June, July, and August), with a coefficient r=0.69 (p<0.001) (Fig. 9a). The correlation is stronger than that with the summer NAO ($r=-\mathrm{0.60}$; p<0.001), suggesting a more direct link with the atmospheric conditions over Greenland as represented by the GBI. The atmospheric conditions during high/low-GBI periods were analyzed by creating composites for 500 hPa geopotential height and surface air temperature anomalies for the five years with the highest/lowest summer GBI (Fig. 9b–e; Table S2). High summer GBI resulted in similar atmospheric conditions as observed during large daily runoff events, which were characterized by an atmospheric ridge over Greenland that resulted in strengthened southerly winds into northern Greenland (Fig. 9b). This is consistent with previously reported atmospheric conditions during the negative NAO phase (Fettweis et al., 2013; McLeod and Mote, 2016). The composite for the surface air temperature anomaly showed pronounced warming in northern Greenland, as well as the Canadian Arctic Archipelago, with a greater magnitude than in the southern region (Fig. 9c). Conversely, the surface temperature in the north cooled during the five years with the lowest summer GBI (Fig. 9e), due to a southward shift of the cold air mass in the polar region (Fig. 9d). Therefore, the significant correlation between annual runoff of Qaanaaq Glacier and the GBI is due to the high sensitivity of air temperatures in the north to anticyclonic conditions over Greenland. The pronounced increase in annual runoff since 2000 coincided with the years with the top-ranked summer GBI (Fig. 9a), indicative of the effects of large-scale atmospheric variability. The results are consistent with recent large-scale circulation changes reported over Greenland (Noël et al., 2019). Qaanaaq Glacier is located in the region with the most pronounced increase in summer cloud cover and subsequent atmospheric warming. Our results highlight the importance of monitoring glacier hydrology in northwestern Greenland, given its direct response to recent changes in atmospheric circulation.

The high sensitivity of the air temperature in northern Greenland to blocking events highlights the importance of synoptic-scale atmospheric conditions over Greenland for accurately predicting future runoff variability. Our results demonstrate that a similar pressure pattern corresponded with increases in runoff on both long and short timescales. However, such atmospheric circulation patterns are not captured by the Earth system models used in CMIP5 and CMIP6, as none of these models reproduced the recent trend of increasing GBI (Hanna et al., 2018; Maddison et al., 2024). Inadequate representation of the atmospheric conditions leads to an underestimation of the glacier runoff and flood occurrence in northern Greenland.

6 Conclusions

To investigate long-term glacier runoff variations in northwestern Greenland, we conducted runoff modelling using a glacier energy–mass balance model and climate reanalysis dataset for 1950–2023. Our results identified an increase in the occurrence of large daily runoff events in recent decades, with all the exceptionally large runoff events (>99.9th percentile) having occurred since 1990. The largest (8.5 m³ s⁻¹ in 2023) and second largest (7.2 m³ s⁻¹ in 2001) daily runoff events coincided with flood events observed in Qaanaaq, which demonstrated substantial impacts on the society of the Greenlandic coast. These rapid increases in runoff were due to intense rainfall (contributions of 79 % and 78 %, respectively) caused by an atmospheric river directed toward Qaanaaq, which in turn was generated by an anticyclone centered over Greenland. The pressure pattern strengthened the advection of warm and humid southerly air to the north. High daily runoff was also caused by intense glacier melting triggered by strong and warm winds, due to southerly winds strengthened by a low-pressure system located over the Canadian Arctic Archipelago.

Long-term annual glacier runoff increased at a rate of 0.7 Mt per decade during 1950–2023, associated with the atmospheric warming with 0.24 °C per decade and 2.2-fold increase of rainfall from 1951–1960 to 2010–2020. The variations in annual glacier runoff were significantly affected by synoptic-scale atmospheric conditions as represented by the GBI (r=0.69; p<0.001). Analysis of the atmospheric conditions in periods of high summer GBI showed that stronger warming occurred in northern Greenland and the Canadian Arctic Archipelago than in the south, suggesting high sensitivity of glacier melt in the northern regions. The warming during high-GBI periods was caused by strengthened warm air advection due to Greenland blocking, highlighting the significant effect of anticyclonic conditions over Greenland on both short- and long-term glacier runoff in northwestern Greenland. Given the inadequate representation of Greenland blocking by climate models, these models are likely to substantially underestimate future glacier runoff and flood occurrence.