Between July 30, 2016 23:42 AKDT (Alaska daylight time) and August 30, 2016 20:14 AKDT, a field experiment monitoring sub-debris glacier melt and temperature at an array of locations was conducted (data gaps will limit different experiments in this study to subsets of this time range). Using a FLIR T620 camera (uncooled microbolometer, 640 × 480 resolution, spectral range 7.5–14.0 μm, accuracy ±2°C, thermal sensitivity <0.04°C at 30°C), surface temperature was monitored at a 15 min interval for 164 h, broken into six observations periods under different meteorological conditions. Due to data gaps, there were 102 h, split between four observation periods (1 August 11:42 AKDT to 3 August 10:57; 3 August 3:45 to 4 August 16:13; 10 August 21:47 to 11 August 8:02; and, 16 August 23:15 to 17 August 19:45) where measurements from all data sources were collected simultaneously. The thermal camera was mounted on a surveying tripod and deployed at the same coordinate location for each of the six periods. Emissivity was set to one and distance set to 0 to facilitate post-processing. Within each oblique thermal image, seven distinct image segments were captured, each with a near-homogeneous debris thickness (Supplementary Figure S1). The segments, ordered in distance from the thermal camera and corresponding to locations shown in Figure 1, are:

1. 38 cm (average) debris thickness (10 m from sensor)

Air temperature and relative humidity are required input to solve for surface temperature from a thermal image. Measurements of 1.5 m air temperature and relative humidity were collected at a 10 min interval throughout the duration of this study at site 1 (on top of structure 8) labeled in Figure 1. The measurements were made using an ONSET HOBO U23 Pro v2 Temperature/Relative Humidity Data Logger housed in a radiation shield. A second, identical structure and sensor configuration was deployed to the same debris-covered location as well as a bare ice location (Figure 1, site 5) in 2012 spanning the same period of time as the data collected in 2016 (31 July 23:42 AKDT to 30 August 20:14 AKDT). These 2012 data are used in this study to amend a 2016 sensor deployment deficiency, providing hourly correction factors to quantify the difference in air temperature over debris cover and bare ice which is notably different (Supplementary Figure S2). To provide context to analysis of debris saturation, precipitation was also measured at site 1 (on top of structure 8) labeled in Figure 1, using a HOBO tipping bucket Rain Gauge Data Logger.

Thermal infrared cameras are often sold along with integrated software packages that facilitate post processing. There is nothing inherently wrong with using proprietary software packages, but given the complexities and scales that glacier research is conducted on, it is helpful to have oversight/control of the equations used to decompose the several entangled signals present in at-sensor radiance values. In contrast to other geophysical applications of thermal imagery, e.g., volcanology where dramatic thermal anomalies (+100s of degrees) are the signal of interest (Spampinato et al., 2011), applications in glaciology mainly use absolute temperature values that have a relatively low deviation (∼10–30°C over varying debris thicknesses), thus further increasing the importance of careful image processing. Here, I present a method to process field-based, oblique time-lapse thermal imagery which poses unique difficulties from, 1) variable atmospheric attenuation that can be present in one frame imaging several surfaces at different distances (in this study, ranging from a few meters to a kilometer away from the sensor); 2) variable atmospheric attenuation through time; and 3), image shift and rotation from both imaging an unstable landscape and, in the case of this study, acquiring the image time-series from an unstable location on the glacier surface. The method uses Thermimage (Tattersall, 2017) to automate the extraction of raw values from FLIR images and (Tỳč and Gohlke, 2015) to co-register images in the series. This processing routine was developed specifically for a glacier setting by including atmospheric corrections accounting for variability in near-surface conditions over bare and debris-covered ice surfaces present along the line of sight distance from the thermal camera sensor to the target surface. Aside from this correction, which can be omitted, the image processing routine is suitable for any Earth science application with a few meters to a few kilometers line of sight surface to sensor distance. The code uses both R and Python packages combined in one open-source Jupyter Notebook available at (https://github.com/samherreid/ThermalTimelapse).

The value assigned to each pixel in a thermal image is a quantization of the net radiant intensity, W _tot (Wm⁻²), received at the sensor within a set, sensor specific, spectral range (7.5–14.0 μm for the camera used in this study). Following the formulation from Usamentiaga et al. (2014), the surface temperature of an object, T _s (K), can be calculated by: T s = W t o t − ( 1 − ε o b j ) * τ a t m * σ * ( T r e f l ) 4 − ( 1 − τ a t m ) * σ * ( T a t m ) 4 ε o b j * τ a t m * σ 4 (1) where ɛ _obj is the emissivity of the object, τ _atm is the transmittance of the atmosphere between the sensor and the object, σ is the Stefan-Boltzmann constant (5.67 × 10^–8 Wm⁻²K⁻⁴), T _refl (K) is the reflected temperature and T _atm (K) is the temperature of the atmosphere.

None of these quantities (apart from σ) are constant in time or space. However, variability in ɛ _obj through time, e.g., from rock surfaces becoming wet, was assumed to be negligible. In this study, ɛ _obj was varied in space based on site specific values extracted from an ASTER image, acquired on August 30th, 2016 at 13:25 AKDT, that was processed to estimate surface emissivity (AST05) (following Gillespie et al., 1998; Abrams et al., 2002). These values averaged to a debris cover emissivity of 0.94 and an ice emissivity of 0.97. T _refl was assigned for each thermal image by extracting the mean temperature of the aluminum poles of the weather station structure visible in the field of view of each image (Figure 1). This loosely follows the reflector method described in Usamentiaga et al. (2014); however, variation from vectors normal to a smooth cylindrical aluminum surface rather than normal to a complex, randomly oriented crumpled and then flattened aluminum surface (a more favorable configuration) was not quantified. Atmospheric temperature (T _atm , Section 3.2) measurements coincident in time with each thermal image were extracted, and a correction (derived from local data collected 4 years earlier in 2012) was applied to account for air temperature variability for the portions of each thermal image that were debris-free. A set of correction factors were computed for each hour of the day by finding the difference between hourly averaged 1.5 m air temperature collected at location Figure 1 1) and the same measurement recorded at location Figure 1 5) (Supplementary Figure S2).

The final parameter needed to solve Eq. 1, τ _atm , can be formulated as the product of the two main quantities that cause signal attenuation received at the thermal camera, τ a t m = τ m τ s , (2) where τ _m is molecular absorption by constituent gases and τ _s is scattering by particles in the atmosphere (Gaussorgues, 1994). Following Gaussorgues (1994), τ _m is simplified to account for the two dominant constituent gases of the atmosphere, water vapor and gaseous carbon dioxide: τ m = τ H 2 O τ C O 2 . (3)

As electromagnetic radiation travels through the atmosphere from the target glacier surface to the infrared sensor, some of the radiation is absorbed by atmospheric water vapor molecules (Gaussorgues, 1994). Discrete volumes of interest requiring a solution for atmospheric water vapor content can be approximated as a solid angle or ellipse-based cone, where the point of the cone is at the sensor and the ellipse base approximates the glacier surface radiation source area surrounding each (rectangular) pixel in a thermal image. The number of water vapor molecules present within this cone is a function of the local partial pressure of water vapor and the presences of gaseous water vapor molecules (Gaussorgues, 1994). These quantities are governed both by predictable factors (e.g., elevation, diurnal and seasonal cycles) and chaotic factors (e.g., wind and weather systems). Considering these factors, τ H 2 O cannot be considered static even at sub-hourly time scales. For this study, the volume of the ellipse-based cone was simplified to a one-dimensional distance, x, from the object surface to the infrared sensor. This is a practical simplification even in an oblique setting because objects that are close to the sensor have a sufficiently disproportionate pixel resolution to the length of x scale ratio of millimeters to meters and objects imaged in the distance have a scale ratio of centimeters to 100s of meters. The quantity of water vapor along x is most commonly expressed as a height of precipitable water, which is the amount of liquid water that would result from the condensation of all of the present water vapor molecules (Gaussorgues, 1994). Precipitable water, h, can be expressed as h = ∫ 0 x q v ρ d x , (4) where q _v (dimensionless) is specific humidity, the mass mixing ratio of water vapor to the total mass of the moist air along x, and ρ (kg m⁻³) is the density of the moist air. h has units of kg m⁻² which is equal to a one-dimensional height of water in mm. h is frequently solved for over a vertical column from the ground to the top of the atmosphere. Many parameterizations exist for this quantity (e.g., Maghrabi and Al Dajani, 2013), but are not applicable to a horizontal x with a variable length and also where pressure and temperature profiles cannot be approximated as simple functions of elevation. h is therefore estimated using measured relative humidity, RH (dimensionless), T _atm (K) and sensor to surface distance, x (m).

The longest x over which h was solved for was 1 km with variable ground-surface types (alternating between debris cover and bare glacier ice, Figure 1C) where RH and T _atm cannot be assumed constant. For this study, the landscape imaged repeatedly was broken into segments where each segment contains area with a near constant debris thickness (Supplementary Figure S1) and a similar ground to sensor distance, x (Figure 1; 1–7). For each segment, h is approximated from the quadrature of Eq. 4 by h = q v i c e ρ i c e x i c e + q v d e b ρ d e b x d e b , (5) where the distance, x, is the sum of two components accounting for variability in the atmosphere above two distinct surfaces, bare ice, x _ice , and debris cover, x _deb (Figure 1). Using a set of standard equations, T _atm and RH can be used to solve for q _v and ρ over bare ice ( q v i c e , ρ _ice ) and over debris cover ( q v d e b , ρ _deb ) (Supplementary Appendix). With a solution for h from Eq. 5, spectral transmittance through the atmosphere considering molecular absorption of water vapor, τ H 2 O , can be estimated as a function of wavelength (λ, μm) (Passman and Larmore, 1956; Gaussorgues, 1994) (Supplementary Appendix). Numerical integration of τ H 2 O ( λ ) over the specific thermal camera spectral range (in this study 7.5–14.0 μm) enables a solution of τ H 2 O .

Solving for spectral transmittance through the atmosphere considering molecular absorption of gaseous carbon dioxide, τ C O 2 and signal attenuation from scattering by particles in the atmosphere, τ _s , can be estimated more simply as functions of λ and sensor to object distance, x (Gaussorgues, 1994) (Supplementary Appendix). Together, these calculations of τ H 2 O , τ C O 2 and τ _s enable an estimate of τ _atm for each image segment that accounts for the specifications of the thermal camera, the physical setting of the experiment and the instantaneous near-surface atmosphere.

The unique setting of imaging a glacier surface with simultaneous measurements of non-zero melt means that bare glacier ice can be expected to be at the pressure melting point. While small impurities are present on even bare ice surfaces with the potential to raise the surface temperature slightly above 0°C, I used this setting of a known, in-frame temperature to identify and correct an assumed linear −2°C sensor bias for all of the images used in this study.

The procedure described here forms the basis of a thermal image processing routine that considers variable surface types, variable line of sight distances and variable near-surface atmosphere along the line of sight distance within a single thermal image. This routine was applied to 684 images acquired over a total span of 164 h. This process was fully automated including automated rotation, translation and scaling data shifts (Tỳč and Gohlke, 2015) to co-register images to match manually defined image segments (black lines in Figure 1B). Locations in the image where unnatural objects were present (e.g., an ablation stake and the corner of a weather station) were removed to not disrupt the mean, median and percentile values computed for each segment.

Contemporaneous with the thermal image time-series, sub-daily melt measurements were made within four of the seven segments defined in Figure 1: Segments 1, 2, 4, and 5 with debris thicknesses 38, 8, 4, and 0 cm, respectively. Melt measurements were made using a “glacier selfie stick” approach where a graduated, rigid aluminum ablation stake was drilled into bare, or sub-debris ice with a visible spectrum time-lapse camera loosely attached and floating on the surface to photograph the progressive exposure of graduations 5 times per day (Figure 2, Supplementary Video S1). For the ablation stake located closest to the thermal camera [Figure 1 (8), 38 cm debris thickness], the visible spectrum time-lapse camera was fixed to the floating meteorological station. These measurements enabled a sub-daily measurement of melt and were validated against periodic, field-based manual measurements. The presence of the camera assembly at the ablation stake also impacted melt, which was most visible at the bare ice location. Instances where the camera assembly caused clear excess lowering were corrected for by adding back the height of an unnatural (few cm) trench, yet these instances were rare and the assembly remained largely flush with the visible surrounding surface (Supplementary Video S1). Each of the four image time-series were post-processed to derive melt rates by manually measuring the graduation exposure rate. While the cameras acquired five images per day, the frequency of melt measurements were also a function of melt. If the melt rate between images was low, e.g., for the locations with a thicker debris cover, measurements intervals were lengthened until a discernible change could be confidently measured. This makes selecting a meaningful common dt for the computed melt rate difficult. For the purposes of the modeling effort described in Section 3.6 below, where high frequency (10 or 15 min) measurements were downsampled to a 1 h frequency, these melt measurements were resampled to 1 h. This resampling produces error in the melt model validation, but preserves the diurnal variation in the melt rate below thin debris cover and at bare ice. This diurnal variation is a unique signal since most melt measurements of bare glacier ice or sub-debris melt are averaged over days to weeks or months. The bare ice ablation stake, with the highest melt rate, needed to be occasionally re-drilled to maintain a continuous record. These disruptions to the melt record and subsequent post processing proved difficult to mitigate, where disruptions to the melt rate only causes obvious momentary unrealistic values, while a large and more subtle error can accrue in the cumulative melt signal. Large errors were manually smoothed by taking an average between the surrounding days melt measurements at the same time of day. Finally, to have all of the ablation records begin at the same time, linear regression was used to solve for melt over the hours needed to be subtracted from earlier emplaced ablation stakes to align with the initial measurement of the last emplaced ablation stake.

Melt measurement error was estimated at 0.12 cm per measurement, which sums to 0.6 cm per day if all five images were suitable for use. This value is the result of watching the landscape shift in the background of each image which is caused by movement of the camera assembly unrelated to melt-driven surface lowering. Figures 2C,D shows an extreme shift event where the camera tilted backwards overnight giving the impression of glacier accumulation. While the distance to the background landscape is unknown as well as the angle of rotation, it is still a clear indication of when the measurement will contain an error, requiring a correction factor that will carry through the remaining melt measurements when solving for cumulative melt. A video (Supplementary Video S1) of the complete visible time-lapse datasets from the four locations shows that while events like the one shown in Figures 2C,D do occur, through much of the time-series, the camera position (and background landscape) remains stable, facilitating meaningful melt measurements.

For the three locations of measured ablation with debris cover [Figure 1(1,2,4)], temperatures were recorded within the debris layer in 4 cm increments from the debris-ice interface to the surface. These measurements were collected using ONSET U23 Pro v2 External Temperature Data Loggers at a 10 min increment (resampled to 1 h for the modeling effort, described in Sections 3.4, 3.6). At the location of the thickest debris cover [Figure 1(1), 38 cm debris thickness], further analysis was conducted to investigate the energy transfer within the debris layer concurrently with the ablation and thermal time-lapse measurements. Thermal diffusivity, in a strictly conductive energy transfer setting, can be solved for using the finite-difference method to approximate the Biot-Fourier diffusion equation, ∂ T ∂ t = κ c ∂ 2 T ∂ h d 2 (6) where T is temperature (°C), h _d is depth (mm), t is time and κ _c is thermal diffusivity (mm ² s ⁻¹) (Hinkel et al., 1990; Conway et al., 2000; Rowan et al., 2021). The subscript, c, denotes the exclusion of any heat transfer mechanism besides conduction.

A second approach was used to solve for thermal diffusivity from diurnal amplitude decay with depth, incorporating all heat transfer mechanisms. In order to automate the selection of diurnal maximum and minimums, a Savitzky-Golay filter was applied (polynomial order 3, window length 51) to smooth the jagged nature of a surface, or near-surface temperature profile as clouds pass overhead. Considering the change in temperature amplitude with depth for each day, amplitude derived thermal diffusivity, κ _a can be calculated, κ a = ω ( h 1 − h 2 ) 2 2 ( ln ( A ( h 1 ) ) − ln ( A ( h 2 ) ) ) 2 (7) where A(h ₁) and A(h ₂) are temperature amplitudes at depths h ₁ and h ₂, respectively Wang et al. (2020). ω is the radial frequency, ω = 2 π ϕ (8) where ϕ is the period of the diurnal cycle set to 86,400 s (24 h).

Using the smoothed temperature profiles, phase shift was calculated at each measurement depth by finding the lag in time between the maximum temperature at depth and the maximum temperature at the surface for each day. The slope of phase shift with depth is also proportional to thermal diffusivity (Anderson, 1998).

The simultaneous hourly measurements of sub-debris melt under variable debris thicknesses and the corresponding surface temperatures offer an ideal framework to evaluate, and possibly improve upon, methods to solve for debris thickness and sub-debris melt.

Exponential scaling is often used to describe the nonlinear relation between surface temperature and debris thickness (e.g., Mihalcea et al., 2008b), yet local variability intrinsic to the relation are embedded within the model parameters, limiting transferability. Here, I propose a method similar to Kraaijenbrink et al. (2017) that draws on information present in the same thermal image used to estimate debris thickness to constrain the exponential scaling. The approach here differs from Kraaijenbrink et al. (2017) by 1) introducing one model parameter that I hypothesize will be stable through time and for other locations on Earth; and 2) looking both within the glacier domain as well as outside the glacier domain, to the local valley walls, to constrain a realistic relation between debris thickness and surface temperature even if there is no thick debris cover present within a thermal image (or no debris cover at all).

I define T s * as the “warmest local temperature” extracted from the thermal image. T s * is used as a reference surface temperature representing the local radiative forcing without heat sinks. For any location with little or no prior knowledge of the setting, for example, a routine iterating through all glaciers on Earth, T s * can be defined as T s * = max ( max ( T s ) , max ( T s b u f f e r ) ) , where T _s is surface temperature across the glacier domain and T s b u f f e r is off-glacier surface temperature within a set buffer distance outside of the glacier domain to capture the local valley walls. The assumption behind looking off-glacier is that the maximum temperature in close proximity, and at a similar elevation, will be from a rock surface with a similar emissivity and surface roughness to a thick debris cover that is decoupled from heat sinks regardless of whether there is thick, thin or no supraglacial debris present at that location. The maximum function is used to discard less ideal, e.g., vegetated, surfaces that will likely be cooler than till or bedrock. Measurements from this study show the surface temperature of thick debris cover can exceed even south facing valley wall temperature, but the two quantities, thick debris T _s and valley wall T _s , are correlated (Figure 3). Since the field experiment presented here is conducted with the knowledge that there is thick debris cover, T s * is set to max(T _s ) which, for this time-series, is equivalent to T _s at the 38 cm debris thickness segment. If using only max(T _s ) or max ( T s b u f f e r ) , it is possible for T _s to exceed T s * which can be accommodated by using a piecewise function. Where available, T s * is taken from the thermal camera data. For the sub-debris melt analysis expanding beyond the thermal camera observation windows, T s * is set to the temperature measured by the contact thermistor at the top of the thermal profile at the 38 cm debris thickness segment (Section 3.5).

With T s * established, a relation to solve for debris thickness, h _deb , can be expressed as, h d e b ( T s ) = min a 9 ⁡ exp T s * T s * − T s , h m a x , 0 ≤ T s < T s * h m a x , T s ≥ T s * (9) where a ₉ is the single model coefficient (subscript referring to the equation number in this paper, subsequently developed model coefficients, a _x , b _x and c _x , will follow the same, subscript to the equation number, notation) that I hypothesize will be stable through space and time because much of the local and regional variability, both in terms of micrometeorology and broad geographical/orographical effects, will be captured in T s * . h _max is a user defined value of “thick” debris cover where associated sub-debris melt rates can be assumed to have asymptotically approached a stable, low value. Here, I set h _max to 40 cm. While an exponential function is generally successful at capturing the relation between surface temperature and debris thickness for debris thicknesses around 40 cm or less, towards thicker debris cover, it is unrealistic that small, incremental increases in surface temperature should return exponentially thicker debris cover. Lacking data to support, for example, a third-order polynomial fit that would smoothly transition to a constant debris thickness while surface temperature continues to increase, I have opted to flat-line debris thickness to make clear that what can be safely resolved is that the debris is thick, but the model is unable to predict an absolute thickness. This method also cannot predict debris thicknesses when temperatures are below freezing.

I compare this approach to two recently proposed/published methods from Rowan et al. (2021) and Rounce et al. (2021). Analysis by Rowan et al. (2021) shows a power-law relation between debris thickness and near-surface temperature which they invoke to suggest solving for sub-debris melt from surface temperature data. Following this framework, I evaluate the inverse of the power-law from (Rowan et al., 2021, Figure 7A), where h d e b ( T s ) = a 10 T s b 10 . (10) I also evaluate the approach used by Rounce et al. (2021) where debris thickness is estimated using the Hill equation, h d e b ( T s ) = T s b 11 c 11 a 11 − T s 1 c 11 . (11)

With h _deb either known or estimated from one of the methods described above, it becomes feasible to estimate sub-debris melt rate, b ̇ . For the method I propose here, I continue to use T s * as defined above, but now use it to first solve for debris cover surface temperature. This may sound redundant, since T s * was derived from distributed surface temperature data. However, the above method was developed to solve for stable (over the timescale considered here) debris thickness. To solve for debris thickness, especially over large spatial scales, there are at best a few suitable (cloud and snow free) satellite-based thermal images acquired per year. To force a method to derive b ̇ , T s * is not discovered, as above, but rather measured (or derived) continuously from a fixed location. For this study, T s * is again set to T _s at the “thick” 38 cm debris thickness segment from thermal camera data where available, and from the contact thermistor at the debris surface for the longer time-scale analysis. In future applications, T s * can be measured from either a contact thermistor at the surface of thick debris (which is an exceptionally easy sensor to deploy in the field) or by taking repeat thermal images below as much of the cloud canopy as possible which could be possible from a vantage point off-glacier where a power source might be more readily available.

Using T s * and one model coefficient, distributed debris surface temperature, T _s , can be estimated at the same temporal interval as T s * where debris thickness is known, T s ( T s * , h d e b ) = T s * − ( T s * ⁡ exp ( − a 12 h d e b ) ) . (12) With distributed T _s derived through time, b ̇ can be estimated as, b ̇ ( T s * , h d e b ) = a 13 ( 1 − b 13 ) ( T s * − T s ( T s * , h d e b ) ) + a 13 b 13 T s * , T s * ≥ 0 0 , T s * < 0 (13) where a 13 = b ̇ 0 / T s * , the ratio of bare ice melt rate, b ̇ 0 , and T s * . This term scales the relation to melt. I hypothesize a ₁₃ will remain relatively stable through space and time because it is a ratio of two physically related quantities. If true, this would mean that b ̇ 0 would not need to be explicitly derived from other melt modeling methods to solve for bare ice and sub-debris melt rates. a ₁₃ will likely need to be recalibrated between geographic regions where the energy balance terms are not similarly weighted and the hypothesis may be rejected where the weighting is locally variable, e.g., where turbulent fluxes are an important component of the energy balance (Steiner et al., 2018). b ₁₃ is a term that accounts for the thermal inertia of a debris cover, where melt below thick debris is forced by energy transfer occurring over a period longer than the 24 h diurnal cycle. Here, I parameterize b ₁₃ as b 13 = median ( b ̇ 38 ) / median ( b ̇ 0 ) , where b ̇ 38 is the melt rate below a “thick” 38 cm debris cover (note b ̇ 38 is a melt rate, while b ₁₃ is a model parameter). This ratio is the empirical middle value of a fairly difficult to obtain set of measurements. It scales the empirical asymptotic convergence of the full (in time) melt rate dataset collected for this study and is possibly stable through space and time. If stable, this would facilitate transferability and reduce the frequency these cumbersome measurements would need to be repeated. With the two model parameters expressed in their expanded ratio form and explicitly solved for (meaning without taking the median values), the last term of Eq. 13 reduces from b ̇ 0 T s * b ̇ 38 b ̇ 0 T s * , to b ̇ 38 . This sets the melt rate for the region of the curve after the exponential decay approaches the asymptote. Because the parameterized form is scaled by T s * , melt will not continue throughout the winter as surface temperature drops to 0°C, or below, from the piecewise case of sub-freezing temperatures.

While Eq. 12 has the same arguments as Eq. 13 and could thus be easily written into Eq. 13, it is worthwhile to consider it as a separate step because surface temperature measurements are easier to obtain than sub-debris melt measurements for the purposes of method validation.

The analysis presented in this study is over a small portion of one glacier where factors like elevation dependent temperature lapse rates can be neglected. In order to derive distributed results from the new methods presented here, T s * will need to be extracted at elevation bins to solve for h _deb , and corrected for elevation to solve for b ̇ .

Finally, I evaluate the method proposed above against the approach used by Rounce et al. (2021) where b ̇ is estimated using a second-order reaction rate equation, b ̇ ( b ̇ 0 , h d e b ) = b ̇ 0 1 + a 14 b ̇ 0 h d e b . (14) While Rounce et al. (2021) used an energy balance model to solve for b ̇ 0 and a ₁₄, for this study, I simply fit Eq. 14 to measured b ̇ 0 , thus evaluating the expression as if bare ice melt was externally modeled perfectly (under the assumption that my bare ice melt measurements are correct). To evaluate the methods used in this study, the coefficient of determination was found between instantaneous measurements (debris thickness, surface temperature or glacier melt rates) and the respective instantaneous modeled values.

Three ASTER thermal images were acquired under clear sky conditions that were coincident in space and time (2 min, 3 min, and 3 h separation) with the terrestrial based, high resolution thermal images collected within this study. Of the two ASTER images with a separation of minutes, one was acquired during the night (August 28, 2016 23:16 AKDT) and the other during the day (August 30, 2016 13:25 AKDT) about 1.5 days later. These images were processed on-demand by NASA to surface kinetic temperature (AST08) which has a spatial resolution of 90 m and a temperature accuracy and precision of about 1.5 K (Gillespie et al., 1998; Abrams et al., 2002).

Due to the high spatial variability in surface characteristics that can occur over very short length scales within a debris cover (Kraaijenbrink et al., 2018), a correction was applied to the ASTER thermal pixels that were coincident with this study. While acknowledging the ASTER processing algorithm is not linear, the most simplistic approach to applying any correction at all is through linear spatial averaging. Using the high resolution, manually digitized bare glacier ice and ice cliff map used as a validation dataset in Herreid and Pellicciotti (2018), the fraction of bare glacier ice within each ASTER pixel was found. Assuming the temperature value assigned to the ASTER pixel is an area weighted average of the true temperature distribution within the pixel area and assuming the area mapped as bare ice or ice cliff has a surface temperature of 273.15 K, the surface temperature of the remaining debris-covered fraction of the pixel can be computed. This debris area only temperature, T s c o r , can be defined as T s c o r = T s s a t L s a t 2 − ( 273.15 β i c e L s a t 2 ) ( 1 − β i c e ) L s a t 2 (15) where T s s a t (K) is the raw satellite temperature for the whole pixel, L _sat is the length of one edge of a satellite pixel (90 m for ASTER) and β _ice is the fraction of a pixel that can be assumed to have a surface temperature of 273.15 K.

For each thermal image segment with a near-homogeneous debris thickness (Figure 1, Supplementary Figure S1), median and percentile statistics were computed (Figure 4). The time-evolution of the median and spread of temperatures within each segment illustrate when it is most feasible to automate the differentiation of these segments and assign a debris thickness. The ideal setting is when there is near-homogeneous temperatures within a segment (narrow individual ranges, Figure 4F), heterogeneous temperatures across segments (high stacked spread, Figure 4G) and even, un-clustered distribution of the spread (swaths equal in magnitude, Figure 4H).

While satellites can only acquire surface temperature data when the surface is not obstructed by clouds, Figure 4 suggests that the two cloud-free days (August 29 and August 30) were the least optimal of the collected time-series. During these clear sky conditions, the percentile spread in individual segments is high, and the spread between segments is uneven to nonexistent. The warmest two segments (with debris thicknesses of 38 and 8 cm) became nearly isothermic during the peak heat of the day. On the contrary, it appears that overcast days without precipitation offer the most ideal conditions for deriving debris thickness from surface temperature (some rain occurred during this time-series but no major rain events were recorded). While this could be useful for small scale studies using terrestrial or airborne/UAV imaging platforms that can operate under a cloud ceiling, it adds to the challenge of using thermal satellite data to derive debris thickness.

There are reasonable arguments to use thermal data from both night and day-time imagery to derive debris thickness. During the night the near-isothermic state described above would be avoided, and, drawing on the stored heat signal of a thick debris cover, or its absence for a thin debris cover, the surface temperature may scale well with debris thickness. On the other hand, during the day, direct short-wave radiation will warm thick debris cover surfaces while thin layers remain cool, even under direct incoming short-wave radiation. Both approaches seem reasonable and possibly a combination will be most successful. The key circumstances to avoid are a night isothermic state when debris cover of all thicknesses cool off, or during peak heat of a clear day when most rock surfaces are heated to a high temperature and decouple from a predictable relation with the ice below. Within the data acquired in this study, there is evidence of all four states: an optimal night setting, an optimal daytime setting, isothermic nights and isothermic days. Both isothermic states (for all segments at night and the two segments with the thickest debris cover during the day) occur during the clear sky days, August 29 and 30. The remaining measurement days have variable degrees of spread but differentiation seems feasible at any time of day. Based solely on the spread of the full thermal histogram over the observation period (Supplementary Video S2), daytime offers a more broad spectrum of temperature to discretize.

Three empirical equations relating surface temperature to debris thickness were evaluated using the same input data over the same 102 h of observation (Figure 5). Each approach was evaluated in four ways: 1) the number of required model coefficients and their stability through time; 2) the optimal performance of the method where a new set of model coefficients were calibrated for each hourly instance in the 102 h time-series; 3) method performance when holding the model coefficients stable through time as the median values of the model coefficient set derived in 2); and 4) using the model coefficients that were best fit one and 2 h prior.

General results from this analysis show that the methods from Rowan et al. (2021) and Rounce et al. (2021) (Eqs. 10, 11, respectively) easily translate variable surface temperatures through time into stable and correct debris thicknesses when calibrated at every timestep (Figures 5A,B). A notable exception for both methods is an under estimation of debris that is 4 cm thick, both methods returned a near-zero debris thickness for every timestep. The method presented here (Eq. 9), is notably less stable, especially at 8 cm debris thickness, where, at its worst, the debris thickness estimate is off by nearly a factor of 2. A visible advantage of this method is the ability to better estimate the 4 cm thick debris cover. Given these visible shortcomings and advantages, when calibrated at every timestep, all three methods have a high (>0.90) coefficient of determination for every timestep (Figure 5F).

The equations from Rowan et al. (2021) and Rounce et al. (2021) produce a nearly identical curve (Supplementary Video S3, frame examples at the top of Figure 5), yet Eq. 11 from Rounce et al. (2021) requires one additional model coefficient (c ₁₁) that appears to remain constant and does not absorb any temporal variability (Figure 5C). The remaining two model coefficients in Eq. 11 (a ₁₁ and b ₁₁) and the two in Eq. 10 from Rowan et al. (2021) (a ₁₀ and b ₁₀) show considerable variability through time with varying adherence to a diurnal pattern (Figures 5C,D). The two parameter sets are most stable during mid-day on August 4th, which for Eq. 11, means the values are constrained within the ±100% deviation from median bounds, which is still considerable variability. For the rest of the series, both approaches deviate outside of these bounds to absorb the surface temperature variability. While these two methods are very accurate when instantaneously calibrated, their coefficients of determination are rarely above zero when the median model coefficient values were used (Figure 5G, Supplementary Table S1). Rounce et al. (2021) implemented a clever methodology to derive, and optimize a parameter set for each surface temperature instance within a dataset covering all mountain glaciers on Earth. In this case, temporal transferability is not necessarily relevant, but even a short dt between neighboring thermal images might cause a depreciation of results. To explore this, I used each method to solve for debris thickness with the current timestep surface temperature, but used the best fit model coefficients from one and 2 h prior (Figure 5H). Results from this experiment show a dramatic loss of performance, where R ² averaged across the time-series decreased from 0.99 (no lag), 0.53 (−1 h) and 0.37 (−2 h) for Eq. 11 (Rounce et al., 2021); and 0.98 (no lag), 0.54 (−1 h) and 0.41 (−2 h) for Eq. 10 (Rowan et al., 2021).

While the method proposed here, Eq. 9, carries more variability than the other two methods when calibrated at every timestep, it has some key advantages when used in less ideal conditions. Nearly all of the local variability in surface temperature is accounted for by T s * , which means the single model coefficient is simply scaling a physical relation and not needing to vary itself to accommodate different thermal distributions. The parameter stability apparent in Figure 5E carries through both experiments, of a constant median value and −1 and −2 h, with a notably high coefficient of determination. For the hours lag experiment, average R ² remains constant at 0.98 (no lag), 0.98 (−1 h) and 0.98 (−2 h). While these results show more variability from Eq. 9 around the 8 cm debris thickness, it is a quantifiable variability. Because evidence suggest the method is stable, an appropriate error can be assigned and assumed also stable enabling confident and constrained estimates of debris thickness through time.

Following the same framework used to evaluate debris thickness methods above, I propose a method to solve for bare ice and sub-debris melt from debris thickness and surface temperature and compare it to the method from Rounce et al. (2021). The method used by Rounce et al. (2021) (Eq. 14) is not forced by surface temperature, but is a similar simple parameterization of sub-debris melt. For the 102 h time-series, where hourly sub-debris melt is known, measured bare ice melt, b ̇ 0 , and constant, known debris thicknesses were used as input to Eq. 14 and a least squares fit was used to optimized a ₁₄, the single model coefficient, for each timestep. Similarly, the 102 h time-series was used to fit the sub-debris melt method presented in this study for each timestep, but instead of melt data, Eq. 12, solving for debris surface temperature, was forced with thermal camera derived surface temperature over thick debris ( T s * ) and the same constant, known debris thicknesses. The model coefficient, a ₁₂, was fit using thermal camera derived surface temperature for the known debris thicknesses. This solution of debris surface temperature was then used in Eq. 13 to solve for sub-debris melt where a ₁₃ was set to the ratio of b ̇ 0 and T s * , adjusted for each timestep, and b ₁₃ was set to the ratio of median melt rate measured below 38 cm of debris and median melt rate of bare ice (both medians computed from the entire time-series in this study, shown in Figure 6). This means that for this limited 102 h time-series, both methods are prescribed measured bare ice melt rates and will fit the (identical) bare ice validation time-series perfectly. This enables the evaluation of how each method resolves sub-debris melt starting from a perfect solution of bare ice melt, and enables an evaluation of model coefficient stability.

The general results of this experiment are shown in Figure 5I and the Supplementary Video S3 (where three frames are shown at the top of Figure 5. Error within the melt measurements make model validation less robust. I do not have complete confidence that instances of notable deviation, for example, measured sub-debris melt being higher below 8 cm of debris than below 4 cm of debris in the morning of August 4th, are correct. Still, both methods preform reasonably well with average R ² values of 0.85 for Eq. 14 (Rounce et al., 2021) and 0.71 for Eq. 13 (Figure 5L). Repeating the experiment with the median model coefficient, average R ² reduces to 0.77 for Eq. 14 and is 0.50, notably less for Eq. 13 (Figure 5M). Using median model coefficients seems to cause Eq. 13 to fail dramatically during nights when melt is nearly shut off for all debris thicknesses. Interestingly, these nighttime failures are not necessarily caused by a failure of the surface temperature estimates from Eq. 12, which remain stable and accurate for most of the time-series (Figures 5L,M). During the night of August 3rd, R ² of estimated surface temperature crashes, but appears independent of when sub-debris melt is unable to be resolve, for example, the evening of August 1st, afternoon of August 2nd (preceding the surface temperature crash) and the morning of August 11th. Model coefficients appear to remain mostly stable for both methods (Figures 5J,K), which is encouraging for the method presented here because this means that the ratios contained within a ₁₃ and b ₁₃ are stable, physical relations that can now be prescribed as constants, eliminating the need for an explicit solution of b ̇ 0 , which Eq. 14 will always need. By prescribing these ratios as constants, I next apply both approaches to a wider time-series of measurements and evaluate the loss of performance from my method forced with surface temperature and debris thickness alone, against Eq. 14 forced with measured b ̇ 0 .

For the time interval August 1st 2:00 AKDT to August 28th 13:00 AKDT 2016, thermal camera images were not continuously collected, so T s * was set to contact thermistor data collected at the surface of the 38 cm debris cover temperature profile (Sections 3.5, 3.6). Using T s * , debris thickness and median values of a ₁₂, a ₁₃, and b ₁₃ derived above (Supplementary Table S1), Eqs. 12, 13 were used to estimate sub-debris melt below 0, 4, 8 and 38 cm of debris, where sub-debris was simultaneously measured for validation (Figures 6F–I; Supplementary Video S4). Forced with measured b ̇ 0 and the median value of a ₁₄ derived above, Eq. 14 was also used to estimate sub-debris melt, but only below 4, 8, and 38 cm of debris since melt at 0 cm ( b ̇ 0 ) was used as model input. Eq. 14 from Rounce et al. (2021) out performs the method presented here (Figure 6J) with averaged R ² values of 0.50 and 0.38 for the two methods, respectively. While the instantaneous fit to Østrem curve measurements have a relatively low R ², estimated melt summed over the 28 days observation period was close to the measured value (Figure 6K; Table 1). The largest error was modeled melt falling 11 cm below the lower error bound for summed measured melt at 4 cm debris thickness. Summed melt at the lower error bound was 111 cm, indicating the least accurate melt measurement from the method presented in this study had a 10% error. The bare ice estimate was off by only 5 cm, indicating 3% error from the middle measured value (the estimate fell within the error bounds). Some instances where one, or both, methods performed poorly appear to be non-random, correlating with rainfall (Figures 6E,J). The next section considers moisture entering the system and how it could disrupt sub-debris melt estimates.

For a statistical glacier melt model framework emphasizing simplicity, minimal input and transferability, there will always be a trade off in performance relative to more sophisticated approaches with more comprehensive input data. Still, it is worthwhile to consider additional data when available to look for patterns of when the simple method fails and if there might be a simple fix, or an approach to infer method confidence as a function of some variable. From the melt modeling experiment conducted in this study, there are numerous instances where the method fails, yet there is an apparent pattern of failure around rain events (Figures 6E,J). Evaporative cooling of the debris and interstitial water modifying the thermal properties of the debris layer (Collier et al., 2014) are two possible explanations for the phase shift gradient variability visible in Figures 6B,C. While at least one rainy day shows a phase shift gradient that is reduced to near zero, it seems more common that introducing moisture into the debris layer causes more thermal resistance, a more gradual phase shift gradient and an increase in the number of hours it takes diurnal surface warming to reach a depth of 38 cm. During days of heavy rain, apparent thermal diffusivity formulated to include all energy transfer mechanisms, κ _a , spikes while apparent thermal diffusivity considering only conduction, κ _c remains essentially static (Figure 6D). Drawing on this observation, the difference between κ _a and κ _c appears to have the potential to be a proxy measure of precipitation, where relative precipitation rate may even be inferred (Figure 5E). Emphasizing again a methodology favoring simplicity and minimal field sensor deployment, it might be advantageous to establish a network of surface and sub-debris thermistor strings rather than one or two costly and cumbersome weather stations.

While I provide evidence in Section 4.1 suggesting clear days are the least favorable for deriving debris thickness from thermal data (Figure 4, August 28 to August 30), it was during these days that two ASTER scenes were acquired concurrently within minutes (a third within hours) of a high resolution, field-based thermal image. The timing of the ASTER images are shown in Figure 4F (red dashed vertical lines) and show that the August 28th acquisition narrowly missed the nighttime isothermic state, and captured the August 30th warming limb of the diurnal cycle. The sub-pixel bare ice correction described in Section 3.7 was applied to each image (Figure 7) and a comparison between ASTER and field-based thermal camera surface temperature data is shown in Figure 8 where the ASTER data is draped over a local DEM [from Herreid and Pellicciotti (2018)] and obliquely angled to mimic the field of view of the field-based thermal image. While acquired under less favorable conditions, the relation between the ASTER pixels contained within the field of view of the field-based thermal camera image agree with R ² values of 0.70 and 0.78 for the nighttime August 28 and daytime August 30 images, respectively. The near isothermic state during the night is apparent, yet a thermal delineation can still be made between the near-field 38 cm thick debris cover and the adjacent 8 cm thick debris (segment boundaries are shown in Figure 1). The ability of these two ASTER images to replicate high resolution field-based data during less than ideal conditions is promising for satellite based methodologies to derive debris thickness.