Estimation of the Antarctic surface mass balance using the regional climate model MAR (1979–2015) and identification of dominant processes

The Antarctic ice sheet mass balance is a major component of the sea level budget and results from the difference of two fluxes of a similar magnitude: ice flow discharging in the ocean and net snow accumulation on the ice sheet surface, i.e. the surface mass balance (SMB). Separately modelling ice dynamics and SMB is the only way to project future trends. In addition, mass balance studies frequently use regional climate models (RCMs) outputs as an alternative to observed fields because SMB observations are particularly scarce on the ice sheet. Here we evaluate new simulations of the polar RCM MAR forced by three reanalyses, ERA-Interim, JRA-55, and MERRA-2, for the period 1979–2015, and we compare MAR results to the last outputs of the RCM RACMO2 forced by ERA-Interim. We show that MAR and RACMO2 perform similarly well in simulating coast-to-plateau SMB gradients, and we find no significant differences in their simulated SMB when integrated over the ice sheet or its major basins. More importantly, we outline and quantify missing or underestimated processes in both RCMs. Along stake transects, we show that both models accumulate too much snow on crests, and not enough snow in valleys, as a result of drifting snow transport fluxes not included in MAR and probably underestimated in RACMO2 by a factor of 3. Our results tend to confirm that drifting snow transport and sublimation fluxes are much larger than previous model-based estimates and need to be better resolved and constrained in climate models. Sublimation of precipitating particles in low-level atmospheric layers is responsible for the significantly lower snowfall rates in MAR than in RACMO2 in katabatic channels at the ice sheet margins. Atmospheric sublimation in MAR represents 363 Gt yr⁻¹ over the grounded ice sheet for the year 2015, which is 16 % of the simulated snowfall loaded at the ground. This estimate is consistent with a recent study based on precipitation radar observations and is more than twice as much as simulated in RACMO2 because of different time residence of precipitating particles in the atmosphere. The remaining spatial differences in snowfall between MAR and RACMO2 are attributed to differences in advection of precipitation with snowfall particles being likely advected too far inland in MAR.

1 Introduction

Mass loss from the Antarctic ice sheet (AIS) and therewith its contribution to the sea level budget results from the difference of two fluxes of a similar magnitude: ice flow discharging in the ocean (D) and net snow accumulation on the ice sheet surface, i.e. the surface mass balance (SMB). The total ice sheet mass balance (SMB minus D) can be assessed using satellite altimetry, gravimetry, or the input–output method (Shepherd et al., 2018), which all request SMB estimates. The input–output method, which consists in separately modelling ice dynamics and SMB, is also the only way to project future trends.

SMB as used in this study is the sum of mass gains (mainly snowfall accumulation and some riming), mass losses (mainly surface and drifting snow sublimation, some liquid water run-off), and drifting snow transport (defined as the horizontal advection of the drifting snow), which can lead to either mass gain or mass loss. Snowfall rates are 1 order of magnitude larger than all of the other SMB fluxes at the continental scale (Lenaerts et al., 2012b), with the largest amounts found along the ice sheet margins due to cyclonic activity in the Southern Ocean and to the orographic lifting of relatively warm and moist air masses (van Wessem et al., 2014; Favier et al., 2017). Accumulation patterns are highly variable at the kilometre scale and from year to year (e.g. Agosta et al., 2012). Consequently, proper observations of SMB require a high spatial coverage (e.g. stake lines, accumulation radars plus ice cores for layer dating and snow density) and a temporal sampling spanning several years (Eisen et al., 2008). Even if efforts have been made to fulfil those requirements, ground-based observations are scarce and carry with them high logistical costs in this cold, windy, and remote environment. Interpolation techniques used to interpolate the scarce SMB observations (Vaughan et al., 1999; Arthern et al., 2006) encounter major caveats (Magand et al., 2008; Genthon et al., 2009; Picard et al., 2009).

This is why many AIS mass balance studies use output of regional climate models (RCMs) to estimate ice sheet SMB for the recent decades (e.g. Rignot et al., 2011; Gardner et al., 2018; Shepherd et al., 2018). In order to obtain a good agreement with observations, atmospheric models require accurate large-scale circulation patterns together with a proper representation of snow surface processes, clouds, and turbulent fluxes and a relatively high horizontal resolution to properly resolve the complex ice sheet topography at the margins.

Here, we present new simulations of the RCM MAR, applied for the first time over the whole AIS, but already widely used for polar studies, e.g. in Greenland (Fettweis et al., 2013, 2017), Svalbard (Lang et al., 2015), Adélie Land (Antarctic coastal area; Gallée et al., 2013; Amory et al., 2015), and Dome C (Antarctic plateau; Gallée et al., 2015). We compare MAR-simulated SMB with the state-of-the-art RCM RACMO2 (van Wessem et al., 2018). We use available SMB observational datasets to show that MAR and RACMO2 perform similarly well in simulating the SMB spatial gradients. In addition, we identify significant processes that still need to be included or improved in both RCMs.

In Sect. 2, we describe MAR and its specific set-up for Antarctica, together with RACMO2, the forcing fields, observational datasets, and methods designed for model evaluation. In Sect. 3, we show that both RCMs share common biases against observed SMB, resulting from drifting snow transport fluxes. Secondly, we analyse SMB differences among models and show that many of the discrepancies can be attributed to low-level sublimation of precipitation in katabatic channels and to the difference in precipitation advection inland. Finally, in Sect. 4, we summarise our main findings and discuss further efforts to be achieved for a better assessment of the AIS SMB.

2 Data and methods

2.1 Regional modelling

2.1.1 Regional atmospheric models

For the first time, the polar-oriented regional atmospheric model MAR is applied for decades-long simulations over the whole AIS. MAR atmospheric dynamics are based on the hydrostatic approximation of the primitive equations, fully described in Gallée and Schayes (1994). Prognostic equations are used to depict five water species: specific humidity, cloud droplets and ice crystals, raindrops, and snow particles (Gallée, 1995). Sublimation of airborne snow particles is a direct contribution to the heat and moisture budget of the atmospheric layer in which these particles are simulated. The radiative transfer through the atmosphere is parametrised as in Morcrette (2002), with snow particles affecting the atmospheric optical depth (Gallée and Gorodetskaya, 2010). The atmospheric component is coupled to the surface scheme SISVAT (soil ice snow vegetation atmosphere transfer; De Ridder and Gallée, 1998) dealing with the energy and mass exchanges among surface, snow, and atmosphere. The snow-ice part of SISVAT is based on the snow model CROCUS (Brun et al., 1992). It is a one-dimensional multilayered energy balance model which simulates meltwater refreezing, snow metamorphism, and snow surface albedo depending on snow properties. We used MAR version 3.6.4, simply called MAR hereafter. In this version the physical settings are the same as in MAR version 3.5.2 used for Greenland (Fettweis et al., 2017), except for the adaptations detailed below.

Grid. Projection is the standard Antarctic polar stereographic method (EPSG:3031). The horizontal resolution is 35 km, an intermediate resolution that results from a computation time compromise in order to run the model with multiple reanalyses and global climate model forcings over the 20th and the 21st centuries. The vertical discretisation is composed of 23 hybrid levels from ∼2 m to ∼17 000 m above the ground.

Boundaries. The topography is derived from the Bedmap2 surface elevation dataset (Fretwell et al., 2013). Because the Antarctic domain is about 4 times larger than the Greenland domain, the circulation has to be more strongly constrained. This is why we use a boundary relaxation of temperature and wind in the upper atmosphere starting from 400 hPa (∼6000 m above the ground) to 50 hPa (upper level), as in van de Berg and Medley (2016), whereas relaxation starts from 200 hPa in Fettweis et al. (2017).

Parameterisations.

• a.

  The surface snow density ρ_s (kg m⁻³) is computed as a function of 10 m wind speed ws₁₀ (m s⁻¹) and surface temperature T_s (K):

  $$\begin{array}{}\text{(1)}& {\mathit{\rho }}_{\mathrm{s}}=\mathrm{149.2}+\mathrm{6.84}{\mathrm{ws}}_{\mathrm{10}}+\mathrm{0.48}{T}_{\mathrm{s}},\end{array}$$

  with minimum–maximum values of 200–400 kg m⁻³. This parameterisation was defined so that the simulated density of the first 50 cm of snow fits observations collected over the AIS (see Fig. S1, with the snow density database detailed in Table S1 in the Supplement).

• b.

  The aerodynamic roughness length z₀ is computed as a function of the air temperature, as proposed in Amory et al. (2017). The parameterisation was tuned so that z₀ fit the observed seasonal variation between high (>1 mm) summer and lower (0.1 mm) winter values in coastal Adélie Land, for air temperatures above −20 ^∘C. For lower temperatures, z₀ is kept constant and set to 0.2 mm, in agreement with observed z₀ values on the Antarctic plateau (e.g. Vignon et al., 2016);

• c.

  As in Fettweis et al. (2017), the MAR drifting snow scheme is not activated because this scheme was sensitive to parameter choices (Amory et al., 2015). An updated version of the drifting snow scheme is currently being developed and evaluated for application at the scale of the whole ice sheet.

We compare MAR results over the AIS to the latest outputs of the regional atmospheric model RACMO2 version 2.3p2 (van Wessem et al., 2018), called RACMO2 hereafter, using a horizontal resolution of 27 km, a vertical resolution of 40 atmospheric levels, and a topography based on the digital elevation model from Bamber et al. (2009). This regional model is developed by the Royal Netherlands Meteorological Institute (KNMI) and has subsequently been adapted for modelling the Antarctic climate and its SMB (van de Berg et al., 2006). It includes a drifting snow scheme (Lenaerts et al., 2012a), an albedo routine with prognostic snow grain size (Kuipers Munneke et al., 2011), and a multilayer snow model computing melt, percolation, refreezing, and run-off (Ettema et al., 2010).

MAR and RACMO2 models were developed independently. We will not detail here the many physical parameterisation differences between both RCMs, but we will later highlight some of them we show having a significant impact on the modelled SMB.

2.1.2 Forcing reanalyses

Regional atmospheric models are forced by atmospheric fields at their lateral boundaries (pressure, wind, temperature, humidity), at the top of the troposphere (temperature, wind), as well as by sea surface conditions (sea ice concentration, sea surface temperature) every 6 h. Consequently, regional atmospheric models add details and physics to the forcing model in the middle and lower troposphere and at the land or iced surface, whereas large-scale circulation patterns are driven by the forcing fields. We forced MAR with three reanalyses over Antarctica in order to evaluate the uncertainty in the simulated surface climate arising from the uncertainty in the assimilation systems: the European Centre for Medium-Range Weather Forecasts “Interim” re-analysis (hereafter ERA-Interim, resolution ∼0.75^∘, i.e. ∼50 km at 70^∘ S; Dee et al., 2011), the Modern-Era Retrospective Analysis for Research and Applications version 2 (hereafter MERRA-2, resolution ∼0.5^∘; Gelaro et al., 2017), and the Japanese 55-year Reanalysis from the Japan Meteorological Agency (hereafter JRA-55, resolution ∼1.25^∘; Kobayashi et al., 2015).

The regional atmospheric model RACMO2 is forced by ERA-Interim. We focus our study to the period 1979–2015, as reanalyses are known to be unreliable before 1979, when satellite sounding data started to be assimilated (Bromwich et al., 2007).

2.2 Observations

2.2.1 SMB observations and sectors of strong SMB gradients

We use SMB observations of the GLACIOCLIM-SAMBA dataset detailed in Favier et al. (2013) and updated by Wang et al. (2016). This dataset is an update of the one assembled by Vaughan et al. (1999) following the quality-control methodology defined by Magand et al. (2007). It includes 3043 reliable SMB values averaged over more than 3 years. We add accumulation estimates from Medley et al. (2014), retrieved over the Amundsen Sea coast (Marie Byrd Land) with an airborne-radar method combined with ice-core glaciochemical analysis.

The first-order feature of the Antarctic SMB is a strong coastal–inland gradient, with mean values ranging from typically greater than 500 kg m⁻² yr⁻¹ at the ice sheet margins to about 30 kg m⁻² yr⁻¹ in the dry interior plateau (Fig. 1a; see also, e.g. Wang et al., 2016). As observations only cover 5 % of MAR grid cells over the ice sheet, we divide that sparse observation dataset into 10 sectors detailed in Table 1 and shown in Fig. 2. Six of them are stake transects with a stake every ∼1.5 km, which have been proven very valuable for evaluating modelled SMB (Agosta et al., 2012; Favier et al., 2013; Wang et al., 2016). The four other sectors are composed of more scattered observations covering large elevation ranges (Victoria Land, Dronning Maud Land, and Ross Ice Shelf–Marie Byrd Land).

Table 1 Sectors extracted from the GLACIOCLIM-SAMBA database.

¹ Medley et al. (2014). ² Clausen et al. (1979). ³ Venteris and Whillans (1998). ⁴ Vaughan et al. (1999). ⁵ Magand et al. (2007). ⁶ Frezzotti et al. (2004). ⁷ Frezzotti et al. (2007). ⁸ Pettré et al. (1986). ⁹ Agosta et al. (2012). ¹⁰ Verfaillie et al. (2012). ¹¹ Goodwin (1988). ¹² Ding et al. (2011). ¹³ Wang et al. (2016). ¹⁴ Higham and Craven (1997). ¹⁵ Wang et al. (2015). ¹⁶ GLACIOCLIM-BELARE. ¹⁷ Picciotto et al. (1968). ¹⁸ Mosley-Thompson et al. (1995). ¹⁹ Mosley-Thompson et al. (1999). ²⁰ Anschütz et al. (2011).

2.2.2 Model–observation comparison method

RACMO2 outputs are bilinearly interpolated to the 35×35 km MAR grid. For each SMB observation, we consider the four surrounding MAR grid cells, from which we eliminate ocean grid cells. We also eliminate surrounding grid cells with an elevation difference with the observation greater than 200 m (missing elevation of observation is set to Bedmap2 elevation at 1 km resolution). Finally, we bilinearly interpolate model values of the remaining grid cells at the observation location (see schematic in Fig. S2).

As we restrict our modelling study to the 1979–2015 period, we only consider observations beginning after 1950. For observations beginning after 1979, we time-average model outputs for the same period as the observation. We keep observations beginning before 1979 only if they cover more than 8 years, and in this case we compare the observed value with the modelled value time-averaged for 1979–2015.

In a last step, we average-out the kilometre-scale variability in the observed SMB (Agosta et al., 2012) by binning point values onto grid cells. For each grid cell containing multiple observations, we average all observations contained in the grid cell weighted by the time span of observations, and in the same way we weight-average the modelled values interpolated to observation locations. This way, we obtain consistent observed and modelled averaged values on grid cells.

We discard 66 observations beginning before 1979 and spanning less than 8 years. We also discard 12 observations for which the four surrounding grid cells fall in ocean and seven observations located at specific topographic features for which none of the four surrounding grid cell has an elevation difference of less than 200 m with respect to the actual location. After this, we retain 559 model–observation comparisons.

Figure 1 MAR SMB for the period 1979–2015: (a) mean annual SMB, with coloured dots showing the observed SMB values (shared colour scale); (b) standard deviation of annual SMB; (c) standard deviation divided by mean annual SMB; (d) difference between MAR and observed SMB on MAR grid cells, following the methodology detailed in Sect. 2.2.2. Magenta dots in panels (b) and (c) show the location of SMB observations. Solid grey lines are contours of surface height every 1000 m a.s.l. Latitude circles are −60, −70, and −80^∘ S, and longitude lines are from 145^∘ W to 145^∘ E by steps of 45^∘.

Figure 2 Modelled vs. observed SMB for sectors and transects as detailed in Table 1. RACMO2 outputs are bilinearly interpolated to the MAR grid. SMB values are first averaged on MAR grid cells (Sect. 2.2.2) then along a chosen grid direction (Fig. S2) or by elevation bins. Distance along the transect starts at the coast. Uncertainty of observed SMB (grey shaded area) is the standard deviation of observations contained in each grid cell (sub-grid variability), estimated as a function of the mean observed SMB (see Fig. S3). Despite SMB values corresponding to grid cell averages, we display one marker for each observation, with the x axis corresponding to the observation location along the transect or elevation. For observed SMB plots, markers with white faces are for bins containing fewer than 10 observations and black faces for bins containing more than 10 observations. Magenta bands mark grid cells in which more than 15 % of precipitation sublimates in the katabatic layers according to Grazioli et al. (2017). The map shows the main Antarctic basins: Antarctic Peninsula in purple, West Antarctic ice sheet in green, and East Antarctic ice sheet in orange. Ice shelves are mapped in blue and grounded islands in red, and the blue line shows the location of the grounding line.

3 Results

3.1 Evaluation of the modelled SMB

The large spatial Antarctic SMB gradients, shown in Fig. 1a as modelled by MAR forced by ERA-Interim for the period 1979–2015, coincide with a strong interannual variability (Fig. 1b), expressed by a standard deviation of ∼22 % of the mean SMB on average over the ice sheet (Fig. 1c). MAR SMB shows no systematic spatial bias (Fig. 1d), with a mean bias of 6 kg m⁻² yr⁻¹ (4 % of the mean observed SMB), as well as a very strong correlation with the observed SMB (R²=0.83, p value <0.01, computed on the logarithm of SMB values, as SMB distributions are log-normal). RACMO2 shows similar performance (mean bias of −3 kg m⁻² yr⁻¹, R²=0.86, computed on the logarithm of SMB as well).

The model–observation comparison by sectors (Fig. 2) reveals a good representation of the coast-to-plateau SMB gradients by both RCMs. MAR and RACMO2 are in good agreement despite MAR not including drifting snow processes whereas RACMO2 does, except in Ross–Marie Byrd Land and in Victoria Land where MAR simulates larger SMB than RACMO2. Another noticeable result is that MAR forced by ERA-Interim, JRA-55, and MERRA-2 gives very similar results for the SMB spatial pattern, not only at the observation locations (Fig. 2) but also at the ice sheet scale (comparisons of MAR SMB for different forcing reanalyses are shown in Fig. S4, with colour map scales 10 times smaller than in Fig. S5 in which MAR is compared to RACMO2). This is why we focus on MAR forced by ERA-Interim in the following.

We find no significant differences in the SMB simulated by MAR and RACMO2 when integrated over the ice sheet or its major basins (Table 2). SMB is driven by snowfall amounts, which are more than 10 times larger than other SMB components. Snow sublimation in RACMO2 is the sum of sublimation at the surface of the snowpack and of drifting snow sublimation and is approximately 50 % larger than in MAR, which only includes surface snow sublimation. However, surface snow sublimation alone is almost 2 times larger in MAR than in RACMO2 (Table 2 and spatial patterns shown in Fig. S6), which we investigate in the next section. Modelled surface melt is less than half of the sublimation amount; however liquid water almost entirely refreezes into the snowpack in both models (maps of MAR- and RACMO2-modelled melt amounts are shown in Fig. S7). Temporal variability in the SMB and its components is fully driven in both RCMs by the forcing reanalyses and are therefore strongly correlated with each other (time series shown in Fig. S8). We do not elaborate on the SMB temporal variability here as this aspect will be further detailed in a forthcoming study.

Table 2 Antarctic integrated SMB on average for 1979–2015 ± 1 standard deviation of annual values, in gigatonnes per year. Antarctic ice sheet (AIS) and basin geometry are based on Rignot basins (Shepherd et al., 2018), shown in Fig. 2. RACMO2 is bilinearly interpolated on the MAR grid and the same mask is applied to both models, with area given for this mask. SMB is computed as follows: MAR SMB = snowfall + rainfall − surface snow sublimation − run-off; RACMO2 SMB = snowfall + rainfall − surface snow sublimation − drifting snow sublimation − drifting snow transport − run-off.

3.2 Drifting snow transport features

Figure 3 Annual mean 10 m wind speed, curvature of elevation, and modelled SMB minus observed SMB for the four long transects: (a) Low Dome–Wilkes Land, (b) Zhongshan–Dome A, (c) Mawson–Lambert Glacier, and (d) Syowa–Dome F. Blue lines and colour shading are for MAR (ERA-Interim) outputs and red lines are for RACMO2 (ERA-Interim) outputs. Values are computed as in Fig. 2.

Figure 4 (a) Difference in SMB by grid cell (ΔSMB) between MAR (ERA-Interim) and observations for four transects (Law Dome–Wilkes Land, Zhongshan–Dome A, Mawson–Lambert Glacier, and Syowa–Dome F) vs. surface curvature on MAR grid. Curvature is shifted by ± 1 grid cell according to the maximum correlation with ΔSMB (Fig. S10). Linear regression through the origin is plotted with a pink dashed line. We excluded the regression of two outliers (dots with black outlines) and seven data for which MAR annual 10 m wind speed is lower than 7 m s⁻¹ (squares with black outlines). (b) Estimate of mean annual drifting snow transport based on a scaling of the curvature: drifting snow transport (kg m⁻² yr⁻¹) =α (10⁶ kg m⁻¹ yr⁻¹) × curvature (10⁻⁶ m⁻¹), with α=0 kg m⁻¹ yr⁻¹ for wind speed lower than 5 m s⁻¹, α=3700 10⁶ kg m⁻¹ yr⁻¹ for wind speed greater than 9 m s⁻¹, and α linearly increasing as a function of wind speed in between, around the 7 m s⁻¹ wind speed threshold. Wind speed is the annual mean of 10 m wind speed modelled by MAR (ERA-Interim). Coloured dots show the difference between MAR SMB and observed SMB with the same colour scale. (c) Mean annual drifting snow transport flux in RACMO2 on average for 1979–2015 (kg m⁻² yr⁻¹). Coloured dots show the difference between MAR SMB and observed SMB with the same colour scale.

Local fluctuations of the observed SMB around the smooth modelled SMB gradients are apparent along the four stake transects covering more than 500 km: Law Dome–Wilkes Land, Zhongshan–Dome A, Mawson–Lambert Glacier, and Syowa–Dome F. We related these fluctuations to drifting snow transport. Indeed, the snow eroded from the snowpack is loaded into the atmosphere, where it can sublimate and be transported by the wind. Katabatic winds blowing on the surface of the ice sheet result from the downslope gravity flow of cold, dense air. As a consequence, the surface wind divergence, which drives the snowdrift mass transport, is strongly related to the curvature of the topography, and both have similar spatial patterns (shown in Fig. S9). This is because slopes becoming steeper (crests, positive curvature) will lead to wind speed acceleration (positive wind divergence), and thus to drifting snow export (mass loss), whereas slopes becoming more gentle (valleys, negative curvature) will lead to wind speed deceleration (negative wind divergence), and thus to drifting snow deposit (mass gain).

To test our hypothesis, we computed the mean curvature of the MAR 35 km×35 km elevation field. In Fig. 3, we notice that both RCMs commonly exhibit an excess of accumulation on crests and a deficit of accumulation in valleys, in the range of ±40 kg m⁻² yr⁻¹. To quantify this curvature effect, we correlate MAR SMB bias (ΔSMB) with the curvature. For each transect, we apply a constant shift of ± one grid cell to the curvature in order to find the maximum correlation with ΔSMB. For three out of the four transects, we find only one shift for which the correlation is significant, and for the remaining transect (Syowa–Dome F) we find no significant correlation (Fig. S10). The sign and the amplitude of those shifts are in line with curvature being used as a proxy for wind divergence, as they are consistent with the Coriolis wind deflection westward of the topography gradient (detailed in Fig. S11). After applying those shifts, we find that the difference between modelled and observed SMB (kg m⁻² yr⁻¹) is scaled to approximately 3700 ± 1100 (10⁶ kg m⁻¹ yr⁻¹) times the curvature (10⁻⁶ m⁻¹), with a significant relationship (R²=0.41, Fig. 4a), when the mean annual 10 m wind speed (ws₁₀) is greater than 7 m s⁻¹. For lower wind speed (ws₁₀ <7 m s⁻¹), we no longer observe any relationship between model bias in SMB and curvature (horizontally aligned squares in Fig. 4a). This is consistent with the drifting snow transport process, which requires the wind speed to reach threshold values for the erosion to be initiated (Amory et al., 2015).

Hence, a large part of the discrepancies between modelled and observed SMB are explained by surface curvature when wind speed is sufficiently high, which we relate to the unresolved drifting snow transport in MAR. We are able to catch the drifting snow transport signal because drifting snow sublimation is negligible for the four studied transects, as they are located at high elevation, above 2000 m above sea level (a.s.l.), where the cold atmosphere has a low capacity to be loaded with moisture (see detailed analysis in Fig. S12). The moisture holding capacity of the atmospheric boundary layer is an upper bound for drifting snow sublimation and quickly tends to zero when the mean air temperature decreases below −30^∘C, which is the case along most of the transects, whereas the amplitude of observed SMB fluctuations around the smooth SMB gradient is independent of the temperature (Fig. S13).

Consequently, we propose that drifting snow transport fluxes (ds_tr) not resolved by MAR can be estimated as a scaling of curvature depending on wind speed: ds${}_{\mathrm{tr}}=\mathit{\alpha }\left({\mathrm{ws}}_{\mathrm{10}}\right)\cdot \mathrm{curvature}$ (Fig. 4b). The scaling factor α(ws₁₀) depends on wind thresholds to simulate the transition between no drifting snow transport for low wind speed (α=0 for ws₁₀<5 m s⁻¹) and drifting snow transport scaled to curvature for high wind speed (α=3700 10⁶ kg m⁻¹ yr⁻¹ for ws₁₀>9 m s⁻¹), with a linearly increasing scaling factor between 5 and 9 m s⁻¹ for a smooth transition around the 7 m s⁻¹ wind threshold defined above. That estimate of drifting snow transport fluxes shows little sensitivity to the choice of the wind thresholds and of the scaling factor (see fluxes summed over the ice sheet for different thresholds and scaling factors in Table S2). The spatial pattern of drifting snow transport we obtain is comparable to the one simulated by RACMO2 (Fig. 4c), except that it gives fluxes more than 3 times larger than in RACMO2 (Table S2, and note the different colour map scales between Fig. 4b and c). The drifting snow transport estimate consists in a redistribution of mass with negligible net mass loss over the AIS (total AIS mass gain of ∼75 Gt yr⁻¹ and total AIS mass loss of ∼80 Gt yr⁻¹; see Table S2).

Our drifting snow transport estimate gives a good constraint for drifting snow fluxes above 2000 m a.s.l., where low temperatures induce negligible atmospheric sublimation. As drifting snow transport is proportional to the amount of snow in suspension in the atmosphere, quantifying this flux also enables us to constrain the amount of snow eroded from the snowpack to the atmosphere, which drives drifting snow sublimation fluxes at lower elevation. This is of importance as drifting snow sublimation is a much larger mass sink than drifting snow transport over the whole ice sheet (Palm et al., 2017; Lenaerts et al., 2012a) but is still poorly constrained because observations are very scarce below 2000 m a.s.l. where it occurs.

Drifting snow sublimation included in RACMO2 and not in MAR moistens the surface atmospheric layers, consequently reducing the sublimation at the surface of the snowpack. This might explain the stronger surface snow sublimation in MAR than in RACMO2 (Table 2 and Fig. S6). However, drifting snow sublimation is a potentially larger mass sink than surface snow sublimation, as drifting snow particles are continuously ventilated and fully exposed to the ambient air. Consequently, by accounting for drifting snow in MAR we expect that the drifting snow sublimation mass sink could be enhanced at the expense of surface snow sublimation at the ice sheet margins.

3.3 Sublimation of precipitation in low-level atmosphere

Figure 5 The four maps show mass fluxes in kg m⁻² yr⁻¹ for the year 2015. (a) Difference in SMB between MAR and RACMO2. Blue lines delimitate areas where the SMB difference is 30 % greater than MAR SMB, with solid lines when MAR is greater than RACMO2 and dashed lines when MAR is lower than RACMO2. (b) Same as (a) but for the snowfall amounts at the ground. (c) Same as (a) but for the sublimation of precipitation in the atmospheric layers. (d) Same as (a) but for the maximum snowfall amount (equal to ground snowfall plus atmospheric sublimation). Locations of transects A1–A2 and B1–B2 extracted in Fig. 6 are shown in panels (b) and (d).

As described above, MAR and RACMO2 RCMs forced with ERA-Interim simulate similar spatial patterns for SMB compared to observations (Fig. 2). However, at the ice sheet scale, MAR and RACMO2 SMB show regional discrepancies (shown in Fig. 5a for 2015, and similar to the 1979–2015 mean shown in Fig. S5a), which are primarily the result of differences in simulated snowfall rates (Figs. 5b and S5b). We notice that areas where MAR snowfall is much lower than RACMO2 snowfall (Fig. 5b, blue dashed lines) coincide almost exactly with the pattern of precipitation that is able to sublimate in the low-level atmosphere according to Grazioli et al. (2017). In that study, the amount of atmospheric sublimation is quantified for the year 2015 using atmospheric modelling constrained with precipitation radar observations. Atmospheric sublimation happens because the katabatic surface air flux, moving from highly elevated inland plateau toward sea level, is subject to adiabatic compression when it moves downslope. This compression induces an increase in air temperature, which reduces relative humidity and drives sublimation rates in the lower troposphere (∼ first 1000 m above the ground), enhanced in the katabatic channels at the ice sheet margins.

To deepen this analysis, we re-ran MAR for the year 2015 in order to save the full atmosphere snowfall fields. From the daily 3-D snowfall amounts, we derived the atmospheric sublimation amount from the difference between the maximum snowfall and the ground snowfall in each atmospheric column, as in Grazioli et al. (2017). The same was done for RACMO2. We find that the atmospheric sublimation simulated by MAR (363 Gt for the year 2015 over the grounded ice sheet) is higher than estimated in Grazioli et al. (2017) (299 Gt after interpolation on the same mask) and much higher than simulated by RACMO2 (128 Gt, Fig. 5c). A major difference between MAR and RACMO2 is the advection of precipitation in the atmosphere: in MAR, precipitating particles are explicitly advected through the atmospheric layers until they reach the surface, while in RACMO2, precipitation is added to the surface without horizontal advection, and is able to interact with the atmosphere only in a single time step (6 min in this simulation). Consequently, atmospheric sublimation is likely to be underestimated in RACMO2.

We conclude, in agreement with Grazioli et al. (2017), that atmospheric sublimation is a major mass sink at the ice sheet margins in MAR, as for the year 2015 it represents 16 % of the snowfall loaded on the grounded ice sheet (12 % in Grazioli et al., 2017) and 26 % for areas below 1000 m a.s.l. (17 % in Grazioli et al., 2017).

It is noticeable that very few SMB observations are available in areas where Grazioli et al. (2017) identify low-level sublimation, marked by magenta bands in Fig. 2. Except for Ross–Marie Byrd Land, the only other areas where low-level sublimation is greater than 15 % of the total precipitation as defined by Grazioli et al. (2017) are close to Dumont d'Urville (coastal Adélie Land) and to Syowa (coastal Dronning Maud Land). In those areas the SMB amount is indeed larger in RACMO2 than in MAR and in observations. Both RCMs overestimate SMB around 2000 m a.s.l. in Dronning Maud Land and in Ross–Marie Byrd Land (Fig. 2), which could indicate katabatic channels not being resolved enough by the topography of the models.

3.4 Precipitation formation and advection

Figure 6 MAR- and RACMO2-simulated fields for the year 2015, extracted with a bilinear interpolation for (left) transect A1–A1 and (right) transect B1–B2 (locations shown in Fig. 5b and d). Each panel shows MAR fields (blue lines) and RACMO2 fields (red lines) for (a) surface height, in metres above sea level; (b) maximum snowfall amounts, equal to ground snowfall plus atmospheric sublimation, in kg m⁻² yr⁻¹; and (c) snowfall amounts at the ground, in kg m⁻² yr⁻¹. In (b) and (c), the thick black line is for the difference in snowfall between MAR and RACMO2 (MAR-RACMO2), with green-filled areas when MAR snowfall is larger than RACMO2 snowfall, and brown-filled areas when MAR snowfall is lower than RACMO2 snowfall (same convention as in Fig. 5); the dotted lines are for the atmospheric sublimation modelled by MAR (blue) and by RACMO2 (red), negative when it induces a decrease in precipitation; light coloured bands show crests (light blue, curvature of MAR topography greater than 0.005 10⁻⁶ m⁻¹) and valleys (light yellow, curvature of MAR topography lower than $-\mathrm{0.005}{\mathrm{10}}^{-\mathrm{6}}$ m⁻¹). The thick black arrows show the main 800 hPa wind direction during cyclonic activity.

Differences between MAR and RACMO2 snowfall fields are strongly reduced when considering the maximum snowfall amounts (before sublimation in the low-level atmosphere) rather than the ground snowfall amounts (Fig. 5b and d). However, MAR snowfall rates generally exceed those simulated by RACMO2, by more than 30 % on the lee side of the West AIS (Marie Byrd Land toward Ross ice shelf), on the lee side of the Transantarctic Mountains (Victoria Land), and close to crests at the ice sheet margins. MAR maximum snowfall rates are lower than simulated by RACMO2 windward of topographic barriers and in valleys at the ice sheet margins. This spatial pattern looks similar to the one obtained in RACMO2 when delaying the conversion of cloud ice–water into snow–rain (Fig. 3a of van Wessem et al., 2018). This change led to both ice and water clouds lasting longer in the atmosphere before precipitating and therefore being advected further towards the ice sheet interior (van Wessem et al., 2018).

For a more in-depth analysis, we extract MAR and RACMO2 snowfall rates on two transects at the ice sheet margins (Fig. 6), following the main wind direction during cyclonic activities (locations shown in Fig. 5b and d). On these transects the observed difference in maximum snowfall between MAR and RACMO2 is largely explained by a phase difference in the snowfall peaks windward of the topographic barriers, with MAR peaking closer to the crests than RACMO2 (Fig. 6b). This induces a wave-like pattern of precipitation difference strongly related to the shape of the topography, with larger snowfall amounts in MAR than in RACMO2 just windward of crests and lower snowfall amounts in MAR than in RACMO2 around windward valleys. At the ground, lower snowfall in MAR than in RACMO2 in valleys is amplified by low-level atmospheric sublimation, which peaks in katabatic channels (Fig. 6c).

Observations do not enable us to definitively discriminate one model against the other, but we observe a general tendency for MAR to overestimate accumulation on Ross–Marie Byrd Land and close to ice sheet summits (Dome C, Dome A, Dome F; see Figs. 1d and 2). Close to summits the wind is low, so a missing drifting snow transport process is an unlikely explanation for a positive bias in SMB modelled by MAR (Fig. 4b). Over the Greenland ice sheet, MAR tends to overestimate ice cores based on accumulation inland (Fettweis et al., 2017) while RACMO2 underestimates it (Noël et al., 2018).

We conclude that the differences in MAR and RACMO2 snowfall patterns are very likely related to differences in the advection of precipitation inland, which may arise from (i) the different advection of precipitating particles to the ground described in Sect. 3.3, (ii) different timing of precipitation formation (cloud–precipitation conversion thresholds), and/or (iii) different dynamical response to the topographic forcing, caused by either different dynamical cores or the different resolutions (the 27 km resolution in RACMO2 better resolves the ice sheet topography than the 35 km resolution in MAR).

4 Discussion and conclusion

In our study, we evaluate new estimates of the Antarctic SMB obtained with the polar RCM MAR ran for the first time for decades-long simulations at the scale of the whole AIS. We use model settings comparable to previous MAR simulations over Greenland (Fettweis et al., 2017) but with a specific upper atmosphere relaxation and new surface snow density and roughness length parameterisations. We present simulations of MAR forced by ERA-Interim, JRA-55, and MERRA-2 for the satellite era (1979–2015) in which we can rely on reanalyses products. Remarkably, MAR forced by those three reanalyses gives similar spatial and temporal SMB patterns. We also compare MAR with the latest simulations of the RCM RACMO2 forced by ERA-Interim (van Wessem et al., 2018). We find no significant differences between MAR and RACMO2 SMB when integrated on the AIS and its major basins (Table 2).

As the dominant feature of the Antarctic SMB is its strong coast-to-plateau gradient, we extract stake transects and sectors with large elevation ranges from the GLACIOCLIM-SAMBA SMB observational dataset. We show that both RCMs show similar performances when compared to observations, with a good representation of the SMB gradient (Fig. 2). But more importantly, we outline and quantify missing or underestimated processes in both RCMs.

Along stake transects, we relate 100 km scale fluctuations of observations around the smooth modelled SMB pattern to the shape of the ice sheet captured on the 35 km×35 km MAR grid. Both RCMs accumulate too much snow on crests, and not enough snow in valleys, as a result of drifting snow transport fluxes not included in MAR and probably underestimated in RACMO2 by a factor of 3 (Fig. 4). In the RACMO2.3p2 version used here, the modified drifting snow routine induced almost halved drifting snow transport and sublimation fluxes compared to the previous RACMO2.3p1 version (Lenaerts and van den Broeke, 2012). In a recent study combining satellite observation of drifting snow events and reanalysis products, Palm et al. (2017) estimate the drifting snow sublimation to be about ∼393 Gt yr⁻¹ over the AIS, vs. 181 Gt yr⁻¹ in RACMO2.3p1 and 102 Gt yr⁻¹ in RACMO2.3p2 (van Wessem et al., 2018). Consequently, observational constraints from our study and from Palm et al. (2017) both tend to confirm that drifting snow transport and sublimation fluxes are likely much larger than previous model-based estimates and need to be (better) resolved and constrained in climate models.

We also point out that MAR generally simulates larger SMB and snowfall amounts than RACMO2 inland, particularly on the lee side of the Transantarctic Mountains and on crests at the ice sheet margins, whereas MAR simulates lower snowfall than RACMO2 windward of mountain ranges and promontories. Sublimation of precipitating particles in low-level atmospheric layers is largely responsible for the significantly lower snowfall rates in MAR than in RACMO2 in valleys at the ice sheet margins. As precipitating snow particles have larger time residence in the atmosphere in MAR than in RACMO2 (Sect. 3.3), amounts of precipitation lost by sublimation in katabatic channels are more than twice as much in MAR as in RACMO2. The remaining spatial differences in snowfall between MAR and RACMO2 are attributed to differences in advection of precipitation, with snowfall particles being likely advected too far inland in MAR.

Atmospheric sublimation represents 429 Gt yr⁻¹ in MAR over the whole AIS (peninsula excluded) for the year 2015, 89 % of which is lost below 2000 m a.s.l. and 61 % below 1000 m a.s.l. This might be of importance for the mass balance of glacier drainage basins (SMB minus discharge; Rignot et al., 2008; Shepherd et al., 2018), as ice streams are typically channel-shaped areas affected by low-level sublimation of precipitation. Consequently, we note the importance of saving precipitation fluxes in models at least 1300 m above the ground for comparison with CloudSat products, but ideally at all model levels below 1500 m above the ground to be able to compute sublimation of precipitation in the low-level atmospheric layers. This will become a standard output in forthcoming MAR simulations.

We expect that accounting for drifting snow in MAR will lead to significant improvements in describing the Antarctic SMB and surface climate, as it will enable (1) a quantification of the drifting snow sublimation mass sink, (2) a more realistic representation of relative humidity and temperature in the boundary layer, and (3) an explicit modelling of the drifting snow transport from crests to valleys. Exploring the impact of horizontal and vertical model resolution on drifting snow estimates and on sublimation of precipitation in katabatic channels will also be of importance as those processes are related to the shape of the ice sheet and to the advection of precipitation in the atmosphere. The accuracy of the topography has to be considered as well, as digital elevation models are in constant improvement over the AIS (e.g. Slater et al., 2018) and should be regularly updated in climate models.