1. Introduction
Perennial land ice masses consist of a broad spectrum of sizes, from perennial ice patches (10–2 km2) to ice sheets (107 km2), and these have a range of responses and sensitivities to climate fluctuations (White and Copland, Reference White and Copland2018). Given the ongoing acceleration of climate change, understanding the spatio-temporal dynamics of the cryosphere is a research priority. However, research efforts have focused on the study of larger systems such as glaciers and ice caps, while smaller ice masses have attracted little attention. Ice patches are therefore among the least understood elements of the polar cryosphere despite their widespread distribution across the landscape (Lewkowicz and Harry, Reference Lewkowicz and Harry1991; Woo and Young, Reference Woo and Young2014; Young and others, Reference Young, Brown and Labine2018). Ice patches are a primary contributor to the hydrological cycle in sustaining meltwater flow in summer (Lewkowicz and Young, Reference Lewkowicz and Young1990; Young and others, Reference Young, Woo and Edlund1997; Assini and Young, Reference Assini and Young2012; Young and others, Reference Young, Assini, Abnizova and Miller2013, Reference Young, Scheffel, Abnizova and Siferd2017) and a driver of landscape and ecosystems development (Christiansen, Reference Christiansen1996; Woo and Young, Reference Woo and Young2003). Hence, the disappearance of ice patches in response to climate change has been identified as a critical issue for the stability of polar geosystems (Abnizova and Young, Reference Abnizova and Young2008; Woo and Young, Reference Woo and Young2014).
Ice patch systems represent an intermediate feature between seasonal snow patches and glaciers. The system includes the perennial ice body made of superimposed ice overlain by firn and a deep seasonal snow layer (Kawashima and others, Reference Kawashima, Yamada and Wakahama1993; Davesne and others, Reference Davesne, Domine and Fortier2022). These dynamics are controlled by the net surface mass balance, i.e. the difference between inputs and losses of mass over a hydrological year. Unlike glaciers, they are static and have no distinct accumulation and ablation zones (Glazirin and others, Reference Glazirin, Kodama and Ohata2004; Meulendyk and others, Reference Meulendyk, Moorman, Andrews and MacKay2012). Inputs of mass to the systems primarily consist of seasonal snow accumulation which depends on the precipitation at the regional scale and redistribution processes driven by wind. Ice patches are typically found in sites with winter snow accumulations of several meters or more, either by wind-driven snowdrift in sheltered gullies, concavities and leeward slopes (Watson and others, Reference Watson, Davison and French1994; Brown and Ward, Reference Brown and Ward1996; Woo and Young, Reference Woo and Young2014) or by avalanching at the bottom of steep slopes (Mott and others, Reference Mott2019). In polar regions, the combination of treeless terrain, strong winds and snow dryness makes wind transport the dominant process of snow redistribution (Eveland and others, Reference Eveland, Gooseff, Lampkin, Barrett and Takacs-Vesbach2013a, Reference Eveland, Gooseff, Lampkin, Barrett and Takacs-Vesbach2013b; Parr and others, Reference Parr, Sturm and Larsen2020). In spring and summer, part of the snow accumulation is converted into superimposed ice when in situ snow meltwater and upslope runoff percolate and refreeze within the snow and underlying firn layer (Ballantyne, Reference Ballantyne1978; Davesne and others, Reference Davesne, Domine and Fortier2022). Ablation, on the other hand, occurs through the sublimation and melting of snow/ice which are a function of the surface energy balance (SEB). The SEB is primarily controlled by the regional climate which depends on the latitudinal position, altitude, aspect, synoptic circulation patterns and land-sea contrasts (Alt, Reference Alt1987; Gardner and Sharp, Reference Gardner and Sharp2007). These regional patterns are modulated at finer scales by topography which affects micrometeorological variables (Mott and others, Reference Mott2011, Reference Mott, Vionnet and Grünewald2018; Schirmer and Pomeroy, Reference Schirmer and Pomeroy2020), and by snow/ice properties (Domine and others, Reference Domine, Gallet, Bock and Morin2012; Jennings and others, Reference Jennings, Kittel and Molotch2018). In mid-latitude mountains, extensive field and modeling work on the ablation of snow/ice patches (e.g. Canadian Rockies: DeBeer and Pomeroy, Reference DeBeer and Pomeroy2017; Schirmer and Pomeroy, Reference Schirmer and Pomeroy2020; Swiss Alps: Grünewald and others, Reference Grünewald, Schirmer, Mott and Lehning2010; Mott and others, Reference Mott2019) has highlighted a set of local effects brought by the topography and the size of the snow/ice patches themselves on microclimate conditions and heat exchange processes. It includes topographic shading (DeBeer and Sharp, Reference DeBeer and Sharp2009; Fujita and others, Reference Fujita, Hiyama, Iida and Ageta2010; Mott and others, Reference Mott2019), development of katabatic flow and turbulent heat fluxes (Mott and others, Reference Mott, Daniels and Lehning2015; Schlögl and others, Reference Schlögl, Lehning and Mott2018), cold-air pooling (Mott and others, Reference Mott, Gromke, Grünewald and Lehning2013) and advective heat fluxes from the adjacent bare ground (Neumann and Marsh, Reference Neumann and Marsh1998; Essery and others, Reference Essery, Granger and Pomeroy2006; Mott and others, Reference Mott2011; Harder and others, Reference Harder, Pomeroy and Helgason2017).
In polar regions, except for a few existing studies, our understanding of both regional and local drivers in ice patch ablation processes requires further research. In a pioneer study on an ice patch on Melville Island in the southern Canadian High Arctic, Young and Lewkowicz (Reference Young and Lewkowicz1990) demonstrated the large dominance of net solar radiation in the summer SEB due to the south-facing position of the ice patch and to the presence of an aeolian sediment cover at its top and bottom that significantly decreased albedo. The authors found that the sensible heat fluxes were dampened by the sheltered position of the ice patch located in a break of slope which reduced wind exposure and fog. More recently, a second study carried out by Leppäranta and others (Reference Leppäranta, Järvinen and Lindgren2013) on a snow patch at Basen Nunatak in western Dronning Maud Land in Antarctica highlighted the importance of sublimation processes which represent a major loss of mass and heat sink because of the aridity and the coldness of the climate. These earlier studies demonstrated the importance of considering specific regional drivers and local synoptic effects to understand ablation processes at the scale of an ice patch system.
To further improve the understanding of ice patches dynamics in polar desert conditions, which are facing major changes due to the ongoing climate warming (Vincent and others, Reference Vincent2023), we address the mass and SEB of a perennial ice patch on the northern coast of Ellesmere Island where semi-permanent and perennial ice patches are ubiquitous (Vincent and others, Reference Vincent2011). The climate of this region is influenced by the Arctic Ocean which maintains cooler and moister conditions in summer than in the interior Arctic (Maxwell, Reference Maxwell1981). We hypothesize that both fog and wind have an important control on energy-exchange processes at the snow/ice surfaces, the former by limiting the amount of incoming solar radiation and the latter by increasing turbulent heat fluxes (Essery and others, Reference Essery, Granger and Pomeroy2006; Grünewald and others, Reference Grünewald, Schirmer, Mott and Lehning2010; Mott and others, Reference Mott, Gromke, Grünewald and Lehning2013, Reference Mott, Daniels and Lehning2015). Wind conditions may also be crucial in winter by controlling snow redistribution processes and snow physical properties, which both impact the mass balance of ice patches (Eveland and others, Reference Eveland, Gooseff, Lampkin, Barrett and Takacs-Vesbach2013a, Reference Eveland, Gooseff, Lampkin, Barrett and Takacs-Vesbach2013b; Parr and others, Reference Parr, Sturm and Larsen2020; Davesne and others, Reference Davesne, Domine and Fortier2021). Our experimental work consisted of measuring the mass and energy balance of an ice patch at Ward Hunt Island (WHI), at the northern tip of Ellesmere Island. We specifically address the following objectives: (1) quantify the winter snow accumulation regime and summer ablation rate; (2) identify the regional and local drivers involved in the ice patch energy and mass balance and (3) assess the sensitivity of the ice patch to climate change.
2. Study site
The northern coast of Ellesmere Island (Figs 1a, b) experiences a polar desert climate characterized by temperatures that drop well below −30°C in winter and rarely rise above 5°C in summer, low annual precipitation (<250 mm a–1), snowpack on the ground 9‒10 months of the year, annual net solar radiation <100 MJ m–2 a–1, and an extreme seasonality with complete darkness in winter and 24 h-daylight during summer (Maxwell, Reference Maxwell1981; Barry and Hall-McKim, Reference Barry and Hall-McKim2018). At the synoptic scale, the north coast of Ellesmere Island is influenced by maritime air masses from the Arctic Ocean, and by the Arctic Cordillera, which act as barriers to the regional airflow and moisture circulation from the west and northwest. In summer, cold and moist air moving from the Arctic Ocean is blocked along the coast, which favors the formation of persistent fog and low stratus clouds, especially in July and August when the sea ice opens. The climate of the coastal fringe contrasts with the interior of Ellesmere such as at Eureka (80°N, 86°W) where warmer and drier summer conditions prevail (Table 1) (Maxwell, Reference Maxwell1981; Edlund and Alt, Reference Edlund and Alt1989; Barry and Hall-McKim, Reference Barry and Hall-McKim2018). This ‘Arctic Ocean effect’ significantly reduces the summer melt along the coast of Ellesmere Island and causes a local lowering of the glacier equilibrium line altitudes (Miller and others, Reference Miller, Bradley and Andrews1975; Wolken and others, Reference Wolken, England and Dyke2008), which favors the existence of low-elevation, coastal glacier and ice caps, ice rises and ice shelves, and the persistence of numerous snowbanks and ice patches (Jeffries, Reference Jeffries1982; Braun and others, Reference Braun, Hardy, Bradley and Sahanatien2004; Vincent and others, Reference Vincent2011, Reference Vincent2023).
(a) Map showing the ~1.2 million km2 of the ice-free land surface characterized as a polar desert in the circumpolar Arctic (annual precipitation <250 mm and mean July temperature <5°C). The black line delineates Ellesmere Island, which is the northernmost large island of the Canadian Arctic Archipelago. The black circles locate weather stations of the Arctic desert detailed in Table 1 ((1) Ward Hunt Island (Ca, NU); (2) Alert (Ca, NU); (3) Purple Valley (Ca, NU); (4) Eureka (Ca, NU); (5) Resolute Bay (Ca, NU); (6) Station Nord (Dk)). (b) Localization of Ward Hunt Island (WHI) 6 km north of the nearest shore of Ellesmere Island and about 760 km from the North Pole. (c) Topographic map of WHI (10 m contours) showing the Ward Hunt Lake (WHL), the Ward Hunt Ice Rise (WHIR), the research station and the SILA weather station managed by the Centre for Northern Studies (CEN). The ice patches (blue polygons) were delineated from GEOEYE satellite imagery taken in late summer 2011, one of the warmest summers in recent decades at WHI (Paquette and others, Reference Paquette, Fortier, Mueller, Sarrazin and Vincent2015). The red square indicates the location of the studied ice patch at the foot of the north face of Walker Hill.
Climate of Ellesmere Island's northern coast (Ward Hunt Island and Alert) relative to other polar desert localities of the Canadian High Arctic, Greenland and Antarctica for the period 1990–2019, except for Ward Hunt Island (2006–19), Purple Valley (2009–19) and Wright Valley (1994–2019)
MAAT, mean annual air temperature; MATs and MATw, mean summer and winter air temperature; RHS, mean annual precipitation and summer relative humidity. Summer averages are calculated for Jun-Jul-Aug for Northern Hemisphere sites and Dec-Jan-Feb for Southern Hemisphere sites, and conversely for winter averages.
Source: 1CEN, 2021; 2Environment Canada, 2021; 3Copland and Mueller, Reference Copland and Mueller2021; 4https://www.dmi.dk/vejrarkiv/normaler-gronland/; 5Doran and Fountain, Reference Doran and Fountain2021.
At WHI (83°05′09″N, 74°06′19″W) (Figs 1b, c), data provided by the SILA (meaning climate in Inuktitut) weather station provide information on the polar climate of the northern coast of Ellesmere (CEN, 2021). The mean annual air temperature and the average summer temperature were −17.1 and 0.4°C, respectively, for the period 2006–2019. Summer temperatures were therefore colder than at other stations in the High-Arctic but not as cold as in the McMurdo dry valleys in Antarctica (Table 1). The melt season is short, generally spanning from mid-June to the end of August and is characterized by high relative humidity (85.5% on average). Precipitation data are not available for WHI, but it is likely similar to the average of 148 mm water equivalent per year (w.e. a−1) measured at Alert, 170 km to the southeast (Environment Canada, 2021). Although WHI does not currently have glaciers, the cryosphere is omnipresent on and around the island (Figs 1b, c). The sea channel between the south coast of the island and Ellesmere is occupied by pack ice and remnants of the Ward Hunt Ice Shelf (Mueller and others, Reference Mueller, Vincent and Jeffries2003) while the north and west coast of WHI is bordered by the Ward Hunt Ice Rise (WHIR) which is an ice cap that was formed by the anchoring of the ice shelf on the sea floor about 1500 years ago (Braun and others, Reference Braun, Hardy, Bradley and Sahanatien2004). The island presents a cold and thick continuous permafrost, with a mean annual temperature of −13.5°C at 3 m depth according to the measurements provided by a thermistor cable for the period July 2015–July 2019 (CEN, 2021).
On the island, perennial ice patches constitute the most prominent cryospheric elements of the landscape. Our study focuses on an ice patch located at the base of the north face of Walker Hill (Fig. 1c). This ice patch (hereinafter IP1) lies between 40 and 60 m above sea level (a.s.l.) and is about 150 m long, 25 m wide and 3.2 m thick in its central section. IP1 is mostly made of superimposed ice whose physical properties vary according to the age of the ice. The coarse crystals that characterize the basement ice reflect the advanced age of the ice patch, which has likely been present for several millennia (Fortier and Davesne, Reference Fortier and Davesne2021a; Davesne and others, Reference Davesne, Domine and Fortier2022).
3. Instrumentation and data acquisition
3.1. Snow regime
At the SILA weather station (Fig. 1c; Table S1), the hourly snow height (SH in m) was recorded by a Campbell Scientific ultrasonic SR50 sensor from August 2006 to July 2019. The sensor was installed on a topographic break produced by a raised-beach ridge where a multi-year snow patch occurs. These data from a preferential accumulation site were used to infer the snow regime on IP1 since both sites have the same north/northeast orientation. To explore in detail the relationships between snowpack development and the wind conditions, we used the hourly wind data, i.e. maximum wind speed (V max in m s–1) and direction (DIR in °), recorded by an RM Young 05103-10 anemometer (Table S1). The maximum hourly wind speed was analyzed to capture wind gusts that have the potential to transport snow. These data were used to calculate the wind index (WI) which describes the sum of hourly wind speeds from September to May above the threshold of 6 m s‒1 for unconsolidated snow transport (Li and Pomeroy, Reference Li and Pomeroy1997; Sturm and others, Reference Sturm, Liston, Benson and Holmgren2001). Deposition/erosion events (i.e. a rapid change in SH of ±0.05 m) for each cold season from 2005 to 2019 were identified and their relation with wind conditions was examined.
3.2. Mass-balance variables
3.2.1. Ice patch thickness and volume
The total volume of IP1 was calculated by measuring the volume of seasonal snow and the ice body. Snow heights (SH in m) were manually probed every 2 m along ten georeferenced transects (Fig. 2) using a 3.50 m graduated rod. Manual observations of SH were interpolated across the ice patch with a kriging interpolation routine in ArcGIS, and integrated to obtain a total volume of snow for each survey. The horizontal resolution of the resulting snow height grids was 0.2 m. Surveys were carried out on 5 June 2016, 8 June 2017 and 12 June 2019 when the snowpack was at its maximum thickness according to the SR50 measurement at the SILA station. In 2019, an additional survey was conducted during the melt season on 5 July.
3D view of the studied ice patch (IP1). Red lines indicate the snow-depth measurement transects and blue dotted lines indicate the ground-penetrating radar (GPR). The black cross locates the automatic weather station installed during field season 2019, the red triangle shows the site where the air temperature was recorded in 2017 and the blue polygons indicate the location of the flumes in the downslope margin of IP1. The background image is an orthomosaic captured with an uncrewed air vehicle (UAV) on 17 July 2019.
The total volume of ice was estimated with ground-penetrating radar (GPR) surveys (Fig. 2) conducted on 7 July 2019 with 100 MHz antennas (Pulse EKKOpro System, Sensors and software®) (detailed in Appendix A). The resulting dataset was imported and interpolated by kriging in ArcGIS, providing a raster of the ice thickness (horizontal grid resolution of 0.5 m) from which a digital terrain model (DTM, horizontal grid resolution of 0.5 m) of the ice patch bed was produced (Ødegård and others, Reference Ødegård2017). The total volume of snow and ice (V total in m3) was then converted into water equivalent (w.e.) using Eqn (1):
where V i and V s are the volumes of ice and snow (m3); ρ i, ρ s and ρ w are the ice, snow and water densities (kg m‒3). We assumed average densities of ~905 kg m–3 for the ice body as reported by Davesne and others (Reference Davesne, Domine and Fortier2022). For the dry snow in early June, we also relied on average values reported by Davesne and others (Reference Davesne, Domine and Fortier2022) for the snowpack on IP1 in 2016 (ρ s = 442 kg m‒3) and 2017 (ρ s = 406 kg m‒3) and we made an additional snowpit excavation on 6 June 2019 following the method detailed in Domine and others (Reference Domine, Gallet, Bock and Morin2012). For melting snow, we used the density (ρ s = 534 kg m‒3) measured in the snowpack on IP1 on 2 July 2017 (Davesne and others, Reference Davesne, Domine and Fortier2022). Snow density can tend to increase as the summer progresses, so using a constant value may lead to a slight underestimation of our snowmelt values.
3.2.2. Daily ablation of the snow over the ice patch
The summer ablation on IP1 was monitored by combining high-resolution measurements of snow/ice surface changes using terrestrial laser scanning (TLS) (VX spatial station Trimble®) with measurements of the ice patch extent based on a continuous record of pictures taken by an automatic time-lapse camera (Fortier and Davesne, Reference Fortier and Davesne2021b; Appendix B). We performed daily repeated TLS surveys of the ice patch surface from 10 June to 12 July 2017 (32 days) and from 10 June to 20 July 2019 (40 days) (details in Appendix C and Fig. S1). The daily height change of the snow surface (ΔSHTLS in m d‒1) on the ice patch system was calculated by subtracting two consecutive scans. The results were then statistically validated by comparison with daily snow height data provided by the SILA station (ΔSHSILA in m) where the melt rate is expected to be similar to that of IP1 since the two sites are close and share the same orientation. ΔSHTLS was then converted into volume of water equivalent in ArcGIS using the Surface Volume tool (3D Analyst) and the aforementioned density value for melting snow, yielding the daily volume change (ΔV TLS in m3 d–1 w.e.). When ΔV is negative, it represents the daily melt rate (M in m3 d–1 w.e.) assuming no precipitation occurred during that day. The DTM of the IP1 bed measured on 7 July 2019 by GPR served as the reference to assess the daily evolution of the total volume of the ice patch.
3.3. Environmental monitoring during the melt season
3.3.1. Meteorological data
The analysis of the ablation processes relied on detailed knowledge of the micrometeorological conditions on IP1. In 2017, only the air temperature (T a in °C) was recorded at hourly intervals near the frontal edge of IP1 by a sensor (U22, Hobo®) installed in a radiation shield 1.5 m above the ground-surface (Fig. 2; Table S1). The wind conditions were extrapolated from the wind data recorded at the SILA station (Appendix D). In 2019, an automatic weather station (AWS) was installed on IP1 from 10 June to 19 July 2019 (Fig. 2). The AWS monitored T a, relative humidity (RH in %), wind speed (WS in ms–1), wind direction (DIR, °), incoming shortwave radiation (SW↓ in W m–2), net radiation (Q* in W m–2) and atmospheric pressure (AP in mbar). The station consisted of a 5 m-high steel mast inserted into a hole drilled into 2 m of snow and 50 cm of underlying ice (Table S1). T a and RH sensors were installed in a radiation shield and data were sampled every 15 min and recorded by a data logger (Micro station, Hobo®). SW↓, Q* and WS/DIR were sampled every 30 s, and the 15 min averages were recorded by a Campbell Scientific CR1000 datalogger. The pyranometer and net radiometer were installed parallel to the surface of the ice patch (~20°). The height of all sensors was adjusted as the melting progressed to keep the sensor's orientations and heights as constant as possible.
In addition to the AWS, the snow temperature (T snow in °C) was monitored by a thermistor cable inserted into the 2.3 m-thick snowpack from 7 June to 19 July (Table S1). It consisted of a series of seven thermistors at 190, 150, 120, 90, 60, 30 and 5 cm above the ice surface connected to a data logger (Systems Smart Reader Plus 8, ACR system) that recorded temperatures at hourly intervals. Liquid precipitation was measured using a metric rain gauge (resolution 0.2 mm) installed near the base camp and read each day at 20:00 h. Snowfall was negligible. For both study seasons 2017 and 2019, the dominant daily sky conditions (fog, overcast, clear sky) were inferred from the visual inspection of time-lapse photographs.
3.3.2. The hydrological regime of the ice patch
The outflow discharge of meltwater was monitored by three flumes (1 m long and 40 cm high) installed in major streams that drained IP1 (Fig. 2). The flumes were equipped with a pressure sensor (Hobo U20; Table S1) that continuously recorded the water level (h u in cm) (Paquette and others, Reference Paquette, Fortier and Vincent2017). h u was then converted to free-flow discharge (Q f in cm3 s–1) (Siddiqui and others, Reference Siddiqui, Lashari and Skogerboe1996). The sums of the flows measured at the three flumes were noted Q f total. This dataset was used as a proxy to infer the fluctuation in the outflow hydrograph during the study period and to correlate it with the melt rate.
3.4 Surface energy balance
As IP1 was covered by a thick seasonal snowpack throughout the study period, we used the standard approach for the computation of the SEB of a melting snow surface (Price and Dunne, Reference Price and Dunne1976; Gray and Male, Reference Gray and Male1981; Harder and others, Reference Harder, Helgason and Pomeroy2018). The SEB is given by Eqn (2):
where Q m is the net energy flux at the snow surface, Q* is the net radiation, Q h and Q e are the sensible and latent turbulent heat exchange respectively, Q c is the conductive heat flux to the underlying snow and ice and Q p is heat flux due to precipitation. All energy terms are expressed in MJ m–2 d–1 and are positive when they contribute energy to the snow surface. They were calculated on an hourly basis and then averaged to obtain the daily energy fluxes. Since the slope of the snow surface was about 20° during the study period and Q* was measured parallel to the surface and not vertically (see section 3.3), the SEB was calculated parallel to the snow surface. Turbulent energy exchange Q h and Q e were estimated by the bulk aerodynamic approach based on the gradient of temperature, vapor pressure and wind speed between the snow surface and the lower atmosphere (e.g. Boike and others, Reference Boike, Roth and Ippisch2003a). The SEB was calculated from hourly mean AWS measurements of Q*, T a, RH, WS. Equations for Q h, Q e and Q c are detailed in Appendix E.
For deep and cold snowpacks, modeling the melt rate in the early warm season is challenging because the surface starts to melt during the ripening phase before the cold content (CC in MJ m‒2) of the whole snow column is exhausted. In order to avoid an underestimation of the surface melt, we considered a top layer of 0.25 m (Garen and Marks, Reference Garen and Marks2005). We assumed that part of Q m is first used to change the energy stored within the top snow layer by raising its temperature (T snow). The melt energy (Q m) then starts to produce the melt of the snow surface (ΔSWESEB in mm d−1 w.e.) when the top snow layer becomes isothermal, i.e at the melting point (T m) (cold content = 0), following Eqn (3) (Oke, Reference Oke1987; Young and others, Reference Young, Assini, Abnizova and Miller2013):
where ρ w is the density of liquid water and L f is the latent heat of fusion (334 kJ kg‒1). ΔSWESEB was then corrected [÷ cos(20°)] to give the surface melt along the vertical plane to be consistent with TLS measurements.
The daily cold content (CC) of the snow layer is calculated from Dingman (Reference Dingman2002):
where c i is the heat capacity of the ice (2102 J kg–1 K–1), ρ s is the density of snow derived from the snowpit dug on 6 June 2019 and h is the height of the snow layer. T snow is the average temperature of the snow layer derived from thermistors and T m is the temperature of melting (i.e. 0°C).
4. Results
4.1. Snow and wind regimes
At the SILA snow gauge, the mean end-of-winter snowpack thickness (SHmax) between 2006 and 2019 was 1.52 m, and was characterized by a low interannual variability (std dev. = 0.25 m) (Fig. 3a). Only 2018 and 2019 stand out with a snowpack close to 2 m. In 2014 and 2018, the snow patch persisted throughout the melt season. On average, cold season snow accumulation events greater than 0.05 m occurred 7.5 ± 2.3 times over the 2006–19 observation period. 100% of the deposition events were associated with maximum hourly winds (V max) above 6 m s‒1, indicating drifting as a key process in the snowpack build-up.
Snow and wind regimes recorded at the SILA station from 2006 to 2019. (a) Time-series of the daily snow height (incomplete data in 2010/11, 2011/12, 2014/15) and wind index (WI) of the cold season (Sept–May) (red crosses: incomplete wind data in 2006/07, 2007/08, 2010/11 and 2011/12). The green line represents the mean maximum snow height from 2006 to 2019. (b) Wind rose for the cold season (Sept–May) constructed from hourly wind data for the 2006–2019 period. (c) Close up on cold seasons 2015/16, 2016/17, 2017/18 and 2018/19 showing the hourly snow height and daily WI.
Deposition events were also usually associated with subsequent erosion events that partially, and sometimes completely, removed fresh snow deposits within hours (e.g. 17 October 2015; Fig. 3c, Fig. S2). The net mass balance between deposition and erosion of each snowdrift event was largely positive in the early cold season (Fig. 4a). On average, maximum gains occurred in September (+0.40 m on average) and they progressively tended to zero as the season progressed and the topographic hollow was filled. On average, ~78% of SHmax was already accumulated by the end of December (Fig. 4b). Nearly 75% of the wind events >6 m s‒1 were associated with snow accumulation events (ΔSH > 5 cm) in October while this fraction dropped to 40% later in the season (Fig. 4b). The most extreme snowdrift events, generally in September, can bring a considerable amount of snow over a short time (e.g. 7 September 2011: +1.16 m; 19 September 2007: +0.9 m; 11 September 2009: +0.75 m). The prevailing winter wind direction for the speed range 5–10 m s‒1 was from the south-southwest and east-southeast, while strong winds (>10 m s‒1) were mainly from the west. This same pattern repeated year after year but the intensity of wind events was subject to substantial interannual variation as revealed by the wind index (WI) (Figs 3a, b). The windiest cold season was recorded in 2016/17 (WI = 21 410) while the calmest was in 2013/14 (WI = 10 600).
Average evolution of the snow regime throughout the cold season in the snow patch at the SILA station for the period 2006–19. (a) Mean monthly snow gains and losses (including erosion and compaction) and the resulting net change in snow height and (b) mean monthly snow height expressed as a proportion of the total end-of-winter snow height (SHmax in m) and the number of snow accumulation events (ΔSH > 0.05 m) as a proportion of the total number of potential blowing snow events (V max > 6 ms–1). Vertical bars indicate the std dev.
4.2. Mass-balance variables at the ice patch IP1
4.2.1. Seasonal snow accumulation
Figure 5 shows the snow conditions on and around IP1 at the beginning of melt seasons 2016, 2017 and 2019. Photographs and snow height measurements show large interannual variability in the snow distributions. Maximum snow height exceeded 2 m on IP1 for the 3 years studied. However, the spatial extent of the deep snow (>1.50 m) was variable between years. In summer 2017, deep snow was restricted to a narrow band along the upper ridge of IP1, while in summer 2019, it evenly covered the entire ice patch.
Time-lapse photographs of the ice patch system (top panels) and maps of the end-of-winter snow height (SHmax in m) in early June of 2016, 2017 and 2019, with corresponding average snow column density (ρ s in kg m–3), average snow height (SH in m) and total snow volume (V snow in m3 w.e.). The black line delimits the contour of the ice body as it was at the end of each summer.
Some differences were also observed in snow density: average early summer snow column densities were 442 kg m‒3 in 2016, 406 kg m‒3 in 2017 and 395 kg m‒3 in 2019 (Fig. S3). Differences in snowpack thickness and, to a lesser extent, in the average density of the snow, explain the significant interannual variability in the total volume of accumulated snow (V snow in m3 w.e., Fig. 5).
4.2.2. Ice thickness
The strongest GPR reflections (red dotted line in Fig. 6a) reveal the contact between the ice and the underlying substrate (rock or sediment). The uppermost reflection (blue dotted line) corresponds to the interface between the seasonal snowpack and the ice surface. IP1 bed topography is concave, with a steep slope in the upper section and a slope break ~45 m from the head of the ice patch. In July 2019, the maximum ice thickness was 3.8 m, and the total volume of ice was estimated at 14 100 m3 (12 500 m3 w.e.). Including the remaining snowpack above as measured on 5 July 2019 (V snow = 10 500 m3, giving 5 500 m3 w.e.; Fig. S4), the total volume of the ice patch system (V total) was 24 600 m3 for a w.e. of 19 600 m3.
(a) 2D visualization of a GPR profile made in the central part of the ice patch in 2019 using 100 MHz antennas (the track is represented by the black line in (b). (b) Interpolated map of the ice thickness based on the series of GPR tracks made in July 2019 (see Fig. 2 for all GPR tracks location). The black line in (b) delimits the contour of the ice body at the end of summer 2016.
4.2.3. Snowmelt monitoring in summers 2017 and 2019
At the scale of IP1, the ablation spatial patterns (ΔSHTLS in m) for 2017 and 2019 are shown in Figure 7. The snow height decreased by 0.85 m between 10 June and 12 July 2017 and by 1.33 m between 10 June and 20 July 2019. The maps revealed strong spatial variability in the surface ablation over the ice patch with a similar pattern in both years. Ablation was strongest just below the upper edge and on the west side of the ice patch (zones 1 and 2) while it was slightly below the average value in the central and east parts (zones 3 and 4).
Maps of the total surface ablation (ΔSHTLS in m) during the study periods 2017 and 2019 based on daily terrestrial laser scan (TLS) surveys of snow/ice surface position. The black boxes represent the averaging area used to highlight the spatial variability of the ablation. Mean (±std. dev.) values for the whole ice patch and each zone are presented below each panel.
Figure 8 shows the daily evolution of snowmelt on IP1 during 2017 and 2019. In 2017, the seasonal snowpack was already discontinuous on 10 June and declined rapidly to disappear almost completely on 25 June. The total volume of the ice patch system, which combines the volume of ice and seasonal snow, was ~25 500 m3 at the beginning of the melt season (Fig. 8a). The melt rate was low between 10 and 16 June 2017 (Fig. 8a) and then increased rapidly to reach a first peak on 23 June (−950 m3 d‒1 w.e.). Meltwater started to flow out at this time at the flumes (Fig. 8b). On 28 June, a major snowstorm re-filled the ice patch system with ~1 850 m3 w.e. of new snow and temporally interrupted the water outflow. The rapid melt of this snow in the following days led to a peak in ablation on 30 June (1 100 m3 d‒1 w.e.). The highest ablation occurred on 6 July (1 200 m3 d‒1 w.e.), leading to a high peak in the meltwater discharge (Qftotal). At the end of the study period on 10 July 2017, the total volume of the ice patch was 9 500 m3 and its areal extent was 6 500 m2 (Figs 8a, c). Analysis of the time-lapse photograph dataset revealed that the minimum extent was reached on 11 August with 3 850 m2 (Fig. S5).
Monitoring of snowmelt from 10 June 2017 to 12 July 2017 (left side) and from 10 June to 20 July 2019 (right side): a and a’; daily evolution of the total ice patch system volume (snow and ice together) (V total in m3) along with the daily volume change in water equivalent (ΔV in m3 d−1 w.e.) measured by terrestrial laser scan (TLS). The red bars (negative ΔV) give an estimation of the melt rate (M in m3 d‒1 w.e.); b and b’. Daily meltwater outflow discharge (Qftotal in m3 d‒1) measured at flumes; c and c’. Picture of IP1 illustrating the condition the system at the end of each study period.
In 2019, the decline of the seasonal snowpack was slower than in 2017 in the first half of June, but it accelerated considerably in the second half. Some seasonal snow patches around IP1 persisted until the end of the study period. The total volume of the IP1 system was ~38 000 m3 in early June, 50% higher than in 2017 (Fig. 8a’). The melt rate exhibited less pronounced daily variations and peak flows were weaker than in 2017 (Fig. 8b’). The highest daily ablation peak over the study period was reached on 9 July with a loss of 720 m3 d‒1 w.e., resulting in a peak flow of 217 m3 d‒1 at the flumes. On 20 July 2019, the ice patch extent was 10 300 m2 for a total volume of 16 200 m3, which represents a loss of 58% of its initial volume in 40 days, (Fig. 8c’).
In both 2017 and 2019, daily changes in snow height at IP1 and SILA weather station were highly correlated (Pearson's correlation test, r = 0.91, p < 0.01 in 2017 and r = 0.83, p < 0.01 in 2019; Fig. S6). The daily fluctuations of Qftotal show a strong correlation ΔSHTLS (m3 d‒1) on IP1 (Pearson's correlation test, r = 0.67, p < 0.01 in 2017 and r = 0.87, p < 0.0001 in 2019 (Fig. S7) except at the beginning of the melt season and during the snowstorm event in 2017 due to the lag in the response of the hydrograph to the surface melt rate.
4.3. Surface energy balance
4.3.1. Meteorological conditions
An overview of the meteorological and environmental variables for summers 2017 and 2019 is given in Figure 9.
Summary of hourly measured meteorological variables during the study periods. (a) 2017 air temperature (T a in m) and mean wind speed (WS in m s–1). (b) 2019 T a and relative humidity (RH in %); (c) 2019 WS and daily precipitation; (d) 2019 incoming solar radiation (SW↓ in W m–2) and net radiation (Q* in W m–2); and (e) contour map snow temperature (T snow in °C) in 2019 with snow height evolution (dashed black line) at the thermistor chain site on IP1 derived from TLS surveys. The background colors indicate the sky conditions; white for fog, grey for overcast and blue for clear sky.
In summer 2017 (Fig. 9a), the mean daily air temperature (T a), recorded in the downslope margin of the ice patch fluctuated from –4 to 9°C, with an average of 1.7°C for the whole period. The mean hourly wind speed (based on SILA measurements, modified with correlations; Appendix D) was 3 m s−1 over the study period with a maximum value of 8 m s‒1, with a prevailing direction from the west. There was fog, cloud cover and clear sky, 10, 65 and 25% of the time, respectively.
In summer 2019 (Fig. 9b), T a at 0.5 and 2 m above the snow surface had an average of 1.7 and 2.4°C, respectively, over the study period. The minimum T a (−2°C) was recorded at 0.5 m during foggy and calm conditions on 2 July while the highest value (10.1°C) was recorded during warm and dry spells associated with strong winds on July 13. The relative humidity (RH) was generally high, reflecting the common foggy conditions, and fluctuated in opposite phase to the air temperature.
In 2019, the mean hourly wind speed was 1.9 m s‒1 with a maximum value of 7 m s‒1, and a prevailing direction from the west (Fig. 9c). There was fog, clouds and clear-sky 30, 30 and 40 of the time, respectively. Precipitation fell essentially as rain and drizzle and reached a cumulative value of 7 mm. During clear-sky conditions, the incoming solar radiation (SW↓) followed a diurnal cycle with maxima >500 W m2 and minima <150 W m2. Net radiation (Q*) had maxima between 110 and 175 W m2 and minima as low as –50 W m2 (Fig. 9d). Finally, the thermistors chain showed that the snow temperature (T snow) in the top layer of the snowpack reached the melting point on 13 June while the temperature was ~−13°C near the snow/ice interface (height of 5 cm). T snow gradually increased at first from the snow surface. After mid-June, latent heat release at or near the snow–ice interface due to liquid water percolation and refreezing also warmed up the base of the snowpack.
4.3.2. Daily mean evolution of the SEB terms in 2019
The daily variability of the SEB is presented in Figure 10a. Both net surface radiation Q* and the sensible heat flux Q h were sources of heat to the snow surface throughout the period. Q* varied from 1.4 to 6 MJ m–2 d−1, increasing with decreasing cloudiness. The sensible heat flux Q h varied from 0 to 6.5 MJ m–2 d−1. The highest values of Q h were observed during windy and warm days (e.g. 9–10 and 17 July). The latent heat flux Q e varied between −1.3 and 2.3 MJ m–2 d−1. Positive values indicate condensation, peaking on 9–10 July, while evaporation/sublimation resulting in surface heat loss was observed on June 10 and 11. Heat loss by conduction (Q c) at the snow surface was small (−0.15 MJ m–2 d−1) and was limited to the first 2 days of the period as the cold content (CC) of the surface layer (0–0.25 m below snow surface) was rapidly eliminated (Fig. 10a). The sum of all SEB terms gives a net gain of heat (Q m) throughout the period, peaking on 9 July with 12 MJ m–2 d−1.
(a) Daily mean values of energy balance terms with net radiation (Q* in MJ m–2 d–1), sensible heat flux (Q h in MJ m–2 d–1), latent heat flux (Q e in MJ m–2 d–1), heat flux by conduction to the underlying ice (Q c in MJ m–2 d–1), the cold content (CC in MJ m–2) in the 0–0.25 m surface layer, and the melt flux (Q m in MJ m–2 d–1). The background colors indicate the sky condition; white for fog, grey for overcast and blue for clear-sky; (b) comparison between daily melt rate calculated from the surface energy balance (ΔSWESEB in mm d–1 w.e.) and measured by TLS (ΔSHTLS in mm d–1 w.e.). RMSE means root mean square error.
On average, Q* provided 74% of total energy gains to the ice patch while 21.5% came from Q h and 4.5% from Q e. However, an analysis of variance indicated that 77% of the variance in the Q m is explained by Q h. The daily TLS melt rate, ΔSHTLS, shows good agreement with the SEB melt rate, ΔSWESEB, with a root mean square error of 3.54 mm d−1 and a determination coefficient of 0.68 (Fig. 10b).
The contribution of Q* is, however, less significant than the value reported by Young and Lewkowicz (Reference Young and Lewkowicz1990) on a semi-perennial snow patch on Melville Island where Q* accounted for 85% of the total energy gains over a period from June 22 to August 10. This may be explained in part by a high surface albedo on IP1 where dust concentration on the snow surface was low during the study period.
Three examples of common summertime synoptic conditions along the northern coast of Ellesmere Island (Fig. 11) were selected: clear sky (28 June), fog (5 July) and overcast (9 July). A summary of measured meteorological and calculated SEB terms for each case is given in Table S2 and an overview of the synoptic conditions (i.e. 500 hPa geopotential heights and 850 hPa air temperature) obtained from the NCEP–NCAR reanalysis is presented in Figure S8:
(a) Case 1, clear sky (June 28; Fig. 11a): A high-pressure system (Fig. S8), 1029 hPa, resulted in clear sky conditions, cool air temperatures (T a−daily = 3.1°C at 2 m) and low winds (Table S2). SW↓ reached 550 W m–2 around 7:00 h (local time), corresponding to a peak in Q* of ~120 W m‒2. The combination of low westerly winds and relatively cold temperatures limited Q h. Q e remained negligible. On average, the ice patch system gained 4.5 MJ m–2 over this day, leading to melting of 12.5 mm d‒1 w.e. Q* was the dominant heat source (79%).
(b) Case 2, fog (5 July; Fig. 11b): Synoptic conditions were characterized by a shallow trough between two low-pressure systems (Fig. S8), giving calm conditions with weak winds and persistent fog. Ta remained near the freezing point (T a−daily = 0.3°C at 2 m respectively) (Table S2). Maximum SW↓ reached ~400 W m‒2 and Q* peaked at only 110 W m‒2 around 8:00 h (local time). Because of the calm and near-freezing conditions, both Q h and Q e remained low. The supply of heat to the ice patch was only 2.8 MJ m–2 d−1 on average, producing a melt of 8.2 mm d‒1 w.e. Q* contributed for 90% while Q h and Q e contributed for 2 and 8%, respectively.
(c) Case 3, overcast (9 July; Fig. 11c): That day was marked by the passage of a warm front linked to a low-pressure system that circulated further south over central Ellesmere Island and that advected warm air (>5°C at 850 hPa) from the southwest (Fig. S8). At WHI, it was cloudy and windy with light rain and drizzle. The 15 min average wind (WS) reached 7.9 m s‒1 with gusts exceeding 10 m s‒1 (Table S2). T a reached slightly above 8°C around 5:00 h (local time). Q h peaked at 180 W m‒2 and remained higher than Q* during most of the day. Q e remained weak with a maximum value of 40 W m‒2. The gain of energy to the ice patch was 11.9 MJ m–2 d−1, resulting in substantial melting of 37.2 mm d‒1 w.e. Q h was the main term of the SEB providing 57% of the total energy. Q* and Q e contributed about 27 and 16%, respectively.
Diurnal cycle of the meteorological variables, including the 15 min average air temperature (T a in °C) at 2 and 0.5 m above the snow surface, wind speed (WS in m s–1) and incoming solar radiation (SW↓ in W m–2), and the hourly evolution of surface energy balance terms, i.e. net radiations (Q* in W m–2), sensible heat fluxes (Q h in W m–2), latent heat fluxes (Q e in W m–2), the melt fluxes (Q m in W m–2) and time-lapse photos for the three different cases that reflected the dominant conditions encountered along the north coast of Ellesmere Island: (a) case 1 was characterized by sunny and calm conditions; (b) case 2 was a foggy day and (c) case 3 was marked by the influence of a low-pressure system with strong winds and warm air advection. The red cross in the picture locates the automatic weather station on the ice patch. Note that the daily solar radiation peak early in the morning due to the northeast orientation of the ice patch.
5. Discussion
5.1. Main mass gain to the ice patch system
Ice mass gain in IP1 occurred primarily through the aggradation of superimposed ice by the refreezing of meltwater upon the old ice surface at the beginning of the warm season (Davesne and others, Reference Davesne, Domine and Fortier2022). Although inflows of water from the upslope occur and contribute to feeding ice growth (Ballantyne, Reference Ballantyne1978; Lewkowicz and Harry, Reference Lewkowicz and Harry1991), seasonal snow accumulation represents the primary mass gain to the system (Lewkowicz and Young, Reference Lewkowicz and Young1990; Davesne and others, Reference Davesne, Domine and Fortier2022).
Along the north coast of Ellesmere Island, precipitation is limited during the cold season (<150 mm w.e. at Alert between September and May; Environment Canada, 2021; Fig. S10a). These low precipitation amounts are unevenly redistributed across the landscape by winds so that the overall snowpack in polar deserts is thin and discontinuous (<0.30 m; Domine and others, Reference Domine2018; Royer and others, Reference Royer2021; Davesne and others, Reference Davesne, Domine and Fortier2021), but locally significant amounts of snow can accumulate by wind action (Figs 3, 5). Since most of the topographic niches where ice patches have developed are relatively limited in volume, we consider them to belong to the class of so-called ‘filling’ snow traps that are filled by snowdrift fluxes by the end of winter, in contrast to ‘nonfilling’ snow traps where the potential topographic volume that can be filled exceeds what can be supplied by precipitation and snowdrift fluxes (Benson and Sturm, Reference Benson and Sturm1993; Parr and others, Reference Parr, Sturm and Larsen2020). As shown in Figure 4, the ‘filling’ snow traps filled up quickly at the beginning of the cold season (80% of SHmax was reached before January on average at SILA) due to frequent winds above the drifting snow threshold of 6 m s–1. At WHI, the fetch providing snow to IP1 is virtually unlimited because of the smooth surface of the surrounding cryosphere (i.e. WHIR, frozen Arctic Ocean), and this probably explains the rapid filling of IP1, as early as September–October (Fig. 3). Later in the winter, the retention capacity of the snow traps rapidly decreases as it reaches its equilibrium-drift surface (Tabler, Reference Tabler1975).
The equilibrium surface of IP1 is primarily controlled by a break of slope at the cliff base (Figs 6, 12). The interannual variability of maximum snow accumulation in topographic niches is little affected by the amount of winter precipitation (Fig. S10b). It is however strongly modulated by wind conditions, which show significant interannual variability, as highlighted by the wind index (Figs 3a, c; Davesne and others, Reference Davesne, Domine and Fortier2021). The contrasted snow conditions during 2016–19 at both SILA (Figs 3a, c) and IP1 (Fig. 5) illustrate the effects of varying winter wind intensities on the maximum volume of snow stored in snow traps. Figure 12 shows that SHmax is positively correlated with the number of wind events in the range of 4–10 m s–1, suitable for snow transport, saltation (Li and Pomeroy, Reference Li and Pomeroy1997) and trapping in sites such as IP1 (e.g. winters 2017/18; 2018/19; Figs 3, 5). Conversely, SHmax is inversely correlated with wind events above 15 m s–1, the threshold for snow transport by suspension (Pomeroy, Reference Pomeroy1989). Such winds are associated with high erosion and sublimation rates and are unfavorable to accumulation, even in sheltered areas, due to high turbulence (Pomeroy, Reference Pomeroy1989). It results in a lowering of the equilibrium surface which leads to a decrease of SHmax (e.g. winter 2016/17).
Scatter plots showing the relationship between the frequency of wind events and the maximum snow height (SHmax in m) observed at the SILA site (left y-axis) and the maximum snow volume (V snow in m3 w.e.) measured on IP1 (right y-axis). (a) For wind events in the range of 4–10 m s–1; (b) for wind events with speed >15 m s–1. The dashed black lines are linear regressions for the SILA station only.
An additional explanation for Figure 12 is the indirect effect of wind on snow redistribution through its influence on snow properties and its ability to be eroded. Under strong winds, snow grains are rapidly fragmented by mechanical metamorphism and can form hard dense wind slabs made of small and mostly rounded grains (Domine and others, Reference Domine, Cabanes and Legagneux2002). The hardness and density of this snow are generally proportional to the wind speed during the snow accumulation. Thus, stronger winds produce a snowpack with a higher threshold for erosion (Jaedicke and Sandvik, Reference Jaedicke and Sandvik2002) so that the snow will be less likely to be eroded during moderate wind events (4–10 m s–1) which limits snow supply. Although denser snow increases SWE, this does not balance out the reduced snow supply to the ice patch system (Fig. 5).
We thus suggest that the high frequency of moderate wind events (4–10 m s–1) represents the best conditions for snow trap filling (e.g. 2018; 2019, Figs 3, 5) while a high frequency of strong wind events (>15 m s–1) considerably limits accumulations even in the most sheltered sites (e.g. 2017; Figs 3, 5).
5.2. Factors controlling summer ablation
Snowmelt becomes significant in early summer as soon as the upper snow layers reach the melting point (Fig. 9) so that all energy inputs to the snow surface are almost exclusively consumed as latent heat by melting or sublimation. This represents the main energy sink in the SEB (Fig. 10). In 2019, the SEB was dominated by net radiation (74%), consistent with other studies on similar topics (e.g. Young and Lewkowicz, Reference Young and Lewkowicz1990; Boike and others, Reference Boike, Roth and Ippisch2003a, Reference Boike2003b; Westermann and others, Reference Westermann, Lüers, Langer, Piel and Boike2009; Leppäranta and others, Reference Leppäranta, Järvinen and Lindgren2013). The daily variability in net melt energy, however, was driven by the evolution of turbulent fluxes and in particular by the sensible heat flux (Q h). Turbulent fluxes are closely related to weather conditions (air temperature, relative humidity, wind speed) and are therefore expected to be highly variable in time and space and to be strongly impacted by climate fluctuation.
5.2.1. Synoptic control on surface energy balance and melt rate
Summer ablation totals depend on synoptic conditions that develop over the Arctic Basin. Three specific synoptic types (Fig. 11) suggest that fog and wind speed are key meteorological controls on ice patch ablation. Their occurrences were examined in the regional context to determine under what synoptic situations they typically occur along the northern coast of the Canadian High Arctic (Alt, Reference Alt1987; Table S3).
Fog (case 2; Figs 9, 11) tends to preserve ice patches since it reduces solar radiation inputs and near-surface air temperatures. Turbulent heat inputs under such conditions are negligible, though latent heat contributions from condensation on the surface and freezing of supercooled fog droplets may be slightly positive. Fog forms under high-pressure systems over the Arctic Ocean or shallow low patterns characterized by low wind and low atmospheric mixing, favoring persistent temperature inversions that trap humidity near the surface (synoptic Type I of Alt, Reference Alt1987) (Fig. S8; Pope and others, Reference Pope, Copland, Alt, Copland and Mueller2017). Along the northern coast of Ellesmere Island, fog typically forms by advection of oceanic air in the lowest layer of the atmosphere that was cooled and moisture-laden as it circulates over the fractured melting sea ice (early July 2019). The frequent fog considerably cools the summer climate of the coastal fringe (Figs S9a, b; Table 1; Maxwell, Reference Maxwell1981) producing the ‘Arctic Ocean effect’ (Edlund and Alt, Reference Edlund and Alt1989; Koerner, Reference Koerner2005; Braun, Reference Braun, Copland and Mueller2017) and also affects the local-scale SEB.
Windy conditions increase ice patch ablation by reducing the occurrence of fog (Fig. S9c), and by limiting the development of stable thermal inversion layers near the snow/ice surface that prevents turbulent heat transfer (Synoptic Type I of Alt, Reference Alt1987; case 1 in Fig. 11a and Table S2; Mott and others, Reference Mott2011). Wind enhances turbulent heat transfers, in particular when strong winds are combined with warm air advection (e.g. 6 July 2017, Fig. 9; case 3 in Fig. 11; Fig. S8). Low 500 hPA geopotential heights centerd on southern Ellesmere Island and high pressures over Greenland direct relatively warm (T a at 850 hPa >5°C) and moist air to northern Ellesmere Island from the southeast (Fig. S8; Type II to Type III classes, Alt, Reference Alt1987). The melt effectiveness of warm and wind conditions over late-lying snow patches and seasonal snow has been frequently reported in alpine (e.g. Mott and others, Reference Mott2011) and Arctic regions (e.g. Lewkowicz and Young, Reference Lewkowicz and Young1990).
5.2.2 Ice patch albedo
The exposure of ice, and its timing, will affect summer ablation totals due to albedo. The optical properties of the ice patch surface can vary greatly throughout the summer and from one year to the next (Fig. 13). When seasonal snow persists for a long time on the ice patch (e.g. 2018 and 2019; Fig. 13), the albedo remains high and favors the preservation of the underlying ice body. In contrast, summers characterized by severe ablation or reduced snow accumulation (e.g. 2017; Fig. S5) are marked by rapid exposure of superimposed ice, resulting in a drop in surface albedo (Gardner and Sharp, Reference Gardner and Sharp2010), which positively feeds back on melting. This shift in surface albedo is even more dramatic when a thin layer of impurities and sediment forms on the exposed ice (Davesne and others, Reference Davesne, Domine and Fortier2022), as it can reduce the surface albedo to 0.2 (Young and Lewkowicz, Reference Young and Lewkowicz1990; Fujita and others, Reference Fujita, Hiyama, Iida and Ageta2010). Summer snowstorms, such as in 2017 (Fig. 8), can temporarily increase surface albedo. However, this snow generally melts within a few days and, according to our field observation over the last 12 years, this kind of event is rather unusual. We thus consider the effect of summer snowfalls on the energy balance negligible during the warm season.
Time-lapse photographs of the ice patch during the melt seasons 2017 and 2018 showing a contrasting situation in terms of surface albedo. In 2017, the pronounced ablation led to progressive exposure of the dirty ice surface giving a dark coloration to the ice patch while in 2018 the seasonal snowpack persisted throughout the summer.
5.2.3. Heat advection processes
Repeated TLS surveys highlighted the small-scale variability in melt rate (Fig. 7) that reflects heat advection (i.e. lateral heat transport), either by wind (turbulent heat) or by running water.
In the case of IP1, the prevailing westerly wind produced local advection of sensible heat from the upwind blocky ground surfaces toward the western edge of the ice patch (Fig. 7). As the season progressed, the advection process became more pronounced due to the reduction of snow-covered areas around the ice patch (Fig. 8), which enhanced the turbulent exchange between the solar-heated ground and the air. As a result, the melt rate increased faster on the upwind edge than on the downwind edge (Fig. 7). Such local-scale heat transfers to patchy snowpacks have also widely been documented in the alpine context (e.g. Essery and others, Reference Essery, Granger and Pomeroy2006; Mott and others, Reference Mott2011, Reference Mott, Schlögl, Dirks and Lehning2017) and in the High-Arctic (Young and Lewkowicz, Reference Young and Lewkowicz1990; Hardy, Reference Hardy1996). Advection of heat by running water can also have a significant impact on local-scale ablation as explained by the higher melt rates measured along the upper edge of the ice patch (Fig. 7). As shown in Figure 9e, water circulation on the snow–ice interface provides a significant amount of heat to the ice patch system early in the summer. Later in the season, water flow greatly influences the spatial variability of local-scale ablation, as evidenced by higher melt rates measured along the upper edge of the ice body (Fig. 7) and thermal erosion (Fortier and others, Reference Fortier, Allard and Shur2007) produced by runoff over the ice surface that forms a series of linear channels, sometimes up to 1 m deep and 0.50 m wide, carved into the ice patch (Davesne and others, Reference Davesne, Domine and Fortier2022).
These local processes are strongly conditioned by the ground-surface thermal regime in the surrounding areas. Any change in the ground-surface properties (e.g. albedo, vegetation, hydrology, temperature) could therefore greatly modify the intensity of heat advection to the ice patch.
5.3. Sensitivity of the ice patches to climate change
5.3.1. Long-term evolution
The significance of ice patches for the polar landscape and ecosystem development (Christiansen, Reference Christiansen1996; Woo and Young, Reference Woo and Young2003; Woo and Young, Reference Woo and Young2014) makes their long-term preservation a key issue for the stability of the polar geosystem. Based on the analysis of ice cores retrieved from IP1, Davesne and others (Reference Davesne, Domine and Fortier2022) suggested that this ice patch is likely at least several centuries old, potentially dating back to the end of the Mid-Holocene Warm Period. As such, despite its limited size, IP1 is less sensitive to climate change than the nearby cryosphere which has experienced an unprecedented decline over the last two decades (Mueller and others, Reference Mueller, Vincent and Jeffries2003, Reference Mueller, Copland, Jeffries, Copland and Mueller2017; Braun and others, Reference Braun, Hardy, Bradley and Sahanatien2004; Paquette and others, Reference Paquette, Fortier, Mueller, Sarrazin and Vincent2015). The current preservation of ice patches until today, despite the recent warmer summers, is indicative of a partial decoupling of their dynamics from the climatic trend due to the primary control by the topographical constraints. This confirms similar findings on ice patches (Lewkowicz and Harry, Reference Lewkowicz and Harry1991; Fujita and others, Reference Fujita, Hiyama, Iida and Ageta2010; Ødegård and others, Reference Ødegård2017) and glacierets (Kuhn, Reference Kuhn1995; Hoffman and others, Reference Hoffman, Fountain and Achuff2007; DeBeer and Sharp, Reference DeBeer and Sharp2009).
We suggest that the key factor in maintaining a long-term mass balance close to zero at WHI is the development of a negative feedback between summer ablation and winter accumulation (Glazirin and others, Reference Glazirin, Kodama and Ohata2004). Indeed, for limited-fill ice patch systems fed by drifting snow, winter accumulations are generally negatively correlated with the end-of-summer ice patch size and volume because the amount of new snow that can accumulate in the system in winter depends on the volume available in the topographic niche. However, the 3 years of data in IP1 show that this relationship is not straightforward since snow availability is modulated by wind conditions (4–10 m s–1 winds) which can favor snow loading after a summer of low ablation (e.g. 2017/18, Fig. 5) or limit it (>15 m s–1 winds) after a summer of strong ablation (e.g. 2017/17, Fig. 12).
This negative feedback mechanism dampens the effects of climate variations on ice patch mass balance but does not suppress them since the summer ablation remains greatly influenced by the micro-meteorological conditions, especially air temperature, wind and fog.
5.3.2. Potential future evolution
According to ensemble simulations from the Coupled Model Intercomparison Project 5 (CMIP5), under the intermediate RCP 4.5 greenhouse gas concentration scenario (Taylor and others, Reference Taylor, Stouffer and Meehl2012), the mean annual temperature is projected to increase by 3°C by 2050 on the northern coast of Ellesmere Island (Climate Change Institute, 2021). Sea ice extent in the Arctic Basin would consequently continue to decline rapidly, greatly increasing turbulent fluxes from the open ocean to the atmosphere. This would enhance the rise in air temperature in fall (+4.5°C) (Fig. S11a), possibly modifying Arctic cyclone characteristics, and could result in more frequent and severe storms, especially in fall (Akperov and others, Reference Akperov2019; Oh and others, Reference Oh, Sushama and Teufel2020), as revealed by projected negative sea-level pressure anomaly on the northern coast of Ellesmere Island by the CMIP5 models (Fig. S11b). This phenomenon would be associated with a substantial increase in precipitation, especially in fall (+28%) and winter (+29%) by 2050 (Fig. S11c). Regional increase in cyclonic activity is also expected to result in more extreme surface winds in the coming decades (Mioduszewski and others, Reference Mioduszewski, Vavrus and Wang2018; Oh and others, Reference Oh, Sushama and Teufel2020).
These expected climate changes will produce contrasting effects on the snow regime at WHI. On the one hand, increased precipitation (Fig. S11c) could bring increased snow availability during the cold season with the potential of higher snowdrifting. On the other hand, the increase in cyclonic activity would potentially be accompanied by more frequent extreme wind events (>15 m s‒1), providing less favorable conditions for a complete snow filling of the topographic niches with situations similar to that observed in 2017 becoming more frequent. Projected increases in air temperature could alter snowpack conditions, particularly in the fall, with potentially wetter snow of lower erodibility. In addition, an increase in rain-on-snow events in winter is also projected across the High Arctic due to more frequent warm spells (Bintanja and Andry, Reference Bintanja and Andry2017). Although rain-on-snow events are unlikely at the latitude of WHI by 2050 during winter, they could occur in fall (September–October). A wetter snowpack with ice layers would then be less likely to be eroded and transported. Warmer conditions in winter and spring would also increase the temperature of the end-of-winter snow and the underlying ice body, reducing the snowpack cold content and therefore the duration of superimposed ice growth in early summer (Davesne and others, Reference Davesne, Domine and Fortier2022). Therefore, despite increased precipitation, we anticipate that climate change will result in a decline in mass inputs to the ice patch system.
Although summer warming is expected to be moderate by 2050 (+1°C) under scenario RCP4.5 (Fig. S11), this would result in longer and more intense melt seasons. In addition, increased air temperatures could potentially trigger several positive feedback loops that increase ice patch ablation. For example, complete meltout of seasonal snow will expose darker ice surfaces (e.g. 2017) and enhance ice patch ablation. Also, higher temperature of the ground-surface around the ice patches, due to warmer air temperature and earlier melting of the seasonal snowpack, could enhance heat advection to ice patches by way of thermal erosion from upslope runoff and lateral airflow, thereby further melting ice patch margins.
How fog and wind conditions will evolve under a warmer Arctic climate is also a key element in assessing the fate of ice patches. These two variables are closely related to each other because they both depend strongly on sea ice conditions in the Arctic Basin. First, an increase in open water surface will likely produce more frequent fog and thus preserve the coastal cryosphere, giving rise to a negative feedback that could counteract, or at least dampen, the effect of temperature warming. However, the expected enhanced cyclonic activity and stronger surface winds in summer would result in unfavorable conditions for fog formation (Fig. S9c). In addition to limiting fog, the more frequent wind in summer would increase turbulent heat transfer to ice patch systems with a potentially higher frequency of strong melting events (e.g. case 3). Finally, the higher humidity in the Arctic Basin could also result in increased cloud cover which, unlike fog, would enhance surface warming by increasing downward long-wave radiation (Huang and others, Reference Huang, Ding, Dong, Xi and Baxter2021).
6. Conclusion
The apparent stability and resilience of ice patches at WHI to the recent warming, despite their limited volume, is the result of a feedback mechanism that operates on winter accumulation (i.e. greater accumulation after strong ablation, and conversely), and on summer microclimatic conditions dominated by recurring fog that reduces the heat input to ice patch systems. These topoclimatic conditions dampen the effects of climate variations on ice patch dynamics but do not eliminate them. The air temperature trends over the coming decades as well as the future frequency of fog episodes and cyclonic activity along the northern coast of Ellesmere Island will therefore be critical to the preservation of ice patches. Based on the regional climate projections, we speculate that most of the ice patches at WHI will lose their perennial status in the next decades, turning into semi-permanent ice patches or even snow patches, as has already happened further south in the High Arctic (Woo and Young, Reference Woo and Young2014). A cycle of ice patch destruction–restoration following short-term fluctuations in summer conditions (e.g. Arctic and North Atlantic Oscillations; Holland, Reference Holland2003) will create a higher interannual variability in the water supply to the polar desert ecosystem, with likely important consequences for the dynamics of hydrological and geomorphic processes, and the evolution of vegetation. Given that the climate of polar deserts encompasses a broad spectrum of summer temperature, humidity conditions and topography, it is likely that the dominant drivers of ice patch mass balance vary significantly among polar sub-regions and cannot be modeled by simplistic large-scale methods such as temperature index models.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2023.44.
Data availability
CT-scan images of IP1 as well as internal temperatures and time-lapse photographs are available for download at the public data repositories of Nordicana D, managed by the Centre for Northern Studies. The complete citations for the dataset are: Fortier D and Davesne G (2021a) Computed tomography (CT) scans of ice cores collected in two perennial ice patches at Ward Hunt Island, Nunavut, Canada, v. 1.0 (2018–2018). Nordicana D85 (doi: 10.5885/45696CE-6CC0B447C7F84E6A) Fortier D and Davesne G (2021b) Monitoring the snow and temperature regime of a perennial ice patch, Ward Hunt Island, Nunavut, Canada, v. 1.0 (2017-2019). Nordicana D91 (doi: 10.5885/45720CE-0F556C84D96948F7).
Acknowledgements
This research has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), including the Discovery Frontiers project Arctic Development and Adaptation to Permafrost in Transition (ADAPT); the Networks of Centres of Excellence program ArcticNet; the Canada Research Chair program; the Northern Scientific Training Program; the Canadian Foundation for Innovation: Canadian Northern Studies Trust; Centre d'études Nordiques (CEN); and Fond de Recherche du Québec-Nature et Technologie (FRQNT). F. D. was supported in part by the French Polar Institute (IPEV). Logistical support was provided by the Polar Continental Shelf Program (PCSP) and Parks Canada graciously granted us the use of their facility. We are also grateful to Derek Mueller (Water and Ice Research Laboratory, Carleton University) for providing comments on this study. Finally, we thank the editors and the two anonymous reviewers for providing valuable feedback on the early versions of the manuscript.
Author contributions
G. D. performed field operations with the assistance of F. D., data analysis and wrote the manuscript. D. F. and F. D. contributed to the sampling design, to the manuscript writing and obtained research funding. C. K. provided inputs for data analysis and contributed to the writing of the manuscript.
Appendix A
Ground penetrating radar (GPR)
For this study, GPR data were sampled in standard reflection mode with a step size of 1 m using a hand-held antenna. Because of the difference between the dielectric properties of the firn/ice and the underlying soil/bedrock, the detection of ice thickness is reliable and requires minimal signal processing. This work was implemented with the software EKKO project V2R3. The low-frequency noise of the signal was removed using the filter ‘dewow’ (Meulendyk and others, Reference Meulendyk, Moorman, Andrews and MacKay2012). The conversion of the electromagnetic signal propagation time to depth was made using a mean velocity (m ns−1) deduced from common midpoint surveys. A topographic correction was applied to the GPR data based on terrestrial laser scan (VX spatial station; Trimble ®) surveys and the ice thickness was extracted at 1 m intervals along each track. The ice thickness estimation was then validated with the ice coring made in 2017 (Davesne and others, Reference Davesne, Domine and Fortier2022).
Appendix B
Time-lapse photograph
From the raw set of time-lapse photographs, only one image per day was selected for analysis. Images captured in the afternoon were preferred because the ice patch is in the shadow of Walker Hill, avoiding light saturation issues. For zone B which focuses on the ice patch system IP1, we relied on a set of ten georeferenced images, including orthomosaics derived from uncrewed aerial vehicle surveys undertaken during the 2019 melt season and aerial photographs taken in summer 2017, to convert the snow/ice cover fraction derived from the time-lapse camera into a horizontal projected areal extent (A in m2). The georeferenced images were analyzed in ArcGIS to manually delineate the perimeter of IP1 and extract the horizontal projected areal extent. The relation between the fraction area and the corresponding value of areal extent is shown in Figure 14.
Relation between the fraction of area of the IP1 system derived from the time-lapse camera and the corresponding value of areal extent derived from detection on ArcGIS.
Appendix C
Terrestrial laser scanning
Measurement of the ice patch surface change was performed by repeated terrestrial laser scanning (TLS) surveys (Fischer and others, Reference Fischer, Huss, Kummert and Hoelzle2016). The TLS consists of a highly accurate method to obtain dense point clouds of a target surface. We used a VX spatial station (Trimble®, single 3D point accuracy 10 mm at ≤250 m). The deviation in the z-direction between various scans, estimated using unchanged reference surfaces such as rock outcrops, was lower than ±15 mm on average. The station remained at the same location through the monitoring periods to reduce the potential error linked to the instrumentation resettlement. The station position was linked to a geodesic landmark recorded using a Global Navigation Satellite System (R8 GNSS; Trimble®, precision x–y ± 8 mm and z ± 15 mm). Elevations obtained were corrected by the Canadian Centre for Remote Sensing and orthometric heights were used. The postprocessing consisted of cleaning the point clouds using the commercial software Trimble RealWork 7.1. They were subsequently used to create high-resolution 3D digital surface models using a triangulated irregular network method in ArcGIS (Esri, version 10.5.1).
Appendix D
Wind speed estimation
The wind data at IP1 were extrapolated from the time series of the SILA weather station recorded by a 10 m-high anemometer. The first step was to extrapolate the wind speed for a height of 2 m, corresponding to the height of the anemometer installed in 2019 at IP1. This extrapolation was based on the logarithmic wind law that has generally been used to model the vertical profile of wind speed over flat terrain (Oke, Reference Oke1987) given by Eqn (A1):
where v is the velocity at a height of 2 m above ground level (z), v ref is the wind speed measured at z ref of 10 m, z 0 is the roughness length. Since the terrain around the SILA weather station consists of a gravelly surface, we assumed a z 0 of 0.01 m. We then used the relationship between the average hourly wind speed (m s–1) recorded by the automatic weather station (AWS) installed on IP1 during the study period 2019 with the data extrapolated for a 2 m-height at SILA station. Figure 15 shows the good relationship between both datasets (R 2 = 0.71). The equation associated with the linear regression line (WSIP1 = 0.6087 × WSSILA + 0.4832) was used to estimate the wind speed at 2 m on IP1 in 2017.
Correlation between the hourly average wind speed derived from SILA measurements adjusted using Eqn (A1) for a 2 m-height level at the SILA station from time series provided by the anemometer at 10 m and by the AWS on IP1 during the study period 2019.
Appendix E
Surface energy balance
The turbulent heat fluxes Q h and Q e (W m–2) were calculated based on the bulk aerodynamic approach, including stability correction, as formulated by Price and Dunne (Reference Price and Dunne1976). This method is based on differences in wind speed, potential temperature and specific humidity between the measurement level and the surface. Equations (A2) and (A3) describe the calculation method of Q h and Q e, respectively, using measurements of wind speed, temperature and humidity at the 0.5 m level (Young and Lewkowicz, Reference Young and Lewkowicz1990; Boike and others, Reference Boike, Roth and Ippisch2003a):
where ρ a is the air density (1.27 kg m–3), C a is the specific heat of air (1.005 × 103 J kg–1 K–1), L is the latent heat of vaporization or sublimation (being 2.48 × 106 and 2.83 × 106 J kg–1, respectively), AP is the atmospheric pressure (hPa) and D is the bulk exchange coefficient. T a and T ss are the respective air and snow surface temperatures (°C), e a is the vapor pressure (hPa) at T a, and e ss is the saturation vapor pressure at the surface of a melting snowpack (6.11 hPa). The snow surface temperature T ss has not been measured but it is assumed to be 0°C when the T a is positive and to be equal to T a when T a is negative. This assumption implies a neutral gradient in sub-zero conditions, leading to a possible underestimation of sensible heat flux. However, given that freezing conditions were very marginal during the study period, we consider that this simplification had negligible consequences on the SEB.
For neutral conditions, D is equal to the momentum D m:
where uz is the wind speed (m s–1) at height z, k is von Karman's constant (0.4) and z 0 is the roughness length. For the snow, z 0 values of 0.001–0.005 m are commonly assumed (Liston and Hall, Reference Liston and Hall1995; Bash and others, Reference Bash and Moorman2020). Here, we calculated the turbulent fluxes with an intermediate value of z 0 = 0.0025. Given the relative importance of the surface roughness, a sensitivity analysis was performed using values of 0.001 and 0.005. It resulted in <10% change to the turbulent fluxes on average at the scale of the study period, considered here as an acceptable error.
Under conditions other than neutral, D is adjusted with the bulk Richarson number (Ri) expressed as:
where g is the gravitational constant.
Under stable conditions Ri < 0:
Under unstable conditions Ri > 0:
The heat flux by conduction through the snowpack (Q c in W m–2) was calculated from the temperature gradient between the temperature of the snow surface (T ss in °C) and the snow temperature at −0.25 m (T snow in °C) and using the average snow thermal conductivity measured during snowpit excavation on 6 June 2019 (Fig. S3). Since the upper thermistor was at 40 cm below the snow surface, we linearly interpolated the snow temperature at 25 cm below the surface from the temperature profile.
Q c was expressed as:
where ΔT is the temperature gradient between the snow surface and the snowpack bottom and R is the thermal insulation effect of the snow. R (m2 K W–1) can be calculated from the thickness h (m) and thermal conductivity k (W m–1 K–1) of the various snow layers (i) that comprise the snowpack as given by Eqn (A9) (Domine and others, Reference Domine, Barrere and Morin2016):
Open access
