The NASA Eulerian Snow on Sea Ice Model (NESOSIM) v1.0: initial model development and analysis

The NASA Eulerian Snow On Sea Ice Model (NESOSIM) is a new, open-source snow budget model that is currently configured to produce daily estimates of the depth and density of snow on sea ice across the Arctic Ocean through the accumulation season. NESOSIM has been developed in a three-dimensional Eulerian framework and includes two (vertical) snow layers and several simple parameterizations (accumulation, wind packing, advection–divergence, blowing snow lost to leads) to represent key sources and sinks of snow on sea ice. The model is forced with daily inputs of snowfall and near-surface winds (from reanalyses), sea ice concentration (from satellite passive microwave data) and sea ice drift (from satellite feature tracking) during the accumulation season (August through April). In this study, we present the NESOSIM formulation, calibration efforts, sensitivity studies and validation efforts across an Arctic Ocean domain (100 km horizontal resolution). The simulated snow depth and density are calibrated with in situ data collected on drifting ice stations during the 1980s. NESOSIM shows strong agreement with the in situ seasonal cycles of snow depth and density, and shows good (moderate) agreement with the regional snow depth (density) distributions. NESOSIM is run for a contemporary period (2000 to 2015), with the results showing strong sensitivity to the reanalysisderived snowfall forcing data, with the Modern-Era Retrospective analysis for Research and Applications (MERRA) and the Japanese Meteorological Agency 55-year reanalysis (JRA-55) forced snow depths generally higher than ERAInterim, and the Arctic System Reanalysis (ASR) generally lower. We also generate and force NESOSIM with a consensus “median” daily snowfall dataset from these reanalyses. The results are compared against snow depth estimates derived from NASA’s Operation IceBridge (OIB) snow radar data from 2009 to 2015, showing moderate–strong correlations and root mean squared errors of ∼ 10 cm depending on the OIB snow depth product analyzed, similar to the comparisons between OIB snow depths and the commonly used modified Warren snow depth climatology. Potential improvements to this initial NESOSIM formulation are discussed in the hopes of improving the accuracy and reliability of these simulated snow depths and densities.

Abstract.The NASA Eulerian Snow On Sea Ice Model (NE-SOSIM) is a new, open-source snow budget model that is currently configured to produce daily estimates of the depth and density of snow on sea ice across the Arctic Ocean through the accumulation season.NESOSIM has been developed in a three-dimensional Eulerian framework and includes two (vertical) snow layers and several simple parameterizations (accumulation, wind packing, advection-divergence, blowing snow lost to leads) to represent key sources and sinks of snow on sea ice.The model is forced with daily inputs of snowfall and near-surface winds (from reanalyses), sea ice concentration (from satellite passive microwave data) and sea ice drift (from satellite feature tracking) during the accumulation season (August through April).In this study, we present the NESOSIM formulation, calibration efforts, sensitivity studies and validation efforts across an Arctic Ocean domain (100 km horizontal resolution).The simulated snow depth and density are calibrated with in situ data collected on drifting ice stations during the 1980s.NESOSIM shows strong agreement with the in situ seasonal cycles of snow depth and density, and shows good (moderate) agreement with the regional snow depth (density) distributions.NE-SOSIM is run for a contemporary period (2000 to 2015), with the results showing strong sensitivity to the reanalysisderived snowfall forcing data, with the Modern-Era Retrospective analysis for Research and Applications (MERRA) and the Japanese Meteorological Agency 55-year reanalysis (JRA-55) forced snow depths generally higher than ERA-Interim, and the Arctic System Reanalysis (ASR) generally lower.We also generate and force NESOSIM with a consensus "median" daily snowfall dataset from these reanalyses.The results are compared against snow depth estimates de-rived from NASA's Operation IceBridge (OIB) snow radar data from 2009 to 2015, showing moderate-strong correlations and root mean squared errors of ∼ 10 cm depending on the OIB snow depth product analyzed, similar to the comparisons between OIB snow depths and the commonly used modified Warren snow depth climatology.Potential improvements to this initial NESOSIM formulation are discussed in the hopes of improving the accuracy and reliability of these simulated snow depths and densities.

Introduction
Snow on sea ice is a crucial component of the polar climate system.Its low thermal conductivity modulates sea ice growth through the cold winter months (e.g., Maykut and Untersteiner, 1971;Sturm et al., 2002), while its high surface albedo limits solar radiation absorption and thus inhibits sea ice melt in spring and summer (e.g., Warren, 1982;Grenfell and Perovich, 1984;Perovich et al., 2002).Conversely, freshwater production from snowmelt facilitates melt pond formation in spring-summer which lowers the surface albedo and promotes sea ice melt (Eicken et al., 2002(Eicken et al., , 2004)).The accumulation of snow on sea ice also modulates the freshwater flux into the ocean, a key component of the freshwater budget of the Arctic (e.g., Serreze et al., 2006).
Estimates of snow depth on sea ice are also a required input for deriving sea ice thickness from satellite altimetry, e.g., from ESA's CryoSat-2 (e.g., Laxon et al., 2013) and NASA's upcoming Ice, Cloud, and land Elevation Satellite (ICESat) mission (Markus et al., 2017).The altimetry technique involves measurements of sea ice freeboard, the extension of Published by Copernicus Publications on behalf of the European Geosciences Union.
sea ice above a local sea level, and estimates of snow depth to derive sea ice thickness, with snow depth being one of the primary sources of uncertainty for both laser and radar altimetry (e.g., Giles et al., 2007).Poor knowledge of snow density provides a further source of uncertainty through its influence on the ice freeboard and radar penetration into the snowpack (e.g., Giles et al., 2007;Kern et al., 2015).
Unfortunately, direct observations of snow depth and density across the polar oceans are very limited, due to difficulties in remotely sensing this relatively thin (O(10 cm)) and heterogeneous medium and logistical challenges associated with in situ data collection.Passive microwave data have been used to estimate snow depth over first-year ice on a basin scale across both poles (e.g., Markus and Cavalieri, 1998;Comiso et al., 2003;Maass et al., 2015), although these data are arguably more relevant for the firstyear-dominated Antarctic sea ice pack and tend to underestimate snow depth in deformed sea ice regimes (e.g., Worby et al., 2008;Brucker and Markus, 2013).Combinations of satellite and/or airborne sensors with variable snow penetration depths are also being explored as a means of producing basin-scale snow depth estimates (e.g., Armitage and Ridout, 2015;Guerreiro et al., 2016;Kwok and Markus, 2017), although this approach is still in its infancy and has limited temporal coverage.NASA's Operation IceBridge has provided airborne measurements of snow depth on sea ice since 2009 (Kurtz et al., 2013).However, the Arctic snow depth data collected are primarily limited to the western Arctic sea ice cover in spring (the spring 2017 campaign also included a flight over the eastern Arctic Ocean), while the Southern Ocean data have only been briefly explored to date (e.g., Kwok and Maksym, 2014).For the Arctic, a climatology of snow depth produced from Soviet drifting station data collected prior to 1991 (Warren et al., 1999) is still commonly used as a basin-scale snow depth product.The Soviet drifting station data also provide the only observationally based basin-scale assessment of snow density currently available.This snow climatology is expected to be outdated due to the rapid changes experienced in the Arctic climate system over the last few decades (Webster et al., 2014), although recent efforts have been made to modify this climatology based on ice type (halving the climatology over first-year ice, e.g., Laxon et al., 2013;Kwok and Cunningham, 2015).
Due to these observational limitations, the sea ice community often utilizes simple models of snow depth forced by reanalyses (primarily snowfall data) (e.g., Maksym and Markus 2008;Kwok and Cunningham, 2008;Blanchard-Wrigglesworth et al., 2018).More sophisticated snow on sea ice models are available, such as SnowModel, a terrestrial snow model recently adapted for sea ice environments (Liston et al., 2018), as well as the prognostic snow layer included in sea ice climate model components, such as the Los Alamos Sea Ice Model (CICE; Hunke and Lipscomb, 2010) and the Louvain-la-Neuve Sea Ice Model (LIM), which have recently undergone various improvements to their snow physics (Holland et al., 2011;Lecomte et al., 2015).
In this study, we present a new model to derive snow depth (and density) across the Arctic Ocean.Our aim is to develop a model of physical and computational simplicity to allow for a detailed assessment of the sensitivity of snow depths to the various input forcing data needed to produce seasonal basin-scale snow depths.The spread in reanalysis-derived snowfall estimates over the Arctic Ocean is high (Boisvert et al., 2018), while the importance and uncertainty of other forcing data (e.g., ice concentration and drift) are still largely unknown.We also wanted a model that could be forced with observed ice concentration and drift to help accurately constrain the seasonal sea ice cycle -a challenge for the more sophisticated sea ice models described above.Our overall goal is that NESOSIM can be used to produce reliable basin-scale daily snow depth and density estimates needed for satellite altimetry calculations of sea ice thickness for both historical analyses and near-real-time operations across the polar oceans.We thus expect the model to increase in complexity with future model developments, e.g., new parameterizations or improvements to existing parameterizations as needed.A secondary utility of the model will be the production of daily/monthly/seasonal snow depths from reanalysis data that can help guide climate modeling research efforts addressing the representation and importance of snow on sea ice in the global climate system.
In the following sections, we present and describe the model configuration/physics, the sensitivity of the model to the input forcing data (e.g., reanalyses snowfall, satellitederived ice drifts) and model calibration/validation efforts.We focus this initial study solely on the Arctic; however, our plan is for the model to be applied and tested in a Southern Ocean framework in the near future.We conclude by looking ahead to potential improvements in the model physics and planned future activities related to our efforts to improve our understanding of snow on sea ice.

Model description
The NASA Eulerian Snow On Sea Ice Model (NESOSIM) is a three-dimensional, two-layer (vertical) Eulerian snow budget model developed with the primary aim of producing daily estimates of snow depth and density across the polar oceans.NESOSIM includes several parameterizations that represent key mechanisms of snow variability through the snow accumulation/growth season and two snow layers to broadly represent the evolution of both old, compacted snow and new, fresh snow.NESOSIM is not currently configured to produce snow depth estimates through the melt season (late spring through summer) due to the lack of surface melt processes included in this initial model formulation.The NESOSIM v1.0 model schematic is shown in Fig. 1.We decided on a Eulerian snow budget approach (as opposed to a Lagrangian approach, e.g., Kwok and Cunningham, 2008) for a number of reasons: (i) it provides a framework flexible to the availability (or lack of) ice drift data, increasing the utility of the model in regions/time periods where ice drift data might be lacking; (ii) it provides a simple assessment of the spatial significance of the parameterized budget terms included in the model, including ice dynamics; and (iii) the parameterizations developed in this framework can be more easily transferred to other Eulerian sea ice models (e.g., the sea ice climate model component CICE) included in general circulation models (GCMs).The following subsections detail the model setup and various parameterizations currently included in NESOSIM.

Model configuration
NESOSIM v1.0 includes two vertical layers on an x-y horizontal grid, with each horizontal grid cell and snow layer featuring a prognostic snow depth and fixed snow density.This two-layer approach was taken to represent the strong differences in properties between dense snow, associated with wind slab, and fresh snow from recent snowfall, while keeping the model computationally efficient and the model physics easily trackable.As stated previously, our plan is that NESOSIM will be used for studying snow on sea ice across the Arctic and Southern oceans; however, for this initial analysis, we run the model on a 100 km × 100 km polar stereographic grid covering the Arctic Ocean and peripheral seas (model domain shown later in Fig. 5).This grid resolution was chosen due to considerations of computational efficiency and the horizontal resolutions of the various input data.The model is forced with daily data of snowfall and near-surface winds from reanalysis data, satellite passive microwave ice concentration and satellite-derived ice drifts.
At each daily time step, an effective snow depth within each grid cell is produced from our various snow budget terms (described in the following subsections) using a forward Euler method as and where t denotes the daily time index, the second index indicates the relevant snow layer (0 = new, 1 = old), and x and y are the horizontal grid indices.NESOSIM v1.0 includes two vertical layers: a "new" layer, ρ n s , which represents recent snowfall, and an "old" layer, ρ o s , which represents snow that has undergone wind compaction and snow grain metamorphism (Colbeck, 1982;Sturm and Massom, 2017).These two fixed snow densities are justified in more detail in the following subsection.We track the evolution of an effective snow depth within each grid cell (the volume of snow per unit grid cell area) for simplicity.The actual snow depth over the ice fraction is calculated by dividing the effective grid-cell snow depth by the grid-cell ice concentration.
We also calculate a bulk snow density, which is the weighted average density across the two snow layers, as Note that the bulk snow density is masked if the respective ice concentration in the given grid cell is less than 15 %, or the effective snow depth is less than 2 cm.While the model tracks the snow budget terms for all grid cell ice concentrations, only grid cells with an ice concentration above 15 % are shown in the analysis, to prevent spurious interpretations in regions of near open water conditions.Each annual model run is initialized in the middle of summer (default of 15 August; rationale discussed in Sect.2.5) and run until the following spring (1 May).This early summer start time was chosen to include the significant snowfall expected across the central Arctic through August (Radionov et al., 1997;Warren et al., 1999;Boisvert et al., 2018) while also avoiding the periods of significant snowmelt in late spring to early summer/midsummer.We acknowledge that this end-of-August time period still likely includes surface melt events that are not captured/included in this model but are hoped to be addressed in future model developments.We also apply a variable initial snow depth (at t = 0) across our model domain, as discussed in Sect.3.4.New ice that subsequently forms in a given grid cell is assumed to be snow free.

Snow accumulation
To accumulate snow on a given grid cell, the snowfall water equivalent from our reanalysis field is converted to snow depth using a representative snow density.Snow pit and density data from the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment and the Soviet drifting ice station data helped guide the parameterization of our seasonal snow density evolution.Initially, snow accumulates into the new, fresh snow layer within a given grid cell as where S f (in units of kg m −2 ) is the gridded daily snowfall across the model domain and A is the gridded daily ice concentration.The density of the new snow layer is fixed at ρ n s = 200 kg m −3 .This value implicitly represents a combination of cold, dry snowfall (∼ 150 kg m −3 ) and wet snowfall (∼ 230 kg m −3 ) based on direct observations over Arctic sea ice (Radionov et al., 1997;Sturm et al., 2002).
Snow can be transferred from the new snow layer to the old snow layer depending on the strength of the near-surface wind forcing.The old snow layer is an implicit combination of two layers that, on average, comprise the majority of the snowpack bulk mass: wind slab and depth hoar (Sturm et al., 2002;Sturm, 2009).The density of wind slab ranges between ∼ 300 kg m −3 and ∼ 410 kg m −3 on average (Colbeck, 1982;Radionov et al., 1997;Warren et al., 1999;Sturm et al., 2002), while depth hoar has an average density of ∼ 150-250 kg m −3 (Colbeck, 1982;Sturm et al., 2002).Based on SHEBA data, the Arctic snow cover consists of slightly more wind slab than depth hoar, comprising ∼ 80 % of it collectively (Sturm et al., 2002).For this reason, we use a weighted average of the higher-end values of wind slab and depth hoar as the density value for the old snow layer, ρ o s .However, we note that the density and ratio of wind slab and depth hoar layers depend on several factors including the atmospheric conditions during precipitation events, sea ice surface roughness, snow depth and the internal snowpack temperature gradient (Sturm et al., 2002).We did experiment with alternative snow densities (e.g., the wider spread of 150 and 400 kg m −3 ) but found this provided worse correspondence with the seasonal snow density evolution compiled from in situ Soviet station data (introduced in Sect.3.4).
When wind speeds are greater than 5 m s −1 , the change in snow depth from wind packing between the two snow layers, respectively, is given as h wp s (t, 0, x, y) = −αT d h s (t, 0, x, y) for U (t, x, y) > ω, (5) where U is the 10 m wind speed, ω is a wind action threshold for wind packing to occur (default of 5 m s −1 ), α is a windpacking coefficient which determines the fraction of the new snow layer that is transferred into the old snow layer (default value of 5.8 × 10 −7 s −1 ), and T d is the number of seconds in our daily time step (equivalent to 86 400 s).The second grid index in Eq. ( 2) (values of 0 and 1) represents the vertical snow layers.The wind action threshold of 5 m s −1 was determined based on observational and modeling studies of blowing snow in the terrestrial Arctic and sea ice environments (Pomeroy et al., 1997;Radionov et al., 1997;Sturm and Stuefer, 2013), while the wind-packing coefficient is a free, unconstrained parameter in the model.

Ice-snow dynamics
Snow within a given grid cell can also evolve due to ice motion.Here, we adapt the ice concentration budget approach used in, e.g., Holland and Kimura (2016) (and more recently in Petty et al., 2018) to snow depth as where u i is the daily gridded ice motion.As in the ice concentration budget studies discussed above, we can expand this into a divergence-convergence term and an advection term as and where h div s is the change in effective snow depth from divergence-convergence, i.e., changes due to spatial gradients in ice motion, and h adv s is the change in snow depth from advection, i.e., changes due to spatial gradients in snow depth (assuming constant drift).Note that this parameter is applied to both "old" and "new" snow layers concurrently.

Blowing snow lost to leads
Snow within a grid cell can also be lost to leads/open water in the ice pack due to the impact of wind forcing, i.e., blowing snow lost to leads.This parameter is expected to be most significant in regions where high lead fractions, wind speeds and snowfall (e.g., the marginal ice zone in the North Atlantic sector of the Arctic) are expected to result in significant wind-blown snow lost to leads/open water (e.g., Leonard and Maksym, 2011).Note that we only apply this wind loss term to the new snow layer as we assume the "old" windpacked snow layer is immune to the impact of wind forcing (e.g., Petrich et al., 2012;Trujillo et al., 2016).The blowing snow to leads term is calculated as where β is a blowing snow coefficient (default value of 2.9 × 10 −7 m −1 ).This is a free, unconstrained parameter in the model, with its default value chosen through our model calibration efforts.We also keep track of snow that enters the ocean through snowfall into the open water fraction and blowing snow lost to leads, a quantity of relevance to those interested in the freshwater budgets of the polar oceans.This is given as For model testing, we also ran NESOSIM with different combinations of the model parameterizations discussed above.When we turn off the wind-packing parameterization, snow remains fixed in the "new" snow layer, despite the strength of the wind forcing, so the model effectively becomes a one-layer model.To account for the low bias in snow density expected by constraining the snow density to the density of fresh, new snow, we forced this snow layer with the daily climatological snow density based on Warren et al. (1999), which we refer to as ρ-W99.
3 Model forcing and calibration/validation data In the following subsections, we describe the forcing data and calibration/validation data used in this study, including atmospheric forcing data (snowfall and winds), satellitederived ice motion, satellite-derived ice concentration, Soviet drifting station snow depths-densities (for model calibration) and Operation IceBridge snow depths (for model validation).

Atmospheric forcing
We use snowfall data provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim (ERA-I) reanalysis.ERA-I is a global reanalysis that utilizes a 4-D variational data assimilation scheme (Dee et al., 2011).We use the 12-hourly ERA-I snowfall data from 15 August 1980 to 1 May 1991 and 15 August 2000 to 1 May 2015.We use the 0.75 • × 0.75 • horizontal resolution data, which are summed to produce daily snowfall estimates across the Arctic.ERA-I snowfall data have been used in previous studies exploring snow accumulation over Arctic sea ice (e.g., Kwok and Cunningham, 2008;Blanchard-Wrigglesworth et al., 2018), while comparisons of reanalysis-derived precipitation data with coastal weather stations suggest ERA-I is one of the better products available for Arctic studies (Serreze and Hurst, 2000;Lindsay et al., 2014).A more detailed comparison of snowfall-precipitation estimates over the Arctic Ocean has recently been carried out alongside this study (Boisvert et al., 2018), which we expect to build on in the future.
We explore the sensitivity of our results to the input snowfall data by forcing the model with snowfall estimates provided by three additional reanalysis-derived snowfall products.Unfortunately, not all reanalyses provide direct estimates of snowfall (and rainfall) and instead provide just total precipitation, e.g., the data from the widely used National Centers for Environmental Prediction (NCEP) -National Center for Atmospheric Research (NCAR) reanalyses 1 and 2, so we focus our analysis on three other commonly used reanalyses that provide direct estimates of snowfall: the Japanese Meteorological Agency 55-year reanalysis (JRA-55); NASA's Modern-Era Retrospective analysis for Research and Applications (MERRA); and the Arctic System Reanalysis, version 1 (ASRv1), as described below and summarized in Table 2.
We MERRA.NASA's MERRA is a global reanalysis that utilizes a 3-D variational data assimilation scheme within the Goddard Earth Observing System (GEOS-5) data assimilation system (Rienecker et al, 2011).We use the daily MERRA snowfall data from 15 August 1980 to 1 May 1991 and 15 August 2000 to 1 May 2015.The data are provided at a horizontal resolution of 0.5 • (latitude) by 0.66 • (longitude).Note that an updated version of MERRA (MERRA-2) is also available, but is known to have a high precipitation bias compared to the other reanalyses (Boisvert et al., 2018), so we exclude this from our study.
ASRv1.ASRv1 is a regional reanalysis based on the Weather Research and Forecasting model (polar WRF) that utilizes a 3-D variational data assimilation scheme and is adapted for the polar regions (Hines and Bromwich, 2008).The ASRv1 data are only available from 2000 to 2012, so we use the daily snowfall data from 15 August 2000 to 1 May 2012, which is provided at a horizontal resolution of 30 km × 30 km.
Considering the expected importance and uncertainty of the reanalysis-derived snowfall for deriving snow depth, we also produce a synthesized snowfall dataset by taking the median snowfall across the gridded snowfall products for each daily grid cell (data referred to as MEDIAN-SF).We use the gridded ERA-I, JRA-55 and MERRA snowfall data, as these products all cover the longer-term  time period.
NESOSIM also requires daily estimates of near-surface winds to drive the wind-packing and wind loss terms, which we take from the ERA-I reanalysis for all reanalysis model runs.Jakobson et al. (2012, Fig. 2) show that ERA-I winds had the lowest biases of several reanalysis-derived nearsurface wind estimates compared to Tara drifting station data.We compute the magnitude of the winds from the 6-hourly u − v vectors before averaging to produce a daily (gridded) wind magnitude dataset.
We linearly interpolate all the daily snowfall (and ERA-I wind magnitude) estimates onto our 100 km × 100 km polar stereographic model domain.Gridding scripts written in Python are included in the NESOSIM GitHub code repository (https://github.com/akpetty/NESOSIM,last update: 8 November 2018).

Satellite-derived ice motion data
We primarily make use of the daily Polar Pathfinder ice motion data, version 3 (Tschudi et al., 2016) made available through the National Snow and Ice Data Center (the product is referred to herein as NSIDCv3).A daily ice motion vector is calculated using a cross-correlation technique applied to sequential daily satellite images acquired by passive microwave satellite sensors (i.e., a 1-day lag in parcel tracking) which are blended via optimal interpolation with estimates from the International Arctic Buoy Programme (IABP) and wind data from the NCEP/NCAR reanalysis.The data are available daily from October 1978 to February 2017 (at the time of writing) at a horizontal resolution of 25 km × 25 km.In this study, we use the daily data from 15 August 1980 to 1 May 1991 and 15 August 2000 to 1 May 2015.We grid the daily ice motion data onto our 100 km model domain (using linear interpolation) and smooth the data using a simple Gaussian filter (as in Holland andKimura, 2016 andPetty et al., 2018).
Recent studies have explored the uncertainty in satellitederived ice motion data (Sumata et al., 2014) and errors introduced by the NSIDC interpolation methodology (Szanyi et al., 2016).We thus also explore the sensitivity of the model results to the input ice motion data by forcing the model with ice motion estimates provided by three additional satellitederived ice motion products, as described below and summarized in Table 3.
Ocean and Sea Ice Satellite Application Facility (OSI SAF).The European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) produce a number of low-resolution sea ice motion products from satellite passive microwave sensors and scatterometry (Lavergne, 2010).Here, we use the merged ice motion product, which increases coverage and reliability over their single-sensor drift products (Lavergne, 2010).The merged drift product uses a 2-day lag in ice parcel tracking and a continuous maximum crosscorrelation (CMMC) method to optimize the drift product and has been available daily (October through April) since 2010 at a horizontal resolution of 62.5 km × 62.5 km.
CERSAT.The Centre ERS d'Archivage et de Traitement (CERSAT), part of the Institut Français de Recherché pour l'Exploitation de la Mer (IFREMER), produces a number of ice motion datasets by merging various combinations of satellite passive microwave and scatterometry data (Girard-Ardhiun and Ezraty, 2012).Here, we use data produced from the merging of Advanced Scatterometer (ASCAT) and the Special Sensor Microwave Imager (SSM/I) data, which have been available daily (September to May) since 2007 at a horizontal resolution of 62.5 km × 62.5 km.Note that CERSAT provides data using both a 3-and a 6-day lag in the tracking of ice displacement, but we use the 3-day lag data as this is closest to the 1-day lag used by the NSIDCv3 product.
Kimura.The Kimura drift data are produced using brightness temperatures obtained by the Advanced Microwave Scanning Radiometer for Earth Observing System (AMSR-E) from January 2003 to September 2011 and the Advanced Microwave Scanning Radiometer 2 (AMSR2) from July 2012 to December 2016 using a cross-correlation approach (see Kimura et al., 2013 for more details).Wintertime (November-December, January-March) ice motion vectors are derived using the 36 GHz channel, while the summertime drifts used in this study (August-October, April) are derived using the 18 GHz channel, to maximize the reliability and coverage of the data.The data are provided at a 60 km × 60 km horizontal resolution.
We use data from these three additional products from 15 August 2010, 2012, 2013 and 2014 to 1 May of the subsequent years, a period of coincident data coverage across the four drift products (including NSIDCv3).We linearly interpolate all the daily drift datasets onto the 100 km × 100 km polar stereographic model domain used in this study.As highlighted above, not all the products produce drift estimates in August, or even September, so for those products we assume no ice motion through this period.To investigate the importance of ice motion, we also run the model assuming no ice motion for the entire model simulation (NODRIFT), as discussed in more detail later.

Sea ice concentration
We use the daily Bootstrap sea ice concentration (SIC) data, version 3 (Comiso, 2000(Comiso, , updated 2017)), which are produced from passive microwave brightness temperature estimates and made available through the NSIDC.We choose to primarily use the Bootstrap over, for example, NASA Team data (Cavalieri et al., 1996(Cavalieri et al., , updated 2017)), another commonly used SIC dataset, as Bootstrap SIC data are less sensitive to surface melt, producing higher concentrations in general.We use the NASA Team data in a sensitivity study to explore the sensitivity of the model to this choice of sea ice concentration data.Due to differences in satellite orbit and sensor characteristics, the SIC data feature a time-varying pole hole depending on the passive microwave sensor used.As we require consistent SIC data across the pole hole, we follow the approach of Petty et al. (2018) and apply a mean SIC calculated in a 0.5 • halo around the variable pole hole to all grid cells within the pole hole.The data are provided at a 25 km×25 km resolution polar stereographic grid from 1978 to 2016, and we use the daily data from 15 August 1980 to 1 May 1991 and 15 August 2000 to 1 May 2015.We linearly interpolate the daily SIC data onto our 100 km × 100 km model domain.Note that a gap in the passive microwave record exists from 3 December 1987 to 13 January 1988, so we do not run the model through the 1987-1988 winter period.

Soviet station data and initial conditions
We use in situ snow data collected on the former Soviet Union's drifting ice stations for initial model calibration and to help guide our choice of initial conditions (Radionov et al., 1997;Warren et al., 1999;Fetterer and Radionov, 2000).The drifting ice stations were in operation in 1937 and 1954-1991, although in this study we use the field observations collected from 1980 to 1991 due to the temporal overlap with the model forcing data.During the drifting ice stations, snow depth data were collected every 10 days in 10 m intervals along a 500 or 1000 m survey line.Snow density measurements were made every ∼ 100 m along the same survey lines, and atmospheric conditions were recorded at near-daily frequencies.Despite their limited spatial coverage, these data provide the most complete record of snow and atmospheric conditions to date over the Arctic sea ice pack.
Initial conditions.We initialize the model on 15 August of each year with a snow depth representing the fraction of snow assumed to have survived the summer melt season and/or accumulated during summer.The August snow depth climatology compiled by Warren et al. (1999, referred to herein as W99) from the Soviet station data suggests significant amounts of snow (up to 10 cm) are present in late summer, especially over the central Arctic sea ice north of Greenland (Radionov et al., 1997).This inclusion of an initial snow depth was also guided by our preliminary model calibration studies that showed that including these initial conditions provided a better match with the seasonal snow depth observations (calibrations presented later).To produce initial mid-August snow depths, we use a near-surface air-temperaturebased scaling of the August W99 snow depth climatology to account for changes in the duration of the summer melt season (e.g., Markus et al., 2009).Briefly, we calculate the annual number of days with continuous, above-freezing air temperatures (taken from the ERA-I reanalysis), which we refer to here as the summer melt duration.To create an initial (August) snow depth estimate for a given year, we linearly scale the W99 August snow depth climatology based on the summer melt duration of the chosen year and the climatological summer melt duration given in Radionov et al. (1997).If the melt duration is longer than the climatological mean in a specific region, the scaled August climatology reflects a reduction in snow depth in August due to the longer melt season.The snow depth is then distributed evenly over the "old" and "new" snow layers based on the climatological observations that some snow persists through summer (Radionov et al., 1997) and the occurrence of summer snowfall events (Radionov et al., 1997;Perovich et al., 2017).While admittedly this is a crude approach for parameterizing an initial snow depth, our sensitivity studies demonstrated that initial conditions were necessary to improve the comparison with the drifting station observations (as presented and discussed in the following section) and indicate that late summer snowfall events might play a significant role in establishing the snow cover on Arctic sea ice prior to the fall-winter season (Warren et al., 1999).The August W99 snow depth climatology and temperature-scaled initial snow depth estimates (for 2012 and 2013) are shown in Fig. 2.

NASA's Operation IceBridge data
We compare our NESOSIM snow depth estimates with spring snow depths collected by NASA's Operation Ice-Bridge (OIB) airborne mission.NASA's OIB mission began collecting airborne observations of the polar regions in 2009, bridging the gap between NASA's ICESat mission, which retired in 2009, and the future ICESat-2 mission scheduled for launch in the summer of 2018 (Markus et al., 2017).The OIB aircraft carry a suite of instruments designed to measure both land and sea ice, including their overlying snow cover.Here, we primarily make use of snow depth estimates derived from the ultra-wideband snow radar (Panzer et al., 2013), which have a footprint size over sea ice of 5-10 m.These snow depths are thought to carry an uncertainty of several centimeters, although this depends strongly on the ice-snow conditions, the particular snow radar system being used and various other factors, e.g., geolocation errors associated with the plane pitch and roll (e.g., Kurtz et al., 2013;Kwok et al., 2017 and references therein).Various algorithms have been developed to produce snow depth estimates from the OIB snow radar data (Kwok et al., 2017), with the products showing broad agreement in the regional snow depth distributions, but significant intraregional and interannual differences, due primarily to changes in the radar configuration and algorithm tuning.To account for these differences, we use the snow depth data from the (i) snow radar layer detection (SRLD) (Koenig et al., 2016), (ii) NASA Goddard Space Flight Center (GSFC) (Kurtz et al., 2013) and (iii) Jet Propulsion Laboratory (JPL) (Kwok and Maksym, 2014;Kwok et al., 2017) at the raw 5-10 m snow radar resolution from 2009 to 2015.We bin the raw OIB snow depth data onto our 100 km model grid and keep only the grid cells that included a significant quantity (> 1000 points) of the raw snow depth data.The OIB data are provided for various days through spring of the relevant campaign (data from mid-March to early May, depending on the campaign year), so we grid the OIB data daily and compare this with coincident (daily) NESOSIM snow depth estimates.The OIB data are collected mainly over the western Arctic sea ice, limiting our validation efforts to this region of the Arctic.

Model calibration and analysis
We carried out model calibration over the period 15 August 1980-1 May 1991 due to the coincident Soviet station data available during this period.As noted previously, this excludes the 1987-1988 winter season due to the lack of complete sea ice concentration data available during this period.
As stated earlier, our initial calibration efforts involved manually tuning NESOSIM to improve the general fit with the mean seasonal snow depth-density cycles shown in the Soviet station data.Specifically, we included the temperaturescaled initial August snow depths and tuned both the windpacking coefficient, α (Eq.5), and blowing snow coefficient, β (Eq.6).We decided against a more optimized calibration effort due to limitations in the calibration data, i.e., its sparse availability in space-time and differences in spatial scales.We instead used the Soviet station data to guide our model choices to achieve a more realistic seasonal cycle in snow depth and density.We also decided against specific model configuration parameter tuning due to these limitations in the calibration data; however, this should be con-sidered when analyzing the model performance, especially with regard to our validation efforts (i.e., more sophisticated and/or configuration-specific tuning could improve the comparisons shown).
In Fig. 3, we show comparisons of our NESOSIM results using the default model configuration (summarized in Table 1) and ERA-I snowfall forcing with the drifting station snow depth and density data.Figure 3 shows both the mean seasonal cycle based on all drifting station data points and coincident model grid cell values over this time period binned monthly and the correlations of snow depth and snow density for all coincident data (described in Sect.3.4).Our calibrated NESOSIM results agree well with the mean seasonal cycle in snow depth (r = 0.96 with a low bias of ∼ 3 to 7 cm) and snow density (r = 0.97, no significant seasonal bias) in the drifting station data.The large spread in the in situ snow density in September-October is due to the survival of snow through the summer melt season (high density) and recent autumn snowfall (low density).The correlations between the raw drifting station data and NESOSIM snow depths are lower, but still strong (r = 0.74), while the snow www.geosci-model-dev.net/11/4577/2018/Geosci.Model Dev., 11, 4577-4602, 2018 NO indicates that the parameterization/initial conditions have been turned off.The different colors then represent a doubling of individual model parameters, with all other settings fixed to the default settings (see Table 1).The black crosses and line represent the default/ERA-I results (as shown in Fig. 3).
density correlation strength is moderate (r = 0.58), suggesting the model may be better capturing regional variability in snow depth over snow density.It should be noted, however, that snow density is highly variable in space and subject to large measurement uncertainties when collected in situ (Sturm, 2009).In general, the moderate-high correlations and seasonal comparisons provide confidence in the utility of NESOSIM for estimating snow depths across the Arctic.In Fig. 4, we highlight the sensitivity of NESOSIM to the chosen model configuration/sophistication, broadly representing the heuristic model tuning that was undertaken.First, we tested the results of NESOSIM with different combinations of the various model parameterizations included.Note that, as discussed at the end of Sect.2, when we turn off the wind-packing parameterization, the model essentially becomes a one-layer model, so we use a fixed Warren et al. (1999) seasonal snow density climatology (constant density value across the Arctic).As this is based on the same drifting station data we compare our results to, it is perhaps unsurprising that this configuration provides a better match with the seasonal drifting station snow depth cycle, including deeper snow depths (and reduced low snow depth bias) from November to April.We chose to develop NESOSIM to allow for the production of snow depths that agree well with the old drifting station snow climatology but are able to also respond to the expected interannual variability and trends in Arctic climate over recent decades -hence the decision to develop and include a simple bulk density parameterization.
Including the blowing snow loss parameterization resulted in slightly lower snow depths (∼ 2 cm), but no significant change in snow density.This parameterization can impact the bulk density implicitly by reducing the amount of fresh snow contributing to the total snow depth-density.As the drifting station data are collected primarily within the central Arctic where ice concentrations are close to 100 %, it was expected that including blowing snow loss would not result in significant differences, as this parameterization is expected to provide more of an impact in lower ice concentration regimes, where unfortunately in situ snow depth data are lacking.Including the initial snow depths resulted in a small increase in snow depth and density, especially earlier in the seasonal cycle, as expected, reducing the low bias compared to the drifting station data.The seasonal correlations were similarly high across these model configurations, highlighting the primary role of the model configuration choices in determining the general bias of the seasonal snow depth-density cycle.
As a simple demonstration of the sensitivity of the model to the poorly constrained-unconstrained model parameters introduced in NESOSIM (the wind-packing threshold, ω, the wind-packing coefficient, α, the blowing snow loss coefficient, β), Fig. 4 also shows results from NESOSIM with these three model parameters individually doubled (based on the default/ERA-I configuration).Doubling the windpacking threshold, ω, (from 5 to 10 m s −1 ) has a large impact on both the snow depth and density.By essentially reducing the likelihood for wind packing to occur, the snow  accumulates and remains in the fresher "new" snow layer for longer, significantly reducing the bulk snow density and increasing the seasonal snow depths.While this does produce snow depths that appear to agree better with the drifting station data, the low bias in the seasonal snow density suggests this is unphysical.Doubling the wind-packing coefficient, α (from 5.8 × 10 −7 to 1.16 × 10 −6 ), has broadly the opposite effect, as expected, reducing the snow depths by increasing the transfer of snow from the fresher "new" snow layer to the denser "old" snow layer.Doubling the blowing snow loss coefficient, β (from 2.9 × 10 −7 to 5.8 × 10 −7 ), has a negligible impact, again likely due to the location of the in situ data away from the lower concentration ice regimes where this process is more significant.
As stated earlier, the differences in spatial scales and data coverage (time and space) make interpreting these comparisons/calibrations challenging.Specific model configurations may be required based on user demands, and our expectation is for these calibrations to evolve as new calibration data are made available and physical parameterizations introduced/updated.Note that we also compared the simulations of NESOSIM forced by the MERRA and JRA-55 snowfall data (Figs.S2 and S3).In general, the seasonal correlations with the drifting station data were similar to the ERA-I re-sults, but the correlations of the raw data were slightly lower for JRA-55 (r = 0.69 for snow depth and r = 0.58 for snow density) and significantly lower for MERRA (r = 0.44 for snow depth and r = 0.57 for snow density).As discussed in Sect.3, it is likely that specific model configuration tuning could improve these comparisons and the later validation efforts, but we decided against a more optimized calibration approach due to the limitations in the Soviet station data.
As discussed in Sect.3.1, we also produced a synthesis snowfall dataset (MEDIAN-SF) using the median snowfall across the gridded ERA-I, JRA-55 and MERRA datasets.The MEDIAN-SF forced results are similar to the ERA-I results (Fig. S4), in general, and show correlations similar to ERA-I and JRA-55 (r = 0.68 for snow depth and r = 0.58 for snow density).The MEDIAN-SF seasonal snow depths have a reduced low bias compared to the ERA-I results, although this difference is small.For the rest of this analysis, we choose to mainly focus on the MEDIAN-SF forced results using the default configuration (Table 1) for simplicity.We provide a further assessment of the impact of the snowfall data in the following regional analysis and when we analyze the regional distributions across the more recent (2000-2015) time period.

Sensitivity studies and model validation
Here, we present and analyze the NESOSIM results from 2000 to 2015, a period broadly defined as the New Arctic, considering the rapid sea ice declines during this time period (e.g., Serreze and Stroeve, 2015).This period also covers the temporal range of NASA's ICESat (2003 to 2008) and ESA's CryoSat-2 (2010 onwards) satellite altimetry missions, meaning the snow depth-density results presented here are planned to be of more relevance for those estimating sea ice thickness from these freeboard measurements.The period also includes temporal overlap with the ASR forcing data and various satellite-derived ice motion products used.We provide examples of the model evaluation figures for the 1980s time period in the Supplement (Figs.S5 and S6).A more detailed study accounting for differences in the input forcing data is likely needed before any conclusions can be made regarding potential trends in seasonal Arctic snow depths, which is beyond the scope of this paper.We hope to explore trends in our simulated snow depths in future work, however.
We focus our analysis on the Arctic Ocean (AO, everything north of 60 • N) and three specific regions that were chosen to represent different components of the Arctic sea ice-climate system: (i) the central Arctic (CA, captures the thicker/multi-year ice over the north of Greenland), (ii) the eastern Arctic (EA, the increasingly first-year-ice-dominated sea ice regime) and (iii) the North Atlantic (NA, a region influenced by the transpolar ice drift and the North Atlantic storm track), as shown in Fig. 5.    Figure 6 shows the seasonal snow depth and density evolution across our four study regions for the 2000-2015 time period using the default/MEDIAN-SF configuration (Table 1).The AO and CA regions especially show strong initial increases in snow depth through fall (August to October) with the snow depth increasing at a slower rate from November to May, which is in good agreement with the W99 climatology.The EA and NA regions show a more uniform increase in snow depth from August to April.The NA region shows more daily snow depth variability, which was expected due to the strong ice drifts and the location of the NA storm track where passing cyclones can deposit large quantities of snow in a short period of time.It is also worth noting the small decline in snow depth through September-October in the EA region which is driven by reduced snowfall and snow densification due to wind packing through this period.By 1 May, the mean snow depths (and interannual variability, calculated as 1 standard deviation of the annual values) are given as 27.8 ± 1.9 cm (AO), 31.8 ± 4.0 cm (CA), 23.2 ± 2.9 cm (EA) and 42.5 ± 8.1 cm (NA).The 1 May snow depth results are summarized in Table 4 to aid comparison with the snow depths produced in the following sensitivity studies.
We see stronger increases in the bulk snow density through fall across all regions (also shown in Fig. 4), with this density increase slowing through winter-spring, especially in the CA region, after December.The AO, CA and NA regions also show an interesting initial decrease in snow density, which is driven by the accumulation of new snow (with a Figure 6.Seasonal snow depth (black) and bulk density (green) evolution across the four study regions (shown in Fig. 5) initiated from 15 August 2000-2014 and run until 1 May of the following year using the MEDIAN-SF/parameter settings (Table 1).The thick lines show the mean values over this time period, while the shaded areas represent the interannual variability (1 standard deviation).
lower density) compared to the equal mix of old and new snow densities included in our initial conditions.The mean bulk snow densities as of 1 May are given as 309 ± 2 kg m −3 (AO), 323 ± 4 kg m −3 (CA), 311 ± 3 kg m −3 (EA) and 318 ± 6 kg m −3 (NA).

Budget analysis
Here, we discuss the relative contributions to the seasonal snow depth evolution from the various snow budget terms currently included in NESOSIM.Results of the various NE-SOSIM budget terms and the total snow depth and bulk density are shown in Fig. 7 across our four study regions for this 2000s time period.The black (green) lines and/or shading that represent the snow depth (bulk density) are the same as the results shown in Fig. 6.Across the AO region, we see that accumulation is higher than snow depth, as expected (higher by ∼ 30 cm by 1 May, around double the 1 May snow depth), with wind-packing (∼ 20 cm) and wind-blowing snow lost to leads (∼ 10 cm), providing significant reductions in snow depth.In the EA and CA regions especially, the blowing snow loss term is negligible, while in the NA region it is more significant (contributes a sink of ∼ 18 cm by 1 May).The NA region also shows a small (∼ 2 cm) increase (decrease) in snow depth driven by snow-ice divergence (ice-snow advection).
To further explore the different budget terms, we also show maps of the various budget terms as of 1 May over the same time period, as shown in Fig. 8.The maps highlight that many of these terms, especially the ice-snow dynamics (advection and convergence), exhibit high spatial variability, which the regional means discussed previously mask.For example, the NA region shows a strong mix of positive snow advection and convergence adjacent to the coast of Svalbard (i.e., snow is advected into the region and is constrained against the coastline) but an advection out of the region further to the north as the ice either drifts down towards Svalbard/the Fram Strait or into the central Arctic.
The ice dynamic behavior around the pole is thought to be spurious considering interpolating issues across the pole hole in the NSIDCv3 drift product (Szanyi et al., 2016), which is one reason why we did not include this region in our regional analysis.In the following section, we assess the sensitivity of our results to the input ice motion dataset, which will provide some further information as to the reliability of these dynamic budget terms.
As stated earlier, we hope to explore these decadal and regional differences more in future work.However, it is worth noting that the regional snow budget results and 1 May budget maps using the same default/MEDIAN-SF configuration but run for the 1980s time period were similar to the 2000s time period results (Figs.S4 and S5).The noteworthy differences in the budget terms include a less significant increase in blowing snow lost to leads in the CA region and less convergent driven snow depth increases in the new period, although accumulation and wind packing still dominate the budget terms for both periods.The NA results also do not show the advection-driven reduction in snow depth in March-April that was present in our 2000s results.

Reanalysis sensitivity study
In Fig. 9, we show the seasonal-regional snow depths from NESOSIM forced by the various reanalysis-derived snowfall estimates (ERA-I, JRA-55, MERRA and MEDIAN-SF) from 2000 to 2015 and the ASRv1 forced results which are only available up to 2012, as described in Sect.3.1.The 1 May results are summarized in   sults show significant differences in the seasonal snow depths across all regions (up to ∼ 10 cm across all regions).The rankings of snow depth between the different products are broadly consistent across the four regions, with JRA-55 and MERRA producing consistently higher snow depths (except in the EA region, where MERRA produces slightly higher snow depths) and ERA-I consistently lower.The MEDIAN-SF snow depths are, in general, slightly higher than the ERA-I forced snow depths.In the CA region, we can see that MERRA, ERA-I and MEDIAN-SF forced results are all broadly similar, with JRA-55 significantly higher (by ∼ 5 cm from October onwards).It is thus expected that the MEDIAN-SF snowfall data will have excluded much of the high JRA-55 snowfall data (the benefits of using a median instead of a mean snowfall).Despite the NA region having the highest snow depths and interannual variability, the intrareanalysis spread is similar to the other regions.The ASRv1 forced snow depths in the AO, CA and EA regions are signifi-cantly lower during the December-April time period, despite showing strong similarities to the other reanalysis-forced results in August-November.The ASRv1 results in the NA region, however, are very similar to the ERA-I forced results.Note that we tested the impact of the different time periods by producing the same figure for the 2000-2012 period (not shown), which showed that the differences between ASRv1 and the other products were similar and not sensitive A. A. Petty et al.: NESOSIM to this time period difference.The results further allude to the need for a consensus (e.g., median) snowfall dataset to force the model with the consideration of the large uncertainty in reanalysis-derived snowfall.
Figure 10 shows maps of the mean snow depths on 1 May over the same 2001-2015 time period, for the model simulations forced by the MEDIAN-SF snowfall, then the differences from this MEDIAN-SF simulation using the four individual snowfall products.
The maps highlight the regional variability across the products but consistency in the MERRA/JRA-55 (ASR) high (low) difference compared to MEDIAN-SF.The JRA-55 and MERRA forced results both show significantly higher (10-20 cm) snow depths through Bering Strait, the NA/Fram Strait region and the southern Labrador Sea.The ERA-I results show slightly lower snow depths over most of the Arctic, small increases around the Canadian Arctic Archipelago and larger decreases in the Fram and Bering Strait region, driven by the larger differences in these regions in the MERRA/JRA-55 forcings.The magnitude of the precipitation events in Fram Strait are often large, but highly variable, due to the active storm track and the resulting difficulties of producing reliable precipitation rates during these events (Boisvert et al., 2018).As discussed earlier, it is challenging to determine from this study any particular reanalysisderived snowfall dataset that might be more appropriate (or an obvious outlier) for producing accurate snow depth estimates across the Arctic.However, the MEDIAN-SF forced results appear to provide a useful synthesis of the available snowfall data.

Ice motion sensitivity
We also explore the sensitivity of NESOSIM to the input satellite-derived ice motion data available during this period.Here, we show results from the default/MEDIAN-SF configuration forced by four different satellite-derived ice drift products: NSIDCv3, Kimura, CERSAT and OSI SAF, as described in Sect.3.2.Due to limitations in the temporal coverage of the different drift datasets, the model is only run for 4 years initialized from 15 August 2011 to 2015 (excluding 2012 initialized runs as Kimura data are not available due to gaps in the AMSR-E/AMSR2 record).The regional snow depth estimates from NESOSIM forced by these four ice drift products are shown in Fig. 11, with the 1 May results summarized in Table 4.In general, the ice drift sensitivity study shows a smaller spread in the mean snow depths across the different products (up to ∼ 3.5 cm), compared with the reanalysis sensitivity study (up to ∼ 13 cm).We also show results of NESOSIM forced with no ice drift (NODRIFT), which demonstrates that including ice drift appears not to be a crucial process for capturing the variability in snow depth at this regional scale; i.e., ice dynamics appear to be a clear second-order term compared to snowfall when analyzed at this regional scale.
The most obvious impact of ice drift is in the EA and NA regions.In the EA region, the inclusion of ice drift reduces the snow depth by 1.4-3.6 cm, with the magnitude depending on the ice drift product chosen (the Kimura forced results show the biggest decrease in this region).In the NA region, the inclusion of ice drift increases the snow depth by 3.8 to 5.3 cm (the NSIDCv3 forced results show the biggest increase in this region).
Figure 12 shows maps of the snow depths averaged on 1 May over the same 2011-2015 time period, for the model simulations assuming no drift (NODRIFT), then the differences from this NODRIFT simulation using the various ice motion products.In general, the results show strong similarity in the spatial impacts of ice motion, including strong decreases in snow depth (up to ∼ 10 cm) in the northeastern sector of the Arctic, and increases (up to ∼ 10-20 cm) in the region directly north and west of Svalbard.There are clear differences between the different ice motion results, though, with the NSIDCv3 and Kimura forced results showing more of an impact on snow depth in the peripheral Arctic regions, e.g., strong decreases in the north and increases in the south Bering Strait, and strong increases in the Labrador Sea.This is thought to be driven primarily by the increased spatial coverage of these data compared to OSI SAF and CERSAT, which may be masking some of the ice motion data in these regions of low ice concentration and uncertain ice drift.The maps also highlight that, at more local scales, the ice dynamic contribution to snow depth variability could be significant.The data around the pole hole are also questionable in some of the products and may be related to interpolation issues across the pole hole.More specifically, the NSIDCv3 and OSI SAF forced simulations show increases in snow depth at the North Pole, which are not apparent in the CERSAT and Kimura simulations, suggesting this increase is likely spurious.
In general, Figs.11 and 12 suggest that the NSIDCv3 (Polar Pathfinder) forced simulations exhibit no obvious biases compared to the results using the other ice motion products, except for the issues of spurious snow depths within the pole hole and issues around the ice edge.

Ice concentration sensitivity
Finally, we present and discuss the snow depth results from NESOSIM driven by two different satellite-derived ice concentration products (Bootstrap and NASA Team), as described in Sect.3.3.The regional snow depth estimates from NESOSIM forced by these two ice concentration products over the 2000-2015 time period are shown in Fig. 13, with the 1 May results summarized in Table 4.In general, the ice concentration sensitivity study demonstrates that the choice of ice concentration product is significant, with differences of several centimeters between the two simulations across the study regions (e.g., ∼ 7 cm differences in the NA snow depths).5), initiated from 15 August 2000 to 2014 and run until 1 May of the following year, forced by the Bootstrap (magenta) and NASA Team (blue) ice concentration datasets.The thick lines show the mean (daily) regional snow depths over this time period, while the shaded areas represent interannual variability (1 standard deviation).All model runs use the default/MEDIAN-SF configuration.
This was somewhat expected given the known low bias in the NASA Team concentration data (e.g., Meier, 2005;Ivanova et al., 2015), reducing the concentration of sea ice for snow to accumulate on.More specifically, the Bootstrap data use daily variable tie points and are thus thought to improve the distinction between surface melt and open water.The lower concentrations also increase the blowing snow lost to leads term (as this is a function of the open water fraction).The snow budget terms using the NASA Team concentration data are shown in the Supplement (Fig. S7) to highlight this further, with all regions showing reduced snow accumulation and blowing snow lost to leads increased, and now significant, across all regions.Again, we believe the Bootstrap data better represent the seasonal ice conditions, although we appreciate that uncertainties still remain regarding the treatment of surface melt-melt ponds and their impact on snow accumulation/depth.

Validation with Operation IceBridge data
Here, we present and discuss comparisons of our NESOSIM snow depth estimates with NASA's Operation IceBridge spring snow depth data from 2009 to 2015, as described in Sect.3.5.We first show the basin-averaged results for the various OIB snow depth products each spring (from 2009 to 2015) and the coincident NESOSIM snow depth estimates, to assess how well NESOSIM captures the mean snow depth and expected interannual snow depth variability across this broad region of the Arctic.As discussed in Sect.3.5, the OIB flights mainly cover the western Arctic sea ice pack, broadly within and to the west of the central Arctic domain used in our earlier regional analyses, although this does vary each year.Maps of the OIB snow depth retrievals across the different products are given in Kwok et al. (2017).Figure 14 highlights the significant and variable spread in the annual mean OIB snow depth estimates (product spread of ∼ 5 to 20 cm depending on the year), with the OIB-JPL snow depth retrievals consistently higher and less variable than the other two OIB products (SRLD and GSFC).The reanalysis-forced NESOSIM snow depths exhibit a more consistent spread of ∼ 5 cm, with the JRA-55 forced results consistently higher than the other reanalyses.This was expected based on our previous analyses (e.g., the central Arctic results shown in Fig. 11b).The large spread in the OIB snow depths makes it challenging to assess the reliability and accuracy of our NESOSIM results.In general, however, there is broad agreement between the NESOSIM and OIB results www.geosci-model-dev.net/11/4577/2018/Geosci.Model Dev., 11, 4577-4602, 2018 in terms of the mean snow depths and the broad pattern of interannual variability.
To assess how well the model captures regional snow depth variability, we show scatter plots in Fig. 15 of the NESOSIM/MEDIAN-SF snow depths and the three OIB snow depth products from 2009 to 2015 (the regressions for each year are given in Fig. S8).A summary of the correlation coefficients (r) and root mean squared errors (RMSEs) across the three OIB products and NESOSIM forced by the three individual (and median) reanalysis products for individual years and for all years of data is given in Table 5, with the regressions shown in Figs.S9-S11.
In general, the comparisons are highly variable and depend mainly on the chosen analysis year and the reanalysis snowfall dataset, rather than the OIB product.The correlations between the OIB snow depth retrievals and the NESOSIM snow depths improve significantly in 2012 (r = 0.63-0.75)compared to the proceeding years (r = −0.15-0.61).The improved correlations in 2012 onwards coincide with increases in the OIB flight coverage, which include more of the central Arctic and Beaufort and/or Chukchi seas, meaning the data better represent the regional variability in snow depths across the western Arctic.The strength of the correlations is highest in 2012 and 2013, while the RMSEs are lowest (< 10 cm) between 2011 and 2013, especially in the ERA-I and MEDIAN-SF forced results.The OIB-SRLD RMSEs are generally lower than the RMSEs calculated with the other OIB products between 2010 and 2015, but significantly higher in 2009 when the signal-to-noise ratio of the earlier version of the snow radar used on OIB was higher (Kwok et al., 2017).The 2009 OIB snow depth results should thus be treated with caution.
The "all years" results in Table 5 provide a summary of the correlations using all the OIB snow depth retrievals from 2009 to 2015.The MERRA forced results produce significantly lower correlations to the OIB snow depths (r = 0.37 to 0.43) compared to the other reanalyses, while the ERA-I forced results show the highest correlations (r = 0.53 to 0.64) and lowest RMSEs (9-10 cm).The MEDIAN-SF results show slightly lower correlations (r = 0.47-0.58)and higher RMSEs (10-11 cm) compared to ERA-I.In general, however, the moderate-to-strong correlations give us confidence that NESOSIM is producing reasonable snow depth estimates across the western Arctic.The RMSEs of ∼ 10 cm imply the expected level of accuracy in our NESOSIM snow depths, although these validations are hindered by uncertainty in the OIB snow depth retrievals (Kwok et al., 2017) and a lack of OIB retrievals in the eastern Arctic Ocean.
In the sensitivity studies presented earlier, we focused primarily on the MEDIAN-SF simulations, due in part to considerations of high snowfall variability in regions of high and uncertain precipitation, e.g., the North Atlantic sector.The OIB data lack coverage in this region, however, making it hard to assess if this synthesized forcing snowfall produces more accurate snow depths in these more challenging regions of the Arctic.Data from the 2017 OIB flights into the eastern Arctic will hopefully provide some assessment of our NE-SOSIM snow depths in this region of the Arctic, however (the data were not available for this study but were made available during the review phase of this paper).Our contemporary (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) NESOSIM results still lack validation of the simulated snow densities due to the lack of basin-scale density data available during this time period.
We can further assess the performance of NESOSIM by comparing these results with comparisons of OIB and the commonly used Warren snow depth climatology (Warren et al., 1999).As stated in the introduction, more recent uses of this climatology tend to apply a scaling factor (usually 50 %) to the snow depths over first-year ice.We follow the same approach here, using the EUMETSAT OSI SAF (www.osi-saf.org)ice-type product which is derived from a combination of passive microwave and scatterometry data at 10 km horizontal resolution (Breivik et al., 2012).We derive daily modified W99 snow depths (referred to herein as mW99) for the same 100 km bins used in Figs. 14 and 15 (where we have OIB data), with the comparisons shown in Fig. 16.
In general, these comparisons are similar, although in some cases the mW99 snow depths compare better with the OIB snow depths, depending on the product analyzed.Note that the bimodal nature of the mW99 data is due to the binary ice-type weighting scheme, which does significantly improve the comparison to the OIB data based on the RMSE and correlation coefficients (comparisons of unmodified W99 data are given in Fig. S12).The low RMSE values in the JPL-OIB comparison are driven by the very good agreement in the mean snow depth, while the GSFC and SRLD products tend to show a slight low bias, increasing the RMSE values.Figure 14 shows that the OIB-JPL product exhibits less interannual variability than the other products, which may provide some explanation for the better correspondence with the W99 climatology.We also carried out OIB comparisons by delineating by ice type (first-year ice and multi-year ice) using the same OSI SAF product discussed above.However, the results were mixed and were also strongly dependent on the OIB product analyzed.Such delineations are also hindered by the lower coverage of first-year ice in the OIB data, despite this becoming an increasingly dominant component of the Arctic sea ice pack.We thus chose to exclude this analysis from our discussion for simplicity.
In general, the modified Warren results are useful for placing the NESOSIM-OIB comparisons in context, which clearly show a higher spread and tend to suffer from positive and negative biases depending on the OIB product chosen.Bias correcting the NESOSIM snow depths could improve these comparisons, but uncertainty still remains regarding which OIB product better represents "truth" at this basin scale.Table 5. Correlation coefficient (r, top rows) and root mean squared error (RMSE, bottom rows) from the correlations between the various reanalysis-forced NESOSIM results and OIB-derived snow depths.The MEDIAN-SF scatter plots for all years of data are shown in Fig. 15, with other reanalysis-forced scatter plots given in the Supplement.model produced very strong agreement with the seasonal cycles of snow depth and density and good (moderate) agreement with the regional snow depth (density) distribution.
The model was run for a contemporary period (2000 to 2015) to produce seasonal snow depth and density estimates representative of the New Arctic climate system.A budget analysis provided insight into the relative processes contributing to our modeled seasonal evolution in snow depth, with snow accumulation driving increases in snow depth and wind packing reducing snow depth (through an increase in the bulk snow density).Blowing snow lost to leads provided a significant sink of snow, but only in the lower ice concentration, high wind-snow depth regime of the North Atlantic sector.
The model showed strong sensitivity to the reanalysisderived snowfall forcing data, with the MERRA/JRA-55 (ASR) derived snow depths generally higher (lower) than ERA-I.We derived a new synthesized snowfall dataset based on the median ERA-I, MERRA and JRA-55 snowfall data.We briefly assessed the sensitivity of NESOSIM to the input concentration data, with our results suggesting that the choice of concentration product (Bootstrap and NASA Team explored in this study) can have a significant impact on the simulations and should not be overlooked.We also explored the sensitivity of NESOSIM to the input ice motion data, where we showed this had a second-order effect compared to the choice of reanalysis snowfall forcing in our regional mean comparisons.The ice motion still appears to be important at smaller spatial scales, e.g., by reducing snow depths in the eastern Arctic and driving higher snow depths north of Svalbard and within the Fram Strait.
We compared our NESOSIM snow depths against spring snow depths derived from data collected by NASA's OIB since 2009 (up to spring of 2015) from three different algorithms.Our comparisons show moderate-strong correlations for the data collected from 2012 to 2015, but weaker correlations before this.The root mean squared differences were around 10 cm, but depend on the year analyzed, snowfall forcing and the OIB product analyzed.The ERA-I and MEDIAN-SF forced results showing the best correspondence with the OIB snow depths.These results were compared with comparisons between OIB and the modified Warren snow depth climatology, which showed similar correlations and root mean squared errors.
Errors in snow depths of around 10 cm are thought to contribute to errors in ice thickness estimates derived from laser (radar) satellite altimetry of ∼ 70 cm (50 cm) assuming typical freeboards of ∼ 30 cm (Giles et al., 2007).We expect that further model development, calibration, and validation are needed to improve accuracy and reliability in the NE-SOSIM snow depths-densities, to improve their utility in ice thickness retrieval analyses.

Future work
This initial formulation of NESOSIM (v1.0) has focused on (i) incorporating several key snow parameterizations needed to capture the regional and seasonal variability in snow depth and density across the Arctic Ocean and (ii) providing a framework simple and computationally efficient enough to run the various sensitivity studies needed to assess the importance of input forcing data.As highlighted throughout the paper, our relatively simple snow model is expected to undergo improvements to its model physics in efforts to increase its potential accuracy and reliability, together with further analyses of the input forcing data, especially snowfall (extending the precipitation comparison of Boisvert et al., 2018).Examples of expected future improvements to NESOSIM include the following, in order of priority: -Incorporation of snow thermodynamics.The temperature evolution of the snowpack and snowmelt-refreeze processes is modeled, allowing us to run NESOSIM year round.Challenges will include accurately modeling or parameterizing the temperature profile through the snow layers and the possible retention of meltwater within the snowpack and its impact on the snow density.
-Increased vertical snow layers.Depth hoar is included as an explicit snow layer as we introduce the snow thermodynamics described above.Model validation will be an obvious challenge.
-Snow-ice formation.NESOSIM is currently run independent of the sea ice state, meaning we include no information regarding the potential for snow-ice formation -the depression of the snow layer below sea level and the conversion of snow to ice.This is thought to be particularly important for running NESOSIM in the Southern Ocean (e.g., Massom et al., 2001), where snow-ice conversion is expected to be more prevalent, but also in our North Atlantic sector (e.g., Granskog et al., 2017).Challenges will involve incorporating observed or simulated sea ice freeboard.
-Increased horizontal resolution.As we look towards the launch of NASA's ICESat-2 and the production of sea ice thickness from the derived freeboard product, we hope to increase the model resolution and conduct assessments of the ability of NESOSIM to capture smaller-scale (< 100 km) snow depth variability.
Snow depth and density information collected during the Norwegian Young Sea Ice (N-ICE2015) expedition (Merkouriadi et al., 2017) and the upcoming Multidisciplinary drifting Observatory for the Study of Arctic Climate (MO-SAiC) will provide crucial insight into the importance of smaller-scale phenomena not currently included in NE-SOSIM, while our model results can hopefully provide useful basin-scale context to the measurements being taken.NESOSIM is being made available as an open-source project (https://github.com/akpetty/NESOSIM) to encourage continued model development and active engagement with the snow on sea ice community.The model code is written in Python, an open-source programming language (Python Software Foundation, https://www.python.org/), to better enable future community development efforts.Our hope is that the model will continue to evolve as additional snow processes are incorporated, especially as new field and remote sensing snow observations are collected and made available for calibration/validation.

Figure 1 .
Figure 1.Schematic of the NASA Eulerian Snow On Sea Ice Model (NESOSIM) v1.0 presented in this study.The red (blue) text indicates processes that result in a loss (gain) of snow."Dynamics" indicates the combination of ice-snow advection and convergence-divergence which can cause either loss or gain of snow.
use the daily JRA-55 snowfall data from 15 August 1980 to 1 May 1991 and 15 August 2000 to 1 May 2015.The data were obtained from the National Center for Atmospheric Research's Research Data Archive at a horizontal resolution of 0.56 • × 0.56 • (∼ 60 km), downscaled from the original 1.25 • × 1.25 • Gaussian grid.The data are being produced on a near-real-time basis (2-to 6-month data latency).

Figure 2 .
Figure 2. (a) Warren climatology of August snow depth, (b, c) the initial conditions used in this study (broadly representing the snow depth as of 15 August) for 2012 and 2013, respectively, calculated using near-surface air temperature scaling.

Figure 3 .
Figure 3.Comparison of NESOSIM snow depth (a, c) and snow density (b, d) data with drifting Soviet station data collected between 1981 and 1991.Panels (a, b) show the mean seasonal evolution of the snow depth and density for the model (blue) and Soviet station data (black), with the data binned into the different months the data were collected.The shaded area represents 1 standard deviation from the annual monthly mean.Panels (c, d) show scatter plots of all points for which there were temporal crossovers.The r values indicate the correlation coefficient, while the colors indicate the different stations that collected the data.The NESOSIM data are from the default/ERA-I model configuration.

Figure 5 .
Figure 5. Map of the Arctic model domain and regions used in this study: AO: Arctic Ocean, CA: central Arctic, EA: eastern Arctic, NA: North Atlantic; BS: Bering Sea, LS: Labrador Sea (peripheral seas also discussed in the paper).

Figure 7 .
Figure7.Seasonal snow budget evolution across the four study regions (shown in Fig.5), initiated from 15 August 2000 to 2014 and run until 1 May of the following year using the default/MEDIAN-SF NESOSIM simulations.The thick lines show the mean daily regional values over this time period, while the shaded areas represent the interannual variability (1 standard deviation).

Figure 8 .
Figure 8. Snow budget terms as of 1 May, averaged over the 2001 to 2015 time period using the default/MEDIAN-SF NESOSIM simulations.The black lines show the four study regions used throughout this study.Panels (a) to (g) are integrals from the 15 August start date.Note the different color bar scales in panels (h) to (k).

Figure 9 .
Figure 9. Seasonal snow depth evolution across the four study regions (shown in Fig. 5) initiated from 15 August 2000 to 2014 and run until 1 May of the following year, forced by five different reanalysis snowfall products.This figure also includes results using the ASRv1 forced simulations (which are limited to 15 August 2000 to 1 May 2012).The thick lines show the mean (daily) regional snow depths over this time period, while the shaded areas represent interannual variability (1 standard deviation).All model runs use the default forcings/parameter settings.

Figure 10 .
Figure 10.Simulated snow depths on 1 May (averaged over 1 May 2001 to 2015), using the MEDIAN-SF snowfall forcing (a) and then the difference to the simulations forced by the four different snowfall products (b-e).The ASRv1 forced results are limited to 1 May 2012.Red (blue) colors indicate the individual reanalysis-forced simulations have higher (lower) snow depth.

Figure 11 .
Figure 11.Seasonal snow depth evolution across the four study regions (shown in Fig. 5), initiated on 15 August 2010, 2012, 2013 and 2014, and run until 1 May of the following year, forced by four different ice motion datasets and assuming no ice motion (NODRIFT).The thick lines show the mean (daily) regional snow depths over this time period, while the shaded areas represent interannual variability (1 standard deviation).All model runs use the default/MEDIAN-SF parameter settings.

Figure 12 .Figure 13 .
Figure 12.Modeled snow depth on 1 May (averaged over 1 May 2011, 2013, 2014 and 2015), assuming no ice motion (NODRIFT, a) and then the difference to the simulations forced by the four different ice motion products and the mean snow depth from the four different forced model runs.

Figure 14 .
Figure 14.Comparisons of the annual mean snow depths from NE-SOSIM (default configuration) forced with different reanalyses and the various Operation IceBridge (OIB) snow depth products.The blue (red) shading represents the annual mean spread across the different NESOSIM results (OIB products).The markers are spread across the shaded areas to improve readability.

Figure 15 .
Figure 15.Scatter plots of the three OIB snow depth products binned onto our 100 km model grid and coincident NESOSIM/MEDIAN-SF snow depth estimates for all years of data from 2009 to 2015, including the correlation coefficient (r) and RMSE.The contours show the kernel density estimate of the distributions.

Figure 16 .
Figure 16.As in Fig. 15 but showing comparisons of modified Warren snow depths (mW99) against the OIB snow depths.

Table 1 .
Default model forcings and parameter settings used by NE-SOSIM.

Table 2 .
Summary of the four different reanalysis datasets used in this study (data availability often subject to change/updates; information is given at the date of submission).
* Different resolutions available in some cases.NRT: near-real time.

Table 3 .
Summary of the different ice motion datasets used in this study based on information obtained at the time of submission.
* , SSM/I Jan 2007-present Public/NRT * These agencies produce drift datasets using different individual and/or combinations of satellite sensors not utilized in this study.NRT: near-real time.

Table 4 .
Mean snow depths as of 1 May across the four study regions (rows, regions given in Fig.5) for NESOSIM using different forcings and time periods (columns).The numbers in brackets represent interannual variability and are calculated as 1 standard deviation of the annual values.The default NESOSIM configuration is MEDIAN-SF snowfall, ERA-I winds, NSIDCv3 ice motion and Bootstrap ice concentration.Note that these 2011-2015 ice motion sensitivity runs exclude the 2012-2013 winter season due to the lack of Kimura drift data. *