Enabling Advanced Snow Physics Within Land Surface Models Through an Interoperable Model‐Physics Coupling Framework

Accurate estimation of snow accumulation and melt is a critical part of decision‐making in snow‐dominated watersheds. In this study, we demonstrate a flexible methodology to couple a detailed snow model, Crocus, separately to two different land surface models (LSMs), Noah‐MP and Noah. The original LSMs and the coupled models (Noah‐MP‐Crocus and Noah‐Crocus) are used to simulate snow depth, snow water equivalent, and other water and energy states and fluxes. The results of simulations are compared against a wide range of independent gridded and point scale reference data sets. Our results show that coupling the detailed snow model, Crocus, with the LSMs improves the snow depth and snow water equivalent relative to independent observations. Overall, larger improvements are obtained with coupling Crocus to the Noah LSM, with the coupled Noah‐Crocus configuration reducing the RMSE and bias of snow depth from 2% to 12% and 57% to 75%, respectively, relative to Snow Data Assimilation System (SNODAS) and snow product from the University of Arizona. On the other hand, smaller improvements are obtained by coupling Crocus with Noah‐MP. The Coupled Noah‐MP‐Crocus reduces the snow depth bias but slightly degrades the RMSE of snow depth and snow water equivalent. The corresponding impacts in other water budget terms such as evapotranspiration, soil moisture, and streamflow, however, are mixed, pointing to the significant need to improve the coupling assumptions of these processes within land models. Overall, the interoperable coupling framework demonstrated here offers the opportunity to include more detailed snow physics and processes, and to advance data assimilation systems through improved exploitation of information from snow remote sensing instruments.


Introduction and Motivation
Meltwater from seasonal snow in snow-dominated basins is the major source of freshwater for billions of people and plays an important role in the socio-economic development of many countries (e.g., Barnett et al., 2005;Sturm et al., 2017;Viviroli et al., 2007).Despite the importance of seasonal snow, there are significant uncertainties in our understanding of the quantity of water stored in snowpacks, spatial and temporal variability of seasonal snow (Margulis et al., 2015), and the temporal evolution of seasonal snowpacks.Though snow measurements from in situ networks are available in many regions of the world, large-scale spatially distributed characterization of seasonal snowpacks remains challenging due to the low sampling density of the in situ observations.Remote sensing data, due to its spatial coverage, can partially mitigate the challenges posed by sparse in situ measurements of the surface layer variables.However, it's important to note that certain surface mass fluxes, such as evaporation and runoff, cannot be directly measured through remote sensing.Numerical snow models can be leveraged to gain insights into the temporal and spatial distribution of snow, and the evolution of snowpack over time.
Snowpack models, based on the level of complexity, can be classified into three main categories (Vionnet et al., 2012): one-layer snow schemes, multi-layer (e.g., 3 to 5 layers) snow schemes, and detailed snow schemes.The first category is commonly used in numerical weather prediction (NWP) models, and regional circulation models (RCMs) or general circulation models (GCMs).These models (e.g., the simplified Simple Biosphere Model (SSiB; Xue et al., 1991)) use a snow scheme to provide boundary conditions (e.g., albedo) for the atmospheric model.The second category of snow schemes is commonly used in land surface models (LSMs) to characterize mass and energy fluxes to and from the snowpack.These snow schemes account for some snow processes such as snow melt-refreeze, snow settlement, and meltwater runoff, among others.Snow schemes such as Snow 17 (Anderson, 1976), the Community Land Model (CLM; Oleson et al., 2010), and the Joint UK Land Environment Simulator (JULES; Best et al., 2011) are some examples of these intermediate complexity snow physics.The recent iteration of intermediate complexity snow schemes uses more snow layers (e.g., the Explicit Snow (ES) scheme (Decharme et al., 2016)).These schemes are recently implemented in NWP systems, for example, at the European Centre for Medium-Range Weather Forecasts (ECMWF; Arduini et al., 2019) and for GCMs (Decharme et al., 2019).Snow schemes in the third category are complex mass and energy balance models that employ snow metamorphism laws to characterize snow microstructures to derive the temporal evolution of snowpack and its mechanical (e.g., snow stratification and shear stresses) and thermal (e.g., albedo and emissivity) properties of snowpack.Well-known detailed snowpack models in this category include Crocus (Brun et al., 1989(Brun et al., , 1992)), the Snow Thermal Model (SNTHERM; Jordan, 1991), the SnowModel (Liston & Elder, 2006), and SNOWPACK (Bartelt & Lehning, 2002).
Though detailed snowpack models have been developed in the past decades, their significant computational cost has prevented the use of these models over large domains.There are limited operational or routine applications of detailed snow models for hydrological applications and avalanche hazard forecasting (e.g., Lafaysse et al., 2013;Langlois et al., 2009;Mortimer et al., 2020).Furthermore, these operational simulations are temporally and spatially limited.To develop large-scale estimates of seasonal snow, many operational systems employ approaches such as data assimilation to ingest observational inputs of snow cover, snow depth (SD) and snow water equivalent (SWE) from ground, airborne, and space-borne platforms (Carroll et al., 2001;Drusch et al., 2004;Takala et al., 2011).However, there is significant uncertainty in SWE and SD estimates from these analysis systems due to the lack of detailed snow processes (e.g., Barlage et al., 2010;Broxton et al., 2016), pointing to the need for detailed snow schemes to better represent the snow processes in large-scale simulations.
Despite the need for reliable snow physics, there have only been a few efforts to implement detailed snow schemes in RCMs and LSMs.For example, Météo-France has implemented all three types of snow schemes in the SURFace EXternalisée (SURFEX) platform (Vionnet et al., 2012).The Modèle Atmosphérique Régional (MAR) RCM (Gallée & Duynkerke, 1997;Gallée & Schayes, 1994) also uses an older version of Crocus to improve snow atmosphere interactions.In this paper, we describe an interoperable implementation of a detailed snow model within a land surface modeling framework.The development of such a system combines the use of more advanced snow physics from the snow model with the water, energy, and carbon cycle representation available within the LSMs.The system also provides detailed information about the temporal evolution of the snow microstructure and the vertical profile of snow density and temperature, essentially replacing the simpler snow physics available within the LSM.The impact of the detailed snow model on mass and energy fluxes can then be extended for land-atmosphere modeling, though not explored in this paper.We used the NASA Land Information System (LIS; https://lis.gsfc.nasa.gov/;Kumar et al., 2006;Peters-Lidard et al., 2007) as the platform to implement the interfaces between the snow model and LSMs.LIS is a widely used open-source land surface modeling and data assimilation framework, which encompasses several land models.In addition, the framework includes support for a large suite of meteorological and remote sensing data for driving the models and assimilation to constrain their estimates.LIS is also used in both research and operational settings, This paper is arranged as follows: Section 2 contains a description of the models and methods used in this work.The experimental design is given in Section 3. The results and evaluation of the proposed methodology are discussed in Section 4. Finally, key conclusions and future research directions are reported in Section 5.

The NASA/GSFC Land Information System Framework (LISF)
LIS is a land surface modeling infrastructure designed to facilitate the efficient utilization of terrestrial hydrological observations.LIS is designed in a flexible manner to support model simulations over user-specified regional or global domains using a suite of LSMs in a computationally scalable manner.A key feature of LIS is the availability of data assimilation and inverse modeling tools that facilitate the integration of information from remote sensing observations with the LSMs (LIS-DA; Kumar et al., 2008).LIS also includes a separate data preprocessing system known as the Land surface Data Toolkit (LDT; Arsenault et al., 2018), and a post-processing data analytics environment known as the Land surface Verification Toolkit (LVT; Kumar et al., 2012).LDT includes the capabilities for processing land parameter data from a variety of data sources, preprocessing support for data assimilation and uncertainty estimation, bias correction, and downscaling of forcing data, among others.LVT enables the comparison of model estimates against a large suite of in situ, remote sensing, and other model analysis products.LIS has been used to develop several custom Land Data Assimilation System (LDAS) instances, including the Global LDAS (GLDAS: Rodell et al., 2004), North American LDAS (NLDAS; Xia et al., 2012), Famine Early Warning Systems Network (FEWSNET) LDAS (FLDAS;McNally et al., 2017), Western U.S. LDAS (WLDAS; Li et al., 2021;Erlingis et al., 2021).Additionally, LIS enables the U.S. Air Force 557th Weather Wing's (557 WW) near real-time, operational global LDAS.The development of Crocus coupling within LIS allows the possibility of transitioning its use to these near real-time and end-user focused LDAS environments.

Noah-3.9 and Noah-MP-4.0.1
The Noah-3.9 LSM (Ek et al., 2003) employs a single-layer bulk snow representation, with the SWE calculated using a bulk snow density assumption (Jonas et al., 2009).The heat exchanges at the snow-atmosphere and snowsoil interfaces and processes of snow melt, sublimation, and accumulation are also represented in the model (Suzuki & Zupanski, 2018).The determination of rain-snow occurrence in the model is based on whether the air temperature is above or below 0°C.Prior studies have documented limitations of the simple Noah snow physics, including early snow melt and underestimation of ground heat flux (Niu et al., 2011).
The Noah-MP-4.0.1 model (Niu et al., 2011;Yang et al., 2011) was developed by the incorporation of multiple and new physics capabilities into the Noah LSM.The physics enhancements in Noah-MP include a multilayer snowpack, and multiple options for surface water infiltration, runoff, and groundwater, including the representation of an unconfined water table depth (Niu et al., 2007), among numerous other options.Noah-MP uses a 3-layer snow structure with minimum thicknesses of 0.045, 0.05, and 0.2 m.The number of active snow layers in the model grows and shrinks based on these minimum thicknesses.Noah-MP also includes a prognostic snow density model (Anderson, 1976) by accounting for factors such as snow compaction and melt metamorphism.An explicit snow interception model (Niu & Yang, 2004) calculates the interception of snowfall by the canopy and sublimation from it.The Jordan scheme (Jordan, 1991) is used to calculate the rain-snowfall partitioning.The other improvements include parameterization schemes for other processes (e.g., snow surface albedo, snow cover fraction, surface turbulent exchange coefficient, liquid water retention and refreezing, and snow surface roughness length) (Niu et al., 2011).Despite the development of a multilayer snowpack module in Noah-MP, a recent study (Cho et al., 2022) shows that Noah-MP significantly underestimates the SD and SWE relative to independent observations.Furthermore, the enhanced snowpack module does not provide information about the microstructure of snowpack.More details of the snow physics of Noah-MP can be found in Niu et al., 2011;Yang et al., 2011.metamorphism module, a snow discretization module, and an integrated surface albedo module.Crocus considers changes in liquid water content due to snowmelt, refreezing, and evaporation during a model time step.Determining the maximum number of layers is a tradeoff between computational cost and accuracy.However, Crocus uses a number of rules to set vertical discretization to provide a more realistic representation of snowpack.In Crocus, snow microstructure is defined by sphericity, dendricity, grain size, and a variable that shows the history of liquid water or faceted crystal in the layer.Implementation of metamorphism law in Crocus significantly improves the simulation of the evolution of snowpack density, permeability, and grain parameters.These detailed characterizations of snowpack properties enhance estimates of snowpack states and fluxes and make Crocus distinct from other snow models.We extracted the Crocus model from SURFEX8.1.1,available online at https:// opensource.umr-cnrm.fr/projects/surfex_git2and we have updated the code with the latest bug fix.The physics and validation of Crocus are detailed in Vionnet et al. (2012).The current version of the Crocus supports the Two-streAm Radiative TransfEr in Snow scheme (TARTES; Libois et al., 2013), the light-absorbing impurities scheme of Tuzet et al. (2017), and the multi-physics version of Crocus (Lafaysse et al., 2017).

Model Development
The primary objective of this study was to develop an interoperable methodology for integrating land surface processes that are currently absent in LSMs.Based on this strategy, we introduced a new methodology to couple a special model such as Crocus with LSMs. Figure 1 shows a schematic of LIS and how Crocus is integrated within LIS to enable interoperability at the model physics level across different LSMs.The model infrastructure in LIS already incorporates a suite of LSMs and hydrological models, with support for a large set of land surface parameters, meteorological boundary conditions, ensemble data assimilation, model calibration, and uncertainty estimation tools around these models.As described in Kumar et al. (2006), these capabilities are implemented in a highly modular fashion, with connectivity and data exchanges between them and facilitated through the features of the Earth System Modeling Framework (ESMF; www.earthsystemmodeling.org).ESMF provides generic functionalities for defining data to be exchanged between model components, by embedding information on the model grid and metadata within the ESMF data structures.This infrastructure has been expanded to allow the notion of specialized models to interact with the existing suite of LSMs in LIS.To show the functionality of the methodology and as a proof-of-concept study, we focused on snow processes that are lacking in LSMs.We hypothesized that the snow modules in the current LSMs are not sufficient to simulate the snow processes.Therefore, there is a need for a more complex snow module to simulate snow processes more accurately.Crocus is implemented as the first instance of demonstrating this flexible coupling between the LSM and a specialized physics model.Here we demonstrate the coupling of Crocus with two different LSMs: Noah-MP-4.0.1 model (Niu et al., 2011;Yang et al., 2011) and Noah-3.9(Ek et al., 2003).Comparison of model simulations with Noah and Noah-MP models integrated with Crocus helps in examining the added value of the Crocus integration as well.The use of these two models is also employed to demonstrate the flexibility of the LIS environment to enable interoperability at the model physics level.
At each timestep, the specialized model (Crocus) interacts with the given LSM through a set of import and export states (which are ESMF objects).In the coupled LSM-Crocus setup, Crocus receives inputs from a given LSM (ground temperature, soil liquid, and frozen water content) at each time step, and at the end of the time step, it returns the updated total SWE and total SD to the LSM.Therefore, even though the original snowpack schemes are running in the background, as part of the LSM physics, the LSM's total column snow states and layers are updated with the Crocus total snow state values.SD and SWE are essentially overwritten by the calculations from Crocus at each timestep.In this way, Crocus simulates the evolution of snowpack and provides detailed information about snow layers.The snow module in the LSM receives total SD and SWE from Crocus and propagates the changes in SD and SWE to the other states and fluxes through the LSM physics.We used the snow re-layering routines within the Noah-MP snow physics to convert total SD and SWE into a three-layer snow structure.Initially, we updated the total SD and the snowpack layer structure to match the SD values provided by Crocus.This is followed by the update of the snow layer ice content, liquid content, and temperature of each snow layer.For Noah, which has a single-layer snow scheme, we update the SD and SWE using total SWE and total SD from Crocus.Then other modules in Noah physics use these modified snow properties to update the other states and fluxes.Since the complexity of the physics between Crocus and the LSM varies significantly, it is impossible to transfer all the details of the snow physics from Crocus to LSM.The update approach used here provides a simple and effective approach to transfer the key outputs from the detailed snow physics calculations back to the LSM.This update scheme is similar to a sequential open water budget data assimilation system where a certain number of relevant state variables are updated in response to observational inputs.However, unlike the sporadic updates in a sequential data assimilation environment, the coupling between Crocus and the LSM provides spatially and temporally continuous updates at every timestep and every grid cell.The more frequent updates of the coupled system help to minimize the impact of water balance closure issues that are often present in sequential data assimilation environments.The coupling also relies on LSM physics to propagate the impact of the updated variables to other states and fluxes while ensuring water and energy balance closures.The loose coupling method introduced in this study currently does not account for the detailed interactions between snow and ground.LSMs compute the ground temperature using the updated SD and SWE and export that to Crocus.LSMs do not benefit from the computation of the heat transfer through the snowpack computed by the advanced physics of Crocus (two-way coupling).Therefore, the impact of the coupling method used in this study on the accuracy of snow microstructure needs further investigation.The snow cover fraction in Crocus within SURFEX approaches one quickly to avoid unrealistic temperatures at the bottom of the snowpack.We reinitialized the snow fraction to zero for snow-free ground.If the ground is not snow-free, we set the snow fraction to 1.We assume the difference in snow cover fraction converges (after a few timesteps) should not significantly impact our results in this study.Note also that the snow cover outputs of the coupled system are generated by the LSM, which recalculates snow cover based on the updated SD and SWE.This design allows for a flexible definition of the import and export fields based on the LSM being used.
The Crocus model itself is incorporated into the LIS framework through a number of interfaces that define the specifications of: (a) initial conditions, which describe the initial state of the land surface; (b) boundary conditions, which describe both the upper (atmospheric) fluxes or states also known as "forcings" and the lower (soil) fluxes or states; and (c) parameters, which are functions of soil, glacier fraction, topography, and other surface properties.A wrapper code is developed to compute the Crocus input parameters such as soil thermal conductivity, solar zenith angle, and air density, and invoke the Crocus physics from LIS.There are two interface subroutines that connect the wrapper code to the core LIS environment.These interfaces are used to provide the initial conditions, forcings, and other surface parameters to run the Crocus model within the LIS framework and to exchange information between the LSM and Crocus.Note that the Crocus physics code itself is incorporated within LIS without any modifications.This allows for the independent, continued maintenance and upgrades of both LIS and Crocus.
For this proof-of-concept study, we used a 50-layer snow configuration.We also did not activate and evaluate the blowing snow sublimation, the advanced parametrizations of the radiative transfer scheme (TARTES), and the light-absorbing impurities.Our implementation of Crocus within the LIS is 1-D, and there is no lateral exchange between grid cells, therefore, we cannot simulate wind-induced snow transport.Also, wind-induced snow transport is not essential on the continental scale simulations with a grid resolution of 12.5 km.Note that the current implementation of Crocus does not account for snow and canopy interaction.We did not modify the forcing data to compensate for the lack of canopy interaction.
The design for coupling specialized models such as Crocus within LIS provides several advantages.The specific set of interfaces designed for the specialized model integration isolates the implementation details.As a result, a model such as Crocus can be integrated without the consideration of all the other capabilities and features (data assimilation, uncertainty estimation, etc.) in LIS.On the other hand, once the model is implemented, it can readily be used with all the existing features in LIS with minimal additional development.For example, model integrations using Crocus can immediately employ all the support for meteorological and parameter data available within LIS.It can also exploit existing capabilities to use routing models, as demonstrated in the manuscript.Though the primary goal of integrating Crocus within LIS is to facilitate coupling with existing LSMs, Crocus can be run by itself as a stand-alone model.Finally, the integration of Crocus enables the use of existing capabilities in LIS such as data assimilation, though it is not the focus of this manuscript.

Experimental and Evaluation Design
We use the North American Land Data Assimilation Phase 2 (NLDAS-2) configuration to conduct all experiments in this study.The NLDAS-2 domain covers the Continental United States (CONUS) (25°-53°N, 125°-67°W ) at 1/8°spatial resolution.The NLDAS-2 near-surface meteorological data (Xia et al., 2012) is used to force the model, which includes gauge-based daily precipitation temporally disaggregated with radar data, bias-corrected shortwave radiation, and surface meteorological analysis.Additional details about the NLDAS-2 configuration can be found in Kumar et al. (2018).The simulations are performed with a 30-min time step for a 40-year time period from 1980 to 2020.The initial conditions are first generated by spinning up the Noah-MP and Noah LSMs from 1980 to 2020 and then reinitializing the Noah-MP and Noah LSMs again in 1980.The Hydrological Modeling and Analysis Platform (HyMAP; Getirana et al., 2012) river routing scheme within the LIS is used to generate routed streamflow estimates from the gridded runoff fields of the Noah-MP and Noah LSMs over CONUS.An offline Crocus simulation (Crocus uncoupled to LSM) was also performed to compare outputs from this simulation with those from coupled and uncoupled systems.
We used four different independent gridded and point scale measurements to evaluate the impact of using Crocus as a snow component of Noah-MP/Noah LSM on SD and SWE.The data sets include: (a) the daily 25-km gridded SD data from the Canadian Meteorological Centre (CMC; Brown & Brasnett, 2010); (b) daily 1-km gridded SD and SWE data from NOAA National Weather Service's National Operational Hydrologic Remote Sensing Center (NOHRSC) Snow Data Assimilation System (SNODAS; Barrett, 2003); (c) daily 4-km gridded SD and SWE data from University of Arizona (UASD, UASWE; Zeng et al., 2018); and (d) the daily Global Historical Climate Network (GHCN; Menne et al., 2012) point scale SD measurements.Following Kumar et al. (2014), a subset of the GHCN stations is used for evaluation in this work.To evaluate the impact of the detailed snow model Crocus on runoff, we use daily streamflow data between 2000 and 2020 from the U.S. Geological Survey over the 572 headwater watersheds with minimal impact from reservoir operations (Kumar et al., 2014(Kumar et al., , 2016(Kumar et al., , 2018)).The impact of coupling a detailed snow model Crocus with LSMs on simulated evapotranspiration (ET) is evaluated using Atmosphere-Land Exchange Inverse (ALEXI; Anderson et al., 2007), a daily gridded 4 km ET product derived from a thermal-infrared sensor (available from 2002 to 2016).The surface soil moisture (SM) and root zone soil moisture (RZSM) from the coupled and uncoupled simulations are compared against in situ SM measurements from the International Soil Moisture Network (ISMN; Dorigo et al., 2011Dorigo et al., , 2021)).The data set consists of quality-controlled in situ SM profiles from different networks.In our comparison, we used hourly surface and weighted (by layer thickness) average RZSM data from 1,251 stations over the CONUS.It is well known that SM is a wetness index and differs from one platform to another.Therefore, the anomaly correlation metric is used to evaluate the impact of the Crocus model coupling on SM.We also computed the normalized bias and normalized RMSE of SM and RZSM for coupled and uncoupled models (Figure 14). Figure 2 shows the locations of different evaluation regions used in this study.

Results
Here we provide a detailed evaluation of the impact of the coupled snow model Crocus with Noah-MP-4.0.1 and Noah-3.9LSMs on key land surface states and fluxes.The primary focus of the evaluation is on snow characteristics such as SD and SWE.In addition, we examine the impact of the Crocus coupling on runoff (Q) and streamflow, SM, and ET.The results from the coupled Noah-MP-Crocus and Noah-Crocus models are compared against modeled estimates from non-coupled Noah-MP and Noah LSM runs and several reference data products.The evaluation provides insights into the performance of the coupled system and helps to quantitatively characterize the improvement or degradation of states and fluxes over the study domain.The post-processing and development of the model evaluations presented here are enabled by LVT.To evaluate the temporal impact of the coupling of the Crocus snow model, the domain averaged mean seasonal cycles of SD, SWE, SM, RZSM, Q, and ET from coupled models (Noah-MP-Crocus and Noah-Crocus) and original models (Noah-MP and Noah) are computed and presented in Figure 4. Using the detailed snow model Crocus increases the Noah-MP (Noah) peak SD and SWE by 44% (50%) and 34% (70%), respectively.The largest impact of the change in snow states is observed on peak Q, which increases by 85% (90%).Smaller impacts on peak SM, RZSM, and ET are also observed, which increase by 7% (4%), 9% (3%), and 6% (<1%), respectively.The time-lagged impact from increases in snow accumulation on different surface variables is presented in Figure 4.The domain averaged peak SD and SWE occurs in February with the coupled models generating more snow.As the snowpack melts from warmer temperatures in the spring, increases in SM, RZSM, and Q are observed.The availability of larger snowpack in the models coupled with Crocus also results in the larger magnitudes of these variables in the springtime period.As ET becomes dominant in the summer months, the increases in SM and RZSM reduce gradually, but the wetter soils contribute to increased evapotranspiration.In the following sections, the quantitative evaluation of the coupled Noah-MP-Crocus and Noah-Crocus are presented through comparisons against Noah-MP and Noah with their native LSM snow physics and against several independent reference data sets.

Snow Depth and SWE
As noted previously, the coupling of Crocus with Noah-MP/Noah has a significant impact on the simulated SD and SWE. Figure 5 shows the spatial pattern of pixel-wise time-averaged SD differences between Noah-MP with its native snow model and four independent SD observations (left panels) and coupled Noah-MP-Crocus and the same set of observations (right panels).Both Noah-MP and coupled Noah-MP-Crocus show more snow accumulation than the CMC data (Figures 5a and 5b).Noah-MP significantly underestimates SD relative to SNODAS, UA, and GHCN stations (Figures 5c, 5e, and 5g), whereas the coupled Noah-MP-Crocus increases the SD  estimates relative to SNODAS, UA gridded products, and GHCN stations (Figures 5d, 5f, and 5h). Figure 5 also indicates that the coupled model overestimates the SD in higher elevation areas and underestimates in lower elevation regions in the SNODAS and UA comparisons.Note that CMC, SNODAS, and UA are also modeled products and therefore are subject to errors related to factors such as representativeness, spatial resolution, and modeling methods.Mass loss due to blowing snow sublimation affects the grid-averaged snowpack so that it can explain some of the overestimations of SD at higher elevations.We observed a similar pattern in Noah and Noah-Crocus comparisons (results not shown).
The spatially distributed comparison provides insights into the performance of the coupled model relative to observations and helps to quantitively characterize the difference between the model and observation over different geographical regions.Domain-averaged seasonal cycle time series over different regions of CONUS (Figure 2) shown in Figure 8 help interpret the regional nature of the seasonal cycles of SD relative to independent observations.The offline Crocus model produces significantly higher SD estimates than both the Noah and Noah-MP models.Nevertheless, the coupling framework has substantially increased the SD estimates, bringing them closer to those of the offline Crocus model.The coupled and uncoupled models overestimate SD relative to CMC, whereas the coupled Noah-MP-Crocus is in good agreement with SNODAS, and the coupled Noah-Crocus is in good agreement with the UA SD product (Figure 8a).Coupled Noah-Crocus is in more agreement with UA SD over the Northern Mountains and Plains, and Southeast (Figures 8c and 8g).Comparing SD from model integrations with SD from UA shows that Noah-MP is in more agreement with the UA values during the accumulation seasons, whereas during melt season, SD from the coupled Noah-MP-Crocus is in closer agreement to the UA SD.While both coupled systems use Crocus and the same forcing data, the differences between the modeled SD and SWE from coupled systems are partially due to differences in ground heat flux (as an input to the Crocus from the LSMs).For example, higher ground heat flux from Noah-MP relative to Noah during the melt season contributes to melting the snowpack faster, and SD in both coupled models converges toward each other at the end of melt seasons.-2017, 2003-2018, 2000-2018, and 2000-2016 for CMC, SNODAS, UA, and GHCN, respectively.The red colors represent a reduction in the RMSE when using the coupled model, and the blue colors represent an increase in the RMSE, relative to independent observations.
In situ measurements and gridded data products used in this study provide valuable information for the results evaluation.However, significant uncertainties are associated with these data products (e.g., Figures 6-8 and Tables 1 and 2) due to heterogeneous availability in space and time, differences in spatial and temporal resolutions, different measurement techniques, equipment error and malfunction, and human error.Limitations in the atmospheric forcing (especially in the precipitation amount and phase) are also a potential reason for some differences in the gridded products (CMC, SNODAS, and UA).Since uncertainty estimates of these data sets are usually unknown.Therefore, caution must be taken in the interpretation of results, especially when the improvements/degradations are marginal.
Figure 9 shows scatterplots of domain-averaged SD estimates from coupled and uncoupled Noah-MP and Noah relative to SD from the SNODAS and UA products.Overall, the Noah-MP estimates show a systematic negative bias relative to SNODAS, which is consistent with the SD map and domain-average seasonal cycle time series.Using the detailed snow model Crocus as a snow module of Noah-MP reduces the overall bias.The coupled Noah-MP-Crocus model moves SD estimates toward 1:1 line.However, the distribution of data points around the regression lines is similar to the uncoupled Noah-MP (Figures 9a and 9b).We also compared the scatterplots of model estimates relative to SD from the UA product.Overall, the improvement over the entire CONUS is smaller and the coupled model tends to overestimate SD (Figure 9b).SD estimates from both Noah and coupled Noah-Crocus underestimates values relative to the SNODAS and UA products.The coupled Noah-Crocus significantly improves SD estimates, reduces the negative bias, and moves the estimates toward the 1:1 line (Figures 9c and 9d). Figure 9d also shows that the data points are distributed more closely around the regression line.
Table 1 shows the root mean square error (RMSE) and bias of various model configurations relative to different independent observations.Over the entire CONUS domain, the coupled Noah-MP/Noah Crocus configurations provide SD estimates that have better agreement with SNODAS and UA data sets (see the bold-faced fonts).Overall, coupled Noah-Crocus configuration improves the RMSE (bias) of SD relative to SNODAS and UA SD by 2% (57%) and 12% (75%), respectively.Coupled Noah-MP-Crocus configuration improves the bias of SD relative to the SNODAS and UA SD by 83% and 60% and degrades the RMSE of SD relative to the SNODAS and UA SD by 10% and 16%, respectively.The RMSE and the bias of the offline Crocus relative to the SNODAS and UA SD estimates are comparable with those from the coupled systems.The coupled systems outperform the offline Crocus relative to the CMC and the GHCN.
Table 2 shows the RMSE and bias of SWE from uncoupled and coupled Noah-MP and Noah relative to the SNODAS and UA SWE products.Overall, coupled Noah-Crocus configuration improves the RMSE (bias) of SWE relative to the SNODAS and UA SWE by 2% (40%) and 13% (73%), respectively.Coupled Noah-MP-Crocus configuration improves the bias of SWE relative to the SNODAS and UA SWE by 42% and 89% and degrades the RMSE of SWE relative to the SNODAS and UA by 12% and 16%, respectively.The RMSE and the bias of the offline Crocus relative to the SNODAS and UA SWE estimates are comparable with those from the coupled systems.
An evaluation of the impact of coupling Crocus is examined for two extreme years in Figure 10.The years 2008 and 2012 are the wettest (highest snow accumulation) and driest (lowest snow accumulation) years over the simulation period and availability of SNODAS and UA data.).The analysis of the wettest and driest years over the simulation period indicates that the coupled models are capable of properly simulating the evolution of SD during extreme climate conditions (wettest and driest).Furthermore, the comparison between uncoupled and coupled models shows the positive impact of using the detailed snow model Crocus in our simulations.

Runoff and Streamflow
Streamflow and its temporal variations are important hydrological processes in snow-dominant water basins.Reliable estimates of timing and magnitude of river discharge are dependent on accurate simulations of winter snow storage and subsequent melt.The improvement in Nash-Sutcliff efficiency (NSE) and RMSE in streamflow computed as a normalized information contribution (NIC; Kumar et al., 2018) is shown in Figure 11.The NIC for NSE and RMSE are computed as a measure of improvement in streamflow using the following equations: where subscripts c and o represent the coupled LSM-Crocus and original LSM with their native snow model, respectively.Positive values of NIC indicate locations where the model coupled with Crocus provides more accurate streamflow, whereas negative NIC values are obtained in places where the model without Crocus performs better.
As shown in Figures 11a and 11b, the coupled Noah-MP-Crocus improves streamflow estimates in the Eastern United States, especially over the upper Mississippi, Ohio, and Mid-Atlantic river basins relative to United States Geological Survey (USGS) streamflow data (the map of river basins is available online at https://water.usgs.gov/wsc/map_index.html).Degradations in streamflow with the coupled Noah-MP-Crocus model are observed over the upper and lower Colorado, Missouri, Great Lakes, and South Atlantic river basins, which are likely due to the overall overestimation of snow accumulation over these basins.Comparatively, the use of Crocus with Noah generally improves streamflow estimates throughout the domain, as shown in Figures 11c and 11d.The Noah-Crocus coupling significantly improves streamflow skill in many river basins, including the upper Mississippi, lower Mississippi, Ohio, Tennessee, Great Lakes, New England, Mid Atlantic, South Atlantic, and Columbia relative to the USGS streamflow data.Over some of the river basins such as the upper and lower Colorado, Arkansas, Missouri, and the Great Basin, some degradations from Noah-Crocus are also observed.
Note that changes in streamflow skill are observed even in some of the non-snowy regions (e.g., southeast U.S.).This is an artifact of the rain-snow determination scheme in the model, which leads to triggering the Crocus model physics when a small amount of snowfall exists.In such areas, the absolute magnitude changes in the skill metric (NSE or RMSE) are indeed small.

Evapotranspiration
The impact of coupling the Crocus detailed snow model with Noah-MP and Noah LSMs on evapotranspiration is evaluated using ALEXI.Overall, the use of the coupled Noah-MP-Crocus slightly increases the ET over CONUS, whereas the use of Crocus with Noah has little overall impact.The improvement can be seen over the west coast as shown in Figure 12h.

Soil Moisture
The SM and RZSM from the coupled and uncoupled simulations are compared against in situ SM measurements from the ISMN.Mean monthly SM for model and in situ measurements are computed first for each grid cell.Then daily anomaly is computed by subtracting the monthly mean values from the daily SM value.The anomaly correlation is defined as the correlation between SM anomaly from modeling (coupled and uncoupled with the Crocus) and in situ measurements.Figure 13 shows the map of the change in anomaly correlation R for SM and RZSM for the simulation period of 2000-2021 for coupled Noah-MP-Crocus relative to Noah-MP (a and b) and Noah-Crocus relative to Noah (c and d).In this map, cool colors indicate locations of improvement, and warm colors indicate locations with degradation.Overall, the coupled models degrade both SM and RZSM.The coupled model slightly improves SM at higher elevations and degrades SM in lower elevations in these comparisons.Slightly larger improvements are obtained in RZSM, especially over the Cascades, and Rocky Mountains.

Discussion
Detailed parametrizations of snow physical processes used in Crocus are the main driver of the differences between the SD and SWE from the coupled and uncoupled systems.Noah model uses a simple single-layer snow scheme that lacks detailed representations of underlying snow processes.Though the Noah-MP snow model is more advanced with a 3-layer snow scheme, it still lacks the representation of processes such as snow grain evolution and metamorphism.The Crocus multilayer energy balance model, on the other hand, consists of detailed parametrization of snow processes as described in Section 2.3.The implementation of detailed snow physics from Crocus has enabled the coupled system to simulate processes such as permeability and grain parameters.It also provides more advanced physics for the simulation of snowpack density, melt and refreeze, and snow layer structure.These detailed characterizations of snowpack properties enhance estimates of snowpack SD and SWE, which in turn impact estimates of SM, RZSM, Q, and ET.
Figure 14  The bottom row of Figure 14 shows the normalized RMSE of each variable.The coupled Noah-Crocus improves the RMSE of SD, SWE, and Q and slightly degrades the RMSE of SM, RZSM, and ET.The coupled Noah-MP-Crocus slightly degrades the RMSE of all variables.While the bias reduction means the coupled model moves the estimates toward the reference data set, the marginal improvement or degradation of the RMSE of some variables indicates that the distance between the estimated and observed values remains high.Therefore, additional enhancements such as the integration of canopy and snow interaction and snow and ground energy exchange into the coupling system may be needed to translate the beneficial impacts on SD and SWE to other variables.
It is interesting to note that though coupling of Crocus improves the snow simulations in Noah (i.e., both bias and RMSE) and Noah-MP (i.e., bias), corresponding impacts on other water budget variables such as ET, streamflow, and SM are mixed.The climatology of these processes in the default model is a significant factor in whether model coupling to Crocus leads to additional improvements or degradations.For example, streamflow simulations are improved when Crocus is introduced with Noah, whereas with Noah-MP, degradations are observed relative to the default Noah-MP model.This is because the snow amounts updated in Noah by Crocus and timing of melt events led to improved agreement in terms of streamflow.Similarly, degradations in SM with the added use of Crocus are likely due to changes and shifts in the SM anomalies from additional snow reducing the level of agreements with in situ observations.Recent studies have also identified the limitations of the inherent coupling assumptions between surface states and fluxes as a key factor in enabling efficient information transfer within LSMs (Crow et al., 2018(Crow et al., , 2019)).Improvements in the realism of the native LSM processes through methods like calibration are needed to ensure consistency in translating improvements in snow simulations from the inclusion of a more detailed snow

Summary
This work presents a new framework that couples Crocus, a detailed snow model, with LSMs within the NASA LIS using the open water budget data assimilation concept.As a demonstration, Crocus is coupled with two LSMs: Noah-MP version 4.0.1 and Noah version 3.9.Through the definition of explicit interfaces that define the processes and states to exchange between the LSM and Crocus, model coupling is done in a flexible manner.In each timestep, Crocus runs the snow physics using the forcing data, snow-ground interface temperature, and soil water content from the LSM and returns the total profile SWE and SD states to the LSM.Then the LSM physics propagate these changes in SD and SWE to relevant water, energy, and carbon cycle fluxes and states within the LSM.This coupling framework demonstrates interoperability at the model physics level, which can be a paradigm for enabling more complex physics within existing land models, instead of hardwiring them within each and every LSM.
The utility of coupling the detailed snow model Crocus with Noah-MP/Noah is evaluated by comparing SD and SWE and other related states and fluxes with the original Noah-MP/Noah (i.e., Noah-MP/Noah with the native snow module) models.Additionally, the impact of coupling to Crocus is examined by comparing to several independent reference data products.The evaluation is performed over the CONUS using NLDAS-2 forcing and domain configuration for 2000 to 2021.The results show that the current medium-complexity snow scheme in Noah-MP/Noah is not sufficient to capture the snow physical processes and causes an overall underestimation of SD and SWE.The use of the more detailed snow scheme from Crocus that employs metamorphism laws to simulate the evolution of snow microstructure, improves the estimation of SD and SWE over the entire domain relative to the CMC, SNODAS, UA, and GHCN data sets.
The impacts of the detailed snow model Crocus on other water cycle variables such as streamflow are, however, mixed.The coupled Noah-Crocus leads to significant improvements in the RMSE and NSE skill metrics on major river basins and mixed results and an overall degradation in a few water basins.However, the coupling of Crocus with Noah-MP leads to an overall degradation of NIC metrics.The impact of using the Crocus snow model on ET was evaluated by comparing the modeled ET against an independent thermal remote sensing-based data set.The impacts of coupling efforts on ET are not statistically significant.The coupled Noah-Crocus generates almost the same ET as the uncoupled Noah and the coupled Noah-MP-Crocus slightly increases the ET.However, the peak ET occurs one month earlier in the coupled Noah-MP-Crocus than found with ALEXI.There is also an overall degradation in SM estimates with reduced anomaly R values with model simulations that employ the coupled Crocus.Some improvements in RZSM are observed in higher elevation areas such as the Cascade and Rocky Mountains.
Though our results showed an improvement in SD and SWE, a corresponding similar improvement in ET and SM is not observed.Uncertainties in the LSM parameterizations and the associated biases of these variables likely contribute to these discrepancies.As documented in recent studies, there are significant limitations of the inherent coupling assumptions used in the LSMs between prognostic states such as SM and snow with water budget estimates of ET and Q.This study points to the need for improving the representation of terrestrial coupling relationships within the LSM to guarantee information transfer from enhancements to processes and states, such as snow.Overall, the results from this coupling effort are promising and show the importance of using a detailed snow model in hydrological simulations.The coupling framework is modular and can be extended to enable the interaction between a given specialized physics model and the general LSM.For example, the framework can be extended to incorporate advanced physics formulations of other processes such as a groundwater scheme or canopy structure, while maintaining the linkages with the water, energy, and carbon cycle processes that are available within the native LSM.Note also that a byproduct of this coupling effort is the development of detailed snow profiles containing optical, thermal, and mechanical properties of snow layers.The availability of these detailed snow properties is particularly important for data assimilation systems that incorporate measurements from air and space-borne platforms.For example, the microstructure and thermal information of the snowpack is needed within radiative transfer models that translate between the raw measurements and geophysical variables.However, the evaluation of snow microstructure was left for future studies.We expect that the development of coupled modeling systems such as the one described here will enable more optimal utilization of snow remote sensing measurements to advance data assimilation efforts.

Figure 1 .
Figure 1.Schematic of the integration of Crocus within LIS as a "specialized" physics model.

Figure 2 .
Figure 2. A detailed map showing the location of different evaluation regions over CONUS; (a) entire domain, (b) central US, (c) Northern Mountains and Plains, (d) Midwest, (e) Northeast, (f) Southeast, (g) Southwest, and (h) west coast.The color bar shows the elevation (in m) of the study domain.

Figure 3 .
Figure 3. Changes in the time-averaged SD, SWE, SM, RZSM, Q, and ET from coupled Noah-MP-Crocus relative to the Noah-MP with its native snow model for 2000 to 2021.

Figure 5 .
Figure 5. Time-averaged differences of modeled SD (m) from Noah-MP with native snow model (left) and Coupled Noah-MP-Crocus (right) relative to four reference data sets.The periods in the mean SD comparisons are2000-2017, 2003-2018, 2000-2018, and 2000-2016  for CMC, SNODAS, UA, and GHCN, respectively.The red colors represent an overestimation of SD, and the blue colors represent an underestimation of SD relative to independent observations.
The coupled Noah-MP-Crocus model slightly overestimates SD relative to SNODAS in January, February, and March and slightly underestimates SD in April, May, and Jun.The coupled Noah-Crocus slightly underestimates SD compared to the UA SD products.Comparing SD from model simulations and SNODAS over different regions shows that the coupled Noah-MP-Crocus is in more agreement with SNODAS over Northern Mountains and Plains, Southeast, and Northeast (Figures 8c-8e, and 8f).

Figure 6 .
Figure 6.Changes in RMSE of SD (m) from coupled Noah-MP-Crocus relative to the original Noah-MP compared to (a) CMC, (b) SNODAS, (c) UA SD product, and (d) GHCN data as the reference.The time periods in the RMSE comparisons are 2000-2017, 2003-2018, 2000-2018, and 2000-2016  for CMC, SNODAS, UA, and GHCN, respectively.The red colors represent a reduction in the RMSE when using the coupled model, and the blue colors represent an increase in the RMSE, relative to independent observations.

Figure 7 .
Figure 7. Similar to Figure 6, but for Noah and coupled Noah-Crocus.

Figure 8 .
Figure 8. Domain-averaged seasonal cycle of SD for different regions over the CONUS.
Figure 10a shows the domain-averaged time series of SD for the water year (WY) of 2008 (1 October 2007, to 30 September 2008) for coupled Noah-MP-Crocus and coupled Noah-Crocus relative to the SNODAS and UA data over the CONUS.The coupled models perform well

Figure 9 .
Figure 9. Scatterplots of domain-averaged SD from model estimates (blue: Noah-MP/Noah; red: Noah-MP-Crocus/Noah-Crocus) versus SNODAS (a, c) and UA (b, d) SD products over the CONUS.The blue and red lines are linear regression lines corresponding to original LSMs and coupled LSMs, respectively.

Figure 10 .
Figure 10.Time series of domain-averaged SD for (a) WY 2008 (one of the wettest years of the simulation period) and (b) WY 2012 (one of the driest years of the simulation period).

Figure 11 .
Figure 11.Improvements in streamflow: (a) NSE for Noah-MP-Crocus, (b) RMSE (both shown as NIC using the USGS daily streamflow observations as the reference for Noah-MP-Crocus), (c) the same as (a) but for Noah-Crocus, and (d) the same as (b) but for Noah-Crocus.Blue colors indicate improvement from the coupled Crocus-LSM and red indicates degradation.

Figure 13 .
Figure 13.Differences in anomaly R values (during 2000-2020) for (left) SM and (right) RZSM from Noah-MP-Crocus relative to the Noah-MP with native snow model (top) and Noah-Crocus relative to the Noah with native snow model (bottom).The cool colors indicate improvements and warm colors represent degradations, from the use of coupling with Crocus.

Figure 14 .
Figure 14.Top row: Normalized mean bias for different states and fluxes (SD, SWE, SM, RZSM, Q, and ET).The columns in each panel represent the models Noah, Noah-Crocus, Noah-MP, Noah-MP-Crocus, in that order, and rows in each panel represent different regions over the CONUS (ALL).Bottom row: similar to the top row but for the Normalized RMSE.

Table 1
RMSE and Bias of the Modeled SD [mm] Relative to SD From Different Independent Data Sets Note.The bold font indicates improvement in the coupled model relative to the original model.

Table 2
RMSE and Bias of the Modeled SWE [mm] Relative to the SNODAS and UA SWE Product Note.The bold font indicates improvement in the coupled model relative to the original model.NAVARI ET AL.
All modeling efforts overestimate ET over the central US, Northern Mountains and Plains, and Southwest (Figures 12b, 12c, and 12g).All models underestimate the ET over Northeast, Southeast, and west coast (Figures 12e, 12f, and 12h).
shows the normalized bias and the normalized RMSE of different states and fluxes for uncoupled and coupled Noah-MP and Noah LSMs for several regions over the CONUS relative to UA products, ISMN SM data, USGS streamflow data, and ALEXI ET product.The top row shows the normalized bias, derived by dividing the original bias value by the maximum absolute bias value of each variable (i.e., SD, SWE, SM, RZSM, Q, and ET).The bottom row shows the normalized RMSE, derived by dividing the original RMSE by the maximum RMSE for each variable.The coupled Noah-Crocus reduces the bias for all variables (i.e., SD, SWE, SM, RZSM, Q, and ET) over the CONUS relative to the Noah estimates.While the coupled Noah-Crocus reduces the SD and SWE bias in all sub-regions (i.e., the central US, Northern Mountains and Plains, Midwest, Northeast, Southeast, Southwest, and west coast.), the results for other variables are mixed, with a general improvement.The coupled Noah-MP-Crocus reduces the SD, SWE, Q, and ET bias and degrades the SM and RZSM bias over the CONUS.