Snowmelt control on spring hydrology declines as the vernal window lengthens

Snow disappearance date and runoff center of volume (R-COV) are modeled with the University of New Hampshire Water Balance Model, WBM [18, 19]. WBM is a process-based, modular, gridded hydrologic model that simulates spatially and temporally varying water volume; it is amongst the earliest developed [20] global hydrologic models (GHMs) and can be scaled to any study domain. WBM represents all major land surface components of the hydrological cycle, and tracks fluxes between the atmosphere, above-ground water storages (e.g. snowpack), soil, vegetation, groundwater, and runoff. A digitized river network connects grid cells, enabling simulation of flow through river and groundwater systems. Full documentation is provided in [18] and [19]; here we provide details on an update to the snow water equivalent simulation methods. Runoff as simulated by WBM has been used and validated extensively in studies of the northeastern US (e.g.[ 21], [22], and references therein) and in global studies that encompass other seasonally-snow covered forested regions from 40–60°N [18, 19, 23, 24]; additional runoff validation is given here in supplementary table S3.

As described in [18] and [19], WBM takes daily precipitation and temperature inputs and determines internally if the precipitation falls as snowfall or rain based on temperature thresholds of −1 °C for snowfall and 1 °C for snowmelt. Snow water equivalent is the balance between snowfall and snowmelt. For this study, a modification has been made to account for the impact of large orographic gradients on snowfall and snowmelt. The elevation distribution of each model grid cell is calculated from a 30 meter digital elevation model, resulting in binned elevation categories of vertical bands (bin size = 250 m elevation change). A temperature lapse rate of −6.4 °C km⁻¹ is applied to the mean daily temperature at the reference elevation for each binned elevation category, resulting in an adjusted mean temperature for the portion of each grid cell in each elevation category.

Here we define snow disappearance as the date by which <0.1 mm WBM-simulated snow water equivalent (SWE) remains, and no additional SWE is accumulated until the next winter season. We applied a mask to each simulated 30 year climatology (early 2010–2039, mid- 2040-2069, and late-century 2070–2099) ensemble mean (i.e. figure 2) that excluded grid cells where more than 50% of years in the specified time slice had fewer than 20 d when snow water equivalent was greater than 30 mm. Modeled SWE is validated against historical (1950–2013) SWE observations at 1034 stations (supplementary figure S1) by driving WBM with the gridded reanalysis climate product used to statistically downscale the LOCA daily temperature and precipitation projections [25]. While the ideal validation metric would be the day of year on which snow disappears, most snow observational records within our study domain are weekly to biweekly, obscuring the exact day on which snow cover was lost. Comparison of historically-simulated SWE disappearance date trends to the few observational SWE sites with daily observations show that, for the overlapping observational and simulation period of 1990–2005, both the model and the observations show no significant trend (see supplementary materials §3.1 for details).

Budburst, which indicates the start of the vernal window closing, is simulated using the thermal time model within the open source phenor modeling suite [26]. We use Phenocam data from [27] to parameterize the phenor thermal time model; this dataset represents a quality-controlled, peer-reviewed subset of all available Phenocam Network (https://phenocam.sr.unh.edu/) observations within the study domain. By comparing the Phenocam plant function type descriptions to the MODIS MCD12Q1-derived IGBP land cover categories [28], we find that these Phenocam observations are representative of the vegetation and land cover of the region, with some small exceptions. Of the 58 available phenocam sites from [27] within the study domain, 43 sites identify deciduous broadleaf as the primary plant functional type, and 16 of these sites have evergreen needleleaf listed as the secondary plant functional type (figure S2). 63.5% of the study domain is classified as deciduous broadleaf or mixed forest by the MODIS MCD12Q1 land cover product [28]; we consider these 43 sites to be representative of the study domain vegetation on average. We used the historical temperature data to which the LOCA climate data were downscaled [25] along with MODIS 3 d transition dates at the 43 sites to optimize the thermal time model of budburst date within phenor. We recognize that this does not represent the other land cover types (representing 10.3% of the domain); further research is needed on urban and developed land and cropland land cover types in this region. The results of model calibration are shown in figure S4.

Vernal window onset is defined as the date after November 1 on which modeled SWE is less than 0.1 mm and does not become positive again until the following snow season. Vernal window closure is the day of year when budburst occurs. Vernal window length is calculated as the number of days between snow disappearance and budburst. Grid cells in which there are fewer than 20 d with SWE >30 mm are considered to have no vernal window, as snow cover is insufficient to calculate a snow disappearance date. For the purpose of calculating regional averages consistently across all metrics (figure 2 and table 1), model grid cells with fewer than 20 d with SWE >30 mm are assigned a vernal window onset date of January 1 in order to estimate region-wide snow disappearance dates. In this way, snow-free grid cells are assigned the longest possible vernal window, rather than assigning a '0', which artificially shortens the regionally averaged vernal window length as snow-free grid cells are dropped from the regional average and it becomes more heavily weighted by colder grid cells. For grid cells that lack seasonal snow cover, vernal window duration is therefore the number of days between January 1 and budburst. Since the vernal window duration is a function of SWE disappearance date and budburst date, we rely on the model validations of those metrics as validation for simulating the vernal window.

To assess how a lengthening and disappearing vernal window impacts spring hydrology, we track the spring runoff center of volume (R-COV; defined as the 50th percentile of cumulative runoff from Nov. 1 through May 31), over time in our simulation results, following the methods of [1]. Similarly, we define the precipitation center of volume as the 50th percentile of cumulative precipitation from Nov. 1 through May 31, and the snowmelt center of volume as the 50th percentile of cumulative snowmelt from Nov. 1 through May 31. These precipitation and snowmelt metrics are used to evaluate the relative importance of snowmelt versus precipitation in controlling runoff, as described below in §3.2.

We used generalized least squares regression to both determine trends in vernal window transition dates and to compare differences in the rates of change in vernal window events. The trend analysis encompassed both the historical period and two climate forcing scenarios. Univariate statistical models evaluated trends over time in snow disappearance, R-COV, and budburst. For multiple regression models, dependent variables were day of year at which two vernal window events occurred (snow disappearance and budburst; snow disappearance and R-COV; or R-COV and budburst). Independent variables were year crossed with a covariate indicating vernal window event type (snow disappearance, budburst, or R-COV). We considered rates of change, i.e. modeled slopes, to be statistically significantly different if p < 0.05. Using the nlme package [29] in R 3.5.2 [30], we compared statistical models that included variety of autocorrelation and variance structures, using the multi-model inference statistic Akaike's Information Criteria to choose the model with the best overall fit [31]. Model fits that did not improve with the addition of autocorrelation or variance structure were assumed to meet the assumptions of linear regression [32].

Multiple linear regression of R-COV to identify the relative importance and explanatory power of snowmelt COV and rainfall COV was done using the R relaimpo package version 2.2–3 [33]. Confidence intervals were estimated using bootstrapping with 1000 iterations, and the relative importance was calculated using the Lindeman, Merenda, and Gold (LMG) method. Results are summarized in table S4.

We find that both the onset (snow disappearance; figures 2(A), (B)) and closing (budburst; figures 2(C), (D)) of the vernal window have shifted to an earlier day of year over the historical period (1980–2005) and will continue this earlier shift through the late 21st century (2070–2099). In addition, the date at which the vernal window opens changes more rapidly than the date at which it closes (table 1). Historically, the snow disappearance date has shifted at a significantly faster rate than budburst (p = 0.0003), advancing by −1.7 d per decade as compared to a shift of only −1.0 d per decade for budburst. In the future, snow disappearance continues to advance significantly more rapidly than budburst date (p < 0.0001 in differences between the rates of change between snow disappearance and budburst for both low and high forcing scenarios). By late century, the snow disappearance date shifts earlier by a total of −15 (low forcing) to −25 (high forcing) days, or an average of −1.6 (low emission) to −3.1 (high forcing) days per decade over the region where a snow cover was maintained through the entire 21st century. It is notable that some portions of the study domain will no longer have any snow cover through the winter by late century; here, we define the presence of snow cover in grid cell as at least 20 simulated days with SWE > 30 mm. Assigning a 'snow disappearance' date of January 1 to this portion of the domain that has no winter snow cover each year, the full regional average of snow disappearance advances by −26 (low forcing) to −42 (high forcing) days. Budburst date shifts earlier by only −9 (low forcing) to −14 (high forcing) days, which is −1.0 (low forcing) to −1.7 (high forcing) days per decade by late century (2070–2099).

Differential rates of change between the two events that mark the start and end of the vernal window results in an overall lengthening of this period of +15 to +28 d from historical to late century time periods (figures 2(E), (F)). We find that historically, from 1980 to 2005, the vernal window lengthened at a rate of +1.6 d per decade. While this historical trend continues at a similar rate under low forcing (+1.5 d per decade from 2010 to 2099), we find an increased rate of change (+3.3 d per decade from 2010 to 2099) in the high forcing scenario. Notably, the vernal window length begins to stabilize at the end of the century in the low forcing scenario, gaining only +0.92 d per decade from 2070–2099. This is in contrast to the high forcing scenario, in which the vernal window gains +4.0 d per decade in the late part of the century.

Since the vernal window onset is defined by snow disappearance, areas that lose winter snow cover due to climate change also lose their vernal window (figures 2(E), (F), grey areas, and figure 3(A)). Historically, 27% (±11%) of the study area had <20 d of simulated snow water equivalent <30 mm from Nov. 1 through May 31; this metric is used here to define seasonally snow-free areas. The snow-free and therefore vernal window-free region of the study area increases to 43% (±14%) under low forcing, and 59% (±16%) under high forcing scenarios (figure 3(A)). This loss of snow cover impacts not only the existence of a defined vernal window, but also notably alters the characteristics of spring hydrology as snowmelt typically drives the spring freshet, a pulse of freshwater that rapidly enters river systems near the opening of the vernal window. To assess how a lengthening and disappearing vernal window impacts spring hydrology, we track the spring runoff center of volume (R-COV. There is no significant trend in the LOCA-projected precipitation center of volume (COV), allowing for a controlled assessment of the impact of the vernal window; all changes in simulated R-COV here are caused by changes in temperature-driven snow dynamics, liquid versus solid precipitation, and evapotranspiration.

Zoom In Zoom Out Reset image size

Like the start and end of the vernal window, we find that the spring R-COV date shifts earlier, with a rate of −1.2 (low forcing) to −2.0 (high forcing) days per decade over the 21st century (table 1). The vernal window-free portion of the study region has no change in R-COV, which is consistent with expectations based on a constant precipitation COV. An unexpected finding, however, is that despite retaining snow cover and a vernal window, spring hydrology in the most northerly and high-elevation parts of our study domain are no longer snowmelt-dominated by the late 21st century in the high forcing scenario (figure 3(C)); rather, the R-COV timing shifts closer to that of the vernal window-free region. Stated another way, even the portion of the domain that retains a snow cover has a change in spring hydrologic character toward a rainfall-dominated regime. Multiple linear regression analysis of region-wide R-COV as a function of snowmelt COV and rainfall COV (with rainfall defined as liquid precipitation) shows that, historically, snowmelt COV was relatively more important to spring freshet timing than rainfall COV (a scaled relative importance of 0.56 compared to 0.44). By late century, the relative importance of these two explanatory variables switches, with rainfall COV importance greater than snowmelt COV under both high and low forcing scenarios (figure 4(A)). Where a vernal window remains present, snowmelt COV retains its position as the more important of the two drivers, yet its relative importance is reduced from 0.68 historically to 0.60 (low forcing) and 0.56 (high forcing), while rainfall COV increases from 0.32 historically to 0.40 (low forcing) to 0.44 (high forcing) (figure 4(B)). Notably, the explanatory power of the combination of these two variables declines from 77% historically to 53% by late century (high forcing) in the vernal window present region, suggesting that additional mechanisms increase in their control over spring hydrology.

Zoom In Zoom Out Reset image size