A Physics-Constrained Neural Differential Equation Framework for Data-Driven Snowpack Simulation

TL;DR

This work addresses the challenge of predicting seasonal snowpack evolution in climate models by learning a physically constrained rate equation for snow depth from hydrometeorological forcings using a neural differential equation framework. A two-component architecture combines a trainable predictor with hard-threshold constraint layers to enforce conservation-like bounds on $dz/dt$, enabling stable integration at multiple time steps without retraining. On 44 SNOTEL sites for training and additional out-of-sample sites, the model achieves a median daily snow-depth error of $8.8\%$ with $NSE > 0.94$ when SWE data are included, and generalizes to unseen sites where traditional parameterizations struggle, while SWE-prediction errors remain near $12\%$. The method shows robust density predictions and can operate at finer temporal resolutions, offering a scalable, physics-consistent alternative for global and regional climate applications and potential applicability to other constrained dynamical systems.

Abstract

This paper presents a physics-constrained neural differential equation framework for parameterization, and employs it to model the time evolution of seasonal snow depth given hydrometeorological forcings. When trained on data from multiple SNOTEL sites, the parameterization predicts daily snow depth with under 9% median error and Nash Sutcliffe Efficiencies over 0.94 across a wide variety of snow climates. The parameterization also generalizes to new sites not seen during training, which is not often true for calibrated snow models. Requiring the parameterization to predict snow water equivalent in addition to snow depth only increases error to ~12%. The structure of the approach guarantees the satisfaction of physical constraints, enables these constraints during model training, and allows modeling at different temporal resolutions without additional retraining of the parameterization. These benefits hold potential in climate modeling, and could extend to other dynamical systems with physical constraints.

Figures (13)

• Figure 1: Structure of the model's predictive component. Blue lines indicate a trainable linear transformation of the input (of $k$ scalar variables), including a bias. Colors indicate the activation function used upon collection at the node (ReLU = Rectified Linear Unit, ELU = Exponential Linear Unit, Id = Identity), as noted in the legend. The hyperparameter $n$ sets the width of the internal mixing layer.
• Figure 2: Architecture of $M$, highlighting the constraint component attached to the predictive component (the grey pentagon) described in Fig \ref{['f:pred']}. The chosen structure enforces increasing snow depth only under precipitation and non-negativity of snowpack height, and is equivalent to a max/min block on the output. Weight colors indicate the constant's sign and activation functions follow the color scheme given in Fig. \ref{['f:pred']}.
• Figure 3: Distribution of SNOTEL sites used for training the network. (a) Training sites as visualized over the United States. (b) Training sites visualized with elevation vs their average nonzero snowpack height, $\Bar{z}_+$. The color bar is a visual indicator of $\Bar{z}_+$ for visualization on the spatial map.
• Figure 4: Performance of $M$ and $\tilde{M}$ against Snow17 with (SN17O) and without (SN17) observational $\mathrm{SWE}$ data for generating $z$ timeseries over the 44 validation sites. RMSE indicates root mean square error, NSE the Nash-Sutcliffe Efficiency, and SPE the average L1 error normalized to the average nonzero depth, to measure percent error. Boxes outline the 25% to 75% quantiles, with the bar at the median, while whiskers mark the extremes and dots indicate outliers, which lie beyond 1.5 times the interquartile range (box width) from the box. Vertical axes limits are chosen to show all data points.
• Figure 5: Performance of the neural parameterization against Snow17 over the 14 testing sites in depth timeseries generation. Labeling convention follows that in Fig. \ref{['f:validbox']}. One outlier for SN17 with an NSE of -2.8 is not shown on the plot to aid in the scaling of the other values.