1. Introduction
Peripheral glaciers and ice caps are often particularly vulnerable to climate warming and exhibit more rapid responses to climatic shifts, making them sensitive indicators of climate change (Bjørk and others, Reference Bjørk2018). Peripheral glaciers and ice caps in Greenland play a critical role in environmental and climatic processes. Although they account for only 4% of Greenland’s glaciated area, they contribute around 11% to the overall ice loss from Greenland (Khan and others, Reference Khan2022). Moreover, because many of these glaciers and ice caps occupy ecologically diverse catchments, variations in their mass balance and freshwater runoff can significantly affect downstream ecosystems (Hopwood and others, Reference Hopwood2020).
Snow accumulation is one of the most challenging mass-balance components to capture in surface mass-balance models. While surface melt is largely governed by atmospheric factors that scale with elevation, snow accumulation remains more difficult to model, especially in rugged terrain, because the spatial distribution of snow accumulation depends on intricate interactions among topography, wind transport and redistribution processes (Hodgkins and others, Reference Hodgkins, Cooper, Wadham and Tranter2005; Machguth and others, Reference Machguth, Eisen, Paul and Hoelzle2006; Dadic and others, Reference Dadic, Mott, Lehning and Burlando2010; Messerli and others, Reference Messerli, Arthur, Langley, How and Abermann2022; Voordendag and others, Reference Voordendag, Goger, Prinz, Sauter, Mölg, Saigger and Kaser2024). Due to the difficulty and large computational requirements of modelling snow drift by wind, most mass balance models use a conceptual approach where accumulation is a function of elevation and topographical parameters such as slope, aspect and curvature (Sold and others, Reference Sold2016). However, assessing the heterogeneity in snow distribution is critical because the amount of winter accumulation strongly influences glacier mass balance as it affects the timing of glacier ice exposure and consequently the albedo feedback mechanisms during the melt season (Réveillet and others, Reference Réveillet, Vincent, Six and Rabatel2017). Sold and others (Reference Sold2016) demonstrated that mass balance estimates of an Alpine glacier drastically improved when spatially extensive snow depth measurements were assimilated into the model.
Despite the importance of winter snow accumulation in determining annual glacier mass balance, spatially extensive measurements are scarce for peripheral glaciers and ice caps in Greenland (Mernild and others, Reference Mernild, Kane, Hansen, Jakobsen, Hasholt and Knudsen2008; Machguth and others, Reference Machguth2016; Tsutaki and others, Reference Tsutaki, Sugiyama, Sakakibara, Aoki and Niwano2017; Hynek and others, Reference Hynek2024). Existing in situ observations—largely from automatic weather stations (AWS), snow pits or manual snow probing—are by nature point-based and may not adequately represent glacier-wide accumulation patterns (Sold and others, Reference Sold2013; Pulwicki and others, Reference Pulwicki, Flowers, Radić and Bingham2018). Ground-penetrating radar (GPR) has become a crucial tool for mapping snow accumulation on glaciers and ice sheets as it can efficiently cover large areas (Kohler and others, Reference Kohler, Moore, Kennett, Engeset and Elvehøy1997; Machguth and others, Reference Machguth, Eisen, Paul and Hoelzle2006; Miège and others, Reference Miège, Forster, Box, Burgess, McConnell, Pasteris and Spikes2013; Hawley and others, Reference Hawley, Courville, Kehrl, Lutz, Osterberg, Overly and Wong2014; McGrath and others, Reference McGrath2015; Lewis and others, Reference Lewis2019; Karlsson and others, Reference Karlsson, Razik, Hörhold, Winter, Steinhage, Binder and Eisen2020). For instance, towed GPR surveys on Alaskan glaciers have revealed strong snow accumulation patterns reflecting climate contrasts, wind-redistribution and microclimate effects (McGrath and others, Reference McGrath2015), but also interannual persistence of snow distribution patterns (McGrath and others, Reference McGrath, Sass, O’Neel, McNeil, Candela, Baker and Marshall2018). Helicopter-borne systems, for example applied on Alpine glaciers, have resolved fine-scale heterogeneities caused by wind redistribution, highlighting accumulation patterns that conventional point measurements often miss (Machguth and others, Reference Machguth, Eisen, Paul and Hoelzle2006; Sold and others, Reference Sold2013, Reference Sold2016). While many GPR-derived accumulation studies exist over the Greenland Ice Sheet (e.g. Miège and others, Reference Miège, Forster, Box, Burgess, McConnell, Pasteris and Spikes2013; Hawley and others, Reference Hawley, Courville, Kehrl, Lutz, Osterberg, Overly and Wong2014; Lewis and others, Reference Lewis2019; Karlsson and others, Reference Karlsson, Razik, Hörhold, Winter, Steinhage, Binder and Eisen2020), the only known GPR application on a smaller ice caps and glaciers in Greenland is on Freya Glacier (a few kilometres from A.P. Olsen Ice Cap). Here, extreme snowfall events and subsequent avalanches can locally enhance snow accumulation by two- to three-fold (Hynek and others, Reference Hynek2024).
Reanalysis data provide an alternative for estimating snow accumulation where direct observations are limited. The Copernicus Arctic Regional Reanalysis (CARRA) dataset offers atmospheric variables—including precipitation and snow water equivalent (SWE, which is the depth of water that would result if the mass of snow melted completely)—at a 2.5 km grid resolution (Schyberg and others, Reference Schyberg2020). Recent studies have started evaluating the performance of CARRA precipitation and SWE outputs over the Arctic, coastal Greenland (van der Schot and others, Reference van der Schot, Abermann, Silva, Rasmussen, Winkler, Langley and Schöner2024) and non-glaciated regions in Svalbard (Maniktala, Reference Maniktala2022). Additionally, some studies have integrated CARRA variables into glacier mass balance models for Svalbard (Schmidt and others, Reference Schmidt, Schuler, Thomas and Westermann2023) and used CARRA precipitation data to examine rain-on-snow events on the Greenland Ice Sheet (Box and others, Reference Box2023). However, these have primarily relied on sparse point measurements, and no systematic, glacier-wide assessment has yet been conducted to evaluate how well CARRA represents the spatial distribution of winter snow accumulation.
In this study we present a multi-year dataset of winter snow accumulation covering mainly the eastern catchment of A.P. Olsen Ice Cap (APO). These data were collected by extensive GPR surveys conducted in 12 years between 2008 and 2024. Building on this dataset, our objectives are to (i) analyse the spatial and temporal variability of snow accumulation over the 16-year period, (ii) evaluate how well existing AWS capture glacier-wide snow accumulation and (iii) evaluate CARRA’s ability to reproduce observed snow accumulation over APO. To this end, we compare the GPR-derived observations to AWS measurements of snow depth and CARRA precipitation. This work offers new insights into the spatial variability of winter snow accumulation on a peripheral ice cap in Northeast Greenland and the suitability of using CARRA precipitation for mass balance modelling in data-sparse regions.
2. A.P. Olsen Ice Cap study site
A.P. Olsen Ice Cap is located in northeast Greenland (Fig. 1) and covers approximately 300 km2, with elevations spanning from ∼200 to 1500 m a.s.l. (above sea level). The eastern catchment drains into the Zackenberg River and has been monitored as part of the Greenland Ecosystem Monitoring (GEM) programme (www.g-e-m.dk) since 2008 (Citterio and others, Reference Citterio2017; Larsen and others, Reference Larsen, Binder, Rutishauser, Hynek, Fausto and Citterio2024). As part of these monitoring efforts, annual GPR surveys and snow density measurements have been conducted regularly from 2008 to 2024, providing a spatially extensive record of winter SWE (Fig. 1). The GPR surveys support three AWS that are installed along an elevation transect on the eastern catchment outlet glacier, continuously recording meteorological variables, snow height, ice ablation and englacial temperatures (Larsen and others, Reference Larsen, Binder, Rutishauser, Hynek, Fausto and Citterio2024).
Figure 1. Overview map of A.P. Olsen Ice Cap (location marked with red dot on the Greenland map) including the GPR survey profiles collected between 2008 and 2024 and the location of the three AWS. The red dotted line represents the average end-of-summer snow line altitude (SLA, 1150 m a.s.l.) contour (Larsen, Reference Larsen2025). Background map is a Sentinel-2 satellite image from 2020 compiled by the Danish Agency for Climate Data, made available at www.dataforsyningen.dk.
Previous investigations found considerable spatial variability in the winter snow accumulation over the Zackenberg River drainage basin (including both glaciated and non-glaciated areas), largely driven by wind redistribution and topographical influences (Mernild and others, Reference Mernild, Liston and Hasholt2007). Wind observations at the AWS in the accumulation zone (ZAC_A) show that winter winds are primary from the northwest, with a secondary mode from southwest (Fig. S1). As in any river catchment with seasonal snow, the winter SWE over APO and the surrounding terrain has a large impact on runoff in the Zackenberg River where discharge usually peaks in early summer due to snowmelt. Discharge from the Zackenberg River impacts the marine environment by transporting sediment into the fjord which significantly influences light availability and thus the marine ecosystem (Sejr and others, Reference Sejr2022).
3. Data and methods
3.1. Radar data collection and processing
GPR surveys were performed in 12 years between 2008 and 2024, before the onset of melt to assess winter snow-depth distribution over APO (Table 1). Surveys could not be carried out in 2009, 2013 and 2019 due to insufficient snow cover preventing snowmobile access to the glacier, and in 2020 and 2021 due to COVID-19 travel restrictions. A MALÅ Geoscience ProEx unit with 500 MHz (survey years 2008–2011) or 800 MHz (2012–2024) shielded antennas was towed on a plastic sled behind a snowmobile at ∼10 km/h. This setup yielded a mean trace spacing of 0.4 m (ranging from 0.1 to 1.3 m for different survey years). Survey lines were slightly different from year to year, but followed the centre line of the APO eastern catchment outlet glacier and, when possible, zig-zag routes across the glacier were conducted to capture cross-glacier variability (Fig. 1).
Table 1. Overview of the GPR survey parameters along with the mean snow pit density and the reconstructed SWE from radar profiles (SWEGPR) and CARRA total precipitation (SWECARRA_TP). Uncertainty estimates in SWEGPR were derived by propagating the density uncertainty estimates. SWECARRA_TP values represent the mean of all grid cells containing at least 200 m of radar survey lines, and SWEGPR represents the mean of all corresponding radar data points within these CARRA grid cells
(x) marks years with manual snow probings.
We used a handheld (single frequency) Global Navigation Satellite System (GNSS) receiver to record the position with an estimated horizontal uncertainty of ∼5 m. GNSS data were filtered to remove outliers, and missing coordinates due to lower sampling rate of the GNSS relative to the radar were interpolated linearly. Because GNSS elevation accuracy is limited, we used ArcticDEM at 32 m resolution (Porter and others, Reference Porter2023) to obtain surface elevations along the GPR profiles.
Radar data were processed in the ReflexW-2D software package (Sandmeier Scientific Software) using the following steps: application of a dewow filter to remove low-frequency noise, a time-zero correction to align traces to the first significant peak (one value was used per survey year), band-pass filtering and background removal to remove noise and enhance signal clarity, and a linear gain function to amplify deeper reflections (Cassidy, Reference Cassidy and Jol2009).
The last summer surface (LSS) was identified using a combination of semi-automated phase-follower and manual picking, depending on the continuity of the reflectors. In regions where no clear reflector was visible, no layer was traced. To maintain consistency, LSS identification was guided by snow pit stratigraphy, manual snow probing and crossover radar profiles. The transition from winter snow to glacier ice in the ablation zone produces a pronounced reflection (yellow lines in Fig. 2a-b), while multiple internal melt or rain layers and older summer surfaces in the accumulation zone can complicate LSS identification (Fig. 2c-d).
Figure 2. Radargrams showing a typical GPR profile collected in 2022 in (a) the ablation zone of APO, with picked LSS marked in (b), and (c) collected in the accumulation zone with picked LSS marked in (d).
Snow depth was calculated by converting two-way travel time using column-averaged radar velocity derived from the empirical relationship by Kovacs and others (Reference Kovacs, Gow and Morey1995), based on measured snow densities (
${\rho _{snow}}$, Table 1), resulting in radar velocities ranging from 0.229 to 0.236 m/ns. Velocity uncertainties, calculated by propagating snow density error estimates, range from 0.005 to 0.015 m/ns, corresponding to snow depth uncertainties of 3–6 cm. Snow water equivalent (
$SWE = {h_{snow}}\ast{\rho _{snow}}/{\rho _{water}}$) was obtained from snow depth (
${h_{snow}}$) also using the measured snow densities. In years without density measurements, we use the mean column-averaged density (344.4 ± 31 kg/m3). Finally, a 2 m moving average filter is applied to the GPR snow depths, followed by down sampling to 2 m equidistant intervals along the profiles. This removes SWE data points where the radar was standing still (typically at the start and end of the profile) and simplifies visualization.
SWE uncertainties were estimated by combining uncertainties from the travel time to snow depth conversion and propagating the snow density error estimates in the SWE calculation. Additional uncertainty may arise from manually picking the LSS, since the exact reflector can be ambiguous or multiple reflectors may be present. Incorrect picking of our dataset by one pixel would translate to an uncertainty of 0.9–2.8 cm. However, the picking was performed by a single interpreter across all survey years, ensuring internal consistency across the dataset.
3.2. Snow probes and snow pits
Manual probing was conducted in five of the survey years to measure the depth from the snow surface to the ice surface or LSS along the GPR profiles. Snow pits were excavated nearly every year. Each pit extended to the LSS, and the stratigraphy (e.g. presence of ice layers) and snow densities were recorded at 10–20 cm intervals. Column-averaged densities from individual snow pits ranged from 261 kg/m3 to 428 kg/m3 (Fig. 3) and showed no clear relationship with elevation, slope, aspect or snow pit depth. Thus, we use the mean of all column-averaged snow pit densities for each year to convert snow depths into SWE (Table 1). Snow density uncertainty is estimated from the standard deviation of all column-averaged densities for each year, with an overall mean uncertainty of 31 kg/m3.
Figure 3. Snow pit depth-density profiles over A.P. Olsen Ice Cap from 2008 to 2023. Black dashed vertical line is the mean density over all profiles for a given year, with the grey dashed lines indicating one standard deviation from the mean, used as uncertainty estimate. Snow pit locations are shown in Fig. S3.
3.3. AWS snow height
The three AWS on APO (ZAC_L, ZAC_U, ZAC_A) measure snow height using sonic ranger sensors (Campbell Scientific SR50). The data used here has been processed following the approach described in Larsen and others (Reference Larsen, Binder, Rutishauser, Hynek, Fausto and Citterio2024). Additionally, a one-day running mean filter was applied to the hourly data in order to remove sensor noise and outliers. Snow depth at each AWS was derived by subtracting the measured distance to the surface from the maximum sensor height observed in July or August each year. To estimate SWE from these point measurements, we used the same mean annual snow density as applied to the GPR data.
3.4. CARRA dataset and evaluation
The Copernicus Arctic Regional Reanalysis (CARRA) dataset provides atmospheric and surface variables over the Arctic at a 2.5 km horizontal grid resolution, with three-hourly outputs extending from 1991 to the present (Schyberg and others, Reference Schyberg2020). CARRA is produced via the HARMONIE-AROME numerical weather prediction system (Bengtsson and others, Reference Bengtsson2017), forced by ERA5 reanalysis (Hersbach and others, Reference Hersbach2020) and a lot of additional surface observations from various programmes such as ASIAQ (Greenland Survey), PROMICE (Programme for Monitoring of the Greenland Ice Sheet [Fausto and others, Reference Fausto2021]) and GC-Net (Greenland Climate Network [Vandecrux and others, Reference Vandecrux2023]). Here, we used CARRA total precipitation to derive winter snow accumulation (hereafter referred to as SWECARRA_TP) and compare it to the GPR observations. SWECARRA_TP is calculated from daily precipitation forecasts (+6 h and + 30 h lead times initiated at 00:00 UTC), summed from the onset of winter to the GPR survey date. We defined the onset of winter as the first day of 5 consecutive days with daily maximum temperature below 0°C, based on CARRA 2 m air temperatures at the ZAC_U AWS location (Table 1). While atmospheric variables such as temperature and relative humidity from weather stations (including the ASIAQ Zackenberg climate station installed on land ∼30 km southeast of APO) are assimilated, precipitation is purely model-derived (Yang and other, Reference Yang2020). We note that CARRA also provides a SWE variable; however, a preliminary evaluation revealed a systematic underestimation of winter SWE over APO (Fig. S9), and we therefore only use the total precipitation SWECARRA_TP.
To compare CARRA outputs with the GPR-derived SWE (SWEGPR), we averaged SWEGPR values within each CARRA grid cell that contained at least 200 m of GPR profiles. We quantified agreement between the two datasets using the Pearson correlation coefficient (r), the root mean square error (RMSE, for absolute differences) and the percent difference (to express relative differences):
$${\text{RMSE }}\left[ {{\text{m w}}.{\text{e}}.} \right]\, = \sqrt {\frac{1}{n}\mathop \sum \nolimits^ {{\left( {SW{E_{CARRA}} \,-\, SW{E_{GPR}}} \right)}^2}} , \kern56.5pt$$
$${\text{Percent difference }}\left[ \% \right]\, = \frac{{\left( {SW{E_{CARRA}} \,-\, SW{E_{GPR}}} \right)}}{{SW{E_{GPR}}}}100\% , \kern36.3pt(2)$$4. Results
4.1. Validation of GPR-derived snow depths with manual snow probes
From 2008 to 2024, 561 km of radar data were collected (on average 47 km per survey year), and snow depths were successfully derived along 92% of this total profile length. Fig. 4 compares GPR-derived SWE (SWEGPR, averaged over the 10 nearest traces) with manual snow probe measurements (SWEProbe).
Figure 4. Comparison of GPR-derived snow water equivalent (SWEGPR) with manual snow probings (SWEProbe) in the ablation (black) and accumulation zones (blue). SWEGPR represents the mean of the 10 closest radar traces to the snow probing location. The black- and blue lines indicate the linear regression fits for the ablation- and accumulation zones, respectively.
In the ablation zone, SWEGPR strongly correlates with SWEProbe (R2 = 0.81, RMSE = 0.05 m w.e.). A few outliers where SWEGPR exceeds SWEProbe likely reflect probes stopping at internal ice layers rather than the true glacier surface (e.g. Fig. 2b, probe at end of the profile).
In the accumulation zone, the correlation is weaker (R2 = 0.4, RMSE = 0.12 m w.e.), and discrepancies show both positive and negative offsets. Multiple internal layers can cause the probes to measure above or below the LSS. In contrast, the LSS picked from the radar layers was traced consistently across the dataset and carefully quality controlled for any misalignments or human errors (e.g. Fig. S2). We therefore ascribe the discrepancies to the inherent uncertainties of probe observations. Potential error sources related to the radar layer picking include uncertainties in the radar velocities, the GPR’s vertical range resolution (∼7–12 cm) and Fresnel zone effects (i.e. GPR snow depths are integrated over a broader footprint compared to the manual snow probing measurement). Nevertheless, the strong correlation in the ablation zone supports the reliability of GPR-derived snow depths, with RMSE values of 0.05 m w.e. in the ablation zone, and 0.12 m w.e. in the accumulation zone serving as additional uncertainty estimates.
4.2. Spatiotemporal distribution of winter snow accumulation
The spatial distribution of GPR-derived winter SWE (SWEGPR) exhibits considerable interannual variability across the APO eastern catchment area and outlet glacier, ranging from 0.3 m w.e. in 2011 to 0.93 m w.e. in 2015 (Fig. 5, Table 1). Although the accumulation zone typically has a thicker snowpack compared to the ablation zone, no clear linear relationship between SWE and elevation is observed (Fig. 6a). Similarly, no significant glacier-wide relationships with aspect or slope (using the ArcticDEM at 32 m resolution [Porter and others, Reference Porter2023]) were detected (Table S1). We note that although the total survey length varied between years, the profiles consistently covered similar regions, and we therefore consider the annual mean SWE values to be representative.
Figure 5. GPR-derived winter SWE (SWEGPR, line-data with white background) on top of CARRA SWE calculated from accumulated total precipitation (SWECARRA_TP) over APO from 2008 to 2024. CARRA grid cells are limited to cells containing at least 200 m of radar survey lines. The black dotted line represents the average end of summer snow line altitude (SLA, 1150 m a.s.l.) contour, and thin black lines indicate elevation contours at 100 m intervals. The white triangles show the AWS locations.
Figure 6. (a) Winter SWE versus elevation. Grey lines represent individual SWEGPR survey years, and the black line indicates the mean SWEGPR across all surveys. Data were binned into 50 m elevation bands repeated every 25 m. For comparison, SWECARRA_TP is averaged over all survey years, and only grid cells with GPR data are considered. Blue dots represent the mean SWECARRA_TP across all years for each grid cell, and the dashed lines indicate linear interpolation between the grid cells. Black triangles are the mean SWE derived from the sonic ranger snow height measurements at the AWS. (b) Along-profile standard deviation of SWEGPR (calculated over 50 m segments) versus elevation.
Despite the low correlation between SWE and surface elevation, the data reveal distinct zones of accumulation consistent across most survey years (Fig. 6a). For instance, the glacier terminus (<675 m a.s.l.) consistently has a localized thicker snowpack—on average 17% greater than the rest of the ablation zone. Moving upslope, snow depths are smaller and relatively uniform between 700 and 900 m a.s.l. and then increase again towards the equilibrium line where the narrow outlet glacier transitions onto the wider firn plateau. In the accumulation zone, snow depths are relatively uniform but decrease slightly above 1400 m a.s.l. Small-scale spatial variability, expressed as the standard deviation of SWE over 50 m intervals along the radar profiles, is higher in the ablation zone than in the accumulation zone (Fig. 6b, Fig. S5). Unbinned SWEGPR versus elevation for each year are shown in Fig. S4, and maps of along-track variability for each survey year are presented in Fig. S5.
4.3. Distribution of winter snow accumulation in the vicinity of the AWS stations
AWS stations ZAC_L and ZAC_U are located where the ablation zone shows high small-scale (quantified over 50 m profile segments) variability in SWE (Fig. 6b). Comparing sonic ranger-derived SWE (SWESR50) with GPR-derived SWE (SWEGPR) from within a 25 m radius shows that SWESR50 underestimate surrounding SWE by 25% (ZAC_L) and 32% (ZAC_U) (Fig. 7a, b). Expanding the radius to 1 km increases the underestimation to 45% at ZAC_L and 40% at ZAC_U. This comparison indicates that the ablation zone AWS, especially ZAC_L, are situated in locally low-accumulation areas, which can lead to significant underestimation of winter SWE when extrapolating this dataset to the broader glacier area. Considering the snow pit measurements instead of the SR50 measurements, a similar pattern of general underestimation with respect to the surrounding 1 km occurs at ZAC_L, whereas at ZAC_U snow pit measurements are in closer agreement with the surrounding SWEGPR.
Figure 7. (a–c) Comparison of the sonic-ranger (black) and snow pit (orange) derived SWE at the three AWSs (ZAC_L, ZAC_U, ZAC_A) with GPR-derived SWE (SWEGPR) within 25 m and 1 km radii from the AWS, (d–f) Mean SWEGPR within 25m versus 1 km of the AWS locations.
At ZAC_A in the accumulation zone, SWESR50 underestimates local SWE by 60% at both 25 m and 1 km radii (Fig. 7c). However, when comparing SWEGPR with snow pit-derived SWE at ZAC_A, we find a good agreement, suggesting that the apparent underestimation in SWESR50 is likely due to measurement limitations or issues when converting sonic ranger heights to snow depths. In particular, there are significant gaps during the winter months for several years, and in the early part of the record, the AWS may have sunk into the firn due to densification, affecting snow height measurements. In contrast to the AWS in the ablation zone, comparing SWEGPR from within a 25 m radius to a 1 km radius around ZAC_A shows negligible difference (0.9%), implying that ZAC_A’s location is more representative of surrounding conditions (Fig. 7d–f).
4.4. Evaluation of CARRA winter SWE
At the glacier-wide scale, CARRA-derived SWE from accumulated total precipitation (SWECARRA_TP) captures the overall pattern observed from SWEGPR well, with generally higher SWE in the accumulation zone than the ablation zone (Fig. 5). When averaged over all grid cells containing GPR data, SWECARRA_TP correlates well with SWEGPR (r = 0.93, RMSE = 0.07 m w.e., Fig. 8). Notably, there is no clear relationship between the extent of over- or underestimation in SWECARRA_TP and the total amount of winter snow accumulation (Fig. 9a).
Figure 8. Correlation between mean SWE from the radar measurements (SWEGPR) and CARRA accumulated total precipitation (SWECARRA_TP). CARRA values represent the mean of all grid cells containing at least 200 m of radar survey lines, and SWEGPR represents the mean of all corresponding radar data points within these CARRA grid cells. The solid black line shows the linear regression lines, and the root mean square error (RMSE) is the error between the CARRA SWE and SWEGPR values.
Figure 9. (a) Comparison of the mean SWE from the GPR measurements and CARRA accumulated total precipitation for all radar survey years. Only CARRA grid cells and GPR points were used where the grid cells contain at least 200 m of radar survey lines. (b) Percentage difference between SWEGPR and SWECARRA_TP against elevation for each year. Each data point represents a single CARRA grid cell with corresponding radar data, with the elevation taken as the mean of the GPR data points within that grid cell. Black dots indicate CARRA grid cells only surveyed in 2018. Black line is the mean over all years (100 m elevation bins repeated every 100 m). A map showing the percentage difference for each grid cell is shown in (Figure S6).
Although SWECARRA_TP provides a reasonable estimate of the total mean winter accumulation over APO, it fails to reproduce the complex, non-linear relationship between snow depth and elevation observed in the GPR data (Fig. 9b). In the ablation zone, CARRA omits much of the thick snowpack near the outlet glacier terminus, leading to an average underestimation of SWE by 16%. In contrast, CARRA tends to overestimate SWE by 17% in the area between the ZAC_L and ZAC_U stations. In the accumulation zone, CARRA underestimates SWE in most years, but significantly overestimates in 2018, skewing the overall mean across all survey years to a 6% overestimation. The time-series of CARRA accumulated precipitation compared to the SR50 measurements suggests that the CARRA SWE overestimation in 2018 stems from an exceptional snowfall event in the CARRA model in February 2018 (Fig. S8).
5. Discussion
5.1. Spatiotemporal variability of winter snow accumulation
Our multi-year GPR dataset highlights both significant interannual variability in winter SWE and relatively consistent large-scale altitudinal patterns. Interannual repeatability of snow distribution patterns have previously been observed in Alaska and the Alps, largely governed by stable topographic features (e.g. ridges) and prevailing wind patterns (Sold and others, Reference Sold2016; McGrath and others, Reference McGrath, Sass, O’Neel, McNeil, Candela, Baker and Marshall2018). Although the accumulation zone of APO typically exhibits higher winter snow accumulation than the ablation zone, there is only a weak correlation with elevation, implying that wind-driven deposition and redistribution likely play a dominant role. This aligns with previous studies suggesting that wind is the dominant factor in snow accumulation patterns on Alpine glaciers (Machguth and others, Reference Machguth, Eisen, Paul and Hoelzle2006; Dadic and others, Reference Dadic, Mott, Lehning and Burlando2010; Godio and Rege, Reference Godio and Rege2016), which has also been observed over the greater Zackenberg region (Mernild and others, Reference Mernild, Liston and Hasholt2007). Wind speed data from the three AWS stations (Larsen and others, Reference Larsen, Binder, Rutishauser, Hynek, Fausto and Citterio2024) suggest that conditions for wind-driven snow redistribution may be common at APO during winter, with wind speeds above 7.5 m/s (Li and Pomeroy, Reference Li and Pomeroy1997) on 32% of winter days at ZAC_A, 21% at ZAC_U, and 18% at ZAC_L. We note, however, that in other glacierized regions with larger elevation gradients, such as the Andes, Himalayas, or Alaska, orographic precipitation often plays a more dominant role in determining accumulation patterns (Houze, Reference Houze2012).
Additionally, we note that there is a lack of correlation between snow depth and aspect or slope in our dataset, which is likely due to a limited range of these terrain properties sampled in measurements (Fig. S7, Table S1), making it difficult to fully evaluate correlations with these terrain characteristics.
We consistently observe a thick winter snowpack at the terminus of the APO eastern catchment outlet glacier, presumably caused by the narrow gully and surrounding valley walls favouring wind drifting and trapping snow in that zone. This thick snow cover at the terminus likely plays an important role in delaying the exposure of glacier ice during the early melt season, thereby influencing local melt rates and glacier runoff. While the inability to capture such localized zones of thicker snow in modelling efforts can lead to an overestimation of melt and underestimation of mass balance, underestimation over this relatively small area at APO is expected to have a minor effect on the glacier-wide mass balance.
Overall, we observe a high spatial variability of winter SWE in the ablation zone. Such pronounced heterogeneity has also been observed on Alaskan glaciers (McGrath and others, Reference McGrath2015), and may result from an underlying rough glacier surface including surface depressions from crevasses, supraglacial streams and moulins that can trap snow locally. Additionally, there could be strong turbulent winds altered by the surrounding valley walls rising 300–350 m above the glacier surface that create uneven snow deposition, wind-scoured areas and snow drifts. In contrast, the accumulation zone appears more uniform in terms of snow cover, likely because it comprises a wider firn plateau with smoother underlying LSS and fewer topographic barriers.
5.2. AWS representation of regional snow cover
The three automatic weather stations (AWSs) at APO offer valuable continuous records of snow height, yet our data show that ZAC_L and ZAC_U systematically underestimate snow accumulation in the surrounding area by 40–45% when comparing their immediate vicinity (25 m) to a broader area (1 km). These stations are located near the glacier centreline where the convex ice surface places them on local topographic highs, while the margins form shallow depressions. Winds are typically stronger at the centreline, exposing the AWS locations to more wind erosion of snow, while the margins are more sheltered and prone to snow deposition. The combined effect of local topography and wind redistribution around these stations likely results in lower-than-average snow depths, introducing errors when extrapolating AWS-based measurements over larger areas. In the accumulation zone, ZAC_A seems more representative of its surroundings (the difference between 25 m and 1 km radii is ∼0.9%), but sonic ranger-derived SWE show large deviations from GPR-based measurements in all years. These deviations likely stem from errors in converting sonic ranger heights to snow depths, for example, due to challenges in pinpointing the last summer surface (LSS) in firn areas.
Overall, this comparison underscores the importance of spatially comprehensive measurements for improving our understanding of snow accumulation processes and subsequently model winter mass balance, as point measurements alone—such as those from weather stations—must be interpreted with caution and do not necessarily represent glacier-wide accumulation.
5.3. CARRA performance in winter SWE
Our results show that accumulated total precipitation from CARRA (SWECARRA_TP) correlates reasonably well (r = 0.93) with the mean GPR-derived winter SWE. In most years, SWECARRA_TP therefore provides a robust first-order estimate of the mean total snow accumulation, suggesting that CARRA can be a useful tool for filling gaps in observational data. However, expected given the relatively coarse resolution of the CARRA grid cells compared to the glacier’s scale, SWECARRA_TP does not fully reproduce the complex finer spatial patterns observed from the GPR data. For example, the consistent zone of thick winter snow accumulation at the glacier terminus is notably absent in the CARRA estimate, leading to significant underestimation at lower elevations. Conversely, SWECARRA_TP tends to overestimate SWE compared to the GPR observations at mid-elevations in the ablation zone. This discrepancy may arise because the 2.5 km CARRA grid cells include higher surrounding terrain (on average 120 m higher than the GPR-sampled locations), artificially increasing precipitation estimates, and because CARRA does not account for wind-driven snow redistribution. Furthermore, due to logistical constraints, the majority of the GPR survey lines cover the glacier’s centre, where snowpacks may be generally thinner and more prone to wind erosion. As a result, the higher snow accumulation at the glacier margins are underrepresented in the GPR data but included the CARRA data estimates.
In the accumulation zone, SWECARRA_TP are generally closer to the GPR measurements, suggesting that smoother topography and less wind redistribution at higher elevations may facilitate a more uniform distribution of winter snow accumulation. Still, notable under- or overestimation in some years in the accumulation zone could result from short-term factors such as synoptic weather patterns or anomalous storms affecting the CARRA precipitation (e.g. in 2018), but this requires future investigation and is beyond the scope of this study.
Some of the observed discrepancies may also be caused by errors in GPR-derived SWE estimates due to misidentification of the LSS, in particular in the accumulation zone. Additionally, missing or sparse snow density observations likely introduce some uncertainties in the GPR-derived SWE. Lastly, not accounting for sublimation in SWECARRA_TP could lead to overestimation of SWE; however, Mernild and others (Reference Mernild, Liston and Hasholt2007) note that sublimation in the Zackenberg region is expected to play a minor role.
Overall, our findings align with previous assessments of CARRA outputs in other regions. Schmidt and others (Reference Schmidt, Schuler, Thomas and Westermann2023) compared CARRA-forced mass balance model with in-situ mass balance stake observations in Svalbard and identified a slight underestimation of precipitation over most glaciers. Similar to our conclusions, they attributed these discrepancies to missing topographic and wind-redistribution processes at the model’s grid scale.
Studies using the CARRA snow water equivalent variable, which is directly available from the reanalysis data, found underestimation of snow depths over non-glaciated low-precipitation sites in Svalbard (Maniktala, Reference Maniktala2022). In contrast, van der Schot and others (Reference van der Schot, Abermann, Silva, Rasmussen, Winkler, Langley and Schöner2024) found that this CARRA SWE variable generally performed well in estimating maximum winter SWE at weather stations in coastal Greenland, including the Zackenberg region. It is important to note that those comparisons were conducted on non-glaciated terrain, whereas our study examines glaciated surfaces, where distinct processes influence the energy balance. We briefly compared the CARRA SWE variable against our GPR dataset over APO; however, we found that it systematically underestimates winter SWE (Fig. S9). This is likely because the underlying snow model schemes lack a glacier ice component (Schyberg and others, Reference Schyberg2020), and therefore do not fully capture the energy balance processes on glacier surfaces.
In contrast to previous studies, we use spatially extensive GPR data rather than point-based observations alone, giving deeper insight into sub-grid scale processes. The absence of explicit wind-redistribution schemes in CARRA is evident in our results at APO, particularly in the ablation zone where spatial variability is high. Overlooking these small-scale processes can lead to significant under-or overestimates of snow depth that propagate into mass balance calculations.
6. Conclusions
This study presents a multi-year dataset of winter snow accumulation derived from GPR surveys at A.P. Olsen Ice Cap (APO) in Northeast Greenland, offering new insights into the spatial and temporal variability of snow accumulation on a peripheral glacier. Our results show that wind-driven deposition and localized topographic effects lead to significant heterogeneity, including notably thick snow at the glacier terminus (17% higher than the rest of the ablation zone), underscoring the inadequacy of using elevation or sparse in situ records to capture glacier-wide conditions. Comparisons with automatic weather station measurements suggest that the two stations in the ablation zone are situated in locally low-accumulation areas likely due to wind-driven snow depletion and are thus unrepresentative of snow accumulation over the surrounding region. These results emphasize the need for higher-resolution observational networks and modelling approaches that incorporate wind redistribution processes and local topographic complexity to properly resolve winter snow cover on peripheral glaciers which are commonly located in complex terrain. Nevertheless, evaluation against CARRA indicates that the accumulated precipitation variable provides a reasonable first-order estimate of the mean snow accumulation over APO, despite the high-sub grid variability in snow accumulation over areas with complex terrain.
7. Datasets
Raw and processed GPR data, along with the derived snow depths can be found on dataverse.geus.dk (https://doi.org/10.22008/FK2/02ALVH). GPR-derived snow depths can also be downloaded from https://data.g-e-m.dk. Snow pit density measurements are available on dataverse.geus.dk (https://doi.org/10.22008/FK2/9QEOWZ) (last accessed 4-July-2024). AWS sonic-ranger and wind speed data from the APO AWS are available through https://data.g-e-m.dk (last accessed 1-Mar-2025). Wind direction data is available on dataverse.geus.dk (https://doi.org/10.22008/FK2/AX8PXW). CARRA data was downloaded from the C3S Climate Data Store (CDS) (Schyberg and others, Reference Schyberg2020) (last accessed 03-Oct-2024). ArcticDEM (Mosaic, Version 4) (Porter and others, Reference Porter2023) was downloaded from https://www.pgc.umn.edu/data/arcticdem/ (last accessed 14-Nov-2024).
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.10080.
Acknowledgements
This research has been supported by the Danish Energy Agency (Climate and Environmental Support) within the framework of GEM Greenland Ecosystem Monitoring. We thank the reviewers and editor for their constructive and thoughtful feedback, which strengthened this manuscript.
Competing interests
None of the authors have any competing interests.
Open access
