Significant Spatial Variability in Radar‐Derived West Antarctic Accumulation Linked to Surface Winds and Topography

Across the Antarctic Ice Sheet, accumulation heavily influences firn compaction and surface height changes. Therefore, accumulation varies over short distances (<25 km), complicating the derivation of ice sheet mass changes from altimetry and reducing how accurately field measurements can be spatially extrapolated. However, current atmospheric reanalyses have grid spacings (>25 km) that are too coarse to resolve this variability. To address this limitation, we construct a fine‐scale accumulation product from airborne snow radar observations by superimposing along‐track fluctuations in accumulation onto an atmospheric reanalysis product. Our resulting airborne product reflects large‐scale (>25 km) orographic precipitation patterns while providing robust and unprecedented insight into Antarctic accumulation variability on subgrid scales. On these smaller scales, we find significant, regionally dependent accumulation variability (σrelative>40%). This variability in accumulation is correlated with variability in topographic surface slope in the wind direction (p<0.01), confirming that subgrid‐scale accumulation variability is driven by snow redistribution by wind.


Introduction
Since 1979, the West Antarctic Ice Sheet (WAIS) and the Antarctic Peninsula Ice Sheet (APIS) have contributed 6.9 ± 0.6 and 2.5 ± 0.4 mm to global sea level change (GSLC), respectively (Rignot et al., 2019). A grounded ice sheet's impact on GSLC is determined by its mass balance or the difference between accumulation onto its surface and flux of ice across its grounding line (Shepherd et al., 2018). Therefore, if we are able to better constrain accumulation over the WAIS and APIS, we can reduce uncertainty in their contributions to GSLC.
Over the WAIS and APIS, virtually all surface meltwater refreezes in place and meltwater runoff is negligible (e.g., Kuipers Trusel et al., 2013). Therefore, we are able to approximate accumulation over these two ice sheets as precipitation minus evaporation, P − E. Evaporation is almost entirely comprised of sublimation since nearly all liquid water produced at the surface of the WAIS and APIS instantaneously percolates and refreezes (Trusel et al., 2013).
These spatial variations in accumulation have direct implications for ice core and altimetry studies. Since ice cores are point measurements, one must consider local spatial variability in accumulation when extrapolating their observations over large areas. Accumulation rates at fine spatiotemporal resolution are also necessary for estimating temporal changes in firn compaction. Accurately modeling these firn compaction rates across the AIS is necessary to derive mass change from temporal changes in surface elevation measured by laser or radar altimetry (e.g., Sandberg Sørensen et al., 2011;Shepherd et al., 2018;Zwally et al., 2015). Small-scale spatial variations in accumulation are not resolved in these estimates of mass change; models of firn compaction are typically forced by data from regional climate models or reanalyses, whose resolutions are too coarse to account for the effects of snow redistribution by wind. An accumulation product that can match the resolution of current altimetry measurements (<5 km) would allow for more accurate estimates of mass change from surface height changes (Csatho et al., 2019).
The National Aeronautics and Space Administration Operation IceBridge (OIB) snow radar is a promising tool to capture small-scale variations in accumulation as it is able to track internal stratigraphy with an along-track sampling frequency of one trace per every 5 to 10 m. The airborne OIB campaigns have been collecting observations over large portions of the WAIS and APIS annually for the past decade (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019). Up until present, analysis of OIB snow radar derived accumulation rates over the AIS has been limited to the middle-to-upper Thwaites and Pine Island glacier catchments (Medley et al., 2013(Medley et al., , 2014, coastal ice domes (Lenaerts et al., 2017), and the Larsen C ice shelf (Kuipers Munneke et al., 2017).
In this study, we present a new snow accumulation product that represents a much more expansive area than that of previous studies. In addition to Thwaites and Pine Island glacier catchments, our product provides coverage of most of coastal West Antarctica and parts of the APIS, Ronne Ice Shelf, and South Pole region. The total length covered by all of our flight lines (17,500 km) is over seven times that of the most extensive previous study of snow radar-derived accumulation over the AIS (Medley et al., 2013).
In order to achieve this degree of coverage, we created an algorithm that is unique in its independence from ice cores and other in situ information. This independence allows us to apply the algorithm to snow radar observations covering any region of the AIS, regardless of its proximity to ice cores. To circumvent the need for external input, we initially assume that P − E from MERRA-2 is accurate for length scales over 25 km, as it performs best among other reanalysis products (Medley & Thomas, 2019). We then generate an accumulation product that reflects this reanalysis product, including all of its errors, on scales over 25 km. However, our accumulation product also resolves smaller-scale variations that we derive from snow radar data. To analyze these small-scale spatial variations in accumulation, we focus on a relative accumulation product that both reflects the subgrid-scale accumulation variation and is shown to be largely independent of MERRA-2's absolute accumulation field. In our analysis, we correlate variations in relative accumulation to variations in topography, wind speed, and wind direction. Lastly, we discuss the implications of these relationships on calculating accumulation over regions of the WAIS and APIS not surveyed by the airborne radar.

Absolute Radar-Derived Accumulation Rates
We estimate accumulation rates with unprecedented coverage of the WAIS and APIS at 100-m along-track spacing from 29 OIB snow radar flights between 2010 and 2017. We focus this analysis on along-track gradients in accumulation and on their regional variations. Therefore, we pick a subset of 29 surveys that spatially covers a large portion of these two ice sheets. Additionally, we exclude flight lines with heavily discontinuous internal reflection horizons (IRHs). These discontinuities result from either (a) rapid changes in aircraft altitude in response to steep surface topography or (b) disruptions in firn stratigraphy caused by reworking of the snow at the surface in low accumulation conditions. Given our data subselection and quality control, we create an expansive product that includes 175,373 individual accumulation measurements.
IRHs within the snow and firn are interpreted as isochrones (i.e., representative of a single time or event); therefore, the mass between horizons (cumulative mass) can be used to calculate accumulation rates (Eisen et al., 2005;MacGregor et al., 2009;Medley et al., 2013Medley et al., , 2014Medley et al., , 2015Sinisalo et al., 2003). In this study, we calculate accumulation rates from snow radar data using two isochrones: the surface and a single IRH.
We divide the snow radar data into 25-km along-track segments, tracking the surface and a continuous and clearly resolved IRH within the first 50 m of snow and firn. We use semiautomated layer picking software to track these horizons, visually verifying and manually retracking as necessary. Our software produces one layer pick per 100 m along track. To obtain a greater signal to noise ratio, we preferentially select deeper (closer to 50 m) IRHs wherever possible, as a horizon's along-track variation tends to amplify with depth.
By tracking the clearest IRH within the first 50 m of snow and firn, we average accumulation over multiple years. For each segment, this time span depends on the accumulation rate, depth of the tracked IRH, and the year of the OIB overflight. Although this methodology results in time spans that vary from segment to segment, tracking the clearest IRHs reduces uncertainty in layer picking. Moreover, we assume that spatial variations in accumulation are persistent in time for length scales between 100 m and 25 km (Medley et al., 2013). Given this assumption, we are able to track the best-resolved IRHs to quantify spatial variability in accumulation by averaging over multiple years.
To calculate accumulation rates, we first compute cumulative mass, or the integral of density with depth, between two horizons. At every 100 m along track, we record the two-way travel time (TWTT) between the surface and the IRH for a total of 250 layer picks per 25-km flight segment. At each of these 100-m-spaced layer picks, we generate a depth-density profile using a semiempirical, steady-state firn compaction model (Herron & Langway, 1980). We use mean P − E and 2-m air temperature from MERRA-2, in addition to initial snow density of 350 kg/m −3 based on ice core observations (see supporting information Figure S1), as input into this model. By calculating the depth-permittivity profile from the modeled depth-density profile (Kovacs et al., 1995), we are able to convert TWTT to depth. We integrate the depth-density profile from the IRH to the surface at each layer pick to find cumulative masses that vary from pick to pick along the segment. This method of calculating cumulative mass is described in detail in Medley et al. (2015).
In this ice sheet-wide study, we determine the age of an IRH using a methodology that differs from previous, region-specific studies (e.g., Medley et al., 2013Medley et al., , 2015, in which a layer's age was determined either from ice cores or by counting annual layers down from the ice sheet surface. Neither of these methods are feasible for our study, as ice core data are not available everywhere and firn stratigraphy cannot always be assumed to be annual. Instead, we tailor a layer's age to suit MERRA-2's conditions, finding a layer age that will force accumulation rates to converge to MERRA-2 P − E on segment-scale (25 km) by arithmetic mean.
We calculate the age of our single IRH since its original deposition, using the assumption that MERRA-2 P−E is accurate for the 25-km length for which we track this layer. We first derive depth-age profiles for every pick along a 25-km segment using the same MERRA-2-forced firn compaction model described above. Since the depth-age profiles are determined by P−E, these profiles do not change within the 25-km segment because a MERRA-2 grid cell is larger than each segment. In reality, the accumulation rates and the depth-age profiles both vary along track, ensuring that the age along the segment remains constant. Instead, the depth-age profiles forced by MERRA-2 are constant, so our calculations result in a range of ages based on our range of depths. Since we assume the layer to be isochronal and intend to force accumulation rates to MERRA-2 P−E on a large scale, we force the mean of this range of estimated ages to equal the MERRA-2-driven layer age.
The quotient of cumulative mass from an IRH to the surface, CM(x), and MERRA-2-driven layer age, age, provides an absolute accumulation rate for each layer pick, x, along each 25-km segment: As in Medley et al. (2015), we refeed this accumulation rate into the firn compaction model to create a new depth-density profile, a new cumulative mass, and a new output accumulation rate. We repeat this process until the output accumulation rate converges with the input accumulation rate (one to two iterations); thus, our cumulative masses and accumulation rates are self-consistent with the depth-density profiles from which they are calculated. This scheme results in very slight changes in along-track accumulation rates (<0.1%).
This process yields an absolute accumulation product in units of meters of water equivalent per year (m w.e. yr −1 ). This airborne product reflects reanalysis on the scale of a segment (25 km) but resolves 100-mto 25-km-scale accumulation variability. Since this absolute accumulation product is forced to MERRA-2 P − E on the length scale of a MERRA-2 grid cell, all errors in MERRA-2 P − E are reflected in our absolute accumulation product. For a more detailed analysis of the relationship between our absolute accumulation product and MERRA-2 P − E, see supporting information Text S1.

Relative Radar-Derived Accumulation
Relative accumulation, RA, is the quotient of the pick's cumulative mass, CM(x), and 25-km mean cumulative mass, CM: Relative accumulation is unitless, centered about one, and largely independent of MERRA-2, as relative accumulation is heavily dependent on variations in TWTT along an IRH. For that reason, errors in relative accumulation are found to be very small, with a spatial mean of 0.002 (See supporting information Text S1 for more details).

Mean Slope in Wind Direction
Mean slope in the (mean) wind direction (MSWD) is defined as the dot product of the topographic surface slope with the annual mean near-surface (10 m) wind direction (Scambos et al., 2012). In this study, topographic surface slopes are derived from the Antarctic CryoSat-2 1-km resolution digital elevation model (DEM; v.1.0) of surface topography (Helm et al., 2014), and wind vector components are based on annual mean (1980-2017) near-surface wind from the MERRA-2 reanalysis that is linearly interpolated onto the DEM grid.
We also quantify the impact of flight direction relative to topography and wind direction on accumulation variability. To do so, we calculate the alignment index, which we define as |û ·̂| * |ŝ ·̂|, where flight direction iŝ, topographic slope direction isŝ, and wind direction is û.

Radar-Derived Accumulation Product
On spatial scales exceeding 25 km, our radar-derived absolute accumulation product is forced by MERRA-2 P − E, which is shown in Figures 1a and 1b for reference. The spatial distribution of absolute accumulation within MERRA-2 grid cells is illustrated in Figures 1c and 1d. The resolution of absolute accumulation is significantly higher than that of MERRA-2 reanalysis; in the Amundsen Sea Embayment sector, several MERRA-2 grid cells contain over 1,000 radar-derived accumulation measurements each. In addition to its fine-scale resolution, our product represents a wide spatial extent of OIB flight lines over the WAIS and APIS, as demonstrated in Figures 1c through 1l. Regions covered range from coastal zones, where absolute accumulation rates are generally 1.0 m w.e. yr −1 or higher, to the high elevation interior around the South Pole, where absolute accumulation rates are 0.03 m w.e. yr −1 or lower (Figures 1e and 1f). These absolute accumulation rates reflect MERRA-2 P − E on the length of a MERRA-2 grid cell while representing spatial variations in accumulation on subgrid scales.
In Figures 1g and 1h, our relative accumulation product illustrates these smaller-scale variations with respect to a 25-km binned mean in the along-track direction. To quantify these variations, we calculate (e, f and g, h) Radar-derived absolute accumulation rates and relative accumulation, respectively. (i, j and k, l) The 25-km standard deviation of radar-derived absolute accumulation rates and relative accumulation, respectively. The second and fourth columns show the same data as the first and third columns, respectively, but centered on the Amundsen Sea Embayment sector, as indicated by the grey boxes in the first and third columns. Red box in (i) outlines the Eights Coast.
standard deviations of accumulation for each 25-km segment. Absolute accumulation rates are highly variable (Figures 1i and 1j) with some regions, such as the Eights Coast, showing standard deviations as high as 0.6 m w.e. yr −1 . Standard deviation in relative accumulation reveal variability irrespective of the coast-to-interior gradient in accumulation (Figures 1k and 1l). These standard deviations reach as high as 50% in the Amundsen Sea Embayment sector. The region around the South Pole also shows high spatial accumulation variability, with areas of standard deviations in relative accumulation reaching 40%.

Relationship Between Wind, Surface Topography, and Accumulation
Our relative accumulation product demonstrates a relationship between accumulation, wind, and surface topography. A 25-km flight segment inland of the Getz Ice Shelf grounding line, shown in Figure 2, illustrates this correlation. This flight line was chosen because of the parallel alignment of the near-surface wind and the OIB flight (<10 • ). The variable depth of the tracked horizon in Figure 2b indicates that accumulation rates are not uniform across this 25-km segment. Relative accumulation rates in Figure 2c reveal this along-track spatial variability. Comparing Figures 2c and 2d, we find that lower accumulation rates colocate with positive topographic surface gradients with respect to wind and higher accumulation rates colocate with negative gradients with respect to wind. Their inverse correlation is shown in Figure 2e (r = −0.69;  p < 0.01; n = 251). Note that all values of p are adjusted for spatial autocorrelation, as detailed in supporting information Table S1. This relationship between wind and topographic surface slope indicates that wind is (1) scouring snow from exposed topography with positive slopes with respect to wind direction and (2) redepositing snow downwind, onto negative slopes with respect to wind direction.
The effects of wind and slope on accumulation can be examined in both horizontal dimensions using gridded MSWD. MSWD, the dot product of wind vectors and topographic surface slopes, is related to small-scale accumulation variability across the ice sheet. In Figures 3a and 3b, high variability in MSWD underlays highly variable accumulation. Likewise, decreased variability in relative accumulation is colocated with near-zero MSWD. Relative accumulation and MSWD do not exhibit a direct correlation with each other, but are related through their spatial patterns in variability.
We quantify this relationship by gridding and comparing standard deviation in MSWD ( MSWD ) to standard deviation in relative accumulation ( RA ) in Figure 4. MSWD and RA are correlated at a 99% confidence interval with correlation coefficients of 0.21 and 0.28 for grid cell spacings of 10 and 25 km, respectively (n = 3,496, Figure 4e; n = 516, Figure 4f). For alignment indices above 0.4, MSWD and RA were correlated at or above a 98% confidence interval, with slightly higher correlation coefficients of 0.24 and 0.44 for grid cell spacings of 10 and 25 km, respectively (n = 1,558, Figure 4e; n = 232, Figure 4f). Additionally, extreme values of MSWD and RA are also related. For both the upper and lower tenth percentiles, MSWD and RA are colocated at a rate of 85% ± 1%. For more detail about our statistical testing, see supporting information Table S1.

Discussion and Conclusions
This study should be regarded as an initial step toward providing a high-resolution, gridded accumulation product over the AIS. We found a relationship between variability in accumulation and surface topography that could be used to inform the downscaling of a reanalysis accumulation product to the grid of a DEM (1 km). Currently, the relationship between relative accumulation and MSWD has limitations. First, we assumed that the near surface wind climatology in MERRA-2 (0.5 • by 0.625 • ) downscales linearly to a 1 km resolution to compute MSWD. Second, by using a wind climatology to derive MSWD, we ignored seasonality, past variability, and trends in AIS wind climate. However, the near-surface winds on Antarctica (especially in the interior) are predominantly of katabatic nature (Parish & Cassano, 2003), very constant in direction (van den Broeke & van Lipzig, 2003), and forced by downslope surface topography. We therefore assume the AIS near-surface wind field and associated MSWD grid to be relatively persistent and expect little changes in the general AIS wind field in the next decades.
We highlight a few avenues to fine-tune the relationship between MSWD and relative accumulation. MSWD could be recomputed using a new higher resolution (31 km) reanalysis product, European ReAnalysis-5 (ERA-5), if the product is shown to compare better with Antarctic observations than MERRA-2. Furthermore, we could obtain a more accurate MSWD grid by downscaling reanalysis wind fields using an algorithm that modifies wind speed and direction to account for the effects of topography (Liston & Sturm, 1998). Moreover, further refinement of the topography can be obtained using higher-resolution elevation models (Howat et al., 2019;Shean et al., 2016). Additionally, the timing of wind events with respect to surface snow density could be used to modulate MSWD as denser snow is less likely to erode. If a strong inverse correlation between MSWD and relative accumulation is found, a high-resolution gridded (1 km) accumulation product across the AIS could be downscaled from reanalysis based on this relationship.
In our study, we constructed an airborne accumulation product across the WAIS and APIS by using both the snow radar's high along-track resolution and OIB flights' wide spatial coverage. The fine-scale resolution of the instrument allowed us to quantify variability on length scales of 100 m to 25 km. The magnitude of this along-track variability in accumulation rates is nonuniform over the ice sheets. Most notably, high spatial variability of accumulation occurs in the Amundsen Sea Embayment sector, with standard deviation of relative accumulation reaching 50%. Across our data set, we found a significant correlation (p < 0.01) between MSWD and RA . The relationship between extreme values of MSWD and RA indicates that MSWD can be used as a threshold index for RA . Overall, our assessment shows that small-scale (<25 km) variations in accumulation are significant, regionally dependent, and driven by wind redistribution of snow.