Downscaled surface mass balance in Antarctica: impacts of subsurface processes and large-scale atmospheric circulation

Antarctic surface mass balance (SMB) is largely determined by precipitation over the continent and subject to regional climate variability related to the Southern Annular Mode (SAM) and other climatic drivers at the large scale. Locally however, firn and snow pack processes are important in determining SMB and the total mass balance of Antarctica and global sea level. Here, we examine factors that influence Antarctic SMB and attempt to reconcile the outcome with estimates for total mass balance determined from the GRACE satellites. This is done by having the regional climate model HIRHAM5 forcing 5 two versions of an offline subsurface model, to estimate Antarctic ice sheet (AIS) SMB from 1980 to 2017. The Lagrangian subsurface model estimates AIS :::::::: Antarctic SMB of 2473.5±114.4 Gt per year ::: yr-1, while the Eulerian subsurface model variant results in slightly higher modelled SMB of 2564.8±113.7 Gt per year ::: yr-1. The majority of this difference in modelled SMB is due to melt and refreezing over ice shelves and demonstrates the importance of firn modelling in areas with substantial melt. Both the Eulerian and the Lagrangian SMB estimates are within uncertainty ranges of each other and within the range of 10 other SMB studies. However, the Lagrangian version has better statistics when modelling the densities. There is a mean bias in modelled density of -24.0±18.4 kg m-3 and -8.2±15.3 kg m-3 for layers less than 550 kg m-3 for the Eulerian and Lagrangian framework, respectively. For layers with a density above 550 kg m-3 the bias is -31.7±23.4 kg m-3 and -35.0±23.7 kg m-3 for the Eulerian and Lagrangian framework, respectively. The mean firn 10 m temperature bias is 0.42-0.52 C. Further, analysis of the relationship between SMB in individual drainage basins and the SAM, is carried out using a bootstrapping approach. This 15 shows a robust relationship between SAM and SMB in half of the basins (13 out of 27). In general, when SAM is positive there is a lower SMB over the Plateau and a higher SMB on the westerly side of the Antarctic Peninsula, and vice versa when the SAM is negative. Finally, we compare the modelled SMB to GRACE data by subtracting the solid ice discharge, and find that there is a good agreement in East Antarctica, but large disagreements over the Antarctic Peninsula.There is a large difference between published estimates of discharge that make it challenging to use mass reconciliation in evaluating SMB models on the 20


20
The Antarctic Ice Sheet (AIS) has the potential to raise global sea level by 58 m (Fretwell et al., 2013) and it is therefore of utmost importance to understand its role in present sea level change in order to project it into the future. At present the AIS contributes 0.3±0.16 mm yr -1 to sea-level rise based on the average ice mass loss of 109±56 Gt yr -1 between 1992 and 2017 (Shepherd et al., 2018). An accelerating mass loss has been observed in West Antarctica and over the Antarctic Peninsula (AP), in the last four decades (Forsberg et al., 2017;Rignot et al., 2019). In the light of this acceleration, climatic changes are of 25 particular interest due to their role in inducing ice-sheet dynamic instability, by changing the mass influx to the ice sheet. The ice sheet mass balance (MB) can be split in an atmospheric and ice dynamic component: where D is the solid ice discharge in the form of iceberg calving, SMB is the surface mass balance composed of precipitation (P, snowfall and rain), sublimation and evaporation (S) from the surface, runoff (RO) of meltwater and erosion of blowing 30 snow. However, blowing snow is not taken into consideration in this study, so the SMB is defined here as: SMB=P-S-RO. Of these components, precipitation is by far the largest contributor (Krinner et al., 2007), and consists primarily of snow at higher altitudes. Melt and runoff of surface melt is largely confined to ice shelves and elevations less than 1400 m above mean sea level (Bell et al., 2018). Sublimation and evaporation are however important across most of the continent due to low humidity and high wind speeds (Palm et al., 2017). If SMB < D the total mass balance is negative and the ice sheet loses mass, and thereby 35 contributes to global sea level rise. Here we focus on the SMB component of the mass balance, to pin-point the immediate forcing to ice sheet dynamic instability. To estimate the SMB we use an atmospheric Regional Climate Model (RCM) to force a subsurface model, which outputs the SMB.
Regional climate models are most often used to downscale coarser global models and reanalysis because they add further detail, due to their higher resolution, e.g in the mountainous areas where the climate can be affected from local orography 40 creating katabatic winds or orographic forced precipitation (Rummukainen, 2010;Feser et al., 2011;Rummukainen, 2016).
Furthermore, RCMs also improve the physical representations of specific processes over polar areas . Mottram et al. (2020) evaluated Antarctic SMB calculated from the outputs from five different RCM simulations driven by ERA-Interim (1987. These five models showed mean annual SMB ranging from 1961±70 Gt yr -1 to 2519 ±118 Gt yr -1 . In the literature, individual evaluations of different RCMs such as COSMO-CLM 2 (Souverijns et al., 2019), MAR v3.6.4: 45  and RACMO2.3p2 ( Van Wessem et al., 2018) are found to be in the same SMB range. The overall model spread in SMB models corresponds to approximately 2 mm of sea level change yr -1 . Mottram et al. (2020) also showed that when compared to in-situ observation from both automatic weather stations and glaciological stake measurements, the data availability proved insufficient to distinguish between better performing model estimates. Fettweis et al. (2020) found similar conclusions for Greenland, where the RCMs displayed different strengths and weaknesses when evaluated both spatially and 50 temporally. Mottram et al. (2020) and Verjans et al. (2021) furthermore showed, that subsurface processes that drive melt and refreezing are extremely important when estimating the SMB. Hence, we here include firn processes by forcing a newly developed full-subsurface SMB model for Antarctica with the RCM HIRHAM5 (Christensen et al., 2007(Christensen et al., ) over 1979(Christensen et al., -2017 assess the effects of firn processes on estimates of ice sheet SMB. This subsurface model accounts for the physical properties of the uppermost part of the AIS, including density and temperature and the SMB.

55
Acknowledging that it might be challenging to judge the performance of the SMB model against in-situ observations , we also compare our modelled SMB results with a GRACE gravimetry estimate of the mass balance to determine any systematic biases. Finally, studies have shown that precipitation is not only the largest contributor to Antarctic SMB (Krinner et al., 2007;Agosta et al., 2019) it also has a spatial heterogeneous distribution varying over time, which affects the SMB (Fyke et al., 2017). Regional scale events like the heavy snowfall in Dronning Maud Land have an important 60 measurable effect on Antarctic SMB (Lenaerts et al., 2013;Turner et al., 2019). Different representations of these may explain differences between modelled SMB e.g.  as well as discrepancies between the GRACE mass balance and SMB-D solutions. Our study therefore also quantifies how regional climate indices affect SMB on basin scale.
Regional circulation patterns including ENSO (El Nino Southern Oscillation), the BAM (Baroclinic Annular Mode) and the Pacific-South American patterns (PSA1 and PSA2) have previously been identified as important determinants on weather and 65 climate variability in Antarctica (Turner, 2004;Irving and Simmonds, 2016;Marshall and Thompson, 2016). However, empirical orthogonal functional analysis of southern hemisphere 500 hPa geopotential height (Marshall et al., 2017), demonstrates that the southern annular mode (SAM) is the most important of these regional circulation indices. Further, Kim et al. (2020) found a multi-decadal relationship between the SAM and variations in the SMB, for these reasons we concentrate on its effects in this study. The SAM is an atmospheric phenomenon found across the extratropical southern hemisphere that influences the 70 climate over and around Antarctica (Fogt and Marshall, 2020). Marshall et al. (2017) found that the phase of the SAM, which describes pressure anomalies and precipitation in the southern hemisphere (Fogt and Bromwich, 2006), strongly affects the precipitation pattern over the AIS. Studies have shown that the phase of SAM, can have great impact on the surface climate in Antarctica, such as the temperature (Thompson and Solomon, 2002;Van Lipzig et al., 2008), sea-ice extent (Hall and Visbeck, 2002), pressure  and especially precipitation (Van Den Broeke and 75 Medley and Thomas, 2019). Other studies (Marshall et al., 2017;Dalaiden et al., 2020) have found that a positive SAM reduces precipitation over the Antarctic plateau and increases it over the western AP and in some coastal areas in East Antarctica.
Finally Vannitsem et al. (2019) found that the Antarctic SMB is influenced by the SAM in most of the coastal areas of East Antarctica and large parts of West Antarctica. Therefore, we also investigate the spatial distribution of SMB over the grounded AIS (GAIS) in relation to the phase of the SAM.

80
The aims of this study are thus; to estimate present-day Antarctic SMB using our subsurface model forced with the RCM HIRHAM5, compare and evaluate two subsurface model versions against each other and in-situ data. Furthermore, we estimate the MB, using our modelled SMB results combined with discharge values, and compare it with GRACE. Finally, we investigate the relationship between the SAM and the SMB. This is done in the following structure; first, the methods are presented, where the RCM HIRHAM5, the two subsurface models and their set-up are described. This is followed by the results, where the 85 modelled SMB results are shown, including evaluation against in-situ measurements of SMB, firn temperature and density.
Finally, the MB is estimated and evaluated against GRACE data, and we discuss the influence of SAM on SMB, followed by the conclusions. The HIRHAM5 RCM is a hydrostatic model with 31 atmospheric layers, developed from the physics scheme of the ECHAM5 global climate model (Roeckner et al., 2003) and the numerical weather forecast model HIRLAM7 (Eerola, 2006). HIRHAM5 has been optimised to model ice sheet surface processes that often are neglected or simplified in global circulation models. For a full description we refer to Christensen et al. (2007) and Lucas-Picher et al. (2012). Here HIRHAM5 is forced at the lateral boundaries at 6-hourly intervals with relative humidity, temperature, wind vectors and pressure from the ERA-interim reanal-95 ysis (Dee et al., 2011), further, daily values for sea ice concentration and sea surface temperature are also used. HIRHAM5 calculates the full surface energy balance at the surface, based on model physics as described in Lucas-Picher et al. (2012); Langen et al. (2015); Mottram et al. (2017). HIRHAM5 also calculates the amount of snowfall, rainfall, water vapor deposition and snow sublimation that occurs at the surface. Finally, for the HIRHAM5 Antarctic simulations, we used the Antarctic domain defined in the Coordinated Regional Climate Downscaling Experiment (CORDEX) (Christensen et al., 2014) and downscaled 100 it further to 0.11 • (≈12.5 km) spatial resolution with a dynamical time step of 90 seconds.

Subsurface model:
The subsurface model was originally built on ECHAM5 physics (Roeckner et al., 2003) but has been updated to include a sophisticated albedo scheme. Following Langen et al. (2015) the shortwave albedo is computed internally and uses a linear ramping of snow albedo between 0.85 below -5 • C and 0.65 at 0 • C for the upper-level temperature. The albedo of bare ice is 105 constant at 0.4. Furthermore, a transition albedo is calculated for thin snow layers on ice, based on Oerlemans and Knap (1998) with an e-folding depth of 3.2 cm for snow. Moreover, the snow and ice scheme is further developed and thereby updates the subsurface snow layers with snowfall, melt, retention of liquid water, refreezing, runoff, sublimation and rain (Langen et al., 2015. Thereby, the subsurface model is forced with the snowfall, rainfall, evaporation, sublimation and surface energy fluxes from HIRHAM5. These include net latent and sensible heat fluxes and downwelling shortwave and longwave radiative 110 fluxes for 6-hourly intervals over the period 1979-2017. To reduce RCM spin-up effects, such as misrepresentation of the physical state of the atmosphere e.g temperature, the first year is removed from the results. Furthermore, the model has been tuned to mimic the average behavior of the ice sheet surface at a 5-12 km scale. It cannot resolve subpixel processes. However, the small-scale features caused by surface melt translate into an increase of water content in the model. The subsurface scheme is updated hourly by interpolating the 6 hourly forcing files to 1 hourly time steps. To ensure a smooth transition between two 115 6 hourly files, a linear interpolation in time between the two nearest 6 hourly files is used. The horizontal resolution of the subsurface model follows the 0.11 • native resolution of HIRHAM5. As the Antarctic SMB may be sensitive to the subsurface model setup we here use two versions of the subsurface model . Common for both model versions is the albedo scheme, their meltwater percolation, firn compaction and heat diffusion schemes. Meltwater in excess of the irreducible water content (Coléou and Lesaffre, 1998) is transferred 120 vertically from one layer to the next using a parameterization of Darcy flow developed by Hirashima et al. (2010), with hydraulic conductivity values calculated from Van Genuchten (1980); Calonne et al. (2012) and coefficients from Hirashima et al. (2010). The impact of ice content on a layer's conductivity is described by the parameterization by Colbeck (1975). When meltwater can infiltrate into a subfreezing layer, it is refrozen and latent heat is released. Firn density is updated at each time step for compaction under each layer's overburden pressure using the parameterization by Vionnet et al. (2012).

125
The two model versions differ in the management of the layers within the model. The first model version developed by Langen et al. (2017), has 32 subsurface layers with a fixed predefined mass, expressed in m water equivalent (w.eq.), given by D N = D 1 λ N −1 , where N is the given layer and D 1 = 0.065 m w.eq. This fixed model implies an Eulerian framework, meaning that when snowfall occurs at the surface, it is added to the first layer and an equal mass from that layer is shifted to the underlying layer. The same goes for each layer in the model column. The same procedure is followed when mass is removed 130 from the top layer due to runoff or sublimation. Then each layer takes from their underlying neighbour an amount of snow/firn equivalent to the mass lost at the surface. The temperature and density of the layers are updated as the average between the snow or firn that is received by the layer, and what remains there. In the following we refer to this model version as the Fixed model.
The second model version uses a Lagrangian framework for the layer evolution developed by Vandecrux et al. (2018Vandecrux et al. ( , 2020a b). Layers evolve through splitting and merging dynamic based on a number of weighted criteria. This dynamical model, henceforth referred to as Dyn model, has 64 subsurface layers, the number of which are fixed during the simulation. When snowfall occurs at the surface, it is first stored in a "fresh snow bucket". When this snow bucket reaches 0.065 m w.eq., its content is added as a new layer at the surface of the subsurface scheme and two layers need to be merged elsewhere in the model column. The layer merging scheme assesses how likely a layer is to be merged with its underlying neighbour based 140 on seven criteria: the layers' difference in temperature, density, grain size, water content, ice content, depth and the thickness of the layers. The first five criteria makes it preferable to merge layers with small differences, the sixth criterion, makes it preferable to merge deep layers rather than shallow layers, in this case the shallow layer limit is set to 5 m w.eq., this criterion carries twice the weight of the first five. The final criterion says that no layer can be thicker than a maximum thickness, in this case 10 m w.eq., this is set to avoid the deepest layers continuing to grow. A weighted average of the criteria, where the 145 first five are weighted equally, while the depth and thickness criteria are weighted double and triple respectively, is used by the model to determine which layers should be merged. When surface sublimation or runoff occur it is taken from the snow bucket and then from top layer. When a layer decreases in thickness and its mass reaches 0.065 m w.eq., then it is merged with the underlying layer and another layer can be split in two elsewhere in the model column. The splitting routine is based on two criteria; thickness of the layer where thick layers are more likely to split; and shallowness where shallow layers are more likely 150 to split. The two criteria are weighted 60/40. However, the minimum thickness of any layer is always 0.065 m w.eq. to avoid numerical instability. The bottom of the lowest model layer is assumed to exchange mass and energy with an infinite layer of ice with a temperature, like in the Fixed model Langen et al. (2015), calculated from climatological mean of the HIRHAM5 2 m temperature.
Another difference between the two model versions is that the dynamic-layer model melts simultaneously the snow and ice 155 content of the top layer while the Fixed-layer model melts the snow content first then the ice content of the top layer. This Melt First snow then ice Snow and ice simultaneously Snow and ice simultaneously update aims at preventing the top layer from becoming only ice and a barrier to meltwater infiltration. Furthermore, the Dyn models runoff is routed downstream using Darcy's law and the local surface slope, whereas the Fixed model follows Zuo and Oerlemans (1996) and excess water in a layer cannot be transferred to the underlying neighbour. Both the Fixed and Dyn versions require a fresh snow density value when adding snowfall at the surface. We here use the Antarctic parameterization 160 from Kaspers et al. (2004) which use local climatological means of skin temperature, 10 m wind speed and accumulation rates, here the means from HIRHAM5 has been used.

Experimental Set-up
The Fixed model was initialized with a firn column with uniform density of 330 kg m -3 , and a temperature at the bottom of the firn pack given by the climatological mean of the HIRHAM5 2 m temperature. Spin-up was performed by repeating a All three model simulations (summarized in Table 1) provide outputs of monthly and yearly means of all 3D variables 175 (density, grain size, firn temperature and ice/water/firn content), daily 2D fields (SMB, runoff, super imposed ice, melt, albedo, ground heat flux, refreezing, diagnosed snow depth (which is an estimate based on the snow concentration in each layer), net shortwave and net longwave radiation) of the surface variables, furthermore daily columns for specified coordinates interpolated to the nearest grid cell, have been retrieved, for comparison of in-situ measurements. For the two simulations with dynamical layer thickness, the daily 3D fields are interpolated into a fixed grid, with the same number of layers, so time averages could 180 be calculated.
6 2.4 Regional drivers and mass balance The SAM is characterized in Fogt and Bromwich (2006) as the zonal pressure anomalies in the high southern latitudes having opposite sign to those of the mid latitudes. The SAM drives the westerly winds around Antarctica but the stream oscillates north-south. The SAM can have three phases: positive, neutral or negative, where positive creates a higher pressure over the 185 mid latitudes and lower pressure over Antarctica, and thus moves the westerly winds closer to Antarctica. A negative SAM creates a lower pressure over the mid latitudes and a higher pressure over Antarctica, moving the westerly winds north. When neutral there is no pressure difference anomaly. To investigate how the phase of SAM affects the SMB, monthly SAM data, as calculated by Marshall (2018) Observing the mass balance can be helpful to assess the spatial patterns of SMB and evaluate the modelled results. Mass  3 Results In the model mean  of the three SMB simulations (Fig. 1a), we see that the majority of the total AIS (ToAIS) has a positive SMB, only a few regions show a negative SMB: Larsen ice shelf, George IV ice shelf, coastal regions of Queen Maud Land, the Transantarctic mountains, near Amery ice shelf and some coastal areas in East Antarctica. Near Vostock in East Antarctica, the SMB is less than 25 mm w.eq. yr -1 . The SMB increases towards the coast due to higher precipitation. The 215 highest SMB is greater than 2000 mm w.eq. yr -1 and is found on the windward (western) side of the AP, whereas the most negative SMB, -500 mm w.eq. yr -1 , is found on the leeward (eastern) side of the AP (Fig. 1a). All the model simulations show nearly identical SMB values over the GAIS, however they differ the most near the coast in West Antarctica and the AP as Fig.   1b shows. Here, we see that the δSMB (model minus mean) show that the Fixed version has a higher SMB of up to 550 mm w.eq. over the Larsen ice shelf relative to the model mean, in Dyn03 the SMB values differ between -350 and 400 mm w.eq. 220 from the model mean, this change occurs over a few grid cells. In Dyn15 the SMB differ up to -650 mm w.eq. compared to the model mean over the Larsen ice shelf. Since the Fixed is above the model mean, over the Larsen ice shelf, and Dyn15 is below the model mean, it looks like that the rapid change from negative to positive δSMB in Dyn03 over Larsen ice shelf is due to lack of spin-up. Below the AP, of the coast of Ellsworth Land and Marie Byrd Land, the Fixed version models a lower (-75 mm w.eq.) SMB than Dyn03 (35 mm w.eq.) and Dyn15 (50 mm w.eq.) all relative to the model mean. Around Alexander

225
Island in the Bellingshausen sea, both the Fixed and Dyn15 have a lower SMB compared to Dyn03. The differences in spatial distribution show that in areas where melt occurs, the SMB is very sensitive to which subsurface scheme is used.  The model differences are seen in the integrated values for precipitation, SMB, melt, refreezing and runoff, for both ToAIS and the GAIS (Fig. 2), and summarized in Table 2. As all model simulations are forced using the same precipitation field ( Fig. 2a) and since the precipitation is the main driver of the SMB, the variability of the modelled SMB closely follows the 230 precipitation variability. The spread in modelled mean melt, refreezing and runoff are respectively 1%, 11% and 8% smaller when including the ice shelves compared to only taking the GAIS, whereas the spread in mean SMB becomes 3% greater. To better compare the melt, refreezing and runoff from the different simulations, the fraction of runoff to melt is shown in figure   2f. Dyn03 has the smallest runoff fraction whereas Dyn15 and Fixed are quite close to each other. This implies, that even

Evaluation against observations
Koenig and Montgomery (2019)  Modelled firn densities are evaluated using the SumUp dataset (Koenig and Montgomery, 2019). When disregarding firn cores shallower than 2 metres, there were 139 density profiles left (Fig. 4), all the references for the firn profiles can be found 255 in the reference list. These profiles vary in depth, from a few metres to 100 metres, but the majority are drilled to 10 metres depths. Knowing the coring date, we compare it to the modelled density of the nearest grid cell the same date. Before the inter-comparison, the modelled and observed density profiles were interpolated to the same vertical resolution (if the model resolution is higher than the core resolution, the model is interpolated to fit the core resolution and vice versa). In the SumUp dataset 96 profiles had the exact date given, 7 SumUp profiles only had year and month given, here the modelled mean density 260 of the given month were compared, finally 36 cores had only the year given, in these cases the modelled mean density of January were compared, as we assume they were most likely collected in the middle of the standard Antarctic summer field season. To evaluate the model performance we calculate Mean Difference (MD) and Standard Deviation (SD) between the modelled and observed firn densities. A statistical comparison of the mean difference and one standard deviation between the firn cores and the modelled densities 265 are given in Table 3 for the three simulations. Summed up over the AIS all simulations overestimate the densities below 550 kg m -3 and underestimate the densities above 550 kg m -3 . It is seen that the Fixed version outperforms the Dyn03 and Dyn15 for densities below than 550 kg m -3 . Whereas, Dyn03 and Dyn15 outperforms the Fixed version for densities above 550 kg m -3 .
All three simulation show the best statistics for higher densities, the agreement with the in-situ cores also vary spatially (Fig.   5). Generally the spatial density bias is consistent between the models.
270 Table 3. Mean difference between the modelled and observed firn densities (model -core) and standard deviation of the modelled densities above and below 550 kg m -3 . In total 139 cores were used, see Fig. 4   (2020) show that the HIRHAM5 model estimates higher precipitation over the Filchner-Ronne ice shelf than other RCMs, and the overestimate in density may therefore relate to overestimated precipitation in this area, which is compliant with our Fig.   3. However, as they also note, the lack of continuous SMB observations makes it difficult to be certain if and by how much precipitation is overestimated in this region. It could also be due to an overestimation for melt and refreezing over the ice shelf. For the Ross ice shelf cores and near the South Pole, the Fixed simulation overestimates most of the cores, some of them by 50 to 100 kg m -3 for densities less than 550 kg m -3 , and more than 100 kg m -3 for densities greater than 550 kg m -3 . However, for Dyn03 and Dyn15 we also observe an overestimation of most cores, but only a six them are overestimated by more than 25 290 kg m -3 . 14 Figure 6 shows four of the 139 firn cores: core BER02C90_02 (Wagenbach et al., 1994b) (Fig. 6a), core DML03C98_09 (Oerter et al., 2000a) (Fig. 6b), core FRI14C90_336 (Graf and Oerter, 2006h) (Fig. 6c) and core Site 11 (Morris et al., 2017) ( Fig. 6d). These four cores are selected because they are located in different regions of the AIS and, furthermore, they show different examples of under/over estimations of modelled densities. The Fixed simulation fit quite well (±20 kg m -3 ) with the 295 core taken on Berkner Island (Fig. 6a), whereas Dyn03 and Dyn15 show a lager bias mainly at the surface and the top three meters of the firnpack. The core from Dronning Maud Land (Fig. 6b) has a high vertical resolution, the deeper the cores goes, the smaller the biases become. Cores FRI14C90_336 and Site 11 are taken on the Ronne ice shelf and in Marie Byrd Land respectively. The model densities in FRI14C90_336 are overestimated below one meter depth, the mean bias is 200 kg m -3 . At Site 11 all simulations underestimate the density, however below two meters depth, the underestimation is nearly constant with 300 a mean bias of -20 kg m -3 .  (Wagenbach et al., 1994b), 6b DML03C98_09 taken in 1998 (Oerter et al., 2000a), 6c FRI14C90_336 taken in 1990 (Graf and Oerter, 2006h) and 6d Site 11 taken in 2013 (Morris et al., 2017).
The modelled subsurface temperatures are evaluated against observed 10 metre firn temperature measurements from 49 boreholes (van den Broeke, 2008) (see Fig. 4 for the locations). Most of the temperatures were taken in the 1980's and 1990's, however only the year or decade is known for when these were taken. Therefore they are compared with the modelled mean 10 m firn temperature from 1980-2000. We evaluated the model performance using the Root Mean Square Difference (RMSD),

305
Mean Difference (MD) and coefficient of determination (R 2 ). Subsurface temperatures are only sparsely available in Antarctica. The measured 10 m firn temperatures are compared with the modelled mean 10 m firn temperature of the nearest grid cell (Fig.   7). The red, blue and green lines are the regression lines of 1 st order, for the Fixed, Dyn03 and Dyn15; they have an R 2 of 0.98, 0.97, 0.98, respectively. It is assumed that the in-situ temperatures are true, so the errors are in the modelled temperatures.
For temperatures below -30 C • the three simulations are in agreement, but in warmer firn temperatures >-30C • , the agreement 310 becomes worse. The mean deviation of the three model simulations are listed Table 4.

Discussion
The annual SMB for the three simulations (Table 2) are of the same magnitude as the previous HIRHAM5 SMB estimate of 2659 Gt yr -1 for the ToAIS . However, we model a lower SMB, with only Fixed and Dyn03 within one standard deviation range of Mottram et al. (2020). The lower SMB estimates are due to the inclusion of the runoff component in 315 the SMB calculation. The initial SMB result from HIRHAM5 in Mottram et al. (2020) were only calculated from precipitation, evaporation and sublimation. Calculating the SMB by including a subsurface model results in a more realistic SMB, due to the fact that it takes surface and subsurface processes like energy fluxes, meltwater percolation and refreezing into account.
The spatial distribution of SMB fits reasonably well compared with the SumUp accumulation measurements, however, more measurements, especially in East Antarctica, are needed to be able to do a complete evaluation. Furthermore, the spatial the geographical distribution of precipitation is uneven between these models, with COSMO-CLM 2 much drier in western Antarctica than other models in the comparison. Even using a common ice mask, Mottram et al. (2020) found that the difference in precipitation is around 500 Gt yr -1 between HIRHAM5 (the wettest model) and COSMO-CLM 2 (the driest model in the intercomparison). The high precipitation in regions of high relief in HIRHAM5 is attributed to a wet bias in the precipitation 330 scheme, also identified in southern Greenland and similarly occurring in the RACMO2.3p2 regional climate model (Hermann et al., 2018). In both models this wet bias in steep topography is related to the precipitation and cloud micro-physics schemes . Areas with a negative SMB can be due to large melt rates, which is what we see in the model over the Larsen ice shelf with melt values between 1200 and increasing toward west to 2300 mm w.eq. yr -1 and SMB values in the range of 300 to 1800 mm w.eq. yr -1 increasing toward west. In general all three simulations display a higher melt compared to 335 other RCM studies e.g 71 Gt yr -1 in RACMO2.3p2 (van Wessem et al., 2018) or 40 Gt yr -1 in MARv3.6.4 .
These two numbers are without the AP, but are nevertheless very low compared to our melt rates. Trusel et al. (2013) derived satellite-based melt rate estimates from 1999 to 2009 and over that period, the Larsen ice shelf experienced the largest melt of around 400 mm w.eq. yr -1 . However, these estimates were derived using RACMO2.1, and the satellite detects melt areas on the Larsen ice shelf that were not simulated in RACMO2.1, most likely due to coarse resolution, so the 400 mm w.eq. yr-1 340 might be on the low end. Nevertheless, Trusel et al. (2013) estimates are still three to six times lower than our simulation. This suggests that the subsurface model may compute a melt rate that is too high in at least some locations.
Negative SMB values can also be due to high sublimation rates in e.g., blue ice areas (Hui et al., 2014). For exampe, Kingslake et al. (2017) found blue ice in Dronning Maud Land and near the Transantarctic Mountains, in these areas our SMB model mean also shows negative SMB between -50 and -400 mm w.eq. yr -1 . A closer investigation (not shown) reveals that the 345 negative SMB values in these areas are driven by the sublimation and thereby consistent with the creation of blue ice areas.
The differences in SMB between the model simulations ( Fig. 1b) are largest near the coast in West AIS and especially on the Larsen ice shelf. This is confirmed in Fig. 2b, where the difference in integrated SMB between the model simulations are greater when the ice shelves are included. We attribute the differences between the Fixed and Dyn models to the following differences in model designs: The increased vertical resolution in the Dyn models, with a higher vertical resolution (the top 350 layers can be 6.5 centimeter w.eq. thick) means that the cold content in the upper layers is depleted faster and it starts to melt while the layer below is potentially still below freezing. Whereas the top layers in the Fixed model get thicker rather quickly, which means it takes longer to be brought to melting point and it starts melting. Furthermore, the two versions of the subsurface model, have different melting schemes. In both versions one layer can contain snow and ice at the same time, described with a fraction, However, in the Fixed model melts snow first and then, if there is more energy left, melts the ice. Whereas, the Dyn 355 melts snow and ice simultaneously. This simultaneous melting of snow and ice was introduced in the Dyn version to prevent the top layer being depleted of its snow content and left only with ice (Vandecrux et al., 2018). A top layer composed of ice would then prevent surface melt to infiltrate below the top layer. By melting snow and ice simultaneously, there is always snow in the top layer for meltwater infiltration to happen. This difference of infiltration may cause the snow-melt first to refreeze less and runoff more water than the snow-and-ice-simultaneous melt. To investigate these differences in melt, refreezing and 360 runoff, the runoff fractions have been plotted in figure 2f and listed in table 2. Here it is seen that even though the difference in melt between Dyn03 and Dyn15 is only around 20 Gt yr -1 , the difference between the runoff to melt fractions are larger. The Fixed model melts around 300 Gt yr -1 less than the dynamical versions but the runoff to melt ratio is the same for the Fixed and Dyn15. This means that the Fixed and Dyn15 have the same relative runoff, leading to the same relative refreezing, indicating that this difference does not cause significant partition of melt between refreezing and runoff.

365
The difference in SMB between the three simulations confirms how complex it is to estimate the SMB. Just by changing the subsurface scheme the final result differs by 90 Gt yr -1 . By keeping the same subsurface scheme and changing the spin-up length the final result differs by 110 Gt yr -1 . These changes in SMB illustrate the consequences of including dynamic firn processes since the layer density and temperature and other firn properties are better conserved, potentially allowing more retention and refreezing where there is capacity or reducing it where there is not. Although these differences are currently only 370 a few percent of the total SMB, as the climate warms and melt becomes more widespread in Antarctica e.g. Kittel et al., 2020), accounting for these processes will become more important. Moreover, on a local and regional scale, the differences are more important when determining mass balance in basins or outlet glacier/ice shelves.
The differences between versions with a different spin-up period suggest that the snowpack is not quite in equilibrium in all locations. Therefore, SMB calculations consequently vary due to the amount of melt calculated during the initialisation period.

375
Retention and refreezing of meltwater during spin-up causes different profiles of temperature and density to develop depending on how long the spin-up lasts for. These results therefore emphasise the importance of adequate spin-up and assessment of the effects of snowpack spin-up in producing and using SMB in Antarctica. Vandecrux et al. (2020b) found that the Fixed version smoothes the firn density profiles, when compared to the dynamical version, this is confirmed by our results. One of the criteria for the dynamical version is that it prefers to merge layers deeper 380 than 5 m of w.eq., meaning that the top 5 m w.eq. have a high vertical resolution, this makes it easier to detect changes in density. In areas such as the AP, Ronne-Filchner ice shelf, Ross ice shelf and in coastal areas of Dronning Maud Land where seasonal melt occurs (Zwally and Fiegles, 1994;Wille et al., 2019), meltwater can percolate into the firn and refreeze, creating ice lenses that change the density, but that cannot be detected if the subsurface scheme have layers with a fixed mass even if the vertical resolution is increased (Vandecrux et al., 2020b). Not only is there a difference between the models when evaluating 385 density profiles, this study also shows the importance of spatial evaluation, here the three simulations follow the same pattern by over/underestimating the densities in the same areas (Fig. 5). This systematic bias may indicate either further tuning of densification routines are necessary or that there are systematic biases in accumulation, leading to these errors. The subsurface scheme does not currently incorporate wind-blowing snow processes that may prove important also in correcting biases in accumulation. On the other hand, although 0.11 degree is a high resolution model in Antarctica and thus better captures 390 topographic variability than lower resolution models, it is still relatively coarse when it comes to capturing steep topography.
Errors in orographic precipitation are difficult to measure even in well instrumented basins and are poorly captured in Antarctica where observations are few and far between. The densification bias becomes especially important when using altimetry data to estimate the total MB, like in Shepherd et al. (2018) and Rignot et al. (2019), here the firn densification rate is needed to correct the altimetry data (Griggs and Bamber, 2011).

395
Since the density cores are primarily taken from West Antarctica and Dronning Maud land, these statistics represents complex areas with high precipitation and melt/refreezing events, whereas, density comparisons from less complex areas (low precipitation and no melt/refreezing) such as East Antarctica are sparse. Nonetheless they are still very important. Based on the statistics from these model set ups, the Dyn version is preferred when modelling densities above 550 kg m -3 .
The simulated 10 m firn temperature depends on the thickness and numbers of the layers above the 10 m point. The thickness 400 of a layer determine how conductive heat fluxes are resolved in the near-surface snow. A thicker layer will have more thermal inertia and will require more energy to be warmed up. A thin layer can respond much quicker to fluctuations in the surface energy balance. Differences in simulated temperatures between models, as we see in Table 4 can therefore be explained by vertical resolution, which affects both their calculation of temperature and how the heat is conducted to a depth of 10 m. Note that the models also use different thermal conductivity parameterisations.  signal for precipitation, which confirm our results since precipitation is the main driver of SMB. Basin 26 shows a negative SMBa and basin 27 has a slightly positive SMBa, this is due to the steep orography on the windward side of the AP creating a shadowing effect on the leeward side of the AP. For SAM-the SMBa signal is opposite and the mean magnitude of the signal is 26% larger in all basins (Fig. 9a). During months of SAM+ the average SAM index is 1.45. Figure 9 show that basins 1-5, From 1980 to 2017 the SAM has become more positive (Fogt and Marshall, 2020), this positive trend in the SAM is attributed to stratospheric ozone depletion (Thompson et al., 2011;Fogt et al., 2017). If this trend continues the basins on the leeward side of the AP will see a smaller mass gain in the future, which could accelerate the collapse of the Larsen ice shelf. This is also seen in basin 9, 10, 11 and 12 surrounding the Amery ice shelf and basin 21 where Thwaites glacier is located. Not all of the above mentioned basins show δSMB signals that are statistically robust (i.e. the signals are within the 5 th or 95 th percentiles), 460 but if the trend in the positive SAM continues, it might become an important factor in the future (Fogt and Marshall, 2020).
It is thus important to take the SAM phases into account when investigating the SMB at regional scale, furthermore it is extremely important that global circulation models resolve the SAM realistically if future climate projections are to be used with confidence to make projections of sea level rise from Antarctica.  show the Lagrangian model set up has the lowest bias and standard deviation in density differences for densities grater than 550 kg m -3 , for densities less than 550 kg m -3 the Eulerian performance best. In general all models overestimate the densities 475 on the Flicher-Ronne ice shelf and underestimate the densities in Marie Byrd Land and around the Ross ice shelf. It is therefore clear that there are regional systematic biases. To evaluate our simulated SMB, we compare our simulations with the SumUp accumulation rates, half of the comparing sites fit with ±13%, moreover we also compare our simulations to MB estimations (SMB minus discharge) from GRACE. We use discharge from two sources: Gardner et al. (2018) and Rignot et al. (2019).
There are large differences between the discharge values over the AP leading to our simulations overestimating MB when 480 using D Gardner and underestimating MB when using D Rignot . Over the East GAIS the MB is underestimated using D Rignot , but fit quite well to the GRACE MB when using D Gardner . These disagreements between the two observational datasets makes it hard to distinguish how well the modelled SMB fits with total mass balance estimates.
Regional precipitation is strongly linked to the phase of the SAM as shown by the bootstrap analysis. By using outputs from HIRHAM5 forced with ERA-interim to resolve the SAM correctly, robust signals are identified in 13 out of 27 basins. It is 485 clear that the phase of the SAM affects the spatial distribution of SMB. When SAM is negative, there is a lower SMB on the windward side of the Antarctica Peninsula and a higher SMB over the plateau and vice versa when SAM is positive. This makes the SAM an important factor to evaluate in global models when downscaling models for projecting future Antarctic climate.
Code availability. A Matlab version of the subsurface model code used in these study is available here: