A Scaling Relation for Cryoconite Holes

Tiny cryoconite holes are commonly found on glacier surfaces. Despite a long history of research on them, their influence on glacier‐scale mass balance and runoff are not well understood. We model the absorption of solar radiation at the bottom of cylindrical cryoconite holes, incorporating the three‐dimensional geometry. The simulated holes achieve a limiting steady‐state depth, where the daily melt rate at the bottom of the holes matches that at the glacier surface. This implies a feedback loop restricting the excess ice melt due to the presence of dark supraglacial impurities. The modeled steady‐state depth scales approximately linearly with the radius, consistent with in situ observations at several glaciers across the world. Given the areal coverage and radius distribution of cryoconite holes on a glacier, this scaling yields first‐order estimates of their melt contribution.

and out of the holes, refrozen ice lids covering the holes as commonly seen on Antarctica (Fountain et al., 2008), partitioning of the available energy between melting and sublimation at the glacier surface, heat production due to microbial activities (Fountain et al., 2004) etc. Possibly because of the complexity involved, a general understanding of the steady-state hole geometry is still lacking.For example, it is still debated if there is any depth-radius, or equivalently, volume-diameter relationship (Fountain et al., 2004) for steady-state holes (Jepsen et al., 2010;Onuma et al., 2023;Podgorny & Grenfell, 1996;Wharton et al., 1985).Such a relation, if present, will help estimate the collective glacier-scale melt contribution of cryoconite holes, bridging the scale gap mentioned before.The existing studies employing models with different levels of complexity (Cook, Hodson, & Irvine-Fynn, 2016;Gribbon, 1979;Hittson, 2010;Onuma et al., 2023;Podgorny & Grenfell, 1996;Wharton et al., 1985) did not address this problem.Most of these studies did not consider the three dimensional geometry, and/or concentrated on the transient evolution of a few observed holes (e.g., Hittson, 2010;Onuma et al., 2023).
The above gap motivates this study of the dominant effects of shortwave radiation on the transient evolution of cylindrical cryoconite holes and their steady-state geometry in three dimensions.We develop an idealized model, with a careful treatment of the hole geometry, and simulate an ensemble of cryoconite holes at two geographical locations-the Himalaya and Greenland.A compilation of published data, along with our field data from Machoi glacier (the Indian Himalaya, 34.27°N, 75.53°E), help validate the simulation results.

Idealized Model
On clear and sunny days, the energy balance of a cryoconite hole is dominated by the shortwave solar radiation.Assuming the hole radius r (m) remains relatively stable, the depth of the hole h (m) measured from the glacier surface evolves depending on the difference between the melt rates at the surface and the hole bottom with, (1) Here,   ′  (W m −2 ) is the rate of shortwave energy absorbed per unit area of the glacier surface; SW dir and SW dif are the corresponding direct and diffuse shortwave components entering through the circular opening at the surface and absorbed at the cryoconite disk; SW wal is the contribution from reflected and transmitted shortwave radiations from the wall; ρ and L m are the density (kg m −3 ) and latent heat of melting of ice (J kg −1 K −1 ).We only model the shortwave radiation to keep the model tractable.The energy lost in the form of longwave radiation, turbulent fluxes, heat advected by meltwater, etc. (Cook, Hodson, & Irvine-Fynn, 2016) (see Text S2 in Supporting Information S1) are approximately accounted for via empirical, dimensionless efficiency factors ϵ and ϵ′.The factor ϵ for the bottom of the hole is assumed to be in the range of 0.6-0.8-aballpark range consistent with previous energy-balance studies (Table S1 in Supporting Information S1).As turbulent heat losses are stronger at the glacier surface than at the bottom, the efficiency factor at the glacier surface ϵ′ was taken to be 1.05 − 1.1 (1.1 − 1.15) times larger than ϵ for the Himalayan (Greenland) site.With these efficiency factors, the present model may produce realistic daily mean melt, but cannot capture the actual diurnal variability of the loss terms.As we focus on the daily averaged steady-state depth averaged, this is not a major limitation.A list of all the variables and symbols used here is provided in Table A1.
The absorbed shortwave solar radiation per unit area (W m −2 ) on the glacier surface with broadband albedo α′, is, Here,  1 and  2 are the direct and diffuse solar insolation on a horizontal surface under clear-sky conditions.
Denoting the intersection between the circular direct solar beam and the cryoconite disk by   (Figure S1c in Supporting Information S1), the absorbed direct shortwave per unit area on the disk with an albedo α is, (3) As the factor  ∕ 2 is purely a function of the aspect ratio of the holes (Text S1 in Supporting Information S1), SW dir = SW dir (h/r).Assuming an isotropic incoming diffuse radiation, the absorbed diffuse radiation per unit area at the disk is, 10.1029/2023GL104942 3 of 9 Here Ω is the mean solid angle subtended by the hole opening, averaged over all points on the cryoconite disk.Again, Ω = Ω(h/r) (Text S1 in Supporting Information S1).The direct and diffuse shortwave radiations that enter the hole but are not incident on the cryoconite disc, fall on the wall, which reflects a fraction α′ of them.
Assuming the surface to be Lambertian, that is, a purely diffuse reflection, the total reflected shortwave energy is . The part of this reflected shortwave, which is absorbed per unit area of the cryoconite disc per unit time, is approximated as (Text S1 in Supporting Information S1), Here Ψ is the mean solid angle subtended by the cryoconite disk at the wall, and again Ψ = Ψ(h/r).The variations in the reflected solar radiation coming from different parts of the wall (e.g., Hittson, 2010) are neglected for simplicity.However, incorporating that effect will not change the fact that Ψ = Ψ(h/r), due to the parallel nature of the sun rays.The shortwave energy transmitted through the wall is absorbed at the disk.We assume a broadband exponential decay of the transmitted shortwave within the ice with depth (Gribbon, 1979), with a decay length of λ (m).The rate of absorbed transmitted energy per unit area of the disk is, The expression for χ involves the decay length so that χ = χ(h/r, λ/r) (Text S1 in Supporting Information S1).
We group together the reflected and transmitted shortwave contributions from the wall, SW wal ≐ SW ref + SW pen .

Ensemble of Cryoconite Holes
The above model was used to simulate the evolution of an ensemble of 1,000 cryoconite holes at two locations: (34.5°N, 75.5°E) in the Indian Himalaya, and (77.5°N, 69.5°W) in Greenland.These locations were selected because of the availability of field data on cryoconite holes on Machoi glacier (34.27°N, 75.53°E; this study) and Qaanaaq ice cap (77.52°N, 69.06°W; Onuma et al., 2023).We did not consider any Antarctic location to avoid complications like the sublimation of surface ice (Lewis et al., 1998) and the presence of ice lids covering the holes (Fountain et al., 2004).
A discretised version of Equation 1 was solved numerically for each hole at hourly time steps, starting with an initial depth of 0.03 m.Gridded (1° × 1°) hourly values of diffuse and direct solar radiations for clear-sky conditions and averaged over all the days of the month of June, were obtained from Clouds and the Earth's Radiant Energy System (CERES) data product "Synoptic TOA surface fluxes and clouds (SYN) 1 deg Level 3 data" (Kato et al., 2015).The corresponding mean diurnal cycles obtained by averaging over the period 2015-2022 were periodically repeated over the 30-day simulation period.For each hole, the five free model parameters, ϵ, r, α′, α, ρ, and λ, were chosen uniformly randomly from ranges 0.6-0.8,0.05-0.6m, 0.3-0.6,0.05-0.2,650-850 kg m −3 , and 0.05-0.2m, respectively (Table S1 in Supporting Information S1).The above ensemble allowed assessing the parameter sensitivity of the steady-state geometry.The hourly solar altitude ϕ for the given location was computed for 15th June using the Mathematica® function SunPosition, which was repeated for the simulation period.

Field Observations
During three field trips to Machoi Glacier in September 2019, June 2022, and September 2022, depth and radius of 82 cryoconite holes on the middle ablation zone (3,860-4,380 m a.s.l.) were measured (Table S2 in Supporting Information S1).The radius was estimated as the geometric mean of the measured radii in two mutually perpendicular directions.Some of the holes were filled up with water, partially or fully, and some were in the areas with flowing surface water.

Transient Evolution of Cryoconite Holes
Initially, the simulated shallow cryoconite holes deepened rapidly (Figure 1a).Subsequently, the rate of the deepening slowed down, leading to a daily averaged steady-state depth.The characteristic (e-folding) timescale for the transient growth increased linearly with increasing hole depth, while for a given depth it systematically grew with decreasing surface albedo (Figure 1b).For a relatively high surface albedo α′ ∼ 0.6 and a relatively small radius r ≲ 0.1 m, the steady-state depth was attained within a few days (Figure 1b).In contrast, holes with radius r ∼ 0.5 m and α′ ∼ 0.6, took a couple of weeks to reach the steady state.The parameters α, ρ, and λ had relatively little systematic influence on the response time (not shown).

Steady-State Hole Geometry and Energy Balance
The above steady state corresponds to a steady daily mean depth (Figure 1a), with daily averaged melt rates at the glacier surface and hole bottom being equal.A relatively small periodic, diurnal variation of the depth was present at subdaily scale (Figure 1a).
The simulated steady-state depth increased approximately linearly with hole radius (Figure 2).The steady-state radius-depth relationship depended on latitude, with shallower steady-state holes at the higher latitude site, when all other parameters remained similar.A tendency of sublinear dependence was noticeable for deeper holes (Figure 2) and in long response-time limit (Figure S5 in Supporting Information S1), where a steady state might not have been achieved over the simulation period due to a slow response.At any given radius, the steady-state depth decreased strongly with the ratio of absorptances of the ice and cryoconite  1− ′ 1− (Figure 2), and with the surface albedo α′ (Figure S3 in Supporting Information S1).The impact of the variation of cryoconite albedo α was relatively less pronounced (Figure S3 in Supporting Information S1).The steady-state depth did not have any apparent dependence on either the ratios of the efficiency factors ϵ/ϵ′ (Equation 1), or ice density ρ, or penetration depth λ (Figures S3 and S4 in Supporting Information S1).

Depth-Radius Scaling
All the shortwave terms in Equation 1, except SW pen , depended solely on the ratio r/h.Therefore, when SW pen is relatively small, scaling both r and h by the same constant factor (geometric scaling) does not change the energy-balance.This implies for given fixed values of the parameters (α, α′, λ, ϵ, and ϵ′) and forcing, the steady-state holes have the same aspect ratio.This is confirmed by the simulation results (Figure 2).The value of this ratio varies with the choice of model parameters and the forcing.The ratio of ice-surface to cryoconite absorptances 1− has a strong control on it (Figure 2).As both solar elevation and the ratio of direct to diffuse shortwave radiation depend on the geographical location, the aspect-ratio is also expected to vary with location and season.
A strong independent evidence in favor of the above linear scaling comes from in situ observations on several glaciers from Antarctica, Greenland, and the Himalaya (Figure 2c).The compiled data cover a wide range of latitude, meteorological condition, and hole morphology, indicating the robustness of the scaling.Despite the idealized theory and a limited number of observations, the consistency among the theory, simulations, and field observations supports the claim of a linear depth-diameter scaling of steady cryoconite holes, which is the key result in this letter.
The observed holes are unlikely to have been in a steady state at the time of measurement.This, and the processes ignored in our theory (Text S1 in Supporting Information S1), including time variations of model parameters, are likely causes of the scatter in the observed depth-radius relationship (Figure 2c).However, the figure also indicates that the scaling relation may not be restricted to the steady states, and is more general.Additional simulations with time-varying model parameters confirms this (Figure S6 in Supporting Information S1).This is not unexpected from a scaling theory perspective (e.g., Bahr et al., 2015).Interestingly, both simulations and observations suggest an increase in scatter with latitude (Figure S7 in Supporting Information S1).Some of the previously reported models (e.g., Podgorny & Grenfell, 1996;Wharton et al., 1985), which ignored the three dimensional geometry of the holes, predicted an uniform steady-state hole depth under a given climatic condition.Both our model results and the global-scale in situ data (Figure 2c) contradict that.The three-dimensional hole geometry appears to be a critical ingredient to model cryoconite holes.(Onuma et al., 2023), Antarctica (Fountain et al., 2004;Weisleitner et al., 2020), and the Himalaya (this study; Takeuchi et al., 2000) suggests a similar linear scaling of the depth of cryoconite holes with their radius, albeit with significant scatter.

Energy Fluxes
Whenever the hole bottom absorbs energy at a higher rate relative to the surface, the holes deepen (Equation 1).This reduces the shortwave energy available and the melt at the bottom continuously, till the melt matches that at the surface, and a steady state is attained.With all other factors remaining the same, a high (low) glacier surface melt implies a shallower (deeper) steady-state hole (Figure 2).This is consistent with the decrease of simulated steady-state hole depth with decreasing surface albedo α′ (Figure S3 in Supporting Information S1), and is supported by observations on glaciers in the Himalaya (Takeuchi et al., 2000) and Greenland (Onuma et al., 2023).Since glacier surface melt increases at lower elevation, shallower steady-state holes are expected toward the glacier terminus.
The available daily mean shortwave energy at the surface (Kato et al., 2015) were similar in both the Himalayan (376 W m −2 ) and Greenland (378 W m −2 ) sites.While the daily maximum was higher in the Himalaya, longer daylight hours largely compensated for it in Greenland (Figure 3; Figure S3 in Supporting Information S1).However, the decay of the available energy at the hole bottom with depth was faster in Greenland (Figures 3a  and 3b), likely due to a lower solar elevation.This resulted in shallower simulated holes in Greenland (Figures 2a  and 2b).In contrast, field data showed relatively shallow holes on Himalayan glaciers (Figure 2c), possibly due to a systematically lower surface albedo on Himalayan glaciers than considered in our simulations causing a relatively higher surface melt (Rashid et al., 2021;Takeuchi et al., 2000).

Effects of Cryoconite Holes at Glacier-Scale
The steady-state depth sets an upper limit on the excess melt at cryoconite holes during their transient growth, compared to the surrounding glacier surface.During a period of steady weather, once the holes become steady, there is no excess daily averaged melt.Thus, the formation of cryoconite holes efficiently limits excess ice loss when dark, low-albedo impurities are present.
While the volume of all the cryoconite holes on a given glacier may be impossible to measure directly, field surveys can obtain the frequency distribution of hole radius P(r) and their fractional areal coverage f c , both of which may vary with elevation.With such data, an upper bound on the glacier-scale excess melt contribution of cryoconite holes per unit area b c (m yr −1 ) can be estimated as, Here h refers to the steady-state depth of the holes.A limited survey can be used to establish the best-fit aspect ratio on a given glacier, obtaining the linear function h(r).Then the above integral can be evaluated to estimate b c .Thus, the depth-radius scaling relation circumvents the scale gap discussed before, and is of practical importance.This highlights the need to study the shape of  () and its spatio-temporal variability.
For an order-of-magnitude estimate of the total excess melt at the cryoconite holes, let's assume an arbitrary surface coverage of f c = 0.1, uniform hole radius of r ∼ 0.02 m), h/r ∼ 2, α′ ∼ 0.5, and α ∼ 0.1 at the end of the ablation season.These parameters are reasonable for Himalayan (e.g., Takeuchi et al., 2001) or Antarctic glaciers (e.g., Fountain et al., 2004).Then the mean albedo of the surface including the holes, when viewed from the top, is ∼10% lower than the clean-ice albedo.However, the excess melt per unit area on the glacier with steady cryoconite holes is only f c × h ∼ 0.004 m over the entire ablation season, which is negligible for all practical purpose compared to a few meters of annual ablation typically observed.With a higher fractional coverage, or deeper and wider holes, the excess melt may be higher, but still not comparable to the clean-ice melt rates.The feedback associated with the formation of cryoconite holes is thus very efficient in restricting excess melt in the presence of dark impurities on glacier surface.

Effect of Simplifying Assumptions
The present three-dimensional model considers the effect of shortwave radiation, ignoring several complicating factors (see discussions in Text S2 of the Supporting Information S1), which may potentially cause departures from the proposed scaling.However, as the depth-diameter scaling predicted by the idealized theory is consistent with available field observations across the globe, the processes ignored here either obey the scale symmetry, or are not large enough to destroy it.A more detailed model study, incorporating some of these processes, can be used to ascertain that.A larger in situ data set on the geometry and energy balance of cryoconite holes, spanning more glaciers from different climatic zones and representing different seasons, will also be helpful.

Conclusions
An idealized model that captures the most dominant physical process often provides useful insights into complex systems.We build an idealized model of the effects of shortwave radiation on a cylindrical cryoconite hole in three-dimension.Both the theory and simulations suggest that the holes achieve a steady-state depth under steady weather conditions, where the daily averages of melt rates at the glacier surface and at the hole bottom are equal.This state is characterized by a linear depth-diameter scaling relation, which is supported by limited, but global-scale observations.This relation allows up-scaling of the effects of cryoconite holes to glacier scale, given the distribution of the hole radius and their fractional coverage.The geometric feedback from a deepening cryoconite hole reduces the shortwave absorption at the low-albedo cryoconite, efficiently limiting excess melt on glacier surface in the presence of dark surface impurities.More field data on the hole geometry and local energy-balance processes can help validate the proposed depth-diameter scaling over a wider range of geographical and climatic settings.

Appendix A: Notations
See Table A1 for a list of the symbols used in this letter.AB conceived the study, in discussion with CS, IR, and SV.AB developed the theory, with inputs from CS. AB did the numerical simulations, analyzed the data, and wrote the manuscript.NAN and AMC performed the field measurements, with help from IR and AB.All authors discussed and edited the manuscript.The people of Matayan village (Ladakh), Danish, Faisal, Imtiyaz, Saquib, and Ulfat helped during the field work.Useful inputs from an anonymous reviewer, reviewer Yukihiko Onuma, and the editor are acknowledged.

Figure 1 .
Figure 1.(a) Initially, the hole depth grows rapidly, but finally tends to a steady diurnally averaged value at both the locations.(b) The e-folding response time needed to reach the steady depth scales approximately linearly with hole depth.The dashed gray line (y = 10x) is a guide to the eye.The symbol colors denote the albedo of the surface ice.A small fraction (∼10%) of the holes with final depths smaller than the initial depth are excluded from the plot.

Figure 2 .
Figure 2. (a and b) Simulated steady-state cryoconite holes are characterized by an approximately linear scaling of the depth with the radius at both the locations.Symbol colors denote the ratio of ice-surface to cryoconite absorptances,  1− ′ 1− .(c) In situ observations at glaciers in Greenland(Onuma et al., 2023), Antarctica(Fountain et al., 2004;Weisleitner et al., 2020), and the Himalaya (this study;Takeuchi et al., 2000) suggests a similar linear scaling of the depth of cryoconite holes with their radius, albeit with significant scatter.

Figure 3 .
Figure 3. (a and b) Diurnal variation of direct (SW dir ) and diffuse (SW dif ) shortwave radiations coming through the hole, and the diffuse radiation coming from the wall (SW wal ) that are absorbed at the cryoconite disk.The sum of the three components is SW tot .While SW dir and SW dif are uniquely determined by the depth to radius ratio, SW wal depends on the penetration depth (see text for details) leading to the scatter seen.(c and d) Diurnal cycles of the shortwave absorbed at the glacier surface  (  ′  ) and at the cryoconite (SW tot ) for the same holes as in this figure.The three components of SW tot are shown.