Temporal and spatial variability in surface roughness and accumulation rate around 88° S from repeat airborne geophysical surveys

We use repeat high-resolution airborne geophysical data consisting of laser altimetry, snow and Ku-band radar and optical imagery acquired in 2014, 2016 and 2017 to analyze the spatial and temporal variability in surface roughness, slope, wind deposition, and snow accumulation at 88° S, an elevation bias validation site for ICESat-2 and potential validation site for CryoSat2. We find significant small–scale variability (< 10 km) in snow accumulation based on the snow radar subsurface stratigraphy, 15 indicating areas of strong wind redistribution are prevalent at 88° S. In general, highs in snow accumulation rate correspond with topographic lows resulting in a negative correlation coefficient of r2 = -0.32 between accumulation rate and MSWD (Mean Slope in the mean Wind Direction). This relationship is strongest in areas where the dominant wind direction is parallel to the survey profile, which is expected as the geophysical surveys only capture a two-dimensional cross section of snow redistribution. Variability in snow accumulation appears to correlate with variability in MSWD. The correlation coefficient between the standard 20 deviations of accumulation rate and MSWD is r2 = 0.48 indicating a stronger link between the standard deviations than the actual parameters. Our analysis shows that there is no simple relationship between surface slope, wind direction and snow accumulation rates for the overall survey area. We find high variability in surface roughness derived from laser altimetry measurements on length–scales smaller than 10 km, sometimes with very distinct and sharp transitions. Some areas also show significant temporal variability over the course of the 3 survey years. Ultimately, there is no statistically significant slope–independent relationship 25 between surface roughness and accumulation rates within our survey area. The observed correspondence between the small–scale temporal and spatial variability in surface roughness and backscatter, as evidenced by Ku-band radar signal strength retrievals, will make it difficult to develop elevation bias corrections for radar altimeter retrieval algorithms.


Introduction
Polar ice sheets play a critical role in Earth's climate system. Measurements from satellites and aircraft reveal that the ice sheets 30 of Greenland and Antarctica are changing at an accelerating rate suggesting increasing rates of global sea-level rise as the ice sheets melt (e.g., Vaughan et al., 2013). To project future rates of sea-level rise, numerical models of an ice sheet's response to climate forcing require input data of the surface mass balance and its spatial and temporal variability. Observing changes in icesurface elevation from satellite and airborne platforms has long been recognized as a powerful tool for assessing and quantifying ice sheet mass balance (e.g., Abdalati et al., 2010;Krabill et al., 2000;Thomas and Investigators, 2001;Zwally et al., 2002). 35 The southern convergence of all Ice, Cloud, and land Elevation Satellite-2 (ICESat-2, (Markus et al., 2017)) and CryoSat-2 (Wingham et al., 2006) ground reference tracks at 88° S (McConnell et al., 1997;Mosley-Thompson et al., 1999;Winski et al., 2019, Helm et al., 2014. Because of the density of tracks, the small impact of surface processes, and the region's relative quiescence, 88° S is the primary ICESat-2 land-ice validation site in the southern hemisphere (Brunt et al., 2019a;Brunt et al., 2019b). Both radar and laser altimeters are potentially prone to elevation biases related to surface roughness and slope. For laser 40 altimeters such as ICESat-2, increased surface roughness causes broadening of the return signal, which can cause elevation biases up to 0.2 m (Smith et al., 2019). When surface roughness changes seasonally the elevation biases will also change with time (Smith et al., 2019). For radar altimeters such as CryoSat-2, smoother surfaces will have larger return signal strength compared to rougher surfaces which also changes the shape of the return waveform potentially causing elevation biases (Kurtz et al., 2014). Since radar altimeters penetrate below the ice surface, volume backscatter from subsurface firn will also impact the return signal waveform 45 (e.g., Nilsson et al., 2016). Radar extinction with depth depends on the dielectric permittivity of firn, which is primarily a function of firn density (Kovacs et al., 1995). Changes in firn density are often related to changes in snow accumulation rates (e.g., Grima et al., 2014) making radar elevation biases potentially a function of spatial changes in accumulation rates. Furthermore, windinduced anisotropic features of firn can introduce azimuth dependent elevation biases (Armitage et al., 2014).
Previous Antarctic studies have reported relationships between surface slope, roughness and snow accumulation rates (e.g., Arcone 50 et al., 2012;Dattler et al., 2019;Fahnestock et al., 2000;Grima et al., 2014;Hamilton, 2004;King et al., 2004). To better understand potential correlations between altimetry elevation biases and geophysical parameters of the ice surface, we are specifically studying the spatial and temporal variability of surface roughness and accumulation rate over the ICESat-2 validation site at 88° S. We use repeat high-resolution airborne data consisting of laser altimetry, snow and Ku-band radar and natural color imagery acquired as part of the National Aeronautics and Space Administration's (NASA) Operation IceBridge (OIB) mission to analyze spatial and 55 temporal variability in surface roughness, slope, accumulation rate, and Ku-band radar backscatter along a 1400 km circle around 88° S (Fig. 1). We start with a description of the survey area and the data sets we use (Section 2). We then focus on surface roughness in Section 3, first describing its relevance to cryospheric research, the methods we use to estimate surface roughness, followed by a description of temporal and spatial variations in surface roughness that we observe in our data. We then describe in Section 4 how we estimate snow accumulation rates, followed by a description of the spatial variability we observe in snow 60 accumulation rates in our survey area in Section 5. Section 6 explores the relationship between surface, slope and wind direction.
The impact of surface roughness on radar backscatter and therefore elevation bias in radar altimeter measurements is discussed in Section 7. Section 8 concludes the paper.

Data sets and survey area
Our survey area is situated on the East Antarctic plateau in the hinterland of the Transantarctic Mountains (Fig. 1). The ice surface elevation along the 1400 km long survey line around 88° S varies between 2450 and 3100 m, with bedrock elevations ranging from -1450 to 1785 m (Morlighem et al., 2020) and Fig. 4. The thinnest ice along 88° S is 1190 m thick and the thickest ice reaches 70 4100 m (Fig. 4). Ice surface velocities in our survey area are generally below 10 m yr -1 (Mouginot et al., 2019). Accumulation measurements from snow pits and shallow firn cores are sparse (Favier et al., 2013;Picciotto et al., 1971). The survey area is in a region of low snow accumulation (7-10 cm annual water equivalent) (Arthern et al., 2006;McConnell et al., 1997;Mosley-Thompson et al., 1999;Winski et al., 2019) and low surface slope (0.11° ± 0.10°, Fig. A1a) (Helm et al., 2014). Together with the low snow accumulation rate and low ice surface velocities this makes it an ideal area for calibration and validation of spaceborne 75 altimeters (Brunt et al., 2019a;Brunt et al., 2019b).
Accumulation of snow on the Antarctic ice sheet is primarily the result of precipitation of snow. The precipitated distribution of accumulated snow is subsequently modified spatially by wind-driven erosion and deposition. Sublimation of accumulated snow, both, in the form of wind-driven sublimation of airborne snow particles and surface sublimation removes accumulated snow and therefore mass from the surface and further modifies the initial deposition pattern resulting from precipitation (e.g., Frezzotti et al., 80 2007, and references therein). For slopes ≥ 0.002 Das et al. (2013) found wind-scoured areas in East Antarctica with negative surface mass balance similar to the wind-glaze area described by Scambos et al. (2012). Using European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data, Casey et al. (2014) estimate that around half of the snow accumulation at the South Pole comes from periodic moisture-bearing storms traversing the Filchner-Ronne and Ross Ice Shelves towards the pole from West Antarctica. It is likely that the percentage of snow accumulation from such cyclonic events is even higher at 88° S near 85 the Transantarctic Mountains compared to the South Pole since the area is closer to the source of moisture, complicating the relationship between surface slope, wind direction and snow accumulation. Therefore, there is no simple relationship between surface slope, wind direction and snow accumulation rates.
We use airborne geophysical data collected during 6 NASA OIB survey flights in 2014, 2016 and 2017. The data consists of highresolution laser altimetry, natural color imagery and snow and Ku-band radar data and is available from the National Snow and Ice 90 Data Center (NSIDC).

Airborne Topographic Mapper (ATM)
The ATM is a conically-scanning laser altimeter that measures the surface topography of a swath beneath the aircraft at a 15° offnadir angle (Krabill et al., 2002). The range from the laser altimeter to the surface is converted to geographic position by integration 95 with platform Global Positioning System (GPS) and attitude/Inertial Measurement Unit (IMU) measurement subsystems. The conical scan geometry results in a near-constant angle of incidence and intersecting laser footprints allow for pointing biases to be determined over any type of surface (Harpold et al., 2016;Martin et al., 2012). The two generations of instruments used in this study, T4 and T6 have a pulse repetition frequency of 3,000 Hz, a wavelength of 532 nm, and a pulse width of 6 ns full width at half maximum (FWHM). The ATM scanner has a swath width of 240 m at a nominal flight elevation of 460 m above ground level 100 (AGL) and a footprint diameter of ~0.8 m. As a result of the conical scan pattern, the density of spot elevation measurements varies across the swath from 0.03 footprints m -2 at the center to 0.37 footprints m -2 at the edge. At a nominal aircraft speed of 130 m s -1 the average spacing between point elevation measurements is ~5 m in the center of the scan and <1 m near the edge. The vertical accuracy of an individual laser spot measurement is estimated to be 7 cm with a vertical shot-to-shot precision of 3 cm (Martin et al., 2012). We use both the L1B and L2 (ICESSN) data products (Studinger, 2013(Studinger, , updated 2018(Studinger, , updated 2018. 105

Riegl LMS-Q240i laser altimeter
The University of Alaska, Fairbanks (UAF) operates a commercially available Riegl LMS-Q240i airborne laser scanner together with IMU and dual-frequency GPS subsystems for attitude and precise position. The system is a near-infrared linear, unidirectional scanner that scans the surface in parallel lines. The system acquires measurements at 10,000 Hz with a footprint size of 1.0 m at 460 m AGL and at a ±30° off-nadir scan angle. The average spacing of laser footprints both along track and across is ~1 m at 460 110 m AGL and a ground speed of 85 m/s (Johnson et al., 2013;Larsen, 2010Larsen, , updated 2018. Results over 20% of our study area at 88° S show that 2 flights from a 2017 UAF laser altimetry survey had a <10 cm bias and a surface measurement precision of <10 cm (Brunt et al., 2019a).

Snow and Ku-band radar
The snow radar is an ultra-wideband microwave radar that operates over the 2-8 GHz frequency range and is developed and 115 operated by the Center for Remote Sensing of Ice Sheets (CReSIS) at the University of Kansas (KU). The system is a frequencymodulated continuous-wave (FMCW) radar that images the stratigraphy in the upper ~40 m of the ice sheet with a bandwidthlimited range resolution of 2.6 cm in firn with a density of 550 kg/m 3 and 3.8 cm in air (Leuschen, 2014(Leuschen, , updated 2018Panzer et al., 2013;Rodriguez-Morales et al., 2014). At these low accumulation sites, preservation of reflection horizons is greatly reduced due to a slow rate of burial. Also, the ambient conditions required to generate seasonal reflections might not always be present 120 (such as depth hoar). At the nominal flight elevation of 460 m AGL the snow radar has a footprint size of approximately 10 m across track and 14.5 m along track (Panzer et al., 2013). A 9-by-4 (range bin-by-trace) median filter is applied to the data to minimize noise. The Ku-band radar altimeter is identical in design but operates over the frequency range of 12-18 GHz for mapping subsurface stratigraphy in the upper 10 m of polar firn (Gomez-Garcia et al., 2012;Paden et al., 2014Paden et al., , updated 2018Rodriguez-Morales et al., 2014). Since both radars have the same bandwidth (6 GHz) the bandwidth-limited range resolution of the Ku-band 125 radar is the same as the snow radar.

Digital Mapping System natural color imagery
The Digital Mapping System (DMS) is a digital camera that acquires natural color, high-resolution images at 10 cm pixel size at the nominal flight elevation of 460 m AGL (Dominguez, 2010(Dominguez, , updated 2018 (Fig. 2). The camera is operated by NASA's Airborne Sensor Facility located at the Ames Research Center. Images are approximately 380 m across swath and 570 m along swath and 130 cover the entire ATM swath width. A combined IMU and GPS system for precise position and attitude information is part of the instrument package. DMS images are acquired with overlap between consecutive images to ensure data continuity. The difference in geolocation between distinct elongated snow surface features (sastrugi) between overlapping, orthorectified images is on the order of several meters. The DMS images are referenced to the RADARSAT 200 m Digital Elevation Model (DEM) (Liu et al., 2015). 135

Survey flights 145
Between 2014 and 2017 six NASA OIB airborne geophysical survey flights were completed to acquire data around 88° S (Table   1). Two flights are necessary to complete the entire small circle around 88° S with the platforms and bases of operations used.
Snow and Ku-band radar data and DMS images are only available for 2014 and 2016 years. The combination of simultaneous laser altimetry, snow radar stratigraphy and natural color imagery on a regional scale provides a unique data set to study small scale deposition and erosional processes and their temporal and spatial variability. 150  3 Ice surface roughness

Background
The surface roughness of polar ice sheets is primarily a result of ice dynamics and surface-atmosphere interactions on varying 155 temporal and spatial scales. In general, ice flow over rugged bedrock topography causes roughness features that can extend from several hundreds of kilometers to a few kilometers depending on ice thickness, flow speed and basal conditions (e.g., , and references therein). These large-scale variations in ice surface topography caused by ice dynamics are not the topic of this analysis. Here, we focus on small-scale surface roughness or surface texture that spans from several meters to hundreds of meters and is primarily the result of ice-atmosphere interactions, such as wind deposition and wind-induced ablation or erosion, 160 the predominant types of surface-atmosphere interactions in the area of 88° S. Sastrugi are the dominant form of small-scale surface roughness in the interior of polar ice sheets and are known to form parallel to the prevailing wind direction.Their orientation can therefore be used to infer time-averaged prevailing wind directions (e.g., Bromwich et al., 1990;Gow, 1965). Lister and Pratt (1959) describe sastrugi on the order of 20-30 meters along the route of the Commonwealth Trans-Antarctic Expedition that also crossed 88° S. Figure   In recent decades, especially in the rapidly warming West Antarctic region, synoptic heat-and moisture-bearing storms have reached the South Pole area (e.g., Harris, 1992;Nicolas and Bromwich, 2011). Such storms and cyclonic events have been less 180 common, though still occur in the interior of East Antarctica (e.g., Gorodetskaya et al., 2014;Hirasawa et al., 2000). Individual sastrugi can be eroded and reform during a single storm (Warren et al., 1998). Most changes, however, occur as a result of seasonal changes between summer and winter months (e.g., Gow, 1965). Gow (1965) shows that sastrugi form during winter months resulting in a rough surface, and subsequently get eroded during the summer by sublimation and deflation at South Pole. This effect mostly results in a flattening of the subsurface layer stratigraphy and therefore does not affect our surface roughness results. 185

Relevance of surface roughness and slope for altimetry and surface mass balance
The surface roughness and slope of ice sheets affect several processes that are relevant for ice sheet mass balance (e.g., Gow, 1965, and references therein; van der Veen et al., 2009). King et al. (2004) describe small scale variations in accumulation rate on the order of 1 km that appear to be associated with wind-borne redistribution as a function of slope. Hamilton (2004) found significant variability in snow accumulation rates due to the interaction of prevailing winds with meter-scale surface topography; where, for 190 example, a concave depression can receive up to 18% more accumulation than adjacent steeper snow surface topography. Similarly, Arcone et al. (2012) mapped accumulation patterns in East Antarctica that are created by wind-blown deposition on windward and leeward slopes. Slope-dependent accumulation is also related to spatial variations in firn density (Grima et al., 2014) which impacts mass balance estimates from altimetry data. Small-scale roughness contributes to noise in firn core records and therefore accumulation rate estimates (van der Veen et al., 1998;van der Veen et al., 2009). Studies by van der Veen et al. (1998;2009) 195 used ATM roughness estimates over Greenland to determine the uncertainty in water equivalent (w.e.q.) accumulation estimates from shallow firn cores.
Surface roughness also affects the albedo and bidirectional reflectance distribution function of ice sheets (Leroux and Fily, 1998;Warren et al., 1998). Nolin et al. (2002) used ATM roughness estimates to calibrate Multi-angle Imaging SpectroRadiometer (MISR) roughness estimates, and Nolin and Payne (2007) derived then relationships between ice surface roughness and near-200 infrared albedo using ATM and MISR data. Ongoing satellite and modelling investigations on radiative impacts of surface roughness and sastrugi continue to illuminate angular relationships and parameterizations that can be key to quantifying bidirectional reflectance distribution function and albedo sensitivities in ice surface studies (e.g., Corbett and Su, 2015;Kokhanovsky and Zege, 2004;Larue et al., 2019). As ice sheet surface roughness mapping and modeling capabilities improve, it will be possible to more accurately include the radiative effects of surface roughness. Surface roughness furthermore affects 205 thermodynamic fluxes because it affects boundary layer processes through the aerodynamic roughness and therefore the surface energy balance (Boisvert et al., 2017;Chambers et al., 2019;Nolin and Mar, 2018;Palm et al., 2017).

Surface roughness estimates -methods
There are several diverse approaches to quantifying topographic irregularity or surface roughness (e.g., Smith, 2014, and references therein). In general, roughness metrics are not only scale and orientation-dependent, but also impacted by the spatial resolution, 210 footprint size and sample spacing of the input data. One commonly used metric for surface roughness is the standard deviation σ of small-scale elevation fluctuations from a mean or de-trended surface in a given area or over a given length (e.g., Das et al., 2013;Smith, 2014, and references therein). In order to minimize potential effects from anisotropy in surface roughness over different length scales and orientations we have calculated surface roughness over an area roughly square in size common to both laser altimeters. Individual spot elevation measurements are binned into 0.06° longitude segments (240 m in length at 88°S). We 215 then fit a 3 rd order polynomial regression model through all spot elevation measurements within a longitude segment. We define the standard deviation σ of the residuals as a metric for surface roughness.
For estimation of surface roughness from snow radar data, we pick an initial surface for each radar trace record (every ~5 meters along track) by finding the maximum slope in the radar return power across 9 discrete time intervals called range bins. The size of a range bin in firn with a density of 550 kg/m 3 is 1.8 cm. Starting at the initial surface pick, we keep sliding the surface pick one 220 range bin deeper (or later in time) while the slope for the range bin remains above 3 standard deviations of the mean, which provides our final surface. Next, surface picks that lie outside of a 15-range bin window from a smoothed surface are discarded and set to the smoothed surface, however, very few data points are discarded (1.8%). Surface roughness is estimated from residuals to the smoothed surface fit. Specifically, surface roughness for a given location is calculated as the standard deviation in surface rangebin residuals for locations within a 250-meter radius. This radius was selected to ensure consistency with laser-altimeter-derived 225 roughness values. Finally, the range-bin roughness is converted to heights by using the radar wave velocity in air (2.998E 8 m/s).
For a closer look at very fine-scale spatial and temporal changes in surface roughness around 88° S we use roughness estimates contained in the ATM Level 2 smoothed ice surface data product, known as ICESSN s. The ICESSN data product includes slope and roughness estimates in overlapping 80 m × 80 m platelets across the swath (Studinger, , updated 2018. The root mean square of the residuals of a plane fit through the platelets is an estimate of the surface roughness. Removing the mean results in the 230 root mean square being equivalent to the standard deviation σ. Figure 4 shows the surface roughness estimates around 88° S latitude from three different instruments and over the course of 3 different years. In general, there is good agreement between the roughness estimated using the ATM laser altimeter and the synchronous roughness estimates from the snow radar (Fig. 4 b, c). Because of the smaller footprint size and higher sampling density the ATM laser-derived roughness estimates are slightly larger than the radar estimates, with few exceptions. The radar 235 estimates also reveal more scatter, probably caused by the much lower range resolution of the radar compared to the ATM laser altimeter. The mean difference between the laser minus snow radar roughness estimates is 0.6 ± 1.9 cm for 2014 and 1.6 ± 2.1 cm for 2016; however, the spatial patterns, which are of main interest for this study, are nearly indistinguishable. There are no obvious spatial patterns in the roughness difference between laser and radar that would reflect a geophysical signal. Because of the higher point density, roughness derived from the UAF Riegl system is larger than the ATM roughness estimates. The mean difference 240 between the 2017 UAF minus 2016 ATM laser estimates is 1.2 ± 2.4 cm.
A 370 km long segment between 150°W and 100°E was repeated in 2017 within 3 days with the same instrument. It is unlikely that the surface was significantly altered within 3 days and therefore the difference between the two estimates can be used as an approximate estimate of the instrument-specific accuracy and precision of the laser-derived roughness estimates. The mean surface roughness for the first east-bound flight is 9.5 ± 1.5 cm and 10.1 ± 2.0 cm for the later west-bound flight. The mean μ of the 245 difference in surface roughness between the two 2017 flights is 0.02 cm with a σ of 0.7 cm. For comparison, a separate study by Das and others (2013) over Dome Argus with a Riegl LMS-Q240i scanner show a similar range of roughness as our measurements around 88° S.

Surface roughness, slope and elevation around 88° S -results
In order to distinguish ice surface roughness features caused by ice dynamics from roughness features that are a result of iceatmosphere interactions, the proximity of the ice surface to bedrock topography and bedrock roughness must be understood (Fig.   4 a). The large ice thickness (Fig. 4) and slow ice flow velocities (< 10 m yr -1 (Mouginot et al., 2019)) in the survey area, combined with the small window size we use to calculate the roughness make it unlikely that any of the roughness features we observe are 260 ice dynamic related. Thus, we interpret the roughness characteristics shown in Figs. 3 and 4 as caused by ice-atmosphere interactions. Figure 3b and 4 show several spatially coherent segments with distinct roughness characteristics around 88° S that appear not to be related to ice dynamics. In general, the surface roughness estimated from snow radar and laser altimetry data varies between 2 and 25 cm. The smoothest surface in 2016 and 2017 is between 175° W and 60° E and includes Titan Dome. The smooth segment 265 also coincides with the highest ice surface elevations and shallowest surface slopes around 88° S (Fig. 4).
The segment between 70° W and 100° W shows a pronounced increase in roughness in 2014 (Fig. 4b). The MODIS Mosaic of Antarctica (MOA, (Haran et al., 2014)) shows that this segment is near the edge of a megadune field that is mostly north of 88° S (Fig. 3b). Megadunes are long-wavelength surface ripples (Fahnestock et al., 2000) with amplitudes on the order of a few meters (peak to trough) and wavelengths of several kilometers (Fahnestock et al., 2000;Scambos and Fahnestock, 1998). The typical 270 elevation pattern of megadunes is not visible in the 88° S laser altimetry data. The likely reasons for this are that the airborne geophysical data was collected at the edge of the dune field and the orientation of the dune crests is subparallel to the airborne geophysical data. The orientation of dune crests between 70° W and 100° W is approximately perpendicular to the prevailing surface wind direction from MERRA-2 (Fig. 3a) consistent with findings from Fahnestock et al. (2000).
A second megadune field can be seen between 130° W -145° W and 150° W -155° W (Fig. 3 and Fig. 4). The surface of this 275 dunefield is less rough in 2017 compared to 2014 and 2016. . Data for MOA was collected between 11/2013 and 03/2014. The temporal stability of megadune fields remains poorly understood. Fahnestock et al. (2000) found 60 m of dune migration over a 34 year period for a megadune field in the vicinity of Vostok Station, approximately 1100km from the dune field discussed here. Therefore dune migration rates could vary significantly between the two sites. Since our survey area is near the edge of the dune field we cannot rule out that over the course of 4 years the edge of the dune field has migrated out of the coverage of the airborne 280 geophysical data.
The slope of the ATM ICESSN nadir platelets shows many distinct peaks that are aligned very well between the 2014 and 2016 data, indicating that these features are stable in location (Fig. 4e). The mean μ and standard deviation σ of the surface slope around 88° S is 0.20° ± 0.16° and never exceeds 1.5°.

Temporal changes in surface roughness -results 285
We use ATM Level 2 ICESSN roughness estimates for a closer look at multiyear temporal changes in surface roughness around 88° S (Fig. 5) because they are calculated over 80 m × 80 m platelets. The 2014 data reveals several areas where surface roughness doubles over very short spatial scales of only a few hundred meters (Fig. 4a). These features, labelled A, B and C in Fig. 4b, are several tens of kilometers wide and appear to be oriented parallel to the main sastrugi direction visible in simultaneously collected ATM spot elevation data and Digital Mapping System (DMS) imagery. Fig. 5 shows a close-up of one of the features (C). The 290 rougher surface features are also present in the simultaneously collected CReSIS snow radar data (Fig. 5 a, c).
These areas of increased surface roughness disappear in 2016 or seem to be significantly reduced in amplitude with the sharpness of the edges significantly blurred. Both the laser derived surface roughness and the roughness estimated from snow radar data seems to be slightly lower than the smooth area. In 2017, the segments labelled A, B and C appear to have no distinct roughness anomaly (Fig. 4d) compared to the surrounding areas. 295 Figure 5: Laser spot elevation measurements and DMS imagery over roughness feature "C" (Fig. 4b)

can be seen in 2014 that is visible in both, the laser derived surface roughness, and the length of the shadows in DMS imagery (b). The roughness doubles over a distance of ~200 m (a). The orientation of the boundary is parallel to the dominant sastrugi orientation (b). In 2016 and 2017 the distinct change in roughness seems to have been smoothened out. 4 Snow accumulation rates derived from snow radar data and MERRA-2 -methods
For snow-radar derived accumulation calculations, we first stack traces to an approximate along-track separation of 100 meters (~ 310 18 traces), which largely reduces noise in the return power especially at depth. These stacked echograms are then combined into segments of ~100 km. A single radar reflection horizon, assumed isochronous, is tracked through the 100-km segment. The actual horizon picked will likely vary from segment to segment because we chose to map the strongest and most continuous reflection within that segment. The horizon is picked semi-automatically. First, the user visually selects the horizon of interest. The range bin with the strongest return power within a 15 bin window is then selected as the horizon "pick" for that trace. That pick is 315 extended laterally across all traces by finding the strongest return power in adjacent traces within the 15 bin window. The user can then modify the picks if they deviate from their visual interpretation. The user can also eliminate portions of a given horizon if visual inspection deems horizon differentiation impossible. The spatial variability in accumulation rates and the varying strength of signal return prevent the calculation of temporally consistent accumulation rate from a single continuous horizon around the entirety of 88°S.. Thus, the accumulation rates estimated for each 100-km segment will span differing time intervals; however, 320 because of the relatively low accumulation rates, the majority span several decades minimizing the impact of interannual variability.
For each 100-km segment, we estimate the spatial variability in snow accumulation using the aforementioned horizon picks.
Typically, radar derived accumulation rates rely on knowledge of the horizon age as well (e.g., Medley et al., 2013;, but a lack of nearby dated ice-core stratigraphy or clearly defined annual horizons restricts our ability to assign an age to our horizon 325 picks. Because our work is focused on evaluating the spatial variability in snow accumulation, we develop a method that approximates the age of a given horizon through combination of horizon depths and MERRA-2 mean annual precipitation-minusevaporation (P-E). In such a manner, our large-scale mean accumulation rates are forced to large-scale MERRA-2 P-E, however, rates are allowed to vary on < 1 km length-scales from our radar horizon picks. We detail the methodology below.
Assuming each horizon pick within a given segment is isochronous, we need to determine a way to approximate the age of that 330 horizon. To do so, we begin by determining the mean accumulation rate and 2-meter air temperature from MERRA-2 over the entire segment, and use those variables to model steady-state firn density and age profiles using Herron & Langway (1980). The two-way travel time (τ) between the surface and the horizon pick is converted into depth (d) assuming where c is the speed of light and ε is the integrated dielectric permittivity of the material above the horizon and x is the distance 335 along flight line. Specifically, we use the modelled depth-density profile to generate depth-dielectric permittivity based on Kovacs et al. (1995), which is then used to relate depth and two-way travel time with a depth-varying radar-wave velocity. Using this model, we interpolate horizon two-way travel time to depth. Depths vary along-track (as our layer pick varies), but our largescale depth-age model from Herron & Langway (1980) does not; thus, we will estimate a variable along-track age of the radar horizon. The use of a single firn density model for the entire 88° S circle is justified because the difference between the minimum 340 and maximum accumulation rates is relatively small. The along-track variability is counter to our initial assumption that the radar horizons are isochronous; however, when we take the average age along-track, we effectively force the overall mean accumulation rate for the segment to the large-scale MERRA-2 mean. We then use this age to calculate spatially varying accumulation rates along the entire segment as outlined by Medley et al. (2015). In such a manner, we force the large-scale mean accumulation rates to those prescribed by MERRA-2 but allow for small-scale variability derived from the snow radar horizon picks in the absence 345 of independent estimates of firn depth-age profiles (Dattler et al., 2019). For comparison we have plotted all existing accumulation measurements of Favier et al. (2013) in Fig. 6c over our MERRA-2 and radar-derived accumulation rates. These snow pit measurements include data from the 1962-1963 South Pole Traverse (Taylor, 1971. While there is general agreement it should be pointed out that Favier et al. {, 2013 #109) applied the quality rating of Magand et al. (2007), which identifies all snow pit data points shown in Fig. 6c as low quality and subsequently excludes these data points from the quality controlled version presented 350 in Favier et al. (2013). Further limitations of the comparison are the long time between the snow pit measurements and airborne data and the large variability in snow accumulation rates on length scales of 10 km that can be seen in the radar-derived snow accumulation rates.

Spatial variability in snow accumulation rates -results
Previous work used the Mean Slope in the mean Wind Direction (MSWD) for studying relationships between surface slope and 355 spatial variability in snow accumulation rates (e.g., Das et al., 2013;Dattler et al., 2019;Scambos et al., 2012). MSWD is defined as the scalar dot product between the surface slope with the mean wind direction (Scambos et al., 2012). Here, we use the timeaveraged zonal and meridional wind components u and v from MERRA-2, transformed in to a Cartesian polar-stereographic projection, to calculate the mean wind direction. Analyzing relationships between surface slope, accumulation rates, and mean wind direction at 88° S is limited by the latitudinal resolution of the MERRA-2 reanalysis model, which is 0.5° or 55 km, as well 360 as the cross-sectional nature of the geophysical surveys (i.e., the data represent a 2-dimensional cross section). Given the narrow swath width of the ATM laser data (240 m) we use the ice surface slope derived from the CryoSat-2 DEM (Helm et al., 2014) at 1 km resolution to calculate the MSWD. The slope south of 88° S is only weakly constrained due to the absence of elevation data imposing further limitations on the analysis. The MERRA-2 26 year 10 m wind field is interpolated to the CryoSat-2 DEM grid cell locations. The difference in spatial resolution between the surface DEM and MERRA-2 will result in MSWD uncertainty from 365 oversampling the wind field. Because of the small slopes in the study area, however, we don't anticipate complex wind fields where actual wind orientation would significantly deviate from the MERRA-2 reanalysis model. Small topographic features, however, are not represented by the 10 m surface wind field as will be discussed later. In general, annual snow accumulation is around 5 cm w.e. yr -1 and is highest near the backside of the Transantarctic Mountains 375 near 150° W, a region that is influenced by precipitation from cyclonic events penetrating the area from the Bellingshausen and Ross Sea sectors (Casey et al., 2014;Nicolas and Bromwich, 2011) (Fig. 6). The highest accumulation rates (up to 20 cm w.e. yr -1 ) near 150° W coincide with a megadune field (Fig. 3b) and appear to be in a local topographic low at the flank of Titan Dome that trends perpendicular to the aerogeophysical survey profile around 88° S (Fig. 6A and Fig. A1a). The surface depression coincides with a 1000 m deep and 25 km wide bedrock low perpendicular to the profile ( (Studinger et al., 2006), and Fig. 3a). The 380 dominant wind direction is near perpendicular to the survey profile and follows the trend of the ice surface low (Fig. 3a). The high accumulation area also shows high surface roughness combined with steep slopes (Fig. 4 and Fig. A1). (2014)Another area with relatively high snow accumulation rates is located between 45° W and 0° longitude and is also exposed to precipitation from the Weddell Sea sector (e.g., Casey et al., 2014).
The snow accumulation rate at 88° S is spatially highly variable over very short length scales of several kilometers (Fig. 6c). Based 385 on visual inspection small-scale variability in snow accumulation rate correlates with small-scale variability in ice surface elevation (Fig. 6a), suggesting that wind-driven erosion and deposition is a primary process of snow accumulation. The relatively constant surface elevation between 60° E and 150° E shows very little variation in snow accumulation. In contrast, the short scale (< 10 km) undulations in ice surface elevation between 75° W and 45° E correspond to a highly variable (0 -20 cm w.e. yr -1 ) snow accumulation pattern with similar length scale (Fig. 6 a, c and Fig. A2). The dominant wind direction in the western segment is 390 subparallel to the profile (Fig. A2). Here, several pronounced peaks in snow accumulation rate correspond to topographic depressions in ice surface elevation (grey dashed lines in Fig. A2) indicating windblown deposition of snow. The eastern part of the profile has wind direction oblique or perpendicular to the profile. However, several peaks in snow accumulation rate still correlate with topographic depressions. Near 90° E the wind direction is parallel to the profile (Fig. 3). A pronounced peak in snow accumulation rate at 90° E correlates with 20 m deep depression in surface topography that is several kilometers wide (Fig. A3). 395 Accumulation decreases on the lee side of the topographic high at the western shoulder of the depression and increases towards the lowest part of the depression where it reaches its highest point (Fig. A3). The general correlation of highs in accumulation rate with lows in topography results in an overall negative correlation coefficient of r 2 = -0.33 between accumulation rate and MSWD ( Fig. 7a). DMS natural color imagery and laser altimetry data shows typically two dominant wind directions (Fig. 2). Our MERRA-2 wind direction is an average over seasonal variations in wind speed and likely reflects a wind direction somewhere between the 400 two dominant orientations of sastrugi. Since we have no knowledge of when a particular layer of snow has been deposited during a year it is not possible to do a more detailed analysis. The upper subset consists of data points above µ + σ and is marked by red circles and red r 2 values. The remaining data points fall within 410 ±σ from µ with r 2 shown in orange (data points are not marked for clarity).
The surface roughness derived from ATM laser altimetry reflects roughness on a scale of <250 m and does not reflect ice surface slope changes on length scales of several km. However, the MSWD (Fig. 6d) shows the same pattern of high variability between 75° W and 45° E and fairly constant values between 60° E and 150° E. To quantify the relationship between variability in surface slope, wind direction and accumulation rates we use the standard deviation σ of the MSWD and snow accumulation rate calculated 415 over a 20 km long moving window (Fig. 6e). Figure 6e shows the standard deviation of the MSWD and accumulation rate. In general, higher σ in accumulation rate generally occurs in areas with higher σ in MSWD. The match is strongest near 90° E where wind orientation is parallel to the profile.
The correlation coefficient between the standard deviations of the accumulation rate and the MSWD is r 2 = 0.50 indicating a stronger link between these variables than the actual parameters (Fig 7b). The magnitude of the correlation coefficient, however, 420 is dependent on the length scale used to calculate the standard deviation. Dattler et al. (2019) find a similar behavior between σ accumulation rate and σ MSWD. Visual inspection of Fig. 7 suggests that the relationship between accumulation rate and MSWD and the σ of accumulation rate and σ of MSWD is more pronounced for larger magnitudes of the variables. A kernel density estimate quantifies the probability density estimate of nearby points and allows visualization of the point density using a color scale for the data points (Fig. 7). We divide the data set into 3 subsets using the mean µ and σ: the lower end is defined by values 425 below µ -σ, while the upper end are values above µ + σ. The remaining points that are within ±σ from µ form the center subset.
We have calculated correlation coefficients r 2 for all subsets. In general, the correlation is strongest for the upper subsets, while the lower subsets show weak correlation. This is different from Dattler's et al. results (2019) who finds that the lower end also shows strong correlation. The upper tenth percentile of our data has an r 2 of 0.85 similar to Dattler's et al. results (2019). Most of Dattler's et al. data (2019) is located over high accumulation areas in West Antarctica. A possible explanation for the weak 430 correlation could be the very low accumulation rates in our area, combined with very small slopes and low wind speed. Noise in the elevation data will have a stronger impact on MSWD calculation and similarly, subtle changes in surface slope are likely below the resolution of MERRA-2, therefore resulting in a noisier and thus uncorrelated lower subset.

Relationship between surface roughness, slope and wind direction -results
Wind-related deposition and ablation processes could cause spatial roughness variations depending on surface slope and wind 435 direction. For example, windblown deposition of snow into concave surface depressions and ablation on up-slope areas could create spatial surface roughness patterns that correlate with slope and wind direction. We use the MSWD to determine if slopes that are exposed to uphill winds have different surface roughness than slopes experiencing primarily downhill winds (Fig. 8a). We calculate r 2 for up-slope winds (MSWD < 0) and down-slope winds (MSWD > 0). Neither the up-slope winds (r 2 = -0.14) nor down-slope winds (r 2 = 0.29) show any statistically significant correlation between surface roughness and slope as can be seen in 440 the scatter plot (Fig. 8a). Similarly, the surface roughness does not seem to be correlated with snow accumulation rates (r 2 = 0.28) indicating that there is also no statistically significant slope-independent relationship between surface roughness and accumulation rates within our survey area (Fig. 8b). Correlations may exist in smaller local areas, but our data shows that there is no consistent relationship between surface roughness, slope and wind direction on a regional scale within our survey area. However, our analysis is constrained by using 2-dimensional high-resolution roughness estimates and correlating it with 3-dimensional wind fields and 445 surface slope with much lower spatial resolution.

Radar backscatter and surface roughness
Surface roughness impacts the return signal of radar altimeters and can therefore cause elevation biases (e.g.,van der Veen et al., 2009, and references herein) similar to slope-dependent errors in altimetry data (Helm et al., 2014;Slater et al., 2018). Radar 455 backscatter in radar altimeters such as CryoSat-2 is a function of surface roughness. Surface roughness at the length-scales of the radar wavelength (2.2 cm) predominantly contributes to radar backscatter which cannot be resolved by our laser data. Changes in surface roughness cause changes in radar backscatter through changes of the echo waveform which can introduce range biases in the retrieval of surface elevation (Arthern et al., 2001;Kurtz et al., 2014). Figure 9 shows the maximum of the Ku-band radar return signal strength over the distinct roughness features identified from laser altimetry data (see Section 3.5). The maximum 460 return energy over rougher surface areas is about 3 dB lower than over the smooth areas in between. The difference in return signal strength is even more pronounced for the snow radar (4 dB, not shown). Stacked Ku-band waveforms shown in Fig. A4 show 3 dB higher surface return power at 3.053 µs over a smooth surface compared to the rough surface. The amplitude of the subsurface backscatter below 3.06 µs, however, is similar in strength over smooth and rough areas (Fig. A4). This observation is consistent with Gow's (1965) finding that heat from radiation causes crystal growth on the flanks of sastrugi, resulting in loosely bonded 465 crystals that are prone to erosion by moderate winds (Gow, 1965). This differential sublimation-deflation driven redistribution of snow which flattens the surface topography at the end of the summer resulting in relatively flat subsurface stratigraphy compared to the surface topography. The difference in waveform shapes between smooth and rough surfaces suggests that radar altimeters are potentially prone to elevation errors when threshold or leading edge trackers are being used for range retrieval. The Ku-band radar's relatively wide bandwidth and small footprint size allows resolution of the surface and sub-surface layers and thus accurate 470 tracking of the surface elevation. However, the reduced bandwidth and significantly larger footprint size of CryoSat-2 LRM returns does not allow for the resolution of individual layers, but instead leads to a pronounced broadening of the return waveform when the total cumulative backscatter of the sub-surface layers is close to, or exceeds, the backscatter from the surface layer.
The relatively small-scale nature and temporal variability of these features would require the use of more sophisticated retrieval techniques to better account for differences caused by the lower relative surface backscatter of rough areas. The elevation biases 475 caused by temporal and spatial variability in surface roughness are in addition to elevation biases caused by wind-induced anisotropy in the firn that have been identified from cross-over analysis (Armitage et al., 2014).  Fig. A4. b) Maximum of relative return signal strength from 2014 snow radar data. Red line is a running mean calculated over 350 radar traces (~2 km), which is similar in size to the 1.65 km CryoSat-2 footprint in low resolution mode (LRM) over smooth surfaces (Scagliola, 2013). The strength of the surface return is around 3 dB weaker over rough areas compared to smooth areas.

Conclusions 485
We have mapped the spatial and temporal variability in surface roughness and snow accumulation rate on a regional scale along a 1400 km circle around 88° S. We find significant small-scale variability (< 10 km) in snow accumulation based on snow radar subsurface stratigraphy, indicating areas of strong wind redistribution are prevalent at 88° S. The observed small-scale variability in snow accumulation rates is not captured by existing reanalysis models such as MERRA-2, which have low spatial resolution.
Our analysis shows that there is no simple relationship between surface slope, wind direction and snow accumulation rates for the 490 entire survey area. Previous studies have primarily focused on smaller regions often showing good correlation between surface slope and accumulation rates and are often used to falsely extrapolate parameters and relationships to larger regions beyond the study area. While we also observe these local correlations between surface slope, wind direction and accumulation rates, our results show that even for a homogenous area like the East Antarctic plateau near the South Pole such simple relationships don't exist on a regional scale. At the same time, we note that our accumulation rate measurements are a simple 2-dimensional view; until we 495 have 3-dimensional mapping of accumulation rates, these relationships might remain elusive. Our results underline the importance of regional-scale studies to derive accurate regional-scale parameterizations and relationships in light of expanding data sets, advances in high-performance computing and sophistication in model development. Similarly, we find high variability in surface roughness derived from laser altimetry measurements on length-scales smaller than 10 km, sometimes with very distinct and sharp transitions. These areas also show significant temporal variability over the course of the 3 survey years. We also find that surface 500 roughness does not seem to be correlated with snow accumulation rates. There seems to be no statistically significant slopeindependent relationship between surface roughness and accumulation rates within our survey area. The observed small-scale temporal and spatial variability in surface roughness will make it difficult to develop elevation bias corrections for radar altimeter retrieval algorithms.
Author contributions. MS led the analysis of the laser altimetry, optical imagery, integration of results and prepared the manuscript with contributions from all co-authors. BM derived surface roughness and accumulation rates from snow radar data and MERRA-2 and wrote the corresponding manuscript sections. KB, KC, and TN contributed to the analysis and interpretation of surface 515 roughness and accumulation rates. NK and TO contributed to the interpretation of radar backscatter and surface roughness. SM contributed to the analysis and interpretation of the ATM laser altimetry data. All authors helped write the paper. Brunt, K. M., Neumann, T. A., and Larsen, C. F.: Assessment of altimetry using ground-based GPS data from the 88S Traverse, Antarctica, in support of ICESat-2, The Cryosphere, 13, 579-590, 10.5194/tc-13-579-2019, 2019a. Brunt, K. M., Neumann, T. A., and.: Assessment of ICESat-2 Ice Sheet Surface Heights, Based on Comparisons Over the Interior of the Antarctic Ice Sheet, Geophysical Research Letters, 46, 13072-13078, 10.1029/2019gl084886, 2019b.