SENSITIVITY AND MODEL REDUCTION OF SIMULATED SNOW PROCESSES: CONTRASTING OBSERVATIONAL AND PARAMETER UNCERTAINTY TO IMPROVE PREDICTION by

Abstract The hydrology of high-elevation, mountainous regions is poorly represented in Earth Systems Models (ESMs), yet these ecosystems play an important role in the storage and land-atmosphere exchange of water. As much of the western United States’ water comes from water stored in the snowpack (snow water equivalent, SWE), model representation of these regions is important. This study assesses how uncertainty in both model parameters and forcing affect simulated snow processes through sensitivity analysis (active subspaces) on model inputs (meteorological forcing and model input parameters) for a widely used snow model. Observations from an AmeriFlux tower at the Niwot Ridge research site are used to force an integrated, single-column hydrologic model, ParFlow-CLM. This study finds that trees can mute the effects of snow albedo causing the evergreen needleleaf scenarios to be sensitive primarily to hydrologic forcing while bare ground simulations are more sensitive to the snow parameters. The bare ground scenarios are most sensitive overall. Both forcing and model input parameters are important for obtaining accurate hydrologic model results.

Approximately one-sixth of the world's population live in snowmelt-dominated regions and rely on snowmelt for their water supply (Barnett et al., 2005). Snow is a critical component of the hydrologic cycle in mountain headwaters systems, with melting snow providing seventy percent of the western United States' water supply (Chang et al., 1987). In fact, eighty-five percent of streamflow for the Colorado River comes from snowmelt in the Colorado headwaters region and the majority of Colorado's water supply is from snow fall and subsequent melt in the mountains (Ikeda et al. 2010, Barnett et al., 2005). Though snow is a vital piece of the hydrologic cycle, hydrologic models have little ability to accurately capture snowmelt in highly heterogeneous mountain systems largely because of lack of available data in mountainous regions (Bales et al., 2006). Most hydrologic models use datasets with coarse spatial resolution and the extrapolation of point data across a large region, which cannot capture the inherent heterogeneity of these systems (Lundquist, 2010;Dettinger, 2004). Yet, accurate characterization of this heterogeneity is invaluable for accurately modeling snowpack in these mountainous regions.
Given that the variability in snow parameters and forcing affect the modeled output of snow process, such as SWE, it is imperative that the model physics are correctly understood. One way to better understand snow model physics in mountain systems is through sensitivity analysis. Sensitivity analysis considers the uncertainty that is introduced during the process of collecting meteorological data or parameterizing models and extrapolating across a heterogeneous region. It is important to understand how these uncertainties affect model output in order to understand why the model is producing a certain result, along with determining the accuracy of that result. With regard to modeling snow water equivalent, some studies have evaluated the differences between models and their snow formulations, but few studies have conducted a sensitivity analysis for these snow formulations and parameters (Chen et al., 2014).
Studies that have conducted a sensitivity analysis on snow parameters in heterogeneous mountain regions find that the topographic characteristics and snow albedo parameterizations play a large role in determining the snow-covered area and SWE model outputs (Engel et al., 2017;Houle et al., 2016). This emphasizes the need for further sensitivity analysis into the snow parameters as well as the forcing parameters-which are driven, in part, by topographic heterogeneity-in order to better understand the model process linked to SWE calculation.

METHODS
This study conducts a sensitivity analysis on 900 model simulations with varying land cover type, meteorological and snow model parameters. This allows for a better understanding of how each of land cover types, forcing variables, and snow parameters contributes to the modeled output of snow water equivalent.

Integrated Hydrologic Model
Using an integrated and coupled model for this work allows for modeling the feedbacks of surface and sub-surface processes simultaneously, including groundwater, soil moisture, and atmospheric variables, all of which are important to snow accumulation and melt (Dai et al, 2003;Dickenson et al., 1993;Domine, 2011).
This study uses ParFlow-CLM (PF-CLM), a fully-integrated, coupled hydrologic model (Figure 2.1a). ParFlow solves the three-dimensional Richards equation for subsurface flow and the shallow water equations for the surface flow (Jones and Woodward, 2001;Kollet and Maxwell, 2006). CLM simulates the surface energy and water balances and passes this information to ParFlow, allowing for communication between all layers of the model (Maxwell and Miller, 2005). ParFlow-CLM has been used to analyze many different hydrologic problems relating to surface water-groundwater interactions, atmospheric interactions, residence times, and land cover-groundwater interactions (Condon et al., 2015;Gilbert et al., 2017;Bearup et al., 2016;Pribulick et al., 2016;Rihani et al., 2010).

ParFlow-CLM Model Physics
ParFlow-CLM contains a robust physical representation of snow processes, the focus of this study. CLM contains up to five snow layers and calculates snow degradation and accumulation based on albedo, snow age, snow metamorphism, and temperature (Dickenson et al., 1993). As precipitation falls, the canopy intercepts some of the precipitation until the default canopy capacity of 0.1mm is reached, at which point the precipitation falls to the ground. The fraction of precipitation that falls as direct throughfall is determined by: q "# = P * −e )*., * -./ 0 12 ./ 0 (2.1) where P is precipitation, L E AI is exposed leaf area index and S E AI is exposed stem area index. The fraction of precipitation that falls as rain varies linearly with air temperature as: where the rest of the precipitation falls as snow, with T frz being 273.16K. The depth of new snow is based on the density of fresh snow calculated from the Alta relationship (Anderson, 1976;Giddings and LaChapelle, 1961). When precipitation falls to the ground as snow, it is added to the already existing snow layers increasing snow depth and SWE. If no snow is present on the ground, a new snow layer is formed when the snow depth exceeds 0.01m.
Shortwave radiation is split between direct and diffuse shortwave radiation (SW t,dir , SW t,dif ) and the corresponding albedos are used to calculate the shortwave radiation absorbed by the surface by: Snow albedo is calculated separately for direct and diffuse solar radiation for either the visible or near-infrared waveband. Snow albedo also decays with age due to the accumulation of dust and dirt along with increasing grain size (Warren and Wiscombe, 1980;Dickenson et al., 1993). α X = α X,GYI * 1 − C X *  et al., 1993). Non-dimensional snow age (t snow ) is estimated by accounting for three processes that will alter snow albedo: 1) change in grain size with vapor diffusion as the snowpack ripens, approximating the vapor pressure using surface temperature (T g ) as: 2) increase in effective grain size near freezing approximated as: 3) the darkening of snow due to the deposition of soot or dirt (age3), which is assumed to be 0.3. The change in snow age over a time step of dt, (Dt snow ) combines these three age-effects: ∆τ FGHI = 10 )q * Δt * A N + A n + A t (2.8) which is then added to the current snow age to obtain the snow age for the next time step. Snow albedo is then further corrected for solar zenith angle, which describes the angle of the sun perpendicular to the surface for a given time, latitude, and longitude using equation: where cff is the snow albedo correction factor for a zenith angle greater than sixty degrees, coszen is the cosine solar zenith angle for the next time step, and sl is a factor that helps control albedo zenith dependence.

Study Site
The study site is the Niwot Ridge Long Term Ecological Research site, which lies in the Rocky Mountains near Nederland, CO, 65 km north-west of Denver (Figure 2.2).
Niwot Ridge houses an AmeriFlux eddy covariance (flux) tower that has been collecting data since 1998 and is located at 3050m elevation. This flux tower data is used to drive the model. The surrounding subalpine forest is about 97 years old and in a state of aggradation after recovering from early twentieth century logging (Monson et al., 2002).
The site was established in 1980 and is representative of alpine ecosystems in the southern Rocky Mountains (Blanken, 1998-). Due to the long-term and well-curated data record, this site is an ideal location for sensitivity analysis.

Model Design
While sensitivity analyses have been performed using large-scale, threedimensional models, it is computationally expensive because these models can contain millions of parameters (Condon et al., 2013). Therefore, this study uses a single-column approach to understand the physics and constraints of the snow processes of PF-CLM.
This approach decreases input parameters and computational expense. Modeling one water year at a single-column scale takes about twenty minutes on a laptop and still accounts for all the processes needed to run a sensitivity analysis. The domain represented here is a single-column with dimensions of 1 m 3 and five subsurface layers.
The site vegetation is characterized by evergreen needleleaf. This study uses a onedimensional meteorological forcing comprised of incoming shortwave and longwave radiation, air temperature, wind speed, atmospheric pressure, and specific humidity from the flux tower at Niwot Ridge for water year 2014. The precipitation data used in this study comes from the U.S. Climate Reference Network (USCRN), which is a network of climate monitoring stations across the United States (Diamond et al., 2013).
This data was obtained from the Boulder 14W gage, which is about one kilometer from the AmeriFlux site and has been collecting data since 2003. This precipitation is used rather than the AmeriFlux precipitation data due to a flux tower precipitation gage  This approach is valid if the active subspace is one-dimensional and the output is monotonic with respect to the input. This can be verified by the sufficient summary plot (SSP). The SSP (step 5) reveals whether the relationship between the quantity of interest and the combination of input parameters (the active variable) fits a univariate, scalar-value function fit, such as linear, quadratic, or cubic. Each point on the SSP corresponds to the inputs and output from one model realization. This takes the fourteen-dimensional input space, identifies the direction in which the fourteen parameters create a line, and plots this linear combination against the output, thus collapsing the fourteen-dimensional space to a one-dimensional space (Constantine, 2015).
The weights of the linear combination of the inputs reveal to which parameters the model is most sensitive. Input parameters with weights of larger magnitudes are more important because a change in these parameters causes a greater change in the output than other input parameters. A negative weight means an increase in the input parameter results in a decrease in the output and a positive weight means an increase in the input parameter results in an increase in the output (Jefferson et al., 2015).

Active Subspaces on Model Parameters
The meteorological forcing and the physical snow parameters embedded in the model control SWE. Shortwave radiation (sw) and longwave radiation (lw) in the direction toward Earth's surface, precipitation (precip), air temperature (airt), wind speed from west-to-east (windu) and from south-to-north (windv), atmospheric pressure (press), and specific humidity (hum) form the meteorological forcing. The value for snow decay due to accumulation of particles (age3), visible and near-infrared albedo of new snow (snal0 and snal1, respectively) (Eq. 5a, Eq. 5b), decay constant for visible and near-infrared albedo of snow (cons and conn, respectively) (Eq. 5a, Eq. 5b), and the factor that controls albedo zenith dependence (sl) (Eq. 9) comprise some of the important snow parameter values embedded in the model. By understanding how modifying these input parameters affects the quantities of interest of this study-output of peak SWE, monthly averaged SWE, timing of melt, and timing of peak SWE-the PF-CLM inputs are better constrained leading to more accurate model results that can be scaled to coarser resolution models. Additionally, running a sensitivity analysis on the snow and forcing parameters allows for a sensitivity comparison between values based on field data and values embedded in the model.
Using the accuracy measurements from the manuals for the instruments that collect the meteorological data and the literature surrounding snow parameters in hydrologic models (RM Young, 2000;Vaisla, 2016;Campbell Scientific, 2011;CS105;Hollinger and Richardson, 2005;Warren and Wiscombe, 1980;Warren, 1980;Grenfell et al., 1994;Dickenson et al., 1993;Bonan et al., 2002), we chose a range of values from which to sample each parameter for this sensitivity analysis (

RESULTS
After perturbing both the snow parameters and the forcing and computing the active subspace on 600 simulations, we observed the relationships between the input parameters and quantities of interest: monthly-averaged SWE, peak SWE, time of total melt, and time of peak SWE (Figure 3.1). We then ran an additional set of 300 simulations only changing the snow parameters and keeping the forcing consistent with the baseline and using bare ground land cover to better determine the sensitivity of the model to the snow parameters. The model sensitivities change depending on land cover and on the quantity of interest.

Effect of Land Cover on SWE
SWE is also controlled by land cover type (Figure 3. ground temperature model outputs (Figure 3.3). The ground temperature refers to the top of the snowpack when snow is present. This ground temperature remains frozen longer in the evergreen needleleaf baseline simulation than for the bare ground baseline simulation. The bare ground simulation also reaches greater ground temperatures than the needleleaf simulation, resulting in a warmer snowpack where snow will melt more quickly. Since SWE is controlled by land cover type as well as the input parameters mentioned above, we ran a set of 300 simulations using land cover type evergreen needleleaf and a set of 300 simulations using bare ground to better understand the relationships between SWE, land cover type, and PF-CLM input parameters.   The months most sensitive to changes in input parameters as a whole are April, May, and June, shown by the larger slopes on the sufficient summary plot. This means the linear combination of weights of parameters is greater in these months than the others. Snowmelt occurs during these months leading to an overall greater sensitivity when the weights are linearly combined than in the other months.        Figure 3.7 Active subspace results plot for peak SWE of the evergreen needleleaf simulations. Figure 3.8 shows the sensitivity of peak SWE to the parameters of the bare ground simulations. In these runs, PF-CLM is sensitive to both shortwave and longwave radiation as well as both the near-infrared and visible albedo coefficients (Table 3.2).

Monthly Averaged SWE
Peak SWE for these runs varies from 151.3mm to 471.5mm. The ranges of modeled peak SWE are greater for the bare ground simulations than the evergreen needleleaf simulations revealing a greater sensitivity of the bare ground model to input parameters.
The least-squares regression for the linear model in all simulations is 0.99, which means the variation in peak SWE can be explained by the linear combination of weights of the input parameters.

Active Variable Weight Peak SWE (mm)
Active Variable Weight  1mm. This is a smaller range than the bare ground runs that perturb all input parameters; however, it is not as small as the evergreen needleleaf runs. This indicates the difference in ranges of these simulations is due to the forcing parameters and land cover type.   Table 3.2 Active subspace parameters with a weight of magnitude 0.25 or greater for the peak SWE quantity of interest.

Time of Peak SWE
Time of peak SWE is the last quantity of interest studied. However, neither the evergreen needleleaf nor the bare ground simulations allow us to use active subspaces to study the sensitivity of the timing of peak SWE to the input parameters. Figures 3.13, 3.14, and 3.15 do not show a univariate, scalar-value function fit of the time of peak SWE to the active variable, which is a requirement for computing the active subspaces (Constantine, 2015).

DISCUSSION
The sensitivity of model parameters varies with both quantity of interest and land cover type. This study shows that trees mitigate the effects of sublimation more than they cause snow interception. One might expect the increased sublimation from bare ground to cancel with the increased snowfall from little vegetation interception (Biederman et al., 2014) leading to a bare ground snowpack similar to that with trees.
These results suggest, however, that sublimation causes a lower snowpack than lack of interception causes snow to accumulate. The effect of greater sublimation influence can be seen by the lower snowpack in the baseline bare ground simulation compared to the baseline evergreen needleleaf simulation.
The evergreen needleleaf land cover causes less shortwave radiation and wind to penetrate the snowpack. An increase in these two parameters causes an increase in sublimation (Biederman et al., 2014). However, since the effect of these processes are so minimal with evergreen needleleaf cover, a third parameter important to sublimation takes over-humidity. This is seen by an increased sensitivity to humidity in the months of snow accumulation (winter) and amount of peak SWE; an increase in humidity results in decreased sublimation and increased SWE. However, in the spring and summer months, the snowpack decreases, air temperatures increase, and humidity is no longer a driver of SWE since shortwave radiation is now large enough to take over the melting of snow. Humidity is not as important in the bare ground model because the absence of trees allows for radiation and wind to affect the snowpack.
The trees also dampen the effect of the snow parameters. CLM attributes a fractional vegetation cover of 0.8 to evergreen needleleaf land cover. This causes the model to become insensitive to snow albedo since such little shortwave radiation reaches the ground (Dickenson et al., 1993). However, the bare ground simulations are sensitive to the visible and near-infrared albedo coefficients as shortwave radiation can penetrate the snowpack.
The month with the lowest number of important parameters for both land cover types is September. This is due to September having the lowest average SWE since radiation and air temperatures are melting snow before it can accumulate in the evergreen needleleaf model. The bare ground model is insensitive to snow parameters in September due to lack of snow.
Furthermore, the overall sensitivity of the bare ground model is greater than the evergreen needleleaf model shown by larger ranges of all quantities of interest. The difference in ranges between the bare ground and bare ground with baseline forcing simulations can be attributed to the influence of the forcing parameters.

CONCLUSION
In this study, a single-column, integrated hydrologic model, ParFlow-CLM, is used to evaluate the sensitivity of four quantities of interest-monthly averaged SWE, peak SWE, time of total melt, and time of peak SWE-to the eight forcing parameters and six snow parameters using active subspaces. This allows for a global sensitivity analysis of fourteen parameters, which are then collapsed into a one-dimensional linear combination and compared to modeled output values (Constantine, 2015).
This study found that sublimation is the main driver of snow melt with humidity being the dominant input parameter, in the evergreen needleleaf model, throughout the winter months until radiation is strong enough to overcome humidity. The bare ground SWE is less than the evergreen needleleaf in the baseline runs and more sensitive overall to changes in the input parameters evidenced by greater ranges in the quantities of interest to change of the active variables than the evergreen needleleaf outputs. Due to no vegetation cover, the bare ground simulations are also sensitive to the visible and near-infrared albedo coefficients. Trees mitigate the effects of sublimation-by blocking wind and radiation-and also cause the model to be less sensitive to changes of the input parameters as a whole, with an especially decreased sensitivity to the snow parameters. This study reveals the importance of land cover type on the ParFlow-CLM snow model as well which physical processes the model deems important for evergreen needleleaf and bare ground land cover. This shows the importance of measuring radiation and humidity correctly in the field for all land cover types as well as creating correct albedo parameterizations for bare ground models.
This study focused on two land cover types, however, conducting a similar sensitivity analysis on all ParFlow-CLM land cover types would lead to a better understanding of the role of land cover on snow processes in hydrologic models.
Furthermore, conducting a sensitivity analysis on more pieces of the model (thermal processes, energy process, etc.) would lead to a deeper understanding of model physics and the sensitivity of all model parameters.
Through sensitivity analysis, the importance of the meteorological parameters, snow parameters, and land cover type to the important snow processes of hydrologic models can be analyzed. Mountain snowpack melt accounts for the majority of the Colorado River water supply, the western United States' water supply, and one sixth of the water supply for the world's population (Ikeda et al., 2010;Chang et al., 1987;Barnett et al., 2005). This highlights the importance of modeling these processes correctly, which requires an understanding how model input parameters correspond with modeled snow water equivalent.