https://doi.org/10.5061/dryad.2z34tmpvz

Description of the data and file structure

This research documents changes in riparian ecosystems in the Colorado River headwaters in Rocky Mountain National Park. There are time series data that extend 10-70 years back in time. There are photo time series and data time series. Some files are in GIS format.

Files and variables

File: allwell_dtw.csv

Description: Depth to water table data 

Variables

File: beaverwell_dtw.csv

Description: Depth to water table for wells near historic beaver dams

Variables

File: BeaverPonds1953to2019.zip

Description: GIS files of extent of beaver ponds in the study area from 1953-2019.

File: climate_wb_flow.csv

Description: Water balance data.

Variables

Missing data code: NA

File: pondarea.csv

Description: Area of open water in the study area through time.

Variables

File: willowarea.csv

Description:   Willow area data for time series in the study area.

Variables

File: Vegetation.csv

Description: Vegetation composition data.

Variables

File: willowstems.csv

Description: Time series data on willow height inside and outside exclosures in the study area.

Variables

File: VegMetricsType.csv

Description: Vegetation data analyzed in this paper.

Variables

File: WillowPatchHeight.zip

Description: GIS data of tall and short willow patches in the study area through time.

Code/software

Appendix S1: Study area

The Kawuneeche Valley is bisected by the Colorado River channel that ranges from 5 to 10 m wide with a riparian zone 0.7 to 1.2 m above the river (Woods 2001), and up to several hundred meters wide (Figure 2 in main text). The valley is bordered on the east by the Front Range consisting of Precambrian metamorphic rocks and on the west by the Never Summer Range of upper Oligocene granitic magmas. Valley slopes are covered by Pleistocene era lateral moraines (Braddock and Cole 1990), with alluvial fans along hillslope margins. Organic soils 0.4 to >1.5 m thick are present in fens along valley margins in many areas supported by hillslope ground water, and silt loam. Loamy sand texture mineral soils elsewhere in the valley average 0.9 m thick. The valley is underlain by 3 to 4 m of gravel alluvium, with 15–122 m of Holocene and upper Pleistocene alluvium below that (Braddock and Cole, 1990).

Appendix S2: Methodological and analysis details

All analyses in this paper were conducted using R version 4.3.3 (R Core Team 2024) with libraries and functions as detailed below. All Geographic Information System (GIS) based work was done using ArcGIS Pro version 3.1 (ESRI 2023). Data and select R code are housed on Dryad (DOI 10.5061/dryad.2z34tmpvz).

Water balance

We estimated a daily water balance (WB; Lutz et al. 2010, Thoma et al. 2020) at a central location near vegetation and ground water sites (see Figure 2, main narrative for this location). Our WB model partitions precipitation (from DAYMET version 4; Thorton et al. 2021) into rain or snow following Jennings et al. (2018), the latter of which accumulates until temperatures (also from DAYMET) become warm enough to melt the snow (Hock 2003). We estimate a direct runoff parameter following NRCS (2017). We calculate potential and actual evapotranspiration following Oudin et al. (2005) with slope and aspect estimated from a digital elevation model to account for a localized heat load. We use SSURGO soil data for an initial soil water holding capacity (NRCS 2021). Subsequent soil water estimates account for actual evapotranspiration and an estimate of shade from plot level vegetation data. Climatic water deficit is estimated as the amount of additional water vegetation would use if available, calculated as the difference between actual and potential evapotranspiration (Stephenson, 1998). The model calculates water in excess of soil storage capacity as runoff. Runoff values are used to accumulate a daily cumulative runoff or “drainage” term following Croke et al. (2005) We assume a unit hydrograph to derive coefficients for “quick” or event-based runoff and “slow” or base runoff. These coefficients are optimized with calibration to DTW by iterating parameters singly and as pairs across a range of values, maximizing adjusted r2 with DTW.

General additive models of climate, water balance, hydrology and willow stem height

We restrict most data in this paper to a May – September growing season under the assumption that this is the portion of the water year that is likely most important to vegetation and ungulate use of vegetation. A growing season-based time step also reduces variance within and across years. We summarize most variables as medians across each growing season. Medians are a preferred summary statistic for data with non-normal distributions, outliers and high variance. (i.e., Wilcox et al 2018). Most of our climate and hydrology data are of this form (Steinschneider et al. 2015). For variables where the response of interest is a sum (i.e., precipitation) we instead use a total across each growing season. We plot monthly (or weekly for ground water) point values in Figure 3 and 5 in the main narrative but expect readers to focus on the fitted lines in these Figures which reflect the overall period of record trends (slope, intercept, and point wise confidence intervals, see below).

We use Generalized Additive Models (GAM; Wood 2011, 2017) for most analyses. GAMs are well suited to data that are inherently nonnormal and that can have nonlinear relationships with their predictors. A GAM allows non-linear or “smooth” functions to be applied to all or a subset of predictor(s) as patterns in data indicate. Smooths are regularized nonparametric functions that can fit simple to very complex non-linear relationships. GAM models may also include predictors with a linear relationship with the response and/or with an a priori or defined non-linear relationships, such as a log. GAMs are additive, as the name suggests, allowing interpretation of each variable in the model independent of the others. Each smooth or linear term represents the unique effect of a single predictor on the response variable, controlling for the other predictors. The effect of each variable is assessed without interactions, meaning that the model assumes the effect of one predictor does not change based on the value of another predictor. All GAM models were estimated using the mgcv package in R (Wood 2023).

Models were developed using a model selection process. We first estimated an intercept only null model. Variables were added and retained with four selection criteria: 1) change in Akaike Information Criterion (AIC; Venables and Ripley 2002) with a delta AIC threshold value of 2.0; 2) overall model and individual predictor interpretability; 3) adjusted r2 values; and 4) a candidate predictor’s p value. In general, we use a structural approach following Grace and Irvine (2020) and initially developed and included candidate predictors as those with likely direct or readily explainable indirect connections to a response.

GAMs with sufficient sample size within a growing season used a mixed model with years treated as a random factorial term and all other predictors fixed. All time-based terms in models were centered and normalized. Where exploratory work suggested a linear trend was present and also improved a GAMs quality, we allowed a linear or parametric POR trend estimate. In models that included a month within each growing season effect, these relationships were clearly non-linear and estimated as a smooth. The total number of observations for climate models was 365 (5 months across 73 years from 1950 to 2023) while water balance responses was 215 (data spanned 1980 to 2023) and for flow responses that began in 1953, the total number of observations was 350. Estimates of the period of record trend had effective sample sizes set by the number of years in the model (i.e., 73).

The GAM of DTW included sites as a random term and was improved by additional predictors for Colorado River flow, one year lagged snow water equivalent (SWE; both best fit as a smooth) and a linear term for climatic deficit. The lagged SWE term likely captures lags of water availability in a simpler way than if we had over specified the model with an explicit month term. Finally, the DTW model also included a test of beaver influence as a factorial term. We felt PDO is an at least decadal phenomenon, and we determined its finest scale time step of monthly was too coarse to apply in the DTW model. The total number of observations for DTW models was 942 with the nonlinear period of record trend based on 21 years.

We model the effect of months within each year’s growing season as a nonlinear or smoothed term in each climate, water balance and hydrologic GAM (note that SWE has only one value per year and so we cannot do this with this response). We restrict the complexity of this smooth to avoid overfitting. While a smooth on months does not explicitly model a month effect of test for independence between months, it provides a flexible non-linear fit to seasonal (monthly) trends. This helps model shared seasonal structure by accounting for patterns across months, reducing unexplained variability, and, in effect, isolating the estimate of an overall trend across the period of record (the primary goal of these models) from the within year variability. Likewise, the random effect used for years across the period of record term helps model variability across years non-linearly by treating years as a random grouping factor. This approach has several benefits. It captures unmeasured factors influencing the response that vary across years (e.g., climate anomalies or environmental conditions) and by modeling year-to-year variability, the term accounts for correlations induced by shared year-level factors, which can help approximate independence across years. We present partial residual plots in this Supplemental’s Results section to help visualize individual model terms and their fits.

Models for total and mean pond size in each year could only use a single predictor for years (as a factor). Effective sample size for total and mean pond size was seven as summarized from 2401 patches into total areas and mean pond size.

We report all p-values and interpret as statistically significant any result as meaningful if its p value is below or equal to 0.05. However, test results that are between 0.05 and ~0.10 are still likely meaningful and we interpret these marginal cases accordingly.

We visually examine autocorrelation and partial autocorrelations via the R autocorrelation functions ACF() and PACF() and then test for an effect using a Ljung and Box test (Ljung and Box 1978) with each lag set based on the model's sample size. If this test indicates that residuals were autocorrelated we further examined the ACF and partial ACF plots. We do not expect model results are strongly influenced even with a significant Ljung and Box test if there was not a strong suggestion of autocorrelation in the residuals as shown in these plots. We include partial plots in this Supplemental’s Results section when there is clear evidence of a likely meaningful autocorrelation.

Confidence intervals (95%) were estimated around each linear trend’s coefficient. For non-linear trends, confidence intervals were pointwise estimates to allow the intervals to follow the variability of a smooth.

Nonparametric models of total pond and willow patch area

Differences in median tall and short willow patch size between 1999 and 2019 were assessed with non-parametric Kruskal Wallis tests from the base stats package in R as these variables were strongly non-normal. Sample size for willow patch models ranged from 199 to 454 depending on height class and inside/outside the park. To test for differences in total area across only two periods, we used a permutation test to generate a null distribution of differences in sums expected under the assumption of no true difference between the periods. We then calculated an exact probability as the sum of all differences in sums greater than the observed difference divided by the number of permutations.

Vegetation metrics

Total absolute non-native invasive taxa cover

Nativity describes whether a species is found within its area of evolutionary origin and/or arrived without human intervention. Nonnative species are often introduced through intentional or unintentional human action (Pysek et al. 2004, Fertig 2011). Nonnative species can be invasive and can have undesirable effects on ecosystem function (Byers et al. 2002, Levine et al. 2003, Fridley et al. 2007). They have been linked to reduced overall species diversity (Meiners et al. 2001), altered resource dynamics (Ehrenfeld 2003), and shifted interactions between species (Christian and Wilson 1999). Nonnative species are often the focus of park management of vegetation, including at RMNP. We created an index of the degree of invasion by weighting the average of the relative cover of taxa by coefficients of invasiveness (“I-ranks”) as developed by Morse et al. (2004). I-ranks are scored based on the ability of a species to change ecosystem processes; invade relatively undisturbed ecological communities; disperse to new areas readily; and cause substantial impacts on rare or vulnerable species or ecological communities, or high-quality examples of more common communities.

Conservatism

Conservatism describes a species’ fidelity to a specific habitat or range of environmental conditions absent of human disturbance (Wilhelm and Ladd 1988, Herman et al. 1997, Matthews et al. 2015). Anthropogenic impacts can cause dramatic shifts in ecological processes and habitat conditions and push disturbance regimes outside a natural range of intensity, frequency, and duration. More conservative species are not able to quickly respond to such rapid alterations compared to broad-niche generalists and are often the first to disappear from habitats heavily impacted by human activities. The composition of conservative species at a particular site integrates spatial and temporal impacts and can serve as an indicator of ecological integrity or condition. To assess conservatism, we use metrics based on “coefficients of conservatism” (also known as “C-scores”) assigned to the flora of Colorado by a panel of experts following the methods described by Swink and Wilhelm (1994). Colorado C-scores were originally assigned in 2007 (Rocchio 2007b) and updated in 2020 (Smith et al. 2020). C-scores range from 0 to 10 and represent the estimated probability that a plant shows high fidelity to landscapes relatively unaltered from pre-European settlement conditions. Nonnative species are given C-scores of 0 by default. Low C-values are assigned to species that demonstrate little fidelity to unaltered landscapes or have wide ecological tolerances and may be found almost anywhere. High C-values are assigned to species only found in high quality natural areas and that cannot tolerate habitat degradation. Using the C-scores of all or select species present within a site, a suite of metrics can be calculated that convey different aspects of the site’s condition and disturbance history. We estimated the conservatism of each sample by averaging cover weighted conservatism scores for all taxa.

Wetland affinity

Wetland affinity measures the prevalence of species in a community that have a demonstrated ability because of morphological or physiological adaptations and/or reproductive strategies to achieve maturity and reproduce in wetlands. Higher wetland affinity in a sample indicates that vegetation composition includes a higher proportion of wetland obligates and thus likely higher and/or more stable water tables. Because the mix of species at a site responds over time and space to variation in hydrologic regime, wetland affinity integrates seasonal and annual fluctuations in groundwater levels (Grace et al. 2012). Wetland affinity is likely a useful, integrated proxy for more complex and expensive measures of groundwater hydrologic regime (Loheide and Gorelick 2007). For wetland affinity, we estimated a mean cover weighted score across all taxa in a sample.

Composition and structure

Multivariate vegetation composition analyses were conducted using vegan 2.6-4 (Oskanen et al. 2022) and the R vegan library (Oksanen et al. 2024). Wetland community types were first identified with hierarchical agglomerative cluster analysis using Sorenson distance and flexible beta linkages less than or equal to -0.25. Indicator species analysis (Dufrêne and Legendre 1997) identified species characterizing each cluster and a Monte Carlo analysis was applied to understand the significance of species’ indicator values within each cluster. Indicator species analysis was used to prune the dendrogram and optimize the number of clusters (McCune and Mefford 2018). We averaged p-values across all species for each cluster level using Monte Carlo analysis and chose the cluster with the lowest average P value. Final site type assignment in each year and indicator species were reviewed and adjusted based on our long-term experience in the system. Non-metric multidimensional scaling (NMDS) was used to determine community dissimilarity across years and wetland community types. Analysis of similarities (ANOSIM) were used to test the significance of the difference in species composition among community types across and within years. NMDS paired with vector fitting (function ‘envfit’ in vegan) was used to examine the relationships between overall community composition and total absolute non-native invasive taxa cover, mean conservatism score, and mean wetland affinity score. To visualize and describe community composition changes across years and wetland types, a single NMDS was used. To determine which summary vegetation metrics (total absolute non-native invasive taxa cover, mean conservatism score, and mean wetland affinity score were predictive of community composition in each sample year, two NMDS with two corresponding PERMANOVAs were run, each with 1000 permutations based on Euclidean distances. Total number of observations for the ordination was 784 (number of vegetation records across both 1998 and 2021).