Flood events caused by discharge from Qaanaaq Glacier, northwestern Greenland

Abstract

As a result of climate warming, glacial meltwater discharge has been increasing in Greenland. During the summers of 2015 and 2016, there were rapid increases in discharge from Qaanaaq Glacier in northwestern Greenland. These discharges resulted in floods that destroyed the road linking the settlement of Qaanaaq to Qaanaaq Airport. Field measurements were performed and a numerical model of glacier runoff was developed to quantify these discharges. The high discharge associated with the 2015 flood, estimated at 9.1 m³ s⁻¹ (hourly mean), resulted from intensive glacier melting due to warm air temperature and strong winds, while the high discharge associated with the 2016 flood resulted from heavy rainfall (90 mm d⁻¹) that led to a peak discharge estimated at 19.9 m³ s⁻¹. The developed model, when used to investigate future glacier runoff under warming conditions, revealed a nonlinear increase in glacial melt with increasing temperature. Additionally, the model forecasted a threefold increase in total summer discharge, owing to a 4 °C rise in temperature. Thus, this study quantified the impact of a changing climate on glacier runoff, which gives insight into future risks of flood hazards along the coast of Greenland.

Introduction

From 1972 to 2018, the Greenland ice sheet experienced a decrease in mass, contributing to a sea-level rise of 13.7 ± 1.1 mm (Mouginot and others, 2019). Approximately 60% of the ice-sheet mass loss since 2000 can be attributed to the increasingly negative surface mass balance owing to the warming trend of the atmosphere over Greenland (Enderlin and others, 2014; van den Broeke and others, 2016). Accordingly, since 1992, fresh water runoff from the Greenland ice sheet has increased by 16.9 ± 1.8 km³ a⁻² (Bamber and others, 2012; Ahlstrøm and others, 2017). Glaciers and ice caps located along the coast of Greenland are also affected by atmospheric warming. Reportedly, for the period 2003–08, the mass loss of the peripheral glaciers and ice caps in Greenland was −28 ± 11 Gt a⁻¹ (Bolch and others, 2013). Given that glaciers and ice caps are located at relatively low elevations, they are sensitive to recent atmospheric warming trends in coastal regions (Mernild and others, 2011; Hanna and others, 2012).

An increase in meltwater runoff from the Greenland ice sheet, glaciers and ice caps poses a serious threat to the local communities. In the summer of 2012, there was an increase in the discharge from the Watson River, which is an outlet of the Kangerlussuaq sector of the ice sheet in central western Greenland, resulting in a washout and the destruction of the Watson River Bridge on 11 July 2012 (Mikkelsen and others, 2016; van As and others, 2017). Additionally, on 12 July 2012, the Greenland ice sheet experienced an extreme melt event in which surface melt covered up to 98.6% of the ice sheet (Nghiem and others, 2012). It has been suggested that this melt event was responsible for the disaster that was observed at Watson River (Mikkelsen and others, 2016; van As and others, 2017).

Floods caused by glacial discharges have also occurred in northern parts of Greenland. On 21 July 2015 and 2 August 2016, in the Qaanaaq region in northwestern Greenland, the outlet stream of Qaanaaq Glacier flooded. Qaanaaq Glacier is an outlet glacier of Qaanaaq Ice Cap and has, over the past decade, been studied by Japanese scientists (Sugiyama and others, 2014; Machguth and others, 2016; Tsutaki and others, 2017; Takeuchi and others, 2018). The flooding of this glacial stream resulted in the destruction of a road linking the settlement of Qaanaaq to Qaanaaq Airport (Fig. 1). Since the establishment of Qaanaaq Airport in 1991, three similar incidents have been recorded – in 1997, 2001 and 2009 (personal communication from K. Petersen, 2018). The floods in 2015 and 2016 destroyed the road which had been reinforced after the previous events, thus we presume that the flood magnitude was greater than those in 1997–2009. Based on existing literature, hydrological monitoring has never been performed in this outlet stream and the discharge parameters associated with the flooding events that occurred are still unknown. Therefore, to monitor and/or forecast similar or even more severe events in the future, investigations are necessary.

Fig. 1. Aerial photograph of the destruction of the road linking the settlement of Qaanaaq to Qaanaaq Airport by the outlet stream from Qaanaaq Glacier. (Photograph was taken on 3 August 2016, from a helicopter).

Glacier hydrology is linked to flood prevention, hydroelectric power generation, water resource assessment and glacier mass balance. Considering the various studies on temperate glaciers in alpine settings (e.g. Moore, 1993; Hock and Noetzli, 1997; Verbunt and others, 2003; Zhang and others, 2007; Fujita and Sakai, 2014; Bravo and others, 2017), relatively few have focused on glaciers and ice caps in Greenland (Bøggild and others, 1999; Liston and Mernild, 2012). Given the remoteness of Greenland and its harsh environment, river discharge measurements in this area are scarce (Rennermalm and others, 2012; Hasholt and others, 2013; van As and others, 2017). Previous glacier runoff modeling in Greenland primarily focused on the southwestern part, owing to the limited river discharge time series that are available. Sustained efforts (Hasholt and others, 2013; van As and others, 2017) provide highly accurate long-term discharge data for Watson River and can be used for a better understanding of surface mass balance and ice-sheet hydrology (Mikkelsen and others, 2016; Mernild and others, 2018). Compared to southern Greenland, very few studies have been conducted on glacial runoffs in northern Greenland (Mernild and others, 2008).

Glaciers and ice caps located in Greenland and other Arctic are characterized by ice temperatures that are lower than the pressure melting point (Blatter and Hutter, 1991; Curffey and Paterson, 2010). Given that the ice is impermeable, glacial meltwater runs on the surface of the glacier rather than draining into the glacier bed (Hodgkins, 1997; Irvine-Fynn and others, 2011). Therefore, surface processes play a central role in the seasonal and short-term patterns of glacier discharges in Greenland. Particularly, the retention and refreezing of meltwater within snow layers control the magnitude and timing of the runoff (Hodgkins, 1997). Despite the importance of glacier surface hydrology and glacier discharge in rapidly changing coastal environments in Greenland, studies of the link between the two are scarce.

Motivated by the flood events that occurred in 2015 and 2016, in this study we investigate the discharge of the outlet stream from Qaanaaq Glacier. Based on field observations between 2016 and 2019, we developed a numerical model to quantify the discharges during the 2015 and 2016 floods. The model was also used to forecast glacial discharges under warming atmospheric conditions.

Study site

Qaanaaq Ice Cap is located in northwestern Greenland (77°28′ N, 69°14′ W) and covers an area of 289 km² at an elevation between 30 and 1110 m a.s.l. (above sea level) (Sugiyama and others, 2014), as shown in Figures 2a and b. Qaanaaq Glacier is one of the outlet glaciers of this ice cap that flows southwestward and terminates on land at 190 m a.s.l. (Fig. 2c). The mean equilibrium line altitude of the glacier over the 2012–16 period was 910 m a.s.l. (Tsutaki and others, 2017). The glacial meltwater drains into an outlet stream that intersects a road before emptying into the ocean (Fig. 2c). This road connects Qaanaaq Village and Qaanaaq Airport; thus, it is an important part of the infrastructure for the ~600 residents of this village and its visitors. Based on a digital elevation model provided by the Greenland Ice Mapping Project (Howat and others, 2014), the outlet stream has a catchment area of 11 km², of which 78% is covered with the glacier (Fig. 3). The glacier area below 700 m a.s.l. is relatively small (i.e. 17% of the total glacier area), because it narrows into valley geometry.

Fig. 2. (a) Landsat 8 OLI image (24 July 2017) showing northwestern Greenland, including the Qaanaaq region. The inset shows the location of the Qaanaaq region in Greenland; the red box indicates the area shown in (b). (b) Sentinel-2 image of Qaanaaq Ice Cap (21 August 2018); the red box indicates the area shown in (c). (c) Sentinel-2 image of the study site (21 August 2018), showing the locations of the weather stations (♦), discharge measurement sites (+) and surface melt and snow survey sites (•).

Fig. 3. Hypsometry of the catchment area of the outlet stream from Qaanaaq Glacier partitioned into glacier-covered (black) and ice-free (white) areas.

Since 2012, Japanese research groups have been conducting in situ studies on Qaanaaq Glacier. Sugiyama and others (2014) carried out glaciological observations and reported melt rates, ice velocity and ice thickness along the flowline. Tsutaki and others (2017) reported near-surface ice temperature, as well as surface mass balance and ice velocity measured from 2012 to 2016. Aoki and others (2014) installed an automatic weather station on the ice cap in 2012 (SIGMA-B in Fig. 2c) and it has been operational since then. Qaanaaq glacial microbes and their influence on glacier melt have also been studied (Uetake and others, 2016; Takeuchi and others, 2018). Darkening resulting from microbial activity reduces ice surface albedo and enhances surface melt on the ice cap (Sugiyama and others, 2014). To quantify the surface elevation change of ice caps in this region, satellite-based research was conducted (Saito and others, 2016) and the results indicated that the mass-loss rate of Qaanaaq Ice Cap was 1.8 ± 0.1 m a⁻¹ between 2007 and 2009.

During the summer months (June, July and August) of the 2016–19 period, glaciological observations were performed on Qaanaaq Glacier at six survey sites covering the glacier snout up until the accumulation area, from an elevation of 243–968 m a.s.l. (Fig. 2c) and the discharge of the outlet stream was measured during the summer months of 2017–19.

Methods

Field measurements

Glacier surface melt

Snow and ice melt measurements on Qaanaaq Glacier were conducted by utilizing six aluminum poles installed at 243, 427, 584, 739, 839 and 968 m a.s.l. (Fig. 2c). Since 2012, mass balance and ice velocities are also available at these sites (Sugiyama and others, 2014; Tsutaki and others, 2017). The poles were drilled into the snow or ice surface using a mechanical ice drill (Kovacs Enterprise) attached to an electric drill driver (Makita HR262DRDX). Surface melt and snow accumulation were inferred by frequent recordings of the surface height at temporal intervals of 1–14 days from 12 June to 2 August 2017, 6 July to 19 August 2018 and 29 June to 22 August 2019. Conversion of surface height to melt/accumulation in water equivalent was done by assuming ice density of 917 kg m⁻³ and utilizing measured snow densities that were conducted with a 100 cm³ box-type stainless steel density cutter and an electronic scale.

Proglacial discharge

Outlet stream measurements were performed 1.4 and 2.0 km away from the glacier terminus in 2017 and 2018–19, respectively (Fig. 2c). A pressure sensor (HOBO U20-001-04) was installed in the stream to record water pressure between 21 July and 3 August 2017, 5 July and 18 August 2018, and 27 June and 26 August 2019. The sensor was mounted on an aluminum pole and fixed on the stream floor to measure water pressure at 5 min intervals. After correcting the data for atmospheric pressure variations using data recorded by another sensor exposed to air, the water pressure values were converted to water stage. The accuracy of the pressure sensor was translated to a water depth of ±3 mm.

During the 2017, 2018 and 2019 field observations, the discharge of the stream was measured 6, 15 and 12 times, respectively. For the discharge calculations, water depth and current were measured every 0.5 m across the stream using an electromagnetic current meter (YOKOGAWA ES-7603) at a depth of 60% from the surface, which represents the mean current of a vertical water column. Each measurement was performed for 1 min at each location and the mean current was integrated over the vertical cross section of the stream. The uncertainty of each discharge measurement was 8% based on the method described by Sauer and Meyer (1992).

The water stage observations were used to compute the discharge time series, assuming a power function between the water stage and discharge (Fig. 4). Relationships between water stage and discharge were determined independently for 2017 and 2018/19, given the different locations of the water pressure sensors for the different periods. The discharge uncertainty obtained from the water stage measurements was 9% and 22% in 2017 and 2018/19, respectively (Sauer and Meyer, 1992; Rennermalm and others, 2012).

Fig. 4. Relationships between water stage and discharge of the outlet stream from Qaanaaq Glacier in (a) 2017 and (b) 2018/19. Solid curves, coefficient of determination (R ²) and RMSE represent the results of the power regression analysis of the data. Shaded areas represent the estimated uncertainty ranges.

Meteorological observations

Meteorological observations were recorded from 1 June 2015 to 31 August 2019 by an automatic weather station at SIGMA-B (944 m a.s.l.) (Fig. 2c) (Aoki and others, 2014). Air temperature, downward shortwave and longwave radiation, relative humidity, wind speed, atmospheric pressure and ice/snow surface height were recorded with hourly resolution. In this same period, we recorded air temperature and liquid precipitation at the Qaanaaq Village at 68 m a.s.l. (Fig. 2c). The air temperature sensor (T&D RTR-502L) and the rain gauge (Climatec CEM-TBRG) operated from 1 June 2015 to 27 August 2019, recording their respective parameters at the hourly resolution, with accuracies of ±0.3 °C and ±0.2 mm, respectively.

Numerical modeling

A glacier runoff routing model was developed, based on the 2016–19 field observations. The model consisted of a surface energy balance and linear reservoir schemes and was used to compute snow/ice melt and glacier runoff, respectively. Its purpose was to reconstruct the 2015 and 2016 discharges, to better understand the underlying processes that led to the flood events and also to forecast glacial discharges under the atmospheric warming trend in Greenland.

Surface melt model

The surface energy-balance model used in this study was based on equations that are commonly used for snow and ice melt studies (Kondo, 1994; Oerlemans, 2001). These formulations have been previously employed in the studies of glacier mass balance (e.g. Fujita and Ageta, 2000; Klok and Oerlemans, 2002; Van As, 2011; van Pelt and others, 2012) and glacier runoff (e.g. Liston and Mernild, 2012; Fujita and Sakai, 2014). The energy flux available for melting (H _M) was expressed as the residual of the net shortwave radiation (S _net), net longwave radiation (L _net), sensible heat flux (H _S), latent heat flux (H _L), conductive heat flux due to rainfall (R) and conductive heat flux from subsurface (G), i.e.,

(1)$$\matrix{ {H_{\rm M} = S_{{\rm net}} + L_{{\rm net}} + H_{\rm S} + H_{\rm L} + R + G.} \cr } $$

The energy flux into the glacier was considered to be positive. In this study, given that the contribution of the conductive heat flux from the subsurface and rainfall was relatively small during the melt season, it was ignored (Kondo and Yamazaki, 1990; Van As, 2011). The sensible and latent heat fluxes were estimated using the bulk method and were expressed respectively as follows:

(2)$$\matrix{ {H_{\rm S} = c_{\rm p}\rho C_{\rm S}U\lpar {T_{\rm a}-T_{\rm s}} \rpar \comma \;} \cr } $$

(3)$$\matrix{ {H_{\rm L} = l\rho C_{\rm S}U\lsqb {h_{\rm r}q\lpar T_{\rm a}\rpar -q\lpar {T_{\rm s}} \rpar } \rsqb \comma \;} \cr } $$

where c _p = 1006 J kg⁻¹ K⁻¹ represents the specific heat capacity of air under constant pressure, C _S = 0.002 represents the bulk coefficient for snow and ice surfaces proposed based on field data (Kondo and Yamazawa, 1986; Kondo, 1994) and previously used for glacier runoff modeling (e.g. Fujita and Sakai, 2014), l represents the latent heat of evaporation of water (2.5 × 10⁶ J kg⁻¹) when the ice surface is melting (T _s = 0°C) or the latent heat of sublimation of ice (2.83 × 10⁶ J kg⁻¹) during freezing (T _s < 0°C). U, T _s, T _a and h _r represent wind speed, surface temperature, air temperature and relative humidity, respectively. In Eqn (1), as long as H _M > 0, the surface temperature, T _s was equal to 0°C. Otherwise, Eqn (1) was solved for T _s by substituting H _M = 0. The density of air (ρ) and the saturated specific humidity (q) were given by the following equations (Kondo, 1994):

(4)$$\matrix{ {\rho = 1.293\displaystyle{{273.15} \over {273.15 + T_a}}\left({\displaystyle{\,p \over {\,p_0}}} \right)\left({1-0.378\displaystyle{{h_{\rm r}e} \over p}} \right)\comma \;} \cr } $$

(5)$$\matrix{ {q = \displaystyle{{0.622e} \over {\,p-0.378e}}\comma \;} \cr } $$

where p and p₀ = 1013.2 hPa are atmospheric pressure and its standard value. The saturated vapor water pressure (e) was calculated from air temperature values using Tetens formulation (Tetens, 1930). The net shortwave radiation (S _net) and the net longwave radiation (L _net) were expressed as:

(6)$$\matrix{ {S_{{\rm net}} = \lpar {1-\alpha } \rpar S_{{\rm down}}\comma \;\;} \cr } $$

(7)$$\matrix{ {L_{{\rm net}} = L_{{\rm down}}-\;\varepsilon \sigma {\lpar {T_{\rm s} + 273.15} \rpar }^4\comma \;} \cr } $$

where S _down and L _down represent the downward shortwave and longwave radiation, respectively; σ, which represents the Stefan–Boltzmann constant, was equal to 5.67 × 10⁻⁸ W m⁻² K⁻⁴; and ɛ, which represents longwave emissivity, was assumed to be 1. For the bare ice surface (α _i), the surface albedo, α was optimized and the optimization was performed at each pole measurement site within the range α _i = 0.2–0.8. To identify the best fit between the computed melt and the observations made in the 2016–19 period, RMSE was minimized; α _i = 0.6 was assumed for the upper two poles because of the unavailability of ice surface measurements. The influence of this assumption on melt calculations was insignificant, given that these two poles were mostly snow-covered. The ice surface albedo determined at each pole was linearly interpolated to other elevations, and the albedo of fresh snow (α _s(0)) was given by the following equations (Fujita and Sakai, 2014):

(8)$$\matrix{ {\alpha _{\rm s}\lpar 0 \rpar = 0.9\;\lsqb {T_{\rm a} \lt -1.0^\circ {\rm C}} \rsqb \comma \;\;} \cr } $$

(9)$$\matrix{ {\alpha _{\rm s}\lpar 0 \rpar = \displaystyle{{\lpar {\alpha_{\rm s}\lpar \infty \rpar -0.9} \rpar \lpar {T_{\rm a} + 1.0} \rpar } \over {4.0}} + 0.9\;\lsqb {-1.0^\circ {\rm C} \le T_{\rm a} \lt 2.0^\circ {\rm C}} \rsqb .\;} \cr } $$

The lower limit of snow albedo α _s(∞) was assumed to be 0.53 (Klok and Oerlemans, 2002). Furthermore, snow albedo, which decreases as a function of time (Kondo and Xu, 1997; Oerlemans and Knap, 1998), was expressed as:

(10)$$\matrix{ {\alpha _{\rm s}\lpar t \rpar = \lpar {\alpha_{\rm s}\lpar {t-\Delta t} \rpar -\alpha_{\rm f}} \rpar {\rm e}^{-\lpar {1/t^\ast } \rpar } + \alpha _{\rm f}\comma \;} \cr } $$

where t represents the time after snowfall, Δt = 1 h represents the model time step and t* was considered to be 360 h, following Bougamont and others (2005).

The model input variables were air temperature, relative humidity, wind speed, atmospheric pressure, downward shortwave and longwave radiation measured at SIGMA-B at 1 h intervals. The model calculation was performed for the period of 1 June to 31 August in 2015–19. Initial snow depth on 1 June was taken after the glacier observations made each year. The surface energy balance was computed for every hour and for every 100 m elevation band. The air temperature at each elevation band was linearly interpolated from those observed at SIGMA-B and Qaanaaq Village, while wind speed and relative humidity were assumed to be uniform given the relatively small study area (5.5 km from SIGMA-B to the glacier snout).

The surface melt rate, M was calculated in the water equivalent thickness using Eqn (11).

(11)$$\matrix{ {M = \displaystyle{{H_{\rm M}} \over {\rho _{\rm w}l_m}}\comma \;\;} \cr } $$

where ρ _w and l _m represent the density of water (1000 kg m⁻³) and the latent heat of fusion of ice (3.34 × 10⁵ J kg⁻¹), respectively.

Linear reservoir model

The simulation of the outlet stream discharge was based on a linear reservoir model (e.g. Hock and Jansson, 2005), taking into account the meltwater storage of snowpacks and the refreezing within firn layers or on ice surfaces (Fig. 5). The sum of the glacier surface melt and the liquid precipitation was considered as the input of the reservoir (I).

(12)$$\matrix{ {I = M + P_{\rm r}.} \cr } $$

Fig. 5. Schematic representation of the routing model used in this study. S, storage of a linear reservoir; Q, runoff from a linear reservoir; S ₀, retention capacity; and k _ice and k _snow, storage coefficients for ice- and snow-covered surfaces, respectively.

Solid (P _s) or liquid precipitation (P _r) was distinguished based on the air temperature at each elevation band i.e.,

(13)$$\matrix{ {P_{\rm r} = P\;\lsqb {T_{\rm a} \ge 2^\circ {\rm C}} \rsqb \comma \;} \cr } $$

(14)$$\matrix{ {P_{\rm s} = P\;\lsqb {T_{\rm a} \lt 2^\circ {\rm C}} \rsqb \comma \;} \cr } $$

where P represents the liquid precipitation measured in Qaanaaq Village. The rain-snow threshold temperature was chosen based on previous studies in Greenland (Oerlemans and Vugts, 1993; Bougamont and others, 2005). Outflow from the reservoir (Q) resulted when the volume of the stored water (S) exceeded the retention capacity (S ₀), which represents the amount of water stored or refrozen on the glacier surface until the end of the melt season. The discharge from the reservoir was assumed to be proportional to the volume of the stored water that exceeded the retention capacity, with a storage coefficient k.

(15)$$\matrix{ {Q\lpar t \rpar = 0\;\lsqb S \lt S_0\rsqb \comma \;} \cr } $$

(16)$$\matrix{ {Q\lpar t \rpar = \displaystyle{{S-S_0} \over k}\;\lsqb {S \ge S_0} \rsqb .} \cr } $$

The stored water volume was given by a volume balance equation, i.e.,

(17)$$\matrix{ {\displaystyle{{dS} \over {dt}} = I\lpar t \rpar -Q\lpar t \rpar .} \cr } $$

Eqns (15)–(17) were solved forward in time with a time step of 1 h for every 100 m elevation band. The river discharge (Q _sum), which represents the sum of the discharge from an elevation band Q_i was given by:

(18)$$\matrix{ {Q_{{\rm sum}} = \mathop \sum \limits_{i = 1}^n Q_i.} \cr } $$

Based on a previously reported measurement from Greenland, S ₀ in the ablation area was considered equal to 0.1 m (Cooper and others, 2018). In the accumulation area, it was assumed to increase from 0.3 m in 900‒1000 m a.s.l. to 0.5 m in ⩾1100 m a.s.l. as a firn layer thickens up-glacier, based on the superimposed ice thickness observed on Qaanaaq Glacier (Tsutaki and others, 2017). The determination of storage coefficient k was based on surface properties, i.e., k _land, k _ice and k _snow for ice-free and ice and snow-covered surfaces, respectively. A smooth transition between k _ice and k _snow was achieved by introducing a snow depth function, d in the ice–snow transition zone.

(19)$$\matrix{ {k = k_{{\rm ice}}\;\lsqb {d = 0\;{\rm m\;w}.{\rm e}.} \rsqb \comma \;\;} \cr } $$

(20)$$\matrix{ {k = k_{{\rm ice}} + \lpar {k_{{\rm snow}}-k_{{\rm ice}}} \rpar \times \displaystyle{d \over {0.3}}\;\lsqb {0 \lt d \lt 0.3\;{\rm m\;w}.{\rm e}.} \rsqb \comma \;} \cr } $$

(21)$$\matrix{ {k = k_{{\rm snow}}\;\lsqb {d \ge 0.3\;{\rm m\;w}.{\rm e}.} \rsqb .\;} \cr } $$

To identify the best fit between the modeled and observed discharges in the 2017–19 period, the parameters k _land, k _ice and k _snow were optimized and the optimization was performed within parameter ranges of 5 < k _land  < 50, 1 < k _ice < 15 and 10 < k _snow < 70 (h), respectively based on previous studies of similar-sized glaciers (Hannah and Gurnell, 2001; Hock and Jansson, 2005). For parameter tuning and for the evaluation of the performance of the model, the Nash–Sutcliffe coefficient (NSC) (Nash and Sutcliffe, 1970) was used.

Results

Field measurements

Meteorological conditions

The mean temperatures in summer (June, July and August) at SIGMA-B in 2015, 2016 and 2019 – 1.7, 1.5 and 2.0°C, respectively – indicated remarkably warm atmospheric conditions compared to those observed in 2017 and 2018, which were −0.4 and −1.0°C, respectively (Table 1). As shown in Figures 6a and b, persistent high air temperatures were observed prior to the 2015 and 2016 floods. From 26 June to 5 August 2015, in particular, the daily mean air temperature was consistently above the melting point of ice; over the 5-day period prior to the flood on 21 July, the daily mean temperature was above 5°C and a maximum of 9.5°C was recorded on 18 July, suggesting intensive glacier melt. Similarly, warm conditions were observed before the 2016 flood event. The daily mean temperature was above 0°C from 3 July to 4 August 2016 and a maximum of 7.7°C was observed on 7 July. In addition to the warm conditions, extraordinarily heavy rainfall (90 mm d⁻¹) occurred on 2 August 2016 (Table 1, Fig. 6b), the date of the flood event. The meteorological data showed an even greater amount of precipitation in 2015 (102 mm d⁻¹ on 5 August) as shown in Figure 6a. The mean wind speed in summer ranged from 3.2 to 4.0 m s⁻¹ with the maximum daily wind speed of 14.6 m s⁻¹ observed on 30 July 2017 (Table 1).

Fig. 6. Air temperature at the SIGMA-B station (dotted), liquid precipitation at Qaanaaq Village (blue), stream discharge (solid) with uncertainty (shaded gray) and the melt rate of Qaanaaq Glacier at 839 m a.s.l. (light blue) observed during the 2015–19 period. The vertical purple dashed lines in (a) and (b) represent the 2015 and 2016 flood dates.

Table 1. Maximum and mean daily temperature for the summer months (June, July and August) during the 2015–19 period at the SIGMA-B station on Qaanaaq Glacier. The maximum daily and summer rainfall levels, maximum and mean daily wind speeds are also listed. The dates in brackets represent the event dates.

Melt rate

Significant variations in glacier surface melt were observed for the different years. This is consistent with the abovementioned temperature trends (Figs 6b–e). In 2016, the mean melt rate from 3 July to 30 July, which was one of the most intensive melt periods over the 2016–19 study period, at 839 m a.s.l. was 17.4 mm w.e. d⁻¹ (Fig. 6b). The mean air temperature during this period was 4.2°C. The following year (2017), the melt rate was relatively low (< 4 mm w.e. d⁻¹), except from 21 July to 2 August, when a mean melt rate of 16.6 mm w.e. d⁻¹ was observed under a mean temperature of 2.2°C (Fig. 6c). In 2018, the melt rate was generally lower than that observed in previous years, i.e., <6 mm w.e. d⁻¹ at 0.5°C from 23 July to 19 August as shown in Figure 6d and in 2019, a relatively higher melt rate (15.9 mm w.e. d⁻¹) was observed from 29 June to 22 August (Fig. 6e) and the mean temperature during this period was 4.0°C.

Discharge

During the first week of the observational period in 2017, there was a gradual increase in the discharge, with daily fluctuations (Fig. 6c). However, there was a rapid increase in the discharge on 28 July 2017, which peaked at 9.8 m³ s⁻¹ and was the highest value recorded over the observation period (2017–19) as shown in Figure 6c. In 2018, the discharge was relatively small until the end of July, which was followed by an increase to a peak discharge of 3.4 m³ s⁻¹ on 6 August 2018 (Fig. 6d). In 2019, the discharge was consistently high, and a peak value of 5.5 m³ s⁻¹ was observed on 25 July (Fig. 6e). Additionally, in 2019, the amplitude of the daily fluctuations was significantly greater than those observed in previous years.

It should be noted that our measurements were utilized for validation of daily discharge from all Greenlandic drainage basins estimated from output of two regional climate models (Mankoff and others, 2020). The comparison of the data from 2017 to 2019 indicated that one of the models (MAR) (e.g. Fettweis and others, 2017) overestimated the total discharge by 29%, whereas the other (RACMO) (e.g. Noël and others, 2019) underestimated by 56%. The model performance at the Qaanaaq Glacier basin was relatively poor (NSC = 0.30 for MAR and −0.01 for RACMO). We speculate that the resolution of the climate models was too coarse (7.5 km for MAR and 5.5 km for RACMO) to accurately reproduce discharge from a small basin. We also suspect that a lack of long-term meteorological observations in the Qaanaaq region affected the accuracy of the models.

Numerical modeling

Surface melt

The ice surface albedo was optimized for the observations in the 2016–19 period. For the poles at 243–739 m a.s.l., α _i = 0.37–0.75 was considered, resulting in an RMSE = 0.09 m w.e. and a correlation coefficient R ² = 0.90 between the simulated and observed melt rates (Fig. 7). The optimized ice albedo was particularly low at 739 m a.s.l. (α _i = 0.37), where the ice was covered with dark material (Sugiyama and others, 2014). The tuned model was used to reproduce the cumulative melt at SIGMA-B during summer (June, July and August) in the 2015–19 period (Fig. 8). The total summer melt obtained in 2015, 2016 and 2019 was greater than that obtained in 2017 and 2018. The total melt in 2019 was 0.98 m w.e., which was 3.9-fold higher than that obtained in 2018 (0.25 m w.e.). In 2015 and 2016, the cumulative melts prior to the flood events were 0.63 and 0.62 m w.e., respectively.

Fig. 7. Scatter plot for observed and simulated melt at all observation sites in 2016 (crosses), 2017 (triangles), 2018 (circles) and 2019 (boxes). The dashed line represents the linear regression of the data. Coefficient of determination (R ²) and RMSE between simulated and observed melt rates are given.

Fig. 8. Cumulative melt at 950 m a.s.l. simulated using the surface energy-balance model for 2015 (red), 2016 (blue), 2017 (green), 2018 (orange), 2019 (brown); the mean of the 2015–19 data (dashed). Vertical lines represent the dates of the flood events in 2015 (orange) and 2016 (cyan).

The simulated daily melt was compared to the discharge observations shown in Figure 9. The modeled melt rates were substantially greater than the observed discharges for most of the observational period, which was likely due to meltwater retention and refreezing on the glacier. These results indicated the necessity of a linear reservoir model to accurately compute the discharge time series.

Fig. 9. Daily discharge measured at the outlet stream from Qaanaaq Glacier (black), daily melt rate simulated using the surface energy balance model (snowmelt, light blue; and ice melt, dark blue) and liquid precipitation observed in Qaanaaq Village (orange).

Stream discharge

Agreement between modeled and observed discharge was greatly improved by introducing the linear reservoir model (Fig. 10). After the optimization of the storage constants (k _ice = 4 h, k _snow = 55 h and k _land = 10 h), the NSC between the runoff routing model and the discharge measurements within the 2017–19 period was 0.71. This value was found to be well above the threshold values proposed for sufficiently accurate hydrological models (Moriasi and others, 2007). The performance of the model was particularly good in 2017 and 2018 (Figs 10a and b). The RMSEs between the model and the observed measurements for these 2 years were 0.6 and 0.3 m³ s⁻¹, respectively, and in 2019, a higher RMSE was observed (0.9 m³ s⁻¹). The modeled discharge in 2019 underestimated and overestimated the measurements in early July and mid-August, respectively, as shown in Figures 10c and d.

Fig. 10. (a)–(c) Simulated (red) and observed (blue) hourly discharge with uncertainty (shaded) within the 2017–19 period. (d) Scatter plot of simulated and observed hourly discharge in 2017 (triangles), 2018 (circles) and 2019 (boxes). The dashed line represents the linear regression of the data. Coefficient of determination (R ²) and RMSE between simulated and observed discharge are given.

The model was applied to the summer months of 2015 and 2016 (Fig. 11). A significant increase in discharge was simulated in the first week of July 2015, owing to the influence of rain (Fig. 11a); for the rest of the month, it was maintained at a high level, i.e. the daily peak discharge never dropped below 2.5 m³ s⁻¹ until the end of the month. The total discharge obtained in July (8.9 × 10⁶ m³) accounted for 73% of the total discharge in the summer of 2015. The daily peak discharge increased progressively, reaching 9.1 m³ s⁻¹ on 21 July, the day of the flood event. According to the model, 75% of the discharge on 21 July originated from ice melt, whereas the rest originated from snowmelt. Additionally, a higher peak discharge of 13.3 m³ s⁻¹ was obtained on 5 August 2015. This event could be attributed to rainwater, which accounted for 62% of the discharge on that day.

Fig. 11. Simulated discharge of the outlet stream from Qaanaaq Glacier in 2015 and 2016. The discharge is partitioned into rain (orange), snowmelt (light blue) and ice melt (dark blue). The shaded gray areas represent the dates of the flood events.

A large amount of discharge was simulated for 2016 as well (Fig. 11b). In July 2016, the total discharge was 7.1 × 10⁶ m³, which corresponded to 81% of that obtained in July 2015 and 58% of that obtained during the summer months of 2016. On 2 August 2016, owing to rainfall on the day of the flood event, the discharge increased sharply to 19.9 m³ s⁻¹, which was the greatest simulated discharge over the study period (2015–19). The contribution of rainwater to the discharge on that day was 73%, while rainwater from 2 to 4 August accounted for 12.5% of the summer discharge obtained in 2016.

Discussion

Drivers of the floods in 2015 and 2016

The modeled discharge in 2015 and 2016 revealed details of the 21 July 2015 and 2 August 2016 flood events (Fig. 11). On the day of the 2015 flood, the mean melt rate over the glacier was 51 mm w.e. d⁻¹. This episodic melt event was due to a considerable heat flux elevated by northeasterly strong winds under warm conditions, i.e., a daily mean wind speed of 7.4 m s⁻¹ and a temperature of 5.4°C at SIGMA-B (944 m a.s.l.). Similar meteorological conditions were associated with the rapid increase in discharge observed on 28 July 2017 and 6 August 2018 (Figs 6c and d). The strong winds (11.5 m s⁻¹ in 2017 and 17.8 m s⁻¹ in 2018) enhanced the glacier surface melt to 53 mm w.e. d⁻¹ and 25 mm w.e. d⁻¹ in 2017 and 2018, respectively, leading to peak discharges. Thus, both the simulated and the observed meteorology revealed that wind speed, as well as temperature, play an important role in melt-induced discharge events.

The discharge during the 2016 flood resulted primarily from rainwater (90 mm d⁻¹ on 2 August 2016) as shown in Figures 6b and 11b. Previous studies have reported an increase in precipitation in the Arctic over the past decades (Kattsov and Walsh, 2000; Walls and others, 2020) and owing to warming atmospheric conditions, as well as retreating sea ice that favor a greater moisture transport from the ocean, future increases in rainfall in the Arctic have been projected (Bintanja and Selten, 2014; Bintanja and Andry, 2017). Thus, in the future, a higher risk of flooding owing to an increase in heavy rainfall in the Arctic should be expected.

In addition to the obvious influence of glacier melt and rainfall on the 2015 and 2016 floods, the effect of snowpack disappearance on glaciers has also been suggested as another contributing factor to the flood events. Given that Qaanaaq Glacier has an impermeable cold ice surface, snowpack is the most important meltwater storage mechanism. Therefore, the disappearance of the snowpack facilitates immediate runoff into the proglacial stream. Under the influence of summer melt and winter snow accumulation, the seasonal pattern of the snowline migration on Qaanaaq Glacier is significantly different from year to year. In 2015 and 2016, persistently high air temperatures in June and July resulted in intensive glacier surface melt (Figs 6a, b and 8). Thus, the snowline was above SIGMA-B (944 m a.s.l.) on 31 July 2015 and 2016, and the largely exposed bare ice had small water retention capacity, which led to immediate runoff of water generated by the intensive melt event in 2015 and the rain event in 2016. Glacier surface conditions also affect short-term discharge variations. The diurnal patterns in the hydrograph that were attributed to the exposed bare ice were most pronounced in 2019 (Figs 10a–c). On 24 July 2019, the snowline was at 839 m a.s.l., where the surface was still covered by 0.85 and 1.25 m thick snow layers in the same week in 2017 and 2018, respectively. The extent of bare ice area in Greenland has been increasing (Ryan and others, 2019); hence, significant impact on glacial streams owing to discharge regimes should be expected.

Discharge under warming climate conditions

The summer air temperature in Greenland is projected to rise by 2°C in 2050 and 4°C in 2100, based on estimations made using Representative Concentration Pathway 8.5 (RCP 8.5), which assumes a high greenhouse gas emission rate (Vizcaíno and others, 2014). This warming enhances glacier surface melt, thereby increasing the risk of future floods. To quantify the future proglacial discharges of Qaanaaq Glacier, the glacier runoff routing model was run assuming the abovementioned projected temperature rise.

As an example of the sensitivity experiments, we show discharge simulated by imposing 2°C air temperature increases on the meteorological conditions from 2017 (Fig. 12). The results showed that a 2°C increase in air temperature resulted in a 144% increase in the total summer discharge (from 5.6 × 10⁶ to 13.8 × 10⁶ m³). The magnitude of the increase was more pronounced during the early melt season and the discharge in June was projected to be 4.1-fold greater than that based on 2017 observations. Owing to this increase in temperature, the peak discharges on 28 July and 14 August were higher than the discharges obtained during the 2015 flood event, suggesting more frequent and intense floods in the future. Most of the increase in discharge could be ascribed to two mechanisms. The first is the linear increase of sensible heat flux with air temperature (Eqn (2)). This effect was more significant under windy conditions, a conclusion supported by the peak discharge obtained on 28 July. The importance of turbulent heat flux in intensive melt events was demonstrated by the surface energy budget over the ablation zone of the Greenland ice sheet during the exceptionally large melt episodes in July 2012 (Fausto and others, 2016). Second, liquid precipitation increased – particularly in the high elevation area – following the conditions given by Eqns (13) and (14). For example, in the 400–500 m a.s.l. elevation band, there was an increase in the total summer liquid precipitation by 106% (i.e. 86.2–177.6 mm).

Fig. 12. Simulated discharge for the conditions observed in 2017 (blue) and for a 2°C warming (red). The horizontal dashed line represents the simulated peak discharge obtained during the flood event that occurred on 21 July 2015.

Further experiments were performed by imposing a 0.5–4°C temperature rise on the meteorological conditions observed in the 2015–19 period (Fig. 13). Relative to the mean of the experimental results based on the 2015–19 observations, a threefold increase in the summer discharge was projected (i.e. 9.5 × 10⁶ to 28.1 × 10⁶ m³) under a 4°C warming. Approximately 60% of the discharge increase resulted from ice melt, which increased nonlinearly as air temperature increased (Fig. 13), showing change rates of 2.3 × 10⁶ m³ K⁻¹ and 3.1 × 10⁶ m³ K⁻¹ for a change in temperature from 0 to 2°C and 2 to 4°C, respectively. The variation of ice melt with air temperature was nonlinear because a greater portion of the ice surface was exposed after the up-glacier migration of the snowline. The experiments performed using 2018 observation data show a more pronounced nonlinear increase in the total discharge (3.5 × 10⁶ m³ K⁻¹ and 5.1 × 10⁶ m³ K⁻¹ for a change in temperature from 0 to 2°C and from 2 to 4°C, respectively) as shown in Figure 13. The reason for this is that in 2018, the snowline migrated up-glacier over a large distance from a relatively low elevation. The variations of snowmelt and precipitation with increasing temperature were more linear and increased at rates of 1.6 × 10⁶ and 0.3 × 10⁶ m³ K⁻¹, respectively, for a change in temperature from 0 to 4°C. In 2016, a significant increase in discharge was obtained, owing to an increase in temperature from 0 to 2°C (Fig. 13) given that a greater level of precipitation was observed that year (Table 1).

Fig. 13. Simulated total summer discharge resulting from the imposition of a 0–4°C warming on the observed meteorological conditions in 2015 (+), 2016 (×), 2017 (▴), 2018 (▾), 2019 (■). The closed circles represent the means of the overall results (2015–19). The mean total discharge is partitioned into rain (orange), snowmelt (light blue) and ice melt (dark blue).

Our experiments quantified the increase in peak discharges as well as the total summer discharges on Qaanaaq Glacier, owing to warming projected over the next 100 years. Similar changes are expected for other glaciers and ice caps in Greenland and the Arctic, and the quantified discharges provide information that can be used to put measures in place that mitigate the impact of climate warming on communities in Greenland.

Conclusions

To investigate the underlying mechanisms associated with the glacial flood events that occurred in Qaanaaq, northwestern Greenland, field measurements on Qaanaaq Glacier and its outlet stream were made during the summer months of the 2016–19 period. Based on the field data, a glacier runoff routing model was developed by combining surface energy balance and linear reservoir models.

Based on an NSC of 0.71, the results of the glacier runoff routing model agreed sufficiently with the discharge observed within the 2017–19 period. The peak discharges reproduced using the model for the 2015 and 2016 flood events were 9.1 and 19.9 m³ s⁻¹, respectively. Both the model and the observed meteorological data showed that the 2015 flood resulted from intensive glacier melting (51 mm w.e. d⁻¹) that occurred on 21 July 2015 owing to warmer air temperature and strong wind, while the 2016 flood resulted from heavy rainfall (90 mm d⁻¹) that occurred on 2 August 2016, suggesting that wetter conditions in the Arctic will lead to increased risk of similar disasters in the future.

The developed model was used to quantify future glacier discharges in a warming climate. The numerical experiments showed a nonlinear increase in discharge with increasing temperature, represented by a threefold increase in the total summer discharge, owing to a 4°C increase in temperature, relative to the observed measurements. The projected increase in discharge could be attributed to enhanced melt and change of precipitation from snow to rain. Therefore, in this study, the influence of changing climate on the glaciers and outlet streams along the coast of Greenland was quantified and results obtained contribute to a better understanding of the impact of a changing climate on human society at large.