Meltwater runoff and glacier mass balance in the high Arctic: 1991-2022 simulations for Svalbard

. The Arctic is undergoing increased warming compared to the global mean, which has major implications for fresh-water runoff into the oceans from seasonal snow and glaciers. Here, we present high-resolution (2.5 km) simulations of glacier mass balance, runoff and snow conditions in Svalbard from 1991-2022, one of the fastest warming regions in the Arctic. The simulations are created using the CryoGrid community model forced by both CARRA reanalysis (1991-2021) and AROME-ARCTIC forecasts (2016-2022). Updates to the water percolation and runoff scheme are implemented in the CryoGrid model 5 for the simulations. In-situ observations available for Svalbard are used to carefully evaluate the quality of the simulations and model forcing. The overlap period of 2016-2021, when both CARRA and AROME-ARCTIC data are available, is used to evaluate the consistency between the two forcing datasets. We find a slightly negative climatic mass balance (cmb) over the simulation period of -0.08 m w.e. yr − 1 , but with no statistically significant trend. The average runoff was found to be 41 Gt yr − 1 , with an significant increasing trend of 6.3 Gt 10 decade − 1 . In addition, we find the simulated climatic mass balance and runoff using CARRA and AROME-ARCTIC forcing are similar, and differ by only 0.1 m w . e . in climatic mass balance and by 0.2 m w . e . in glacier runoff when averaged over all of Svalbard. There is, however, a clear difference over Nordenskiöldland, where AROME-ARCTIC simulates significantly higher mass balance and significantly lower runoff. This indicates that AROME-ARCTIC may provide high-quality predictions of the total mass balance of Svalbard


Introduction
Glaciers and ice caps are considered to be good indicators of climate change.During the last decades, glaciers and ice caps worldwide have been responding to a globally warming climate by melting at increasing rates (e.g.Vaughan et al., 2013;Huss and Hock, 2018).The Arctic has experienced greater warming than the global average due to positive feedbacks triggered by changing sea ice cover, the so-called Arctic amplification (e.g.Serreze and Francis, 2006;Graversen et al., 2008;Lind et al., 2018).As sea ice continues to retreat, further warming in the Arctic is expected (e.g.IPCC, 2019).
In particular, the region around the Barents Sea, which includes the archipelagos of Svalbard, Franz Josef Land and Novaya Zemlya, has experienced pronounced warming in recent decades due to disappearing sea ice (e.g.Screen and Simmonds, 2010;Lind et al., 2018).For example, the Svalbard archipelago has the strongest observed warming in Europe since the 1960's, with temperatures increasing at a rate of 0.5 • C decade −1 (Nordli et al., 2014).Even under the moderate RCP4.5 emission scenario, which projects a global temperature increase of 1.1 -2.6 • C by 2100 relative to the 1986-2005 period, temperatures in the Barents Sea region are projected to increase by 5-9 • C (AMAP, 2017).
Although the volumes of ice on Svalbard and the Russian Arctic are only equivalent to a global sea level rise of about 15 mm (Fürst et al., 2018) and 32 mm (Farinotti et al., 2019), respectively, they are estimated to be some of the most important regional contributors to sea level rise in the 21st century (e.g.Meier et al., 2007;Church et al., 2013;Hock et al., 2019).In addition to sea level rise, meltwater from retreating glaciers is important for river hydrology, fjord circulation, and terrestrial and marine ecosystems (e.g.Carroll et al., 2017;Hopwood et al., 2020).
While few glaciological measurements in the Russian Arctic are available, glaciological measurements of the surface mass balance (smb) have been conducted on Svalbard since the 1960's, and show an accelerating mass loss (e.g.Hagen et al., 2003;Schuler et al., 2020).However, these observations only exist in a small area, therefore, dedicated energy/mass balance models are an important tool to determine the mass balance of the whole archipelago.
The change in the climatic mass balance (cmb), which equals the smb plus refreezing in the firn column, has previously been estimated for Svalbard using energy balance models at different spatial and temporal resolutions (Lang et al., 2015;Aas et al., 2016;Østby et al., 2017;Van Pelt et al., 2019).Previous studies of the mass balance evolution primarily use energy balance models forced by either the ERA-Interim reanalysis product (Dee et al., 2011), the ERA-40 reanalysis product (Uppala et al., 2005), or ECMWF operational analysis, all of which are further downscaled to a 1-10 km resolution either using a regional climate model (e.g.Lang et al., 2015;Aas et al., 2016), statistical downscaling (e.g.Østby et al., 2017), or both (e.g.Van Pelt et al., 2019).Since the energy/mass balance models used in these studies are forced by similar global climate products, inter-model disagreements are due to uncertainties from the downscaling of the climate forcing, differences in physical parameterizations, or spatial resolution.The mass balance of the Russian Arctic has not previously been estimated using energy balance models, but has been simulated by various global temperature-index models (e.g.Radić et al., 2014;Hock et al., 2019).
In recent years, high-resolution simulations of the meteorological conditions over the Barents Sea region have become available.The high-resolution Copernicus Arctic Regional ReAnalysis (CARRA) dataset (Schyberg et al., 2020;Yang et al., 2021) https://doi.org/10.5194/egusphere-2022-1409Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.was published in 2021.It is a reanalysis product with a 2.5 km resolution, downscaled from ERA5 (Hersbach et al., 2020) by the state-of-the-art weather prediction model HARMONIE-AROME (Bengtsson et al., 2017).CARRA includes a number of improvements over ERA5, such as the assimilation of a large amount of additional surface observations, extensive use of satellite data, and an improved representation of sea ice.CARRA is likely the best high-resolution estimate of the meteorological parameters available in the Barents Sea region currently available due to the complex physics contained within the model and the large amount of assimilated data.Further downscaling is not required since it already contains such a high spatial resolution, which avoids introducing more uncertainties.In addition, the high-resolution AROME-ARCTIC numerical weather prediction system has provided forecasts over the Barents Sea region since late 2015 at 2.5 km resolution.AROME-ARCTIC is built upon HARMONIE-AROME, the same state-of-the-art weather prediction model as CARRA, but since it is a forecast product, it likely contains larger uncertainties in the estimates of the meteorological variables.
Here, we present a high-resolution (2.5 km) dataset of the glacial mass balance, runoff, refreezing, and seasonal snow amount of Svalbard from 1991-2022.The dataset is created using a full energy balance model forced by the CARRA reanalysis and AROME-ARCTIC forecasts.The produced dataset has the potential to provide near real-time forecasts of the mass balance and runoff in the region.We use Svalbard as a study region, as automatic weather stations and yearly mass balance surveys are available for evaluation of our model performance, while very few observations are available in the Russian Arctic.Using Holocene glacier records, Lubinski et al. (1999) found that glacier mass balance changes in the Russian Arctic closely mimic that of Svalbard, although the magnitude of mass change differs.We, therefore, expect a model that performs well for Svalbard will also provide a high-quality dataset for the Russian Arctic.
The mass balance simulations are conducted using CryoGrid, a model for simulating the terrestrial cryosphere (Westermann et al., 2022).CryoGrid simulates the energy and mass balance of both seasonal snow and glaciers, and estimates permafrost in non-glaciated areas.
For the first time, we use the novel, high-resolution CARRA reanalysis (Schyberg et al., 2020) as well as AROME-Arctic forecasts as forcing for mass balance simulations on Svalbard.Although other products based on HARMONIE-AROME have successfully been used as forcing for mass balance simulations in the Arctic (e.g.Mottram et al., 2017;Schmidt et al., 2018), neither AROME-ARCTIC nor CARRA have previously been validated for use in mass balance simulations.We therefore thoroughly evaluate both forcings against available observations, and compare the results from both simulations for the period they overlap.Lastly, the CryoGrid model results are evaluated against available observations of mass balance, both from in-situ campaigns and geodetic methods.

Study area
Located in the Norwegian Arctic between 75 and 81 • N, the Svalbard archipelago is currently in one of the fastest warming regions in the world (e.g.Nordli et al., 2014;Isaksen et al., 2016).With a land area of ∼60,000 km 2 , of which about 57% is covered by glaciers (Nuth et al., 2013), it contains about 10% of the glacier area in the Arctic, outside of the Greenland ice https://doi.org/10.5194/egusphere-2022-1409Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.sheet.The glacier types vary between small cirque and valley glaciers to large ice fields and ice caps, with more than 1000 individual mapped glaciers across the archipelago.Around 60% of the glacier area belong to tidewater glaciers (Błaszczyk et al., 2009) which introduce freshwater into the oceans through discharge from subglacial channels or calving at the glacier front.The highest elevations on Svalbard reach 1700 m a.s.l., but the hypsometry of glaciers peaks at 450 m a.s.l (Noël et al., 2020).
While the western side of the Archipelago is kept warm and humid by the Norwegian current, which brings warm Atlantic currents northwards along the western coast (Walczowski and Piechura, 2011), and warm and moist air from southerly flows, the eastern side is colder and drier, dominated by cold Arctic ocean current and dry and moist air masses originating in the north-east (Käsmacher and Schneider, 2011).Precipitation varies wildly across the archipelago, with the highest precipitation rates along in the south and along the west coast (Førland and Hanssen-Bauer, 2003;Winther et al., 2003;Førland et al., 2020).
These patterns in temperature and moisture are reflected in the distribution of glaciers, with the largest glaciers found in the north-east and less extensive glacier coverage along the western side of Svalbard and in central Spitsbergen.
3 Methods and Data

Methods
The simulations presented in this paper were created using the full energy balance model CryoGrid, which is forced by both CARRA reanalysis and AROME-ARCTIC forecasts.The workflow used is described below.
First, both the CARRA reanalysis and AROME-ARCTIC forecasts are evaluated against available observations from automatic weather stations.In addition, the consistency between the two forcings are evaluated for the overlap period (2016)(2017)(2018)(2019)(2020)(2021).
The results of this analysis is described in Supplement S2.
We then perform a 30-year spin-up of the CryoGrid model (described in section 4)for the glaciated gridpoints by repeating the 1991-2000 CARRA forcing to initialise the snow/ice temperature, density, and water content.Initially, the entire domain consists of temperate, pure glacier ice, i.e. the ice temperature of the entire column is 0 • C. Tests were conducted with lower initial temperatures (-5 • C), but it did not affect the temperature profile at the end of the spin-up.At the end of the spin-up period, the runoff, refreezing, subsurface temperatures and climatic mass balance reached stable values.For the non-glaciated land points, only a 2-year spin-up was used.
The energy/mass balance model CryoGrid is then used to simulate the mass balance components of both glaciers and seasonal snow from 1991-2021 using the CARRA reanalysis as forcing.
A second simulation with CryoGrid, this time forced by AROME-ARCTIC, is then conducted from 2016 to present.From 2016 until the summer of 2019, the AROME-ARCTIC model was initialized with too little snow over some glacier points in the ablation area, thus leading to unrealistically high surface and 2m temperatures.To counter this effect, we use the 10m temperature for the AROME-ARCTIC-forced simulation when unrealistically high surface temperatures occur.The AROME-ARCTIC-forced simulation is initialized from the CARRA-forced simulation at the end of 2015.Thus, the initial conditions https://doi.org/10.5194/egusphere-2022-1409Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.
for the 2016-2021 period is identical for the two simulations.Lastly, the output from both CryoGrid simulations is evaluated against in-situ mass balance observations and geodetic estimates.
The AROME-ARCTIC-forced CryoGrid simulation is automatically updated on a daily basis.In this study, we present the simulations spanning until October 1st, 2022.
For the CryoGrid simulations, a fractional glacier mask is created by computing the percentage of glacier coverage in each grid point.The glacier coverage is based on the extent in the 2000s, based on the inventory of Nuth et al. (2013).Any points which have a fractional coverage between 10 and 90% are calculated with both the glaciated and non-glaciated land scheme.
The final results are then weighted based on the fractional glacier coverage.
CARRA is based on the non-hydrostatic numerical weather prediction model HARMONIE-AROME (Bengtsson et al., 2017).It uses ERA5 reanalysis (Hersbach et al., 2020) as boundary conditions and downscales it to a 2.5 x 2.5 km resolution over the European Arctic.The simulations are divided into two domains.Here we use the east domain, which contains Svalbard, Franz Josef Land, Novaya Zemlya and Northern Norway.CARRA currently spans the time frame from September 1990 to December 2021.
Similar to CARRA, AROME-ARCTIC (e.g.Müller et al., 2017) is also based on HARMONIE-AROME, and provides operational forecasts at a 2.5 x 2.5 km resolution over the Barents sea region.It uses ECMWF HRES forecasts as a lateral boundary conditions.The model has been operated by the Norwegian Meteorological Office since October 2015 and provides 66-hour forecasts with hourly resolution every 6-hours.Since this is a real-time forecast product, there are occasionally gaps in the forecast.When possible, we use the forecast initialized at 18UTC, as most data is assimilated at this time.We use a 6-hour lead time and extract data for 24 hours at a time, thus using forecast timesteps 6-30 for the simulations.This is chosen to optimise the prediction quality as well as to avoid spin-up effects.When the 18UTC forecast is not available, we use longer lead times of previous forecasts to find the most recent available estimate at a given hour.In the rare case that no forecast is available for the desired period, we simply interpolate between the previous and following available timestep.
Since CARRA and AROME-ARCTIC are on slightly different grids, we bilinearly interpolate the AROME-ARCTIC fields onto the CARRA grid in order for the final dataset to be consistent.

In-situ data
For evaluation of the AROME-ARCTIC and CARRA forcing data, observations from 26 Automatic Weather Stations (AWSs) are used: 6 stations on glaciers and 20 stations on non-glaciated land (see Fig. 1 and Table 1).The 20 stations on non-glaciated land are all operated by the Meteorological Office in Norway (MET-Norway) and have been assimilated into the AROME-  In addition, observations from mass balance stakes are used for evaluation of the CryoGrid products.A total of 52 measurement points are used, spread over 8 glaciers and ice caps (Table 1 and Fig. 1).The stake heights are recorded once or twice a year (typically in April and September), and are converted into summer and winter mass balance estimates using snow density and snow depth data.Stake data on Austre Brøggerbreen (BRG), Midtre Lovénbreen (MLB), Kongsvegen (KNG), Holtedahlfonna (HDF), and Linnébreen (LNB) have been collected by the Norwegian Polar Institute (e.g.Hagen et al., 1999), with the oldest record dating back to 1967.The Polish Academy of Sciences have measured mass balance stakes on Hansbreen (HBR) since 1989 (Grabiec et al., 2012).University of Oslo and the Norwegian Polar institute started mass balance measurements on Austfonna (AUS) in 2004 (e.g.Aas et al., 2016), while Uppsala and Utrecht universities initiated stake measurements on Nordenskiöldbreen (NSB) in 2006 (e.g.Van Pelt et al., 2012).

Satellite observations
In addition to the in-situ measurements of mass balance performed by stake measurements, we use estimates of the geodetic mass balance for validation of the CryoGrid product.The geodetic mass balance is found by taking the difference between elevation data at different dates to find the change in volume.This volumetric change is then converted into mass balance by assuming a value for the bulk density.Unlike the climatic mass balance estimates provided in this study, geodetic mass balance includes frontal ablation from marine terminating glaciers.Therefore, we only compare our results to the geodetic balance of land-terminating glaciers.
Several studies have provided estimates of the geodetic mass balance of glaciers in Svalbard (e.g.Moholdt et al., 2010;Nuth et al., 2010), but here we use the estimate by Hugonnet et al. (2021) which used Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) imagery for determining the geodetic mass balance of all glaciers on Earth from 2000-2019.
The results are available for all glaciers in the Randolf Glacier Inventory (Pfeffer et al., 2014)  The model has a modular structure, with many different modules that can be added together in various combinations to represent a wide range of surface and subsurface conditions.Information about the different functionalities and structures are described in detail in Westermann et al. (2022).
Three modules are used to determine the stratigraphy of glaciers in Svalbard: a glacier (ice) module, a firn module and a snow module (see Fig. 2).The main components of each module are described below.All modules use the surface energy balance as an upper boundary condition.For simulations of seasonal snow, a simple ground module and a snow module is used (see Westermann et al. (2022) for details on the ground module).

Glacier module
The glacier module consists of layers of pure-ice with a user-defined constant ice thickness.For this study, 47 layers with a thickness between 0.1 and 1 m were used, totaling 20 m of ice.We assume an ice albedo of 0.4 for all of Svalbard.When mass is removed from the model by runoff, evaporation or sublimation, mass is shifted up from an infinite ice reservoir below into the lowest model layer.This is done to prevent the glacier from disappearing during long spin-ups due to the lack of ice flow.
The infinite reservoir is assumed to have the same temperature as the lowest model layer.If there is no snow on the surface, any excess water from rain or melt runs off instantaneously.

Snow and firn module
If snowfall is added to the model, a snow module is added on top of the glacier ice or firn (see Fig. 2).If the snow survives on the glacier surface for more than one year, the snow layer is moved to a firn scheme.The snow and firn scheme have almost the same model physics, but new snow will never be mixed with a firn layer.The snow and firn modules follow a slightly altered CROCUS (Vionnet et al., 2012) snow scheme as described in Westermann et al. (2022).For this study, we mostly follow this scheme but with updated runoff and water percolation schemes.These modifications are described below, while a brief description of the albedo, temperature diffusion, and densification are given in Supplement S1.

Water percolation and runoff
The water in a grid cell is either immobile and bound to the snow or firn, or it flows downwards driven by gravity.The limit between the two regimes is the field capacity θ f c , in this study chosen as 0.05.The vertical water flux q w [m s −1 ] is therefore given by where K is the hydraulic conductivity [m s −1 ] and θ w is the water saturation.The hydraulic conductivity is parameterized in terms of the snow grain diameter d [m], the snow density ρ s , and the effective liquid saturation The hydraulic conductivity of snow is the product of the unsaturated conductivity, K r , and saturated conductivity, K s , i.e.
The saturated hydraulic conductivity (Shimizu, 1970) is given by where g is the gravitational acceleration [m s −2 ] and v w = 1.787 • 10 −6 m 2 s −1 is the kinematic viscosity of water.The unsaturated hydraulic conductivity (van Genuchten, 1980) is given by where the parameter n is given by n = 15.68 exp(−460d g ) + 1. (4) Water is not allowed to flow into an impermeable layer, here defined as layers with a density higher than 830 kg m −3 (Cuffey and Paterson, 2010), or a layer that already has its entire pore space filled.Water which would have otherwise flowed into an impermeable layer becomes available to run off.Runoff does not occur immediately, but depends on a characteristic local runoff scale τ R [days] which increases with surface slope S [m m −1 ] as follows, where c 1 = 0.33 day, c 2 = 25 days, and c 3 = 140 (Lefebre et al., 2003).The runoff per timestep R [m w.e.] is then calculated from the water in excess of the field capacity W ex [m], as where ∆t is the time step in days.This delay in runoff means that water in excess of irreducible saturation may linger in a layer until it either refreezes or runs off.

Vertical discretization
To avoid very thin snow layers, a simplified gridding scheme is used (Zweigel et al., 2021).During each timestep, new snow is added to the uppermost grid cell by calculating a weighted average between all variables describing the new and old snow (density, snow age, snow grain size etc).The water equivalent volume of snow is used as the weighting factor.When the top grid cell exceeds a target snow water equivalent (here 0.02 m) by more than 50%, it is split in two.If the top grid cell is smaller than 50% of the target snow water equivalent, it is merged with the below cell.The grid size of the top snow cell is therefore in the order of 0.01-0.03m w.e.For deeper snow layers, the layer size doubles every 10 layers by splitting/merging of layers.
For the firn modules, the top layer has a maximum snow water equivalent thickness of 0.1 m, and the layer size doubles every 10 layers by merging/splitting of layers.Freshly fallen snow will always fall on top of the firn and never be mixed in with the top layer.

Forcing evaluation
The evaluation of both forcing products against available AWS observations shows that both models generally fit well with observations, but that the bias and root-mean-square-error of the CARRA product is generally smaller than for AROME-ARCTIC.For detailed evaluation of the model forcing against available AWS observations for both CARRA and AROME-ARCTIC, in addition to a discussion on the inter-comparison, we refer to Supplement S2,S3.

Mass balance evaluation 260
Mass balance, b, from stakes on eight glaciers in Svalbard is compared to the CryoGrid simulations of the surface mass balance in Table 2. Here, the surface mass balance is defined as the mass balance in the annual layer, and thus does not include refreezing in firn.The difference in mass balance, ∆b, is defined as the modeled value minus the observed value at a given location.4 compare the stake observations and the nearest model gridpoint value for the CARRA-forced CRYOGRID simulations.Overall, there is a good agreement between the model and the observations, with biases and root-mean-square errors similar to those found in other modeling studies (e.g.Østby et al., 2017;Van Pelt et al., 2019).The largest difference in the winter mass balance occurs at Hansbreen (Fig. 4e) where the model has a large negative bias at all stake locations except at the lowest and highest elevations.There is also a negative bias in the summer mass balance at the low elevation stations, but there is good agreement at glacier stations at higher elevations.The largest average difference in summer occurs at Austre Brøggerbreen, where the CARRA-forced simulation underestimates the mass balance by 0.22 m w.e. on average.Note, however, that both Hansbreen and Austre Brøggerbreen are small glaciers in complex topography, and thus may not be well represented by the the 2.5 x 2.5 resolution of these simulations.
Table 2 also contains the comparison between the stake observations and the nearest model gridpoint for the AROME-ARCTIC-forced simulations.Only the 2016/17 and 2017/18 glaciological years were used for this evaluation.Overall, the CARRA and AROME-ARCTIC forced simulations perform almost equally well over these two years, with similar biases and root-mean-square errors for both the summer and winter balance.There is, however, a larger difference between the estimates of the annual mass balance, primarily due to large differences in the simulations for Nordenskiöldbreen when using the different forcings.For the CARRA forced runs for 2016/17-2017/18, the overestimation in the mass balance of Nordenskiöldbreen during the winter is balanced by excess melt during the summer, leading to only a small bias in the annual comparison.Using AROME-ARCTIC, the mass balance of Nordenskiöldbreen is underestimated both in summer and winter, leading to a large bias and rmse in the annual comparison.
In addition to the in-situ mass balance, we use estimates of the geodetic mass balance of land-terminating glaciers by Hugonnet et al. (2021) to validate the mass balance results.Since the geodetic estimates include refreezing below the annual layer, we here use the climatic mass balance for the comparison.Figure 5 compares the cmb from CryoGrid for five year periods between 2000 and 2020 against estimates from Hugonnet et al. (2021).The simulated cmb is within the uncertainty estimate of the geodetic data for the whole period, except in 2005-2009 where the cmb is slightly higher than the uncertainty estimate (by 0.02 m w.e.yr −1 ).The AROME-ARCTIC-forced simulations are within the uncertainties of the geodetic estimate, but have a slightly lower mass balance than the CARRA-forced simulations for the same period (-0.29 m w.e.yr −1 using CARRA vs. -0.34 m w.e.yr −1 using AROME-ARCTIC.)

Climatic mass balance
The area-averaged climatic mass balance of all Svalbard glaciers for the whole CARRA simulation period is found to be -0.08 m w.e.Fig. 6a-c show the annual, winter, and summer climatic mass balance over Svalbard.The results are shown for each mass balance year, here defined as September to August.The most negative values are found at low-elevation areas in S and SW Spitsbergen (Fig. 6a), with the cmb reaching down to -3.23 m w.e.yr −1 , while the most positive values are found at highelevation areas in central Spitsbergen, reaching a maximum cmb of 1.16 m w.e.yr −1 .The winter mass balance (Fig. 6b) is on average positive in all points, while the summer mass balance (Fig. 6c) is negative except in some high-elevation points in NE and NW Spitsbergen.

e.).
There is no significant trend in winter, summer, or annual cmb over the investigated period.
Figure 6d shows the results using both CARRA and AROME-ARCTIC forcing.Generally, the AROME-ARCTIC simulations contain slightly lower winter cmb, with an average difference of -0.04 m w.e..The summer cmb is more variable, but 305 generally the values in the AROME-ARCTIC simulations are more negative than CARRA (by -0.05 m w.e. on average).Overall, the annual cmb is similar in the two simulations, but with a more negative cmb when using AROME-ARCTIC of about -0.1 m w.e. from 2016-2017, and a more positive cmb when using AROME-ARCTIC in 2019-2021 of approximately 0.1 m w.e..
Based on AROME-ARCTIC, 2021/22 appears to be a record negative mass balance year for Svalbard.
Figure 7 shows the winter (blue bars), summer (red bars), and annual cmb (green line) for eight different regions of Svalbard.
The glaciers in Nordenskiöldland have the most negative annual cmb (-0.73 m w.e.yr −1 ), with glaciers losing mass during all years except 2007/08.The most positive average cmb is in NE Spitsbergen (0.11 m w.e.yr −1 ).For all areas, 2012/13 was a year with a strongly negative summer cmb.In 2019/20, NW-Spitsbergen experienced a record amount of melt (-1.37 m w.e.) in combination with a record low winter cmb (0.18 m w.e.).Most other regions also experienced strong summer melt, with the exception of Edgeøya and Barentsøya where the summer cmb is close to the average over the simulation period.

Seasonal snow
During the summer season, the non-glaciated land areas on Svalbard experience snow-free conditions.The length of the snowfree season depends on the snow disappearance day and the snow onset day, defined as when the snow thickness first drops below (snow disappearance) or above (snow onset) a certain threshold.Here, we set the threshold to 1 cm.During the simulation period, no significant trend is found for the snow disappearance date.There is, however, a significant trend in the snow onset date of 4.0 days decade −1 (p=0.05) and the length of the snow-free season of 6.8 days decade −1 (p=0.02) (Fig. 8).
When considering the different regions as defined in Fig. 7, we find that some areas show significant trends in the snow onset or disappearance date.Nordenskiöldland has a significant trend in both the snow disappearance date (-5.7 days decade −1 , p < 0.01), and the snow free period (11.0 days decade −1 , p < 0.01).NW Spitsbergen, there is a significant trend in the snow onset date (4.3 days decade −1 , p=0.02) and the snow free period (7.2 days decade −1 , p < 0.01).S. Spitsbergen there is a trend in the snow disappearance date (-6.2 days decade −1 , p < 0.01) and the snow free period (10.6 days decade −1 , p < 0.01).
The simulations using AROME-ARCTIC forcing give similar results as those using CARRA for the snow free period, snow onset date, and snow disappearance date.The snow free period is shorter by on average 8.4 days when using AROME-ARCTIC, which is due to a slightly later snow disappearance date (on average 4.2 days later) and a slightly earlier snow onset date (on average 4.2 days earlier).

Refreezing
Refreezing is defined as all liquid water that refreezes within snow and firn, without taking into account that this may melt again.The average annual refreezing for glacier-covered and land areas is 0.24 m w.e. and 0.11 m w.e., respectively.The lowest annual refreezing is simulated at low elevations, where only thin seasonal snowpacks are present, thus limiting the amount of refreezing.The largest refreezing is found in the accumulation zones (Fig. 9), where average values up to 0.32 m w.e.yr −1 are simulated.In these areas, water can percolate down into firn layers and refreeze over the winter season.The spatial distribution of this internal accumulation (defined as refreezing beneath the annual layer) is shown in Fig. 9b, showing that a significant fraction of the refreezing at higher elevations occurs in deeper layers.The average annual internal accumulation is 0.11 m w.e., and thus accounts for almost half of the total refreezing (Fig. 9c).There is a significant negative trend in the refreezing within glaciers of -13 mm w.e.decade −1 (p < 0.01), which is primarily due to a decrease in internal accumulation.
For glacier-covered areas, annual refreezing of melt and rainwater accounts for 25% of the total accumulation, varying between 19% and 32% over the simulation period.There is a significant negative trend in the contribution of refreezing to the total accumulation of -2.1% decade −1 (p < 0.01).
There is a higher amount of refreezing in the AROME-ARCTIC-forced simulations between 0.01-0.05m w.e.yr.This is most likely due to the higher temperatures in AROME-ARCTIC over some glaciers, leading to higher amounts of liquid precipitation in the spring and autumn.The internal accumulation is almost equal in the two simulations (the mean difference between the two simulation is <0.006 m w.e.).The amount of refreezing in seasonal snow is also similar overall, with an average difference of 0.01 m w.e.

Runoff
The simulated runoff from both glacier-covered and non-glaciated land points are shown in Fig. 10.The average runoff for glaciers has a similar pattern as the cmb (Fig. 6), with the highest runoff in low-elevation regions (up to 3.0 m w.e.yr −1 ) and lowest runoff in high-elevation areas in central Spitsbergen (down to 0.02 m w.e.yr −1 ).
The average total runoff from glacier-covered regions is 29 Gt yr −1 and for land regions it is 12 Gt yr −1 .While there is a large variation in runoff from glacier-covered regions, the runoff from land areas is relatively stable throughout the whole period (Fig. 10c).The minimum and maximum runoff from seasonal snow occurred in 1996 (9.5 Gt) and 2016 (19 Gt), respectively.For glacier covered areas, the minimum runoff occurred in 2008 (14.4 Gt), where the discharge was almost equal to that coming from seasonal snow.The runoff during this year was low for the entire Svalbard area, with particularly low rates along the Western coast.The largest runoff from glaciers occurred in 2013 (50.4 Gt), closely followed by 2020 (49.0 Gt).In 2013 there was generally high runoff over the entire peninsula compared to the average values, with especially large runoff rates in southern Spitsbergen and Barentsøya.
There is a significant, positive trend in both the glacier runoff and the runoff from seasonal snow of 5.2 Gt decade −1 (p=0.01) and 1.1 Gt decade −1 (p<0.01),respectively.Comparison to the runoff simulated using AROME-ARCTIC is shown in Fig. 10b,c.The glacier runoff is generally higher 375 in the AROME-ARCTIC-forced simulations for SW Spitsbergen and Barentsøya/Edgeøya, and lower for Kvitøya and Nordenskiöldland.The average difference from 2016-2021 is 0.4 Gt yr −1 , ranging between 2-6 Gt for individual years.The runoff from seasonal snow on non-glaciated land is also similar overall, with the average value from AROME-ARCTIC only slightly larger than CARRA by 0.6 Gt yr −1 .Interestingly, although the glacier runoff is lower in the AROME-ARCTIC simulations for Nordenskiöldland, the runoff from seasonal snow is not.This indicates that the difference between the CARRA and AROME-380 ARCTIC runoff estimates is not due to differences in precipitation.In the first half of the 2010s, the calving rate is up to -0.43 m yr −1 , while in the latter half it is up to -0.5 m yr −1 .The large range in calving rates reflects the uncertainty in the geodetic estimate.

Almost real-time estimates with AROME-ARCTIC
As shown above, AROME-ARCTIC forcing generally produces results for the mass balance and runoff of Svalbard that is similar, within 0.2 m w.e.yr −1 for both variables, to that simulated using CARRA forcing.This indicates that AROME-ARCTIC can be used to create near real-time estimates of the climatic mass balance of Svalbard, although the uncertainties may be larger than generated by the CARRA-forced simulations.Since AROME-ARCTIC output is made available daily on the Norwegian Meteorological Institute server, it is possible to create an automatically updating product of mass balance, runoff, etc. for Svalbard.However, since this is an evolving product, technical difficulties may occur if data formats or naming conventions change, or if the forecast files are missing for longer than a few days.

Uncertainty
Several sources of uncertainty are introduced through the creation of the glacier mass balance dataset in this study.The sources of these uncertainties are comprised of the model physics, the initial model state, atmospheric forcing, glacier extent and topographic simplification.It is, however, difficult to quantify the contribution of each individual source.This section discusses these sources of uncertainty.

Model physics and initialisation
Although the snow and firn scheme is based on the CROCUS model, the physics of which has been used and validated in a number of glacier mass balance and snow studies (Cullather et al., 2016;Schmidt et al., 2017;Verjans et al., 2019), there may still be uncertainties connected to using this model in Svalbard.For example, previous studies have shown that CROCUS does not always perform well under Arctic conditions, therefore, we have made a number of changes to the original model as e.g.
suggested by Royer et al. (2021).However, most of the model variables are based on recommendations from previous studies and have not been tuned for the conditions of Svalbard.Although the model does well when compared to observations, potential biases may arise in other regions of Svalbard.To initialize the sub-surface conditions, a 30-year spin-up was performed.This was done by repeating the forcing from 1991-2000, until the model output was approximately in balance with the applied climate forcing.This could introduce some biases in both the extent and depth of the firn area, as the glaciers may not have been in balance with the 1991-2000 climate in reality.

Model forcing
We find good agreement between CARRA forcing and the meteorological variables and the incoming radiation over glaciers and non-glaciated land (see Supplement S2,S3 for details), although it should be noted that the non-glaciated land based AWSs are assimilated into the CARRA product.Comparison against winter mass balance stake observations shows CARRA precipitation has a low rmse overall (0.21 m w.e.), but is slightly underestimated over most of the glaciers (Table 2).This could be partly due to the spatial resolution of CARRA, as at a 2.5 km resolution the model might miss some of the impact of the terrain on the precipitation distribution, particularly in areas with complex topography.In addition, the simulated mass balance representing a 2.5 x 2.5 km cell may not be directly comparable with point-observations, as heterogenities in the energy and mass balance occur at spatial scales less than 2.5 km.For example, in areas with high wind, redistribution of snow by the wind may have a large effect on the winter mass balance (e.g.Winther et al., 2003).Furthermore, since the stake observations are mainly taken along the glacier centerline, the observations do not reflect the horizontal distribution of the mass balance along In addition, using AROME-ARCTIC to generate a real time dataset adds additional uncertainties as this is a forecast and not a reanalysis product.From 2016 until the summer of 2019, the model was initialized with too little snow over some glacier points in the ablation area, thus leading to unrealistically high surface and 2m temperatures.To try to counter this effect, we use the 10 m temperature for the AROME-ARCTIC-forced simulations when unrealistically high surface temperatures occur, but some biases may still persist.In addition, as previously discussed, due to missing data it is not always possible to use AROME-ARCTIC forecasts with a 6-hour lead time.Using an earlier forecast with longer lead times introduces higher forecast errors, and therefore using forecasts at different time steps may give different results.As an example, Fig. 12 shows the differences in incoming radiation between various forecast lead times during August 2019.between forecast lead times of 6, 12, 18 and 24 hours is shown in the Fig. 12b.The grey area shows the maximum and minimum differences between the incoming radiation with a 6 hour lead time and 12, 18, and 24 hour lead times.The mean differences between the 6-hour and 12/18/24-hour lead times over the whole month is small, ranging between -8.2 to 8.4 W m −2 , but at any given time step large differences (> 100 W m −2 ) occur.Similar effects can be seen in the other meteorological variables, like precipitation, wind speed, and temperature.We expect the effect of the lead time used to be small over monthly or yearly timescales, but it can introduce large errors for specific days or areas.
Furthermore, the re-gridding of the AROME-ARCTIC product to the CARRA grid using linear interpolation may introduce additional errors.

Glacier extent and topography
Throughout the simulation period, we assume the elevation and glacier mask is fixed thus neglecting the effect of ice flow and elevation changes on the mass balance.Both the elevation and the glacier mask are based on observations collected between 2000-2010 and should therefore be representative for most of the investigated period.
The added error of a fixed glacier mask has previously been investigated by Østby et al. (2017).The authors found that the error in the climatic mass balance associated with using a fixed glacier mask (based on observations from the 2000s) as opposed to a time-varying mask was on average 0.02 m w.e.−1 for the period 1957-2014.Since the period investigated in this study is smaller and more closely matches the time period the glacier mask was created, we expect the error due to a fixed glacier mask in our simulations to be equal or smaller than the value found by Østby et al. (2017).
Using a fixed elevation mask may introduce a negative bias in the beginning of our study period, as the elevation may be too low, and a positive bias towards the end of the study period where the used elevation mask may be too high.On average for Svalbard, the glaciers elevation decreased at a rate of 0.36 m yr −1 from 2000-2020 (Hugonnet et al., 2021), while between the mid-1960's and 2005 the glacier elevation outside Austfonna and Kvitøya decreased, on average, at a rate of 0.49 m yr −1 (Nuth et al., 2010).Considering the elevation map used in this study is based on observations from the 2000s, we expect the maximum average deviation to be 10 m.Assuming a change of mass balance with elevation of 3 • 10 −3 m w.e.m −1 , we expect the error associated with the constant glacier mask to be less than 0.03 m w.e.yr −1 .

Comparison with other studies
Several other studies have previously quantified the Svalbard-wide mass balance and runoff.Direct comparison between our results and other studies is in some cases hampered by differences in time period, areal coverage, and the type of mass balance calculated (e.g.estimates from gravimetry or geodetic methods will estimate the total mass balance, including frontal ablation).
Here, we only compare against studies which calculate either the climatic or surface mass balance, and who have published results within our simulation period of 1991-2020.When available, we compare simulated average 2m temperature, yearly precipitation, climatic mass balance, and runoff (Table 3 and Fig. 13).
The climatic mass balance simulated in this study is similar to estimates from other studies.The cmb in this study is slightly less negative than that simulated by Lang et al. (2015), which could partly be due to the higher precipitation in this study.There  are larger differences between our cmb and the estimates by Aas et al. ( 2016) and Østby et al. (2017), where our simulated cmb is lower by 0.22 and 0.05 m w.e.yr −1 , respectively.In the case of Østby et al. (2017), this is partly due to a big difference between the estimates for 2013, where the cmb estimated by Østby et al. (2017) is strongly negative.A more negative mass balance is consistent with the higher average 2m temperatures used for their simulations.On the other hand, our simulations have a more negative cmb than those simulated by Van Pelt et al. (2019).This is likely related to the much lower precipitation in our simulations.When comparing the estimates of specific runoff, our results are consistent with Van Pelt et al. (2019).The precipitation, cmb, and runoff values in this study are consistent with those of Noël et al. (2020).
The temporal evolution of each of the cmb studies are plotted against our results in Fig. 13.The temporal pattern is similar for all the estimates, with high inter-model correlations between 0.8 and 0.9.
https://doi.org/10.5194/egusphere-2022-1409Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.ARCTIC and CARRA products.The 6 glacier stations, on the other hand, have not been assimilated and thus provide independent reference.The glacier stations are located on: Etonbreen, operated by University of Oslo and the Norwegian Polar institute since 2004, (e.g.Schuler et al., 2014), Kongsvegen, operated since 2007 by the Norwegian Polar Institute (Kohler et al., 2017), Vestfonna, operated for two years between 2007-2009 by Uppsala University (Jonsell et al., 2010), and Nordenskiöldbreen and Ulvebreen, operated by Utrecht University since 2009 and 2015, respectively.The measurement interval was between 1-2 minutes, depending on the station.Daily mean observations of the 2m temperature, 2m relative humidity, 10m wind speed, and incoming and outgoing longwave and shortwave radiation is used for the evaluation, when available.

Figure 1 .
Figure 1.The location of the surface mass balance stakes and automatic weather stations used.The blue shaded areas are used for comparison with geodetic mass balance estimates.
at a temporal resolution of 1, 2, 4, 5, 10, and 20 years.Here, we use the 5 year mass balance estimate for all land-terminating glaciers in Svalbard for model comparison (see Fig. 1).https://doi.org/10.5194/egusphere-2022-1409Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.In this study we use and further develop the CryoGrid community model for simulations of the climatic mass balance and meltwater runoff.CryoGrid is an open-source model developed for climate-driven simulations in the terrestrial Cryosphere.

Figure 2 .
Figure2.Evolution of CryoGrid stratigraphy for simulation of (a) glacier mass balance and (b) seasonal snow in this study.For glacier grid points (a), the model stratigraphy consists of up to three modules, which can be combined to represent the following four situations: glacial ice, glacial ice covered with snow, glacial ice covered by firn and snow, and glacial ice covered by firn.Between each module there is an interaction (IA) class, which determines the transfer of heat, water, and mass.A trigger function determines when modules can be added or removed from the stratigraphy.For seasonal snow (b), there is a ground module which can be coupled to a snow module.

5
Figure3shows the average yearly temperature and precipitation in CARRA over 1991-2021, as well as significant trends 245

Figure 3 .
Figure 3. Average a) 2m-temperature and c) precipitation over 1991-2021 in CARRA.Significant (p < 0.05) b) temperature and d) precipitation trend in each points.Stippled areas have no significant trend.

Figure 5 .
Figure 5. Geodetic mass balance of land terminating glaciers (from Hugonnet et al. (2021)) compared to the climatic mass balance simulated in CryoGrid.

Figure
Figure 6d shows the temporal evolution of the summer, winter, and annual cmb.The most positive cmb (0.43 m w.e.) was found in the 2007/08 mass balance year, while the most negative cmb (-0.68 m w.e.) was found in 2019/20.The winter cmb is on average 0.44 m w.e.yr −1 , with a maximum in 2015/16 (0.65 m w.e.) and a minimum in 2001/02 (0.28

Figure 6 .
Figure 6.Climatic mass balance of Svalbard from 1991/92 to 2020/21.The top row contains maps of the average (a) annual, (b) winter, and (c) summer cmb, while (d) shows the temporal evolution of the summer, winter and annual cmb for each mass balance year (September -August).Dark red and blue show the results using CARRA forcing, while lighter colored bars show the results using AROME-ARCTIC (denoted AA).

Figure 7 .
Figure 7. Climatic mass balance for different regions of Svalbard.Dark red bars show the summer balance and dark blue bars show the winter balance calculated with CARRA forcing.The lighter colored bars from 2021/22 is the climatic mass balance calculated with AROME-ARCTIC forcing.The green line is the yearly cmb simulated using CARRA forcing, while the black line is the yearly cmb using AROME-ARCTIC.

Figure 7
Figure 7 also show the yearly cmb of each region calculated with AROME-ARCTIC (black line).For most regions, the cmb simulated using AROME-ARCTIC closely matches that of CARRA, although large deviations for Nordenskiöldland exist (Fig.7h).For the other regions, the average difference between the CARRA and AROME-ARCTIC estimates is < 0.10 m w.e., while for Nordenskiöldland it is 0.41 m w.e..For 2021/22, the AROME-ARCTIC forced simulations show a highly negative 325

Figure 8 .
Figure 8. a) The length of the snow free period and b) the snow disappearance and snow onset date.In both plots, the thick lines shows the average for all non-glaciated land on Svalbard using CARRA forcing, while the stippled line shows the average using AROME-ARCTIC.The colored areas show the span in values for the different regions of Svalbard (see Fig. 7), i.e. the top line shows the region with the maximum average disappearance/onset date and the lower line shows the region with the minimum average disappearance/onset date.

Figure 9 .
Figure 9. (a) Average annual refreezing for glaciers and seasonal snow on non-glaciated land areas.(b) Average annual internal accumulation from glacier-covered areas.(c) The temporal variation in refreezing and internal accumulation for glaciers and seasonal snow.Solid lines show the CARRA forcing and dashed lines show AROME-ARCTIC forcing.

Figure 10 .
Figure 10.a) Average runoff over the whole CARRA simulation period (1991/92-2020/21).b) The difference in runoff between the AROME-ARCTIC and CARRA forced simulations from 2016-2021.c) Timeseries of runoff from glaciers and seasonal snow on nonglaciated land.Dashed lines show the results from the AROME-ARCTIC-forced simulations.

Figure 11
Figure 11 shows the daily accumulated mass balance from CARRA and AROME-ARCTIC forced simulations.The mean of the CARRA-forced simulations from 1991-2021 is shown with a dashed black line, while the minimum and maximum years are shown in grey.The 2020/21 and 2021/22 mass balance years calculated with AROME-ARCTIC are shown in red and blue.In order to better compare the CARRA and AROME-ARCTIC forced simulations, the AROME-ARCTIC-forced estimates are shown with an uncertainty, given as two standard deviations of the differences between the CARRA and AROME-ARCTIC forced simulations from 2016/17-2020/21.In other words, this uncertainty indicates what the cmb would most likely be if CARRA had been used as forcing as opposed to AROME-ARCTIC for these simulations.During the winter months, there is generally little difference between the simulations with the difference forcings, and we can therefore have a high confidence in the AROME-ARCTIC results.During the summer, however, larger differences arise between the different products, which accumulate over the melt season.Based on the 2016/17-2020/21 simulations, the estimated standard error in the cmb is 0.12 m w.e. by the end of the mass balance year.

Figure 11 .
Figure 11.The span in daily accumulated cmb from 1991/92-2021/22 simulated using CARRA and AROME-ARCTIC climate forcing.The 2020/21 and 2021/22 mass balance years are simulated using AROME-ARCTIC (shown in red and blue).The CARRA result from 2020/21 is additionally shown in yellow.The uncertainty of the AROME-ARCTIC estimate is defined as two standard deviations of the differences between the CARRA and AROME-ARCTIC forced simulations from 2016/17-2020/21.

Figure 12 .
Figure 12.(a) Mean difference in August incoming radiation (R↓) between 6-hour and 24-hour forecast lead times.(b) The temporal difference in incoming radiation between 12, 18, and 24 hour lead times versus a 6 hour lead time (grey area).Location of point is shown with a circle in (a).

Fig
Fig.12ashows the mean difference between 6 hour and 24 hour lead times.Overall, the mean absolute difference is small (1.5 W m −2 ) with a maximum average deviation of 7.0 W m −2 .However, at any given location and time step, the difference in incoming radiation between the different lead times may be a large as 335 W m −2 .An example of the temporal difference ] [m w.e.yr −1 ] [m w.e.yr −1 ] [m w.e.yr −1

Figure 13 .
Figure 13.Time series of yearly cmb from different model studies from 1990-2021.

Table 1 .
In-situ data used for evaluation of the forcing and mass balance products.UiO: University of Oslo, NPI: Norwegian Polar Institute, IMAU: Institute for Marine and Atmospheric research Utrecht, PAN: Polish Academy of Sciences, UU: Uppsala University.

Table 2 .
Evaluation of modelled results against observations from mass balance stakes in m w.e.yr −1 .The values from 1991/92-2017/18 shows the comparison with CARRA-forced model simulation, while for 2016/17-2017/18 the observations are compared to simulations using both CARRA and AROME-ARCTIC forcing.The results are given as CARRA forcing / AROME-ARCTIC forcing.Subscripts w, s, and a, respectively, refer to values calculated for winter months, summer months and annually.https://doi.org/10.5194/egusphere-2022-1409Preprint.Discussion started: 8 February 2023 c Author(s) 2023.CC BY 4.0 License.Table 2 and Fig.

Table 3 .
Svalbard climatic variables (precipitation and 2m temperature), climatic mass balance, and glacier runoff from different modeling studies. 1 only precipitation over glaciers.