Modelling the development and decay of cryoconite holes in northwestern Greenland
Cryoconite holes (CHs) are water-filled cylindrical holes with cryoconite (dark-coloured sediment) deposited at their bottoms, forming on ablating ice surfaces of glaciers and ice sheets worldwide. Because the collapse of CHs may disperse cryoconite on the ice surface, thereby decreasing the ice surface albedo, accurate simulation of the temporal changes in CH depth is essential for understanding ice surface melt. We established a novel model that simulates the temporal changes in CH depth using heat budgets calculated independently at the ice surface and CH bottom based on hole-shaped geometry. We evaluated the model with in situ observations of the CH depths on the Qaanaaq ice cap in northwestern Greenland during the 2012, 2014, and 2017 melt seasons. The model reproduced the observed depth changes and timing of CH collapse well. Although earlier models have shown that CH depth tends to be deeper when downward shortwave radiation is intense, our sensitivity tests suggest that deeper CH tends to form when the diffuse component of downward shortwave radiation is dominant, whereas CHs tend to be shallower when the direct component is dominant. In addition, the total heat flux to the CH bottom is dominated by shortwave radiation transmitted through ice rather than that directly from the CH mouths when the CH is deeper than 0.01 m. Because the shortwave radiation transmitted through ice can reach the CH bottom regardless of CH diameter, CH depth is unlikely to be correlated with CH diameter. The relationship is consistent with previous observational studies. Furthermore, the simulations highlighted that the difference in albedo between ice surface and CH bottom was a key factor for reproducing the timing of CH collapse. It implies that lower ice surface albedo could induce CH collapse and thus cause further lowering of the albedo. Heat component analysis suggests that CH depth is governed by the balance between the intensity of the diffuse component of downward shortwave radiation and the turbulent heat transfer. Therefore, these meteorological conditions may be important factors contributing to the recent surface darkening of the Greenland ice sheet and other glaciers via the redistribution of CHs.
Cryoconite holes (CHs), cylindrical water-filled holes formed on ablating ice surfaces of glaciers and ice sheets worldwide, influence the ablation process of glacial ice and function as the habitats of glacial microbes (e.g. Hodson et al., 2010a; Cook et al., 2016; Uetake et al., 2016; Zawierucha et al., 2021). In CHs, dark-coloured sediment, referred to as cryoconite, is deposited at the bottom; it consists mostly of spherical aggregates (cryoconite granules) formed by the entanglement of filamentous cyanobacteria with mineral and organic particles (Hodson et al., 2010b; Langford et al., 2010; Takeuchi et al., 2010, 2014; Wientjes et al., 2011; Uetake et al., 2019). Because such dark-coloured sediment on bare ice absorbs solar radiation more effectively than the surrounding clear ice, CH formation occurs often.
The development and decay of CHs depend on the meteorological conditions near the ice surface. CHs are likely to form well under sunny conditions (McIntyre, 1984; Fountain et al., 2008; Jepsen et al., 2010; Cook et al., 2016; Takeuchi et al., 2018). Although many studies have investigated which meteorological conditions are favourable for the deepening (development) of CHs, there is limited information on their shallowing (decay) dynamics. Takeuchi et al. (2018) reported that CHs collapse under cloudy and windy conditions in summer in northwestern Greenland. They suggested that CHs collapse when the ice surface surrounding the CHs experiences enhanced melt caused by turbulent heat flux exchange. CH collapse events have also been reported in Svalbard and southwestern Greenland (Hodson et al., 2007, 2008; Stibal et al., 2008; Irvine-Fynn et al., 2011; Chandler et al., 2015), while it has been indicated that a higher melt rate at the ice surface induces CH collapse (Hodson et al., 2007).
The development and decay of CHs modulate the ice surface albedo. On the one hand, the development of CHs increases the area-averaged ice surface albedo compared to the case where the surface is uniformly covered with cryoconite granules because cryoconite granules in holes are shielded from sunlight by the hole walls (Bøggild et al., 2010; Takeuchi et al., 2014; Chandler et al., 2015). On the other hand, the collapse of CHs reduces the area-averaged surface albedo by approximately 0.1–0.2 due to the redistributed cryoconite on the ice surface (Irvine-Fynn et al., 2011; Takeuchi et al., 2018). Topologically heterogeneous ice surfaces of the Greenland ice sheet (GrIS) can be classified into four types: clean bare ice surfaces, dirty bare ice surfaces, surfaces with CHs, and meltwater streams (Irvine-Fynn and Edwards, 2014; Chandler et al., 2015; Holland et al., 2019; Tedstone et al., 2020). Bare ice surfaces and surfaces with CHs are estimated to account for approximately 80 %–90 % and 5 %–15 %, respectively, of the southwestern GrIS ablation area (Chandler et al., 2015). Therefore, both the development and decay of CHs are important for the surface albedo dynamics of the GrIS.
Numerical models have been used for reproducing the temporal changes in CH depth; however, only few models consider the decay of CHs. Gribbon (1979) proposed a conceptual model based on the Beer–Lambert law to simulate the development of CHs. By using solar radiation and CH depth, this model calculates the shortwave radiation reaching the cryoconite at the CH bottoms using the absorption coefficient of ice. Although several studies have proposed numerical models following Gribbon (Podgorny and Grenfell, 1996; Jepsen et al., 2010; Cook, 2012), these models focus only on the development of CHs. Therefore, to reproduce CH decay processes, such as CH collapse events, a model that simulates melt not only at the CH bottom but also at the ice surface surrounding the CHs is necessary.
In this study, we established a numerical model to simulate the temporal changes in CH depth, accounting for both the development and decay of CHs, based on a concept different from those in previous studies (e.g. Gribbon, 1979). We evaluated the model using observed temporal changes in CH depth from the Qaanaaq ice cap in northwestern Greenland. We further performed numerical model sensitivity analyses to reveal how the CH depth is controlled by the meteorological variables, hole geometry, albedos of the ice surface surrounding the CHs and the cryoconite therein, and optical parameters of ice. Finally, we assessed the meteorological factors controlling the CH depth changes based on heat component analysis.
The cryoconite hole model (hereafter referred to as CryHo) developed in this study calculates the heat balances and resultant melt rates at the ice surface and CH bottom independently (Fig. 1). All variables and constants described in the following equations are listed in Table 1.
2.1 Heat balance at the ice surface
The heat balance at the ice surface (Qi, W m−2) is described as follows:
where αi is the ice surface albedo (dimensionless); RS s the downward shortwave radiation (W m−2); RLni is the net longwave radiation (W m−2); HSi is the sensible heat flux (W m−2),; HLi is the latent heat flux (W m−2); ε is the emissivity of the snow and ice surface, which is assumed to be 1.0 (dimensionless); RL is the downward longwave radiation (W m−2); σ is the Stefan–Boltzmann constant (5.67 × 10−8 ); and Ti is the surface temperature (K). Subscript i used in the variables refers to the ice surface. All the downward components have positive signs. Note that HLi is restricted to the latent heat of evaporation; it does not include the latent heat of melting. Heat conduction from or toward the glacier ice and heat supplied by rainwater are assumed to be negligible. The turbulent heat fluxes (i.e. HSi and HLi) were estimated using the following bulk formulations:
where cp is the specific heat of air (1006 ), ρa is the air density (kg m−3), C is the bulk coefficient for the snow and ice surface (0.0025, Kondo, 1994), U is the wind speed (m s−1), Ta is the air temperature (K), lE is the latent heat of water evaporation (2.50 × 106 J kg−1), hr is the relative humidity (dimensionless), and q(T) is the saturated specific humidity (kg kg−1). The excess energy (Qi > 0) is used to melt the ice at the surface, once the surface temperature reaches the melting point (273.15 K), and its hourly amount is obtained as follows:
where QMi is the excess heat energy for ice melt (W m−2), Mi is the melting ice thickness (m h−1), th is the length of an hour in seconds (3600 s), lM is the latent heat of ice melt (3.33 × 105 J kg−1), and ρi is the ice density assumed equal to 900 kg m−3.
2.2 Heat balance at the cryoconite hole bottom
The heat balance and melt rate at the CH bottom (denoted by subscript c) are defined as follows:
where Qc is the heat balance at the bottom (W m−2), αc is the albedo (dimensionless), RSdc and RSfc are the direct and diffuse components of shortwave radiation (W m−2), RStdc and RStfc are the direct and diffuse components of shortwave radiation transmitted through ice (W m−2), RLnc is the net longwave radiation (W m−2), QMc is the excess heat energy for ice melt (W m−2), and Mc is the melting ice thickness (m h−1). Because we assume that the CH is partially filled with water, the turbulent heat fluxes at the bottom of the CH (i.e. HSc and HLc) are assumed to be 0 in Eq. (6). The change in CH depth (Dt, m) at a given time (t) is defined as follows:
where Dt−1 is the CH depth at one time step before (m). If the melt rate at the CH bottom is greater than that at the ice surface (Mc>Mi), the CH depth deepens (and vice versa). The initial depth D0 at t=0 in CryHo is a prescribed constant initial condition. For simplicity, we assume the heat balance at the centre of the CH bottom can represent that of the entire bottom.
2.3 Geometry effect on the radiation components
In earlier CH models (Gribbon, 1979; Jepsen et al., 2010; Cook, 2012), downward shortwave radiation is reduced exponentially with increasing CH depth. The novel aspect of our model is that it considers not only the CH depth but also the CH shape geometry. The direct and diffuse components of downward shortwave radiation and downward longwave radiation are also considered independently. Regarding downward shortwave radiation, we first separate the observed downward shortwave radiation into direct and diffuse components (RSd and RSf, W m−2) as follows:
where rdif is the diffuse ratio of downward shortwave radiation (dimensionless; 0–1), which is calculated following Goudriaan (1977):
where rze is the ratio based on the zenith angle (dimensionless; 0–1) and rcld is the ratio based on cloudiness (dimensionless; 0–1). Calculation of rze was performed following Goudriaan (1977), and calculation of rcld was performed following Broeke et al. (2004) and Niwano et al. (2015):
where θz is the solar zenith angle (radian).
The CH geometry is defined as the zenith angle of the edge from the centre of the CH bottom (θc, rad) as follows:
where ϕ is the diameter of CH (m).
Compared with the solar zenith angle, the direct component of shortwave radiation from the sky looking up from the CH bottom (RSdc; W m−2) is calculated as follows:
Meanwhile, the diffuse component of shortwave radiation from the sky looking up from the CH bottom (RSfc, W m−2) is obtained as follows:
Shortwave radiation components transmitted through ice are described in the following section.
Longwave radiation from the sky looking up from the CH bottom (RLc; W m−2) can be obtained in a manner similar to RSfc as follows:
Longwave radiation from the CH wall (RLw; W m−2) is subsequently described as follows:
Here, we assume that the CH wall or bottom temperatures (Tc; K) are equal to the melting point. Because the CryHo does not calculate water level in the CH, the longwave radiation emitted from the CH wall is calculated using snow and ice surface emissivity, which is similar to the water emissivity in Eq. (18). Regarding longwave radiation emitted from the CH bottom, the net longwave radiation (RLnc; W m−2) can be summarized as follows:
2.4 Shortwave radiation transmitted through ice
Depending on whether direct solar radiation illuminates the CH bottom (θz≤θc) or not (θz>θc), the direct and diffuse components of shortwave radiation transmitted through the ice (RStdc and RStfc, W m−2) are described as follows:
where κd and κf are the broadband asymptotic volume flux extinction coefficients (hereafter broadband flux extinction coefficients) of ice (m−1) for the direct and diffuse components, respectively. For the direct component, the path length of the transmitted light is considered, whereas for the diffuse component, it is assumed to be the CH depth.
The broadband flux extinction coefficient for the diffuse component (κf) can be calculated from the spectral flux extinction coefficients. However, this is not a simple integration because the spectral distribution of shortwave radiation transmitted through the ice varies depending on both the solar illumination spectral distribution and the path length in the ice. Thus, we first calculate the insolation spectral distribution in the ice from both the spectral flux extinction coefficients experimentally determined by Cooper et al. (2021) for ice on the GrIS (the term “irradiance attenuation coefficient” was used in their study) and the spectral solar radiation at the ice surface assumed for clear and cloudy skies, computed using a radiative transfer model (Aoki et al., 1999; 2000).
2.4.1 Theoretical basis for calculating bare-ice broadband flux extinction coefficient
Section 2.4.1 presents an expression for the broadband asymptotic volume flux extinction coefficient, κ (m−1), from the spectral asymptotic volume flux extinction coefficient ke(λ) (m−1) of bare ice, where λ is the wavelength (µm). In this study, we simply refer to them as the “broadband (or spectral) flux extinction coefficient.” When the bare ice surface is illuminated by the spectral solar radiation Fs(λ), the attenuated solar radiation Fb(λ) at ice thickness x (m) in bare ice is expressed by ke(λ) as follows:
Using Eq. (22), the spectral flux transmittance (T(λ,x); dimensionless) as a function of x is described as follows:
The spectrally integrated flux transmittance is calculated from Fs(λ) and T(λ,x) as follows:
Using Eq. (24), the broadband flux extinction coefficient (, m−1) is described as follows:
2.4.2 Parameterization of broadband flux extinction coefficient
Using Eq. (25), we calculate the for bare ice in Greenland from the following two datasets for ke(λ) and Fs(λ). The first dataset ke(λ) is employed based on measurements by Cooper et al. (2021) for the bare ice “Layer A” and “Layer B”, which are defined as the ice layers at 12–77 and 53–124 cm depth, respectively. The second dataset Fs(λ) is theoretically simulated at λ = 0.24.0 µm with a spectral resolution of 0.025 µm using a radiative transfer model for the atmosphere–snow-and-ice system (Aoki et al., 1999, 2000), assuming clear- and cloudy-sky conditions in the subarctic model atmosphere over the bare ice surface. As the spectral range of ke(λ) reported by Cooper et al. (2021) is limited around the visible wavelengths (λ = 0.35–0.75 µm), with their Table A1 data (λ = 0.35–0.60 µm) and the calculated values from their transmittance data (λ = 0.60–0.75 µm) being employed in this study, we assume ke(λ) in the other spectral regions λ = 0.20.35 µm and λ = 0.754.0 µm by scaling the pure-ice spectral absorption coefficient ka(λ) (Warren and Brandt, 2008) to connect to the values of Cooper et al. (2021), as shown in Fig. 2. It is inferred that the difference between the scaled ke(λ) (pink curve) and pure-ice ka(λ) (blue curve and red x symbols) is attributed to air bubbles in bare ice for λ ≥ 0.55 µm and impurities and air bubbles for λ < 0.55 µm from the discussion by Cooper et al. (2021). There are large uncertainties in ka(λ) in λ = 0.35–0.75 µm, as shown by the blue curve (Warren and Brandt, 2008) and the red x symbols (Cooper et al., 2021) in Fig. 2. Regarding the scaled ke(λ) validity in other spectral regions, we did not have verifiable data. In fact, the contribution to from ke(λ) in the other spectral regions is very low. The simulated Fs(λ) for clear and cloudy skies for the solar zenith angle θz = 63.1∘ is shown in Fig. S1 in the Supplement.
Figure S2 demonstrates as a function of θz for x = 0.1 m calculated from ke(λ) and Fs(λ) described above using Eq. (25). There is almost no difference in when we use ke(λ) for layers A and B, as observed by Cooper et al. (2021). We also find that the θz dependence of is small in the range of θz from 50–86∘, and thus we determine the θz-independent values by averaging it over this range of θz. By plotting the θz-independent values of as a function of x for clear and cloudy skies (Fig. 3), the values rapidly decrease in the shallower ice layer and change to a gentle slope in the deeper ice layer when x increases. This is because the spectral distribution of Fb(λ) in bare ice changes from a wider distribution to a narrower distribution centred around the minimum ke(λ) wavelength, where the light absorption by bare ice is weakest as x increases. We can then obtain the approximation equations by fitting the calculated values of with a power regression equation, as shown in Table S1 in the Supplement. In this study, we use the terms “κclr” and “κcld” as the values for clear and cloudy skies, respectively. Finally, we determine the following approximation equations (Eqs. 26 and 27) of the broadband flux extinction coefficients as a function of ice thickness for clear (κclr; m−1) and cloudy (κcld; m−1) skies (Eq. 25 and Fig. 3).
where κclr and κcld are calculated using the Layer A values shown in Table S1, as we focus on the vertical change in cryoconite holes below 0.20 m. The broadband flux extinction coefficient for the diffuse component (κf) is determined in a similar manner for the estimation of RSd and RSf as follows:
As there are no available data for estimating the extinction coefficient for the direct component (κd), we parameterize it by introducing a coefficient (rd) as follows:
The coefficient rd was assumed to be . The value 1.66 is used as an approximation value of the effective optical path length when the transmittance of the diffuse component of shortwave radiation is obtained from that of the direct component of shortwave radiation (Liou, 1980). Because we first estimate the diffuse component coefficient, κd was obtained by dividing κf by 1.66.
3.1 Site location
To evaluate the model, we investigated changes in CH depth at the Qaanaaq ice cap in northwestern Greenland (Fig. 4) during the 2012, 2014, and 2017 summer seasons. The ice cap covers an area of 286 km2 and has an elevation range of 0–1200 (Sugiyama et al., 2014). We selected five study sites at different elevations (sites 1 to 5) on the Qaanaaq Glacier, which is an outlet glacier of the ice cap that is easily accessible from Qaanaaq village. Sites 1 to 5 are located in the middle of the glacier at elevations of 247, 441, 672, 772, and 944 , respectively. The equilibrium line elevation of the glacier ranged between 862 and 1,001 in early August 2012, 2014, and 2017 (Tsutaki et al., 2017); hence, the ice surfaces at sites 1 to 5 were exposed during the ablation period from July to August.
3.2 Observations of meteorological variables and surface reflectance
The meteorological data that were used as input in CryHo were collected with an automatic weather station (AWS), which was established at the SIGMA-B site (77.518∘ N, 69.0619∘ W; 944 ), i.e. at Site 5 in this study, on 19 July 2012 (Aoki et al., 2014; Nishimura et al., 2021). Hourly air temperature (Ta; K), relative humidity (hr; dimensionless), downward shortwave radiation (RS; W m−2), downward longwave radiation (RL; W m−2), upward longwave radiation (RLui; W m−2), wind speed (U; m s−1), and air pressure (Pa; hPa) were collected for the model simulations described in Sect. 4. The seasonal changes in the measured air and surface temperatures and the calculated daily surface energy balance during the three studied summer seasons are shown in Fig. 5. The air and surface temperatures support the assumption that all study sites below the SIGMA-B site (i.e. Site 5) were mostly in ablation conditions during the three studied summer seasons. The heat balance was similar across the study years.
To determine the ice surface albedo used as input in CryHo, we measured the spectral reflectance from the visible to near-infrared band (0.350–1.050 µm) on the ice surface using a spectrometer (MS-720, Eiko Seiki Co., Japan) at the study sites during the 2014 and 2017 summer seasons (Table 2). For details on the observation method, we refer to Takeuchi et al. (2015). Although the mean reflectance differs optically from broadband albedo (0.300–3.000 µm), the reflectance measured at Site 3 in 2014 agrees well with the broadband albedo at the same site measured in 2012 (0.4 versus approximately 0.4; Takeuchi et al., 2018). Therefore, we assume that the reflectance shown in Table 2 can be used as a proxy for broadband albedo at the study sites for numerical simulations.
3.3 Monitoring cryoconite hole depth
To collect in situ data that can be used for the evaluation of CryHo, the temporal changes in CH depth were observed using the monitoring device, which consisted of two time-lapse cameras as well as a plastic stick positioned in the centre of a CH, supported by metal angles buried to approximately 1.5 m depth in the ice body (Fig. 6a). The temporal changes in CH depth were derived from variations in the heights of the plastic stick and white metal angle (Fig. 6b–d). To measure the variations, we marked both the stick and angle with clear and black colours at 5 cm intervals. The stick and angle were used to estimate the melt rates at the CH bottoms and ice surface surrounding the CHs, respectively. The variations were measured from hourly images with a resolution of 4608 pixel by 2592 pixel captured by a time-lapse camera (PENTAX WG-10, RICOH Co., Japan) installed at 3 m from the marked stick and angle. The time-lapse camera was placed 1 m above the ice surface and orthogonally oriented to the ice surface. To maintain the camera at the same height from the ice surface, we adjusted the height regularly (i.e. every 5–10 d). To capture the CH collapse event and check whether the plastic stick remained positioned in the CH, we also monitored the CHs using hourly images with a resolution of 1280 pixel by 1022 pixel taken by another time-lapse camera (GardenWatch Cam, Brinno Co., Taiwan, Fig. 6b). The location of the device was selected to have a suite of characteristics broadly identical to those of the monitoring site in a previous study (Takeuchi et al., 2018) in terms of cryoconite loading, ice structure, and ice surface slope of ∼ 5∘. Monitoring was conducted from day-of-the-year (DOY) 194 to 214 in 2014.
To estimate the CH depths, the images were processed using ImageJ (version 1.48, National Institutes of Health, USA; Schneider et al., 2012). The CH depths were measured as relative variations using the length of marks on the stick for the CH bottom (Fig. 6c) and on the angle for the ice surface (Fig. 6d). One pixel corresponds to approximately 1 mm. Because the monitored CH collapsed between DOY 202 and 204, the sticks and angles were re-installed into another CH on DOY 205 (Fig. 6e).
Camera-based observations were conducted at Site 2, near the site where CH collapse occurred in 2012 (Takeuchi et al., 2018). We also conducted manual measurements of the depth and diameter of CHs two to four times during the 2014 and 2017 summer seasons. In addition, we used the dates of the CH collapse events reported by Takeuchi et al. (2018) for model evaluation. CH depths and diameters were manually measured from 20 to 50 randomly selected CHs at each site. Averages and standard deviations of the observed depths and diameters are summarized in Table 2.
To evaluate CryHo, numerical simulations of the CH depth were conducted at several sites of the Qaanaaq ice cap for the three summer seasons (Site−exp). The number of study sites in 2012, 2014, and 2017 were 1 (Site 3), 4 (sites 1 to 4), and 3 (sites 1 to 3), respectively. The depths and diameters observed on the first observational date were used to set the model constants at each site (D0 and ϕ in Table 2). The meteorological data measured at Site 5 were used as input for the simulations. At each site, Ta was corrected based on the air temperature measured at Site 5 and a temperature lapse rate of 7.80 × 10−3 K m−1, which had been previously estimated from field observations on the same glacier in 2012 (Sugiyama et al., 2014). The other meteorological data at each site were assumed to be the same as those measured at Site 5. The observed surface temperature at Site 5 was almost the melting point during the three observational periods; thus, it was safely also used as input for the other sites, since they were located at lower elevations than that of Site 5 (Fig. 5a). In 2014, αi at sites 1, 2, 3, and 4 was assumed to be 0.68, 0.57, 0.40, and 0.41, respectively (Table 2), as mentioned in Sect. 3.2. Similarly, in 2017, αi at sites 1, 2, and 3 was assumed to be 0.56, 0.42, and 0.24, respectively (Table 2). αi in 2012 for Site 3 was assumed to be 0.4, as observed by Takeuchi et al. (2018). We assumed the albedo at the CH bottom (αc) to be 0.1 based on assumptions and observations in previous studies (e.g. Takeuchi et al., 2001; Jepsen et al., 2010). The calculation periods for the simulations in 2014 and 2017 were taken from the first date when CH depths were measured in June or July. In the 2012 simulation, the period was constrained by both CH and AWS observations.
a Change from the control experiment. b Shortwave radiation assumed to be only direct or diffuse component.
We conducted sensitivity tests to assess the sensitivity of the CH depth to input data and model constants, such as air temperature (Ta-exp), radiation components (RS-exp), initial depth (D0-exp), hole diameter (ϕ-exp), albedo at the ice surface (αi-exp), albedo at the CH bottom (αc-exp), extinction coefficients of direct (κd-exp) and diffuse (κf-exp) radiation, solar zenith angle (θz-exp), and zenith angle of the edge from the centre of the CH bottom (θc-exp) (Table 3). Site−exp, i.e. Site 2 in 2014, was used as the control experiment for the sensitivity tests (Ctl−exp). The ranges of the changing parameters, which are summarized in Table 3, were determined based on field measurements (Table 2). The extinction coefficients for κd-exp and κf-exp were obtained from multiplying the original values by factors of 0.25–4.00. The factor range was assumed by referring to the difference between the spectral flux extinction coefficient and absorption coefficient calculated from the imaginary refractive index of pure ice (Fig. 2). In RS-exp, we assumed rdif of Eqs. (9) and (10) to be 0 and 1 in Sd and Sf cases shown in Table 3, respectively. Because we did not consider the light refraction at the air–water surface in the CH, we evaluated the refraction effect via sensitivity analysis (θz-exp) using various incident angles (θz). Although the heat balance at the CH bottom should vary across the width of the hole, particularly between the northern and southern edges, we assumed heat balance only at the centre of the CH bottom for the model simulation. The effect of different edges is evaluated by changing θc (θc-exp). In θz-exp and θc-exp, θz and θc calculated in the model were replaced with the values shown in Table 3, respectively.
To quantify the four components of shortwave radiation reaching the CH bottom (RSdc, RSfc, RStdc, and RStfc) for different CH geometries (i.e. depths and diameters), we conducted two sensitivity tests (RScD-exp and RScϕ-exp, respectively) at 13:00 and 01:00 LT (local time) on DOY 172 (i.e. the time of meridian transit and midnight at the summer solstice) in 2014. The CH depth and diameter ranged from 0.01 to 0.35 m in RScD-exp and from 0.01 to 0.15 m in RScϕ-exp, respectively. The CH diameter and depth were assumed to be 0.05 m in RScD-exp and 0.01 m in RScϕ-exp, respectively. The other model constants were the same as those in Ctl−exp.
To discuss the meteorological factors controlling the CH dynamics, besides the sensitivity tests described above, we also investigated the relationships between the CH deepening or shallowing rates against the input meteorological conditions, and the simulated heat properties at the ice surface at the CH bottom. These relationships were obtained from the input and output data in Ctl−exp. The deepening (shallowing) rate (m h−1) was calculated when a CH deepened (shallowed) compared to its state 1 h earlier.
5.1 Evaluation of the cryoconite hole depth
The CH depth simulations at the different elevation sites in the 3 years (Site−exp) showed that the simulated temporal changes in CH depth, including CH collapse, agreed well with those observed at all sites for the 3 years (Figs. 7 and 8). The CH depth observed by both the camera and manual measurements at Site 2 gradually shallowed from DOY 186 to 202, developed until DOY 212, and shallowed again until DOY 215 (Fig. 7b). The modelled and camera-based CH depths at this site were in good agreement with each other (R2 = 0.78; root-mean-square error of 0.04 m). In addition, the model also simulated the CH collapse event between DOY 202 and 204, during which the camera-based CH depth was not recorded (line gap in Fig. 7b). The model simulations revealed the good performance of the model also at the other sites (Fig. 7a, c, and d). Although no CH depth measurement was conducted in 2012, the model captured the CH collapses observed at Site 3 on DOY 206 and 208 (Fig. 8a, dashed grey lines). The manual measurements in 2017 showed similar temporal changes in CH depth to those in 2014 (Fig. 8b–d), which were also well represented by the model, especially the timing of drastic CH decay around DOY 210 at all sites. Because CryHo diagnoses the CH depth based on the balance between melt at the ice surface (Mi) and at the CH bottom (Mc), the model can capture the CH decay dynamics, which have not been simulated by previous numerical models considering only Mc (Gribbon, 1979; Jepsen et al., 2010; Cook, 2012). This indicates that the concept of CryHo, in which both heat balances (i.e. Qi and Qc) control the temporal change in the development and decay of CHs, is essentially correct. In the next section, we discuss the meteorological and/or ice physical factors that are important for simulating temporal changes in detail.
5.2 Sensitivity of the simulated hole depth to the model parameters
The sensitivity experiment regarding the air temperature (Ta-exp) suggests that CHs tend to be deeper and more stable under cooler conditions (and vice versa) (Fig. 9a). This is likely because the turbulent heat fluxes on the ice surface decrease under cooler conditions; thus, net radiation becomes the main component for ice melting, which is constrained by the albedo at the ice surface and CH bottom. Indeed, field observations indicated that the turbulent heat fluxes during warmer conditions (DOY 198–205) were increased compared to the case during cooler conditions (before DOY 198) at Site 2 in 2014 (Fig. 5b). Under such cooler conditions, the melt rate at the ice surface surrounding the CHs was lower, resulting in deeper and therefore more stable CHs. Indeed, Gribbon (1979) reported that the CH depth increased with the elevation of the glacier in western Greenland. Furthermore, higher melt rates (and hence warmer conditions) at the ice surface are likely to induce CH decay, as suggested by Hodson et al. (2007). Therefore, global warming may cause the collapse of CHs earlier in the melt season compared to the present case, thereby inducing darkening of the ice surface by cryoconite granule redistribution in the future. Nevertheless, a snow–ice–atmosphere coupling simulation using a climate model coupled with CryHo is necessary to evaluate the contribution of CH dynamics to surface darkening under climate change.
The sensitivity experiment regarding the shortwave radiation components (RS-exp) suggests that CHs tend to decay and develop when the direct and diffuse components, respectively, are dominant (Fig. 9b). We then compared the contributions of each radiation component reaching the CH bottom during the experimental period (Fig. 10, red and blue lines). The CH bottom is more accessible by the diffuse component (RSfc+RStfc) rather than the direct component (RSdc+RStdc), except for shallowing depth case. In the model, the direct component of shortwave radiation can reach the CH bottom only when the solar zenith angle θz is smaller than the CH edge θc (Fig. 1 and Eq. 15), and it is transmitted through the ice in the opposite case. Because θz and θc ranged from 56 to 85∘ and 8 to 90∘ during the simulation period at the studied glacier, respectively, the direct component reaching the CH bottom from the hole mouth was very limited. Figure 11a shows that there is no significant difference between the direct and diffuse components reaching the CH bottom even at the time of meridian transit in the summer solstice, suggesting that the diffuse component reaches the CH bottom more frequently than the direct component in the studied glacier. To investigate heat flux to the CH bottom by shortwave radiation from the hole mouth or through the ice, we additionally compared the contributions of each radiation component reaching the CH bottom (Fig. 10, grey lines). The figure shows that the radiation components transmitted through ice to the CH bottom (RStfc+RStdc) were greater than the radiation components reaching the CH bottom from the hole mouth (RSfc+RSdc) when CH developed from DOY 208 to 213, meaning that radiation components transmitted through ice are also important for the heat balance at the CH bottom (i.e. Qc). Further discussion regarding shortwave radiation transmitted throughout the ice is described later.
The sensitivity experiments regarding the CH geometry (D0-exp and ϕ-exp) suggest that differences in the geometrical parameters have little influence on the CH depth. The simulated CH depth starting from different initial depths converged within approximately 2 weeks (Fig. 9c). This is probably because the CH bottom is relatively accessible by shortwave radiation in the case of shallower depths (and vice versa) (Fig. 11), which is quantified by additional experiments (RScD-exp and RScϕ-exp). The diameter shows no significant effect on the CH depth (Fig. 9d), although it strongly affects whether shortwave radiation reaches the CH bottom. This is likely because transmitted shortwave radiation reaches the CH bottom irrespective of the CH diameter. Figure 11b suggests no significant difference in the total shortwave radiation reaching the CH bottom among the different diameters, thereby supporting the ϕ-exp result. Indeed, no significant correlation between the CH depth and diameter has been found (Gribbon, 1979; Cook, 2012). Because an increase in CH diameter allows more direct shortwave radiation to reach the CH bottom, positive feedback of the CH development is possible. However, such feedback is unlikely to occur because Fig. 11b and d suggest that most of the diffuse component of the transmitted shortwave radiation reaches the CH bottom at 0.01 m depth in CHs over 0.03 m in diameter. Furthermore, the observed CH depths and diameters significantly varied among the sites and years (Table 2), suggesting that CH depth is mainly controlled by factors other than the CH diameter.
The sensitivity experiments regarding the albedos of the ice surface and CH bottom (αi-exp and αc-exp) suggest that the difference between ice surface (αi) and CH bottom albedo (αc) is an important factor for reproducing the CH dynamics, especially the timing of CH collapse. The CH depth increases with an increase in αi owing to decreasing Mi, whereas it decreases with an increase in αc owing to decreasing Mc (Fig. 9e and f). Notably, the sensitivity of the CH depth to αi was greater than that to αc. This is probably because shortwave radiation at the ice surface was greater than that at the CH bottom. In addition, αi-exp suggests that a 0.1 decrease in αi induces CH collapse 1 d earlier. Although αc is known to be a key parameter for simulating the CH deepening rate (Podgorny and Grenfell, 1996), there is little information and discussion regarding αi. Because we used constant values of αc and αi in αi-exp and αc-exp, respectively, the CH depth sensitivity to the albedo is equivalent to the sensitivity difference between αi and αc. Therefore, our results first highlight the importance of the difference between αi and αc in simulating vertical CH variations in which CH tends to develop with greater albedo difference and vice versa. Furthermore, the collapse of CHs is likely to reduce the ice surface albedo, owing to the redistribution of cryoconite on the ice surface (Irvine-Fynn et al., 2011; Takeuchi et al., 2018). Considering studies that reported a darkening trend of the GrIS surface over the last 20 years (e.g. Shimada et al., 2016; Tedstone et al., 2017), a positive albedo feedback, during which a decrease in ice surface albedo causes frequent CH collapse events, may occur. While our assumed reflectance was based on the αi observations in 2014, the ice surface albedo should change temporally. The uncertainty in the model simulations, such as that regarding the results around DOY 204 at Site 3 in 2017 (Fig. 8d), may be attributed to temporal changes in the ice surface albedo. To simulate the development and decay of CHs more accurately, temporal changes in the ice surface albedo, including the effects of CH collapse and algal blooms, should be incorporated into CryHo.
The sensitivity experiments regarding the broadband extinction coefficients of shortwave radiation transmitted throughout ice for the direct and diffuse components (κd-exp and κf-exp) suggest that κf is more effective at the CH depth than κd. Both experiments showed that the CH depth with a higher coefficient was shallower, owing to a decrease in Mc (Fig. 9g and h); however, there was a significant difference in the sensitivity to the coefficients. One reason is probably that κd is lower than κf, as in Eq. (29). Figure 10 showed that the diffuse component of shortwave radiation reaching the CH bottom was greater than the direct component of that even though κf is higher than κd, suggesting that the shortwave radiation diffuse component reaches the CH bottom more easily than the direct component. Indeed, Fig. 11 indicates that the reaching fraction of direct component of shortwave radiation transmitted throughout the ice to the CH bottom decreases with θz. This is probably the cause of the higher sensitivity to κf, suggesting that κf is likely a relatively more important parameter compared to κd. However, there was no significant difference in sensitivity between κf and κd in the case of θz≤θc, which is the case for depths shallower than approximately 0.50 m in Fig. 9g and h. In such a case, most of the direct and diffuse components reach the CH bottom without being transmitted through the ice (Fig. 11). Figure 11 also suggests that shortwave radiation transmitted throughout the ice is dominant in the opposite case. Since CH depth ranged from 0 to 0.20 m at Site 2 in 2014, the result supports our argument that transmitted shortwave radiation mainly contributes to CH development as described for RS-exp. In the Antarctic, a lot of CHs of over 0.30 m depth were reported even though the CH mouths were covered with frozen ice (Fountain et al., 2004), suggesting that the contribution of the transmitted shortwave radiation to CH development is substantial. The results highlight the importance of separating the direct and diffuse components of shortwave radiation in CH depth simulations. Further studies of the optical characteristics of ice are necessary because there is little information on the broadband flux extinction coefficients. Furthermore, previous studies reported the development of a porous ice layer, known as the weathering crust, on glaciers worldwide during summer (Irvine-Fynn and Edwards, 2014; Stevens et al., 2022). The weathering crust density is sometimes less than 500 kg m−3 (Müller and Keeler, 1969). Because the shortwave radiation transmittance would be larger in a low-density case than that in the ice density case assumed in this study (900 kg m−3), CHs may also develop under the weathering crust. To accurately simulate temporal change in CH depth, more detailed surface ice properties, such as the weathering crust layer, should be considered for future CH modelling.
The sensitivity experiments regarding the zenith angle (θz-exp and θc-exp) suggest that differences in the zenith angles have little influence on CH depth, except for the instance in which the downward shortwave radiation always reaches the CH bottom from the hole mouth (θz = 0∘). θz-exp showed that CH depth with a higher θz was shallower, owing to a decrease in Mc (Fig. 9i), and θc-exp showed that CH depth with a lower θc was also shallow (Fig. 9j). Notably, the experiments suggest that both θz and θc hardly affect CH depth at over 15∘ and below 60∘, respectively. Snell's law states that the direct component of incident radiation is refracted through the air–water surface. According to this law, the refraction angle is approximately 20∘ smaller than the incident angle θz; therefore, the direct component of the downward shortwave radiation more easily reaches the CH bottom from the hole mouth. However, such refraction is unlikely to affect CH depth because θz is always greater than 55∘ at the studied latitude. Although CryHo calculates Mc at the centre of the CH bottom using θc, Mc may differ at the northern and southern edges of the CH bottom because the zenith angles of the CH bottom edges would differ from θc. θc-exp suggests that CH depth is temporarily non-uniform on the northern and southern edges of the CH bottom. However, CH depth likely becomes uniform again over time according to D0-exp. Indeed, the simulated CH depth using a different θc converged within approximately 2 weeks (Fig. 9j). A previous study suggests that the surface of the CH bottom on steep north-sloping ice is non-uniform due to heterogeneous radiation reaching the CH bottom (Cook et al., 2018). The slope in our studied sites is approximately 5∘. The aspect in our studied sites is southwest. Thus, the surface of the CH bottom might remain uniform on such gentle south-sloping ice.
5.3 Meteorological factors controlling the cryoconite hole dynamics
The sensitivity tests in the previous section showed that CH depth did not equilibrate in any experimental cases at the studied glacier, where meteorological conditions change before the depth reaches equilibrium (Fig. 9). It suggests that meteorological conditions significantly affect CH depth. For example, downward shortwave radiation governs the CH deepening processes (Hodson et al., 2010a; Takeuchi et al., 2018). To better understand the factors controlling the CH deepening and shallowing processes, we compared the heat components with both the CH deepening and shallowing rates. In the shallowing phase (red marks in Fig. 12), CH decay is mainly controlled by the sensible and latent heat fluxes at the ice surface (HSi and HLi) via wind speed (U). Air temperature (Ta) shows no significant correlation with the shallowing rate; conversely, U shows a strong correlation with the shallowing rate (Fig. 12a and d). However, there is no clear relationship between the input downward shortwave radiation, input downward longwave radiation, relative humidity, and direct and diffuse components of downward shortwave radiation at the ice surface (RS, RL, hr, RSd, and RSf shown in Fig. 12b, c, e, g, and h, respectively). The results suggest that strong winds induce CH decay by increasing the surface melt via the turbulent heat fluxes (Figs. 12f–j). Takeuchi et al. (2018) suggested that CHs tend to be shallower under both cloudy and windy conditions. Our analyses also suggest that windy conditions are important meteorological conditions governing the CH decay (Fig. 13a). Because the total downward shortwave is reduced through clouds, although the diffuse component of downward shortwave radiation generally increases under cloudy conditions, turbulent fluxes are likely to be dominant in the total energy budget, resulting in an ice surface that melts faster than the CH bottom. This is probably the reason why the CH collapse events were observed under cloudy conditions.
In the CH deepening phase (blue marks in Fig. 12), although the input shortwave radiation (RS) was weakly correlated with the CH deepening rate, the diffuse component reaching the CH bottom was strongly correlated with the CH deepening rate (RSfc+RStfc), whereas no correlation was found with the direct component reaching the CH bottom (RSdc+RStdc) (Fig. 12l–m). Furthermore, Fig. 12k, which shows the relationship between ice melt at the CH bottom (Qc) and the CH deepening rate, is similar to Fig. 12m. The results suggest an essential contribution of the diffuse component to CH development, which was also revealed by sensitivity analysis. However, the direct component strongly depends on the solar zenith angle and CH geometry (θz and θc); thus, the CH development is likely to be unstable during sunny conditions. Other meteorological variables are unlikely to contribute to the CH development. Although previous studies suggested that sunny conditions induce CH development (Hodson et al., 2010b; Irvine-Fynn et al., 2011; Takeuchi et al., 2018), our results highlight that the diffuse component of shortwave radiation is a key factor governing CH development, rather than merely sunny conditions. Because the CH deepening rate is likely to depend on wind speed and the intensity of the diffuse component of downward shortwave radiation, as shown in Fig. 13b, CHs could develop under sunny conditions depending on wind speed and other meteorological conditions. CH depths are mainly determined by the balance between the intensity of the diffuse component of downward shortwave radiation and the strength of the turbulent heat transfer.
5.4 Other possible factors affecting cryoconite hole dynamics
Our model does not include the effect of water lingering in CHs on the heat balance at the CH bottom because a quantitative understanding of the mechanism of convective heat transport or the buffering effect in the lingering water is insufficient. Such lingering water in CHs may affect the heat exchange between the atmosphere and CH bottom. Heat exchange should not be negligible in the case of large water surfaces in CHs. Although the water level in CHs is not estimated in the model, the refraction through the air–water surface in CHs is unlikely to affect CH depth in the studied glacier as discussed regarding θz-exp. However, such refraction might contribute to CH development in lower-latitude regions such as Asia, where the solar zenith angle is significantly smaller than that in the polar region. In addition to the refraction, reflectance at the water surface would reduce the amount of shortwave radiation reaching the CH bottom. To simulate CH depths globally, such an effect may need to be better incorporated into CryHo. Besides lingering water in CHs, the thickness of cryoconite at the CH bottom, which is not considered by CryHo, is also likely to be a key factor in determining the CH diameter and shape because a portion of the absorbed radiation in the cryoconite at the CH bottom could be transferred laterally and then melt the CH wall (Cook et al., 2010; Cook, 2012). Furthermore, the surface slope of the glacier affects the sediment thickness and thus the CH shape geometry via water flow and changes the shortwave radiation reaching the CHs (Takeuchi et al., 2000; Cook et al., 2018). Although such factors should be incorporated into CryHo in the future, CryHo is likely to reproduce the temporal change in CH depth well, even for relatively large CH on flat surfaces with diameters up to around 0.1 m, which falls within the range of CH shape geometries documented in this study. Compared to the effects of CH shape geometry and topological conditions, the effects of lingering water and sediment thickness in CHs on the CH depth may be small. The CH diameter range of 0.01 to 0.11 m assumed for ϕ-exp is similar to that observed in glaciers in western Greenland, Svalbard, and the Himalayas (Gribbon, 1979; Takeuchi et al., 2000; Cook et al., 2010), suggesting that CryHo can be applied also to these regions.
We established a numerical cryoconite hole model (CryHo) to reproduce the vertical cryoconite hole (CH) dynamics. We evaluated the model using field observations from an ice cap in northwestern Greenland collected during the 2012, 2014, and 2017 summer seasons. CryHo simulates the temporal changes in CH depth based on heat balances calculated independently at the ice surface and the CH bottom. This is a novel concept that is different from previous models simulating CH development. CryHo also considers the CH geometry, which affects the direct and diffuse components of shortwave radiation, while including the optical process of transmitted radiation throughout ice. Field observations revealed that CHs repeatedly developed and decayed over time. CryHo reasonably reproduced the observed temporal changes in CH depth, including the CH shallowing phase at four sites during three melt seasons. Sensitivity tests indicated that the CH bottom is more accessible by the diffuse component – rather than the direct component – of shortwave radiation, thereby controlling the CH development. In addition, both components of shortwave radiation transmitted throughout ice also significantly contribute to CH development, except for shallowing CH or large-diameter cases. Besides wind speed, which has been indicated by earlier studies, we revealed that the difference in albedo between the ice surface and the CH bottom is also an important factor affecting the timing of the CH collapse. Although many studies in Greenland have reported the bio-albedo effect contributing to ice surface darkening, CH collapse could also reduce the ice surface albedo not only by redistributing the cryoconite granules but also by increasing the phototrophic (i.e. algae, cyanobacteria, diatoms, etc.) and heterotrophic (i.e. tardigrade, rotifer, etc.) communities. Therefore, CHs dynamics are interactively connected to physical and microbiological processes on the ice surface. For improved projection of spatiotemporal changes in the bio-albedo effect, CryHo will be useful when coupled with regional or global climate models, although additional observations and model evaluation in various glaciers and ice sheets are necessary.
The codes of the model and the analysis and observational data used in this study are available at https://doi.org/10.5281/zenodo.8189377 (Onuma et al., 2023). Meteorological observations from the SIGMA-B site are available at https://ads.nipr.ac.jp/data/meta/A20220413-006/ (Nishimura et al., 2021).
The supplement related to this article is available online at: https://doi.org/10.5194/tc-17-3309-2023-supplement.
KF designed the study. YO, NT, and TA designed the field observations. KF developed the model with the support of MN and TA. TA parameterized the extinction coefficients of ice. All authors performed field observations. YO, KF, and TA analysed the data and wrote the manuscript. All authors contributed to the discussion.
At least one of the (co-)authors is a member of the editorial board of The Cryosphere. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work has been supported by the Arctic Challenge for Sustainability (ArCS), Arctic Challenge for Sustainability II (ArCS II), and Advanced Studies of Climate Change Projection (SENTAN) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. This study has also been partly supported by the Japan Society for the Promotion of Science (JSPS). We thank Jun Uetake, Naoko Nagatsuka, Rigen Shimada, Sota Tanaka, Ryutaro Sakaki, Koki Ishiwatari, and Akane Watanabe for their support with the field observations.
This research has been supported by ArCS (grant no. JPMXD130000000), ArCS II (grant no. JPMXD1420318865), and SENTAN (grant no. PMXD0722680395) from MEXT, Japan. This study has also been partly supported by JSPS (grant nos. 23221004, 26247078, 26241020, 16H01772, 16H06291, 19H01143, 20K19955, 21H03582, and 23K17036).
This paper was edited by Benjamin Smith and reviewed by David Chandler and one anonymous referee.
Aoki, T., Aoki, T., Fukabori, M., and Uchiyama, A.: Numerical simulation of the atmospheric effects on snow albedo with a multiple scattering radiative transfer model for the atmosphere-snow system, J. Meteorol. Soc. Jpn., 77, 595–614, https://doi.org/10.2151/jmsj1965.77.2_595, 1999.
Aoki, T., Aoki, T., Fukabori, M., Hachikubo, A., Tachibana, Y., and Nishio, F.: Effects of snow physical parameters on spectral albedo and bidirectional reflectance of snow surface, J. Geophys. Res., 105, 10219–10236, https://doi.org/10.1029/1999JD901122, 2000.
Aoki, T., Matoba, S., Uetake, J., Takeuchi, N., and Motoyama, H.: Field activities of the “Snow Impurity and Glacial Microbe effects on abrupt warming in the Arctic” (SIGMA) project in Greenland in 2011–2013, Bull. Glaciol. Res., 32, 3–20, https://doi.org/10.5331/bgr.32.3, 2014.
Bøggild, C. E., Brandt, R. E., Brown, K. J., and Warren, S. G.: The ablation zone in northeast Greenland: Ice types, albedos and impurities, J. Glaciol., 56, 101–113, https://doi.org/10.3189/002214310791190776, 2010.
Chandler, D. M., Alcock, J. D., Wadham, J. L., Mackie, S. L., and Telling, J.: Seasonal changes of ice surface characteristics and productivity in the ablation zone of the Greenland Ice Sheet, The Cryosphere, 9, 487–504, https://doi.org/10.5194/tc-9-487-2015, 2015.
Cook, J.: Microbially mediated carbon fluxes on the surface of glaciers and ice sheets, PhD thesis, University of Sheffield, UK, uk.bl.ethos.559180, 2012.
Cook, J., Hodson, A., Telling, J., Anesio, A., Irvine-Fynn, T., and Bellas, C.: The mass–area relationship within cryoconite holes and its implications for primary production, Ann. Glaciol., 51, 106–110, https://doi.org/10.3189/172756411795932038, 2010.
Cook, J., Edwards, A., Takeuchi, N., and Irvine-Fynn, T.: Cryoconite: The dark biological secret of the cryosphere, Prog. Phys. Geogr., 40, 66–111, https://doi.org/10.1177/0309133315616574, 2016.
Cook, J. M., Sweet, M., Cavalli, O., Taggart, A., and Edwards, A.: Topographic shading influences cryoconite morphodynamics and carbon exchange, Arct. Antarct. Alp. Res., 50, S100014, https://doi.org/10.1080/15230430.2017.1414463, 2018.
Cooper, M. G., Smith, L. C., Rennermalm, A. K., Tedesco, M., Muthyala, R., Leidman, S. Z., Moustafa, S. E., and Fayne, J. V.: Spectral attenuation coefficients from measurements of light transmission in bare ice on the Greenland Ice Sheet, The Cryosphere, 15, 1931–1953, https://doi.org/10.5194/tc-15-1931-2021, 2021.
Fountain, A. G., Tranter, M., Nylen, T. H., Lewis K. J., and Mueller, D. R.: Evolution of cryoconite holes and their contribution to meltwater runoff from glaciers in the McMurdo Dry Valleys, Antarctica, J. Glaciol., 50, 35–45, https://doi.org/10.3189/172756504781830312, 2004.
Fountain, A. G., Nylen, T. H., Tranter, M., and Bagshaw, E.: Temporal variations in physical and chemical features of cryoconite holes on Canada Glacier, McMurdo Dry Valleys, Antarctica, J. Geophys. Res.-Biogeo., 113, G01S92, https://doi.org/10.1029/2007JG000430, 2008.
Goudriaan, J.: Crop micrometeorology: A simulation study, Pudoc, Wageningen, the Netherlands, ISBN 902200614X, 1977.
Gribbon, P. W. F.: Cryoconite holes on Sermikavsak, West Greenland, J. Glaciol., 22, 177–181, https://doi.org/10.3189/S0022143000014167, 1979.
Hodson, A., Anesio, A. M., Ng, F., Watson, R., Quirk, J., Irvine-Fynn, T., Dye, A., Clark, C., McCloy, P., Kohler, J., Sattler, B.: A glacier respires: Quantifying the distribution and respiration CO2 flux of cryoconite across an entire Arctic supraglacial ecosystem, J. Geophys. Res., 112, G04S36, https://doi.org/10.1029/2007JG000452, 2007.
Hodson, A., Anesio, A. M., Tranter, M., Fountain, A., Osborn, M., Priscu, J., Laybourn-Parry, J., Sattler, B.: Glacial ecosystems, Ecol. Monogr., 78, 41–67, https://doi.org/10.1890/07-0187.1, 2008.
Hodson, A., Bøggild, C., Hanna, E., Huybrechts, P., Langford, H., Cameron, K., and Houldsworth, A.: The cryoconite ecosystem on the Greenland ice sheet, Ann. Glaciol., 51, 123–129, https://doi.org/10.3189/172756411795931985, 2010a.
Hodson, A., Cameron, K., Bøggild, C., Irvine-Fynn, T., Langford, H., Pearce, D., and Banwart, S.: The structure, biological activity and biogeochemistry of cryoconite aggregates upon an arctic valley glacier: Longyearbreen, Svalbard, J. Glaciol., 56, 349–362, https://doi.org/10.3189/002214310791968403, 2010b.
Holland, A. T., Williamson, C. J., Sgouridis, F., Tedstone, A. J., McCutcheon, J., Cook, J. M., Poniecka, E., Yallop, M. L., Tranter, M., Anesio, A. M., and The Black & Bloom Group: Dissolved organic nutrients dominate melting surface ice of the Dark Zone (Greenland Ice Sheet), Biogeosciences, 16, 3283–3296, https://doi.org/10.5194/bg-16-3283-2019, 2019.
Irvine-Fynn, T. D. and Edwards, A.: A frozen asset: The potential of flow cytometry in constraining the glacial biome, Cytometry A, 85, 3–7, https://doi.org/10.1002/cyto.a.22411, 2014.
Irvine-Fynn, T. D. L., Bridge, J. W., and Hodson, A. J.: In situ quantification of supraglacial cryoconite morphodynamics using time-lapse imaging: An example from Svalbard, J. Glaciol., 57, 651–657, https://doi.org/10.3189/002214311797409695, 2011.
Jepsen, S. M., Adams, E. E., and Priscu, J. C.: Sediment melt-migration dynamics in perennial Antarctic lake ice, Arct. Antarct. Alp. Res., 42, 57–66, https://doi.org/10.1657/1938-4246-42.1.57, 2010.
Kondo, J.: Meteorology of water environment, Asakura Publishing, Tokyo, Japan, ISBN 4254161107, 1994.
Langford, H., Hodson, A., Banwart, S., and Bøggild, C.: The microstructure and biogeochemistry of arctic cryoconite granules, Ann. Glaciol., 51, 87–94, https://doi.org/10.3189/172756411795932083, 2010.
Liou, K. N.: An Introduction to Atmospheric Radiation, Academic Press, ISBN 0124514510, 1980.
McIntyre, N. F.: Cryoconite hole thermodynamics, Can. J. Earth Sci., 21, 152–156, https://doi.org/10.1139/e84-016, 1984.
Müller, F. and Keeler, M.: Errors in Short-Term Ablation Measurements on Melting Ice Surfaces, J. Glaciol., 8, 91–105, https://doi.org/10.3189/S0022143000020785, 1969.
Nishimura, M., Aoki, T., Niwano, M., Matoba, S., Tanikawa, T., Yamaguchi, S., Yamasaki, T., and Fujita, K.: Quality-controlled datasets of automatic weather station (AWS) at sigma-B site from 2012 to 2020: Level 1.3, Arctic Data Archive System (ADS), Japan [data set], https://ads.nipr.ac.jp/data/meta/A20220413-006/ (last access: 1 December 2022), 2021.
Niwano, M., Aoki, T., Matoba, S., Yamaguchi, S., Tanikawa, T., Kuchiki, K., and Motoyama, H.: Numerical simulation of extreme snowmelt observed at the SIGMA-A site, northwest Greenland, during summer 2012, The Cryosphere, 9, 971–988, https://doi.org/10.5194/tc-9-971-2015, 2015.
Onuma, Y., Fujita, K., Takeuchi, N., Niwano, M., and Aoki, T.: Codes and data set for Cryoconite hole model (CryHo): Version 1.01, Zenodo [code and data set], https://doi.org/10.5281/zenodo.8189377, 2023.
Podgorny, I. A. and Grenfell, T. C.: Absorption of solar energy in a cryoconite hole, Geophys. Res. Lett., 23, 2465–2468, https://doi.org/10.1029/96GL02229, 1996.
Schneider, C. A., Rasband, W. S., and Eliceiri, K. W.: NIH Image to ImageJ: 25 years of image analysis, Nat. Methods, 9, 671–675, https://doi.org/10.1038/nmeth.2089, 2012.
Shimada, R., Takeuchi, N., and Aoki, T.: Inter-annual and geographical variations in the extent of bare ice and dark ice on the Greenland ice sheet derived from MODIS satellite images, Front. Earth Sci., 4, 43, https://doi.org/10.3389/feart.2016.00043, 2016.
Stevens, I. T., Irvine-Fynn, T. D. L., Edwards, A., Mitchell, A. C., Cook, J. M., Porter, P. R., Holt, T. O., Huss, M., Fettweis, X., Moorman, B. J., Sattler, B., and Hodson, A. J.: Spatially consistent microbial biomass and future cellular carbon release from melting Northern Hemisphere glacier surfaces, Commun. Earth Environ., 3, 275, https://doi.org/10.1038/s43247-022-00609-0, 2022.
Stibal, M., Tranter, M., Benning, L. G., and Řehák, J.: Microbial primary production on an Arctic glacier is insignificant in comparison with allochthonous organic carbon input, Environ. Microbiol., 10, 2172–2178, https://doi.org/10.1111/j.1462-2920.2008.01620.x, 2008.
Sugiyama, S., Sakakibara, D., Matsuno, S., Yamaguchi, S., Matoba, S., and Aoki, T.: Initial field observations on Qaanaaq ice cap, northwestern Greenland, Ann. Glaciol., 55, 25–33, https://doi.org/10.3189/2014AoG66A102, 2014.
Takeuchi, N., Kohshima, S., Yoshimura, Y., Seko, K., and Fujita, K.: Characteristics of cryoconite holes on a Himalayan glacier, Yala Glacier Central Nepal, Bull. Glaciol. Res., 17, 51–59, 2000.
Takeuchi, N., Kohshima, S., and Seko, K.: Structure, formation, and darkening process of albedo-reducing material (cryoconite) on a Himalayan glacier: A granular algal mat growing on the glacier, Arct. Antarct. Alp. Res., 33, 115–122, https://doi.org/10.1080/15230430.2001.12003413, 2001.
Takeuchi, N., Nishiyama, H., and Li, Z.: Structure and formation process of cryoconite granules on Ürümqi glacier No. 1, Tien shan, China, Ann. Glaciol., 51, 9–14, https://doi.org/10.3189/172756411795932010, 2010.
Takeuchi, N., Nagatsuka, N., Uetake, J., and Shimada, R.: Spatial variations in impurities (cryoconite) on glaciers in northwest Greenland, Bull. Glaciol. Res., 32, 85–94, https://doi.org/10.5331/bgr.32.85, 2014.
Takeuchi, N., Fujisawa, Y., Kadota, T., Tanaka, S., Miyairi, M., Shirakawa, T., Kusaka, R., Fedorov, A. N., Konstantinov, P., and Ohata, T.: The effect of impurities on the surface melt of a glacier in the Suntar–Khayata mountain range, Russian Siberia, Front. Earth Sci., 3, 82, https://doi.org/10.3389/feart.2015.00082, 2015.
Takeuchi, N., Sakaki, R., Uetake, J., Nagatsuka, N., Shimada, R., Niwano, M., and Aoki, T.: Temporal variations of cryoconite holes and cryoconite coverage on the ablation ice surface of Qaanaaq Glacier in northwest Greenland, Ann. Glaciol., 59, 21–30, https://doi.org/10.1017/aog.2018.19, 2018.
Tedstone, A. J., Bamber, J. L., Cook, J. M., Williamson, C. J., Fettweis, X., Hodson, A. J., and Tranter, M.: Dark ice dynamics of the south-west Greenland Ice Sheet, The Cryosphere, 11, 2491–2506, https://doi.org/10.5194/tc-11-2491-2017, 2017.
Tedstone, A. J., Cook, J. M., Williamson, C. J., Hofer, S., McCutcheon, J., Irvine-Fynn, T., Gribbin, T., and Tranter, M.: Algal growth and weathering crust state drive variability in western Greenland Ice Sheet ice albedo, The Cryosphere, 14, 521–538, https://doi.org/10.5194/tc-14-521-2020, 2020.
Tsutaki, S., Sugiyama, S., Sakakibara, D., Aoki, T., and Niwano, M.: Surface mass balance, ice velocity and near-surface ice temperature on Qaanaaq Ice Cap, northwestern Greenland, from 2012 to 2016, Ann. Glaciol., 58, 181–192, https://doi.org/10.1017/aog.2017.7, 2017.
Uetake, J., Tanaka, S., Segawa, T., Takeuchi, N., Nagatsuka, N., Motoyama, H., and Aoki, T.: Microbial community variation in cryoconite granules on Qaanaaq Glacier, NW Greenland, FEMS Microbiol. Ecol., 92, fiw127, https://doi.org/10.1093/femsec/fiw127, 2016.
Uetake, J., Nagatsuka, N., Onuma, Y., Takeuchi, N., Motoyama, H., and Aoki, T.: Bacterial community changes with granule size in cryoconite and their susceptibility to exogenous nutrients on NW Greenland glaciers, FEMS Microbiol. Ecol., 95, fiz075, https://doi.org/10.1093/femsec/fiz075, 2019.
van den Broeke, M., Reijmer, C., and van de Wal, R.: Surface radiation balance in Antarctica as measured with automatic weather stations, J. Geophys. Res., 109, D09103, https://doi.org/10.1029/2003JD004394, 2004.
Warren, S. G. and Brandt, R. E.: Optical constants of ice from the ultraviolet to the microwave: A revised compilation, J. Geophys. Res., 113, D14220, https://doi.org/10.1029/2007JD009744, 2008.
Wientjes, I. G. M., Van de Wal, R. S. W., Reichart, G. J., Sluijs, A., and Oerlemans, J.: Dust from the dark region in the western ablation zone of the Greenland ice sheet, The Cryosphere, 5, 589–601, https://doi.org/10.5194/tc-5-589-2011, 2011.
Zawierucha, K., Porazinska, D. L., Ficetola, G. F., Ambrosini, R., Baccolo, G., Buda, J., Ceballos, J. L., Devetter, M., Dial, R., Franzetti, A., Fuglewicz, U., Gielly, L., Łokas, E., Janko, K., Novotna Jaromerska, T., Kościński, A., Kozłowska, A., Ono, M., Parnikoza, I., Pittino, F., Poniecka, E., Sommers, P., Schmidt, S. K., Shain, D., Sikorska, S., Uetake, J., and Takeuchi, N.: A hole in the nematosphere: tardigrades and rotifers dominate the cryoconite hole environment, whereas nematodes are missing, J. Zool., 313, 18–36, https://doi.org/10.1111/jzo.12832, 2021.