Comparative analysis of snowfall accumulation over Antarctica in light of Ice discharge and gravity observations from space

The remote and cold Antarctic continent presents unique challenges to quantify precipitation rates from space and in situ observations. This has resulted in large uncertainties in current estimates. In this study, we quantify annual precipitation rates over seven Antarctic basins using a novel mass budget (MB) approach, by building on the recent Landsat based estimate of ice discharge and changes in total water storage from GRACE. The MB precipitation rates are compared with those from CloudSat, GPCP, the Arthern precipitation climatology, the GPM constellation sensors, a few popular reanalysis products, and a regional climate model for two periods: 2007–2010 and 2013–2015. The new estimates are bounded by CloudSat precipitation rates with and without adjustment for the unmeasured near surface precipitation. GPM products significantly underestimate Antarctic precipitation rate, but capture spatial variability that is valuable for bias-adjustment. We find variable performance between products at basin scale, suggesting that an in-depth regional study of precipitation rates is necessary.


Introduction
Accurate estimation of precipitation rates over Antarctica is key to refining our understanding of past and present changes in the mass of the Antarctic ice sheet, sea level change, and the Earth's water and energy budget. In a warming climate, Antarctic precipitation, which is currently almost entirely snowfall, is expected to increase (Frieler et al 2015;Palerme et al 2017). Despite the near consensus in the expected sign of the trend in precipitation, the rate of increase is debated and climate models differ with measurements (Genthon and Krinner 2001) and between models for present-day precipitation rates over Antarctica (Behrangi and Richardson 2018). In addition to changes in precipitation rates, current observational evidence shows an acceleration in ice loss for the West Antarctic Ice Sheet with enhanced losses from those glaciers that flow into the Amundsen Sea Embayment (Rignot et al 2008, Velicogna 2009, Shepherd et al 2012. It is therefore expected that accelerated rates of ice loss resulting from enhanced glacier flow will be offset by increases in snowfall, somewhat moderating the future contribution of the Antarctic ice sheet to sea level rise (Seroussi et al 2020).
Due to the harsh climate and little logistical infrastructure, in situ measurements of precipitation in Antarctica are sparse and the global precipitation products that rely on these observations (i.e. the Global Precipitation Climatology Centre (GPCC)) (Schneider et al 2017) have no coverage over this region. As a result, other observational products, such as the Global Precipitation Climatology Project (GPCP) (Huffman et al 2001;Huffman et al 2009;Adler et al 2017), that uses GPCC for bias correction over land, has relied on satellite-only precipitation estimates. Note that even those available in situ data often have large uncertainties. High wind speeds, typical of the Antarctic, result in significant gauge undercatch and current correction factors seem to be highly uncertain . Satellite products also face large uncertainties over cold regions such as Antarctica due to insufficient sensitivity of sensors to detect and estimate precipitation signals, complex surface emissivities (Ferraro et al 2013), and poor understanding of precipitation microphysics, among others. Nonetheless, Antarctic precipitation has been estimated from various sources. Building on some previous analysis (e.g. Vaughan et al 1999, Glovinetto andZwally 2000, Arthern et al 2006 produced a climatology map of Antarctic snowfall accumulation from in situ glaciological measurements covering the period 1950-2000. They used Satellite observations from the Advanced Microwave Scanning Radiometer (AMSR-E) and snow temperatures derived from satellite observations of thermal infrared emission to guide their interpolation method.
More recently, CloudSat (Stephens and Ellis 2008) was used to estimate Antarctic snowfall accumulation (Palerme et al 2014, Behrangi et al 2016 given its unprecedented sensitivity (~− 28 dBZ) to detect light precipitation and snowfall. These estimates compared reasonably with GPCP, reanalyses, and the Arthern et al climatology Palerme et al 2014), but concerns remained about uncertainties in ice crystal habit, ice particle size distribution, and missing near-surface precipitation. CloudSat precipitation is estimated from the sixth range bin (~1300 m) above the surface to avoid contamination of the reflectivity profile by surface returns (Tanelli et al 2008). (Grazioli et al 2017) used a combination of new in situ data collected over one year on the coast of Adelie Land and showed that dry air in lower troposphere can sublimate about 17% of falling snow before it reaches the surface. Knowing that the CloudSat snowfall estimate is made at 1.00-1.50 km above the surface, CloudSat may not see the reduced snowfall rate due to the sublimation. This may offset some missing precipitation by CloudSat or lead to over estimation by CloudSat. Furthermore, (Palerme et al 2019) concluded that ground clutter can cause CloudSat to significantly overestimate snowfall rates in parts of Antarctica. In 2011 Cloud-Sat faced battery issues that restricted its sampling to daytime only and in 2018 it exited the A-Train.
Gravity Recovery and Climate Experiment (GRACE) satellites, launched in 2002 with science measurements extending through 2017, measured small changes in Earth's gravity caused by nearsurface mass movement (i.e. water, ice, and solid earth). The value of the GRACE observations to estimate snowfall accumulation has been shown by a few studies over Eurasian pan-Arctic (Swenson 2010, Landerer et al 2010, Tibetan Plateau (Behrangi et al 2017), and cold regions in the Northern Hemisphere . GRACE estimates of snowfall accumulation are based on measuring the total water storage anomaly (TWSA) and using the mass conservation principle (Behrangi et al 2017). This approach enables an independent estimate of precipitation accumulation (i.e. gravimetry versus radiometry and in situ observations) and is not affected by several issues affecting other space-based estimates, several of which are discussed above. GRACE has also been used to study the magnitude of recent extreme precipitation events along the East Antarctic coast (Boening et al 2012). However, no comprehensive GRACE-based precipitation analysis has been performed over the entire Antarctic continent. This is partly because estimating snowfall accumulation over Antarctic ice sheets is complicated by ice divergence, the continual export of mass from the interior to the oceans via ice flow .
In a recent study,  reconstructed ice discharge and changes in ice discharge over the Antarctic ice sheet by combining a comprehensive record of recent changes in Antarctic-wide ice flow, calculated from novel feature tracking methods to hundreds of thousands of Landsat image, with optimized flux gate definitions and an earlier mapping of surface velocity . They calculated ice discharge values and their uncertainties for two years (2008 and 2015) and for 27 basins defined by Zwally et al (2002), covering almost the entire Antarctic ice sheet. These new estimates complement improved GRACE-based mass change observation (Watkins et al 2015) and provide an unprecedented opportunity to estimate basin-scale precipitation accumulation using the water budget approach.
In this study we build on this unique opportunity and calculate basin-wide precipitation rate estimates over the Antarctic ice sheet. The new estimates are then compared with existing estimates from Cloud-Sat, the precipitation climatology map of (Arthern et al 2006), the Global Precipitation Measurement (GPM) mission (Skofronick-Jackson et al 2017) constellation sensors, GPCP V2.3, a few popular reanalysis products, and a regional climate model.

Method and data
Precipitation accumulation is calculated by solving the mass conservation equation (e.g. Dingman 2008) for individual Antarctic basins. If t 1 and t 2 represent the start and end of the accumulation time for basin B, precipitation accumulation (P) between time t 1 and t 2 can be calculated as: Where S is the change in storage (from GRACE observations) between time t 1 and t 2 and Q net (t), ET(t), and Sub(t) are net lateral flux, evapotranspiration, and sublimation rates for basin B at time t, respectively. Most of Antarctica experiences below freezing temperature throughout the year, and therefore runoff from surface melting is negligible (Van den Broeke et al 2011, Kuipers Munneke et al 2012, hence Q net is essentially the ice discharge plus basal melt (melting at the bed of the ice sheet due to pressure, geothermal heat flux and friction). ET is also negligible, but sublimation needs to be considered.  figure 1(a)) to 7 larger basins to reduce leakage errors across drainage basin boundaries when analyzing the GRACE data. Ice flow (governs ice discharge) shown with purple lines  Our analysis focuses on two distinct periods: 2007-2010 and 2013-2015. These were chosen to be as closely coincident with available annual ice discharge data (2008 and 2015), and CloudSat data (2007)(2008)(2009)(2010). Ice discharge changes little from year to year in Antarctica , and can be approximated well with a linear trend between the two periods; however, using timespans near the epoch of the ice discharge data availability reduces any interpolation errors. It is advantageous to use multiyear time spans (rather than monthly or annual) to reduce random errors in both GRACE and Cloud-Sat data, as well as other precipitation data products examined.
Note that precipitation over Antarctica, and especially over the basins studied here, is entirely snowfall (Palerme et al 2014). Therefore, hereafter, the terms precipitation and snowfall are used interchangeably.

Mass change estimate from GRACE
Mass change estimates were computed using the JPL RL06M Version 2 GRACE and GRACE-FO mascon solution (Watkins et al 2015, Wiese et al 2019, which parameterizes the Earth's gravity field in terms of equal-area 3 • spherical cap mass concentration elements. The solutions have monthly temporal resolution and a spatial resolution of approximately 333 km. Previous analysis has shown each individual 3 • mascon element to be uncorrelated with neighboring elements at high latitudes (Schlegel et al 2016). To reduce leakage errors across drainage basin boundaries, when analyzing the GRACE data, we merged the 27 Antarctic basins into 7 major basins (figure 1(a)). Land-ocean leakage errors were minimized by applying the Coastline Resolution Improvement (CRI) filter  to the data. Leakage errors between the seven regions used in this study were further corrected by applying a set of scaling factors to the GRACE data to effectively downscale/redistribute mass within each individual mascon element according to the spatial pattern in the altimetry-derived trend map from CryoSat-2 , consistent with the methods described in  for downscaling hydrology-related signals. Additional standard corrections to the GRACE data were applied in this analysis, including replacing the C 20 coefficient with that obtained from satellite laser ranging using TN-14 (Loomis et al 2019), using a geocenter estimate computed using the methods of (Sun et al 2016), and applying a correction for glacial isostatic adjustment (GIA) processes using the regional model of (Ivins et al 2013). The trend estimates presented here were obtained via a simultaneous fit of an annual, semiannual, trend, and a tidal alias signal from the S 2 constituent to the GRACE data. Uncertainties on the estimates of mass change were computed as the root sum of squares of the GRACE measurement error (reported with the JPL RL06M Version 2 mascon solution), leakage correction error, and GIA model error. The leakage correction error is conservatively assigned to be 50% of the magnitude of the full leakage correction (including that applied by the CRI filter). GIA model error is taken to be the 1-sigma spread between four competing GIA models (Whitehouse et al 2012, Ivins et al 2013Peltier et al 2018;Purcell et al 2016): two regional models and two global models.

Ice discharge
Ice discharge (flow of solid ice into the ocean) is calculated for two time periods in 2008 and 2015 reported in . Ice flow is determined from the radar mapping of (Rignot et al 2011) for 2008 and from featuring tracking of Landsat imagery for the 2015. The ice flux is then calculated as the flow of ice through the cross-section of a flux gate that is optimized to follow flight paths of airborne campaigns equipped with radar echo sounders that measure ice thickness and that were flown in close proximity to the ice sheet grounding line (where ice transitions from grounded to floating). The fluxgate was modified for temporal changes in ice thickness using CryoSat-2 satellite altimetry data processed according to (Nilsson et al 2016). This approach assumes that the surface velocity is equal to depth-averaged velocity (i.e. ice motion is entirely due to sliding at its base). This assumption introduces little uncertainty into the analysis . Total discharge estimates also take into account the additional mass flux due to the accumulation of snowfall (net combination of sublimation and deposition as well as snow drift) that occurs between the defined flux gate and the grounding line. This additional mass flux was small and was accounted for using output from a regional climate model (Van Wessem et al 2014;Van Wessem et al 2017). In addition to ice flux, ice loss also occurs from melting at the base of the ice sheet that then finds its way into the ocean. This melting is due to the decreased melting point of ice with pressure (ice thickness), geothermal heat flux and friction at the bed and the ice sheet. We account for this using the estimate of (Van Liefferinge and Pattyn 2013). Ice fluxes were calculated along a flux gate with an average spacing between nodes of~200 m and surface velocities, from which fluxes were derived, were resolved at approximately annual resolution.

Sublimation rate
Sublimation rate is obtained from the new version of the Regional Atmospheric Climate Model V2.3p2 (RACMO2) (van Wessem et al 2018), providing monthly high-resolution (0.25 • by 0.25 • ) simulation of sublimation, precipitation climatology, and other meteorological variables over the ice sheets of Greenland and Antarctica. RACMO2, combines the dynamical core of the High Resolution Limited Area Model (HIRLAM) numerical weather prediction model with ECMWF Integrated Forecast System (ISF) physics. As a regional model, RACMO needs external information at the lateral boundaries and sea surface that are obtained from ERA-I. Temperature, specific humidity, zonal and meridional wind components, surface pressure, sea surface temperature, and sea ice concentration are relaxed towards the fields of a global model at the lateral boundary zone of the model. The new version, RACMO2.3p2, incorporates updated glacier outlines, topography and ice albedo fields.
RACMO2 provides a net combination of sublimation and deposition as well as snow drift. Sublimation over the studied basin is about 7% of total precipitation out of which surface sublimation accounts for about 2.5% and drifting snow sublimation accounts for about 4.3% of total precipitation. These numbers are consistent with those reported by (Van Wessem et al 2018). Note that drifting snow is limited to below 2 m above the surface, and blowing snow particles are above that level . However, the contribution of blowing snow processes compared to the total precipitation is negligible (Wang et al 2016).

Other precipitation data sets
Several other precipitation data sets are used for comparison including CloudSat snowfall, Arthern climatology map, GPCP, several reanalysis, a regional climate model, and GPM products. CloudSat 2C-SNOW-PROFILEP1-R05 (Wood et al 2014) was obtained from http://cloudsat.atmos.colostate.edu. CloudSat snowfall data comes in orbits, retrieved from Cloud Profiling Radar, and is reported for each CPR footprint (~1.7 km along track and~1.4 km cross track). CloudSat revisit time over Antarctica is zonal dependent and can vary between two to five days (Souverijns et al 2018). Therefore, CloudSat product may not be great to study individual snowfall events over Antarctica that often span several hours. However, when averaged over large areas (e.g. 1 • x 2 • grid box) it has been shown capable to capture stable statistics for annual comparison with other products including reanalysis (Palerme et al 2014;Souverijns et al 2018). Given that our studied basins cover very large areas, we do not expect issues due to the temporal sampling of CloudSat in capturing mean annual snowfall rates. Arthern climatology map, compiled using data for years 1950-2000(Arthern et al 2006, has an effective resolution of approximately 100 km by 100 km. The map was acquired from the British Antarctic Survey (https://legacy.bas.acuk/). While the Arthern climatology map does not overlap with our study period, we included it in our comparison as it is a popular product that is sometimes used as an observational reference for comparison with more recent data sets precipitation products at 0.25 • by 0.25 • resolution was obtained, for all of the GPM constellation sensors, from the Goddard Earth Sciences Data and Information Services Center (GES DISC). Note that in high latitudes, including Antarctica, these products still face large uncertainties as discussed in the introduction.  1(b)) are derived from equation (1) using ice discharge plus basal melt, GRACE total water storage anomaly, and sublimation rates (see table 1 for individual components and uncertainties).

Results
The reported uncertainties are based on uncertainties calculated for ice discharge and basal melt (Gardner et al 2018) and GRACE estimate (section 2.1). Actual uncertainties are expected to be larger, partly because uncertainties in sublimation rate and snowdrift are not considered. Two estimates are provided for CloudSat, the original estimate prior to adjustment ( figure 1(c)) and the one after multiplying an adjustment factor to account for the missed precipitation or changes in precipitation profile in the lowest 1.2 km profile of CloudSat due to the ground clutter issue ( figure 1(d)). Similar to (Grazioli et al 2017), the adjustment factor is calculated using the ECMWF Integrated Forecast System (ECMWF IFS) and by dividing the cumulative precipitation at near surface by precipitation accumulation at 1.2 km above surface. Furthermore, we used ERA-I pattern, to extend CloudSat estimates from 81 • S to 90 • S, as CloudSat coverage is within 81 • S/N. GPCP V2.3 (figure 1(e)) estimate over Antarctica is mainly based on precipitation rate estimated from Infrared sounders (Adler et al 2003), in particular the Atmospheric Infrared Sounder (AIRS) data after 2003. The Arthern precipition map (figure 1(f)) is based on compilation of in situ glaciological measurements for 1950-2000, thus may not necessarily represent mean rate over the 2007 to 2010 period. The rest of the panels (figures 1(g)-(l)) display precipitation rates from different reanalysis products, and RACMO2. Figure 1 shows that P MB ( figure 1(b)) falls between the snowfall rates estimated by CloudSat ( figure 1(c)) and CloudSat-adj. (figure 1(d)). P MB is consistent with all of the products, except GPCP, in assigning the highest and lowest precipitation rates to Basins 1 and 2, respectively. GPCP V2.3 (figure 1(e)) overestimates in Basins 2 and 7, compared to P MB and most of the studied products. Arthern climatology map is also different from other products,  including P MB , as it assigns much larger precipitation rate to Basin 2 and much smaller precipitation rate to Basin 6. A quantitative comparison of the precipitation rates at each basin is provided in table 2, with bold fonts representing estimates that fall inside the uncertainty range of P MB . For the 2007-2010 period, P MB is consistent within its uncertainty range with those from RACMO2 in Basin 1, CloudSat-adj in Basin 2, NCEP/DOE in Basin 3, RACMO2, ERA-5, MERRA, and MERRA-2 in Basin 4, GPCP, Arthern, MERRA, and NCEP/DOE in Basin 5, NCEP/DOE in Basin 6, and RACMO2 and NCEP/DOE in Basin 7. For all basins combined, the mean P MB is most similar to RACMO2 and NCEP/DOE, although this agreement can be due to cancelation of over and under estimated precipitation rates across different basins. This suggests that besides comparison of annual precipitation rates (e.g. Behrangi et al 2016, Palerme et al 2017, regional analysis of precipitation intensity is necessary to advance precipitation analysis over Antarctica. P MB does not favor any specific product based on regional analysis, but RACMO2 and NCEP/DOE agree with P MB in 3 and 4 (out of 7) basins within the P MB uncertainty range.
Repeating the analysis for 2013-2015 (figure 2 and table 1) suggests that the outcomes are fairly consistent in the two periods. Precipitation rate from RACMO2 agrees with P MB in 4 (out of 7) basins followed by NCEP/DOE and ERA-5 that agree with P MB in 3 basins within the uncertainty range of P MB . It should be noted that since April of 2011, CloudSat operates only during daytime and misses most of the Antarctic ice sheet during austral winter. Therefore, we could not use CloudSat to produce annual precipitation rates over Antarctica for 2013-2015. Instead, in figure 2 and table 2 we added mean annual precipitation rates from four space born sensors (i.e. SSNIS on DMSP F16, MHS on NOAA18, ATMS on SNPP, and AMSR2 on GCOM-W1), each represent a type of sensor flying as part of the GPM constellation sensors. The GPM Core Observatory (i.e. GMI and DPR) are not included in the analysis as their coverage is limited to ± 67 • . P MB and other products suggest that GPM constellation sensors significantly underestimate the annual mean rate across all seven basins (figure 2 and table 2). The GPM products, however, capture precipitation patterns fairly well and show high correlations with P MB estimates (figure 3). The large systematic errors might partly relate to insufficient sensitivity of sensors to capture snowfall signals, complex surface emissivities, generally cold and dry atmosphere over the region, and poor understanding of precipitation microphysics and particle size distributions. However, if precipitation spatial patterns are captured, systematic errors can be reduced. The red Asterisks in figure 3 display the GPM and GPCP estimates after applying adjusting ratios that are calculated by dividing P MB values by their corresponding estimates from GPM and GPCP during 2007-2010. The results suggest that such bias correction scheme can be effective, at least at basin scale.

Concluding remarks
In this study we calculated annual precipitation (i.e. snowfall) over seven Antarctic basins via a mass budget (MB) approach. This is the first study to take this approach to assessing precipitation rates in Antarctic and is now possible thanks to recent estimates of ice discharge rates (and their uncertainties), changes in total water storage anomaly observed by GRACE, and estimated sublimation from RACMO2. The ice discharge rates are obtained using novel feature tracking methods applied to hundreds of thousands of Landsat images and through optimized flux gate definitions . P MB estimates are then compared with other estimates from Cloud-Sat, the precipitation climatology map of (Arthern et al 2006), the GPM constellation sensors, GPCP V2.3, ERA-Interim, ERA-5, MERRA, MERRA-2, NCEP-DOE, and RACMO2 for 2007-2010 and 2013-2015 periods.
The results show that precipitation estimates from MB fall between CloudSat estimates with and without adjustment. The P MB 's rankings at individual basins show variable performance between products at basin scale, suggesting that regional analysis of precipitation intensity is necessary to advance precipitation analysis over Antarctica. For all basins combined, the mean P MB is most similar to that of RACMO2, MERRA-2, and NCEP/DOE for both 2007-2010 and 2013-2015 periods. During 2007-2010 mean annual precipitation estimates from CloudSat (113.2 mm y −1 ) and ERA-Interim (115.8 mm y −1 ) averaged over all 7 basins were similar, but about 30 mm y −1 less than P MB . Given the total area of the seven basins (11 894 700 km 2 ) and assuming that 360 Gt of ice is equivalent to 1 mm in sea level (Alley et al 2005), this difference is equivalent to~1 mm yr −1 change in global sea level contribution from the Antarctic Ice Sheet when employing the glacier mass budget approach . Similarly, GPCP shows an overestimation of about 18 mm yr −1 compared to P MB , equivalent to~0.6 mm yr −1 of sea level.
Based on P MB , GPM products significantly underestimate annual snowfall accumulation across all seven basins, yet present high spatial correlations with P MB estimates. We showed that by applying adjusting ratios, calculated by dividing P MB values by their corresponding estimates from GPM and GPCP for 2007-2010, systematic errors of GPM products in 2013-2015 can be reduced significantly.
Our analysis shows that CloudSat has better agreement with P MB (and most other products) than GPM sensors to estimate basin-scale annual precipitation rate. This is partly because of the high sensitivity of CloudSat radar to detect snow particles, mitigating some of the issues that GPM sensors face. CloudSat left the A-train in early 2018, but with the recent launch of GRACE-FO and continued production of ice discharge from Landsat  we can produce basin-scale precipitation accumulation using the MB approach for years ahead. The launch of EarthCARE (Illingworth et al 2015), currently planned for 2021, will also restore some of the great capabilities offered by CloudSat. The combination of recent advances in retrieving precipitation over polar regions and new techniques to constrain the estimates creates a great opportunity to further improve quantification of precipitation and advance science over such remote but important regions.