A Lagrangian Snow Evolution System for Sea Ice Applications (SnowModel ‐ LG): Part II — Analyses

Sea ice thickness is a critical variable, both as a climate indicator and for forecasting sea ice conditions on seasonal and longer time scales. The lack of snow depth and density information is a major source of uncertainty in current thickness retrievals from laser and radar altimetry. In response to this data gap, a new Lagrangian snow evolution model (SnowModel ‐ LG) was developed to simulate snow depth, density, and grain size on a pan ‐ Arctic scale, daily from August 1980 through July 2018. In this study, we evaluate the results from this effort against various data sets, including those from Operation IceBridge, ice mass balance buoys, snow buoys, MagnaProbes, and rulers. We further compare modeled snow depths forced by two reanalysis products (Modern Era Retrospective ‐ Analysis for Research and Applications, Version 2 and European Centre for Medium ‐ Range Weather Forecasts Reanalysis, 5th Generation) with those from two historical climatologies, as well as estimates over ﬁ rst ‐ year and multiyear ice from satellite passive microwave observations. Our results highlight the ability of our SnowModel ‐ LG implementation to capture observed spatial and seasonal variability in

S e e h t t p://o r c a .cf. a c. u k/ p olici e s. h t ml fo r u s a g e p olici e s.Co py ri g h t a n d m o r al ri g h t s fo r p u blic a tio n s m a d e a v ail a bl e in ORCA a r e r e t ai n e d by t h e c o py ri g h t h ol d e r s .

Introduction
Sea ice plays a fundamental role in the Earth's climate system, influencing the amount of heat and moisture exchanged between the ocean and the atmosphere, and plays an important role in deepwater formation, a major driver of the global thermohaline circulation.Sea ice additionally protects coastal regions from wind waves and storm surges and influences coastal and marine ecosystems, as well as a number of economic and societal activities.Once the melt season ends, new snow will start to accumulate on the ice, increasing its insulative properties and impacting thermodynamic ice growth.If the snow is thick enough, it can depress the sea ice below the sea level surface, forming snow-ice (Jeffries et al., 1997;Rösel et al., 2018).The snow cover also limits the amount of photosynthetic light reaching the bottom of the ice, impacting underice biota (Mundy et al., 2007).When the melt season starts, the snowmelts provides water for melt pond development (Eicken et al., 2004;Iacozza & Barber, 2001;Perovich et al., 2002).Melt ponds in turn influence the sea ice energy balance and play an important role in how much ice melts each summer (e.g., Schroeder et al., 2014;Stroeve et al., 2014).When the melt water runs into the ocean, it impacts stratification of the upper ocean by providing a fresher surface layer (e.g., Toole et al., 2010).Snow on sea ice additionally increases the skin friction on a ship's hull making ice breaking more difficult (Canadian Government-Coast Guard, 2013).
Besides playing a key role in all aspects of the ice-ocean system, snow is a critical variable influencing the accuracy of thickness retrievals from satellite altimetry.Altimeters do not measure ice thickness directly but instead measure ice or snow freeboard, the height of the ice or snow surface above sea level.Snow depth and density are required to convert this freeboard to ice thickness (e.g., Kurtz et al., 2009;Laxon et al., 2003Laxon et al., , 2013)).Giles et al. (2007) estimated that snow depth and snow density uncertainties contribute 48% and 14%, respectively, to the total error in sea ice thickness retrievals using radar altimetry.Zygmuntowska et al. (2014) evaluated the corresponding impact of snow depth and density on laser altimetry, finding that snow depth is the main driver in thickness uncertainty, contributing up to 70% of the total uncertainty, with density playing a lesser role (30-35%).The difference in uncertainties is because radar measures the freeboard of the ice-snow interface, whereas lidar measures the total freeboard of air-ice or air-snow interface.
To date, estimates of snow depth and density on sea ice have primarily been based on a climatology derived from observations collected decades ago over multiyear ice (MYI) (Warren et al., 1999; hereafter referred to as W99) or on snow accumulation driven by output from atmospheric reanalyses in either a Eulerian or Lagrangian framework (e.g., Blanchard-Wrigglesworth et al., 2018;Kwok et al., 2009;Petty et al., 2018).
In the W99 climatology, monthly mean snow depth north of 80°N increases from a minimum of 4.2 cm in August to about 35.1 cm in May, with a maximum snow depth of 46 cm within the Lincoln Sea.
While observations collected during April 2017 in the Lincoln Sea were in agreement with the W99 climatology (Haas et al., 2017), it is unlikely that this climatology holds in every year or in every region.
Interannual variability in precipitation can be large (Barrett et al., 2020;Boisvert et al., 2018), and a climatology based on snow observations over MYI may no longer be valid in today's Arctic, which is dominated by first-year ice (FYI) (Maslanik et al., 2007(Maslanik et al., , 2011;;Stroeve et al., 2012;Stroeve & Notz, 2018).Furthermore, while sea ice extent trends have been largest in summer, ice conditions in recent years have become more anomalous during spring and autumn (Stroeve & Notz, 2018).This represents a distinct shift in the seasonality of the Arctic as the melt season begins earlier and lasts longer each year (Stroeve et al., 2014), reducing the time over which snow can accumulate.On the other hand, expanses of open water have led to a warmer and moister autumn and winter atmosphere (Boisvert & Stroeve, 2015;Serreze et al., 2012), which may increase winter precipitation amounts compared to 30 years ago.Since the use of an unrepresentative snow climatology can result in substantial biases in satellite-derived sea ice thickness if the snow depth and density depart strongly from this climatology (Bunzel et al., 2018), it is imperative to provide an operational product that captures spatial and temporal variability in snow depth over multiple altimeter missions.
Modeling approaches provide an attractive alternative to using a climatology, yet they are limited by biases and structural errors in reanalysis systems (Chevallier et al., 2017).Accurate representation of physics and interactions between the atmosphere, snow, and ice may also be important shortcomings in modeling approaches.Past approaches used either one atmospheric reanalysis data set and performed a simple accumulation of snow based on a temperature threshold in a Lagrangian framework (Blanchard-Wrigglesworth et al., 2018;Kwok & Cunningham, 2008) or implemented a simple two-layer snow budget model in a Eulerian framework with multiple reanalysis (Petty et al., 2018).These approaches generally neglected biases in reanalysis precipitation, which may differ by more than 30% over the Arctic Basin between different reanalysis products (Barrett et al., 2020), did not account for the full physical processes impacting snow depth accumulation or varying snow density, which in turn influence the snow water equivalent (SWE) (e.g., Sturm et al., 2010), and ignored summer melt processes.Other approaches to map snow depth (not density) have either relied on satellite passive microwave brightness temperatures to map snow depth over FYI (e.g., Markus et al., 2011) and MYI (Rostosky et al., 2018) or used a combination of Ka and Ku band satellite altimeters (e.g., Guerreiro et al., 2016;Lawrence et al., 2018).None of these currently available products provide snow depth and density through a full annual cycle (see Table S1 in the supporting information).
To fill these gaps, SnowModel, originally developed by Liston et al. (2018) to run over sea ice surfaces, was modified to run on a pan-Arctic scale using satellite-derived ice-motion vectors from passive microwave observations, buoy movement, and surface wind fields (Tschudi et al., 2019).This new Lagrangian snow evolution model (SnowModel-LG) accumulates snow in a Lagrangian framework (Liston et al., 2020), producing pan-Arctic daily snow depth, density, and grain size distributions, now available from August 1980 to 31 July 2018 at 25-km spatial resolution.Key improvements in SnowModel-LG compared to SnowModel include blowing snow sublimation, creation of superimposed ice and incorporating ice dynamics.This effort is described in Part I of this series of papers (Liston et al., 2020).This study forms Part II, an evaluation of SnowModel-LG results against snow depth estimates from Operation IceBridge (OIB) snow radar data, ice mass balance buoys (IMBs), snow buoys, MagnaProbes, and rulers.Further intercomparisons are made against retrievals from passive microwave (Markus et al., 2011;Rostosky et al., 2018).We additionally evaluate SnowModel-LG results (depth and density) against two climatologies (Shalina & Sandven, 2018;Warren et al., 1999).Finally, we evaluate regional seasonality in snow depth and density, their interannual variability, and temporal trends since 1980.

Snow Mass Balance Model
This study relies on the newly developed Lagrangian snow evolution model (SnowModel-LG) described in the Part I companion paper (Liston et al., 2020).The model accumulates snow on sea ice using atmospheric inputs of air temperature, humidity, precipitation, wind speed, and direction.SnowModel-LG runs in a Lagrangian framework, where satellite-derived ice motion vectors redistribute individual ice parcels around the Arctic Ocean.The ice parcels and their regridding to a 25-km resolution product provide higher resolution than previous studies (Blanchard-Wrigglesworth et al., 2018;Kwok et al., 2009;Petty et al., 2018).
Compared to earlier methods, SnowModel-LG also has the advantage of including the physics required to simulate most processes that drive the seasonal evolution of snow depth and snow density over sea ice, such as blowing snow redistribution and sublimation, density evolution, snow pack metamorphosis, and formation of superimposed ice.SnowModel-LG is based on components developed previously for terrestrial snow applications with some new features included for sea ice.
A new model component, SnowDunes, was developed to simulate snow redistribution over level ice (Liston et al., 2018).These components make SnowModel-LG unique compared to previous methods since wind and ice roughness play important roles in the snow distribution over sea ice.To accumulate snow on 25-km grid cells and the results that follow, SnowTran-3D was turned off as the scale of snow bedforms (1-10 m) and local-scale topographic variability (1-500 m) are negligible at such a coarse spatial resolution.
The model can be forced with atmospheric information from either atmospheric reanalysis or automatic weather stations using MicroMet (Liston & Elder, 2006).For more details on SnowModel-LG and its implementation for this paper, the reader is directed to the Part I companion paper Liston et al. (2020).
Both reanalyses extend from 1980 onward.Barrett et al. (2020) provides an assessment of the performance of several reanalysis precipitation estimates, including MERRA-2 and ERA5, against gauge-corrected precipitation observations from Russian North Pole drifting stations.While all reanalysis capture the basic spatial pattern and interannual variability of precipitation across the Arctic Ocean, they differ with respect to the magnitude of the precipitation and the number of wet days each month.MERRA-2 is overall wetter than ERA5.The use of two reanalysis helps to better understand the impact of reanalysis product on overall snowpack accumulation and density evolution.

Snow Depth Data Sets
To evaluate SnowModel-LG output against other snow depth estimates, we make use of several snow depth products, each of which provides a different representation of snow depth at different locations, at different times of year, and with different levels of uncertainty.Given disparate levels of uncertainties between various measurement techniques and the large spatial variability of point measurements compared to a 25-km grid cell, we focus on how the different sources of snow depth information represent mean snow depth and standard deviation as well as snow depth distribution (i.e., probability distribution functions (PDFs)), rather than the root-mean-square error.Several algorithms have been applied to this data to retrieve snow depth, which differ in magnitude based on how the air-snow and snow-ice interfaces are detected in the radar returns, as well as how noise and resolution are addressed (see Kwok et al., 2017).Considering the uncertainty in the OIB-derived snow depths themselves, we make use of four different OIB snow depth products, all of which have been previously evaluated in Kwok et al. (2017): (1) the standard product, originally described in Kurtz et al. (2013Kurtz et al. ( , 2015) ) and available from the National Snow and Ice Data Center (NSIDC), hereon referred to as the Goddard Space Flight Center (GSFC) product; (2) the Quicklook product, which acts as an interim to the standard product and uses waveform fitting described in Kurtz et al. (2014), also available from NSIDC and here referred to as Quicklook (GSFC-NK in Kwok et al., 2017); (3) the revised version of the Kwok and Maksym (2014) algorithm, referred to as the Jet Propulsion Laboratory (JPL) product; and (4) the Snow Radar Layer Detection (SRLD) algorithm, described originally in Koenig et al. (2016) to detect the air-snow and snow-ice interfaces over ice sheets but adapted for sea ice as described in Kwok et al. (2017).Details in processing of each snow depth product are found in the references indicated.Data for each product were obtained from 2009 to 2015, QuickLook through 2016.

In Situ Data
In situ snow depth comparisons come from several sources, including IMBs, MagnaProbes, and rulers.Figure 1 summarizes the locations of all in situ and buoy data used, and individual listing of dates for each buoy is provided in Table S2.Here we use IMBs from the U.S. Army Cold Regions Research and Engineering Laboratory (CRREL) and snow buoys from the Alfred Wegener Institute (AWI) (Nicolaus et al., 2017).CRREL originally deployed IMBs in MYI near the North Pole in 1993, with some drifting through Fram Strait and others circulating in the central Arctic, including the Beaufort Sea (Perovich et al., 2017).Snow buoys used in this paper span 2013 to 2016.In total we used 11 AWI and 26 CRREL buoys.Measurements were daily averaged.
MagnaProbe (Sturm & Holmgren, 2018) and/or ruler measurements come from those made during the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment in 1997 to 1998 and from the Norwegian young sea ICE (N-ICE2015) expedition that took place during winter and spring of 2015.SHEBA snow depths were collected using MagnaProbes and/or a graduated ski pole along survey lines throughout the entire year.A total of 21,169 locations in April and May was sampled, giving a mean snow depth of 33.7 ± 19.3 cm (Sturm et al., 2002); thinner snow was observed over flat terrain, and snow depths increased for heavily rubbled ice.Here we have examined all snow lines, which ranged from 200 to 500 m long, and from flat terrain to ridged ice.N-ICE observations came from the R/V Lance, which was originally frozen into the ice near Svalbard at 83°N/22°E and drifted toward Fram Strait.Since the ship drifted toward the ice edge too fast, the ship was relocated three times between 15 January and the end of June.Both FYI and second-year ice types were encountered.Snow thickness surveys were conducted along 800 m long near-weekly transects, with some 2-4 km nonrepeated transects from January to June (Merkouriadi et al., 2017;Rösel et al., 2018).

Satellite Passive Microwave Data
Passive microwave-derived snow depths offer the possibility for a pan-Arctic perspective, and over the full 1979 to present time period.Markus et al. (2011) developed the first method to retrieve snow depth over FYI during the cold season, using brightness temperatures (Tbs) from the Scanning Multichannel Microwave Radiometer, Special Sensor Microwave Imager, and the Special Sensor Microwave Imager Sounder, hereon referred to as PM fyi .Snow depth was only computed over FYI because Tbs over MYI and deep snow were found to be ambiguous.Further, retrievals were limited to premelt conditions as the passive microwave signal is dominated by surface scattering when the snow surface is wet, increasing the emissivity and therefore the Tbs, resulting in unrealistically large snow depths.While the available data product attempts to mask out these wet snow areas to avoid erroneous snow depths, some artifacts remain.Further, the upper snow depth limit is set at 100 cm, and snow depths are only given at 5-cm increments.However, while values up to 100 cm are present, Markus et al. (2006) state the upper limit is about 50 cm.Mean biases with respect to OIB data over FYI were found to be close to 0 cm, but locally the snow depth can be significantly biased (Brucker & Markus, 2013).Data are provided at 25-km spatial resolution in the polar stereographic grid.
The second passive microwave-based snow depth data product comes from a new algorithm developed by Rostosky et al. (2018), which extends the Markus et al. (2011) methodology to cover MYI regions using lower frequency (6.9 GHz) Tbs from the NASA Advanced Microwave Scanning Radiometer AMSR-E and the JAXA Global Change Observation Mission-Water AMSR2 instruments.The approach uses OSI-SAF ice type information to separate out FYI from MYI, and then during March and April, snow depth over MYI is determined by calibration against OIB-estimated snow depths.Thus, while the product extends over MYI, it only does so during the springtime when OIB data were available to constrain the retrievals.It is unclear how representative the results are outside of the OIB time period or regions sampled, yet the regression against OIB was found to be stable in time.This product is referred to as the PM all and is available during the AMSR-E/AMSR2 time period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) at 25-km spatial resolution in the polar stereographic grid.

Climatologies
Finally, we compare our pan-Arctic results to those from Warren et al. (1999) as well as the updated climatology from Shalina and Sandven (2018).SS18 improve upon the W99 climatology by including snow depth data from Soviet airborne (Sever) expeditions collected over 28 years (1959 to 1988) and combining those data with the North Pole drifting stations used in the W99 climatology.The Sever data set contains snow depth and density measurements on the runway and surrounding areas during aircraft landings.Additional statistics regarding snow distribution over hummocks, ridges, and sastrugi are provided.Since aircraft land on level (and often thinner) sea ice, measurements carried out near the aircraft preferentially sampled level ice, which may not be representative of the broader snow depth distribution over more deformed older ice.Further, considering a pan-Arctic perspective, the landings in the 1970s are most representative of the Arctic Ocean, whereas the measurements in the 1960s and 1980s were limited to the coastal regions.This results in an updated snow climatology for March-April-May that covers both MYI and FYI.Details on the merging of the North Pole drifting station data and the Sever data can be found in Shalina and Sandven (2018).

Processing Steps
Because the SnowModel-LG grid cells are larger than any of the in situ and OIB data, there are multiple snow depth values within each SnowModel-LG grid cell.The latitude and longitude of each in situ or OIB value are first matched to the closest SnowModel-LG grid cell using center coordinates.All locations falling within the grid cell are averaged together, and standard deviation and number of values are recorded.
For comparisons with buoy data, as well as the SHEBA and N-ICE data, we start with the initial in situ location (i.e., buoy initial latitude and longitude) and follow the parcel through time.In other words, the SnowModel-LG time series is from a representative parcel located at the position of the earliest date of in situ snow depth observation and then that parcel is followed as it moves around the Arctic Ocean.For comparisons with larger scale (i.e., satellite and climatologies), we first grid all the SnowModel-LG parcel data to the 25-km Equal Area Scalable Earth (EASE) grid and compute monthly means.Comparisons with PM fyi and  PM all are made for October through April: SS18 as an average for March, April, and May since individual snow depths for each month did not vary substantially (not shown), and for W99, we compare all months of the year.

Results-Snow Depth Assessments
Below, we show results using two simulations from SnowModel-LG, one forced using ERA5 and the other with MERRA-2 atmospheric reanalysis.As described in Part I, annual mean precipitation from each of the reanalysis was adjusted using mean calibration scaling factors averaged over 8 years (2009 through 2016) determined by comparing SnowModel-LG snow depths with OIB snow depths from the QuickLook product.In short, for each year, OIB snow depths collected the during spring campaign are assumed to represent snow depth on 1 April.These aggregated OIB snow depths are compared with SnowModel-LG snow depths for 1 April, and the reanalysis precipitation is adjusted by a scaling factor so that the difference between OIB and SnowModel-LG is minimized.Mean precipitation scaling factors for the 2009 through 2016 time period are applied across all years .Mean scaling factors are 1.37 for MERRA-2 and 1.58 for ERA5.More details are provided in Part I.Note that since this is a mean calibration adjustment designed to set the mean difference averaged along all OIB tracks to 0 cm, spatial pattern differences between SnowModel-LG and the OIB data at individual grid points are preserved.It also places the two reanalysis simulations on common levels, so we can look at spatial pattern differences without being distracted by differences in the means.
To compare spatial patterns of snow depths using ERA5 and MERRA-2, we show monthly mean snow depths from two accumulation seasons spaced three decades apart: August 1985 to July 1986 and August 2015 to July 2016.While these individual years are arbitrarily chosen, they represent years in the early and latter part of the 39-year record and highlight the importance of interannual variability and long-term trends in snow accumulation.For reference, mean monthly snow depth for each SnowModel-LG run during these years are summarized in Figures 2a and 2b.Spatial patterns of snow depth follow the expected pattern of accumulation within the central Arctic Ocean already in September and within the Atlantic sector in winter as a result of the North Atlantic storm tracks, with redistribution of snow within the Beaufort Gyre and the Transpolar Drift Stream.Seasonally snow begins accumulating already in August, with maximum snow depths reached in either March or April.While spatial patterns between ERA5 and MERRA-2 are broadly similar, a key difference is that MERRA-2 keeps its snow longer in June (see also Part I; Liston et al., 2020).Snow depth is seen to differ substantially between the 2 years, both in spatial pattern distribution and magnitude.Interannual changes will be discussed further in section 5.

Comparisons With Buoys and Field Campaigns
As mentioned in section 2.3, to compare SnowModel-LG results with buoys and other field campaign data, we extract daily modeled snow depths from grid cells that correspond with locations and date of in situ data and follow that parcel through time (i.e., we start with the initial position of a buoy and find the same parcel at that position and then follow that parcel along the buoy drift track).The data from SHEBA provide a nearly yearlong evolution of snow depth from the 1990s.In our SHEBA comparisons (Figure 3), we aggregate all in situ observations for each day over all transects and show the corresponding SnowModel-LG results from a representative parcel located during the start of the SHEBA (Uttal et al., 2002) snow depth observations.We used the fifth nearest parcel to that start date/location because it was the only parcel old enough to have accumulated enough snow to compare reasonably well with the SHEBA observations; parcels 1-4 were too young to have much snow accumulated even though they were subjected to the same general synoptic meteorological forcing.This highlights the fact that in addition to biases in atmospheric forcing and SnowModel-LG snow-physics representations, model results are impacted by errors in parcel concentrations and trajectories (Liston et al., 2020).In general, SnowModel-LG tracks the SHEBA snow depths well, showing increases from August to May and declines thereafter.During July and August, SnowModel-LG loses all its snow, earlier with ERA5 than with MERRA-2, whereas the MagnaProbe data show between 2 and 8 cm of snow remains.Faster loss of snow in ERA5 may be in part a result of a warm bias in ERA5 (e.g., Graham et al., 2019) and also because ERA5 accumulates less snow in winter.Further, Part I (Liston et al., 2020) shows that the summer melt period under ERA5 forcing is 3 weeks longer than from MERRA-2, which leads to earlier removal of snow.During other months, results using MERRA-2 and ERA5 fall within the range of in situ snow depths, while the SHEBA observations snow more daily variability.
A few tens of thousands of snow depth measurements are available during the N-ICE field campaign.For each leg of the cruise, we plot the starting position of each N-ICE snow depth observation over the entire year from 1 August 2014 to 31 July 2015 together with the observed mean snow depth (red X) and ±1σ (gray shading) (Figure 4).N-ICE snow depths show significant small-scale spatial variability, reflected by the large standard deviation.SnowModel-LG results span the range of variability observed during N-ICE.The divergence between the MERRA-2 and ERA5 in mid-November is the result of a storm that was represented differently in the two atmospheric forcing data sets.Interestingly, the location of the N-ICE ship corresponds to an area where snow depths from ERA5 tend to exceed those from MERRA-2, except during the second leg.
Finally, in contrast to the SHEBA and N-ICE comparisons, where hundreds of snow depth observations are made over transects of varying lengths each week, the IMBs and snow buoys represent point measurements at one location as the ice moves around the Arctic.Given the variable nature of snow over sea ice, especially over deformed sea ice, the IMB and snow buoy snow depths will not be representative of the coarse 25-km SnowModel-LG grid cell.Even over nondeformed areas, or over smooth landfast ice, there can be several meters' length scale related to the formation of snow dunes (Moon et al., 2019).Considering that the dominant errors between buoy data and model outputs are the initial conditions, rather than the snow depth evolution, we set the first buoy observation equal to the average of the MERRA-2 and ERA5 snow depth values on that date and time and let the buoy data evolve with that initial offset.This effectively removes any initial condition error so that we can analyze the snow depth evolution over time.On eight of the buoy data sets, unrealistic jumps in the data (i.e., snow depth increasing by more than 30 cm in 1 day) were removed from the rest of the time series.Unrealistic jumps happened six times (or 0.16% of the time out of 3,712 time steps).Results are presented in Figure 5, with results separated for June to July (red) and August to May (black), representing summer and winter, respectively.Individual buoy plots for each of the 37 buoys have been uploaded as supplemental data.Overall, winter snow depths correlate well, with r 2 values of 0.72 and 0.71 for MERRA-2 and ERA5, respectively.Once the snow cover is melting, however, large disagreements are found, especially with the ERA5 results, which quickly melts away the snow cover, while the buoys maintain a deeper snowpack.This was also observed in the SHEBA comparisons and suggests that SnowModel-LG forced with ERA5 may melt the snow too quickly and that further study of SnowModel-LG's representation of summer processes, or the appropriateness of the ERA5 summer forcing, is warranted.

Snow Depth During OIB
Since snow depths from other OIB products differ from the QuickLook product used to calibrate the reanalysis precipitation, it is helpful to compare results against other OIB data products.Also the calibration only used one mean value for the spatial domain across all years, so the interannual and spatial variability of the model is independent from the product; those are the features we are analyzing here (not the mean in space and time that we have prescribed).Table 1 summarizes means and plus/minus one standard deviation (1σ) of the different OIB snow depth products averaged along all OIB tracks for each year, together with corresponding values for the SnowModel-LG runs forced by MERRA-2 and ERA5.Corresponding scatterplots for each year are provided in the supporting information (Figures S1 and S2), while histograms are summarized for each year in Figure 6.Agreement among SnowModel-LG and the various representations of snow depth from OIB varies from year to year, with no clear OIB product better matching the modeled snow depths than another.A key concern in this analysis is the large difference in both magnitude and interannual variability of mean snow depths among the different OIB-derived values.For example, in 2009 the SRLD mean snow depth is 13 cm compared to 30 cm from the JPL product.In this year, both the MERRA-2 and ERA5 simulations suggest even deeper snow packs at 33 and 31 cm, respectively.Conversely, in 2014, the SRLD product results in a mean snow depth of 32.9 cm compared to only 20.8 cm from QuickLook, with the MERRA-2 and ERA5 simulations falling somewhere in between at 25.6 and 27.5 cm, respectively.
SnowModel-LG mean snow depths from MERRA-2 and ERA5 are generally within 3 cm of each other, with broadly similar PDFs.In some years, the PDF between the modeled and OIB snow depths matches best with the GSFC product (i.e., 2009, 2010, and 2015), while in others it is better with the QuickLook (2012) or the JPL values (2011).For all years combined, SnowModel-LG PDFs most closely follow those from the standard GSFC product.
Considering that there are differences among the OIB snow depth products, we additionally compare the spatial pattern of mean snow depths from averaging the four OIB snow depth products together to those from SnowModel-LG using MERRA-2 (Figure 7); results for ERA5 are similar (Figure S3).Individual spatial plots for each OIB product are additionally provided in the supporting information (e.g., Figures S4-S7).Spatially, SnowModel-LG tends to have shallower snow depths north of Greenland and the Canadian Archipelago and thicker snow packs closer to the pole, though the magnitudes of the difference vary from year to year.While differences are generally within 10 cm of each other, at times they can reach 30 cm.The spatial pattern of differences is at times consistent with the individual OIB products, especially in 2013 when SnowModel-LG systematically simulates nearly 20 to 30 cm less snow than any of the OIB products directly north of Greenland.However, other years such as 2009, 2012, and 2014 show large regional snow depth differences among the OIB data products, resulting in different spatial pattern differences with SnowModel-LG depending on which OIB data product is used.Time series from 2000 to 2018 further reveal large systematic differences in FYI snow depths at the start of the freeze-up period (Figure 9).In some years the PM fyi matches closer to SnowModel-LG than in others, though consistently higher PM fyi snow depths are found early in the freeze-up period.Later in winter and into spring, the PM fyi snow depths converge with those from SnowModel-LG (Figure 9).Differences are noticeably less over MYI; PM myi in March and April is mostly less than that from SnowModel-LG, though in some years they match well (i.e., 2008, 2013, and 2015).

Comparisons With Passive Microwave
From a climatological perspective, we find that the 1981 to 2010 climatological mean FYI snow depth averaged for each month from October to April from PM fyi ranges from 14.5 cm (October) to 17.0 cm (April), whereas snow depths from SnowModel-LG forced with ERA5 and MERRA-2 have larger seasonal variability, ranging from 6.3 to 18.6 cm and 6.9 to 21.0 cm, in October and April, respectively (Table S3).This highlights the result that the passive microwave data likely have a positive snow depth bias at the start of the freeze-up, though by springtime, results are more in line with snow depths expected over FYI.

Comparisons Against Climatology
Finally, SnowModel-LG results are compared against the W99 and SS18 climatologies.Ideally, we should recreate the climatology by extracting snow depth values from all points that entered into the climatology.Since this is not feasible given the larger uncertainties in atmospheric forcing data prior to the satellite    data record, we chose again to compare to two individual years three decades apart to highlight the change from climatology in more recent years.Further, since it is well known that the W99 climatology is not applicable outside of the Arctic Basin we exclude Baffin Bay/Davis Strait, the Barents and Kara Seas, the Bering Sea and the Sea of Okhotsk in these comparisons.While many CryoSat-2 sea ice thickness products rely on the W99 climatology, with snow depths halved over FYI regions, we highlight here differences with the climatologies as they are.For reference, corresponding results from halving the W99 climatological snow depths over FYI based on OSI-SAF first-year/MYI identification for August 2015 to July 2016 are provided in Figure S8 from the MERRA-2 simulations.
With respect to the W99, SnowModel-LG has considerably less snow in 2015/2016 almost everywhere the Arctic Basin (Figure 10), with the exception of a small region of up to 30 cm positive snow depth differences that grows in magnitude from November to May 2016 north of Svalbard.Higher snow depths in this region from SnowModel-LG may reflect increased precipitation from the North Atlantic storm "Frank" that occurred in December 2015/January 2016.Since this region is FYI, halving the Warren et al. (1999) climatology results in a larger area where SnowModel-LG snow depths exceed those from W99 (Figure S8).In contrast, in 1985/1986, SnowModel-LG has overall deeper snow packs that compare more favorably with W99, showing both positive and negative differences on the order of ±10 cm, though regions of larger differences do occur.In particular, the MERRA-2 simulations show a broad region of deeper snow than W99, which grows from localized small positive differences north of Svalbard in October and spreads across the Arctic Ocean in February through May, stretching from north of Greenland toward the coast of the East Siberian and Laptev seas.At times the difference from W99 exceeds 20 cm.Again, this likely reflects the impact of winter storms within the North Atlantic, which bring precipitation to the Greenland, Barents, and Kara seas (e.g., Koyama et al., 2017;Stroeve et al., 2011).In June, SnowModel-LG consistently has less snow than W99, with differences exceeding 30 cm in 2016 and also from the ERA5 simulation in 1986.The negative departure from W99 in June may reflect slower melting relative to climatology in these years.The use of a FYI region flag to halve W99 snow depths results in artificial boundaries in snow depth that are unlikely to be representative of true snow depth variability (Figure S8).
Differences between SnowModel-LG and the new Shalina and Sandven (2018) climatology, shown as a March, April, and May average, reveal similar spatial patterns of positive and negative differences as in W99, but the mean differences are less (typically within 20 cm) (Figure 11).This is mostly a result of larger overall snow depths in W99 compared to SS18 (Figure S9).Results between ERA5 and MERRA-2 are quite similar in 2016, showing thinner snow packs within the Beaufort, Chukchi, East Siberian, and northern Barents Sea and thicker snow packs north of Greenland.In 1986 however, MERRA-2 has more snow accumulated than ERA5 throughout the Arctic Ocean, leading to broad areas of deeper snow packs than the SS18 climatology, which can exceed 25 cm (e.g., East Siberian coastal region).
Finally, in terms of how well SnowModel-LG captures the climatological seasonal cycle, monthly mean snow depths for the central Arctic (CA) from W99 and from SnowModel-LG as averaged from 1981 to 2010 are summarized in Figure S10.Results for other regions are only shown for SnowModel-LG.While the magnitude of snow depth is less in SnowModel-LG, the seasonal cycle is well represented.The key difference is that the snow depth decreases more quickly in SnowModel-LG after the maximum snow depth is reached, whereas W99 snow depths remain relatively high between May and June.

Results-Snow Density
While snow depth is more routinely measured than snow density, snow density is essential for conversion of snow water equivalent (SWE, unit of snowfall) to snow depth, and it can range from as little as 10 kg m −3 to more than 350 kg m −3 , depending on the air temperature and wind conditions (Doesken & Judson, 2000).This means that 1.0 cm SWE of snowfall can in theory (with no vertical compaction) result in 100 to 3 cm of snow (see conversion in Part I).Snow density develops with crystal structure and relative humidity.As soon as the snow accumulates on the sea ice, snowpack densification begins, reducing the snow porosity (fractional volume of air or liquid portion of the snow pack) (e.g., Mizukami & Perica, 2008).At the same time, snow metamorphism results in snow grain modifications.How this impacts the snow density varies depending on the form of metamorphism.In general, the rounder the grain size, the higher the density, which is typically the case for warmer snowpacks.When the snowpack is melting, liquid water fills the 10.1029/2019JC015900 Journal of Geophysical Research: Oceans pore space, further increasing the snowpack density while decreasing snow depth.On the other hand, depth hoar undergoes slower densification (Mizukami & Perica, 2008).
The advantage of SnowModel-LG is that it simultaneously evolves SWE and snow density, from which the snow depth can be determined.Snow depths measured at Soviet drifting stations showed that snow density increases from 250 kg m −3 in September to 320 kg m −3 in May. Warren et al. (1999) found regional variations were small (25 kg m −3 in May).Merkouriadi et al. (2017) report an average snow bulk density of 345 kg m −3 during the N-ICE2015 expedition, though it was found to vary with ice type: over FYI the mean density was 300 ± 54 kg m −3 , while over second-year ice it was 381 ± 237 kg m −3 .During SHEBA, bulk snow density was 340 kg m −3 (Sturm et al., 2002), while the same area from W99 has a mean of 330 kg m −3 .Snow densities observed during the Multidisciplinary drifting Observatory for the Study of the Arctic Climate (MOSAiC) expedition in autumn/winter, snow density was found to range from 218 to 500 kg m −3 .
These in situ observations are useful for evaluating the representativeness of seasonal evolution of SnowModel-LG density results.Figure 12 once again focuses on the years 1985/1986 to 2015/2016 to highlight both the seasonal and interannual variability between decades.Snow density increases from low values in August (~200 to 290 kg m −3 ) to values around 300 to 350 kg m −3 over the central Arctic Ocean in November, up to densities as large as 450 kg m −3 by May over parts of the Arctic Ocean.These densities increase further in June, when the snow is melting, to values approaching 500 kg m −3 , widespread in MERRA-2 in June 1986 and less so in June 2016.Very dense snow (>500 kg m −3 ) remains in July 1986 through most of the central Arctic Ocean from the MERRA-2 simulations.In stark contrast, snow density from the MERRA-2 simulations in July 2016 drop to values near 200 kg m −3 over parts of the Arctic Ocean, because the most of the winter snow has melted and new snow is accumulating, with just a few locations remaining with dense snow from the previous winter.The contrast between ERA5 and MERRA-2 snow densities is most notable in the months of June and July, because of the large differences in snowmelt (see also Part I; Liston et al., 2020), whereas during the winter months they are more consistent.Overall, SnowModel-LG densities are consistent with low values for damp new snow along the ice margin (e.g., 100-200 kg m −3 ), higher values in the interior (~325 kg m −3 ), and increasing over time as snow settles and compacts such that maximum snow density is reached during snow melt in June.This is consistent with seasonal increases from September to May observed by W99, though our seasonal maximum occurs in June rather than May (see also Figure S10).Results are also consistent with values recently observed during the MOSAiC Arctic drift experiment from October to May.
While the seasonal variability in modeled densities are broadly consistent with climatology (e.g., Figure S10), spatial differences can be large in individual years/months.For example, SnowModel-LG density in 1985/1986 is larger than that from W99 in most months and regions, with the exception of the coastal regions, and north of Greenland in August from MERRA-2 and also in the Laptev Sea in August from ERA5 (Figure 13).The largest departures from climatology are found in June (MERRA-2) and also in July (MERRA-2 and ERA5).In June, the MERRA-2 simulations can have up to 180 kg m −3 denser snow, whereas in July, the modeled simulations show at least 150 kg m −3 less dense snow than W99, with an exception in July 1986 when the MERRA-2 results have similar magnitude differences but in the opposite direction.Interestingly, the regions where W99 snow density in 1985/1986 is larger than SnowModel-LG in winter expands throughout large parts of the Beaufort, Chukchi, East Siberian, and Laptev seas in 2015/2016, whereas the W99 snow density remains less than SnowModel-LG north of Greenland and the Canadian Archipelago.Again, this highlights that how well SnowModel-LG matches climatology depends strongly on year and month evaluated.Nevertheless, differences over the central Arctic Ocean are usually within 50 kg m −3 in winter, the region for which the W99 is most valid.Note further that there are some noticeable artifacts in the W99 snow density fields, reflecting the use of the 2-D quadratic function that was fitted to all available SWE observations (e.g., see Figure 2 in Tilling et al., 2018).This can result in artificial boundaries in snow densities, as observed particularly in the August and July difference plots.

Regional and Pan-Arctic Interannual Variability
Finally, we assess how snow depth and density are changing over time.Trends in snow depth and density for each month are shown in Figures 14a and 14b, while time series of regional mean values averaged for winter (October to March) and summer (April to September) are summarized in Figures S11a and S11b (see Meier et al., 2007, for the regional descriptions).Individual monthly snow depth trends are listed in Table 2, along with 1981-2010 mean averages from SnowModel-LG and corresponding regional W99 regional estimates.Climatological mean seasonal cycles are provided in Figure S10, while regional annual snow depths are further shown in Figure S12.Here we include both the Kara and Barents seas in the regional means in Table 2, though the W99 values should be treated with caution.
Regionally, the central Arctic has the deepest snow packs in winter and also in summer from MERRA-2; summer ERA5 simulations have individual years when the snow within the East Greenland Sea exceeds that over the central Arctic.The shallowest snow packs are found in the Barents Sea, followed by the Laptev, Chukchi and Kara seas.Noteworthy is that SnowModel-LG predominantly simulates an overall decline in snow depth since 1980 in all regions and months except for the East Greenland Sea and also within the central Arctic Ocean from ERA5.These trends are for the most part statistically significant at the 95% confidence interval, though none of the positive trends are significant except directly north of Greenland from ERA5 during winter.Journal of Geophysical Research: Oceans The largest negative trends in snow depth are generally found in April and May throughout the Barents, Kara, Laptev, East Siberian, Chukchi, and Beaufort seas, larger for MERRA-2 than for ERA5.One region where MERRA-2 and ERA5 diverge is the Beaufort Sea, where MERRA-2 trends are at least twice as large as those from ERA5, and trends from ERA5 are mostly not statistically significant as averaged over the entire Beaufort Sea domain (e.g., Table 2; individual data points do show significant declines).Negative trends in these regions are consistent with later autumn freeze-up in recent years (e.g., Stroeve & Notz, 2018), reducing the time over which snow can accumulate.They are also consistent with an overall transition in these regions to FYI as the MYI has been reduced by 90%.On the other hand, over the central Arctic, snow depth trends are smaller (or positive from ERA5) during winter, increasing in magnitude during summer (−2.0 to −3.0 cm per decade).The largest negative snow depth trends occur in June over the East Siberian Sea (−4.1 cm per decade from MERRA-2) and in the Kara Sea (−3.3 cm per decade from ERA5).While the Note.Corresponding values from the Warren et al. (1999) climatology are also given.Note that regional snow depths are only computed over ice-covered pixels.
Positive trends are in red.Bold indicates statistically significant trends at the 95% confidence level. 10.1029/2019JC015900 Journal of Geophysical Research: Oceans trends may not seem large, over the last 40 years these trends represent an 8-to 16-cm decline in total snow depth, respectively.These results are consistent with those previously reported in Webster et al. (2014).
Compared to W99, which was based on snow depth measurements over MYI from 1954MYI from to 1991MYI from , our corresponding 30-year (1981MYI from -2010) ) climatological snow depth values are less.This is true also in the central Arctic where large amounts of MYI remain, yet in this region, SnowModel-LG results from MERRA-2 are reasonably consistent with those from W99 (within 5 cm for most months).In the central Arctic, the seasonal cycle from SnowModel-LG is also broadly consistent with that from W99, showing ~30 cm of snow accumulation.In other regions, the accumulation of seasonal snow is on average 10 cm less than from W99, though can be as much as 30 cm less in the Barents Sea, a region where the W99 climatology is not applicable.The only regions where the 1981-2010 mean SnowModel-LG snow depths exceed those from W99 occur in the Kara Sea in May through January and in the Barents and Laptev seas in July (MERRA-2 only).However, caution is needed in these regions because W99 is only valid for the central Arctic.
Finally, in regard to snow density, there is quite a bit of interannual variability (Figure S11b), with slightly negative trends in the marginal seas around the central Arctic Ocean in winter, and no trend over the central Arctic Ocean and the Lincoln Sea nor in the East Greenland Sea (Figure 14b).Negative snow density trends are statistically significant mostly within areas of the Beaufort, Chukchi, East Siberian, Laptev, Kara, and Barents seas during freeze-up until around April/May, when density trends weaken.Trends are larger and more widespread for MERRA-2 than for ERA5.Negative trends are also statistically significant within the East Siberian and Laptev seas in June and over the central Arctic in July from MERRA-2.The July snow density trends in Figure 14b, however, have little meaning, because the snow in July is minimal (Figure 2b) and sensitive to periodic melt and new snowfall.Regionally, monthly mean snow density is largest within the central Arctic and East Greenland Sea areas and generally lowest in the Laptev and Chukchi seas (Figure S10).

Discussion
Our implementation of SnowModel-LG results in different spatial representations of snow compared to previous studies that have also relied on atmospheric reanalysis.Figure 15a  As expected, there is also increased spatial structure reflected in the snow density simulations compared to previous studies as well as differences in magnitude.SnowModel-LG densities are larger than those from NESOSIM within the northern Barents Sea, Fram Strait, and near the pole in November (ranging from 350 to 400 kg m −3 in SnowModel-LG compared to 270 to 300 kg m −3 in NESOSIM) but generally smaller over FYI regions outside of the Barents Sea.In April, SnowModel-LG snow density is larger, with values from MERRA-2 reaching as high as 400 kg m −3 near the pole and within Fram Strait, whereas NESOSIM densities range from 300 to 320 kg m −3 over most of the Arctic Ocean.However, in Baffin Bay, snow density from SnowModel-LG is around 250 kg m −3 , compared to ~300 kg m −3 from NESOSIM.NESOSIM snow density is a weighted average of the two snow layers, with snow density fixed at 200 kg m −3 for new snow and 330 kg m −3 for old snow.
We find the mean length scale of SnowModel-LG is similar to the mean length scale of coincident OIB measurements (Figures S13 and S14).Further, we find that the scale of the modeled snow depth anomalies reflects the scale of the governing processes.If the heterogeneity were a result of noise introduced by the grid resolution or divergent parcel creation, we would not see the same scaling of the spatial heterogeneity.
SnowModel-LG improves upon the use of a climatology by evolving the snow depth and density in response to changing atmospheric and snowpack conditions.Our results lend confidence that this daily snow depth and density product is able to represent the observed spatial and seasonal variability in snow distribution as it captures many of the key features shown in W99 in regard to regions of high/low accumulation and seasonal cycle.However, there are some key differences compared to W99: maximum snow depth is generally reached in April or May, 2 months later than W99, and snow depth is overall declining over time.This is not surprising as the time period over which snow can accumulate on the ice has declined as the melt season length has lengthened.For the Arctic as a whole, the annual snow volume, which integrates snow depth and sea ice extent, shows a steep decline since 1980, at a rate of −8.5% and −5.1% per decade for MERRA-2 and ERA5, respectively.
In addition to the long-term trend in snow accumulation, there is also a large amount of interannual variability, which indicates that it is inappropriate to use snow depth climatologies for assessing interannual variability in ice thickness from satellites such as CryoSat-2 (e.g., Bunzel et al., 2018;Kern et al., 2015).Interannual variability in the modeled snow depths between the two reanalyses is strongly correlated (r=0.86),leading to confidence that our interannual snow depth variations accurately reflect true precipitation variability.Differences in precipitation from 1 year to the next can lead to snow depth differences between the W99 climatology and SnowModel-LG on the order of 40 cm for individual grid cells.A 40-cm difference in snow depth, from, say, 20 to 60 cm would translate into an 81-cm difference in CryoSat-2-derived ice thickness, assuming a freeboard of 30 cm and densities of sea water, snow, and ice of 1023.9, 287.5, and 882.0 kg m −3 , respectively.This uncertainty propagation assumes a perfect knowledge of the sea ice freeboard.Further, given that SnowModel-LG simulates negative snow depth trends over time, with some regions showing a 16-cm decline since 1980, this can result in misinterpretation of long-term thickness changes when applying the W99 climatology across multiple satellite altimeter systems.
Several data sets exist for which to evaluate SnowModel-LG results, yet a key challenge remains their limited spatial and temporal coverage, as well as the spatial representativeness of the 25-km SnowModel-LG grid cell with point-based measurements.An important gap remains the lack of seasonally varying snow depths and densities over all ice types, something which the MOSAiC campaign can help remedy.Further, while the original version of SnowModel adapted for sea ice was rigorously tested during N-ICE (Liston et al., 2018), it has yet to be validated during summer melt.However, even with more in situ data, modeled snow depth and density at 25-km grid cells are not expected to match the highly variable snow depths in point measurements.SHEBA, N-ICE, and IMBs, while the most accurate snow depth measurements available, measure snow depth at meter to kilometer spatial scales compared to the 25-km grid cell of SnowModel-LG.OIB features larger spatial coverage, but it is clear that large differences between OIB snow products exist, and the spatial coverage is still small compared to 25-km grid cells.Satellite observations of snow depth offer pan-Arctic and longer time series perspective at similar spatial resolution, an improvement over the spatial and temporal limitations of field observations, yet have mostly been limited to FYI regions and do not cover the summer melt season.Using radar systems at two different frequencies may provide a way forward for pan-Arctic snow mapping during winter (e.g., Guerreiro et al., 2016;Lawrence et al., 2018), but this strongly depends on scattering horizon assumptions, which may not be valid for seasonally evolving snowpacks, especially over saline snowpacks (e.g., Nandan et al., 2017).With the launch of ICESat-2, a new opportunity for obtaining pan-Arctic snow depth products in conjunction with CryoSat-2 freeboards may also be possible during the cold season (Kwok & Markus, 2018).

Conclusions
This study evaluates the first full seasonal-cycle pan-Arctic, daily varying snow depth, and snow density estimates over Arctic sea ice from 1980 to 2018 using a sophisticated new Lagrangian snow pack model developed for sea ice applications.We test our methodology using two different reanalysis systems, each with different estimates of total precipitation.While SnowModel-LG was tested and bias-corrected within an 10.1029/2019JC015900 Journal of Geophysical Research: Oceans assimilation framework using OIB data, here we also tested outside of the assimilation framework against several OIB snow depth products, IMBs, SHEBA, and N-ICE in situ measurements, and passive microwave data.Each of these data sets provide one representation of snow; SnowModel-LG provides yet another.Our results show reasonable agreement with IMBs, SHEBA, and N-ICE measurements, especially considering representation issues between point measurements and the 25-km grid cells we consider here.On the other hand, our results are not in agreement with those from passive microwave or from climatology.Negative snow depth trends are observed since 1981 throughout most of the Arctic during the cold season months, rendering climatologies based on data collected several decades ago no longer valid in today's Arctic.
Although snowfall rates relate to snow depth, the accumulation of snow over sea ice using SnowModel-LG combines the effects of snowfall, snow redistribution, snow metamorphosis, and ice motion.By using this advanced snow modeling tool, we have seen two things in these simulations: there are important physical processes at work over sea ice that previously published studies have never accounted for before (i.e., blowing snow sublimation), and there are spatial structures in snow accumulation over sea ice that previous studies were not able to reproduce.In particular, compared to other new reanalysis-based snow depth products, the implementation of the higher spatial resolution Lagrangian tracking of individual ice parcels together with improved snow evolution physics representations results in more spatial and temporal variability.This allows for SnowModel-LG to capture known spatial variability in snow accumulation and snowpack evolution when sea ice dynamical redistribution of the snow cover in response to the ice motion is accounted for.The ice dynamics is crucial to capture the observed snow depth spatial variability, serving to increase the spatial variability of snow on sea ice in the Arctic.This is true using either MERRA-2 or ERA5, and despite differences in total annual precipitation, resulting snow depths using either MERRA-2 or ERA5 are within 5 cm of each other for most regions of the Arctic, when the calibration scaling is applied.Nevertheless, ERA5 precipitation was found to perform slightly better overall compared with North Pole drifting station data (Barrett et al., 2020), and its native spatial resolution is higher, suggesting it is a suitable product for producing pan-Arctic snow depth estimates going forward.
Finally, while we have compiled as much snow information as possible, the Arctic Ocean remains poorly observed; this is an unfortunate reality.The scale of our problem, using 25-km grid cells, allows smoothing of high-resolution noise when measurements are available, such as that from OIB, though grid cell representativeness error warrants further study.Further, all observations have their limitations and uncertainties, and different spatial and temporal scale representations.Ideally, a distributed snow depth network would allow for an optimal data assimilation process to aggregate biases from several locations, increasing confidence in basin-wide snow depth retrievals from models.
Pl e a s e n o t e: C h a n g e s m a d e a s a r e s ul t of p u blis hi n g p r o c e s s e s s u c h a s c o py-e di ti n g, fo r m a t ti n g a n d p a g e n u m b e r s m a y n o t b e r efl e c t e d in t his ve r sio n.Fo r t h e d efi nitiv e ve r sio n of t hi s p u blic a tio n, pl e a s e r ef e r t o t h e p u blis h e d s o u r c e.You a r e a d vis e d t o c o n s ul t t h e p u blis h e r's v e r sio n if yo u wi s h t o cit e t hi s p a p er. Thi s v e r sio n is b ei n g m a d e a v ail a bl e in a c c o r d a n c e wit h p u blis h e r p olici e s.

Figure 1 .
Figure 1.Location of buoys from the U.S. Army Cold Regions Research and Engineering Laboratory and snow buoys from the Alfred Wegener Institute, referred to here as buoys, as well as the location of the Norwegian young sea ICE (N-ICE) and the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment.

Figure 4 .
Figure 4. SnowModel-LG simulated snow depth evolution during N-ICE2015, using MERRA-2 (blue) and ERA-5 (green) atmospheric forcing.The Lagrangian parcel was identified that corresponded to the starting position of each N-ICE2015 observing leg (each leg was on a different parcel), and the entire years' snow depth time series was extracted and plotted for each parcel.The observed snow depth average (red X) and plus/minus one standard deviation (vertical black line with red dots) are shown for data collected during the gray observation periods during N-ICE2015 (a) Leg-1 and (b) Leg-2 (Liston et al., 2018; Merkouriadi et al., 2017) and (c) Leg-3 and (d) Leg-4 (Gallet et al., 2017).

Figure 5 .
Figure 5. Scatterplot of buoy snow depth versus SnowModel-LG for (a) MERRA-2 and (b) ERA-5.Black circles correspond to winter values, whereas red circles represent melting snow depths.R 2 value is given for the winter (black circle) snow depths.
Spatial pattern differences between modeled and passive microwave-derived snow depths over FYI in the mid-1980s (October 1985 through April 1986) and the mid-2010s over both FYI (October 2015 through April 2016) and MYI (March, April 2016) are summarized in Figure 8. Spatially, snow depths differ markedly between the two years shown and for individual months, with generally 10 to 20 cm less snow over the FYI

Figure 6 .
Figure6.OIB-derived and SnowModel-LG snow depth relative frequency histograms of 5 cm depth bins over the OIB flight lines from 2009 to 2015, gridded to the 25-km EASE grid.Mean of all OIB products is given by solid gray line.

Figure 7 .
Figure 7. Difference in mean OIB snow depth, as represented by an average over four data products (GSFC, QuickLook, JPL, and SLDR) to that from SnowModel-LG runs using the MERRA-2 reanalysis.Results are shown for each year from 2009 to 2015.

Figure 9 .
Figure 9.Time series from January 2000 to April 2018 of (top) first-year and (bottom) multiyear ice snow depths from SnowModel-LG and passive microwave.MERRA-2 is shown in green, ERA-5 in red, and passive microwave (PMW) in brown.

Figure 14 .
Figure 14.(a) Trends in snow depth from August 1980 to July from SnowModel-LG forced with (top) ERA-5 and (bottom) MERRA-2.(b) Time series of regional mean winter (October through March) and summer (April to September) snow density from SnowModel-LG forced with (top) ERA-5 and (bottom) MERRA-2.

Table 1
Mean Snow Depth and One Standard Deviation (cm) of All OIB Snow Depths From Each Data Product Averaged Along All OIB Tracks for Each Year Note.Corresponding results from SnowModel-LG forced by ERA5 and MERRA-2 are listed, as well as the average from all four OIB snow depth products (OIB all).

Table 2
Regional Mean 1981-2010 Averaged Snow Depths (cm) From SnowModel Runs Forced by ERA5 and MERRA-2 Reanalysis Together With Trends in Parentheses From 1980 to 2016 (cm decade −1 ) Blanchard-Wrigglesworth et al. (2018)bothBlanchard-Wrigglesworth et al. (2018)and the NASA Eulerian Snow on Sea Ice Model (NESOSIM) fromPetty et al. (2018)for April 2014; density results from MERRA-2 in November and April compared with NESOSIM are also provided in Figure15b.While all simulations show deeper snow north of Greenland and stretching across toward the East Siberian Sea, our simulations show considerably more spatial variability than previous approaches.Petty et al. (2018)do not follow the ice parcels through time but instead use a concentration budget approach followingHolland and Kimura (2016).This results in a good match with the SnowModel-NoMo runs, which have no ice motion included (see Part I).Blanchard-Wrigglesworth et al. (2018)rely on NSIDC weekly ice motion vectors in their Lagrangian tracking of ice parcels but do not interpolate the vectors to daily time steps and instead rely on weekly averaged snowfall from ERA5.Further, they run their model at much coarser spatial resolution (75 versus 25 km).Interestingly, both our snow depth and those fromPetty et al. (2018)show a region of less snow northeast of Greenland and less snow in the Laptev Sea compared toBlanchard-Wrigglesworth et al. (2018).