## Abstract

In the remote and challenging terrain of the Himalayan region, accurate measurement of cyclic snow accumulation and depletion is a significant challenge. To overcome this, an attempt has been made in the present study by applying a statistical analysis of MODIS snow time series data with the Seasonal Autoregressive Integrated Moving Average (SARIMA) model from 2003 to 2018 over the Beas river basin. The Box–Jenkins methodology of forecasting is based on the identification using seasonality, stationarity, ACF, and PACF plots; and estimation based on maximum likelihood techniques; and the last diagnostic checking based on the residual and error values have been used. Later, forecasting models have been proposed separately for the snow accumulation period (October–February) as (1,1,1) (0,1,3)^{19} and for the snow depletion period (March–September) as (1,1,1) (1,1,2)^{27} after calibration of the data (2003–2015) and the same were then validated using data (2016–2018). The accuracy assessment of the models has been checked using performance criteria like AIC, MSE, and RSS. The comparison of the forecasting models with the observed data showed a good agreement with *R*^{2} of 0.83 and 0.89 for snow accumulation and snow depletion, respectively. This research highlights the potential of utilizing satellite data and statistical modeling to address the challenges of monitoring snow cover in remote and inaccessible regions.

## HIGHLIGHTS

A methodology has been proposed to obtain the snow accumulation and depletion of the Beas river basin in the western Himalayas.

A time-series data of MODIS terra and aqua optical satellite data has been used from 2003 to 2018 to find out the snow cover area.

The study showed that the models have good agreement between the simulated and observed data in snow accumulation and snow depletion.

## INTRODUCTION

A snow-covered area (SCA) is a hydrological feature that applies to a variety of hydrological applications on a basin scale, which could benefit both local development and water resource management (Şensoy & Uysal 2012). There are usually practical applications such as assessment of snow reserve, simulation of snowmelt runoff through snowmelt runoff models, forecasting floods, flood control, irrigation, electricity production, and making operational decisions regarding reservoir management. In snowmelt-dominated regions, SCAs are useful for understanding, modeling, and predicting atmospheric, hydrological, and ecological processes (Zhou *et al.* 2005; Richer 2009). The region contributes almost 60% of annual streamflow received by melting snow and glaciers to the Beas river as fresh storage of water (Ahluwalia *et al.* 2013). The Beas basin receives heavy snow during the winter (October–February) and received plenty of muddy water during the monsoon/summer (March–September) (Ahluwalia *et al.* 2013). This melted water fulfils the demands of millions of people in various forms such as irrigation, drinking water, fishing, hydroelectric demands, etc. (Colbeck 1978; Kulkarni *et al.* 2010; Lutz *et al.* 2016). On the other hand, the outbreak of glacier lakes and inaccurate assessment of seasonal snowmelt in the form of avalanches may cause natural hazards that damage the life of mankind and infrastructures. Despite potential contribution to the development of water resources management and damages in the form of lives and infrastructure of society, a few studies have been explored related to the accurate assessment of SCA (Singh *et al.* 1995, 1997; Kumar *et al.* 2005; Prasad & Roy 2005; Haritashya *et al.* 2006; Hall *et al.* 2010; Singh & Jain 2013; Han *et al.* 2019; Di Marco *et al.* 2021). Thus, the long-term monitoring and forecasting of glaciers and snowfields are pre-requisite for proper water management and disaster preparedness (Prasad & Roy 2005). Various studies have been performed for rainfall–runoff modelling by calibration and validation of the Curve Number (CN) values for watershed hydrologic response (Rautela *et al.* 2023; Rizul *et al.* 2023). In the Himalayan region, spatial snow variability and accurate assessment of snow accumulation and depletion is a challenging task because of its highly complex topography, rugged terrains, and inaccessibility (Immerzeel *et al.* 2009). To overcome this, the advanced Remote Sensing (RS) and Geographic Information System (GIS) techniques is the only possible solution to tackle such challenges without contact with the actual snow cover region (Dhanju 1983; Matson *et al.* 1986; Dozier & Marks 1987; Hall *et al.* 2002; Dozier & Painter 2004; Jain *et al.* 2010; Kulkarni *et al.* 2010; Sichangi *et al.* 2016; Steele *et al.* 2017; Xiang *et al.* 2019). There are numerous difficulties in data acquisition techniques of optical RS like cloud persistence and hill shadows. In the present study, MODIS 8-day average cloud-Free snow cover composite products of MOD10A2 (Terra) and MYD10A2 (Aqua) consisting of a single tile of 1,200 km × 1,200 km area on earth having a spatial resolution of 500 m from 2003 to 2018 were downloaded and converted it to a useful product. Then, the same time series data were used to apply statistical techniques to forecast a model for the prediction of snow accumulation and depletion (Raicharoen *et al.* 2003). The time series data of MODIS-NDVI have also been utilized to analyse the impact of climate change on crop productivity (Jabal *et al.* 2022). Moreover, the time series data of monthly rainfall have been used to analyse the temporal change in the rainfall and temperature trend over the western ghat using EI-Nino and La-Nina techniques (Mann & Gupta 2022). Box (1976) proposed simple techniques such as the ARIMA (AutoRegressive Integrated Moving Average) model which follows the normal statistical distribution, and Seasonal AutoRegressive Integrated Moving Average (SARIMA) model which accounts for seasonality and periodic fluctuations (Box & Jenkins 1976; Adhikari & Agrawal 2013; Box *et al.* 2015; Khandelwal *et al.* 2015; Astuti & Jamaludin 2018; Adams *et al.* 2019; Farsi *et al.* 2020). The ACF and PACF plots of satellite snow data clearly show the cyclic seasonal variation as maximum during the peak winter and minimum at the end of the monsoon period. Thus, the SARIMA model has been proposed as the most efficient model to predict and forecast future snow accumulation and depletion. Furthermore, snow cover accumulation and depletion estimations based on both satellite data and numerical weather prediction (NWP) data, which are extremely rare for mountainous regions, will also be used in real-time runoff forecasting.

## MATERIALS AND METHODS

### Study area and data used

*et al.*2018). The Hierarchical Data Format of the Earth Observation System (HDF-EOS) format was converted into the original geographic latitude–longitude projection by the nearest neighbor resampling method with the application of the ERDAS Imagine-14 software. Model maker and interpreter tool are used to classify geometrically corrected MODIS data and finally, it has been superimposed on the classified DEM of ASTER 12m. For reducing the underestimation of SCA, mainly caused by cloud cover, the seasonal, temporal, and spatial filters have been applied; and for reducing the overestimation, the combined Terra and Aqua product of snow cover considers snow only if a pixel representing snow in both the products has been used (Muhammad & Thappa 2020).

### SARIMA methodology

*k*; and the Partial Autocorrelation Function (PACF) to decide which autoregressive terms should be involved to expose the time lags based on plotting residuals. In the present study, the SARIMA model based on non-seasonal model ARIMA (

*p, d, q*) and SARIMA (

*P, D, Q*)

*was considered for the prediction and forecasting of snow accumulation and depletion on the seasonal and non-stationary time series. The structure of the SARIMA model is denoted as (*

^{S}*p, d, q*) (

*P, D, Q*)

^{S}and presented in the form of an equation as Equation (1) (Hazarika

*et al.*2020).where represents Gaussian white noise at period

*t*,

*X*as time series data at period

_{t}*t*, and is the constant. ∇

*is the non-seasonal difference component as (1 −*

^{d}*B*)

*, is the seasonal difference component as (1 −*

^{d}*B*)

^{S}*,*

^{D}*B*as the backward shift operator as

*S*represents as seasonal order. The non-seasonal autoregressive operator and moving average operator in Equation (1) represent in form of polynomials as (

*B*) and

*θ*(

*B*) as order

*p*and

*q*, respectively, in Equations (2) and (3):

By considering the relationship within the data, the different SARIMA models in (*p, d, q*) (*P, D, Q*)* ^{S}* were successfully applied. The

*d*and

*D*indicate the order of the non-seasonal and seasonal differencing more than 1 and 2 total seasonal differences. The period value of time series

*S*(seasonality) is based on the dataset as weekly, monthly, yearly, and other seasonal data, respectively.

### Model selection criteria

*et al.*2020). The model identification involves differencing of the time series along with the behavior of the ACF and the PACF to recognize the temporal relations of the transformed data. The ACF function checks the relationships of the previous values with the latest values and the PACF function checks the value of the correlation coefficient between the variable and its time lag (Rebala

*et al.*2019). The parameter estimation of the fitted model involves finding the parameters that maximize the probability of observations using maximum likelihood techniques. For checking the adequacy of the model, the method applied as Akaike's Information Criterion (

*AIC*) (Akaike 1974; Eni & Adeyeye 2015; Chen

*et al.*2018) and mathematical formulation are given in Equation (6) aswhere

*n*refers to the size of the series, and

*p, q, P,*are the parameters of the selected model. The optimal forecast model has been selected as the one having the minimum AIC value of the group (Box

*and*Q*et al.*2015). In the present study, the criteria of mean square error (

*MSE*) and the residual sum of square error (

*RMSE*), and

*AIC*were applied.

## RESULTS AND DISCUSSION

### Model identification

### Stationarity and seasonality test

### Model estimation

Based on the observation in a section of model identification, the final model estimation has been done by computing the seasonal and non-seasonal variations of ACF and PACF plots which were followed by combining these models as the following.

### Seasonal and non-seasonal variation (snow accumulation)

^{19}and put into Equation (1) and the final model has been proposed as Equation (7):

### Seasonal and non-seasonal variation (snow depletion)

^{27}and the same has been put into Equation (1) and the final model has been drawn as Equation (8):

The next step is to find the appropriate values of parameters by searching the models around the above-proposed model based on the selection criterion as lower values of RSSE, MSE, and AIC. The number of models has been searched with different combinations of *p, d, q, P, D*, and *Q* for their performance using selection criteria (Tables 1 and 2). The best-fit models were found as (1,1,1) (0,1,3)^{19} for snow accumulation and (1,1,1) (1,1,2)^{27} snow depletion and the values of parameters are shown in Table 2 accordingly. Table 3 represents the estimated values of the parameters based on the selection criteria.

Model (SARIMA) . | RSSE . | MSE . | AIC . |
---|---|---|---|

(0,1,1) (1,1,0)^{19} | 54,270.3 | 242.278 | 242.28 |

(0,1,1) (0,1,1)^{19} | 32,702.4 | 145.993 | 145.99 |

(0,1,1) (1,1,1)^{19} | 31,251.3 | 140.14 | 140.14 |

(0,1,1) (1,1,2)^{19} | 31,492.6 | 141.859 | 141.86 |

(0,1,1) (0,1,2)^{19} | 31,074.4 | 139.347 | 139.35 |

(0,1,1) (1,1,3)^{19} | 29,868.7 | 135.153 | 135.15 |

(0,1,1) (0,1,3)^{19} | 29,531.7 | 133.026 | 133.03 |

(0,1,1) (1,1,4)^{19} | 28,591 | 129.959 | 129.96 |

(0,1,1) (0,1,4)^{19} | 30,250.5 | 136.88 | 136.88 |

(1,1,1) (1,1,0)^{19} | 44,624.8 | 200.111 | 200.11 |

(1,1,1) (0,1,1)^{19} | 29,641.5 | 132.921 | 132.92 |

(1,1,1) (1,1,1)^{19} | 27,810 | 125.27 | 125.27 |

(1,1,1) (0,1,2)^{19} | 27,198.1 | 122.514 | 122.51 |

(1,1,1) (1,1,3)^{19} | 26,893.6 | 122.244 | 122.24 |

(1,1,1) (0,1,3)^{19} | 26,315.3 | 119.074 | 119.07 |

(1,1,1) (0,1,4)^{19} | 27,032.7 | 122.876 | 122.88 |

(1,1,1) (1,1,4)^{19} | 26,864.7 | 122.67 | 122.67 |

(2,1,1) (1,1,0)^{19} | 44,625.1 | 201.014 | 201.01 |

(2,1,1) (0,1,1)^{19} | 29,624.7 | 133.445 | 133.45 |

(2,1,1) (1,1,1)^{19} | 27,369.1 | 123.842 | 123.84 |

(2,1,1) (1,1,3)^{19} | 26,835.2 | 122.535 | 122.54 |

(2,1,1) (1,1,4)^{19} | 26,874.1 | 123.276 | 123.28 |

(2,1,1) (0,1,2)^{19} | 27,207 | 123.108 | 123.11 |

(2,1,1) (0,1,3)^{19} | 26,273.8 | 119.427 | 119.43 |

(2,1,1) (0,1,4)^{19} | 27,046.4 | 123.5 | 123.50 |

Model (SARIMA) . | RSSE . | MSE . | AIC . |
---|---|---|---|

(0,1,1) (1,1,0)^{19} | 54,270.3 | 242.278 | 242.28 |

(0,1,1) (0,1,1)^{19} | 32,702.4 | 145.993 | 145.99 |

(0,1,1) (1,1,1)^{19} | 31,251.3 | 140.14 | 140.14 |

(0,1,1) (1,1,2)^{19} | 31,492.6 | 141.859 | 141.86 |

(0,1,1) (0,1,2)^{19} | 31,074.4 | 139.347 | 139.35 |

(0,1,1) (1,1,3)^{19} | 29,868.7 | 135.153 | 135.15 |

(0,1,1) (0,1,3)^{19} | 29,531.7 | 133.026 | 133.03 |

(0,1,1) (1,1,4)^{19} | 28,591 | 129.959 | 129.96 |

(0,1,1) (0,1,4)^{19} | 30,250.5 | 136.88 | 136.88 |

(1,1,1) (1,1,0)^{19} | 44,624.8 | 200.111 | 200.11 |

(1,1,1) (0,1,1)^{19} | 29,641.5 | 132.921 | 132.92 |

(1,1,1) (1,1,1)^{19} | 27,810 | 125.27 | 125.27 |

(1,1,1) (0,1,2)^{19} | 27,198.1 | 122.514 | 122.51 |

(1,1,1) (1,1,3)^{19} | 26,893.6 | 122.244 | 122.24 |

(1,1,1) (0,1,3)^{19} | 26,315.3 | 119.074 | 119.07 |

(1,1,1) (0,1,4)^{19} | 27,032.7 | 122.876 | 122.88 |

(1,1,1) (1,1,4)^{19} | 26,864.7 | 122.67 | 122.67 |

(2,1,1) (1,1,0)^{19} | 44,625.1 | 201.014 | 201.01 |

(2,1,1) (0,1,1)^{19} | 29,624.7 | 133.445 | 133.45 |

(2,1,1) (1,1,1)^{19} | 27,369.1 | 123.842 | 123.84 |

(2,1,1) (1,1,3)^{19} | 26,835.2 | 122.535 | 122.54 |

(2,1,1) (1,1,4)^{19} | 26,874.1 | 123.276 | 123.28 |

(2,1,1) (0,1,2)^{19} | 27,207 | 123.108 | 123.11 |

(2,1,1) (0,1,3)^{19} | 26,273.8 | 119.427 | 119.43 |

(2,1,1) (0,1,4)^{19} | 27,046.4 | 123.5 | 123.50 |

Model (SARIMA) . | RSSE . | MSE . | AIC . |
---|---|---|---|

(0,1,1) (1,1,0)^{27} | 27,605.2 | 86.27 | 1,536.11 |

(0,1,1) (1,1,1)^{27} | 18,914.2 | 59.29 | 1,405.40 |

(0,1,1) (1,1,2)^{27} | 18,641 | 58.62 | 1,402.29 |

(0,1,1) (1,1,3)^{27} | 19,518.4 | 61.57 | 1,420.43 |

(0,1,1) (0,1,0)^{27} | 33,912.8 | 105.65 | 1,606.34 |

(0,1,1) (0,1,1)^{27} | 1,9019.3 | 59.44 | 1,405.34 |

(0,1,1) (0,1,2)^{27} | 19,004.4 | 59.58 | 1,407.07 |

(0,1,1) (0,1,3)^{27} | 20,185.4 | 63.48 | 1,430.23 |

(0,1,0) (1,1,1)^{27} | 28,929.2 | 90.40 | 1,552.55 |

(0,1,0) (1,1,2)^{27} | 28,421.2 | 89.09 | 1,548.33 |

(0,1,0) (1,1,3)^{27} | 30,131.6 | 94.75 | 1,570.84 |

(0,1,0) (0,1,1)^{27} | 29,507.5 | 91.92 | 1,557.50 |

(0,1,0) (0,1,2)^{27} | 29,489.9 | 92.16 | 1,559.29 |

(0,1,0) (0,1,3)^{27} | 30,641.9 | 96.06 | 1,574.74 |

(0,1,0) (1,1,0)^{27} | 43,813.2 | 136.49 | 1,696.24 |

(1,1,1) (1,1,1)^{27} | 18,819.4 | 59.18 | 1,405.63 |

(1,1,1) (1,1,2)^{27} | 18,509.4 | 58.39 | 1,401.80 |

(1,1,1) (1,1,3)^{27} | 19,347.7 | 61.23 | 1,419.35 |

(1,1,1) (0,1,0)^{27} | 33,912.5 | 105.98 | 1,608.33 |

(1,1,1) (0,1,1)^{27} | 18,952.4 | 59.41 | 1,404.10 |

(1,1,1) (0,1,2)^{27} | 18,957.7 | 59.62 | 1,408.20 |

(1,1,1) (0,1,3)^{27} | 20,010.4 | 63.12 | 1,429.17 |

(1,1,1) (1,1,0)^{27} | 27,472.3 | 86.12 | 1,536.41 |

Model (SARIMA) . | RSSE . | MSE . | AIC . |
---|---|---|---|

(0,1,1) (1,1,0)^{27} | 27,605.2 | 86.27 | 1,536.11 |

(0,1,1) (1,1,1)^{27} | 18,914.2 | 59.29 | 1,405.40 |

(0,1,1) (1,1,2)^{27} | 18,641 | 58.62 | 1,402.29 |

(0,1,1) (1,1,3)^{27} | 19,518.4 | 61.57 | 1,420.43 |

(0,1,1) (0,1,0)^{27} | 33,912.8 | 105.65 | 1,606.34 |

(0,1,1) (0,1,1)^{27} | 1,9019.3 | 59.44 | 1,405.34 |

(0,1,1) (0,1,2)^{27} | 19,004.4 | 59.58 | 1,407.07 |

(0,1,1) (0,1,3)^{27} | 20,185.4 | 63.48 | 1,430.23 |

(0,1,0) (1,1,1)^{27} | 28,929.2 | 90.40 | 1,552.55 |

(0,1,0) (1,1,2)^{27} | 28,421.2 | 89.09 | 1,548.33 |

(0,1,0) (1,1,3)^{27} | 30,131.6 | 94.75 | 1,570.84 |

(0,1,0) (0,1,1)^{27} | 29,507.5 | 91.92 | 1,557.50 |

(0,1,0) (0,1,2)^{27} | 29,489.9 | 92.16 | 1,559.29 |

(0,1,0) (0,1,3)^{27} | 30,641.9 | 96.06 | 1,574.74 |

(0,1,0) (1,1,0)^{27} | 43,813.2 | 136.49 | 1,696.24 |

(1,1,1) (1,1,1)^{27} | 18,819.4 | 59.18 | 1,405.63 |

(1,1,1) (1,1,2)^{27} | 18,509.4 | 58.39 | 1,401.80 |

(1,1,1) (1,1,3)^{27} | 19,347.7 | 61.23 | 1,419.35 |

(1,1,1) (0,1,0)^{27} | 33,912.5 | 105.98 | 1,608.33 |

(1,1,1) (0,1,1)^{27} | 18,952.4 | 59.41 | 1,404.10 |

(1,1,1) (0,1,2)^{27} | 18,957.7 | 59.62 | 1,408.20 |

(1,1,1) (0,1,3)^{27} | 20,010.4 | 63.12 | 1,429.17 |

(1,1,1) (1,1,0)^{27} | 27,472.3 | 86.12 | 1,536.41 |

Season . | Estimates . | P-value
. | Estimated model . |
---|---|---|---|

Snow Accumulation (1,1,1) (0,1,3)^{19} | 0 | ||

0 | |||

0 | |||

0.009 | |||

0.159 | |||

0.358 | |||

Snow Depletion (1,1,1) (1,1,2)^{27} | 0.18 | ||

0 | |||

0 | |||

0.73 | |||

0 | |||

0.62 |

Season . | Estimates . | P-value
. | Estimated model . |
---|---|---|---|

Snow Accumulation (1,1,1) (0,1,3)^{19} | 0 | ||

0 | |||

0 | |||

0.009 | |||

0.159 | |||

0.358 | |||

Snow Depletion (1,1,1) (1,1,2)^{27} | 0.18 | ||

0 | |||

0 | |||

0.73 | |||

0 | |||

0.62 |

### Diagnostic checking

### Model calibration and forecasting

^{2}) achieved as 0.83 in snow accumulation (Figure 11) while a value of R

^{2}was 0.89 for snow depletion (Figure 12). The below comparison increases confidence in the SARIMA models to represent the 8-day average snow accumulation and snow depletion data at the Beas river and both can be used for forecasting future snow accumulation and snow depletion.

The regional heterogeneity of the underlying processes and inputs causes spatial variability in snow accumulation and melt. Due to this variability, the snow cover melts in patches and has a non-uniform geographical distribution of snow water equivalent (SWE) (Luce & Tarboton 2004). From the distribution function of SWE at peak accumulation, Luce *et al.* (1999) demonstrated how to infer the depletion curve's form. Liston (1999) developed the interrelationship in between the distribution of snow, melt, and snow cover depletion. The study's findings point to the development of the currently covered area and a sub-model that effectively calculates snowmelt. Alternatively, the information on the snow-covered region might be replaced by knowledge of the sub-grid SWE distribution. Shamir & Georgakakos (2007) estimate the SWE and SCA to analyze the SDC for the American river basin using empirical methods and show a significant correlation in between SCA and SWE. However, in the Sutlej river basin, Shukla *et al.* (2017) estimated the variability in the SCA is about 44–56% during 2001–2014. On the other hand, Kour *et al.* (2016) analyzed the effect of slope classes in the variability of SCA of the Chenab basin. Also, the expanse at 3,600–4,300 m elevation, and in the depletion period, 4,300–5,000 m of elevation found a maximum rate of change in the SCA% per 100 m rise in elevation.

The results clearly show that the snow accumulation in the Beas catchment starts in the month of October with a percentage of ∼45%, and rise up to a value of ∼70%, in the month of December, and this magnitude will reach up to a value of ∼93% in the month of January and February because of the western disturbance in all three forecasted years 2016, 2017, and 2018, respectively. Similarly, snow depletion starts in the month of April with a value of ∼80% to a value of ∼23% in the month of September because of a rise in temperature and ablation.

## CONCLUSIONS

In the present study, the SARIMA model is applied to the time series of 8-day average now data of MODIS satellite for snow-dominant area of the Beas catchment (2003–2018). The Box–Jenkins methodology of model forecasting has been applied after splitting the data into two main seasons: snow accumulation (October–February) and snow depletion (March–September). The models predict the next 8 days of weekly snow accumulation and snow depletion for 3 years, based on the past 13 years of data. The performance of the resulting models as (1,1,1) (0,1,3)^{19} for snow accumulation and (1,1,1) (1,1,2)^{27} for snow depletion was also validated by comparing the forecasted values with observed data (2016–2018). The forecasted data showed good agreement with the actual data in terms of coefficient of determination as 0.83 in snow accumulation and a value of 0.89 in snow depletion. Moreover, the results achieved in the study gave increasing confidence in using proposed models to forecast the weekly snow accumulation and snow depletion for at least the next 3 years (2019–2021) for the estimation of hydraulic events such as runoff, flood forecasting, and hydropower potential in the Beas region.

## ACKNOWLEDGEMENTS

The first and the third author would like to give special thanks to Dr D K Singh, presently working as a scientist at the Indian Institute of Technology, Bombay in the field of remote sensing, snow avalanche climatology, and snow cover properties, for the MODIS data preparation in the present study. The second author would like to thank the Department of Civil Engineering, Indian Institute of Technology, Ropar, and the fourth author would like to thank the Department of Civil Engineering, G.B. Pant Institute of Engineering and Technology, Pauri (Garhwal) for providing facilities.

## FUNDING INFORMATION

The author(s) received no specific funding for this work.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

*An Introductory Study on Time Series Modeling and Forecasting*. LAP Lambert Academic Publishing, Germany.

*Journal of the American Statistical Association*

**71**(356), 791–799.

*Time Series Analysis: Forecasting and Control*. John Wiley & Sons, Hoboken, NJ.

*Journal of Water Supply: Research and Technology-Aqua*

**72**(2), 139–159.

*An Introduction to Machine Learning.*Springer, Germany.

*Snowmelt Runoff Analysis and Modeling for the Upper Cache La Poudre River Basin, Colorado*

*MsS Thesis*

*International Journal of Applied Earth Observation and Geoinformation*

**63**, 234–243.