Evaluating and Enhancing Snow Compaction Process in the Noah‐MP Land Surface Model

The accuracy of snow density in land surface model (LSM) simulations impacts the accuracy of simulated terrestrial water and energy budgets. However, there has been little research that has focused on enhancing snow compaction in operationally used LSMs. A baseline snow simulation with the widely used Noah‐MP LSM systematically overestimates snow depth by 55 mm even after removing daily snow water equivalent (SWE) biases. To reduce uncertainties associated with snow compaction, we enhance the most sensitive Noah‐MP snow compaction parameter—the empirical parameter for compaction due to overburden (Cbd)—such that Cbd is calculated as a function of surface air temperature as opposed to a fixed value in the baseline simulation. This enhancement improves accuracy in simulated snow compaction across the majority of western U.S. (WUS) SNOTEL test sites (biases reduced at 88% of test sites), with modest bias reductions in cooler accumulation periods (biases reduced at 70% of test sites) and substantial improvements during warmer ablation periods (biases reduced at 99% of test sites). Relatively larger improvements during warm conditions are attributable to the default Cbd value being reasonable for cold temperatures (≤−5°C). Improvements in simulated snow depth and density with observations outside of the training sites and optimization periods support that the snow compaction enhancement is transferable in space and time. Differences between enhanced and baseline gridded simulations across the total WUS support that the enhancement can have important impacts on snowpack evolution, snow albedo feedback, and snow hydrology.


Introduction
Snow is an important component of global, regional and local climate and weather systems due to its high albedo which modulates the amount of radiation absorbed by the land surface (Barlage et al., 2010;Cohen & Entekhabi, 2001;Gao et al., 2017).Snow is also a critical component of water supply, with more than one sixth of the global population receiving water from snow-dominated regions (Barnett et al., 2005) and approximately 70% of total western United States (WUS) runoff originating as snowmelt (Kapnick et al., 2018;D. Li et al., 2017).Yet, commonly used land-surface-models (LSMs) which inform snow water management show substantial differences in simulated snowpack properties (Chen et al., 2014).Snow density-the ratio of snow water equivalent (SWE) to snow depth-is a fundamental property of snow because it relates SWE to snow depth and affects the snow-surface albedo.Snow compaction modeling in physically based LSMs plays important roles in modeling the impacts of snow on the terrestrial energy budget, particularly in some widely used LSMs (e.g., Noah-MP or the Community Land Model (CLM)) which compute snow cover area and snow albedo as functions of snow depth and density (e.g., Lawrence et al., 2019;Niu et al., 2011).Simulating snow compaction accurately is also important for modeling SWE in simulations that assimilate snow depth observations because the LSM compaction scheme directly relates ingested snow depth observations to SWE (Smyth et al., 2019(Smyth et al., , 2020)).This is critical because in situ snow depth measurements largely outnumber SWE measurements worldwide (Jonas et al., 2009;Sturm et al., 2010).Despite the importance of accurately simulating snow compaction, there is little research focused on enhancing snow compaction schemes in widely used LSMs.
Model uncertainties in snow compaction are exacerbated by the snow albedo feedback (SAF) in coupled systems, which provide climate, weather and water forecasts (Abatzoglou et al., 2021;Liu et al., 2017;Livneh & Badger, 2020;Meng et al., 2018;Siirila-Woodburn et al., 2021;Williams et al., 2022) that inform management of water resources, wildfire, drought response, and agriculture for millions of people in the western United States (WUS) (Abolafia-Rosenzweig, He, & Chen, 2022; Abolafia-Rosenzweig, He, Chen, Ikeda, et al., 2023;Gochis et al., 2015Gochis et al., , 2020;;Kapnick et al., 2018;Livneh & Badger, 2020;Mote et al., 2018;Musselman et al., 2021;Qin et al., 2020).In particular importance to this study, the SAF enhances uncertainties associated with snow compaction biases in coupled modeling systems that employ the widely used Noah-MP or CLM LSMs for operational and research applications, as these models calculate ground snow cover fraction (SCF) as a non-linear function of snow depth and density such that uncertainties in snow compaction propagate to uncertainties in SCF.When a model pixel is partially covered by snow (i.e., SCF < 1.0), the pixel-scale ground albedo is calculated as an average of the albedo associated with the snow cover and the albedo associated with the uncovered bare soil (He et al., 2023a;Lawrence et al., 2019;Niu et al., 2011).Because snow albedo is frequently 2-3 times larger than that of the underlying land surface, uncertainties in SCF can have large impacts on the modeled land surface albedo.Albedo uncertainties propagate to uncertainties in absorbed solar radiation by the land surface, which propagates to uncertainties in sensible heat flux from the land surface to the lower atmosphere, which propagates to uncertainties in near-surface temperature (e.g., cold or warm biases).Overall, these uncertainties propagate to greater ablation uncertainty and even greater uncertainties in snow depth.These feedbacks are consistent with previous analyses which found that coupled modeling systems using the Noah-MP LSM overestimates SCF and surface albedo (He et al., 2019) and underestimates lower atmospheric temperatures (i.e., cold-bias) in the WUS (Liu et al., 2017;Rasmussen et al., 2023).
LSM snow depth accuracy is dependent on the accuracy of SWE and snow compaction (Dawson et al., 2017;Smyth et al., 2019Smyth et al., , 2020; Figure 1).In recent years, a large body of research has worked toward improving model SWE accuracy by accounting for previously unmodeled processes, parameter optimization, reducing biases in meteorological boundary conditions and assimilation of in situ and satellite-observed snow states (Abolafia-Rosenzweig et al., 2021;Abolafia-Rosenzweig, He, McKenzie Skiles, et al., 2022;Broxton et al., 2016;Chen et al., 2014;He et al., 2019;Lim et al., 2022;Liu et al., 2017;Livneh et al., 2013;Oaida et al., 2019;Rasmussen et al., 2023;Smyth et al., 2019Smyth et al., , 2020;;Wang et al., 2019) and a large amount of on-going research continues to target improving LSM SWE accuracy.However, there is little research focused on improving model snow depth accuracy directly through enhancing snow compaction schemes.Previous research quantified deficiencies in the snow density formulation used by the Noah LSM (which employs a single-layer snow model) and proposed an enhanced formulation informed by Snow Telemetry (SNOTEL) snowpack observations (Dawson et al., 2017).However, little to no research has attempted to enhance LSM snow compaction in more comprehensive LSMs which use multi-layer snow models (e.g., Noah-MP).This gap in research is important to address because widely used operational and research application modeling systems leverage such LSMs (Powers et al., 2017;Gochis et al., 2020;C. Li et al., 2017;Reick et al., 2021;Kumar et al., 2006;Peters-Lidard et al., 2007).
Therefore, the primary focus of this research is to enhance the snow compaction scheme in the widely used Noah-MP LSM.Noah-MP computes changes in snow density using empirically based parameters which were derived over four decades ago from an analysis that used relatively small records of data (2 locations with <5-year of data) (Anderson, 1976;Kojima, 1967).This study evaluates and enhances the skill of Noah-MP snow compaction by comparing LSM simulations with multiple decades of data from >700 snow observing stations in the WUS (USDA Natural Resources Conservation Service, 2022).Our overarching hypothesis is that the long-standing empirically based snow densification parameterization used by Noah-MP can be enhanced using the large number (>2,000,000) of observations that are available from the SNOTEL network.To address this hypothesis, our study answers three primary science questions.(a) How accurate is the current snow compaction scheme in Noah-MP?(b) Can the Noah-MP snow compaction parameterization be enhanced to improve agreement with observed snow depth and density?(c) Are the enhancements transferable in space and time?

Noah-MP Simulations
The Noah-MP LSM, developed based on Noah-v3.0 (Chen & Dudhia, 2001;Chen et al., 1996Chen et al., , 1997;;Ek et al., 2003), contains multiple physics options to calculate water and energy processes (He et al., 2023b;Niu et al., 2011;Yang et al., 2011).Noah-MP considers a multi-layer snowpack depending on snow depth, which allows it to simulate snow variables more accurately than Noah.Noah-MP accounts for explicit snow layers when snow depth exceeds 25 mm.Before this threshold, no explicit snow layer exists and the very thin snowpack is combined with the top soil layer.The SCF on the ground is parameterized as a function of snow depth, ground roughness length and snow density (Niu & Yang, 2007).The scheme results in a higher SCF for the same snow depth during accumulation periods when snow densities are relatively low, compared to ablation periods when snow is relatively dense.SCF is used to compute a weighted average of albedos of snow and soil surface to derive the pixel-wise ground albedo and subsequently the ground absorbed solar radiation (Jiang et al., 2020).The Noah-MP snow compaction process is described in detail in Section 2.2 below.Other details of Noah-MP snowpack treatment and water and energy balances can be found in Niu et al. (2011) and He et al. (2023a).
Simulations in this study use the model-physics options consistent with the Weather Research and Forecasting (WRF)-Noah-MP options used in NCAR's (National Center for Atmospheric Research) continental-scale convection-permitting regional climate simulations (He et al., 2019(He et al., , 2021;;Liu et al., 2017;Rasmussen et al., 2023).Noah-MP snow-related parameters follow the values used in the latest release of Noah-MP version 4.5 (https://github.com/NCAR/noahmp/tree/release-v4.5-WRF),where the SCF-formulation parameters have been updated with land-cover dependency to improve simulated surface albedo and temperature (He et al., 2021).Leaf and stem area indices (LAI and SAI) are classified by vegetation type based on the Moderate Resolution Imaging Spectroradiometer (MODIS) monthly climatology from 2000 to 2008 (Yang et al., 2011).The focus of this study is to enhance the snow compaction scheme.Simulations using the default snow compaction parameterization are termed "baseline" simulations, and simulations using the improved compaction scheme are termed "enhanced" simulations (see Sections 2.2 and 2.4 for details on default and enhanced schemes, respectively).
The first suite of simulations is generated over 804 SNOTEL stations across the WUS over the 2009-2018 water years (WYs) to enhance, optimize and evaluate the Noah-MP snow compaction parameterization (detailed in Section 2.4).Additional simulations are conducted to test the transferability of compaction enhancements in time over the same spatial domain but for the prior decade (WYs 1999(WYs -2008)).Both sets of simulations are constrained with observed SWE to reduce uncertainties associated with the snow water budget (see Section 2.3 for details) and use a 10-year spin-up loop.These simulations are forced with humidity, wind velocity, pressure, and downward solar and longwave radiation from downscaled hourly North American Land Data Assimilation system version 2 (NLDAS-2) (Xia, Ek, et al., 2012;Xia, Mitchell, et al., 2012) data to 90-m resolution with topographic adjustment following Liston and Elder (2006) and Gupta and Tarboton (2016).Precipitation and temperature are constrained to equal daily observations at SNOTEL sites following He et al. (2021) to reduce uncertainties associated with meteorological forcing.Gridded open loop simulations (i.e., no data assimilation) across the entire WUS are generated to quantify broadscale impacts of the enhancements compared to the baseline parameterization.These gridded simulations are forced by hourly meteorological fields from a recent USGS-NCAR 4-km long-term hydroclimate reanalysis covering the continental US (Rasmussen et al., 2023).In all simulations, vegetation cover and elevation are based on 30-m data from the National Land Cover Database (NLCD) (Wickham et al., 2021) and the Shuttle Radar Topography Mission (https://doi.org/10.5066/F7K072R7),respectively.

Snow Compaction Formulation in Noah-MP
The Noah-MP snow compaction formulation is based on Anderson (1976) and described in detail in He et al. (2023a).The total fractional change in snow depth due to compaction (∆F DZ,sno ) is computed as the summation of three factors for each snow layer when total snow depth ≥25 mm (i.e., explicit snow layers form): (a) destructive or equi-temperature metamorphism (D DZ,dm ; s 1 ), (b) compaction due to the weight of the overlying layers of snow (i.e., overburden; D DZ,bd ; s 1 ), and (c) melt metamorphism (D DZ,melt ; s 1 ) as following: where ∆t is the model timestep and i is the snow layer.The 0.5 constraint assumes it is unreasonable for snow density to more than double in a single model timestep.D DZ,dm for snow layer i is computed as: where T frz and TSNO i are the melting point of snow/ice (273.16K) and snow temperature for layer i (units of K), respectively.C dm,1 and C dm,2 are empirically based parameters with default values of 2.5 × 10 6 s 1 and 0.04 K 1 , respectively.An upper limit of snow density for destructive metamorphism compaction (D M,max , kg/ m 3 ) is applied so that D DZ,dm,i is updated using Equation 3 in the case that the partial density of ice to total water (W ice,sno,i /D zsno,i ) exceeds D M,max : where D DZ,dm on the right-hand-side of the equation comes from Equation 1, W ice,sno,i [mm] is the ice content in snow layer i, and D zsno,i is the depth of snow layer i.In the case that the depth equivalent for snowpack liquid water (using ρ water = 1,000 kg/m 3 ) exceeds 1% of the snow depth for layer i, D DZ,dm,i is updated using: where C dm,3 is an empirically based parameter (≥1.0), which accounts for increased compaction under conditions of high-water content in the snowpack.
D DZ,bd is computed as: where W abv is the pressure of overlying snow (kg/m 2 ) and W sno,i is the cumulative ice and liquid water mass for layer i. C bd is an empirically based parameter assigned a default value of 0.021 m 3 /kg.η 0 is a viscosity coefficient at 0°C which had a default value of 0.8 × 10 6 kg s m 2 , which was recently increased to 1.33 × 10 6 kg s m 2 in He et al. (2021) to better match SNOTEL observed snow depth.
D DZ,melt is computed as: where F ice,old,i is the fraction of ice to total water in the snow from the previous time step (kg/kg) and F ice,i is this fraction for the current time step for snow layer i.
The thickness for snow layer i is updated using: The most up-to-date default values prior to this study for compaction parameters are summarized in Table 1.

Implementation of SWE Direct Insertion
83% of the variability in biases of baseline simulated snow depth are explained by biases in simulated SWE, relative to SNOTEL observations (Figure 1).Therefore, we consider it necessary to constrain model SWE with observations prior to optimizing and evaluating the compaction scheme to reduce uncertainties, which may affect snow compaction optimization that are not pertaining to compaction processes (e.g., reducing compensatory errors).To do this, we employ the direct insertion data assimilation technique, which updates model SWE to equal observed SWE when valid daily SNOTEL observations are available.In this procedure, first the difference between simulated (SWE mod ) and observed SWE (SWE obs ) are calculated at the time step when SWE is observed: Then, the snowfall rate for this timestep (Q snow ) is set equal to SWE increment (i.e., Q snow = SWE increment /DT).
When SWE increment is positive, Q snow is assumed to have the bulk density of fresh snowfall (BD fall ) and is imposed on the top snow layer as additional ice content: where h sno,1 is the depth of snow layer 1.When SWE increment is negative, snowpack reductions are made from top to bottom imposing non-negativity constraints for each snow layer and translating the negative SWE increment downwards through the snow column by removing the existing SWE from each layer.In this case, the original (i.e., pre-data assimilation) snow densities and partial ice and water ratios for each layer are maintained post-hoc.
In the case that SWE increment removes the entire snowpack (i.e., SWE increment < 0 and |SWE increment | > SWE mod ), a non-negativity constraint is applied and used in the Noah-MP water budget check, such that the Noah-MP water budget error check does not need to be relaxed.All model simulations over SNOTEL stations evaluated herein between the open loop simulated SWE biases and snow depth biases from the direct insertion simulation, indicating the original model SWE biases do not exert control on snow depth biases in the direct insertion simulation.However, maintaining the pre-direct-insertion snow density for each snow layer post-hoc in the direct insertion procedure maintains a portion of the snow depth and density uncertainties that are attributable to inaccurate pre-directinsertion SWE.In an ideal data assimilation experiment, snow density errors attributable to inaccurate SWE would be updated in the data assimilation procedure.However, in practice it is difficult to determine how snow density should be updated based on the SWE analysis increment.For example, if snow is added to the terrestrial system during direct insertion of observed SWE due to a model SWE underestimate, the original model underestimate could be attributable to underestimated snowfall, overestimated vegetation interception and associated sublimation, unmodeled wind redistribution of snow, overestimated ground snow ablation or a combination of factors.Determining the root causes for the model errors is required to accurately update the model snow density in SWE data assimilation, but determining these root causes and including appropriate snow density updates in the data assimilation process is a remaining challenge for snow data assimilation research that is out of the scope of this study.Therefore, we proceed accepting that model SWE errors can impact modeled snow depth and density even in the direct insertion simulations.Thus, although our methodology successfully reduces the propagation of SWE biases to modeled snow depth biases, improvements from the enhanced model simulations presented herein may be partially compensating for snow density errors induced by uncertain model SWE at the sub-daily scale.

Parameter Sensitivity and Optimization of Snow Compaction Parameters
To determine the most impactful snow compaction parameter(s) to be optimized, we first quantify the sensitivity of all Noah-MP compaction parameters (Table 1) by conducting two simulations for each parameter: (a) increasing the default value by 50% and (b) decreasing the default value by 50%.In each of these cases, all compaction parameters are maintained at default values, except respective parameters of interest for each case.
After generating these simulations for each compaction parameter, the range of snow depth from perturbed simulations is used to quantify the sensitivity of snow compaction in Noah-MP to each parameter.This sensitivity procedure assumes that the six tunable compaction parameters have similar uncertainties up to 50%; however, if parameter uncertainties are variable, then this sensitivity test may underestimate or overestimate the sensitivity of associated compaction parameters.
We found the baseline simulation has spatially and temporally varying biases, so it is required to calculate the most sensitive snow compaction parameter (C bd ) as a function of spatiotemporally varying conditions.We aim to relate C bd to a spatiotemporally varying variable already used by the LSM to reduce complexity of enhancements (e.g., adding new variables to the LSM).We choose to only optimize the most sensitive compaction parameter for the following reasons: (a) some other parameters (e.g., η 0 ) have a physical meaning and thus have physical limits and are less empirical; (b) adjusting multiple parameters simultaneously may result in compensatory errors (i.e., the right answer for the wrong reason); and (c) the model is highly sensitive to C bd so adjusting other parameters does not provide substantial additional improvements.Additionally, the compaction enhancement is designed such that minor code updates are required to implement the enhancement in operational modeling systems.
In the optimization procedure, we first generate 40 simulations with unique C bd values ranging from 0.0105 to 0.0315 m 3 /kg ( ±50% perturbations of the default value), while maintaining all other parameters at default values.The 50% range for perturbations is selected such that the simulated snow depth range from perturbed simulations is larger than model biases, which suggests that optimizing in this range is suitable to achieve a minimum model bias (Figures S1 and S2 in Supporting Information S1).We then select two optimal static C bd values at each SNOTEL station: (a) minimizing biases relative to observed snow depth during accumulation periods (prior to April 1), and (b) minimizing biases during ablation periods (April 1 and after).We hypothesized that these optimized C bd values would have a significant relationship with surface air temperature based on the strong relationship between snow density and surface air temperature (Dawson et al., 2017;De Michele et al., 2013;Valt et al., 2018).Figure 2 supports this hypothesis, showing a strong correlation between mean surface air temperature (averaged over accumulation and ablation periods) with seasonally optimized C bd (r = 0.67, p < 0.001) over corresponding periods for instances when seasonally averaged surface air temperature is below freezing.However, there is a very weak relationship for instances when seasonal temperature is above freezing (r = 0.05).Preliminary analyses also compared simulated snow temperature with optimized static C bd values, because in nature snow densification likely responds more directly to the temperature of snow rather than the surface air temperature.However, there is large uncertainty in snow temperature simulated by Noah-MP, resulting in weaker coupling between snow temperature and optimal static C bd values (r = 0.20).Therefore, our snow compaction enhancement leverages the relationship between surface air temperature with snow densification.
We quantify relationships between optimized static C bd values and surface air temperature using seasonal averaging in Figure 2 because determining finer resolution (e.g., daily) optimized C bd values is highly uncertain due to memory in associated state variables.However, Noah-MP can only relate instantaneous C bd with surface air temperature without imposing complex model updates that would be difficult to implement in operational systems.Because the best-fit linear relationship between seasonally optimized C bd and surface air temperature in Figure 2  Given the very weak relationship between C bd with surface air temperature for instances where seasonal temperature is above freezing, we apply variable lower constraints to these evaluated linear models based on evaluations of the enhanced compaction parameterization which revealed a robust statistical relationship between seasonally optimized C bd with surface air pressure for warm instances (r = 0.64, p < 0.01; Figure S3 in Supporting Information S1).
Preliminary analyses revealed binning C bd lower constraints on the basis of surface air pressure provides substantially improved performance relative to the originally optimized simulation which uses the lower constraint of C bd,min = C bd (temperature = 273.16K) (see Text S2 in Supporting Information S1 for details).Therefore our optimization and evaluation of the enhanced compaction scheme employs surface air pressure based binning of C bd,min (Equation 12) which is intended to represent elevation-based effects on snow compaction.T air and psfc are the surface air temperature and pressure, respectively.We use surface air pressure as a proxy for elevation controls on snow compaction in this formulation because elevation is not functionally used in Noah-MP, and thus using surface air pressure (which is already used in Noah-MP) instead of elevation allows simplistic code updates which can be easily implemented in operational systems.Finally, we impose an upper constraint on C bd (0.0315) supported by Figure 4 which indicates minimal temperature control on snow density for instances when surface air temperature is less than 251 K which corresponds with C bd = 0.0315 in Equation 11.Model code updates associated with snow compaction enhancement are publicly available: https://github.com/RAbolafiaRosenzweig/HRLDAS_Snow_Compaction_enhancement/tree/main.

Evaluation Metrics
Metrics evaluating the accuracy of simulated snow depth and density in site-wise evaluations, relative to SNOTEL observations, are calculated only when observed snow depth exceeds 30 mm because the Noah-MP snow compaction scheme is not used for very thin snowpacks.The percent bias (pbias) considers the averaged relative difference between simulated and observed states.
O and P represent the time series mean of observed and simulated states, respectively.The unbiased root mean square error (ubRMSE) is calculated as: where O and P are predicted and observed state time series from respective sites.The Pearson correlation coefficient (r) (Pearson & Henrici, 1896), ranging from 1 to 1, is a widely used benchmark for performance calculated as: where P t and O t are predicted and observed states at timestep t.

Journal of Advances in Modeling Earth Systems
10.1029/2023MS003869

SNOTEL Snowpack Validation Data
In situ observations of SWE and snow depth used in this study are from the United States Department of Agriculture (USDA) Natural Resources Conservations Service (NRCS) SNOTEL network which provide data from high snow accumulation regions throughout the WUS (Serreze et al., 1999).This analysis uses biascorrected and quality controlled (QCed) daily SNOTEL observations of SWE from the Pacific Northwest National Laboratory (PNNL) (Sun et al., 2019;Yan et al., 2018; https://www.pnnl.gov/data-products).Because the non-bias-corrected snow depth measurements can be noisy at times (Dawson et al., 2017;Hill et al., 2019), we implement a quality-control strategy to only consider observations of snow depth with realistic values.Specifically, snow depth observations are screened if: (a) SWE exceeds snow depth in the same unit of mm, (b) snow depth is less than 0 mm or exceeds 50,000 mm, (c) snow depth increases by greater than 1,000 mm in a single day (removing spikiness in the data), (d) snow depth decreases by at least 100 mm in a single day but the corresponding QCed SWE observations do not decrease by at least 100th of the snow depth decrease, (e) data have a greater than 2,000 mm displacement from the best fit line comparing SWE and snow depth.Screening step (v) is implemented assuming a linear relationship between SWE and snow depth, and large outliers may be unrealistic (Hill et al., 2019).After performing this screening procedure, 8.53% of the original snow depth observations were removed leaving 2,040,166 QCed snow depth observations in the analysis period (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018).After screening, there are 781 and 622 SNOTEL stations which we use in our analyses which have more than 300 valid daily snow depth and SWE observations during optimization (2009-2018) and temporal transferability evaluation (1999-2008) periods, respectively.We recognize that this procedure is very conservative and may remove some valid observations, but we consider removing valid observations less detrimental than maintaining invalid observations, particularly because the large number of remaining observations available after screening.

Parameter Optimization
The Noah-MP snow compaction scheme is most sensitive to parameters used to calculate snow compaction from overburden (i.e., C bd and η 0 from Equation 5), with the largest sensitivity to C bd across all vegetation classifications and elevation bins (Figures S1 and S2 in Supporting Information S1).Noah-MP snow depth simulations are relatively insensitive to the remaining four snow compaction parameters used to calculate snow compaction from metamorphic changes.It is possible that these insensitive parameters may show higher sensitivity outside of the western U.S. domain where testing was performed, but these results may also indicate that the model compaction formulation may be able to be considerably simplified (i.e., by reducing the number of tunable parameters).Results from this sensitivity test are consistent with the results in He et al. (2021), where the Noah-MP snow depth is most sensitive to the overburden-induced snow compaction process.Therefore, we proceed to enhance and optimize C bd given the high sensitivity of the Noah-MP snow compaction process to this empirical parameter (Section 2.4).
Noah-MP snow compaction biases using the baseline simulation are variable in space and time (Figures 3 and 5).Therefore, we consider it more valuable and physically reasonable to update C bd as a function of spatially and temporally varying conditions rather than updating C bd to a new spatiotemporally static and uniform value.As discussed in Section 2.4, we relate C bd to surface air temperature using a linear model, and optimize the corresponding linear model coefficients.This decision is supported by a strong relationship between the seasonal mean surface air temperature with seasonally optimized C bd (r = 0.67 ∼ 0.93, p < 0.001) for instances when seasonally averaged temperature is below freezing (Figure 2; Figure S3 in Supporting Information S1).We then leverage the relationship between seasonally optimized C bd with surface air pressure for warm instances (r = 0.64 to 0.69, p < 0.01; Figure S3 in Supporting Information S1) via binning C bd lower constraints to provide further improved performance (see Section 2.4 and Text S2 in Supporting Information S1 for details).Enhanced simulations use the optimized relationships between surface air temperature and pressure with C bd as shown in Figure 2 (i.e., Equations 11 and 12).

Accuracy of Snow Compaction From Enhanced and Baseline Simulations
A dominant source of error in simulating the depth of snowpack is attributable to errors in simulating SWE (i.e., water budget errors), with the large majority of the variability in snow depth biases being explained by biases in simulated SWE (r 2 = 0.83) in the baseline open loop simulation using default parameters (Figure 1a).However, even when errors in daily SWE are removed from the Noah-MP simulation by directly inserting observed SWE (Section 2.3), biases in snow depth are still prominent in the baseline simulation: a mean systematic bias in snow depth of 55 mm (6% of observed mean snow depth) with an interquartile range (IQR) of 45-127 mm (Figure 1b).Some of these remaining biases in snow depth may be attributable to sub-daily errors in model SWE or observational uncertainty, but we proceed under the assumption that the majority of these remaining errors are attributable to uncertainties in the Noah-MP snow compaction scheme which we attempt to reduce through our optimization strategy (Section 2.4).Temporal variability in simulated snow depth from the baseline and enhanced direct insertion simulations strongly agrees with observations averaged across all SNOTEL test sites (r > 0.99) and across each of the three evaluated elevation (Z) bins (r > 0.99) (Figure 3).Across each elevation bin, the enhanced simulation provides slightly higher r (by <0.003) relative to the baseline simulation, and the enhanced simulation has notably lower biases (by 4%-11%) than the baseline simulation for low-and mid-elevation bins and with a slightly lower mean pbias (by 1%) for the high-elevation bin.Snow depth from the enhanced simulation is particularly more accurate during warm times (e.g., ablation periods) and places (e.g., low-elevation sites where the baseline simulation has a mean pbias of 13.2% and the enhanced simulation has a mean pbias of 2%) (Figure 3).These results from the test sites are similar to results from the training sites which also show consistent bias reductions (by 1%, 5%, and 11% for high-, middle-and low-elevation sites, respectively).Similar improvements between training and testing sites support the model enhancement is spatially transferable.
Figure 4 supports simulated snow density from the enhanced scheme is slightly more accurate in representing mean observed snow density during cold periods (surface air temperature <266K); whereas at warmer temperatures, the enhanced simulation is substantially more accurate at representing observed snow density.Similar results from the training and testing sites further support the successful spatial transferability of the model enhancement.Relatively larger improvements from the enhanced simulation for warm instances occurs because the default C bd parameter value is generally suitable for the relatively cool temperatures when compaction occurs slowly.This is consistent with the formulation of our enhanced C bd parameter (Equation 11), which equals the default value at ∼266.3 K ( 6.9°C).Notable biases in snow density for above-freezing instances, even in the enhanced simulation, (Figure 4) may be attributable to propagation of sub-daily SWE errors to snow density (see Section 2.3 for discussion), snow-albedo feedback uncertainty associated with declining snow coverage, the insensitivity of the Noah-MP snow compaction parameterization to melting processes and wet snow compaction which plays a larger role when surface air temperature is above freezing, and other model parameterizations external from the compaction procedure.This indicates that additional model enhancements beyond the scope of this research (e.g., to the SCF process and a higher sensitivity wet snow compaction) may be necessary to further improve Noah-MP simulated snow density.
The baseline simulation generally tends to overestimate snow depth at individual sites with a mean site r, ubRMSE, and pbias of 0.95, 135 mm and 8.3%, respectively.These baseline simulation skill scores are similar for training (0.95, 137 mm, 8.4%) and test sites (0.95, 131 mm, 8.3%), respectively.Baseline overestimates are particularly large in the Pacific-Northwest (PNW) (Figures 5a and 5b).There is also a tendency for the baseline simulation to show decreasing overestimates of snow depth moving from the western to eastern portion of the WUS with an IQR of site pbiases from 0.1% to 13.4% (Figure 5b).The enhanced simulation provides improved skill scores across the majority of sites: 81% of test sites (80% of training sites) have reduced ubRMSE and 88% of test sites (89% of training sites) have reduced pbias.The enhanced simulation tends to provide slightly higher r, with increases at 72% of test sites (69% of training sites) and a mean increase of 0.003 across all sites.The enhanced simulation reduces the mean site pbias to 4.9% across test sites (4.8% across training sites) with a smaller IQR across all sites ( 2.0%-8.7%).The mean-site snow depth ubRMSE from the enhanced simulation (125 mm) is 7% smaller than ubRMSE from the baseline simulation across all sites.Bias reductions are largest in the PNW where the baseline simulation originally had the largest errors.This is because SNOTEL sites in the PNW have relatively warm snow season temperatures (mean November-April temp.= 0.14°C) compared to the other WUS SNOTEL stations (mean November-April temp.= 2.92°C), and the baseline value for C bd is more applicable to cooler sites with relatively slow snow compaction.The maximum degradation in site pbias is an increase by 2.7% and the largest improvement is a pbias decrease by 18.7%.
Improvements associated with the enhanced snow compaction scheme are most prominent during ablation periods when temperatures are relatively warm and snow densifies relatively fast (Figure 6).This is of particular importance to water managers who rely on accurate simulations of spring and summer snowmelt.Namely, there are improvements in ubRMSE and pbias at 87% and 99% of test sites during ablation periods (i.e., post-March) (88% and 98% for training sites).The baseline simulation has a mean site pbias of 24% with an IQR of 11.8%-30.1% across all sites, and the enhanced simulation has a mean site pbias of 16.9% with an IQR of 4.4%-22.6%.(c) high-elevation sites (Z > 2,500 m) (r = 0.76-0.78),with relatively higher r from the enhanced simulation across each elevation bin (Figure 9).The enhanced simulation has a lower pbias in snow density relative to the baseline simulation across mid-and low-elevation bins (by 2.9% and 8.8%, respectively) and a slightly increased pbias at high-elevation bins (by 1.7%).Bias reductions across all test sites (by 2.0%) are mainly attributable to the enhanced simulation predicting higher snow densities in spring and summer relative to the baseline simulation, which tends to reduce biases compared to the relatively high observed snow densities during warm periods.
To quantify broadscale impacts of snow compaction enhancements, we compare spatiotemporally continuous enhanced and baseline open loop simulations across the total WUS (Figures 10 and 11).We select two representative water years in the 2009-2018 record to show comparisons for: 2011 (drought area = 27%) and 2013 (drought area = 82%) which have the smallest and largest western US drought areas (based on the United State Drought Monitor; Svoboda et al., 2002), respectively.Mean monthly SCF between enhanced and baseline simulations differ by at least 0.01 during these months at 3%-20% (varying by month) of all WUS snowy pixels (i.e., pixels where baseline SCF > 0) and up to 0.75 locally.Daily differences in SCF between baseline and enhanced simulations become more skewed in spring months, relative to the accumulation period.For instance, the 95% range in daily SCF differences (enhanced minus baseline) from December-April of 2011 is 0.04 to 0.04; whereas the 95% range in May is 0.02 to 0.00.This indicates that during the accumulation period there is spatial and temporal heterogeneity in the directional effect of snow compaction enhancements on SCF; however, in the ablation period the SCF is typically reduced because the enhanced compaction scheme simulates a relatively denser snowpack when temperatures exceed 6.9°C, which leads to lower snow depth and hence smaller SCF in Noah-MP.This is consistent with analyses presented in Sections 4.2 and 4.3.Differences between baseline and enhanced SCF propagate to differences in pixel-scale albedo (Figure 10).Mean monthly surface albedo   between baseline and enhanced simulations differ by at least 0.01 during these months at 3%-18% (varying by month) of all WUS snowy pixels and up to 0.30 locally.Considering the typical downward solar radiation of a few hundred W/m 2 , the preceding albedo difference would make substantial changes in land surface energy balance.
The aforementioned differences between baseline and enhanced SCF and albedo are attributable to differences in simulated snow depth displayed in Figure 11.Mean monthly snow depth differences between enhanced and baseline simulations reach at least 10 mm in December-May of 2011 and 2013 at 36%-97% (varying by month) of all WUS pixels where snow depth ≥100 mm.The 95% range of daily snow depth differences (enhanced minus baseline) is most evenly distributed from October-January ( 89-82 mm) and becomes increasingly skewed through spring months (May 95% range = 332-2 mm).Mean monthly SWE differences between enhanced and baseline simulations reach at least 10 mm in December-May of 2011 and 2013 at 0%-20% (varying by month) of all WUS pixels where SWE ≥100 mm.The 95% range of daily SWE differences (enhanced minus baseline) consistently skews toward lesser SWE in the enhanced simulation with the most even distributions occurring from October-January ( 9-2 mm).Differences in daily SWE are increasingly skewed through spring (May 95% range = 24-4 mm) and maintain a skewed nature in summer months (June 95% range = 25-7 mm and July 95% range = 10-0 mm).Consistency between snow depth, SCF, albedo and SWE differences support that the enhanced compaction scheme can have important impacts on snowpack evolution, hydrology, and snow albedo feedbacks, which could play a bigger role in coupled land-atmosphere systems.

Discussions
There is a suite of other methodologies we considered to enhance C bd in the Noah-MP compaction scheme.In the results shown above, we optimized m and b parameters that linearly relate C bd to surface air temperature to minimize simulated snow depth errors in direct insertion simulations that constrain simulated SWE to equal observed SWE (discussed in Section 2.4).In addition, we also considered optimizing the parameters to minimize simulated snow density errors, which provides a different set of m and b parameters compared to our selection presented in Equation 11.However, these two optimization techniques result in similar performance in terms of snow density simulations: difference in RMSE across all points in space and time is less than 0.001 (mm SWE/ mm h sno ) and the snow depth optimized solution (used herein) provides a lower snow density pbias.We also considered calculating C bd as a function of other meteorological forcing variables used by Noah-MP.We use only surface air temperature to calculate C bd for cool temperatures (below freezing) because correlations between seasonally optimized C bd with other meteorological conditions were significantly lower (r ≤ 0.48) than those with the surface air temperature for below-freezing instances (Figure 2; Figure S3 in Supporting Information S1).The two meteorological conditions that showed modest correlations with optimized C bd for below-freezing instances are downwards longwave radiation and surface pressure (below-freezing seasonal r = 0.48 and 0.38, respectively).Seasonally optimized C bd has a robust relationship with surface pressure when seasonal temperature exceeds freezing (r = 0.64), and we leverage this empirical relationship in the enhanced C bd calculation via binning lower constraints on the basis of surface air pressure (as discussed in Section 2.4).The statistical snow density-meteorological relationships used to update C bd may differ in parts of the world outside of the western U.S. study domain, and thus further testing of the model enhancement should be evaluated in regions that have distinct vegetation classifications, elevations or climates from those represented by western U.S. SNOTEL sites.
The implemented snow compaction enhancement in this study improves the accuracy of simulated snow depth and snow density across the majority of WUS SNOTEL stations.However, there are still a portion of stationsparticularly colder stations-where the baseline simulation outperforms the enhanced simulation during accumulation periods.This is because the default C bd value (0.021) tends to be well suited for cold temperatures which often coincides with high elevation locations during snow accumulation periods.Notwithstanding, enhanced simulations more frequently improve than degrade skill even at cold places and times.Instances of degraded skill from the enhanced scheme are likely partially attributable to calculating C bd as a function of surface air temperature forcing instead of simulated snow temperature which has high uncertainty (as discussed in Section 2.4).This motivates future research to evaluate and enhance Noah-MP snow temperature to reduce associated uncertainties, which could in turn improve model snow compactions processes.Although the enhanced snow compaction scheme provides the largest improvements in warm instances, there are key remaining biases in snow density from the enhanced simulation for above-freezing temperatures (Figure 4).These remaining biases are partially attributable to model parameters outside of the snow compaction procedure.For instance, we find that doubling the parameter defining the liquid water holding capacity of snowpack (from 0.03 to 0.06) tends to further improve snow density for warm instances, while providing similar performance for below-freezing instances (Figure S11 in Supporting Information S1).Because the enhanced compaction scheme provides the largest improvements at warm places and periods, we hypothesize that snow compaction enhancements developed in this study are likely to provide larger improvements than those reported herein when used in snow simulations over regions that have warmer climates than SNOTEL sites which have a mean November-May 2009-2018 surface air temperature of 269.9 K.

Conclusions
Our study supports that the long-standing empirically based snow densification parameterization used by Noah-MP can be enhanced using the large number (>2,000,000) of observations that are available from the SNOTEL network.We first quantified the uncertainty associated with snow compaction using the baseline scheme with default parameters.In the baseline simulation with direct insertion of observed SWE, which removes daily uncertainties associated with modeled SWE, snow compaction errors are still prominent with a mean systematic bias for snow depth of 55 mm (i.e., snow is less compact than observations) and an IQR of biases of 45-127 mm.The baseline simulation has reasonable accuracy in site-averaged comparisons across all SNOTEL sites (r = 0.95 and pbias = 8.3%), but site-by-site analyses reveal the largest errors in the baseline simulation occurs at sites located in the relatively warmer western portion of the domain, particularly in the PNW.
In the enhancement procedure, we identified and optimized a robust linear equation relating the most sensitive compaction parameter (C bd ) to surface air temperature for below-freezing conditions.For above-freezing conditions, we binned the lower constraints of C bd derived from its relationship with temperature based on surface air pressure, supported by the robust relationship between surface air pressure and optimized C bd for above-freezing instances (r = 0.64; p < 0.01).This enhancement-designed to more accurately reflect controls of temperature and elevation on snow compaction-resulted in improved simulation of snow depth and snow density (reduced pbias and ubRMSE) at the majority of western U.S. SNOTEL test sites (81%-88% of sites) when evaluated across the total optimization period (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018).The enhanced simulation tends to be similarly skilled with the baseline simulation during cooler snow accumulation periods while providing substantially higher accuracy during warmer ablation periods (improved pbias and ubRMSE at 99% and 87% of test sites, respectively).Consistently higher accuracy from the enhanced simulation during ablation periods is attributable to the tendency of the default C bd value being reasonable only for cold temperatures (< 5°C).Similar results for comparisons over training sites and test sites across the 2009-2018 water years (when the scheme was optimized) compared to analyses over the 1999-2008 water years (independent validation years) support that the enhancement presented herein is both spatially and temporally transferable.Further comparisons between the baseline and enhanced broadscale simulations support that the enhanced compaction scheme can have important impacts on snowpack evolution, snow hydrology, and snow albedo feedbacks.Remaining snow density biases during warm periods may indicate that errors are exacerbated by snow-albedo-feedbacks attributable to uncertain snow cover fraction.The enhanced snow compaction scheme has the potential to improve snow albedo with the integration of a physics-based snowpack radiative transfer scheme in Noah-MP which explicitly considers the relationship between snow depth and albedo (Flanner et al., 2007(Flanner et al., , 2021)).

Figure 1 .
Figure 1.Simulated biases in snow water equivalent (SWE) compared with snow depth biases from (a) the baseline open loop simulation and (b) the baseline direct insertion simulation relative to SNOTEL observations.
is sub-optimal for relating instantaneous C bd to surface air temperature we evaluate a suit of 160 Noah-MP simulations with unique m and b parameters (i.e., C bd = T air × m + b) that relate instantaneous C bd to surface air temperature (T air ).The 160 sets of unique m and b parameters are generated by perturbing the m and b parameters from the best-fit the relationship in Figure 2, while constraining m and b to provide reasonable ranges of C bd (0.0105-0.0315) across a normal range of below-freezing temperatures (263-273.15K).

Figure 2 .
Figure 2. Scatter plot comparing the multiyear (2009-2018 water years) seasonally averaged surface air temperature with optimized C bd parameters from Noah-MP simulations.Two dots for each station correspond to accumulation and ablation periods, respectively.Data corresponding to below and above freezing temperatures are represented by blue and red dots, respectively.160 gray lines represent perturbations of the best-fit linear equation for the blue data points with lower constraints for above freezing instances which were evaluated in the optimization procedure.Bold orange, purple and black lines represent the optimized equations for C bd across surface pressure bins (Equations 11 and 12).Correlations (r) for the below-and above-freezing data points are in blue and red, respectively.The horizontal black dashed line represents the default C bd value (0.021).

Figure 3 .
Figure 3. Mean multiyear (2009-2018 water years) time series for snow depth averaged across (a) all 236 test sites, (b) low-, (c) mid-and (d) high-elevation test sites.Shading represents one standard deviation, calculated from the spread in daily snow depth across SNOTEL validation sites (i.e., spatial heterogeneity).

Figure 4 .
Figure 4. Mean snow density from observations (black), enhanced (red) and baseline (blue) simulations averaged by 1K surface air temperature bins at all days with valid observations at western U.S. (a) training and (b) testing SNOTEL sites in the 2009-2018 years.Shading represents respective interquartile ranges.The reference vertical black line marks the freezing temperature.

Figure 5 .
Figure 5. Skill scores (a) ubRMSE, (b) pbias and (c) r from the baseline simulation at SNOTEL sites on the top row and (d-f) changes in skill scores (enhanced-baseline) on the bottom row.Skill scores are evaluated across the entire evaluation period: water years 2009-2018.(e) and (f) present changes in the absolute values of pbias and r at each site.Small (large) dots correspond with training (testing) sites.

Figure 6 .
Figure 6.Same as Figure 5, but for ablation periods only.Skill scores (a) ubRMSE, (b) pbias and (c) r from the baseline simulation at SNOTEL evaluation sites on the top row and (d-f) changes in skill scores (enhanced-baseline) on the bottom row.

Figure 8 .
Figure 8. Skill scores (a) ubRMSE, (b) pbias and (c) r from the baseline simulation at SNOTEL evaluation sites on the top row and (d-f) changes in skill scores (enhanced-baseline) on the bottom row.Skill scores are evaluated across the entire temporal transferability evaluation period: water years 1999-2008.(e) and (f) present changes in the absolute values of pbias and r at each site.Small (large) dots correspond with training (testing) sites.

Figure 9 .
Figure 9. Mean multiyear (2009-2018 water years) time series for snow density from observations and open loop simulations averaged across (a) all 236 test sites with validation data, (b) low-, (c) mid-and (d) high-elevation test sites.Shading represents one standard deviation, calculated from the spread in daily snow depth across SNOTEL validation sites (i.e., spatial heterogeneity).

Figure 10 .
Figure 10.Differences in the mean monthly snow cover fraction (SCF) and pixel-scale albedo between baseline and enhanced 4-km simulations during representative 2011 and 2013 water years.

Figure 11 .
Figure 11.Differences in the mean monthly snow depth and snow water equivalent (SWE) between baseline and enhanced simulations during representative 2011 and 2013 water years.

Table 1
Summary of Compaction Parameters in the Noah-MP LSM

Journal of Advances in Modeling Earth Systems 10
.1029/2023MS003869 apply this direct insertion of SWE unless stated otherwise where simulations that do not apply data assimilation are referred to as open loop simulations.Model code updates associated with direct insertion are publicly available: https://github.com/RAbolafiaRosenzweig/NoahMP_SWE_DirectInsertion_Updates.Constraining model SWE with direct insertion substantially reduces snow depth biases (mean |bias| for open loop and direct insertion baseline simulations = 249 and 124 mm, respectively) and results in a decoupling