Temporal Trend Analysis of Meteorological Variables and Reference Evapotranspiration in the Inter-Mountain Region of Wyoming

: Long-term trends in reference evapotranspiration ( ET ref ) and its controlling factors are critical pieces of information in understanding how agricultural water requirements and water resources respond to a variable and changing climate. In this study, ET ref , along with climate variables that directly and indirectly impact it, such as air temperature ( T ), incoming solar radiation ( R s ), wind speed ( u ), relative humidity ( RH ), and precipitation ( P ), are discussed. All variables are analyzed for four weather stations located in irrigated agricultural regions of inter-mountain Wyoming: Pinedale, Torrington, Powell, and Worland. Non-parametric Mann − Kendall (MK) trend test and Theil–Sen’s slope estimator were used to determine the statistical signiﬁcance of positive or negative trends in climate variables and ET ref . Three non-parametric methods—(i) Pettitt Test (PT), (ii) Alexandersson’s Standard Normal Homogeneity Test (SNHT), and (iii) Buishand’s Range Test (BRT)—were used to check the data homogeneity and to detect any signiﬁcant Trend Change Point (TCP) in the measured data time-series. For the data inﬂuenced by serial correlation, a modiﬁed version of the MK test (pre-whitening) were applied. Over the study duration, a statistically signiﬁcant positive trend in maximum, minimum, and average annual temperature ( T max , T min , and T avg , respectively) was observed at all stations, except for Torrington in the southeast part of Wyoming, where these temperature measures had negative trends. The study indicated that the recent warming trends are much more pronounced than during the 1930s Dust Bowl Era. For all the stations, no TCPs were observed for P ; however, signiﬁcant changes in trends were observed for T max and T min on both annual and seasonal timescales. Both grass and alfalfa reference evapotranspiration ( ET o and ET r ) had statistically signiﬁcant positive trends in at least one season (in particular, the spring months of March, April, and May (MAM) or summer months of June, July, and August (JJA) at all stations, except the station located in southeast Wyoming (Torrington) where no statistically signiﬁcant positive trends were observed. Torrington instead experienced statistically signiﬁcant negative trends in ET o and ET r , in the fall months of SON and winter months of DJF. Over the period-of-record, 26, 31, − 48, and + 34 mm ET o and 28, 40, 80, and + 39 mm in an annual time-step reveal that ET o and ET r were signiﬁcantly and positively correlated to T avg , T max , R s , R n , and VPD , as well as signiﬁcantly and negatively correlated to RH . ( RH ), vapor pressure deﬁcit ( VPD ), incoming solar radiation ( R s ), net radiation ( R n ), and wind speed at 2 m height ( u ); (2) investigate the trend change point (TCP) in the long-term climate data series; (3) investigate variation, trends, and magnitude of change in ET ref for grass and alfalfa ( ET o and ET r , respectively) at annual and seasonal time-steps; and (4) investigate relationships between ET ref and associated climate variables. negative trends in relative humidity ( RH ). Our study indicates that minimum T min is rising at a faster rate than the maximum T max . These increases suggest that cooler days and nights, along with frosts, have become less frequent, while hotter days and nights, and heat waves, have become more frequent. The other two study sites—Pinedale and Torrington—exhibited no statistically signiﬁcant trend in annual ET ref for grass ( ET o ) or alfalfa ( ET r ) over the period-of-record. However, a statistically signiﬁcant negative trend was observed for

ET ref has been extensively used to understand the effects of climate change and variability on irrigated agriculture at point and regional scales [10,11]. The general expectation is that temperature increases result in increases in ET ref ; however, many studies report a decreasing trend in ET ref despite increases in air temperature [12][13][14][15]. For example, Chattopadhyay and Hulme [16] studied pan evaporation (E pan ), and Penman [17] estimated ET ref trends in India for 15 years and reported decreasing trends in both E pan and ET ref . Moonen et al. [14] used a temperature-based Hargreaves and Samani [18] model in Pisa, Italy, over a 122-year period and found a decreasing trend in ET ref . Aside from these studies that show regional declines in ET ref , other studies have shown an increase. Tabari et al. [19] reported a positive trend in Penman−Monteith-ET ref over 20 meteorological stations during 1966-2005 in Iran and demonstrated that an increase in air temperature was the main contributing factor. Contrary to an increase in air temperature, Burn and Hesch [20] in their analysis of E pan at 48 stations along the Canadian Prairies found that wind speed has more influence on trends in E pan compared to air temperature. Irmak et al. [21] suggest that one reason for the inconsistent trends reported in the literature is the use of temperature or radiation-based models such as Thornthwaite [22], Blaney and Criddle [23], Priestley and Taylor [24], and Hargreaves and Samani [25]. Since ET ref can be affected by many climatic variables, models that rely on single climate variables, such as air temperature or radiation, can cause significant bias in ET ref . Physically based energy balance methods such as Penman [17], Monteith [26] or Standardized Penman-Monteith [27], which account for different climate variables, can generally provide more robust assessments of the integrated effect of change in meteorological parameters on ET ref . In general, ET ref calculations based on the Standardized Penman-Monteith method required maximum and minimum temperature (T max , and T min ), maximum and minimum relative humidity (RH max and RH min ), incoming solar radiation (R s ), and wind speed at 2 m (u) data. However, in many regions of the world, including the U.S., some of these variables are often missing or not available from the weather station, especially for long-term analysis. Therefore, it is important to drive the methodology to solve a combination-based energy balance equation with fewer inputs to study the span of a diverse array of conditions thought to influence ET ref .
In the semi-arid state of Wyoming, irrigated agriculture accounts for 80-85% of total consumptive water use [28]. Irrigated cropland acreage varies from 0.60 to 0.65 million hectares annually depending upon available surface water supplies. As Wyoming's agricultural, municipal, industrial, and recreational/environmental water uses increase, the State's limited water supplies have come under increased scrutiny. Although the Wyoming State Constitution decrees that all water inside its boundaries is property of the state, several interstate water compacts and court decrees limit the amount of streamflow that Wyoming can deplete [29]. Wyoming serves as a headwater state to the Missouri-Mississippi Basin, Green-Colorado River Basin, Snake-Columbia River Basin, and Great Salt Lake Basin. Wyoming's apportionment of water in the Colorado River Basin (via the Green River and Little Snake River) is also affected by an international treaty between the U.S. and Mexico, as the Colorado River ultimately flows into Mexico [29]. A robust and accurate estimate of ET ref trends can aid in understanding the long-term impacts of climate change and variability on the future of Wyoming's water availability, crop production, and irrigation management. However, few studies quantified ET ref in the inter-mountainous region of Wyoming; this study is one of the first to explore the long-term trends in climate variables and ET ref . Furthermore, there is a need to identify the critical trend change point (TCP) in the climate dataset, which has not been investigated by many researchers. Therefore, the objectives of this study were to (1) investigate long-term trends in maximum, minimum, and average temperature (T max , T min , and T avg ), relative humidity (RH), vapor pressure deficit (VPD), incoming solar radiation (R s ), net radiation (R n ), and wind speed at 2 m height (u); (2) investigate the trend change point (TCP) in the long-term climate data series; (3) investigate variation, trends, and magnitude of change in ET ref for grass and alfalfa (ET o and ET r , respectively) at annual and seasonal time-steps; and (4) investigate relationships between ET ref and associated climate variables.

Study Area and Climate Data
Located in the Rocky Mountain region of the western U.S., Wyoming covers an area of 253,600 km 2 and is dominated by mountain ranges and high elevation plains with an average elevation of~2042 m AMSL ( Figure 1). The central and western parts of the state comprise mountain ranges and basins, while the eastern portion lies within the Great Plains. The state's climate is characterized as semi-arid to arid with cold winters and mild to warm summers, depending on elevation. Annual average precipitation ranges from 915 mm in high elevation mountains to 150 mm in lower elevation plains, with a statewide annual average of 405 mm. Wyoming has approximately 12 million hectares of land in agricultural production, of which 90% is grazing land and approximately 9% is cropland. Most crop production is only possible with irrigation (surface or sprinkler) because rainfall does not typically supply the required amount of water for reproductive or forage growth. Irrigated agriculture accounts for 85% of total consumptive water use in Wyoming. Each year, approximately 8.5 billion cubic meters of surface water diversion and approximately 0.6 billion cubic meters of groundwater withdrawals support the state's 1.4 to 1.6 million hectares of irrigated lands, which generate $1.8 billion of income annually [30].
Water 2020, 12, x FOR PEER REVIEW 4 of 33 variation, trends, and magnitude of change in ETref for grass and alfalfa (ETo and ETr, respectively) at annual and seasonal time-steps; and (4) investigate relationships between ETref and associated climate variables.

Study Area and Climate Data
Located in the Rocky Mountain region of the western U.S., Wyoming covers an area of 253,600 km 2 and is dominated by mountain ranges and high elevation plains with an average elevation of ~2042 m AMSL ( Figure 1). The central and western parts of the state comprise mountain ranges and basins, while the eastern portion lies within the Great Plains. The state's climate is characterized as semi-arid to arid with cold winters and mild to warm summers, depending on elevation. Annual average precipitation ranges from 915 mm in high elevation mountains to 150 mm in lower elevation plains, with a statewide annual average of 405 mm. Wyoming has approximately 12 million hectares of land in agricultural production, of which 90% is grazing land and approximately 9% is cropland. Most crop production is only possible with irrigation (surface or sprinkler) because rainfall does not typically supply the required amount of water for reproductive or forage growth. Irrigated agriculture accounts for 85% of total consumptive water use in Wyoming. Each year, approximately 8.5 billion cubic meters of surface water diversion and approximately 0.6 billion cubic meters of groundwater withdrawals support the state's 1.4 to 1.6 million hectares of irrigated lands, which generate $1.8 billion of income annually [30].  Climate data for the analysis were obtained from four weather stations near Powell, Worland, Pinedale, and Torrington, Wyoming (Table 1; Figure 1). Stations were selected based on data availability, period-of-record (>50 years), and proximity to agricultural lands. The different study period for different stations was selected to quantify the long-term trend and to compare the difference between the current situation and the Dust Bowl Era drought (1930−1935). In addition, the long-term trend for the common period (1963-2017) when data for each station is available were also discussed for potential differences among selected stations. These stations are operated by the National Weather Service (NWS) in manual mode, with observations made and recorded daily by an observer. All of these NWS stations have daily data for T max , T min , T avg , and P. It is important to note, however, that the NWS site at Worland has undergone a few moves during its history and, even more so, was initially installed in an industrial setting. Automated stations were established near the NWS site in Powell in 2010 and near the NWS site in Worland in 2013. For the NWS site in Pinedale, the nearest automated station was installed in 2011, 15 km southeast in Boulder, WY. For the NWS site in Torrington, WY, an automated weather station was installed nearby in September 1996 and then moved to Lingle, WY, (15 km northwest) in May 2008. The automated weather stations are operated and maintained by the Wyoming State Engineer's Office and the University of Wyoming Water Resources Data System's (WRDS) Wyoming Agriculture Climate Network (WACNet). WACNet is an automated climate network that records hourly and daily climate variables [31]. The data from WACNet stations near the four NWS sites were used to model wind speeds. For the NWS climate data, there are data gaps in the observed T max , T min , and P, which were infilled according to the procedures outlined in Cropper and Cropper [32]. These include the following: (i) if data were missing for one day, the missing values were interpolated between the previous and following day; (ii) if temperature data were missing for more than one day and corresponding T min /T max values were available, then missing values were calculated based on the 11-year average of the T min /T max difference equal to the daily diurnal temperature range (DTR) of the missing day; and (iii) if both the T max and T min values were missing for more than one day, the values were calculated based on the 11-year daily average of the missing day. The 11-year average was used such that interpolated values equal the average of the decadal climate.
where ET re f is standardized grass and alfalfa reference evapotranspiration (ET o and ET r ) (mm d −1 ), R n is net radiation at the reference surface (MJ m −2 d −1 ), G is the soil heat flux density (MJ m −2 d −1 ), T is the mean daily air temperature ( • C), u is the mean daily wind speed at 2 m height (m s −1 ), e s is the saturation vapor pressure (kPa), e a is the mean actual vapor pressure (kPa), (e s − e a ) is the vapor pressure deficit (kPa), ∆ is the slope of the saturation vapor pressure curve (kPa • C −1 ), γ is the psychrometric constant (kPa • C −1 ), C n and C d are numerator and denominator constants that change with a reference surface and calculation time-step. C d accounts for the time-step, bulk surface resistance, and aerodynamic resistance of the reference surface. C n and C d coefficients are 900 and 0.34 for ET o , and 1600 and 0.38 for ET r , respectively, over a daily time-step. In this study, measured R s were available beginning in 2011 for the Pinedale station, 2010 for Powell, 1997 for Torrington, and 2013 for Worland. Earlier R s data were calculated on a daily time-step following the Hargreaves and Samani [25] and Samani [33] approach using measured T max and T min data as: where KT is an empirical coefficient, R s is daily extraterrestrial radiation (MJ m −2 d −1 ), and TD is the mean maximum and minimum temperature difference ( • C). Daily R a and other associated variables and parameters were calculated as a function of day of the year, solar constant and declination, and latitude using the procedures presented by Duffie and Beckman [34]: where G sc is the solar constant of 0.0820 MJ m −2 min −1 , ϕ is the latitude (rad), d r is the inverse relative distance from Earth to Sun, which is calculated as a function of day of the year (J) as: ω s is the sunset hour angle (rad), which was calculated as: and δ is the solar declination angle (rad), which was calculated as: The R n was calculated as the radiation balance between incoming net shortwave radiation and outgoing net longwave radiation [27] as: The R ns is calculated as the radiation balance between incoming shortwave radiation and reflected solar radiation as: where α is surface albedo. In this study, a constant value of 0.23 is used for the analysis. Net longwave radiation is calculated as: where σ is Stefan−Boltzmann constant of 4.903E −9 MJ K −4 m −2 d −1 , T max and T min are the daily maximum and minimum absolute temperature in Kelvin, e a is the actual vapor pressure (kPa), R s is incoming shortwave radiation, and R so is the clear sky solar radiation (MJ m −2 d −1 ), which is calculated as: The daily soil heat flux (G) beneath the reference grass and alfalfa surface is generally small, so it was neglected for 24-h time-steps. Actual vapor pressure (e a ) was calculated as: e a = 0.6108 exp 17.27T dew T dew + 237.3 (11) T dew is the dew point temperature ( • C); however, in the absence of T dew data, many studies assume minimum daily temperatures are in equilibrium with T dew . This assumption is an acceptable approximation under clear sky conditions in sub-humid to humid climates. However, in Wyoming's semi-arid climate, this assumption might not always be valid. In this study, T dew was calculated by subtracting a temperature-dependent term from daily minimum temperature (T min ), following the procedure developed by [35]: where T day is diurnal temperature computed from daily T max and T avg as: The coefficient k is a calibration coefficient that accounts for site-specific tendency for RH below 100%. Values for k were computed on a monthly basis using measured RH data for each automated weather station ( Figure 2). Due to a lack of RH data at the Worland station, k values from the Powell station were used to calculate T dew . Missing RH data were calculated as the ratio of actual vapor pressure (e a ) to saturation vapor pressure (e s ) at the same air temperature as follows: Water 2020, 12, x FOR PEER REVIEW 7 of 33 = (0.75 + 2 × 10 −5 ) The daily soil heat flux (G) beneath the reference grass and alfalfa surface is generally small, so it was neglected for 24-hour time-steps. Actual vapor pressure ( ) was calculated as: is the dew point temperature (°C); however, in the absence of data, many studies assume minimum daily temperatures are in equilibrium with . This assumption is an acceptable approximation under clear sky conditions in sub-humid to humid climates. However, in Wyoming's semiarid climate, this assumption might not always be valid. In this study, was calculated by subtracting a temperature-dependent term from daily minimum temperature ( ), following the procedure developed by [35]: where is diurnal temperature computed from daily and as: The coefficient k is a calibration coefficient that accounts for site-specific tendency for RH below 100%. Values for k were computed on a monthly basis using measured RH data for each automated weather station ( Figure 2). Due to a lack of RH data at the Worland station, k values from the Powell station were used to calculate Tdew. Missing RH data were calculated as the ratio of actual vapor pressure ( ) to saturation vapor pressure ( ) at the same air temperature as follows: The was calculated from measured values of and . The slope of saturation vapor pressure versus air temperature (Δ) was calculated as a function of air temperature following [27].  The e s was calculated from measured values of T max and T min . The slope of saturation vapor pressure versus air temperature (∆) was calculated as a function of air temperature following [27].
Wind speed measurement at 2 m height (u, m s −1 ) is an essential component of the standardized Penman−Monteith equation. Measured wind speed measurements at 3 m height, beginning in 2011 for Pinedale, 2010 for Powell, 1997 for Torrington, and 2013 for Worland were converted to u using the logarithmic function [27]. When data were not available, u was estimated following Curtis and Eagleson [36]. Details of this estimating procedure are provided in Appendix A.
To further understand the impact of climate variables on ET ref , linear regression models were developed to evaluate the relationship between ET o and ET r with T max , T min , T avg , P, R s , R n , u, VPD, and RH. For each climate station, the coefficient of determination (R 2 ) was used to measure the proportion of variability in ET o and ET r explained by the climate variables. Furthermore, variations in ET ref can be caused by both radiometric and aerodynamic variables [27]. Changes in the radiative component can be explained by contributions from R n and ∆, whereas changes in the aerodynamic component can be explained by VPD, u, and ∆ contributions. We divided the ET o and ET r equations into their aerodynamic and radiative components, as follows, to evaluate the contribution of each to total ET o and ET r :

Homogeneity Analysis
Three non-parametric methods (i) Pettitt Test (PT) [37], (ii) Alexandersson's Standard Normal Homogeneity Test (SNHT) [38], and (iii) Buishand's Range Test (BRT) [39] were employed to detect any significant TCP in the measured data time-series. All of the three tests used in this study assume the null hypothesis of data being homogeneous; i.e., no change point exists in the time series. These tests were selected because of their different characteristics; for example, the SNHT test is more sensitive to inhomogeneities in the beginning and later end of the time series, whereas PT and BRT tests work better in identifying the TCP near the middle of time series. In this study, the TCP was investigated for only the measured climate dataset i.e., T min , T max , and P time series.

Trend Analysis
In this study, two non-parametric methods, the Mann−Kendall (MK) trend test [40,41] and Theil−Sen's slope estimator [42,43], were used to detect long-term trends in climate variables. The MK test's null hypothesis states that datasets are independent and identically distributed random variables (i.e., there is no trend or serial correlation among the data series). The alternative hypothesis states that the data follow a monotonic trend. However, in many cases, the observed dataset is auto-correlated (e.g., T max and T min ), which may result in misinterpretation of the trend analysis [44,45]. For example, a positive serial correlation results in wrongful rejection of the null hypothesis of no trend when it is true (type I error). Similarly, a negative serial correlation results in accepting the null hypothesis of no trend when it is false (type II error). To eliminate potential issues in time series data, a modified MK test was used by removing the lag-1 serial correlation component from the time series. This process is called pre-whitening. The lag-1 serial correlation of sample data (k = 1) was computed as [44,45]. The "modifiedmk" open-source library package in R-language was used to perform all the non-parametric MK tests, Theil-Sen's Slope test, and all modified versions of MK test mentioned in this study.

Trends in Climate Variables
3.1.1. Maximum, Minimum, and Average Temperatures (T max , T min , and T avg ) Long-term fluctuations in air temperature can cause distinct changes in ET ref , influencing agricultural water consumption patterns and production outcomes. Because of its high elevation, Wyoming has a relatively cool climate. Table 2 presents descriptive statistics of T max , T min , and T avg for Pinedale, Powell, Torrington, and Worland, WY over the study period, and from 1963-2017. Over the period-of-record, a maximum average annual T avg of 8.7 • C (range: −0.01 to 17.4 • C) was observed for Torrington, WY, and a minimum average annual T avg of 2.3 • C (range: −6.4 to 10.9 • C) was observed for Pinedale, WY. At all the stations, January and July were the coldest and warmest months (data not shown). From 1963-2017, T vag was 73%, 14%, and 13% and T max was 36%, 12%, and 9% higher for Torrington compared to Pinedale, Powell, and Worland, respectively. For the same period, minimum T min was 96%, 98%, and 87% lower for Pinedale than forPowell, Torrington, and Worland, respectively. Table 2. Average values with standard deviation for climate variables used in this study from four weather stations over the period-of-record and from 1963-2017. T min , T max , and T avg = minimum, maximum, and average air annual temperature, respectively ( • C); P = average annual precipitation (mm); R s and R n = average annual incoming shortwave radiation and net radiation, respectively (MJ m −2 d −1 ); u = average annual wind speed at 2 m height (m s −1 ), VPD = vapor pressure deficit (kPa); RH = relative humidity (%); MinT min = minimum annual minimum temperature ( • C); MaxT min = maximum annual minimum temperature ( • C); MinT max = minimum annual maximum temperature ( • C); and MaxT max = maximum annual maximum temperature ( • C).  Annual distribution in T max , T min , and T avg for the period-of record and common period (1963-2017) along with trend lines is presented in Figure 3. Over the period-of-record, significant positive long-term trends (p−values < 0.05) in T max , T min , and T avg were observed at all stations, except Torrington, where there were no significant positive trends, and instead a significant negative trend in T max (Table 3). Total change in T max for the period-of-record was 1.2 • C, 0.9 • C, −1.2 • C, and 1.3 • C at Pinedale, Powell, Torrington, and Worland, respectively. Although there is a small decreasing trend in T max at Powell (−0.01 • C year −1 ) from 1963-2017, the recent warming trends (T max ) at Powell and Worland are more pronounced than during the Dust Bowl Era in the 1930s. For example, from 1988 to 2017, a total change of 0.9 • C and 0.3 • C in T max was at Powell and Worland, respectively (Table 4). In addition, the maximum decadal changes in T max occurred between 1981-1990 for Powell (1.8 • C decade −1 ) and Worland (2.5 • C decade −1 ) ( Table 4). This is due to an increasing frequency of the number of very hot days (days with T max > 35 • C) since 1980 at all stations (data not shown). On a seasonal time-step, significant positive trends in T max are observed for March to August, with a maximum change (p-values < 0.05) of +1.5 • C, +1.21 • C, and +1.6 • C for Pinedale, Powell, and Worland, respectively. In contrast, a significant negative T max trend was observed at Torrington for September to February, with a maximum decrease of −1.7 • C (Table 3).     Table 3. Theil-Sen's slope (Qi) of the meteorological variables and grass and alfalfa reference evapotranspiration at annual and seasonal time-step for four weather stations in the state of Wyoming.

Time-
Step   The significant negative T max trend in Torrington on annual and seasonal scales (especially from 1963-2017 (0.05 • C year −1 ) and 1988-2017 (−0.1 • C year −1 ) (Figure 3; Tables 3 and 4) might be due to the cooling effect of irrigation in southeast WY. For example, southeast Wyoming, where most of the agriculture is dependent on the groundwater resources from the Ogallala Aquifer, had a notable increase in groundwater withdrawal for irrigation by 180 million gallons per day from 2005 to 2015 [46]. However, the relationship between irrigation expansion and its effect on local and regional temperature trends is beyond the scope of this study. A similar negative trend in T max was observed by Lobell and Bonfils [47] at the irrigated sites in California during June-August, owing to irrigation expansion in the region. Analysis of long-term climate observations and climate and land surface modeling efforts by Mahmood et al. [48] and Kueppers et al. [49] have shown that irrigation can consistently reduce T max . A similar negative trend in T max was also observed by Sharma et al. [10] in Panhandle, NE, which is in the same climate zone as Torrington, WY. They reported that the negative trend in T max is due to the expansion of irrigated acres in western NE, which may have contributed to increased atmospheric moisture and cloudiness, resulting in a decrease in T max .
Over the period-of-record, statistically significant increases in T min (p-value < 0.05) at the annual time-step were observed for all stations, except Torrington ( Table 3). The degree of observed changes for T min over the period-of-record was 2.20 • C, 0.66 • C, −0.19 • C, and 1.84 • C at Pinedale, Powell, Torrington, and Worland, respectively. Similar to T max , a negative trend in T min from 1963-2017 at Powell is mainly due to a drop in T min from 1977 to 1986, which led to the overall decreasing trend from that period and thereafter , and a significant positive trend was observed (Table 4). For example, the total change in T min from 1988 to 2017 was 1.2 • C, 1.8 • C, −0.09 • C, and 0.6 • C at Pinedale, Powell, Torrington, and Worland, respectively. In addition, for all stations, the largest decadal changes in T min also occurred in the same period, with a total change of +1.50 • C, +1.80 • C, and +2.  (Table 3). Over the period-of-record, T avg showed a significant positive trend at Pinedale, Powell, and Worland and a significant negative trend at Torrington (Table 3). Although a significant increasing trend was observed in T avg in recent years , the decadal trend in T avg revealed that for both Powell and Worland, the warming trend is much more pronounced than in the Dust Bowl Era in the 1930s ( Table 4). The largest change in T avg is +0.04 • C at Pinedale and +0.01 • C year −1 at Powell, both observed in MAM, and −0.02 • C at Torrington and +0.02 • C year −1 at Worland, both observed in DJF.
It is important to note that climate change may not always be characterized by changes in the means of climate variables, but rather by changes in parameters such as standard deviation, skewness, or median and mode [21]. To further understand positive and negative trends in T max and T min , annual minimum and maximum values of T min and T max (MinT min , MaxT min , MinT max , and MaxT max ) were analyzed (Tables 2 and 3). Increases in T max and T min values at Pinedale, Powell, and Worland are accompanied by significant increases in MinT min and MaxT max both over the study period as well as in recent years . This study indicated that MinT min is rising at a faster rate than MaxT max . These increases suggest that cooler days and nights, along with frosts, have become less frequent, while hotter days and nights, along with heat waves, have become more frequent. Significant negative trends in annual T avg and T max at Torrington are due to significant decreases in MinT max , which overshadows the effect of significant positive seasonal trends in MinT min in December through May. In addition, in recent years (2000 to 2017), approximately 59%, 59%, 65%, and 76% of years at Pinedale; 76%, 53%, 65%, and 65% of years at Powell; and 49%, 29%, 59%, and 88% of years at Worland have higher values of MinT min , MaxT min , MinT max , and MaxT max than their long-term average values, illustrating warming trends in recent years.
The climatic forces responsible for warming in WY are complex and are not fully understood. Many factors could lead to the recent warming (e.g., the elevation dependency of warming, which is most likely due to the combined effects of cloud radiation and snow-albedo feedback, changes in atmospheric circulation, and land-use changes). The significant increase in winter and summer air temperature observed in this study is likely to result in increased winter runoff, reduced peak snow-water equivalent, earlier peak streamflow, and reduced summer streamflow volumes [50]. Agriculture is one of the important economic activities in the state. Crops such as alfalfa, sugarbeet, dry bean, and barley grown in Wyoming are highly dependent on the irrigation practices, which is highly dependent on the winter snowpack accumulation. Continued increase in T max and T min in recent years especially in the northwestern and western parts of the state will have a significant effect on the water supply and soil moisture availability in the region, which can further affect the irrigate agriculture.
Similar 21st century warming trends on annual and seasonal time-steps are reported by Frankson et al. [51] in Wyoming with net warming of 1.4 • C since the beginning of the 20th century. They concluded that since 2000, the state has experienced above-average frequency of extremely hot days and below-average frequency of extremely cold days. These results also align with those observed by Gaffen and Ross [52], who reported an increase of over 0.4 • C decade −1 from 1961 to 1995 in the High Plains. Similar results were obtained by Hamlet and Lettenmaier [53] over the western U.S. They used the Historical Climatology Network (HCN) and Historical Canadian Climate Database (HCCD) data, which showed greater increases in minimum temperatures compared to maximum temperatures. A similar observation was made by Mote [54] in the Pacific Northwest U.S. and reported an increase of 0.7 • C to 0.9 • C in annual averaged air temperature from global averages. In the same regions, Abatzoglou et al. [55] observed an accelerated warming rate with linear trends for the 1970-2012 and 1980-2012 time periods of approximately 0.2 • C decade −1 . Temperature trends observed in this study are also consistent with Santos et al. [56]. They used T max and T min temperatures data taken from 35 meteorological (weather) stations across Idaho, U.S. and observed that the difference between T max and T min is decreasing, indicating that the T min is increasing faster than the maximum air temperature. The increasing trend in air temperature observed in the study is also consistent with the findings of Karl et al. [57] and Meehl et al. [58]. The results observed in this study at Pinedale are also consistent with Hall et al. [59], who showed a pronounced statistically significant increasing trend in T max and T min in the spring and summer months at Pinedale, Wyoming.

Precipitation (P)
Changes in P regimes indirectly have a significant effect on other climate variables and thus ET ref . Figure 4 represents the annual variation in P for Pinedale, Powell, Torrington, and Worland. Like other locations in the western U.S., annual P varies greatly by location and from year to year. Being a headwaters state, Wyoming serves as a major source of water for several downstream states, so changes in P can have broad implications beyond the state's boundaries. The long-term average annual P over the study duration was 292 ± 72, 160 ± 46, 351 ± 83, and 199 ± 62 mm at Pinedale, Powell, Torrington, and Worland, respectively, with 51%, 55%, 51%, and 56% of the years having below-average P. Powell and Worland stations in the Big Horn Basin experience less P than the other stations, likely due to the effect of multiple mountain ranges blocking the flow of moist air from both the east and west.
Over the study duration, all stations except Torrington show a positive trend in either annual or seasonal P. There was approximately 12.5, 31.6, −3.9, and 4.1 mm of change in annual P over the study duration at Pinedale, Powell, Torrington, and Worland, respectively. At Powell and Worland, there was a significant positive trend in P during the spring (MAM). Pinedale instead experienced a significant positive trend in P during the fall (SON), as well as a significant negative trend during the summer (JJA). The positive trend in spring P is due to the rising spring temperatures that cause an increase in the snow line (the average lowest elevation at which snowfalls). This further increases the likelihood that some of the snow will fall as rain, which can reduce water storage in the form of the snowpack, particularly at lower elevations. Heavy spring month P, along with increased spring temperatures, can cause severe flooding by causing the rapid melting of the snowpack, further decreasing water availability during the drier summer months [53].

Incoming Shortwave Radiation (Rs) and Net Radiation (Rn)
Incoming solar radiation (Rs) and net radiation (Rn) are excellent estimators of the radiative budget for ETref. However, in many cases, these datasets are not available for long-term analysis. Over the years, many models have been developed to estimate Rs ranging from simple models that are based on sunshine duration, cloud cover, and temperature extremes to more complex models that require detailed meteorological data such as u, RH

Incoming Shortwave Radiation (R s ) and Net Radiation (R n )
Incoming solar radiation (R s ) and net radiation (R n ) are excellent estimators of the radiative budget for ET ref . However, in many cases, these datasets are not available for long-term analysis. Over the years, many models have been developed to estimate R s ranging from simple models that are based on sunshine duration, cloud cover, and temperature extremes to more complex models that require detailed meteorological data such as u, RH, and T dew [62][63][64][65][66]. Sanchez et al. [63] and Sanchez-Lorenzo and Wild [64] estimated R s based on a homogeneous sunshine duration series to study the long-term (1885−2010) trend in R s in Spain and Switzerland, respectively. Biggs et al. [65] estimated R s using a combination of the Angstrom equation and the modified Hargreaves equation to study the trend in R s in southern India from 1985−1997. Similarly, Bandyopadhyay et al. [66] compared different models to estimate R s and observed that the original Hargreaves method performs better than other methods. Therefore, in this study, due to the limitation of R s and R n data availability, R s was modeled on a daily time-step using the Hargreaves and Samani [18] approach. R n was calculated as the radiation balance between incoming net shortwave radiation and outgoing net longwave radiation. It is important to note the modeling of R s from air temperature and extraterrestrial radiation has some error that can impact the trend and magnitude of changes in  Figure S1). The relationship between the estimated and measured R s improved when it is accounted for at monthly and yearly time−step. It was observed that a 10% increase in R s yielded between 2.8% and 3.9% increase in ET o and between 6.5% and 6.7% increase in R n at all sites. A similar observation was made by Irmak et al. [21], who observed that overestimation in R s would impact the R n value portion of the ET ref estimates. They found that an 8% overestimation of R s would result in a 2% and 6.4% overestimation in ET ref and R n , respectively, with all other climate variables held constant. A similar observation was made by Bandyopadhyay et al. [66], who found a relatively small effect of modeled R s on ET o estimates over 29 stations distributed throughout India. While the overestimation in modeled R s may have a minor impact on the magnitude and long-term trend of ET ref , it is not large enough to impact the annual average trend of R n and ET ref in this study.
The annual average R s over the study duration at Pinedale, Powell, Torrington, and Worland was 18.4 ± 0.5, 16.7 ± 0.4, 18.4 ± 0.6, and 17.7 ± 0.4 MJ m 2 d −1 , respectively (Table 2). Overall, annual R s follows a significant positive trend at Powell, and a negative trend at Pinedale, Torrington, and Worland (Table 3). Total change in R s over the period-of-record was −0.5, +0.3, −0.4, and −0.1 MJ m −2 d −1 at Pinedale, Powell, Torrington, and Worland, respectively. On a seasonal scale, statistically significant increases in R s occurred mainly during summer and fall (JJA and SON) ( Table 3). Time series analysis of R s from 1963 to 2017 revealed mixed results with an increasing R s trend in Powell and Worland and a negative trend in Pinedale and Torrington, with a total change of −0.5, +0.2, −1.7, and +0.01 MJ m −2 , d −1 , respectively (Table 4). However, all the positive trends turn negative from 1988 to 2017 time period. R n follows a similar trend as R s , with the exception of Torrington where R s and R n have inverse trends, particularly in the months of SON and DJF. The negative trend in R s and R n at the Pinedale station is characterized by cloudier skies when compared to the Powell and Worland stations, where clear skies increase R s and R n , resulting in higher T max and T min values.

Wind Speed (u)
Wind speed at 2 m height (u) remains relatively constant throughout the period−of−record for all stations (Table 4), with minor seasonal fluctuations (Table 3, Figure 3). It is important to note that although we presented the trends in u, it is not anticipated because in modeling u (based on the availability of station data), the mean and variance remain relatively constant based on time of day, which preserves skewness and distribution [36]. Because the model uses simulated u, independent of other climate variables, the trend in modeled u does not exhibit any bias when compared to measured station data at all locations, indicating the accuracy of the modeling algorithm used. In addition, the estimated u values are comparable to the measured wind speed data at all stations ( Figure 5). In all cases, mean, variance, lag−1 value of autocorrelation, and skewness of simulated u are in agreement with the measured data ( Figure 5). Similarly, positive skewness in the frequency distribution of u was observed by Ivanov et al. [67], with mean observed and simulated u of 3. 49 (Table 2). While our analyses did not consider the annual trends in u, Xie and Zhu [68] observed a decreasing trend in u at higher altitude state of Tibetan Plateau from 1970 to 2009, with u decreasing at a rate of 0.21 m s −1 decade −1 . A similar decreasing trend in u was observed by Yin et al. [69] with u decreasing at a rate of 0.09 m s −1 per decade in China. In the Central High Plain Region of Nebraska, Sharma and Irmak [10] reported an overall reduction of 0.94 and 0.67 m s −1 on annual (January-December) and seasonal (May to September) time-steps, respectively, from 1986 to 2009. In theory, a decrease in u would cause ET ref to decline, with all other climate variables remaining constant; however, this rarely occurs. It was observed with 1% decrease in daily average u yielded between 0.17-0.3% decrease in ET ref . A similar observation was made by Irmak et al. [21], who reported that the effect o u on ET ref can easily be obscured by changes in other climate variables, especially changes in R s . theory, a decrease in u would cause ETref to decline, with all other climate variables remaining constant; however, this rarely occurs. It was observed with 1% decrease in daily average u yielded between 0.17-0.3% decrease in ETref. A similar observation was made by Irmak et al. [21], who reported that the effect o u on ETref can easily be obscured by changes in other climate variables, especially changes in Rs.

Relative Humidity (RH) and Vapor Pressure Deficit (VPD)
The annual change of average relative humidity (RH) and vapor pressure deficit (VPD) over the study duration and for the common period (1963-2017) are presented in Table 3. In this study, both RH and VPD are calculated from Tdew, which in turn is estimated using Equations (3) and (4) and is based on the availability of measured RH data from each station. On a daily time-step, the relationship between estimated and measured Tdew is strong (R 2 = 0.86; RMSD = 4.2 °C for Powell weather station (Figure S2)). In general, relatively low RH values were observed at all stations ( Table  2). On an annual scale, RH values increase over the period-of-record for Pinedale, Torrington, and Worland stations and decrease at the Powell station (Table 3). This trend remains the same from 1963-2017 for all the stations. At Powell, there is a significant decline in RH of 7.3% from 1922−1930 and 4.5% from 1971−1980, followed by an increase of 3.6% from 1988-2017 (Table 4). At the Worland station, RH increased by 5.9% from 1951-1960 but then decreased by 9.4% from 1981−1990. In theory, VPD decreases with increasing precipitation; however, under the semi−arid conditions of Wyoming, a small increase in precipitation has very little effect on VPD relative to the influence of

Relative Humidity (RH) and Vapor Pressure Deficit (VPD)
The annual change of average relative humidity (RH) and vapor pressure deficit (VPD) over the study duration and for the common period (1963-2017) are presented in Table 3. In this study, both RH and VPD are calculated from T dew , which in turn is estimated using Equations (3) and (4) and is based on the availability of measured RH data from each station. On a daily time-step, the relationship between estimated and measured T dew is strong (R 2 = 0.86; RMSD = 4.2 • C for Powell weather station ( Figure S2)). In general, relatively low RH values were observed at all stations ( Table 2). On an annual scale, RH values increase over the period-of-record for Pinedale, Torrington, and Worland stations and decrease at the Powell station (Table 3). This trend remains the same from 1963-2017 for all the stations. At Powell, there is a significant decline in RH of 7.3% from 1922−1930 and 4.5% from 1971−1980, followed by an increase of 3.6% from 1988-2017 (Table 4). At the Worland station, RH increased by 5.9% from 1951-1960 but then decreased by 9.4% from 1981−1990.
In theory, VPD decreases with increasing precipitation; however, under the semi−arid conditions of Wyoming, a small increase in precipitation has very little effect on VPD relative to the influence of and variation in u, T max , and T min . For all stations except Torrington, VPD follows a significant (p-value < 0.1) positive trend on an annual scale over the study duration as well as in recent years (Table 3) (Table 4).

Trend Change Point (TCP) Analysis
In long-term climatic trend analysis, identifying the trend along with the time of trend changes plays an important role. In this study, three non-parametric tests PT, SNHT, and BRT [37][38][39] were employed to detect any significant TCPs in the measured data time-series (T min , T max , and P) on both annual and seasonal scales. On average, similar results were observed from the three tests used in this study. Test results revealed significant TCPs in T max and T min time-series data; however, no TCPs were observed for the P time series at all stations on annual and seasonal scales (Table 5). For Pinedale, significant TCP occurred in 1993 for the annual T min data series. On a seasonal scale, TCP in the T min data series occurred in 1993 for winter months (DJF) and 1980 for summer and fall months (MAM and JJA). At Powell, significant TCP occurred in 1937 and 1951 for T min and T max , respectively, on an annual scale. On a seasonal scale, TCP was mainly observed during the 1950s and 1980s (Table 5). For Torrington, TCP was observed only for T max , with a change point that occurred in the late-1970s and mid-2000s. For Worland, the TCP occurred in the late-1970's for both T max and T min on annual and season scales. As stated by Alexandersson and Moberg [39], there are many reasons for trend change points, such as relocation of the station, instrument/sensor change, land-use, and management changes, or sudden changes in the environmental condition. For example, in this study, the site at Worland has undergone a few moves during its history, especially in late 1970s, which is also reflected in the change point analysis for both T max and T min . Similarly, the recent trend change point at Torrington (mid-2000s) may also be due to a significant increase in water use for irrigation. For example, southeast Wyoming, where most of the agriculture is dependent on the groundwater resources from the Ogallala Aquifer, had a notable increase in groundwater withdrawal for irrigation by 180 million gallons per day (43% increase) from 2005 to 2015 [46].

Trends in Reference Evapotranspiration (ET ref )
Estimates of atmospheric evaporative demand of reference grass and alfalfa surfaces (ET o and ET r ), based on climate station data over the study duration at annual and seasonal scales, are presented in Table 6. It is important to note that estimated climate variables have some impact on the trend analysis and magnitude of changes of ET ref .
To check the accuracy of the ET ref values calculated using the modeled variables, the estimated ET ref values were compared to the ET ref values when calculated using the measured full climate dataset at all four sites ( Figure 6). Overall, a good relationship was observed between two approaches, with R 2 values varying between 0.86 to 0.97 and RMSD ranging from 0.37 to 0.91 mm d −1 . Maximum ET o values were observed in the years 1988, 1931, 1939, and 1981 at Pinedale, Powell, Torrington, and Worland, respectively. Their minimum ET o values were observed in the years 1993, 1924, 2009, and 1932, respectively. Comparing ET r and ET o at each station, average annual ET r values were 26%, 25%, 30%, and 25% higher than ET o at Pinedale, Powell, Torrington, and Worland, respectively. Greater fluctuation in both ET o and ET r was observed at each station, with 49%, 55%, 49%, and 53% of years having above−average ET o for Pinedale, Powell, Torrington, and Worland, respectively, and 51%, 56%, 45%, and 52% of the years had above-average ET r . For each station except Pinedale, maximum and minimum ET o and ET r values occurred in the same years that maximum and minimum R s values occurred.     On a seasonal scale, maximum ET o and ET r values were observed in JJA, and minimum values were observed in DJF (Table 6). At both seasonal and annual time-steps, comparison among the stations for the common time interval (1963 to 2017) reveals larger ET o and ET r values at Torrington, followed by Worland, Powell, and Pinedale. On average, Torrington's annual ET o values were 23%, 22%, and 21% higher than those at Pinedale, Powell, and Worland. Similarly, Torrington's annual ET r values were 27%, 26%, and 26% higher than those of the other stations. Higher ET o and ET r values at Torrington are due to higher values of T max , T avg , R s , u, and VPD compared to the other stations for the same time interval (1963−2017). For example, annual T avg at Torrington was approximately 73%, 14%, and 13% higher than T avg values at Pinedale, Powell, and Worland, respectively. Similarly, average annual VPD was approximately 32%, 13%, and 4% higher at Torrington compared to other stations. On daily time−steps, histograms of ET o and ET r for each station (Figure 7) are positively skewed with the peak (i.e., mode) values between 0 to 1 mm d −1 for ET o and 1 to 2 mm d −1 for ET r . However, at Torrington, peak (i.e., mode) ET o and ET r values are observed between 1 to 2 mm d −1 and 2 to 3 mm d −1 (Figure 7). Average daily values of 3.0, 3.9, 3.0, and 3.1 mm for ET o and 4.1, 4.2, 5.6, and 4.1 mm for ET r were observed at Pinedale, Powell, Torrington, and Worland, respectively.
A time-series graph of grass− and alfalfa-based reference evapotranspiration (ET o and ET r ) on an annual and daily time-step ( Figure 7) suggested positive trends in both ET o and ET r at Pinedale, Powell, and Worland but negative trends at Torrington. Over the study duration, an overall change of +26, +31, −48, and +34 mm in ET o was observed at Pinedale, Powell, Torrington, and Worland. Likewise, an overall change of +28, +40, −80, and +39 mm in ET r was observed. However, long-term trends at the annual scale were statistically significant only for ET o , and only at Powell and Worland (p-value < 0.1). No annual trends were statistically significant at the other stations, nor for ET r . A similar non-significant slight upward trend in ET r was observed by Clifford and Doesken [70] in Colorado from 1993 to 2007. Positive trends in annual ET o at Powell and Worland are associated with significant positive trends in annual P, T max , T min , and T avg , R s , and VPD and a significant negative trend in annual RH.
On a seasonal scale, Powell and Worland showed significant positive trends for both ET o and ET r during MAM, as well as significant positive trends during JJA for ET o (Table 3). Pinedale shared this significant positive trend during JJA, for both ET o and ET r . The total increase during JJA over the study duration was 26, 1.8, and 22 mm in ET o at Pinedale, Powell, and Worland, respectively. These positive trends in ET o and ET r during spring or summer were all associated with significant increases in T max and T avg during those seasons. Other studies have detected significant increases in T avg as the main cause of increasing ET o (e.g., [16]). A similar increasing trend in spring and summer ET o was observed by Hamlet et al. [71] with a median trend of 7-13 mm decade −1 . They associated this increase with a rising vapor pressure deficit from July to September. Additional contributing factors at some of the stations included significant seasonal increases in T min , R s , R n , u, and VPD (Table 3). These trends imply that crop water requirements will likely increase in these regions because of increases in ET o and ET r during the growing season. If accompanied by decreases in summer precipitation, as observed in Pinedale for JJA (Table 3), these two combined effects could have adverse consequences for agricultural production. These increase trends in ET o and ET r , and its adverse consequences to agriculture have also been reported by Jaksa and Sridhar [72]. They used the Noah land surface model to study the long−term trends of ET c within the irrigated region of Snake River basin, Idaho for the past 30 years spanning between 1980 and 2010 and reported that warming climate and boundary-layer meteorological variables indicated that cooling induced by irrigation are the main cause of the increasing trend in ET c (increase at a rate of 0.78 mm year −1 ) in the region.
Torrington, in contrast to all the other stations, showed significant negative trends for both ET o and ET r , specifically during SON and DJF. These long-term decreases at Torrington may be the result of significant decreases in seasonal T max because of the cooling effect of irrigation [47][48][49]. In addition, a decreasing trend in ET ref may be due to the decreasing trend in P. Trend analysis from 1963 to 2017 (red and blue lines in Figure 7) revealed that both ET o and ET r are decreasing at a much faster rate at Powell and Torrington but increasing at Pindale and Worland. A total change of +26, −34, −187, and +35 in ET o and +28, −43, −291, +42 mm in ET r was observed at Pinedale, Powell, Torrington, and Worland, respectively, from 1963-2017. Torrington's trends are consistent with the findings of Sharma et al. [10], who reported similar negative trends from 1986−2009 in Panhandle, NE, which is in the same climate zone as Torrington, WY. A similar negative ET o trend was observed in China [73], India [16], and the U.S. [21]. All these studies reported that the negative trends in T max are the prominent cause of reduction in ET o .  (Table 4). This decade also saw significant declines in u and VPD and a significant increase in RH. Powell, which is just 120 km northwest of Worland, saw similar conditions during this decade, but the reductions in ETo and ETr were not statistically significant.  (Table 4). This decade also saw significant declines in u and VPD and a significant increase in RH. Powell, which is just 120 km northwest of Worland, saw similar conditions during this decade, but the reductions in ET o and ET r were not statistically significant.
Furthermore, the trend and contribution of aerodynamic and radiative components of ET o and ET r revealed that, for both ET o and ET r , the aerodynamic component contributes more at all stations, ranging between 45% to 59% of the total ET o and 59% to 79% of the total ET r ( Table 7). The maximum aerodynamic contribution to ET o and ET r was observed at Torrington. Trend analysis of these components illustrates a positive trend in aerodynamic and radiation components at all stations, except Torrington, where both components show negative trends. The positive trends of aerodynamic and radiative components at Pinedale, Powell, and Worland are driven by the increasing trend of radiative components (i.e., solar radiation and air temperature) and positive trend of aerodynamic components (i.e., u and VPD). Negative trends in both radiative and aerodynamic components of ET o and ET r at Torrington are mainly driven by negative trends in both temperature and u. Moreover, the negative trend of u at the Torrington station was more evident than at the other three stations, which have a positive ET o trend. Daily estimates of ET o and ET r were analyzed for their potential to inform future estimates of actual crop ET c based on the two-step approach of FAO-56 for calculating crop coefficients (K c ). Values for the grass crop coefficient (K co ) and alfalfa crop coefficient (K cr ) are typically obtained by dividing ET c by ET o and ET r , respectively. Since K co and K cr coefficients are not interchangeable, this study analyzed the relationships between daily ET o and ET r ( Figure S3). The results indicated that, in semi-arid WY, the ratio of ET o /ET r ranged from 0.70 to 0.75. The minimum value was observed for Torrington, while the other stations' ET o /ET r values were very close to 0.75. These values are consistent with those of Irmak et al. [74], who evaluated ET r /ET o (reciprocal to that presented in this study) in other climate regions. For example, the values at Pinedale, Powell, and Worland closely match the values Irmak et al. [74] reported for Phoenix, AZ, which is categorized as arid, while the Torrington ratio resembled values reported for semi-arid Bushland, TX.

Relationship between ET r and ET o and Climate Variables
Small changes in meteorological variables can result in significant changes in ET ref ; therefore, these relationships should be further evaluated. Linear regression models were developed to evaluate the relationship between ET o and ET r with T max , T min , T avg , P, R s , R n , u, VPD, and RH on an annual time-step over the period-of-record (Table 8). Positive and negative signs in Table 8 indicate whether the relationship is positive or negative between a climate variable and ET o and ET r . Results show that, for all stations, nearly all climate variables have a significant correlation with ET o and ET r (p−values < 0.01), with the exception of T min , u, and P (where they are grayed out in Table 8). Most variables are positively correlated with ET o and ET r , except P and RH, which show a negative correlation. T max , R s , and VPD explain the maximum amount of variation in ET o and ET r values across all four stations. Results for Pinedale, for example, indicate that VPD explains approximately 87% of variability in ET o and 83% of variability in ET r , followed by T max (0.70 in ET o and 0.67 in ET r ) and R s (0.66 in ET o and 0.67 in ET r ). T min and u are poorly correlated to ET o and ET r (Table 8). Similar observations occur at Torrington, which is at a similar latitude as Pinedale, with VPD, R s , and T max accounting for~80% of the variability in ET o and ET r . Negative relationships between ET ref and P at Pinedale and Torrington are probably due to higher cloud cover there and decreases in T max , R s , R n , and VPD, all of which are generally associated with increasing P.
At Powell and Worland, in the northwestern portion of the state, the effects of VPD, R s , and T max on ET o and ET r are somewhat less than at the other stations (Table 8). VPD, T max, and R s explain approximately 81%, 59%, and 51% of the variability in ET o at Powell, and 73%, 54%, and 50% of the variability in ET r at Powell. They explain 75%, 71%, and 54% of the variability in ET o at Worland and 68%, 66%, and 59% of the variability in ET r at Worland. Another difference observed at Powell and Worland is the significant correlation of T min and u with ET o and ET r , which explain 10−16% of the variability in ET o and ET r . The negative relationship between ET ref and P is probably due to an increase in cloud cover and decrease in T max , R s , R n , and VPD, which is generally associated with increasing P. Significant relationships between u and ET ref at Powell and Worland, but not Pinedale and Torrington, are likely due to the more arid climate of the northwestern region. In general, u influences ET ref through its effects on moisture content of the air. For example, because of the arid climate of northwestern Wyoming, the saturated air above the surface is replaced by drier air, which increases VPD and hence ET ref . Table 8.
Coefficient of determination (R 2 ) observed between grass and alfalfa reference evapotranspiration (annual ET o and ET r ) and climate variables. T min , T max , and T avg = minimum, maximum and average air annual temperature ( • C); P = average annual precipitation (mm); R s and R n = average annual incoming shortwave radiation and net radiation (MJ m −2 d −1 ); RH = relative humidity (%); u = average annual wind speed at 2 m height (m s −1 ); VPD = vapor pressure deficit (kPa); (+) and (−) signs indicate a positive or negative relationship between climate variables and ET o and ET r ). These findings are consistent with results reported by Maček et al. [75], who examined the trends and relationships between climate and ET o from 18 meteorological stations located in the temperate continental to mountain climate of Slovenia. The authors reported that R s had the largest impact on ET o , followed by T avg and VPD. These results are also in agreement with Irmak et al. [76], who evaluated different climate regions and reported that ET o is highly sensitive to VPD. It should be noted that changes in magnitude or trend of one climate variable cannot solely describe the variation in ET ref . The results indicated that ET ref did not respond well to T min ; however, a significant increase in T min can indirectly increase trends in T avg , which in turn affects VPD and ET ref . Therefore, it is important to account for all climate variables collectively when analyzing variation in ET ref by using combination-based energy balance models as opposed to models based solely on simplified temperature or radiation components. In addition, this also makes it important to understand and quantify the individual contribution of the climatic variables to the changes in ET o and ET r . Therefore, future studies will aim at the sensitivity analysis of ET o and ET r to the changes in climate variables. Future studies will also aim to evaluate the spatio−temporal variation in ET ref and ET a in the inter-mountain region of Wyoming.

Conclusions
Information on long−term trends in climatic variables, reference evapotranspiration (ET ref ), and their relationship to each other is vital in the management of agricultural, water resources, and hydrological systems under variable and changing climate conditions. In this study, long−term trends for reference evapotranspiration (ET ref ) and climate variables that directly and indirectly impact ET ref were evaluated over four sites in the intermountain state of Wyoming: Pinedale (1963-2017), Torrington (1922-2017), Powell (1908-2017), and Worland . In addition, the long-term trend for the common period (1963-2017) when data for each station areavailable was also discussed for potential difference among selected stations. ET ref was calculated on a daily time-step using the combination-based energy balance Standardized Penman-Monteith equation. For all sites, the full climate dataset (T max , T min , RH, u, and R s ) to estimate ET ref was available only for recent years. Earlier climate data were estimated using a previously developed practical and robust procedure using the observed T max and T min data on a long-term basis. It is important to note that estimated climate variables have some impact on the trend analysis and magnitude of changes; however, the impact is not expected to be large enough to influence the results on this study, especially when annual scales are considered. Non-parametric Mann-Kendall trend test and Theil-Sen's slope estimator were used to determine the statistical significance of positive or negative trends in climate variables and ET ref . Three non-parametric methods, (i) Pettitt Test (PT), (ii) Alexandersson's Standard Normal Homogeneity Test (SNHT), and (iii) Buishand's Range Test (BRT), were employed to detect any significant TCP in the measured data time-series. For the data influenced by serial correlation, a modified version of the MK test (pre-whitening) was applied.
A statistically significant positive trend in maximum, minimum, and average annual temperature (T max , T min , and T avg , respectively) was observed at all stations, except for Torrington in the southeast part of Wyoming, where these temperature measures had negative trends. The study indicated that the recent warming trends are much more pronounced than during the 1930s Dust Bowl Era. For all the stations, no TCPs were observed for P; however, significant changes in trends were observed for T max and T min on both annual and seasonal timescale (Table 5). Over the study duration, an overall change of +26, +31, −48, and +34 mm in ET o and +28, +40, −80, and +39 mm in ET r was observed at Pinedale, Powell, Torrington, and Worland, respectively. Two of the four study sites-Powell and Worland-exhibited a statistically significant positive trend in annual ET ref for grass (ET o ), but not for alfalfa (ET r ), over the period−of−record (Table 3). These increases were associated with significant positive trends in annual T max , T min , and T avg and significant negative trends in relative humidity (RH). Our study indicates that minimum T min is rising at a faster rate than the maximum T max . These increases suggest that cooler days and nights, along with frosts, have become less frequent, while hotter days and nights, and heat waves, have become more frequent. The other two study sites-Pinedale and Torrington-exhibited no statistically significant trend in annual ET ref for grass (ET o ) or alfalfa (ET r ) over the period-of-record. However, a statistically significant negative trend was observed for both ET o and ET r at Torrington from 1963-2017, where ET o and ET r are decreasing at a much faster rate of 3.4 mm year −1 and 5.3 mm year −1 compared to other stations. That said, three stations-Pinedale, Powell, and Worland-exhibited statistically significant positive trends in seasonal ET ref during the spring (MAM) or summer (JJA). Torrington showed no statistically significant positive seasonal trends, instead experiencing statistically significant negative trends during the fall (SON) and winter (DJF) seasons. These long−term decreases in fall and winter ET o and ET r are associated with significant decreases in T max , T avg , R s , and in some instances decreased u and VPD and increased RH.
On an annual time−step, ET o and ET r are significantly and positively correlated with T avg , T max , R s , R n , and VPD and significantly negatively correlated with RH (Table 8). Such increases in ET o and ET r , if combined with a decrease in summer (JJA) precipitation (e.g., at Pinedale), may have serious adverse consequences for agricultural production. Increases in ET o and ET r trends in MAM and JJA months could be attributed to increases in crop water requirements for the production of various crops and early snowmelt in mountainous regions that determine seasonal water availability. Therefore, it is important for state agencies, irrigation districts, and water managers to consider these rising trends when discussing future developments for water resources and agriculture in Wyoming, particularly given long-term warming trends in the state. Author Contributions: This paper is a joint effort by several authors. V.S. conceptualized the idea, obtained funding, performed the data analysis, and drafted the manuscript. C.N. and A.B. assisted V.S. in data collection, data analysis, and editing the manuscript. S.I. and D.P. contributed to the intellectual content of the manuscript and provided edits. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
The authors express their appreciation to USDA−NIFA for funding. We are also grateful to University of Wyoming Water Resources Data System for making the long-term climate data and station information available for the study. The mention of trade, firm, or product names is for descriptive purposes only and does not imply endorsement, suitability, or recommendation for use by the U.S.−government−represented universities or the authors.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Methods for Estimating Wind Speed (u)
When wind speed data were unavailable, u was estimated following Curtis and Eagleson [36]. This procedure assumes that the frequency distribution of u is positively skewed, and u is treated as lag−1 Markov process, independent from other climate variables. Furthermore, the mean and variance of the observed u values are not allowed to vary with time of day when estimating u, which ensures that the generated values approximate the property of atmospheric stability following a characteristic diurnal curve. This procedure is outlined as (i) determining statistical parameters (mean, variance, skewness) of u from the measured dataset, (ii) analyzing the frequency distribution of u from the measured dataset, and (iii) generating random numbers of the same statistical characteristics and frequency distribution, then constituting the deterministic part and combining it with the random distribution. To generate skewed daily u data, Curtis and Eagleson [36] suggest that the random term should force skewness on the estimates of the autoregressive model, leading to an approximately gamma distribution of u: u (t−1) = u + ρ s u (t) − u + ε (t−1) σ s 1 − ρ 2 s (A1)