Observed and Predicted Trends in Icelandic Snow Conditions for the period 1930-2100

. This study presents an estimate of historical snow conditions in Iceland and a projection of these conditions, given different emission scenarios. Historical snow conditions were estimated using in situ observations from manned meteorological stations over the period 1930-2021 and by remotely sensed observations from the MODIS instruments over the period 2001-2021. Historical and future climate conditions, as described by each of the 21 General Circulation Models (GCM’s) from the 10 5 th iteration of the Coupled Model Intercomparison Project (CMIP5) as contained in the NASA Earth Exchange (NEX) Global Daily Downscaled Projections (GDDP) dataset, were used to simulate snow conditions in Iceland over the period 1950-2100 under the Representative Concentration Pathways (RCP) RCP4.5 and RCP8.5 with the SNOW-17 model. The results show an increase in the average annual Snow Cover Frequency (SCF) over the historical record detected both in the in-situ (1930-2021) and remotely sensed data (2001-2021). Average annual snow depth measurements also revealed an increasing trend over the 15 historical record. Simulated snow conditions show a substantial decrease in both Snow Water Equivalent (SWE) and SCF over the period 1950-2100, a trend more pronounced under RCP8.5 as compared to RCP4.5.


Introduction
Icelandic climate is categorized as maritime, with mild winters, cold summers, strong winds, frequent precipitation and large spatio-temporal variations in weather and micro-climate Ólafsson et al., 2007). It is strongly influenced 20 by ocean conditions in the Northern Atlantic (e.g. Massé et al., 2008) and mass balance trends of Icelandic glaciers are highly correlated with changes in large-scale ocean circulations (Eythorsson et al., 2018). Since the last glacial maximum the average annual air temperature in Iceland has increased about 4°C (Geirsdóttir et al., 2013;Knudsen et al., 2008;Langdon et al., 2011;Larsen et al., 2011;Sicre et al., 2011). The average air temperature in Iceland has risen by 0.08°C/decade since the 1850's, comparable to the global average, and by 0.5 °C/decade over the period 1980-2016 (Bjornsson et al., 2018). Since 1890 the 25 Icelandic glaciers have lost about 16% of their mass and 18% of their surface area, contributing about 1.5 mm of global sea level rise (Aðalgeirsdóttir et al., 2020;Bjornsson et al., 2013) and are expected to lose most of their remaining mass over the next two centuries at current pace (Adalgeirsdottir et al., 2006;Bjornsson & Palsson, 2008;Jóhannesson et al., 2004;Schmidt et al., 2020). Runoff in Iceland is generally expected to increase in winter as less water is stored in snowpack and runoff from which is expected to vary depending on ocean conditions in the North Atlantic, where recent cooling has led to a slowdown in mass loss of Icelandic glaciers (Noël et al., 2021). Spring melt is generally predicted to begin earlier and autumn snow cover to occur later (Johannesson et al., 2007). Analysis of a recently developed gap filled MODIS snow cover product suggests that the snow cover duration has increased during the period 2000-2018 for all months expect October and November (Gunnarsson et al., 2019). Understanding of future expected changes to snow in Iceland is important for water resources management as it 35 constitutes a significant portion of the regional hydrological cycle, especially in the interior highlands where the majority of the country's energy production occurs, in hydropower plants fed by glacial rivers.
Snow cover monitoring by satellite remote sensing has been studied since the 1960s and several global snow cover products have been produced based on these observations. (Dong, 2018;Frei et al., 2012;Robinson et al., 1993). The MODIS instruments on the Terra and Aqua satellites (Dietz et al., 2012) provide a good balance of spatial and temporal resolution, 40 with two daily observations and 500m x 500m pixels (Aalstad et al., 2020). An imporant variable for snow remote sensing is the Snow Cover Frequency (SCF), the number of days with snow cover divided by the number of valid observations per year (Nolin et al., 2021) which is related to e.g. growing season length and habitability (Callaghan et al., 2011). SCF is a key parameter in the Earth's energy balance (Cohen, 1994) and can be used to analyze the impacts of climate change on the cryosphere (Brown & Mote, 2009). 45 Snow condition estimates by both General Circulation Models (GCM) and Regional Circulation Models (RCM) capture the main traits of annual snow cycle but are know to contain biases due to their relatively simple snow schemes (Frei et al., 2018;Matiu & Hanzer, 2022). In general GCM-RCM pair ensembles continuation of the ongoing reduction in average snow conditions until the middle of the 21 st century (Verfaillie et al., 2018). Improved estimates of snow conditions have been achieved e.g. using various re-analysis (e.g. Fiddes et al., 2019), downscaling (Fiddes et al., 2022;Smiatek et al., 2016) and 50 data assimilation methods applied either to GCM-RCM snow projections or projections of snow conditions by different snow models forced with downscaled and/or bias corrected GCM-RCM data (e.g. Hanzer et al., 2018).
Many snow models have been developed and described in the literature (e.g. Krinner et al., 2018;Magnusson et al., 2015).
The SNOW-17 model was developed for the US national Water Service where it has been used for operational snow forecasting for the past several decades (Anderson, 2006). The SNOW-17 model has been applied to several regional climate change 55 studies Notaro et al., 2014) and has shown good correlation to MODIS Snow-Covered Area (SCA) observations (Franz & Karsten, 2013) . A key advantage of the SNOW-17 model is that it is a conceptual model which simulates snow pack conditions based on a temperature index which is both more computationally efficient compared to full energy balance models and requires fewer and simpler forcing data variables.
The objective of this study was to analyze observed trends and predict the development of snow conditions in Iceland under 60 different plausible climate scenarios. It presents an analysis of historical and future trends in Icelandic climate and snow conditions. Improved understanding of how local snow resources are likely to respond to changing climate conditions is important as these changes are expected to impact local communities and ecosystems as well as changing the challenges and opportunities for exploiting natural resources in cold areas (Eliasson et al., 2017). In this study changes to historical snow cover properties were estimated based on both in-situ and remotely sensed observations. Future snow conditions were projected 65 by modelling based on a globally downscaled and bias corrected ensemble of General Circulation Models (GCM) from the 5 th iteration of the Coupled Model Intercomparison Project (CMIP5). The novelty of this study is the analysis of an extended dataset of in-situ records of snow conditions in Iceland combined with reliable remotely sensed dataset of snow conditions in the area and the comparison of these observations with snow conditions simulated using a trusted snow model run with downscaled and bias corrected temperature and precipitation estimates from an ensemble of 21 the CMIP climate models 70 ensemble on a freely available, cloud based, parallel computing platform.

In Situ Snow Observations
Data on in situ snow measurements at manned monitoring stations were acquired from the Icelandic Meteorological Office  The number of stations reporting snow data is below 10 until 1950 and rapidly increases thereafter, the number of stations recording snow cover status increases prior to those recording snow depth, from the 1960's snow depth has been recorded at more than 60 stations. Figure 1 (right panel) shows that the average annual snow depth from all stations has remained similar throughout the study period. 90

Remote Sensing and Geospatial Data
The MOD10A1.006 and MYD10A1.006 daily snow cover products from the MODIS instruments on NASA's Aqua and Terra satellites (Hall et al., 2016) were used to estimate spatial changes in snow cover over the period of the 2001-2021 water years.
The 'NDSI_Snow_Cover' band was used to estimate the presence of snow in each pixel. The band contains the value of the Normalized Difference Snow Index (NDSI) which leverages the fact that snow is highly reflective within the visible spectrum 100 but not the infrared (Painter et al., 2009). The NDSI values in each pixel are given in a range of 0-100 where a value of NDSI > 0 indicates the presence of some snow in the pixel and a value of 100 that the pixel is likely fully snow covered. The NDSI is not to be confused with Fractional Snow Covered Area (FSCA) which is a measure of the fractional snow coverage of a pixel although there exist commonly used transformations between MODIS FSCA and NDSI, including linear (Salomonson & Appel, 2006), inverse linear (Fiddes et al., 2019) and other regression methods (Alonso-González et al., 2021). Results 105 from the Sentinel-2 mission have also shown good performance of other NDSI-FSCA transformations (Gascoin et al., 2020).
Regression methods have shown good performance on medium resolution observations such as from MODIS, whereas higher resolution observations have been shown to benefit from spectral unmixing (Aalstad et al., 2020). The 'NDSI_Snow_Cover_Basic_QA' band was used to select observations by quality estimate. A 10 x 10m DEM (IslandDEM) was used for topographical information (National Land Survey of Iceland, 2016). 110

Climate Data
The NASA Earth Exchange (NEX) Global Daily Downscaled Projections (GDDP) dataset (Thrasher et al., 2012;Thrasher et al., 2006) was used as an estimate of historical and future climate. The dataset contains global minimum and maximum near surface air temperatures and surface precipitation rates, as estimated by 21 globally downscaled and bias-corrected CMIP5 dataset was published after the conclusion of the present study (Thrasher et al., 2022). Daily average temperature was calculated as the mean of daily minimum and maximum temperatures and the ensemble mean was used to represent future climate. Climate and land cover data from the Global Land Data Assimilation System V 2.0 (GLDAS-2) (Rodell et al., 2004) Dataset was used for parameter estimation and permafrost extent data from the Arctic Permafrost Map (WGS43261) were also used for parameter estimation (Brown et al., 2002). 120

Remotely sensed observations
Binary snow cover classification was derived from the MOD10A1.006 and MYD10A1.006 snow cover products (Hall et al., 2006) (Hall et al., 2016). Data from the 'NDSI_Snow_Cover' band was selected for observations with the highest quality estimate ('NDSI_Snow_Cover_Basic_QA' = 0). The daily mean of 'NDSI_Snow_Cover' band was calculated from both snow cover products. Pixels with 'NDSI_Snow_Cover' > 0 were classified as snow cover (1), and other as no snow (0). The average 135 annual SCF was calculated by counting the number of snow-covered days and dividing by the number of days with valid observations in each pixel, per hydrological year. SCF was calculated based on the highest quality observations, thus excluding lower quality observations as well as missing data due to cloud cover. The availability of MODIS data during polar darkness is temporal limitation for the data set. It limits the application to identify rain on snow events that can cause flooding and deplete areas of snowpack 140

Snow Modelling
Daily snowpack in Iceland was simulated for each hydrological year in the period 1950-2100 using the SNOW-17 model (Anderson, 2006). The model was run in a 0.2-degree resolution with daily average precipitation and temperature data from each of the 21 downscaled and bias corrected CMIP5 GCM's in the NASA NEX GDDP dataset (Thrasher et al., 2012). The model was initialized at the start of each hydrological year in the study period to prevent snow accumulation between years. 145 The model was applied to each of the 21 CMIP5 GCMs in the NASA NEX GDDP dataset and to both the RCP4.5 and RCP8.5 scenarios. These scenarios were chosen to represent both a "business-as-usual" scenario (RCP8.5) and a stabilization scenario (RCP4.5) where anthropogenic climate forcing is assumed to be stabilized by the end of the century. The SNOW-17 algorithm was coded in Google Earth Engine (GEE) as described in Anderson (2006), the simulations were performed in GEE and the input data were accessed through the GEE data catalog. 150 SNOW-17 uses 10 model parameters must be specified by the user for each location. In this study the SNOW-17 parameters were determined at the model resolution across Iceland based on local topology, ecology and hydrology. The recommendations provided by the author of the model (Anderson, 2006) were followed for all model parameters except the melt factors MFMAX and MFMIN which are key model parameters that describe the relation between surface air temperature and snow melt.   (2006)  165 The yearly April 1 st SWE was extracted for each ensemble member. The annual SCF was also calculated for each ensemble member as the number times model grid cell contained snow per year divided by the number of days in that year. The April 1 st SWE was used as it has a long history of use as a snow metric for streamflow forecasting and the SCF was used as it has been suggested as a more appropriate snow metric for a changing climate (Nolin et al., 2021).

Data Analysis 170
This study analyzed a large amount of data on Icelandic snow conditions with the purpose of studying long term trends, to prevent further sources of uncertainty the data products were used as published without additional manipulation. The presence of a statistically significant trend in the time series of in situ observed mean annual SCF, SND was estimated using the Mann-Kendall trend test and by using Sens's estimator of slope method for the MODIS observations. Both of these tests have often been applied to trend analysis in snow cover studies (Notarnicola, 2020;Yilmaz et al., 2019). The Sen's slope method was 175 applied to the remotely sensed observations as it is tolerant to outliers (e.g. Nguyen et al., 2022). The trend test p-values were calculated for the annual SCF and SND timeseries and the trend was assumed statistically significant if p < 0.05, indicating a presence of a monotonic trend within a 95% confidence level. The average annual snow rain ratio in Iceland was estimated from ensemble average of air temperature and surface precipitation data from NEX-GDDP dataset by applying a simple rain/snow partitioning scheme, where precipitation is classified as snow under a set temperature threshold (0°C). Google Earth 180 Engine (GEE) (Gorelick et al., 2016) was used to access data, perform simulations and analyze results. Statistical analysis was performed using GEE and the SciPy toolbox (Oliphant, 2007). ArcMap 10.7.1 was used to produce maps showing the results.

200
The results in Figure 2 show that on average both SND and SCF in Iceland have a positive trend over the period 1930-2021.
The trend is more apparent when considering both fully and patchy snow cover status, (SNC or SNCM ≥ 2) and the data reveal considerable natural climate variability. The MODIS estimates of SCF below and above 500 m a.s.l. are comparable to the insitu estimates of local and mountain SCF, respectively. Iceland is expected to decrease continually, from around 0.6 to around 0.2 and 0.1 for RCP4.5 and RCP8.5 respectively. This 210 trend will be apparent sooner at lower elevations where air temperatures are closer to the snow/rain partitioning threshold. At higher elevations the observed increase in precipitation will result in a temporarily thicker snowpack overall, as air temperatures are further from reaching the threshold, which would offset the increased winter snowmelt and shortening of the snow cover duration associated with temperature rise until the threshold is reached, Iceland. A few areas showed significant decreases in the SCF and most of those were located at the termini of the country's 235 major outlet glaciers, whose retreat has been well documented (Aðalgeirsdóttir et al., 2020;Hannesdóttir et al., 2019;Hauser & Schmitt, 2021) or in coastal areas. All manned observations sites where a decrease in SCF or SND had occurred over the period were all located at low elevation in coastal areas except for one.  Figure 5 shows the results of the simulation of daily snow conditions in Iceland for the period 1950-2100 for both Representative Concentration Pathways (RCP) RCP4.5 and RCP8.5. 5a shows the average winter SWE across Iceland and 5b

Projected Seasonal Snow Conditions 245
shows the simulated average annual SCF along with in situ and MODIS derived SCF estimates.  Figure 5 shows that both SWE and SCF are expected to decrease in Iceland over the course of the 21 st century and that this 255 decrease has been ongoing throughout the study period. Figure 5b shows that MODIS derived SCF estimates over the period of 2001-2021 fit well with the simulated values. In situ observations of local and mountain snow cover status (SNC or SNCM > 2) fall below and above the simulated averages, respectively, as expected. However, although observed and simulated SCF estimates fit well with the magnitude and variability of each other, their trend is opposite. With the observational data consistently trending with increasing SCF while the simulations show a decreasing trend over the historical period, 1950-2020. 260 This pattern of opposing trends is also observed in terms of snow magnitude as the simulated SWE estimates show a decrease in SWE, whereas the observed snow depth measurements (shown in Figure 2) show a significant increase (p = 1.54 * 10 -5 ) over the period 1930-2021. The results also illustrate the substantial natural climate variability in Icelandic snow conditions.
The results in Figure 2 show a positive trend for temperature and precipitation in Iceland over the period 1950-2021. Increasing 265 temperatures result in enhanced snow melt, which is apparent in a flat or decreasing SCF in coastal regions (shown in Figure   4), whereas at higher elevation the increased precipitation enhances winter snow accumulation leading to higher SCF despite the enhanced snow melt during summer, leading to and countrywide increase in average SCF. This effect of increased snow cover at high elevations can be expected to persist until temperatures have risen above freezing for a significant portion of the winter at the highest elevations as well, after which snow-cover is expected to decrease at all elevations. Due to variability in 270 the Icelandic landscape and topography this effect should be more apparent when simulated at a finer spatial resolution.
Recent studies have suggested a regional cooling in the ocean temperatures surrounding Iceland due to changes in the thermohaline circulation in the North Atlantic Ocean (Caesar et al., 2018). This regional cooling, which has been connected to a temporary slowing of glacial ablation in Iceland (Noël et al., 2022), would explain the opposing trends observed in Figure  275 5 as this regional cooling is poorly represented in the CMIP5 models used for the snow simulations in this study. The North-Atlantic cooling trend is projected to halt around 2050 given the results of the Community Earth System Model version 2 (CESM2) (Danabasoglu et al., 2020).

Discussion and Conclusion
The analysis of snow observations showed a significant increase in both Snow Cover Frequency (SCF) and 1 st April SWE, 280 both as estimated from in situ observations over the period 1930-2021 and from observations from the MODIS instruments on NASA's Terra and Aqua satellites over the period 2001-2021. The MODIS observations were comparable with in-situ observations of both local and mountain snow cover status. The results also revealed a large natural variability in snow conditions, which was expected due to the sensitivity of the Icelandic climate to fluctuations in large scale atmospheric and ocean circulations in the north Atlantic region (e.g. Hanna et al., 2004;Massé et al., 2008). The results showed a significant 285 increase in average annual snow depth over all stations for the period 1930-2021.
Simulated Snow Cover Frequency, SCF, was consistent with SCF estimates from both MODIS and in situ observations for the historical period, although the simulated trend was opposite to the trends in both observational datasets. The simulations show that SCF is expected to decrease significantly over the projected period, 2006-2100 especially below 500 m a.s.l. where snow cover is expected to become a rare occurrence by the end of the period, given the RCP8.5 emission scenario. The simulated 290 SWE shows a significant decrease in SWE over the period 1950-2021 whereas average annual SND from all IMO stations has a positive trend over the same period. The results show that the water storage in Icelandic winter snowpack could decrease by about half or 3/4 under the RCP4.5 and RCP8.5 emission scenarios, respectively, over the period 1950-2100.
The results of this study suggest that the increased SCF in Iceland, observed both from remotely sensed and in situ data, is associated with increased precipitation causing a more frequent and thicker snowpack which persist longer, despite enhanced 295 melt rates. This is consistent with Bjornsson et al. (2018) who found annual precipitation to have increased by about 10% during the period 1980-2015. This increasing trend was also observed by Gunnarsson et al. (2019) which used multi-source satellite remote sensing data to show that there had been an increase in snow cover in Iceland for all months except October and November over the period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). The simulated snow conditions are also in agreement with previous projections of a decrease in snow cover frequency and snow mass across Iceland, as rising average temperature causes spring melts to 300 begin earlier and autumn snow cover to occur later (e.g. Johannesson et al., 2007).
The results presented in this study deserve further investigation. Observations of snow conditions reveal a large natural variability which may be affected by large scale circulations in atmospheric and ocean circulations in the northern Atlantic as well as global temperature changes. The observations of both snow cover and snow depth indicate an increasing trend in these variables over the historical period whereas simulated snow conditions predict a decrease in both over the course of the present 305 century, the extent of which is dependent on future emission scenarios. The observed increases in SCF and SWE could be part of natural climate variability induced by low frequency cyclical climate patterns, or by a small amount of extreme weather events. The causes and the impacts of these changes to Icelandic ecology and society should be better understood as future changes to snow conditions will impact the hydrological cycle, which will further affect the local ecology, hazard assessments, water resources management, and hydropower production in the country. 310

Author Contribution
DE and SMG designed the experiments. DE developed the code and performed the analysis and prepared the manuscript. DE and AG gathered, assessed, and prepared the data. SMG, AG and OGBS reviewed the manuscripts and provided significant consultation and contributions throughout the work.

Competing Interests 315
The authors declare that they have no competing interests.

Code/Data Availability
All data used for the analysis in this study are freely available and were accessed either through the Google Earth Engine database or by direct correspondence with the data provider. The datasets used in this study and their source literature and links are provided in Table 3. The code for the snow model and/or the remote sensing analysis can be made available upon 320 request.