Statistical framework to assess long-term spatio-temporal climate changes: East River mountainous watershed case study

Evaluation of long-term temporal and spatial climatic change in mountainous regions is a critical challenge because of the interactive effects of multiple land and climatic factors and processes. Here we present the application of the statistical framework to the assessment of changes of climatic conditions, using data from 17 meteorological stations across the East River watershed near Crested Butte, Colorado, USA, and spanning the period from 1966 to 2021. The framework is developed based on (1) a time-series analysis of daily, monthly, and yearly averaged meteorological parameters (temperature, relative humidity, precipitation, wind speed, etc.), (2) evaluation and time series analysis of potential evapotranspiration (ETo), actual evapotranspiration (ET), aridity index (AI), standard precipitation index (SPI) and standard precipitation-evapotranspiration index (SPEI), and (3) a temporal-spatial climatic zonation of the studied area based on the hierarchical clustering and PCA analysis of the SPEI, because the SPEI can be considered an integrative characteristic of the changes of climatic conditions. The Budyko model, with the application of the Penman–Monteith equation for the estimation of ETo, was used to determine the ET. The time series analysis of the AI is used to identify the periods with energy limited and water limited conditions. Hierarchical clustering of site locations for the three temporal segments of the SPEI showed a significant temporal-spatial shifts, indicating that dynamic climatic processes drive zonation patterns. Therefore, the watershed climatic zonation requires periodic re-evaluation based on the structural time series analysis of meteorological and water balance data.


Introduction
Evaluation of long-term temporal and spatial features of climatic change in mountainous regions is particularly challenging due to strong gradients in energy budgets, topographic complexity and heterogeneous soils and vegetation. This limits our ability to predict the consequences of global change for the ecology and ecosystem services of mountainous regions. Various meteorological (such as temperature and precipitation) and hydrological variables (such as evapotranspiration), as well as various drought indices, which are derived from meteorological and hydrological variables, are amongst the metrics used to characterize climatic changes. For example, to assess spatiotemporal drought variability, researchers commonly use such indices, as the Standardized Precipitation Index (SPI), the Standard Precipitation-Evapotranspiration Index (SPEI), the Compound Index (CI), the Z Index, the Palmer Drought Severity Index (PDSI), etc. (e.g., Bonaccorso et al. 2003;Vicente-Serrano et al. 2010;Beguería et al. 2014;Jiang et al. 2015;WHO 2016). These index values are indicators of the frequency of specific conditions across extended time-series. In this paper, we performed a statistical analysis of SPI and SPEI, but the developed approach is applicable for other climatic indices.
The SPI is a drought index that is used for estimating wet or dry conditions based on a single precipitation time series (McKee et al. 1993(McKee et al. , 1995Edwards and McKee 1997;Guttman 1998Guttman , 1999Guenang and Kamga 2014 observed precipitation would deviate from the long-term mean. The application of SPI relies on two assumptions: (a) the variability of precipitation is much greater than that of other variables, such as temperature or ET o , and (b) the other variables are effectively stationary, i.e., they have no pronounced temporal trend (Mckee et al. 1993). The SPI indicates the frequency of the precipitation value in the corresponding month as estimated from the entire precipitation record. The World Meteorological Organization (WMO) has recommended SPI as a starting point for meteorological drought monitoring (https://www.drought management.info/indices/). However, SPI cannot reflect how other meteorological parameters, such as temperature or relative humidity, interact with precipitation to influence drought characteristics. The SPEI is generally an extension of the SPI, and has been formulated to capture the effect of the meteorological water balance given as the difference between precipitation  -Serrano et al. 2010a,b). The Global SPEI database, SPEIbase, offers long-term, robust information about drought conditions at global scales, with a 0.5-degree spatial resolution and a monthly time resolution (Global SPEI database 2020). The SPEI can be used for determining the onset, duration and magnitude of drought conditions with respect to normal conditions across a variety of natural and managed systems, with time-scales between 1 and 48 months. Like other climatic indices, SPEI can be used to assess drought severity according to its intensity and duration, and can identify the onset and end of drought episodes. Statistical analysis showed that SPEI is more sensitive to drought assessment than SPI, and could more accurately reflect dry/wet alternations in complex regions (e.g., Stagge et al. 2014;Liu et al. 2021). Taking into account the need to analyze the long-term temporal and spatial climatic changes in the mountainous area, the goal of this paper is to present a statistical framework to: (1) conduct a time-series analysis of meteorological parameters (air temperature, dewpoint temperature, relative humidity, precipitation, and wind speed), using the data for 17 locations of meteorological stations at the East River watershed, Colorado, USA, as a case study, (2) calculate and analyze time series of water balance parameters-potential evapotranspiration, aridity index, actual evapotranspiration, SPI and SPEI, and (3) perform a climatic zonation of the watershed area for different segments of the temporal trends SPEI, by means of the hierarchical k-means clustering and Principal Component Analysis.
The 17 locations of meteorological stations represent significant variance in energy (elevation/aspect), geologic and topographic complexity, subsurface properties, and vegetation, which influence the interactions between precipitation and ET and therefore provide a valuable test case to track the ability of these indices to define climatic zonation across a mountainous watershed. We analyzed the meteorological datasets for the period from 1966 to 2021.
The structure of the paper is as follows: Sect. 2 presents a description of the study area, the sources of meteorological data, and data used for calculations. Section 3 presents a flow chart of the statistical approach, the methods and equations used for calculations of the ET o , ET, AI, SPI and SPEI, as well as the approach used for the hierarchical cluster analysis and ordination. Section 4 includes (a) the results of the time-series analysis to assess the timings of breakthroughs, i.e., abrupt changes or shifts of meteorological and water balance parameters, (b) calculations of SPI and SPEI, and hierarchical clustering, i.e., subdivision/partitioning of meteorological stations. Section 5 includes conclusions and directions of future research.
2 Study area characterization and data used for calculations 2.1 Study area of the East River watershed

Location and topography
The study area is the East River watershed, which is a representative headwater basin in the Upper Colorado River Basin, Colorado, USA (Fig. 1). The watershed is located northeast of the town of Crested Butte, Colorado, and covers an area of * 300 km 2 at an average elevation of 3266 m. The East River watershed is a typical example of watersheds in mountainous areas, with pronounced gradients in hydrology, geomorphology, biome type or life zone (montane, subalpine, alpine). A summary of geographic coordinates, types of the land cover at these locations, and elevations is given in  Figure 2 shows a barplot of the elevations of meteorological stations along with the photographs of three meteorological stations, Almont, representing the lowest elevation, Billy Barr-median elevation, and Mexican Cuthighest elevation. These three locations represent different climatic conditions within the East River watershed in terms of the elevation, a landcover, as well as temperature, precipitation, and other meteorological and water balance parameters, and were used to present typical results of calculations in the main body of the paper. The results of calculations for other locations are given in the attached Supplemental Information (SI).

Climate and biome zones
The East River watershed area is defined as having a continental, subarctic climate with long, cold winters and short, cool summers. Excursions in river discharge are driven primarily by snowmelt in late spring to early summer, with mid-to late-summer monsoonal rainfall inducing rapid but punctuated increases in flow. The watershed comprises several biome zones-montane, subalpine, and alpine life zones, which collectively include aspen, meadow, mixed conifer, sagebrush, willow, grasses, sedges, and a diversity of forbs. The area is characterized by variable climatic conditions across the drainage watersheds, with quite variable rainfall and temperature. The watershed has a mean annual temperature of *0°C, with average minimum and maximum temperatures of -9.2 and 9.8°C, respectively; winter and growing seasons are distinct and greatly influence the hydrology and biogeochemistry of the watershed. An average precipitation is 1200 mm yr -1 , the majority of which falls as snow, which is typical for the Upper Colorado River Basin. A typical feature of the East River watershed is that air reaches the mountains, and then moves up the windward side of a mountain and cools. As a result, humidity increases and orographic clouds and precipitation can develop.

Data used for calculations
Meteorological datasets for the period from January 1966 to December 2021 were downloaded from two major climatic databases: (1) PRISM database (PRISM 2022), which is part of the Northwest Alliance for Computational Science and Engineering at the Oregon State University (the abbreviation PRISM stands for ''Parameter-elevation Regressions on Independent Slopes Model), and (2) NOOA Reanalysis database (2022). These databases provide the most complete meteorological datasets, which were calculated using modern weather forecasting models and data assimilation.
The PRISM database was used to provide time series datasets of monthly precipitation, air temperature (minimum, mean, and maximum), vapor pressure deficit (minimum and maximum), and dewpoint temperature (mean) for the period from January 1966 to December 2021. The wind dataset for the same period (with a 3-h frequency) was downloaded from the NCEP/NCAR Reanalysis database with the application of the function ''NCEP.gather'' in the R package 'RNCEP'' (Kemp 2020). A function ''NCEP.aggregate'' of the package 'RNCEP'' was used to convert the 3-h time series of wind velocity into the monthly wind velocity. (The citations to the R packages that are central to this study are given at the end of the Reference list.) Monthly data of relative humidity were calculated using the air temperature (T a ) and dewpoint temperature (T d ), based on the formula given by (Alduchov and Eskridge 1996;Buck 1981;Lawrence 2005) where air temperature and dewpoint temperature are given in°C, and RH is in percent.
3 Methods of data analysis

Workflow flowchart
The developed statistical framework of the data analysis includes a sequence of three main phases, which are summarized in the flowchart in Fig. 3. Phase I of the data analysis involves an exploration data analysis: plotting time trends of temperature, precipitation, and relative humidity, estimation of the time series structural breakpoints, basic statistics, and cumulative probability distribution functions; Phase II is concerned with calculations of ET o , ET, SPI, and SPEI; and Phase III includes hierarchical k-means cluster analysis of SPEI, averaging of the SPEI time trends for each cluster, and plotting hierarchical clustering trees-dendrograms and 2D PCA (PCA stands for Principal Component Analysis) maps, and zonation maps for periods before and after the time series breaks. Preparation of the datasets for calculations was conducted using the methods given in the paper by Faybishenko et al. (2022). Monthly data of all variables were converted to yearly averaged time series data using the function ''apply.monthly'' of the ''xts'' library in R (Ryan et al. 2020). This function is applied to non-overlapping time periods, such as monthly time series, and it is different from a rolling function in that it used to subset the data based on the specified time period (in this case it is the year), and returns a vector of values for each period in the original data. Other R libraries applied for data manipulation are: zoo (Zeileis and Grothendieck 2005) and dplyr (Wickham et al. 2021).

Detecting breakpoints and segmentation of time series
In statistical studies, time series structural breaks (also termed change-points, or breakthrough points, or simply breakpoints) are commonly defined as unexpected (or sudden) changes over time, i.e., time series shifts (Antoch  (Zeileis et al. 2002). The algorithms is based on the determination of the time series of structural changes, i.e., breakpoints or shifts, of time series based on the evaluation of the F-statistics (Chow test statistics), by means of assessing deviations in the classical linear regression model. In general, the algorithm for computing the optimal breakpoints given the number of breaks is based on a dynamic programming approach of Bellman's ''principle of optimality.'' Details about the underlying theory can be found in Bai and Perron (2003),  3.3 Calculations of ET o , ET, SPI, and SPEI.
The Thornthwaite formula (1948) is given by where T mi is the average monthly temperature ( o C); if the monthly averaged T mi \ 0, then ET o = 0; N is the number of days in the month being calculated, L is the average sunshine hours of the month (which was determined from the site latitude given in Table 1), the exponent a is determined from a ¼ 6:75 Â 10 À7 I 3 À 7:711 Â 10 À5 I 2 þ 1:791 Â 10 À2 I þ 0:4924 with I being the Annual Heat Index given by Thornthwaite equation is an empirical model, and is based on the assumption that temperature and radiation are correlated. However, over the course of a year, air temperature tends to lag behind radiation.
The Hargreaves formula is given by (Hargreaves and Allen 2003) ET o ¼ 0:0135R s T m þ 17:8 ð Þ where T m is monthly temperature ( o C), and R s is the global solar radiation at the land surface, given in the units of water evaporation, mm/day. The Penman-Monteith (PM) equation is based on a combination of turbulent transfer and energy-balance models, and is given by where ET o is the reference evapotranspiration (mm day -1 ), R n is the net radiation (MJ m -2 day -1 ), G is the soil heat flux density (MJ m -2 day -1 ), T is air temperature at 2 m height (°C), u 2 is wind speed at 2 m height (m s -1 ), e s is saturation vapor pressure (kPa), e a is actual vapor pressure (kPa), (e s -e a ) is saturation vapor pressure deficit (kPa), D is slope of the vapor pressure curve (kPa°C -1 ), and c is psychrometric constant (kPa°C -1 ). Note that the Thornthwaite and Hargreaves equations are derived to calculate the potential evapotranspiration, which is the cumulative amount of evaporation and transpiration that would occur if a sufficient water source is available. The Penman-Monteith equation is used to calculate the reference evapotranspiration, which is the amount of evaporation and transpiration from a ''reference surface,'' i.e., a hypothetical grass surface. Calculations of ET o were provided using the library ''SPEI'' in R. The day length is determined using the function ''DayLength'' from the library ''solrad'' (Seyednasrollah 2018).

Calculations of ET using the Budyko model
Budyko (1974) hypothesized that a functional relationship between the actual evapotranspiration, ET, and the two climate variables-precipitation, P, and ET o is given by where AI is the Aridity Index given by The Budyko model is subject to the limiting conditions (Wang et al. 2016;Sposito 2017): under wet conditions, ET is energy limited by ET o , and under dry conditions, ET is water limited by P.
The Budyko equation provides a concise and accurate representation of the relationship between annual evapotranspiration and long-term-average water and energy balance at a catchment scale, and has achieved iconic status in soil and hydrological studies (Sposito 2017).

Calculations of SPI and SPEI
The first step in calculations of SPI and SPEI is to choose a particular probability distribution among several typical distributions for the precipitation P (for calculations of SPI) and the difference (P-ET o ) used for calculations of SPEI. We tested the application of the gamma, Weibull, lognormal, and logistic distributions, and found that the gamma distribution fits best the long-term time series of both precipitation and the difference (P-ET o ).
The SPI and SPEI are usually computed for different time scales, and expressed in units of the standard deviation. Positive SPI and SPEI values indicate wet conditions, and negative values indicate dry conditions (Lloyd-Hughes

Hierarchical k-means clustering analysis and mapping
A clustering analysis is generally provided to classify studied locations on the basis of a set of measured variables into a number of different groups such that similar locations are placed in the same group (Kaufman and Rousseeuw 1990;Maechler et al. 2021). Hierarchical clustering, which is one of the most popular modern clustering techniques, is aimed to build a hierarchy of clusters. The application of this technique is useful as it provides flexibility in selecting from the hierarchical tree (i.e., a dendrogram) how many clusters can be selected and/or how many elements of the studied system can be included in each cluster. In this paper, we applied a univariate clustering of the SPEI time series, assuming that the SPEI is a comprehensive measure of the climatic water balance of the studied area. Hierarchical clustering is provided using the Hartigan and Wong algorithm (Hartigan and Wong 1979), with the application of the function 'hkmeans' in the R library 'factoextra' (Kassambara and Mundt 2020). The results are Fig. 8 Graphs of yearly averaged AI for Almont, Billy Barr, and Mexican Cut. Red color: AI [ 1, which corresponds to water limited conditions, and green color: AI \ 1, which corresponds to energy limited conditions. Vertical dashed lines are 95% CI around the breakpoints. Horizontal dashed blue lines are mean values of AI for each segment Fig. 9 The Budyko curve, i.e., the relationship of ET/P vs AI, with the points representing the yearly averaged ET/P vs AI for Almont, Billy Barr, and Mexican Cut presented in the forms of the cluster dendrogram and the PCA plot for the three segments of SPEI time series, which are shown in Sect. 4.2. The optimal number of clusters was determined using the ''elbow'' method, using the function ''kmeans'' of the R package ''stats.'' Figure SI-2 shows that the optimal number of clusters is 5. The function ''fviz_cluster'' of the R package ''factoextra'' was used to visualize the results. Mapping was conducted with the application of the R libraries ''ggmap,'' ''ggrepel,'' and ''ggplot2'' (Kahle and Wickham 2013;Wickham 2016;Slowikowski 2021).

Meteorological parameters
The temporal warming trends, structural breakpoints and their 95% CIs for Almont, Billy Barr and Mexican Cut locations are shown in Fig. 4, and for all other locations in SI-3. The corresponding boxplots are shown in Fig. 5, showing the minimum, lower-hinge (first quartile), median, upper-hinge (third quartile), and maximum of the temperature before and after the breaks. These results show that the median temperature of the 3rd segment in comparison with the 1st segment increased by 1.8°C at Almont (from 3.3 to 5.1°C), 1.6°C at Billy Barr (from 0.5 to 2.5°C), and 1.5°C at Mexican Cut (from -0.8 to 0.7°C). Figures 6 and 7 show that the precipitation increased in the 2nd segment following the decrease to the same (for Almont) or lower levels (Billy Barr and Mexican Cut) in the 3rd segment. The median precipitation dropped from the 2nd to the 3rd segment at Billy Barr by 306 mm/yr (from 1295 to 989 mm/yr) and at Mexican Cut even more significantly-by 530 mm/yr from 1816 to 1286 mm/yr. The same pattern of increasing and then decreasing precipitation was observed for all stations, except NRCS 1141, where precipitation dropped in the 2nd segment, and then again in the 3rd segment ( Figure SI-4).
Based on the results of the two-sample KS test, the null hypothesis H o that the temperature and precipitation trends for each pair of segments (i.e., for the 1st and 2nd segments, and the 2nd and 3rd segments) are from the same  (Tables 2 and  3). The graphs of the dewpoint temperature are the same for all locations: T d dropped in the 2nd segment compared to the 1st, and then increased to either the same or higher levels in the 3rd segment ( Figures SI-5). Relative humidity decreased in the 3rd segment compared to the 1st segment at all locations ( Figure SI-6).

Time series of ET o , AI, and ET
The ET o was calculated using the Thornthwaite, Hargreaves, and Penman-Monteith equations. For further calculations, we used the Penman-Monteith model that provides more reliable E o , comparable with those determined from field observations by Tokunaga et al. (2016Tokunaga et al. ( , 2019 and numerical CLM (Community Land Model) modeling by Tran et al. (2019). Graphs of the time series of ET o , which were calculated from the Penman-Monteith equation, depict the pattern of decreasing ET o in the 2nd segment, followed by the increase in the 3rd segment to the level above that in the 1st segment for all stations ( Figure SI-7). Fig. 11 Graphs of yearly averaged SPEI (12-month averaged). Blue: wetted condtions, and red: drought conditions. Vertical dashed lines are 95% CI around the breakpoints. Horizontal dashed black horizontal lines are mean values for each segment Fig. 12 Regression analysis of the relationship between SPEI and SPI, showing 95% confidence and 95% prediction intervals Figure 8 shows the time series graphs of AI for Almont, Billy Barr and Mexican Cut. On these plots, the horizontal lines of AI = 1 are separating the periods of limited energy (AI \ 1) and periods of limited water (AI [ 1). For example, at Almont, climatic conditions are characterized as water limited for the period from 1966 to 2021, at Billy Barr-as intermittent energy and water limited until 2011, and then mostly water limited, and Mexican Cut-as energy limited for the whole period from 1966 to 2021, with a few short periods of water limited. Figure  The Budyko curve, i.e., the relationship of ET/P vs AI, with the points representing the yearly averaged ET/P vs AI for Almont, Billy Barr, and Mexican Cut, is shown in Fig. 9. The Almont points are located entirely within the water limited part of the Budyko curve; the Billy Barr points are within the energy limited part before the temporal break, and in the water limited part after the break; and the Mexican Cut points are mostly within the energy limited part with a few points in the water limited part of the Budyko curve. Figure 10 shows the time series graphs of ET for Almont, Billy Barr, and Mexican Cut, and Figure SI-9, the ET graphs for all locations. The graphs show that the temporal patterns of ET are the same for all stations: an increase in the 2nd segment compared to the 1st, followed by the decrease of ET in the 3rd segment practically to the same level as that in the 1st segment. Figure 11 shows the time series of yearly averaged SPEI for Almont, Billy Barr and Mexican Cut, showing the shifts from wetted to drought climatic conditions. Figure SI-10 demonstrates the SPEI graphs for all locations, also indicating the shift from wetted conditions in the 2nd segment to drought conditions in the 3rd segment.

SPEI and SPI time series
Time series plots of SPI are shown in Figure SI-11. A linear regression between SPEI and SPI for all locations is shown in Fig. 12, which illustrates the 95% confidence and 95% prediction intervals for the relationship SPEI = f(SPI). As expected, a prediction interval is wider than a confidence interval. Although the majority of the data points are beyond the 95% confidence interval, most of the points are within the 95% prediction interval. A prediction interval accounts for both the uncertainty in estimating the mean, plus random variations of the individual values. A prediction interval is an estimate of an interval in which a future observation will fall, with a certain probability, and

Hierarchical clustering and zonation mapping
To reduce the effect of the high frequency fluctuations of the graphs of SPEI, hierarchical clustering was performed based on the cumulative distribution curves of SPEI for each of the three segments shown in Fig. 13. Figure 14 shows the results of the hierarchical clustering (left) and corresponding PCA plots for the three segments of SPEI time series for all locations. Based on the results of clustering, Fig. 15 shows the maps of the spatial zonation of meteorological stations. Figures 14 and 15 clearly indicate different clustering and zonation patterns for the three temporal segments. Table 4 lists the clusters for the three temporal segments of SPEI. This table indicates that 7 out 17 locations remained in the same cluster for all three segments, but combinations of the meteorological stations into the clusters are different for the three segments. Thus, the watershed climatic zonation requires periodic re-evaluation based on the results of the structural time series analysis.

Conclusions
A statistical framework has been developed to assess longterm temporal and spatial variability of meteorological parameters (temperature, dewpoint, precipitation, relative humidity, and wind speed), as well as calculated time series of ET o , ET, SPI, and SPEI. Calculations were conducted for the period from 1966 to 2021 for 17 locations of meteorological stations located within the East River Watershed, Colorado, which is a typical watershed in the Upper Colorado River Basin. The statistical analysis included the segmentation of meteorological and calculated time series datasets time series into three segments, based on the evaluation of  structural breakpoints (i.e., shifts) of, and zonation of the studied locations by means of hierarchical and PCA clustering of the SPEI time series segments. Time series segmentation analysis and zonation demonstrate considerable changes in climatic conditions with time and space. The response to changing climate forcing is non-uniform across the watershed. A significant shift in the cluster arrangements for the temporal segments indicate that zonation patterns are driven by dynamic climatic processes, which are variable through time and space. Therefore, the watershed climatic zonation requires periodic re-evaluation based on the structural time series analysis of meteorological data. Examples of directions of future research include the assessment of the redistribution of the water and energy limited regions across the East River watershed, a comparison of the time series breakpoints for different meteorological and water balance parameters, application of other programming languages, such as Python or MATLAB, and the application of the statistical framework to other watershed areas.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.