Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 2: Antarctica (1979–2016)

. We evaluate modelled Antarctic ice sheet (AIS) near-surface climate, surface mass balance (SMB) and surface energy balance (SEB) from the updated polar version of the regional atmospheric climate model, RACMO2 (1979– 2016). The updated model, referred to as RACMO2.3p2, in-corporates upper-air relaxation, a revised topography, tuned parameters in the cloud scheme to generate more precipitation towards the AIS interior and modiﬁed snow properties reducing drifting snow sublimation and increasing surface snowmelt.Comparisons of RACMO2 model output with several in-dependent observational data show that the existing biases in AIS temperature, radiative ﬂuxes and SMB components are further reduced with respect to the previous model version.


Introduction
Before being able to accurately predict future changes in the climate and surface mass balance (SMB) of the Antarctic ice sheet (AIS), it is required that its contemporary climate and SMB are realistically modelled. The interaction of the ice sheet with its atmospheric environment can be studied with regional climate models (RCMs) that are specifically adapted for use over the polar ice sheets of Greenland (Noël et al., 2016;Langen et al., 2017;Fettweis et al., 2017) and Antarctica Van Wessem et al., 2014b). Such polar RCMs explicitly calculate the individual components of the ice sheet SMB, such as precipitation, snow sublimation and surface meltwater run-off into the ocean. Combining integrated SMB with estimates of glacial discharge (D) then gives the mass balance of the ice sheet (MB) and its asso-ciated contribution to global sea level change (Rignot et al., 2011c;Shepherd et al., 2012).
Over Antarctica, precipitation is the dominating component of the SMB, contributing 91 % to the total (the sum of the absolute fluxes) mass budget (Van Wessem et al., 2014b). Snowfall rates are expected to rise in the future as a result of the increased moisture-holding capacity of a warming atmosphere (Ligtenberg et al., 2013;Palerme et al., 2017), which potentially compensates (Barrand et al., 2013) or amplifies (Winkelmann et al., 2012) changes in glacial discharge rates.
Surface melt rates along the coastal margins of the ice sheet can be significant and are further amplified by local climate features , possibly leading to ice shelf hydrofracturing (Van den Broeke et al., 2005a;Scambos et al., 2009;Kuipers Munneke et al., 2014), and the acceleration of grounded glaciers and related sea level rise (Rignot, 2004).
Polar RCMs are required to model the above processes realistically, and, in this study, the regional atmospheric climate model (RACMO2) is used at 27 km spatial resolution to simulate the climate and SMB of Antarctica. RACMO2 is able to accurately simulate climate variables over Antarctica (Van Lipzig et al., 2002;Reijmer et al., 2005;Van de Berg et al., 2006;Lenaerts et al., 2012b). However, substantial challenges remain. For instance, RACMO2 systematically underestimates snowfall, and hence the SMB, over the East Antarctic plateau. This negative bias was first identified in Van de Berg et al. (2006) but persisted in subsequent model versions at ∼ 10 % (Van Wessem et al., 2014b). In addition, biases in the net longwave radiation and sensible heat fluxes (Van Wessem et al., 2014a) resulted in biases in the energy available for melt (King et al., 2015) or in the surface energy budget (SEB), the sum of all energy fluxes at the surface. As a result, the previous model versions were systematically too cold over the ice sheet ( . In regions of melt, too little meltwater percolates into and refreezes in the snow column, resulting in too-low snow temperatures . Although this bias was substantially reduced, it still remains at −1.3 K in the latest model version (Van Wessem et al., 2014a).
One of the main causes of these discrepancies is shortcomings in the cloud microphysics. Biases in SMB and SEB are caused by too-thin clouds simulated over the AIS, resulting in too little snowfall, too much downwelling shortwave radiation and too little downwelling longwave radiation (Van Wessem et al., 2014a;King et al., 2015). In addition, these biases are potentially related to an unrealistic fractionation of ice and water content in these clouds, which significantly affects the sensitivity of the (cloud) radiative fluxes to changes in cloud content (King et al., 2015).
Another model issue is the relatively coarse horizontal resolution of RACMO2. A high spatial resolution is important for resolving the impact of the topography on the atmospheric motion in detail, realistically simulating topographyrelated processes such as katabatic winds and orographically forced precipitation (Genthon and Krinner, 2001). Since Van de Berg et al. (2006) and Van Lipzig et al. (2002) the horizontal resolution has been refined from 55 to 27 km in Lenaerts et al. (2012b), and recently climate simulations at 5.5 km resolution have been performed, focusing on specific regions such as the Antarctic Peninsula (AP) (Van Wessem et al., 2016), Dronning Maud Land (Lenaerts et al., 2016a), the West Antarctic coast (Lenaerts et al., 2017a) and Adèlie Land (Lenaerts et al., 2012a). Limitations in the model topography also result from biases in the source data sets used for aggregation (Liu et al., 2001;Bamber et al., 2009). These data sets are typically based on observational data from remote sensing techniques such as airborne (DiMarzio et al., 2007) and satellite measurements (Rignot et al., 2008), but frequent improvements are applied to these data sets as well Borsa et al., 2014).
Here, we discuss the effects of an update from RACMO2.3 version p1 to version p2, which addresses the model challenges presented above. To assess whether this update improves the modelled climate of the AIS in terms of SEB and SMB, we revisit several of the evaluations in previous studies. First, Sect. 2 discusses the model, the updates included and the observational data sets used for model evaluation. Section 3 presents the changes in terms of simulated SEB and SMB for the full ice sheet, as well as the evaluation by comparison with the observations. Section 4 presents the specific changes and evaluation for the 5.5 km model version applied to the AP. Section 5 then discusses remaining challenges and model limitations, followed by a summary and conclusions in Sect. 6. This paper is part of a tandem model evaluation over the Greenland (part 1) and Antarctic (this study) ice sheets.  . This version of the model is specifically applied to the polar regions by interactively coupling it to a multilayer snow model that calculates melt, refreezing, percolation and run-off of meltwater (Greuell and Konzelmann, 1994;Ettema et al., 2010). In addition, snow albedo is calculated through a prognostic scheme for snow grain size (Kuipers Munneke et al., 2011) while a drifting snow scheme simulates the interaction of the near-surface air with drifting snow (Déry and Yau, 1999;Lenaerts et al., 2010). Throughout this study, two model domains are addressed: RACMO2.3/ANT simulates the climate of the full Antarctic ice sheet (AIS), while RACMO2.3/AP is applied to the AP region specifically (Fig. 1).

Surface energy budget and surface mass balance
RACMO2.3 explicitly resolves surface melt by solving the surface energy budget (SEB; W m −2 ), defined as where fluxes directed towards the surface are defined positive, M is melt energy, SW and LW are the (upward-and downward) shortwave and longwave radiative fluxes, SHF and LHF are the sensible and latent turbulent heat fluxes and G s is the subsurface conductive heat flux. Excess energy at the surface is used to produce meltwater which can percolate through the snow column, where it is refrozen or retained. Ultimately, when the snowpack is saturated, the meltwater runs off to the ocean, representing a negative component in the surface mass balance (SMB) of the ice sheet. The SMB, in mm w.e. y −1 , is defined as where P tot represents total precipitation (snowfall, SN, plus rain, RA), SU surface (SU s ) plus drifting snow (SU ds ) sublimation, ER ds drifting snow erosion and RU meltwater runoff, the amount of liquid water (melt and rain) that is not retained or refrozen (RF) in the snowpack. Note that by defining SMB in this way, we include processes in the snowpack such as refreezing; this is formally referred to as the climatic SMB (Cogley et al., 2011).

Model updates
The model update includes small bug fixes and tuning of atmosphere and snow parameterisations, as summarised below. Special effort is made to synchronise the model updates to the Greenland and Antarctica model versions. Most of the updates presented below are therefore also implemented for Greenland, unless noted otherwise. Throughout this paper we will refer to the new RACMO2.3 version as RACMO2.3p2 and to the old version as RACMO2.3p1, where p stands for polar.
a. Previous studies found that the previous model version systematically underestimates snowfall and downwelling longwave radiation in the interior of the ice sheets of both Greenland  and Antarctica (Van Wessem et al., 2014b). Therefore, the critical cloud water and cloud ice content (l crit ) threshold, which governs the onset of effective precipitation formation T ra n s a n ta rc ti c M tn s . for mixed-phase and ice clouds, is increased in the following equation (adapted from ECMWF-IFS, 2008): (3) Here, A is a scaling value which is the cloud fraction for stratiform clouds and the updraught strength for convective clouds; c 0 is the coefficient for autoconversion of cloud ice/water into snow/rain; l cld is the total cloud ice and water content and BF e an enhancement factor for stratiform mixed phase clouds (ECMWF-IFS, 2008). The value of l crit is increased by a factor 2 for convective clouds and stratiform water/mixed phase clouds and by a factor 5 for stratiform ice clouds. This increase leads to both ice and water clouds lasting longer in the atmosphere before precipitating and therefore being advected further towards the ice sheet interior. As a result, we expect increased cloud cover for colder conditions and higher elevations and precipitation simulated further inland. Consequently, we expect to further decrease SEB and SMB biases in the AIS interior, in a way similar to Van Wessem et al. (2014a) and Van Wessem et al. (2014b).
The values of l crit adopted in RACMO2 were obtained after conducting a series of sensitivity experiments, i.e.
1-year simulations, over Greenland to test the dependence of precipitation formation efficiency, spatial distribution and cloud moisture content on l crit and other cloud-tuning parameters. From these experiments, we found a linear relationship between l crit for mixed and ice clouds, which is the vertical integrated cloud content. These new settings were then tested over Greenland and Antarctica for a longer period and proved to raise the accumulation and downwelling radiation fluxes in the interior ice sheets of Greenland (Noël et al., 2018) and Antarctica (see Sects. 3.1, 3.2 and 3.4.3). The new value of l crit remains well within the range of values previously presented in the literature (Lin et al., 1983).
b. The linear saltation snow load parameter (Eq. 24, Déry and Yau, 1999) used in the drifting snow sublimation scheme is halved, i.e. from 0.385 to 0.190 (Lenaerts et al., 2012b), effectively halving the horizontal drifting snow mass and sublimation fluxes, without changing the length or frequency of drifting snow events, which were well simulated (Lenaerts et al., 2012b). This makes the simulated drifting snow fluxes more in line with the limited drifting snow observations over Greenland (Lenaerts et al., 2014), and we include it for the Antarctic simulations as well.
c. The albedo of superimposed ice layers now follows Kuipers Munneke et al. (2011). Previously, the albedo used for superimposed ice was prescribed at 0.55, underestimating surface albedo and overestimating melt. Although the effect of this update is likely minor for Antarctica, where melting and the formation of superimposed ice layers in the current climate are rather limited in occurrence, it can result in better modelled surface albedo and melt locally and in future simulations. d. A number of minor model bugs, mostly related to the snow model and the resulting snowmelt fluxes, were fixed. Improvements were made in the snow grain evolution scheme, leading to faster snow grain metamorphism in the uppermost snow layer. As a result, albedo decay is enhanced, which in turn enhances snowmelt. Two changes applied to the RACMO2 version used over Greenland have not been included in RACMO2/ANT: (1) the size of refrozen snow grains is reduced from 2 to 1 mm (this value was already used in RACMO2/ANT), and (2) the model soot concentration, consisting of dust and black carbon impurities deposited on snow, has been halved from 0.1 to 0.05 ppmv. Over Antarctica, likely no soot is present and model soot content is kept at zero.

Initialisation and set-up
RACMO2.3 uses a vertical mesh of 40 hybrid model layers and a horizontal resolution of 27 km over the AIS and of 5.5 km over the AP. Both model versions are forced at the lateral boundaries ( Fig. 1) by ERA-Interim reanalysis data every 6 h (January 1979-December 2016, Dee et al., 2011). RACMO2.3 is not coupled to an ocean model, and sea surface temperature and sea ice cover are prescribed from the reanalysis. The lateral boundaries are now taken from higher resolution (0.75 • ) ERA-Interim fields than before (1.5 • ), mostly affecting sea ice near the ice sheet margin. RACMO2.3 is coupled to a snow model, a single column time-dependent model that describes the evolution of the firn layer based on a firn densification model (IMAU-FDM). It calculates firn density, temperature and liquid water content evolution based on forcing at the surface by surface temperature, accumulation and wind speed. Surface meltwater percolates into the model firn layer, where it can refreeze, be stored or percolate further down. The retention of meltwater is based on the tipping-bucket method (i.e. liquid water is stored in the first available layer and transported downwards only when it exceeds the maximum capillary retention). Liquid water that reaches the bottom of the firn layer is removed as run-off. More details on the snow model can be found in Ettema et al. (2010), Ligtenberg et al. (2011) andKuipers Munneke et al. (2015). Both RACMO2.3/ANT and RACMO2.3/AP are initialised on 1 January 1979; the initial firn pack is taken from a simulation with an offline IMAU-FDM (Ligtenberg et al., 2011), which has the same physics but a higher vertical resolution than the internal snow model and is driven by the previous RACMO2.3 climatologies for both model versions (Van Wessem et al., 2014b. 2.4.1 Specific RACMO2.3p2/ANT set-up changes a. The model topography in earlier model versions was aggregated from the digital elevation model (DEM) from Liu et al. (2001). The updated topography is based on the (1 km spatial resolution) DEM from Bamber et al. (2009), which is derived from a combination of satellite radar and laser data. As a result, the topography in the coastal regions is now better resolved , likely resulting in better modelled topography-related variables, such as melt, katabatic winds and orographic precipitation.
b. To better simulate SMB interannual variability in RACMO2.3p2, upper-atmosphere relaxation (UAR or nudging) of temperature and wind fields is applied every 6 h for the model atmospheric levels above 600 hPa (Van de Berg and Medley, 2016). UAR is not applied to atmospheric humidity fields in order not to alter clouds and precipitation formation in RACMO2. Additional details  Fig. 1).
c. Upper-air relaxation in RACMO2.3/AP is not used, as the model domain is smaller therefore better constrained at the lateral boundaries, and because the wind relaxation would smooth out orographic precipitation over the detailed AP topography.

Observational data
We evaluate the updated model by revisiting the comparisons of modelled SMB components and SEB in previous studies (Van Wessem et al., 2014a. These observational data include SEB data from automatic weather stations (AWSs), 10 m snow temperature, in situ SMB observations, accumulation radar profiles and solid ice discharge estimates. Detailed descriptions of these methods are provided in their respective papers. Here, we describe updates to the observational data sets only.

Automatic weather stations
We evaluate modelled near-surface wind, temperature and SEB components using nine AWSs located in Dronning Maud Land (DML) in different climate regimes. These AWSs measure all four radiation components as well as humidity and snow temperature, which are used to force an energy balance model (EBM), effectively closing the SEB at AWS locations. A summary of the location and data records of the AWSs is provided in Table 1 and locations are shown in Fig. 1. We use the same time period  for evaluation as in Van Wessem et al. (2014a) and did not extend the few AWSs that have longer records. For more details on the EBM see Van den Broeke et al. (2005a, b) and Reijmer and Oerlemans (2002). All data from AWSs are monthly averaged and compared with data from the same months from RACMO2.3. In addition to these AWSs, in order to evaluate the modelled melt fluxes, we have included SEB data from Neumayer station in DML that span a longer time period  ( Van den Broeke et al., 2009). We use these data to force the EBM and test the sensitivity of the calculated melt to two different surface roughness lengths: z 0 = 28 (a commonly used setting) mm and z 0 = 1 mm (standard setting in the EBM).

In situ SMB and temperature observations
We compare the updated RACMO2.3/ANT modelled SMB with 3234 in situ SMB observations ( Fig. 1) described by Favier et al. (2013). For modelled surface temperature (T s ) evaluation, we use the 64 10 m snow temperature observations ( Fig. 1) used in Van Wessem et al. (2014a), which at this depth are representative of annual average surface temperature.

Accumulation radar
In order to evaluate model performance in capturing the temporal variability in SMB, we use radar-derived annual accumulation rates from 1980 to 2009 over much of the Thwaites glacier catchment area generated by Medley et al. (2013). The accumulation rates were derived from the Center for Remote Sensing of Ice Sheets snow radar system as part of NASA Operation IceBridge (Leuschen, 2014). Similarly, we calculate accumulation rates using snow radar data from two additional surveys ( Fig. 1): the western Getz ice shelf (17 October 2011) and the Ronne ice shelf (15 Novem-ber 2015) near the Institute Ice Stream grounding zone. The areas covered by these new surveys are much smaller than the Thwaites glacier survey (∼ 1600 km) but are larger than the RACMO2.3 grid cell width (Getz: ∼ 75 km, Ronne: ∼ 50 km), averaging out much of the spatial variability. We use a method described in Medley et al. (2015) which allows for spatial variations in the firn density profile in both radar-depth estimation as well as conversion to water equivalence. The method iteratively solves for a depth-density profile for each measurement that is consistent with the radarderived accumulation rate as well as long-term modelled 2 m air temperature from MERRA-2 and an initial density of 350 kg m −3 . For the Ronne ice shelf survey, we are able to generate a 30-year time series ; however, we are only able to record 10 years of accumulation over the Getz ice shelf, which is due to its high accumulation rate (∼ 900 mm w.e. y −1 ). Thus, the model evaluation over the Getz is less robust than over the Ronne and Thwaites surveys.

GRACE
The Gravity Recovery and Climate Experiment (GRACE) mission provides monthly observations of mass changes across the AIS at a resolution of ∼ 300 km. To assess temporal simulated SMB variability, we regionally compare cumulative SMB with the GRACE observations on a regional scale. The GRACE data (CSR RL05) are processed as in Van Wessem et al. (2014b). Mass anomalies are assigned to predefined drainage basins . To better capture the higher SMB variability in the marginal zones, these basins are split up into coastal and interior sub-basins. The assigned anomalies are then converted to pseudo-GRACE observations and adjusted until the differences with the actual GRACE observations are minimised in a least-squares sense. GRACE measures the integral sum of mass anomalies along the satellite orbit. Atmospheric and oceanic mass variability is removed by the processing centres using numerical models. However, residual signals may remain due to model deficiencies, mostly at high frequencies. To reduce the aliasing effects of these signals, both the GRACE and RACMO2.3 data are smoothed with a 7-month moving average filter. As GRACE measures the sum of SMB and ice dynamics, which generally act at slower timescales, all time series are detrended and a quadratic acceleration term is removed.

Drifting snow in Adélie Land
Observations of horizontal drifting snow transport fluxes (TR ds ) were performed in Adélie Land, East Antarctica ( Fig. 1), where surface atmospheric conditions are well monitored at the permanent French Dumont d'Urville station (Favier et al., 2011). The coastal region is characterised by frequent strong katabatic winds that are regularly associated with aeolian snow transport events (Trouvilliez et al., 2014). The evaluation data set includes snowdrift acoustic measurements collected using second-generation FlowCapt ™ devices during 2013 at a coastal location of Adélie Land, as described in Amory et al. (2017). The measurement system consists of two 1 m long tubes superimposed vertically to sample the first two metres above the snow surface, which largely represent the total snow mass flux (Mann et al., 2000). Monthly values of horizontal snowdrift transport from 0 to 2 m are computed from half-hourly estimates of the snow mass flux and are compared with snowdrift transport fluxes from the drifting snow routine of RACMO2.3. The measured flux represents a minimum estimate of the total column flux.

QuikSCAT melt fluxes
We use estimates of surface meltwater production from the satellite radar backscatter time series from QuikSCAT (QS-CAT), which is calibrated with SEB-derived melt flux observations (Trusel et al., 2013). We compare these with modelled melt fluxes for the period 2000-2009 for the ice shelf regions denoted in Fig. 1.

CloudSat-CALIPSO
To evaluate the modelled downwelling radiative fluxes we use a modified version of release 04 (R04) of the CloudSat Level 2B Fluxes and Heating Rates (2B-FLXHR-LIDAR) product (Henderson et al., 2013) with specific adaptations for the polar atmosphere (Van Tricht et al., 2016b), combining cloud properties retrieved by CloudSat-CALIPSO (C-C), reanalysis data from ERA-Interim and surface properties, which drives a broadband radiative transfer model. Data from 2007 to 2010 are available at a horizontal resolution of 2 • × 1 • as described in Lenaerts et al. (2017b) and Matus and L'Ecuyer (2017). The performance of this product over polar regions has been evaluated by statistical comparison with AWS measurements both in the Arctic as well as in the Antarctic (Van Tricht et al., 2016a, b). Root mean square deviations (RMSDs) of the computed fluxes and the AWS measurements from the four Antarctic AWSs are used as uncertainties in the fluxes used in this study.

Upper-air profiles
To evaluate modelled (upper) atmospheric conditions we use radiosonde data from Kohnen Station, close to AWS 9, maintained by the Alfred Wegener Institute (AWI). As part of austral summer campaigns at Kohnen Station (75 • S, 0 • E, 2892 m a.s.l.), in 2005/2006 and 2013/2014, radiosonde measurements were performed four times a day (00:00, 06:00, 12:00 and 18:00 UTC). Radiosondes of type RS92-SGP manufactured by Vaisala, Finland, were used. The sondes carried sensors to measure pressure, temperature and humidity and a GPS receiver to capture wind speed and wind direction. At 06:00 and 18:00 UTC a 200 g balloon was launched.

Discharge estimates
To evaluate RACMO2.3/AP modelled SMB we use the solid ice discharge estimates from Wuite et al. (2015), as also used in Van Wessem et al. (2016). We added similar discharge estimates for the George VI ice shelf (Hogg et al., 2017).
3 Results: full ice sheet at 27 km 3.1 Changes in modelled SEB Figure 2 shows the difference in modelled annual mean (1979-2014) SEB and near-surface variables between RACMO2.3p2 and RACMO2.3p1. Figure 2a shows a relatively uniform increase in LW d of up to 10 W m −2 , a result of the changes in the model cloud scheme, leading to more clouds inland. The increase is largest over the sea ice and along the coastal margins and the escarpment zone of the ice sheet, where cloud coverage peaks. Figure 2b illustrates that surface temperature T s increased uniformly as a result of the increase in LW d . At the East Antarctic plateau changes are smaller. Differences in V 10 m are small and most pronounced in regions of steep topography, where changes in the model surface topography are largest but are also due to changes in the near-surface temperature inversion T inv = T s − T 2 m as a result of the changes in T s . Finally, Fig. 2d shows the resulting changes in the sensible heat flux (SHF), which are closely related to the changes in LW net , causing SHF to drop as the surface temperature inversion is weaker. Figure 2 also shows sharp local changes over polynyas between the iceshelf edge and the sea ice. These differences are mostly related to changes in the ice mask, resulting from the updated surface topography but are also related to the higherresolution ocean and sea ice boundary forcing. The latter change caused sea ice to be better resolved along the coastal margins, with associated changes in the surface fluxes. Figure 3 shows the related changes in the SMB components. First, the patterns in the change in total precipitation (snowfall and rainfall, Fig. 3a) are comparable to the changes in LW d and T s : there is a migration of precipitation from the ocean and coastal margin towards the interior of the ice sheet due to the updated cloud scheme. However, in  Figure 3b illustrates that sublimation fluxes are smaller locally, removing up to 200 mm w.e. y −1 less mass. This is a result of the lowered snow drift saltation parameter, causing drifting snow sublimation SU ds to be less efficient in the coastal and escarpment region where wind speeds and SU ds are the largest. Over the flat ice shelves, as well as some regions where the new topography has resulted in lower elevations, warmer conditions have resulted in a small increase in (surface) sublimation. Surface snowmelt fluxes (Fig. 3c) increased by ∼ 100 mm w.e. y −1 over the coastal region and ice shelves. The changes are mainly caused by the faster snow grain growth, which lowers surface albedo and enhances summer snowmelt. Surface melt has no influence on the SMB, as nearly all meltwater refreezes in the firn. Ultimately, changes in SMB (Fig. 3d) are dominated by changes in P tot , which are the largest. Table 2 summarises the RACMO2.3 integrated values of the SMB components, calculated for the model ice mask including ice shelves, but excluding the AP. In RACMO2.3p2/ANT, integrated SMB amounts to 2229 Gt y −1 with an interannual variability of σ = 109 Gt y −1 , which is an increase of 69 Gt y −1 (3.2 %). Changes in precipitation are mainly caused by a redistribution over the ice sheet; integrated changes over the total ice sheet are small: P tot increased slightly by 14 Gt y −1 (0.5 %) to 2400 ± 109 Gt y −1 . The increase in SMB is mostly caused by a reduction in drifting snow sublimation SU ds , which dropped by 79 Gt y −1 , from 181 ± 9 to 102 ± 5 Gt y −1 . An increase of 22 Gt y −1 in surface sublimation partly compensates for this decrease and total SU decreased by 56 to 161 ± 7 Gt y −1 . Snowmelt increased significantly from 36 ± 17 to 71 ± 28 Gt y −1 , an increase of 97 %. The increase does not affect the SMB, as all extra meltwater (and rain) is modelled to refreeze in the snow (71 ± 28 Gt y −1 ). The remaining SMB components (rainfall, ER ds and run-off) did not change significantly.

Changes in integrated SMB
Integrated SMB for the grounded ice sheet (GIS), which is based on the Rignot et al. (2011a, b) drainage basin definitions excluding the AP, increased by 103 to 1885 Gt y −1 , which is a bigger increase than when including the ice shelves. This difference is caused by a redistribution of snowfall towards the interior of the ice sheet (see Fig. 3a, d) and because changes in SU ds are mainly present over grounded ice (Fig. 3b). This effect is most profound over the East Antarctic ice sheet (EAIS), where the SMB increased by  79 Gt y −1 (7.5 %), compared to an increase of only 17 Gt y −1 (2.7 %) over the West Antarctic ice sheet (WAIS).

Automatic weather stations
Figure 4 compares monthly modelled SEB components with AWS observations. All components show consistent improvements. For net longwave radiation (LW net ), the determination coefficient (r 2 ) increased from 0.69 to 0.77, the bias decreased from −6.9 to −4.5 W m −2 and the RMSD decreased from 7.8 to 7.0 W m −2 . Similar but less pronounced improvements are found for the other components, most notably SHF. Improved correlation is mostly caused by the inclusion of upper-air relaxation in the model: temporal variability is now better constrained (Van de Berg and Medley, 2016). Improvements in bias and RMSD are mainly caused by changes in the cloud scheme. Figure 4e shows small improvements for V 10 m . Temporal variability has improved in a similar way to the SEB, but improvements in slope, bias and RMSD are mostly unrelated to the improvements in the SEB components. These changes are related to the new model topography: the regression slope of V 10 m went up considerably from 0.43 to 0.56. This poor fit is generally caused by the relatively coarse horizontal resolution, resulting in an underestimation of surface slope in steep areas and an overestimation in flat areas (Reijmer and Oerlemans, 2002;Van Wessem et al., 2014a). Therefore, katabatic winds that are strongly related to the surface slope are too weak/strong in sloped/flat regions. Even though the horizontal resolution in RACMO2.3/ANT is unchanged, surface slope is likely better represented in the new model topogra- phy. The improvement in T s (see below) also partly leads to a better representation of the surface temperature inversion T inv , which affects the strength of the katabatic wind forcing. Figure 4f shows modelled T s , as a function of observed T s calculated by closing the SEB. A general shift to higher temperatures is seen for the entire temperature range, but temperature is still underestimated in most months. This warming results in better statistics overall, although at some locations temperature now is/remains overestimated. This overestimation overcompensates some of the biases, which is reflected in the change of the RMSD: while the absolute bias suggests a near-perfect match of modelled and observed T s (−0.14 K), RMSD shows a more modest improvement by 0.26 K.

10 m firn temperature
To evaluate changes in T s over a larger area, Fig. 5 shows the comparison of RACMO2.3/ANT modelled 10 m firn temperature with 64 observations. Slope and spatial correlation remained similar, but the bias and RMSD decreased significantly, by 1 and 0.35 K respectively. The relatively small decrease in the RMSD is caused by the persistent positive temperature bias for the relatively cold East Antarctic plateau. At these locations the updated model now overestimates temperature by up to 3 • , likely related to the slight overestimation of downwelling longwave radiation with the updated cloud scheme. Figure 6a shows the LW d and SW d bias and the observations (right axis) for all model ice sheet grid points averaged in nine surface elevation bins. For SW d the bias decreased by about ∼ 5 W m −2 for all elevation bins. Generally, over the AIS SW d remains overestimated, but the new model version simulates these fluxes better. Improvements in LW d are also evident: here the old model version underestimated LW d for all elevation bins, but more so at lower elevations (< 1500 m). In the updated model the bias for these lower elevation bins, mostly representing the ice shelves and the WAIS, is reduced. However, for the elevation bins >2000 m, which represent the East Antarctic plateau, LW d is now overestimated by ∼ 4 W m −2 . Overall, the downwelling radiative fluxes are likely better represented in the new model version, but errors in these fluxes are large. Only for intermediate elevations, the third to fifth bins, the changes appear larger than the uncertainty in the fluxes, which is based on a comparison with ground-based AWS observations (Van Tricht et al., 2016a).

CloudSat-CALIPSO
Unfortunately, Fig. 6b also shows that a systematic overestimation of both cloud ice and water is present in the new model version. Apparently, there is an error that compensates the biases in the downwelling radiative fluxes: by simulating clouds that are optically too thick, biases in the fluxes and also in P tot have decreased. The general biases in these cloud variables appear relatively large, especially for the coastal bins. However, the biggest fraction of grid points is located in the high-elevation bins as shown in the bar chart in Fig. 6a, where biases are the lowest.

Kohnen Station radiosondes
To assess upper-atmosphere conditions, Fig. 7 compares modelled temperature, relative humidity and wind speed profiles with radiosonde measurements conducted at Kohnen Station (75 • S, 0 • E, 2892 m a.s.l.). Here, we only show the average daily values for January 2014; the comparisons for other years give similar results. Figure 7 shows a consistent improvement in modelled temperature, relative humidity and wind speed. RMSD decreased significantly for all three variables and at all height levels in the atmosphere, a direct result of the upper-air relaxation included in RACMO2.3p2. The improvement is largest for wind speed; for temperature the changes are more variable near the surface (0-1000 m). It is interesting to see the significant improvement in simulated wind speed and temperature at the tropopause inversion (∼ 8000 m). Here, RMSD is largest, which is related to the maximum in wind speed and temperature at this level and the difficulty in simulating the accurate vertical location of the inversion. Figure 8 shows the SMB bias (model-observation) for both model versions averaged in nine surface elevation bins. The uniform increase in cloud cover and P tot , accompanied by the decrease in drifting snow sublimation, increased SMB in all elevation bins. Only in the lowest bin, which represents the ice shelves, SMB decreased because precipitation now falls at higher elevations. The increase in SMB leads to a relative SMB bias that is generally below 5 % for all bins above 250 m elevation. For some bins, SMB is even slightly overestimated. This change is consistent with those in temperature (Figs. 4f and 5b): more clouds in the interior result in higher temperatures and snowfall rates. The fifth and sixth elevation bins (1750-2750 m) show the largest improvements: negative biases decreased from −9 and −16 % to 1 and −6 %. In these bins the improvements are larger than the uncertainty margins and therefore significant. In the coastal and escarpment zone of mainly West Antarctica (the 250-1750 bins), a small overestimation of SMB is now suggested, but the bias is not significantly different from zero. Table 3 compares simulated annual SMB with annual accumulation rates derived from the snow radar in three regions in West Antarctica (Fig. 1). This comparison includes two new but small (unpublished) surveys over Getz and Ronne ice shelves. Correlations for annual values between the measured and simulated SMB have significantly improved for the three surveys, clearly resulting from the better-constrained interannual variability (Van de Berg and Medley, 2016). For the largest of the survey areas, Thwaites glacier, correlation r has improved from 0.68 to 0.90, comparable to the reanalyses in this region. Modelled mean accumulation rates over Thwaites remain slightly underestimated, but changes are small.

Snow accumulation radar
For the ice shelf comparisons correlation (r) increased significantly from ∼ 0.5 to ∼ 0.8. These values are lower than the Thwaites comparison because (1) the accumulation rate over the Ronne ice shelf survey is lower (∼ 190 mm w.e. y −1 ), resulting in a lower signal-to-noise ratio in the radar time series, and (2) the Getz record is shorter.  Figure 9 shows detrended and decelerated mass anomalies from GRACE, RACMO2.3p2 and RACMO2.3p1 for the West and East Antarctica ice sheets, the AP and for the AIS. Changes from RACMO2.3p1 to RACMO2.3p2 are small and likely insignificant as changes fall within the uncertainties of both products, but generally there is a slight improvement. For West Antarctica (Fig. 9a) and East Antarctica (Fig. 9c) the RMSD and correlation both improved slightly. For the AP changes are smaller and not significant. When these regions are combined, we find a slight improvement in the mass anomalies for the whole of the Antarctic ice sheet, but our main conclusion is that both model versions realistically simulate the Antarctic mass anomalies.  Figure 10a shows monthly cumulative horizontal drifting snow transport fluxes (TR ds ) for 2013 from Amory et al. (2017), RACMO2.3p2 and RACMO2.3p1. Modelled TR ds is computed for the full drifting snow layer, while the observed flux is measured for the first 2 m above the surface. Therefore, although the bulk of the flux is contained within the vicinity of the surface (Mann et al., 2000;Lenaerts et al., 2010), the observations represent a lower limit of TR ds . Figure 10a illustrates a significant decrease in modelled TR ds for RACMO2.3p2, as a result of the lowered linear saltation snow load parameter (Déry and Yau, 1999). There is a good match, albeit an average underestimation by RACMO2.3p2 of 5.4 tons, with observed TR ds . However, the underestimation of TR ds partly results from an underestimation of V 10 m (Fig. 10b) by ∼ 2.3 m s −1 on average. Correcting for this bias in wind speed by linearly regressing modelled wind speed and TR ds reveals a likely improved match with the observations. Even though this is a simple first-order correction, it suggests that RACMO2.3p1 TR ds fluxes were too high.

Surface melt
To assess modelled surface melt fluxes, Fig. 11 correlates annual average (2000-2009) melt fluxes from RACMO2.3 and the QSCAT satellite for the whole AIS and six regions. The increase in snowmelt for the total AIS (inset) is evident and the negative bias decreased from −61 to −15 Gt y −1 , while the determination coefficient r 2 increased from 0.75 to 0.81. For the six regions, melt is also better represented, although a few years remain in which melt is underestimated. The largest underestimation is shown for Wilkins ice shelf. Here, observed melt fluxes may be overestimated due to extensive melt ponding and/or saturated firn conditions in this region (Trusel et al., 2013;Välisuo et al., 2014), a feature which negatively affects the QSCAT retrievals. Figure 12 compares cumulative melt for 1990-2015 by RACMO2.3/ANT, and the melt calculated by an EBM (with two roughness length settings) forced with meteorological data from Neumayer station (Reijmer and Hock, 2008;König-Langlo, 2013). Updated model melt correlates best with the melt calculated by the EBM, but the total melt remains somewhat underestimated, consistent with the comparisons with QSCAT. The timing of the melt events is captured well: many of the strong melt summers, e.g. 1994/1995, 2000/2001 and 2013/2014, are now accurately modelled. Year Year Year Year 4 Results: the Antarctic Peninsula at 5.5 km 4.1 Changes in modelled SMB and SEB For RACMO2.3/AP at 5.5 km resolution, changes in SMB and SEB components are mostly consistent with those at a coarser resolution (Fig. 3). This means that for the AP we see an increase in P tot (Fig. 3a) and hence in the SMB (Fig. 3d) over the western slopes, a small decrease in P tot over the eastern ice shelves, a decrease in sublimation (Fig. 3b) and a considerable increase in snowmelt (Fig. 3d) where the RACMO2.3p2/ANT snowmelt is comparable to that simulated by RACMO2.3p2/AP (Fig. 11). There are only slight local differences in topography-related variables such as precipitation due to the different interpolation setting used to ag-gregate the model topography, shifting topographic features approximately one grid box northwards. However, as we will show in the next section, this hardly affects integrated SMB estimates. Table 4 shows the integrated changes for RACMO2.3/AP. Total integrated SMB (1979SMB ( -2014 increased by 16 Gt y −1 to 366 ± 58 Gt y −1 . This increase is mainly caused by changes in precipitation (+14 Gt y −1 ). Even though there are large local changes due to reprojection of the RACMO2.3/AP topography, it hardly affects the total integrated SMB. Changes in SU ds and the resulting SU are similar to those at a coarser resolution. While SU ds dropped by 4 Gt y −1 (44 %), SU s increased by 1 Gt y −1 and total sublimation SU now amounts to 8 ± 2 Gt y −1 , a decrease of 3 Gt y −1 . Snowmelt increased significantly (+10 Gt y −1 ) but most of the meltwater increase Bias Bias Figure 11. Modelled RACMO2.3p2 (blue), RACMO2.3p1 (red) and RACMO2.3p2/AP (green) as a function of observed QSCAT (Trusel et al., 2013)  is refrozen in the snowpack (+9 Gt y −1 ). As a result runoff fluxes remain similar to the fluxes simulated in the old version. Relatively, most of the SMB increase comes from an increase in the eastern Antarctic Peninsula (EAP) SMB (+5.3 %), while the SMB in the western Antarctic Peninsula (WAP) increased by 3.9 %.

Evaluation of modelled SMB
Changes in RACMO2.3/AP are generally small and most changes are also visible at a coarser resolution, and only one previous evaluation with new observational data is revisited. Figure 13 shows RACMO2.3/AP average SMB (1979SMB ( -2016 as a function of glaciers draining in Larsen B (a, Wuite et al., 2015) and George VI ice shelves (b, Hogg et al., 2017) for 1995/1996, when these were assumed to be in balance. Changes due to the model update are small yet different for the two regions. For Larsen B a consistent improvement in modelled SMB is seen: r 2 , bias and RMSD all marginally improved. For George VI SMB values are an order of magnitude larger. The bias improved for the  glaciers where the SMB was formerly underestimated, but worsened for those where it was overestimated. RMSD worsened slightly, but the slope of the fit and the determination coefficient r 2 improved. Overall changes are small and fall well within the uncertainty/interannual variability of the discharge estimates/modelled SMB.  Figure 13. Modelled RACMO2.3p2/AP (blue) and RACMO2.3p1/AP (red) integrated average (1979-2016) SMB as a function of glacial discharge estimates from (a) the Larsen B embayment (Wuite et al., 2015) and (b) the George VI embayment (Hogg et al., 2017). Horizontal error bars represent the uncertainty in the discharge estimate, vertical bars the interannual variability of RACMO2.3.

Clouds
In this update of RACMO2.3 the near-surface climate in terms of the SEB and SMB has generally improved, but several limitations and challenges remain. For instance, changes in cloud scheme parameters governing the conversion of clouds into precipitation have altered the distribution and quantity of clouds over the ice sheet, which in turn has led to improvements in the near-surface climate, as presented in this study. Several of the variables related to the distribution of clouds over the East Antarctic plateau (LW d , precipitation and surface temperature T s ) are now overestimated in the interior, albeit slightly so. Therefore, further improving these biases will require a different approach.
For instance, when the model is used at increasingly higher spatial resolution, it is important whether precipitation is handled prognostically by the model. While it is now assumed that precipitation reaches the surface instantaneously, ideally precipitation should be modelled as a prognostic variable in order to capture its fall time and horizontal displacement. This could significantly affect local precipitation patterns and hence SMB in regions of relatively high wind speeds and snowfall rates, such as the AP and the WAIS.
Several studies suggest that it is not necessarily the amount of cloud but the fractionation of cloud ice and cloud water content that is important for the radiative fluxes: King et al. (2015) showed that significant biases exist in relation to the ice/water content of clouds, and that a realistic simulation mainly requires accurate modelling of the presence or absence of clouds at a particular level in the atmosphere. Matus and L'Ecuyer (2017) and Lenaerts et al. (2017b) show that climate models typically exhibit significant biases in the simulation of clouds and radiation in the polar regions, and that future model developments should focus on the microphysi-cal properties of clouds and their radiative impact. We have also shown in Fig. 6 that, even though the simulation of the radiative fluxes has improved at the surface, RACMO2.3p2 considerably overestimates cloud ice and cloud water, suggesting that compensating errors exist due to a cloud radiation scheme that is not active enough. Further model improvement efforts should also focus on the surface layer turbulent mixing scheme ( Van de Berg et al., 2006).

Cold bias, meltwater and hydrostatic assumption
The updated RACMO2.3/ANT and RACMO2.3/AP remain somewhat too cold at the surface, especially at lower elevations. This is partly due to biases in the SEB and cloud cover, but also likely due to the limited horizontal resolution. A limited horizontal resolution affects the near-surface winds, where biases influence the simulation of SHF and the associated mixing of warmer air to the surface. This limitation is addressed by the higher resolution in RACMO2.3/AP, but a further increase in horizontal resolution moves the model towards the limits of the hydrostatic assumption, and it would be necessary to resolve vertical accelerations and internal and gravity waves. A non-hydrostatic model formulation should be used if these features are to be explored further.
The cold bias in the model 10 m firn temperature mostly affects more temperate regions of the ice sheet (Fig. 5), which suggests that the bias could be related to the refreezing of meltwater. Meltwater percolates into the snow column, refreezes and raises snow temperature. Underestimating surface melt or neglecting meltwater penetration features (Harper et al., 2012;Ligtenberg et al., 2014) could result in too-low snow temperature. Another explanation is improper modelling of the albedo-melt feedback in regions of sustained local melt, such as the melt induced by foehn winds in the AP (Luckman et al., 2014;Van Wessem et al., 2016), or the wind-albedo interaction over blue-ice areas in DML . Further efforts to improve modelled melt fluxes should therefore focus on these features.

Summary and conclusions
In this study we evaluated the modelled Antarctic ice sheet (AIS) climate, surface mass balance (SMB) and surface energy balance (SEB) in RACMO2.3p2, the latest polar version of the regional atmospheric climate model RACMO2 . This model version is applied at 27 km horizontal resolution to the full AIS and at 5.5 km resolution to the Antarctic Peninsula (AP). Updates to the previous model version include additional upper-air relaxation by ERA-Interim reanalyses data, a revised topography, tuned parameters in the cloud scheme that generate more precipitation towards the AIS interior and modified snow properties reducing drifting snow sublimation and increasing simulated surface snowmelt. The updates in the 5.5 km AP simulation are similar but do not include a new topography or upper-air relaxation.
The updated model simulates more clouds towards the interior of the ice sheet, increasing downwelling longwave radiation and snowfall rates. Drifting snow sublimation rates decreased by 43 % and surface snowmelt doubled. In total, these changes led to an integrated SMB for the ice sheet, including ice shelves and excluding the AP, of 2229 Gt y −1 , with an interannual variability σ = 109 Gt y −1 . We evaluated the model with various observational data, including in situ SMB observations , radarderived snow accumulation (Medley et al., 2013), CloudSat-CALIPSO satellite observations (Van Tricht et al., 2016b), GRACE mass changes, horizontal drifting snow transport fluxes , QSCAT-satellite-derived melt fluxes (Trusel et al., 2013), AP glacial discharge estimates (Wuite et al., 2015;Hogg et al., 2017) and SEB data from AWSs. We find a significant improvement in simulated SEB and SMB over the ice sheet, in particular at lower elevations. The largest improvement is found in modelled surface snowmelt, which now compares considerably better with the QSCAT and Neumayer (indirect) melt observations. No significant changes are found for the 5.5 km model version: here, model results are comparable to earlier versions.
This study shows that the latest version of RACMO2 can be used for high-resolution future projections of AIS SMB and SEB. However, limitations remain, mostly related to the cloud microphysics, the horizontal resolution and partly to the meltwater percolation scheme in the snow model. -GRACE mass anomalies (this study). Contact: b.wouters@uu.nl.
All other data used in this study are freely available by contacting the corresponding authors.
Author contributions. JMW, BPYN, WJB and MB conceived this study, decided on the new model setting and performed the analysis and synthesis of the data sets. JMW performed the model simulations and led the writing of the manuscript. CA, GB, CLJ, KK SL, BM CHR, KT, LDT, BW and JW processed and provided observational data sets. JMW, BPYN, WJB, MB, JTML, SRML, CLJ, CHR, EM and LU contributed to the development of the model. All authors contributed to discussions on writing the manuscript. financial and logistical support of the French Polar Institute IPEV (programme CALVA-1013) Graphics and calculations were made using the NCAR Command Language (version 6.3.0, 2017).
Edited by: Xavier Fettweis Reviewed by: Michael Lehning and one anonymous referee