2.1. Ice-penetrating radar

2.1.2. Groundhog radar

Developed at the University of Arizona, Groundhog is a snowmobile-towed, ice-penetrating radar. On the transmit end, Groundhog employs a Kentech impulse generator which is operated at a pulse repetition rate of either 1 or 5 kHz. An Ettus X310 software defined radio was used to digitize radar echoes at 20 MHz. The transmitter and receiver were each mounted within sleds at the center of two resistively loaded 50 m half-dipole antennas. The transmit and receive assemblies are separated by ~15 m, making the total system ~215 m in length (Fig. 2). Each antenna element is encased within hollow braid polyethylene rope such that the rope, and not the resistively loaded antenna element, takes on the brunt of the tension when the system is in tow. To record system positioning information, both the transmit and receive sled were instrumented with Emlid Reach RS2+ GNSS receivers logging at 1 s intervals. Kinematic precise point positioning solutions were generated through the Canadian Spatial Reference System Precise Point Positioning web service (Tétreault and others, Reference Tétreault, Kouba, Héroux and Legree2005), resulting in positioning accurate to the decimeter level.

Figure 2. Groundhog ice-penetrating radar common-offset setup. The receiver (rx) and transmitter (tx) each sat in sleds (red) with Emlid GNSS receivers and were connected to two half-wavelength dipole antenna elements, each measuring 50 m.

Data acquisition was controlled by an operator riding within the receiver sled. Groundhog was towed across the amphitheater and the Great Gorge of Ruth Glacier between 1st May and 9th May 2022, acquiring a total of 143 line-km of common offset radar profile data. We also acquired several profiles by leaving the transmitter stationary and conducting a moveout survey. None of these provided useful data, so they will not be discussed here.

Groundhog data were uniformly processed by removing the mean trace from each profile and applying a 500 kHz–9 MHz bandpass filter. Data were migrated following the Stolt method with an assumed wavespeed in ice of 169 m μs⁻¹ (Stolt, Reference Stolt1978), collapsing diffractions and repositioning dipping reflectors to their true subsurface locations. Profiles were vertically shifted along the fast-time axis to account for the delay between the time of transmission and the airwave-triggered start of each record.

2.1.3. Radar-derived ice thickness and bed elevation

Both ARES and Groundhog data were investigated for the presence of glacier bed returns using the open-source Radar Analysis Graphical Utility software (RAGU; Tober and Christoffersen, Reference Tober and Christoffersen2020). Horizontally continuous and high-amplitude radar returns were interpreted as echoes from the glacier bed (following verification through clutter simulations for airborne data; see Section 2.1.1). We manually digitize (‘pick’) the arrival time of the glacier bed return in ARES and Groundhog data (Figs 3, 4). In ARES data, which has been pulse-compressed to contain approximately zero-phase wavelets, this corresponds to the peak of the return, while in Groundhog data, which contain minimum phase (impulse) wavelets, this is the ‘first-break’ (onset of the first arrival; Wilson, Reference Wilson2021). Digitized bed returns provide point measurements of the two-way travel time to the glacier bed.

Figure 3. ARES example radar profile. (a) Profile A–A′ location over the lowermost ~20 km of Ruth Glacier. Context of (a) shown by red box atop white RGI glacier outline at upper right. (b) ARES radar sounding profile A–A′. Top panel shows uninterpreted profile, second panel shows clutter simulation, third panel shows interpreted profile with altimetry-derived surface elevation in blue and manually digitized glacier bed in orange, and bottom panel shows cross-sectional elevation profile with the glacier surface in blue and the bed in orange.

Figure 4. Groundhog example radar profile. (a) Profile B–B′ location in the Don Sheldon Amphitheater. Context of (a) shown by red box atop white RGI glacier outline at upper right. (b) Groundhog radar sounding profile B–B′. Top panel shows uninterpreted profile, second panel shows interpreted profile with the manually digitized glacier bed in orange and bottom panel shows cross-sectional elevation profile with the glacier surface in blue and the bed in orange.

Point-measurement spacing along radar profiles is a function of the radar pulse repetition rate (1–5 kHz), the system groundspeed, and the trace stacking interval, and is typically ~5 m for both ARES and Groundhog data. Following Armstrong and others (Reference Armstrong2022), the two-way travel time in ice t for a given bed return is converted to ice depth h:

(1)$$h = \left[{{\left({\displaystyle{v \over 2}\left({t + \displaystyle{{x_{{\rm sep}}} \over c}} \right)} \right)}^2-{\left({\displaystyle{{x_{{\rm sep}}} \over 2}} \right)}^2} \right]^{1/2}$$

where v = 169 m μs⁻¹ is the radio wave speed in temperate ice (relative permittivity of 3.15; Evans, Reference Evans1965), c = 300 m μs⁻¹ is the approximate radio wave speed in air, and x _sep is the separation distance between the centers of the receiver and transmitter antenna elements. The transmitter–receiver separation distance for Groundhog data ranged from 105 to 120 m during the survey, varying as a function of antenna stretching and the system meandering across the surface. We use the precise separation distance provided by post-processed transmitter and receiver sled GNSS positions to calculate the depths of digitized bed returns along a given profile using Eqn (1). The midpoint between the transmit and receive sleds was assigned as the ground positioning for digitized Groundhog bed returns. ARES uses a single end-fed antenna, so the antenna separation distance x _sep in Eqn (1) is zero. Accounting for the radio wave travel time from the aircraft to the glacier surface using coincident OIB-AK laser altimetry data (Larsen, Reference Larsen2020), the ice thickness is

(2)$$h = \left({t-\displaystyle{{2 \times z} \over c}} \right)\displaystyle{v \over 2}\;$$

where z is the flight altitude above ground level.

The quantization interval of bed echo interpretations, as well as the radar system timing and positioning together impact the precision of bed echo interpretations. Following Lapazaran and others (Reference Lapazaran, Otero, Martín-Español and Navarro2016), we assess the combined effect of these factors by analyzing the disagreement in ice thickness at profile intersections. Tober and others (Reference Tober2023) demonstrated from extensive crossover analysis that the precision of radar-derived bed returns for OIB-AK airborne radar data is on the order of ~10 m. We similarly assess the precision of Groundhog radar-derived ice thickness measurements in Section 3.1. Laboratory analyses have demonstrated that the dielectric permittivity in temperate ice can range from 3.1 to 3.2 (Evans, Reference Evans1965), translating to wavespeed range of 168–170 m μs⁻¹. We therefore acknowledge an additional uncertainty in the computed ice thickness of 2%.

2.2. Estimating depth of the Great Gorge

As radar-derived measurements of ice thickness were not achievable within the Great Gorge via airborne or surface-based surveys due to the glacier's geometry (Section 4.2), we employed a mass conservation approach to estimate the glacier's thickness using radar-derived measurements of ice thickness combined with satellite-derived glacier surface velocities. Following the principle of mass conservation, the climatic-basal mass balance rate $\dot{b}$ is balanced by the ice flux divergence $\nabla \cdot \vec{q}\;$ and the resulting change in surface elevation ∂h/∂t:

(3)$$\displaystyle{{\partial h} \over {\partial t}} = \dot{b}-\nabla \cdot \vec{q}.$$

The climatic-basal mass balance rate is typically dominated by, and can thus be replaced with, the surface mass balance rate $\dot{b}_{{\rm sfc}}$ (Cuffey and Paterson, Reference Cuffey and Paterson2010), which in contrast to the climatic-basal balance rate is possible to determine through traditional field measurements.

Applying Gauss's theorem and integrating Eqn (3) across the surface area S between two transverse profiles (‘gates’), we arrive at

(4)$$q_{{\rm out}}-q_{{\rm in}} = \mathop \smallint \limits_S^{} \left({{\dot{b}}_{{\rm sfc}}-\displaystyle{{\partial h} \over {\partial t}}} \right){\rm d}S$$

where q _in is the ice influx from the up-glacier gate, and q _out is the ice efflux from the down-glacier gate.

We calculate the ice flux q through segment (dx, dy) of a given gate as

(5)$$q = \bar{h} \times \gamma ( {\overline {vy} \times {\rm d}x-\overline {vx} \times {\rm d}y} ) $$

where between consecutive points $\bar{h}$ is the average ice thickness dx is the x-distance, dy is the y-distance, $\overline {vx}$ is the average x-component of the surface velocity, $\overline {vy}$ is the average y-component of the surface velocity, and γ is the factor relating the observable surface velocity to the depth-averaged ice velocity.

For a glacier under simple shear with no sliding, it can be shown that γ =  (n + 1)/(n + 2), where n is the exponent in the flow law of ice (e.g. Cuffey and Paterson, Reference Cuffey and Paterson2010). For the commonly used flow law exponent of n = 3, this would result in γ = 0.8. This should be seen as a minimum estimate, since a temperate glacier such as the Ruth is expected to have a significant component of basal sliding and γ is more realistically expected to fall between 0.9 and 1 (e.g. Cuffey and Paterson, Reference Cuffey and Paterson2010). Given a sensible value for γ, summing the ice flux through each segment from Eqn (5), we arrive at the total flux through a given gate.

The glacier's mean 1985–2018 surface velocity field was generated using auto-RIFT (Gardner and others, Reference Gardner2018) and provided by the NASA MEaSUREs ITS_LIVE project (Gardner and others, Reference Gardner, Fahnestock and Scambos2023; hereafter referred to as ITS_LIVE). Ice thickness data provided by the Groundhog ice-penetrating radar across the amphitheater were extrapolated to the valley edges and combined with ITS_LIVE surface velocities to constrain the ice influx to the Great Gorge (Fig. S1). To constrain the ice flux down-glacier, we take measurements of ice thickness provided by ARES ~17–20 km up-glacier from the terminus, along a profile which follows the glacier's flow direction and construct a flux gate perpendicular to ice flow. Surface speeds through this gate indicated the potential for an asymmetric cross-sectional bed shape, so we estimate the ice flux there by considering both a symmetric (parabolic) and an asymmetric bed (Fig. S2). Although another set of ARES ice thickness measurements were obtained between ~4 and 11 km up-glacier from the terminus (Fig. 1), these measurements are from debris-covered ice and are therefore less useful in constraining the surface mass balance rate.

Provided the average annual change in surface elevation for the period from 2000 to 2019 by Hugonnet and others (Reference Hugonnet2021; Fig. S3), and with the ice thickness constrained across two flux gates, the only unknown variables remaining in Eqns (4) and (5) are the surface mass balance rate and the γ-factor.

We assume that the surface mass balance rate increases linearly with altitude (e.g. Beusekom and others, Reference Beusekom, O'Nell, March, Sass and Cox2010), and can thus be expressed as

(6)$$\dot{b}_{{\rm sfc}} = \displaystyle{{\;( {z_{{\rm sfc}}-z_{{\rm ela}}} ) } \over {\rho _{{\rm ice}}}}\displaystyle{{{\rm d}{\dot{b}}_{{\rm sfc}}} \over {{\rm d}z}}$$

where z _sfc is the surface elevation provided by a 2012 Interferometric Synthetic Aperture Radar (IFSAR) DEM, z _ela is the equilibrium line altitude where the surface mass balance rate is zero, ρ _ice is the density of ice, and ${\rm d}\dot{b}_{{\rm sfc}}/{\rm d}z$ is the surface mass balance gradient.

While no known records of surface mass balance exist for Ruth Glacier within the published literature, measurements for Gulkana Glacier – one of five USGS benchmark glaciers in Alaska – date back to 1966 (March and O'Neel, Reference March and O'Neel2011). Gulkana Glacier is located on the southern flank of the eastern Alaska Range, at a latitude and continental climatic setting similar to that of Ruth Glacier. While we estimate the equilibrium line altitude for Ruth Glacier based on of summer snowlines observations to be 1500–1550 m (Fig. S4) – 200–300 higher than that reported by March and O'Neel (Reference March and O'Neel2011) for Gulkana – we suspect that the annual mass-balance gradients of both glaciers are comparable. Applying a piecewise-linear fit to the annual point balance data reported by O'Neel and others (Reference O'Neel2019, their Fig. 3) provides a surface mass balance gradient of ~10 mm w.e. m⁻¹ a⁻¹ within the glacier's ablation zone, which we use as a reference point for the expected surface mass balance gradient of Ruth.

In common with many inverse problems in glaciology (e.g. Brinkerhoff and others, Reference Brinkerhoff, Aschwanden and Truffer2016; Rounce and others, Reference Rounce2020), we would expect that our mass conservation model is overparameterized: despite having some a priori information on the expected range of our unknown model parameters (γ, z _ela, ${\rm d}\dot{b}_{{\rm sfc}}/{\rm d}z$), an infinite number of parameter combinations will likely conserve mass between our two radar-constrained flux gates. To locate all sensible parameter combinations, we perform a grid search and calculate the model misfit at our downstream flux gate relative to the median ice flux estimated through the gate from OIB-AK radar-derived ice thickness measurements, considering both a symmetric and asymmetric bed shape (Table 1). While we expect Ruth Glacier to have a significant component of basal sliding, we consider all sliding cases between 0.825 and 0.975 (0.025 step size), given that a γ-factor of 0.8 would indicate no sliding, and a value of 1 would indicate a fully decoupled glacier bed and thus ‘plug-like’ flow (Nye, Reference Nye1951). While we estimate the equilibrium line altitude to be ~1550 m based on summer snowline observations, we consider all values between 1450 and 1650 m in our grid search (5 m step size). All surface mass balance gradients between 1 and 15 mm w.e. m⁻¹ a⁻¹ (1 mm w.e. m⁻¹ a⁻¹ step size) are considered, encapsulating the reference value of 10 mm w.e. m⁻¹ a⁻¹ found at Gulkana Glacier by O'Neel and others (Reference O'Neel2019). We calculate the ice flux at the down-glacier extent of the model domain for all parameter permutations and consider all parameter sets which lead to a relative ice flux error of <10% as sensible (Fig. 5). A misfit threshold of 10% relative to the median value constrained by OIB-AK ice-penetrating radar was chosen such that both end-member ice fluxes for a symmetric and asymmetric bed were captured along with some margin of uncertainty. Additional analysis was performed to assess the resulting discrepancy in modeled ice thickness for a greater misfit threshold (Section 3.2).

Table 1. Mass conservation model grid search a priori estimates

Figure 5. Mass conservation model a priori parameter value grid search. Each panel (a–g) represents a different γ-factor, indicated at the lower left. Contours show the misfit in the modeled ice flux at the down-glacier extent of our model domain, relative to the ice flux estimated there by OIB-AK radar-derived ice thickness measurements. Only relative misfits less than a factor of 0.5 are shown.

For each sensible parameter set, we take the ice thickness at our upstream flux gate, the surface mass balance rate, and the ITS_LIVE surface velocity, and calculate the ice flux across a series of transects provided by the Open Global Glacier Model (OGGM; Maussion and others, Reference Maussion2019) between our up-glacier and down-glacier radar-constrained gates (Fig. S5). Assuming a parabolic glacier cross section for each gate (Nye, Reference Nye1965), we then invert for the bed position while imposing the constraint that the thickness reaches zero at the edges of each gate. Ice thickness along each gate is determined by subtracting the inverted bed position from the 2012 surface elevations provided by IFSAR. Bed elevation and ice thickness across the model domain were linearly interpolated at a spatial resolution of 100 m to provide gridded results. This inversion approach yields a mean ice thickness as well as a marginal standard deviation. We choose to present two times the marginal standard deviation as our model uncertainty, which implies an ~95% probability that the true ice thickness lies within these bounds.