<jats:title>Abstract</jats:title><jats:p>Correctly modeling snow is critical for climate models and for hydrologic applications. Snowpack simulated by six land surface models (LSM: Noah, Variable Infiltration Capacity, snow‐atmosphere‐soil transfer, Land Ecosystem‐Atmosphere Feedback, Noah with Multiparameterization, and Community Land Model) were evaluated against 1 year snow water equivalent (SWE) data at 112 Snow Telemetry (SNOTEL) sites in the Colorado River Headwaters region and 4 year flux tower data at two AmeriFlux sites. All models captured the main characteristics of the seasonal SWE evolution fairly well at 112 SNOTEL sites. No single model performed the best to capture the combined features of the peak SWE, the timing of peak SWE, and the length of snow season. Evaluating only simulated SWE is deceiving and does not reveal critical deficiencies in models, because the models could produce similar SWE for starkly different reasons. Sensitivity experiments revealed that the models responded differently to variations of forest coverage. The treatment of snow albedo and its cascading effects on surface energy deficit, surface temperature, stability correction, and turbulent fluxes was a major intermodel discrepancy. Six LSMs substantially overestimated (underestimated) radiative flux (heat flux), a crucial deficiency in representing winter land‐atmosphere feedback in coupled weather and climate models. Results showed significant intermodel differences in snowmelt efficiency and sublimation efficiency, and models with high rate of snow accumulation and melt were able to reproduce the observed seasonal evolution of SWE. This study highlights that the parameterization of cascading effects of snow albedo and below‐canopy turbulence and radiation transfer is critical not only for SWE simulation but also for correctly capturing the winter land‐atmosphere interactions.</jats:p>