1. Introduction
Quantifying snow water equivalent (SWE) is crucial for assessing the annual mass balance of glaciers, water resources and potential glacier contributions to sea-level rise. Ground-penetrating radar (GPR) offers a valuable quantitative, nondestructive method to estimate SWE along profiles tens to hundreds of kilometers long within a single field season, accounting for variability of snow distribution across larger areas (e.g., Bradford and others, Reference Bradford, Nichols, Mikesell and Harper2009; McGrath and others, Reference McGrath2015; McGrath and others, Reference McGrath, Sass, O’Neel, McNeil, Candela, Baker and Marshall2018). Often, studies conduct surveys early in the year (spring), avoiding the impact of liquid water content (LWC), simplifying assumptions of radar wave velocity when determining snow depth. However, research suggests that the sub-Arctic and Arctic are already rapidly warming (Winski and others, Reference Winski2018; Ballinger and others, Reference Ballinger2023) and are projected to continue doing so, leading to earlier spring snowmelt onset (Zheng and others, Reference Zheng, Jia and Xu2022), more incidences of rain-on-snow events (Pan and others, Reference Pan, Kirchner, Kimball, Kim and Du2018), and greater contributions from glaciers to global sea-level rise (Oerlemans and others, Reference Oerlemans2005; Rounce and others, Reference Rounce2023). Given this trend, assessments of wet snowpacks for SWE will likely become more common. Therefore, our aim is to determine the variability of wet snowpack properties and quantify their impact on SWE estimates using GPR and ground-truth methods.
Most data related to our understanding of mass balance and SWE on Alaskan glaciers come from smaller U.S. Geological Survey (USGS) Benchmark glaciers (areas < 20 km2) (McGrath and others, Reference McGrath, Sass, O’Neel, McNeil, Candela, Baker and Marshall2018; O’Neel and others, Reference O’Neel2019), with few results from larger glaciers (McGrath and others, Reference McGrath2015). McGrath and others (Reference McGrath2015) used GPR surveys across seven glaciers and five climate regions and found that SWE across Alaskan glaciers is highly variable and heavily influenced by wind-shading, elevation, and storm track trajectory. Further, they find that accumulation rates derived from radar exceed those estimated by glaciological methods and that there is significant SWE variability along those centerlines. This included a small survey (several kilometers) on Taku Glacier, suggesting that there is limited spatial variability in end-of-winter SWE on Taku Glacier and greater variability on other Alaskan glaciers in their study.
Determining snow depth and SWE with GPR requires knowledge of the radar wave velocity traveling through the snowpack. Typically, this velocity is assumed to be spatially or temporally homogeneous. Many studies have addressed this limitation by applying a density-dependent bulk radar wave velocity across the study area (Gascon and others, Reference Gascon, Sharp, Burgess, Bezeau and Bush2013; Miège and others, Reference Miège2016; Sold and others, Reference Sold, Huss, Eichler, Schwikowski and Hoelzle2015; McGrath and others, Reference McGrath, Sass, O’Neel, McNeil, Candela, Baker and Marshall2018). However, since radar wave velocity is affected by snowpack density and strongly affected by LWC, both properties that can vary in space and time, this assumption can lead to uncertainty in SWE estimates. Understanding how density and LWC vary spatially is therefore crucial to interpretations of snow depth and calculating SWE from GPR observations. Several recent studies focused on evaluating the effects of density and LWC on radar wave velocity. For example, St. Clair and Holbrook (Reference St. Clair and Holbrook2017) used migration velocity analysis with common-offset radar to determine the most appropriate velocity for a dry snowpack in Wyoming. Recent work in seasonal snowpacks used GPR in combination with destructive (e.g. snow pits and dye tracing) and nondestructive methods (e.g. terrestrial Light Detection and Randing, airborne imagery, and Structure from Motion) for assessing and monitoring spatial variations in seasonal snowpack depth, density, SWE, and volumetric LWC throughout the springtime (Webb and others, Reference Webb, Jennings, Fend and Molotch2018; Webb and others, Reference Webb, Wigmore, Jennings, Fend and Molotch2019; Bonnell and others, Reference Bonnell, McGrath, Williams, Webb, Fassnacht and Marshall2021; McGrath and others, Reference McGrath, Bonnell, Zeller, Olsen-Mikitowicz, Bump, Webb and Marshall2022; Webb and others, Reference Webb, Musselman, Ciafone, Hale and Molotch2022; Bonnell and others, Reference Bonnell2023). These studies find that LWC can vary spatially (Webb and others, Reference Webb, Jennings, Fend and Molotch2018; Bonnell and others, Reference Bonnell2023), diurnally (Bonnell and others, Reference Bonnell, McGrath, Williams, Webb, Fassnacht and Marshall2021), and seasonally (Schmid and others, Reference Schmid, Koch, Heilig, Prasch, Eisen, Mauser and Schweiser2015). Other approaches have implemented multiple antenna configurations to determine radar wave velocities of snow and firn layers. Examples include common mid-point surveys on annual accumulation (McGrath and others, Reference McGrath2015) and multi-offset radar on seasonal snowpacks (Greissinger and others, Reference Greissinger, Mohr and Jonas2018) and on polar firn (Meehan and others, Reference Meehan2021). Our work applies an approach most similar to Bonnell and others (Reference Bonnell, McGrath, Williams, Webb, Fassnacht and Marshall2021) using in situ snow depths to calculate snow LWC.
With continued warming, we expect the window of dry snow to shorten and further work is necessary to better constrain LWC and SWE variations in wet snowpacks to expand the time to capture winter accumulation across large spatial areas. Here, we utilize 400 MHz common-offset GPR surveys in conjunction with in situ snow depth and density measurements to further constrain and assess the impact of wet snow properties on SWE estimates across the Juneau Icefield (JIF) in Southeast Alaska (Fig. 1; Supplementary Fig. S1). We calculate bulk radar wave velocity in 2012 and 2021 at a total of 13 locations with GPR surveys and ground-truth snow depths from snow pits and cores. We then combine these calculations with dry (spring) snow density measurements to quantify the variability in LWC across the study site. Density measurements collected from 31 summer snow pits on JIF in 2012 and 2021 are combined with snow depth interpretations from repeat 400 MHz common-offset GPR surveys to calculate bounds of the range of SWE across JIF in July 2012 and 2021. We estimate variability in relative permittivity and LWC within the snowpack and the resulting uncertainties in SWE estimates, constrained by in situ snow density and depth measurements.
Southern portion of the Juneau Icefield with individual glacier outlines (colored) from the Randolph Glacier Inventory (RGI 7.0 Consortium, 2023). GPR transects from 2012 and 2021 (dotted grey, solid black) are shown with ground-truth locations: snow pits (squares) and cores (dotted circles), superimposed with surface elevation contours at 250 m spacing. Inset map indicates the study area (red) relative to Alaska with base imagery from ESRI. Source DEMs courtesy of the Polar Geospatial Center.
2. Study site
2.1. Juneau Icefield (JIF), Alaska
JIF (Fig. 1), one of Alaska’s largest icefields, is currently losing mass with increased warming (McNeil and others, Reference McNeil, O’Neel, Loso, Pelto, Sass, Baker and Campbell2020), averaging 5.9 ± 0.8 km3 a−1 in volume loss between 2010 and 2020 (Davies and others, Reference Davies2024). Taku Glacier (725 km2), traditionally named T’aakú Kwáan Sít’i (the T’aakú peoples’ glacier; Thornton, Reference Thornton2010), is the main outlet drainage of the southern portion of JIF (McNeil and others, Reference McNeil, Amundson, O’Neel, Motyka, Sass, Truffer, Zechmann and Campbell2021). Since 1946, rates of annual mass loss have increased in correlation with increasing mean summer temperatures (McNeil and others, Reference McNeil, O’Neel, Loso, Pelto, Sass, Baker and Campbell2020). ERA5-Land 2 m temperature data (Muñoz Sabater, Reference Muñoz Sabater2019) indicate a warming of 0.28°C per decade over the JIF region (58–59.3°N, 133.2–135°W) between 1950 and 2024 and with average yearly temperatures more frequently approaching 0°C (Supplementary S.2.). JIF receives between 2.5 and 4.4 m w.e. of winter (October–March) precipitation each year (Roth and others, Reference Roth, Hock, Schuler, Bieniek, Pelto and Aschwanden2018) and experiences prevalent melt during summer months (JJA) producing supraglacial melt ponds and streams in both the accumulation and ablation zones. The established field stations and annual field campaigns of the Juneau Icefield Research Program (JIRP) provide ideal access and logistical support to enable extensive GPR surveys and ground-truth efforts during peak warming (and increased snow water content) months.
3. Methods
3.1. Ground-penetrating radar
To assess the impact of LWC on SWE from GPR, we analyzed approximately 150 km of GPR surveys collected across JIF in July 2012 and in July 2021, following the 2012 GPS tracks for repeatability (Campbell, Reference Campbell2024). Surveys covered eight catchment basins and elevations from ∼1000 to 1800 m, providing data to encompass the spatial variability of snowpack depth, LWC, and SWE. We conducted these surveys by towing a Geophysical Survey Systems Incorporated (GSSI) control unit (SIR-3000 in 2012; SIR-4000 in 2021) and Garmin GPSMap handheld units for georeferencing (see Supplementary S.3. for detailed processing) coupled with a GSSI model 5103 400 MHz antenna on a plastic sled at 3–7 km h−1 via snow machine. Surveys were conducted along the centerline of several branches of Taku Glacier and along three cross-sections (Fig. 1).
GPR data were collected with a time window of 250–400 ns, 2048 samples per trace, 24 traces per second, and were not stacked during collection. The sampling rate and towing speeds of 3–7 km h−1 yield traces at 3.5–8.1 cm intervals across the profiles. After collection, all GPR profiles were time-zero corrected, distance normalized, and stacked 10- to 30-fold using GSSI’s proprietary software, RADAN v.7. Some profiles were filtered with low- and high-pass filters (800 MHz and 100 MHz, respectively) if noise was apparent in the radargrams. However, the data quality was generally high, so no other filtering or processing was required. We identified the annual accumulation–firn boundary as the most spatially continuous horizontal reflection horizon across both the ablation and accumulation zones (Fig. 2) and converted the GPR two-way travel time (TWTT) data to annual snow accumulation depth (z) with:
\begin{equation}z = \,\frac{{TWTT\, \cdot \,v}}{2}\end{equation}where the bulk average velocity v of the radar wave through the snowpack is given by
\begin{equation}v = \,\frac{c}{{\sqrt {{{{\varepsilon }}_s}} }}\end{equation}where c is the speed of light (0.3 m ns−1) and ԑ s is the relative permittivity of snow.
(a) Location of 2012 GPR profiles and snow pit TKG3, with a close-up indicating the location of TKG3 relative to the same-day radar survey. (b) Radargram along the centerline of Matthes and Taku Glaciers with the profile extent in (a) denoted by A-A’. (c) Close-up of radargram at the extent of the red box in (b) at closest point to TKG3 (red arrow).
3.2. Ground-truth snow depth and relative permittivity
In this study, 
${\varepsilon _{\text{s}}}$ is determined by
\begin{equation}{{{\varepsilon }}_s} = {\left( {\frac{{c \cdot {\text{TWT}}{{\text{T}}_{gt}}}}{{2{z_{gt}}}}} \right)^2}\end{equation} using ground-truth measurements of depth (zgt), measured by snow pit and snow core depths located within 25 m of the GPR survey and TWTT of the accumulation layer at the ground-truth site (TWTTgt). The latter are derived from the stacked profile, thus utilizing the average of the closest 10–30 original scans, where GPR transects passed at a distance of 3–24 m from each ground-truth site to both capture the local accumulation and also minimize potential reflections from the snow pit wall itself (Fig. 2). Sensitivity tests were performed on these single TWTTgt picks compared to the average picks over 50 m and the average of 25 traces nearest to the ground-truth site to validate our use of TWTTgt for 
${{{\varepsilon }}_s}$ calculations (Supplementary S.4., Table S1).
The snow pits used for ground-truth snow depth were selected from a dataset collected by JIRP during summer mass balance surveys (Campbell, Reference Campbell2024). Additional ground-truth snow depths were collected in 2021 with a Kovacs Mark II Coring System. GPR surveys and ground-truth snow depth measurements were obtained concurrently whenever possible, but, due to the time required to dig snow pits, this was not possible for three sites in 2012. For these cases, we applied a surface correction using ablation stake measurements (Table 1; see Supplementary S.5.). Additionally, we assessed the relationship between 
${{{\varepsilon }}_s}$ and lapse-rate-derived air temperature (see Supplementary S.6.).
Dates and types of ground-truth data and GPR surveys.
Surface correction was applied to the sites in 2012 where GPR and GT measurements were not collected on the same day.
3.3. Liquid water content
We calculated the volumetric LWC of the snowpack at each ground-truth location using an empirical relationship between dry snow density (ρ d), volumetric LWC (W v), and relative permittivity of dry, wet, and measured snow (
$\varepsilon $d, 
$\varepsilon $w, and 
$\varepsilon $s) (Tuiri and others, Reference Tiuri, Sihvola, Nyfors and Hallikaiken1984; Sihvola and Tiuri, Reference Sihvola and Tiuri1986):
\begin{equation}{\varepsilon _{\text{s}}} = \left( {0.10{W_{\text{v}}} + 0.80W_{\text{v}}^2} \right){\varepsilon _{\text{w}}} + {\varepsilon _{\text{d}}}\end{equation}where the permittivity of dry snow is defined by Tuiri and others (Reference Tiuri, Sihvola, Nyfors and Hallikaiken1984) as follows:
\begin{equation}{\varepsilon _{\text{d}}} = 1 + 1.7{\rho _{\text{d}}} + 0.7{\rho _{\text{d}}}^2\end{equation} For all LWC calculations, we used a single value for ρ d of 0.444 g cm−3, estimated by averaging spring snow densities collected across JIF in 2013–2015 (ranging between 0.400 g cm−3 and 0.550 g cm−3; McNeil and others, Reference McNeil, Sass, Baker and O’Neel2017; O’Neel and others, Reference O’Neel2017), and used a value of 88 for 
$\varepsilon $w for water at the freezing point. Although these densities were collected outside of the study years, we assume these measurements are representative of dry snow density for this snowpack because of the documented consistency in the summer wet snow density across JIF (LaChapelle, Reference LaChapelle1954; Pelto and Miller, Reference Pelto and Miller1990; McNeil and others, Reference McNeil, O’Neel, Loso, Pelto, Sass, Baker and Campbell2020).
3.4. Snow water equivalent
We calculated SWE across JIF using GPR and all available density measurements from JIRP’s mass balance program in 2012 and 2021 (Pelto and others, Reference Pelto, Kavanaugh and McNeil2013; McNeil and others, Reference McNeil, O’Neel, Loso, Pelto, Sass, Baker and Campbell2020), including data from nearby or adjoining glaciers such as Lemon Creek Glacier and Tulsequah Glacier (Fig. 1). These snow depth and wet snow density measurements were from a total of 18 and 13 snow pits in 2012 and 2021, respectively; typically, snow pit depths reached between 300 and 500 cm. Density measurements were sampled from the north-facing wall of each pit at 10 cm intervals using a 500 cm3 sampling cylinder. Any observable ice lenses or ice pipes over 0.2 cm through the snowpack were recorded. For each study year, density data from all snow pits were aggregated and averaged within each 10 cm depth interval to yield the average and SD of wet density (
$\rho $w) at each depth interval across JIF, and an exponential equation was fit to all the density data for each study year (
$\rho $w(z)) (Supplementary Fig. S4).
Because using snow depth-dependent density calculations reduce errors in calculating SWE (Lundberg and others, Reference Lundberg and Richardson-Näslund2006), we calculated SWE along the JIF profiles by integrating the exponential relation derived above, 
$\rho $w(z), from the surface (z = 0) to the snowpack depth (z = d) to capture the density differences throughout the snow column with:
\begin{equation}SWE = \,\mathop {\smallint }\limits_0^d {\rho _w}\left( z \right)dz\end{equation} For the two study years, SWE was calculated using snow depths from mean (
$\overline {{\varepsilon _s}} \,$± 1 SD), minimum (
${\varepsilon _{s,min}}$) and maximum (
${\varepsilon _{s,max}}$) relative permittivity in each year to provide bounds on SWE for the range of observed conditions.
4. Results
4.1. Relative permittivity and LWC
Relative permittivity and LWC (% volume) estimates show variability across the study region, with no consistent trend between 2012 and 2021 (Figs. 3, 4). Using the 13 ground-truth sites, estimates of 
${\varepsilon _s}$ range from 2.06 to 2.62 in 2012 (
$\overline {{\varepsilon _s}} = \,$ 2.48 ± 0.20) and from 2.11 to 5.11 in 2021 (
$\overline {{\varepsilon _s}} = \,$ 3.56 ± 1.07) (Fig. 3, Table 1). Uncertainty in relative permittivity is 0.106, calculated using mean TWTTgt (50.93 ns) and mean snow depth (447 cm) (see Supplementary S.8.). Using the calculated 
${\varepsilon _s}$, LWC estimates range between 1.7% and 5.7% in 2012 (avg. 4.8%), and from 2.1% to 16% in 2021 (avg. 9.6%). While higher elevations exhibited higher relative permittivity values, no dominant spatial trend was observed (Fig. 3). Additionally, LWC values were not consistent between years, nor is there a consistent trend at the same locations (Fig. 4; Supplementary Fig. S5).
Calculated relative permittivity (εs) values in (a) 2012 and (b) 2021 for each snow pit used for ground-truth. Basemap is Landsat 8 Band 8 imagery from 24 April 2021.
Calculated liquid water content (Wv, %) versus elevation (m) for 2012 (orange) and 2021 (blue) with ground truth (GT) points distinguished with different symbols. LWC was calculated using εs (Fig. 3) and the mean spring snow density (0.444 g cm-3).
4.2. GPR-derived snow depths
The average TWTT for the accumulation layer across all radar transects was 51.1 (±16.3) ns in 2012 and 44.5 (±17.8) ns in 2021. For 2012 depth calculations, we used 
$\overline {{\varepsilon _s}} $ of 2.48 (
${\rho _d}$ = 0.444 g cm−3, W v = 4.8%), resulting in an average annual accumulation depth of 524 cm ± 153 cm. In 2021, the 
$\overline {{\varepsilon _s}} $ of 3.56 (
${\rho _d}$ = 0.444 g cm−3, W v = 9.6%) applied across the dataset resulted in an average annual accumulation depth of 351 ± 138 cm.
4.3. Snow pit density and snow water equivalence estimates
To calculate SWE, we established the average wet snow density curves using all available snow pits collected in 2012 and 2021, consisting of 18 snow pits in 2012 and 13 snow pits in 2021 (Fig. 1; Supplementary Fig. S4). Snow pit densities were relatively consistent between years, averaging 0.556 ± 0.030 g cm−3 in 2012 (avg. depth 552 ± 94 cm) and 0.562 ± 0.025 g cm−3 in 2021 (avg. depth 422 ± 88 cm). Snow density follows the expected trend of increasing density with depth (Supplementary Fig. S4). Equations describing the relationship between density and depth for each year are ρ w,2012(z) = 0.403z 0.06 (R 2 = 0.86) and ρ w, 2021(z) = 0.450z 0.04 (R 2 = 0.66).
Using GPR-derived snow depths and density data, we calculated average SWE to be 286 ± 89 cm w.e. (
$\overline {{\varepsilon _s}} $ = 2.48) in 2012 and the average SWE in 2021 to be 193 ± 80 cm w.e. (
$\overline {{\varepsilon _s}} $ = 3.56) (Fig. 5), representing a −32% change between 2012 and 2021 (see Supplementary S.7. for ice lens considerations). However, there is high spatial variability in 
${\varepsilon _s}$ in 2021, which we infer to be due to the variability in LWC. We therefore calculated the range of SWE possible in each year based upon the minimum and maximum calculated relative permittivity for each year (
${\varepsilon _{s,min}},\,\,{\varepsilon _{s,max}}$) and found a potential 13% difference in mean SWE in 2012 and a 45% difference in mean SWE in 2021 (Table 2). Acknowledging those are the maximum differences, we also calculated the potential difference in mean SWE for 
$\overline {{\varepsilon _s}} $ ± 1 SD to be 8% and 33% in 2012 and 2021, respectively.
Snow water equivalent (SWE) estimates for (a) 2012 and (b) 2021, based upon the mean calculated relative permittivity (
$\overline{\varepsilon_s}$) of 2.48 (2012) and 3.56 (2021) (Landsat 8 Imagery, Band 8).
Comparison of the range of annual accumulation snow water equivalent (SWE, cm w.e.) for 2012 and 2021 based on the maximum, average (±1 SD) and minimum values for calculated relative permittivity (
${\varepsilon _s}$) for the shared extent of radar surveys. 
$\varepsilon_s$ used are provided beneath each SWE range (italicized).
5. Discussion
5.1. Spatial and temporal variability in relative permittivity and LWC
Spatial variability in 
${\varepsilon _s}$ across JIF was lower in 2012 relative to 2021 (Fig. 3). Previous studies have documented spatial variability in relative permittivity, for example in seasonal snowpacks in Colorado and Montana, where measured permittivity values range between ∼1.3 and 2.4 (Webb and others, Reference Webb, Musselman, Ciafone, Hale and Molotch2022) and were calculated to be ∼1.4–1.6 (Bradford and others, Reference Bradford, Nichols, Mikesell and Harper2009) and ∼1.5–2.5 (Bonnell and others, Reference Bonnell2023). The range of 
${\varepsilon _s}$ in 2012 calculated here (2.06–2.62) is slightly greater than but more similar to 
${\varepsilon _s}$ measured in those seasonal snowpacks than the range in 2021 (2.11–5.11). In this glaciological setting, differences in relative permittivity are primarily due to differences in density and LWC. In this way, we attribute the wider range of 
${\varepsilon _s}$ reported here to, in part, the greater wet snow density observed on JIF (0.460–0.600 g cm−3) when compared to these seasonal snow studies (∼0.230–0.500 g cm−3), but more strongly attribute this to differences in LWC.
Observed density patterns across JIF (Section 4.3) are similar to those previously reported at this location (LaChapelle, Reference LaChapelle1954; Pelto and Miller, Reference Pelto and Miller1990). Since densification takes weeks to months, while meltwater production and percolation occur on timescales of hours to days, and our results include permittivity values greater than that of glacier ice (3.15), we attribute the variability of calculated relative permittivity in 2012 and 2021 to LWC within the snowpack. LWC values calculated here are spatially and temporally variable, ranging between 1.7% and 16% (by volume), similar to the range presented by Webb and others (Reference Webb, Jennings, Fend and Molotch2018) for seasonal snowpacks (0–19%). On a smaller spatial scale, Techel and Pielmeier (Reference Techel and Pielmeier2011) carried out snow fork surveys in the Swiss Alps to demonstrate spatial variability in LWC of 0–10% (by volume) in a 5 m wide by 80 cm deep snow pit. LWC results presented here agree with the findings of spatial (Webb and others, Reference Webb, Jennings, Fend and Molotch2018) and temporal variability (Techel and Pielmeier, Reference Techel and Pielmeier2011) in LWC.
In addition to spatial variability in 
${\varepsilon _s}$ and LWC, there is also likely substantial temporal variability of meltwater in this snowpack during the summer months, a time of high ablation and rain-on-snow events. Colbeck (Reference Colbeck1975) modeled how different snowpacks respond to rain-on-snow events and finds water moves fastest though an isothermal snowpack, as compared to fresh snow or refrozen snow. Temporal variability in 
${\varepsilon _s}$ and LWC in snow is found to vary diurnally, between days and between seasons in seasonal snowpacks (Techel and Pielmeier, Reference Techel and Pielmeier2011; Schmid and others, Reference Schmid, Koch, Heilig, Prasch, Eisen, Mauser and Schweiser2015; Bonnell and others, Reference Bonnell, McGrath, Williams, Webb, Fassnacht and Marshall2021), and on a glacier in Western Norway (Hart and others, Reference Hart, Rose and Martinez2011). Bonnell and others (Reference Bonnell, McGrath, Williams, Webb, Fassnacht and Marshall2021) similarly observed spatial and temporal variability in LWC in a seasonal snowpack between April and June in Colorado with LWC values remaining low in the spring and increasing during the summer with some areas reaching >10% LWC. In our study, there was not a consistent trend of an increase or decrease in 
${\varepsilon _s}$ at sites repeated in both years, agreeing with this temporal variability observed elsewhere (Figs. 3, 4).
Values of 
${\varepsilon _s}$ in this study are higher than most other ground-truthed GPR studies on glaciers, which are typically conducted in spring, whereas our data were collected in mid-July during a time when temperatures are consistently above freezing and precipitation may fall as rain. Clayton (Reference Clayton2022) conducted a snowpack meltwater percolation study on JIF, near the Matthes–Llewellyn ice divide (∼1800 m), in 2018 and 2019 and observed the movement of meltwater within this snowpack on sub-daily timescales and trends indicative of diurnal melting and refreezing. In a dense and isothermal snowpack such as on JIF, we expect surface melt or rain to move through the snowpack, but this can be impeded by ice lenses. Evidence of ice lensing on JIF supports the idea that there are meltwater or rain events where water percolates deeper into the snowpack before refreezing into an ice lens. Previous work has demonstrated horizontal hydrological connectivity in snowpacks and the effect of snowpack layering and ice lenses on vertical meltwater movement (Webb and others, Reference Webb, Wigmore, Jennings, Fend and Molotch2019; Webb and others, Reference Webb, Musselman, Ciafone, Hale and Molotch2022). It is likely that the ice lenses within the snow impact this water movement and result in increased snow permittivity. Here, the heterogeneity of ice lensing and water movement through snowpacks may additionally help explain such varied and high values of ԑ s on JIF.
Spatial variability in the LWC varies depending on the time of day of radar transects (Webb and others, Reference Webb, Jennings, Fend and Molotch2018; Bonnell and others, Reference Bonnell, McGrath, Williams, Webb, Fassnacht and Marshall2021) and the prior meteorological conditions (Hart and others, Reference Hart, Rose and Martinez2011). In both 2012 and 2021, the average temperature during the days of data collection was above freezing (8.17 ± 1.43°C and 7.21 ± 0.91°C, respectively) at 1196 m on JIF from lapse-rate derived temperatures (Supplementary S.6., Supplementary Fig. S3). However, there is no strong correlation between relative permittivity and air temperature in either year. Some of the high variability in the 2021 data and low variability in the 2012 data may be attributed to precipitation (or lack thereof), but there are no in situ continuous precipitation data available for this study site, and the meteorological data from the Juneau International Airport are too localized to represent precipitation on JIF. It is possible that differences in precipitation and solar insolation contribute to this difference in water content between years, but this cannot be ascertained with the present data. Therefore, this is an area of study that should be investigated in future work.
5.2. Snow water equivalent
The GPR surveys still show patterns that are a result of climatological and topographical factors such as elevation and proximity to the moisture source. We observe SWE gradients with elevation of 30 cm/100 m (R 2 = 0.46) in 2012 and 21 cm/100 m (R 2 = 0.26) in 2021, computed with 
$\overline {{\varepsilon _s}} $ for each year (Supplementary Fig. S6(a, c)). These values are within the range of accumulation gradients presented in McNeil and others (Reference McNeil, O’Neel, Loso, Pelto, Sass, Baker and Campbell2020) from over 20 years of JIF data, and similar to values found by McGrath and others (Reference McGrath2015) ranging between 11.5 and 40.0 cm/100 m for glaciers near the Gulf of Alaska. The impact of elevation may be amplified on marine-proximal glaciers, where the high-elevation areas are also those closest to the coast and therefore receive more precipitation than interior regions of JIF. We observe an inverse relationship between distance from the coastline and snow accumulation depth and SWE of −7 cm km−1 (R 2 = 0.24) in 2012 and −9 cm km−1 (R 2 = 0.46) in 2021 (Supplementary Fig. S6(b, d)). These trends are most notable for more coastal branches of JIF, in agreement with a comparison between Lemon Creek Glacier and Taku Glacier in McNeil and others (Reference McNeil, O’Neel, Loso, Pelto, Sass, Baker and Campbell2020). These results also support the modeling presented in Roth and others (Reference Roth, Hock, Schuler, Bieniek, Pelto and Aschwanden2018) who find that the area of JIF closest to the coastline has the highest winter precipitation and aligns with the wind direction most dominant between October and March.
While these results align within expected values and trend direction, the magnitude of these relationships are impacted by the applied 
${\varepsilon _s}$, either being amplified or diminished (Supplementary Tables S2 and S3). Both the conservative (8–33%) and extreme (13–45%) differences in mean SWE from GPR due to different 
${\varepsilon _s}$ indicate the potential for over- or underestimation in mass balance and water resources. Due to the lack of spatial trends in 
${\varepsilon _s}$, it is not possible at this time to establish these relationships more accurately. The implications of this may either steepen or lessen mass balance gradients and support the need for resolving 
${\varepsilon _s}$ and LWC variations at higher spatial and temporal resolution.
5.3. Survey implications in a warming climate
The onset of snowmelt across Alaska is occurring earlier in the spring due to increasing atmospheric temperatures (Zheng and others, Reference Zheng, Jia and Xu2022). Ballinger and others (Reference Ballinger2023) showed that between 1957 and 2021, the Central Panhandle of Alaska, which includes JIF, experienced a −3.13% change in snowfall equivalent per decade and increases in 2 m air temperature in every season (Ballinger and others, Reference Ballinger2023). Under modeled future warming scenarios, southeastern Alaska is predicted to experience warmer temperatures and less precipitation (Bigalke and Walsh, Reference Bigalke and Walsh2022). Alaskan glacier systems are amongst those predicted to lose mass and contribute the most to global sea-level rise in the next century (Oerlemans and others, Reference Oerlemans2005; Rounce and others, Reference Rounce2023). With warmer seasons throughout the year, JIF could experience an earlier onset of melt or more predominant spatial and temporal variability in LWC. In this way, perhaps even spring surveys of SWE using GPR should consider increasing the spatial and temporal sampling density of ground-truth snow depths for accurate interpretations.
These trends extend beyond Southeast Alaska, with broader implications for other areas in Alaska and for areas of the Arctic that typically have dry snow. Temperatures are rising faster at the poles than at mid-latitudes, especially in the Northern Hemisphere (Holland and Bitz, Reference Holland and Bitz2003). Rain-on-snow events have recently occurred at Summit Station, Greenland (Xu and others, Reference Xu, Yang, Hu, Liang and Vihma2022) and are predicted to occur under future warming scenarios in Svalbard (Hansen and others, Reference Hansen2014). Such changes may require accounting for LWC variability considerations in areas that historically have dry snowpacks.
The spatial variability in 
${\varepsilon _s}$ and, by implication, LWC presented here supports the need for an increased sampling density of ground-truth points. This study demonstrated a large spatial variability in LWC (2.1–16% in 2021), even within 7–10 days, which negates the assumption that conditions are homogenous across a large (>500 km2) study site. Variations of LWC occur on annual, daily, or even shorter timescales, and the degree of variability between the survey years at the same location demonstrates that we cannot assume the same conditions between years even if surveys occur at the same time as years prior. Understanding more about the spatial and temporal variability of meltwater in snow helps us potentially expand the window for data collection of SWE with GPR if we can properly account for the water.
To remedy this situation, we recommend ground-truth points via snow pits or coring at least at the end-point of GPR surveys and potentially in the middle of survey lines if there are substantial changes in topography along the transect, such as more than 300 m of elevation change, rolling topography, or sudden increases or decreases in topographic wind shading. Examples of factors that may spatially and temporally influence snowpack relative permittivity include (but are not limited to) elevation change, aspect, and solar barriers such as nunataks (which may shield the sun and reduce melt). Newer techniques including multi-offset radar may mitigate some of these limitations where the antenna configuration consisting of multiple transmitters and receivers can continuously derive the layer radar wave velocity across the survey (Bradford and others, Reference Bradford, Nichols, Mikesell and Harper2009; Greissinger and others, Reference Greissinger, Mohr and Jonas2018; Meehan and others, Reference Meehan2021). Future research should consider annual surveys, sub-seasonal surveys, and the use of multi-offset radar to better quantify these trends and variability.
6. Conclusions
In this study, we analyzed 2 years of 150 km GPR surveys with ground-truth points from snow pits and snow cores in 2012 and 2021 to assess the variability of summer snowpack 
${\varepsilon _s}$ on JIF. We calculated LWC at each ground-truth location using an empirical relation to infer snowpack LWC ranging from 1.7% to 16%. Finally, we provided bounds on the SWE across JIF in 2012 and 2021 using density data from 31 summer snow pits and the calculated 
${\varepsilon _{s,max}}$, 
${\varepsilon _{s,avg}}$, and 
$\varepsilon_{s,min}$ from 13 ground-truth points across JIF. We find that depth-density trends across JIF within and between years are relatively consistent and that, across the range of calculated 
${\varepsilon _s}$, differences in interpreted depths due to 
${\varepsilon _s}$ introduce substantial uncertainty to SWE calculations (8–33%, with extreme cases of 13–45%) by GPR.
GPR is an effective, nondestructive tool to rapidly extend point-based measurements of SWE across large spatial extents (10s to 100s of km). However, in wet snow environments, GPR assumptions about relative permittivity must be made with caution. Our results show that SWE estimates from GPR in such conditions require extensive ground-truth campaigns due to the spatial and temporal variability in snowpack LWC. As airborne and ground-based radar surveys become more common in high-latitude glacier and ice sheet studies, the expected increase in LWC due to rising temperatures may complicate GPR interpretation, especially in Alaska, western Canada, and polar areas. To ensure accurate interpretations of depth and SWE, ground-truth campaigns should be conducted as close in time as possible to the time of surveying and across the broadest spatial scale feasible. Future work should aim to resolve this LWC variability at higher spatial and temporal resolution to improve SWE estimates.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/aog.2025.10014.
Data availability statement
GPR and Summer Snow Pit Data are available through the Arctic Data Center (Campbell, Reference Campbell2024). Spring Snow Density Data are available through USGS data release (McNeil and others, Reference McNeil, Sass, Baker and O’Neel2017). NWS NOAA NOWData Meteorological Data are publicly available at https://www.weather.gov/wrh/Climate?wfo=ajk.
Acknowledgements
We are grateful for the input from our reviewers, Daniel McGrath, Jack Kohler, and Andrea Creighton, whose feedback greatly improved the manuscript, as well as Peter Koons for his advice on data analysis and Jonathan Maurer for preliminary data processing. We thank the Foundation for Glacier and Environmental Research and the many students, staff, and faculty of the JIRP for contributions to this research in data collection, logistical support, and preliminary data analysis. Funding from the National Science Foundation under award no. 2239668 also supported data analysis and interpretation in this research. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
Author contributions
M.M. conducted data processing, data analysis, and initial writing and revision of the manuscript. S.B. conducted both years of data collection, preliminary data processing, and revisions to the manuscript. S.C. provided funding acquisition, supervision, idea conceptualization, both years of data collection, software resources, and extensive revisions. E.S. conducted the second year of data collection and provided feedback on revisions. K.S. contributed guidance on data analysis and revisions. C.M. provided advice and interpretations of snow pit ground-truth data.
Competing Interests
The authors declare no competing interests.
Open access
