Listen

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 R2 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.

Listen

  • 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.

ACF and PACF plots, AIC, MODIS time series data, MSE, RSSE.4, SARIMA model
Listen

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.

Listen

Study area and data used

Listen
The Beas river emerging from the Beas Kund in the Himalayas is one of the major rivers of north-western India, and stretches up to 470 km to finally meet the river Sutlej in Punjab. The Beas river has a length of 116 kilometers up to the Pandoh Dam has been selected in the present study. The region of interest lies between 31°70′ N – 32°42′ N and 77°5′ E – 77°73′ E. The Parbati and Sainj Khad (glacier-fed) are the main tributaries of Beas that meets at Bhuntar and Larji. Some of the other tributaries like Tirthan, Bakhli khad, and Luni meet the Beas river near Pandoh dam. The drainage network map of the region of interest is shown in Figure 1. The MODIS 8-day average cloud-Free snow cover composite products of MOD10A2 (Terra) and MYD10A2 (Aqua) consisting of a spatial resolution of 500-m gridded geo-referenced cells have been downloaded and then converted to useful product (http://ladsweb.nascom.nasa.gov/ website). The study reported that the MODIS satellite data have higher accuracy and temporal resolution than NOAA and IRS (Jain 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).
Figure 1
Location map of the Beas basin with ASTER DEM showing the lowest and highest points.
View largeDownload slide

Location map of the Beas basin with ASTER DEM showing the lowest and highest points.

Figure 1
Location map of the Beas basin with ASTER DEM showing the lowest and highest points.
View largeDownload slide

Location map of the Beas basin with ASTER DEM showing the lowest and highest points.

Close modal

SARIMA methodology

Listen
In 1976, Box and Jenkins gave a methodology to forecast the best-fit model from the historical time series data. The methodology involved four steps: model identification, parameter estimation, diagnostic checking for appropriateness, and application of the model for forecasting. The important analytical tools for diagnostic checking are Autocorrelation Function (ACF) to check the linear dependence between the observations separated by lag 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)S 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 (p, d, q) (P, D, Q)S and presented in the form of an equation as Equation (1) (Hazarika et al. 2020).
(1)
where represents Gaussian white noise at period t, Xt as time series data at period t, and is the constant. ∇d is the non-seasonal difference component as (1 − B)d, is the seasonal difference component as (1 − BS)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):
(2)
(3)
Similarly, the seasonal autoregressive operator and moving average operator are represented as and ⁠, respectively, in the form of polynomials as order P and Q (Equations (4) and (5):
(4)
(5)

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

Listen
The model selection involved three main steps model identification, parameter estimation, and diagnostic checking in performing the statistical analysis (Eze 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) as
(6)
where n refers to the size of the series, and p, q, P, and Q 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 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.
Listen

Model identification

Listen
The MODIS 8-day cloud-free satellite data of snow from duration (2003–2018) initially has been divided into two seasons: snow accumulation (October–February) and snow depletion (March–September). The trend analysis of original data for snow accumulation shows an increasing trend while snow depletion data has little increasing trend as shown in Figure 2. Moreover, the ACF and PACF were plotted (Figures 3 and 4).
Figure 2
Trend analysis plot for (a) snow accumulation and (b) snow depletion (2003–2015).
View largeDownload slide

Trend analysis plot for (a) snow accumulation and (b) snow depletion (2003–2015).

Figure 2
Trend analysis plot for (a) snow accumulation and (b) snow depletion (2003–2015).
View largeDownload slide

Trend analysis plot for (a) snow accumulation and (b) snow depletion (2003–2015).

Close modal
Figure 3
ACF and PACF plots of snow accumulation (2003–2015).
View largeDownload slide

ACF and PACF plots of snow accumulation (2003–2015).

Figure 3
ACF and PACF plots of snow accumulation (2003–2015).
View largeDownload slide

ACF and PACF plots of snow accumulation (2003–2015).

Close modal
Figure 4
ACF and PACF plots of snow depletion (2003–2015).
View largeDownload slide

ACF and PACF plots of snow depletion (2003–2015).

Figure 4
ACF and PACF plots of snow depletion (2003–2015).
View largeDownload slide

ACF and PACF plots of snow depletion (2003–2015).

Close modal

Stationarity and seasonality test

Listen
The first step of model building is to check whether the available time series to be forecasted has stationarity and seasonality. A check for stationarity means to check any kind of trend over time means the time series data should have a constant mean while a seasonal variation is a regularly repeating pattern over a given duration of the period. The ACF and PACF plot of weekly snow accumulation indicate seasonality as a slowly declining sinusoidal wave and significant peaks at seasonal lag values 19, 38, and 57 with alternating signs (Figure 3). Similarly, the ACF and PACF plots of weekly snow depletion show similar slow decline sinusoidal waves and strong significant peaks at seasonal lag values 27, 54, and 81 with alternating signs (Figure 4). To obtain a stationary time series, i.e., a series having constant mean and no trend over time, first differencing of original data using Augmented Dickey-Fuller (ADF) test along with seasonal period of 19 in case of snow accumulation and a seasoning period of 27 in snow depletion were adopted. After the transformation in the original data, the trend analysis, and ACF and PACF plots of snow accumulation and depletion seem to be stationarity (Figures 5 and 6).
Figure 5
Trend analysis of snow accumulation after taking differencing of lag 1 along with ACF and PACF plots (2003–2015).
View largeDownload slide

Trend analysis of snow accumulation after taking differencing of lag 1 along with ACF and PACF plots (2003–2015).

Figure 5
Trend analysis of snow accumulation after taking differencing of lag 1 along with ACF and PACF plots (2003–2015).
View largeDownload slide

Trend analysis of snow accumulation after taking differencing of lag 1 along with ACF and PACF plots (2003–2015).

Close modal
Figure 6
Trend analysis of snow depletion after taking differencing of lag 1 along with ACF and PACF plots (2003–2015).
View largeDownload slide

Trend analysis of snow depletion after taking differencing of lag 1 along with ACF and PACF plots (2003–2015).

Figure 6
Trend analysis of snow depletion after taking differencing of lag 1 along with ACF and PACF plots (2003–2015).
View largeDownload slide

Trend analysis of snow depletion after taking differencing of lag 1 along with ACF and PACF plots (2003–2015).

Close modal

Model estimation

Listen

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)

Listen
From plotting in Figure 5 (non-seasonal), the ACF cuts off at lag 1 while PACF cuts off at lag 2, so the possibility of MA (1) and AR (2), and for seasonal variation, the ACF shows significant spikes at lag 19, 38, and 57 while PACF plots show the significant seasonal lag at 19, so the possibility of SMA (3) and SAR (1). So, a rough idea of best-fit model estimation in terms of (⁠ has been acquired nearby the above observation by considering the combined behavior as (2,1,1) (1,1,3)19 and put into Equation (1) and the final model has been proposed as Equation (7):
(7)

Seasonal and non-seasonal variation (snow depletion)

Listen
Similarly, from plotting in Figure 6 (non-seasonal), both ACF and PACF cuts off at lag 1 so the possibility of MA (1) and AR (1), and for seasonal variation, the ACF shows significant spikes at lag 27, 54, and 81 while PACF shows the significant seasonal lag at 27, so the possibility of SMA (3) and SAR (1). So, a rough idea of the best-fit model in terms of (⁠ has been acquiring nearby the above observation by considering the combined behavior as (1,1,1) (1,1,3)27 and the same has been put into Equation (1) and the final model has been drawn as Equation (8):
(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.

Table 1

Postulate models with evaluation criterion (snow accumulation)

Model (SARIMA)RSSEMSEAIC
(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)RSSEMSEAIC
(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 
View Large
Table 2

Postulate models with evaluation criterion (snow depletion)

Model (SARIMA)RSSEMSEAIC
(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)RSSEMSEAIC
(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 
View Large
Table 3

Estimates of SARIMA models for snow accumulation and snow depletion

SeasonEstimatesP-valueEstimated 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 
SeasonEstimatesP-valueEstimated 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 
View Large

Diagnostic checking

Listen
The fitted models were examined for adequacy by testing the normal probability and histogram plots of residuals. Moreover, the ACF and PACF plots of residuals at different lags should be within 95% confidence limits. A model with no significant residuals should be considered adequate. The normal probability plot of the residuals for both snow accumulation and snow depletion clearly shows that the residuals were on a straight line (top left of Figure 7 for snow accumulation and top left of Figure 8 for snow depletion). The histogram plots of the residuals were also shown good normality of the residuals for both (bottom left of Figure 7 for snow accumulation and bottom left of Figure 8 for snow depletion). The third measure as the plots of residuals against fitted values also clearly showed that the variance of the error terms is constant with a mean of zero (top right of Figure 7 for snow accumulation and top right of Figure 8 for snow depletion). The fourth/last measure as the plot of residuals against the fitted order of the data does not follow any symmetric pattern with the run order value (bottom right of Figure 7 for snow accumulation and bottom right of Figure 8 for snow depletion). The random behavior of residuals indicates the model is the best fit for both seasons. The plots of the residuals are also within the acceptable limits with a mean of 0.03 and a standard deviation of 10.8 for snow accumulation (Figure 9) and a mean of 0.24 and a standard deviation of 7.58 for snow depletion (Figure 10).
Figure 7
Residual plots of the SARIMA model for snow accumulation (1,1,1) (0,1,3)19.
View largeDownload slide

Residual plots of the SARIMA model for snow accumulation (1,1,1) (0,1,3)19.

Figure 7
Residual plots of the SARIMA model for snow accumulation (1,1,1) (0,1,3)19.
View largeDownload slide

Residual plots of the SARIMA model for snow accumulation (1,1,1) (0,1,3)19.

Close modal
Figure 8
Residual plots of the SARIMA model for snow depletion (1,1,1) (1,1,2)27.
View largeDownload slide

Residual plots of the SARIMA model for snow depletion (1,1,1) (1,1,2)27.

Figure 8
Residual plots of the SARIMA model for snow depletion (1,1,1) (1,1,2)27.
View largeDownload slide

Residual plots of the SARIMA model for snow depletion (1,1,1) (1,1,2)27.

Close modal
Figure 9
ACF and PACF plots of residuals for snow accumulation (1,1,1) (0,1,3)19.
View largeDownload slide

ACF and PACF plots of residuals for snow accumulation (1,1,1) (0,1,3)19.

Figure 9
ACF and PACF plots of residuals for snow accumulation (1,1,1) (0,1,3)19.
View largeDownload slide

ACF and PACF plots of residuals for snow accumulation (1,1,1) (0,1,3)19.

Close modal
Figure 10
ACF and PACF plots of residuals for snow depletion (1,1,1) (1,1,2).
View largeDownload slide

ACF and PACF plots of residuals for snow depletion (1,1,1) (1,1,2).

Figure 10
ACF and PACF plots of residuals for snow depletion (1,1,1) (1,1,2).
View largeDownload slide

ACF and PACF plots of residuals for snow depletion (1,1,1) (1,1,2).

Close modal

Model calibration and forecasting

Listen
The performance of the models was evaluated by forecasting the data for the years (2016–2018) and the same was validated with the actual data. Both the forecasted and the actual values of the 8-day average snow accumulation and depletion data (2016–2018) were fitted well as the coefficient of determination (R2) achieved as 0.83 in snow accumulation (Figure 11) while a value of R2 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.
Figure 11
Comparison of actual and forecasted snow accumulation (2016–2018).
View largeDownload slide

Comparison of actual and forecasted snow accumulation (2016–2018).

Figure 11
Comparison of actual and forecasted snow accumulation (2016–2018).
View largeDownload slide

Comparison of actual and forecasted snow accumulation (2016–2018).

Close modal
Figure 12
Comparison of actual and forecasted snow depletion (2016–2018).
View largeDownload slide

Comparison of actual and forecasted snow depletion (2016–2018).

Figure 12
Comparison of actual and forecasted snow depletion (2016–2018).
View largeDownload slide

Comparison of actual and forecasted snow depletion (2016–2018).

Close modal

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.

Listen

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.

Listen

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.

Listen

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

Listen

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

Listen

The authors declare there is no conflict.

Adams
S. O.
,
Mustapha
B.
 &
Alumbugu
A. I.
2019
Seasonal Autoregressive Integrated Moving Average (SARIMA) model for the analysis of the frequency of monthly rainfall in osun state, Nigeria
.
Physical Science International Journal
1
–
9
.
https://doi.org/10.9734/psij/ 2019/v22i43013
.
Google Scholar
 
Adhikari, R. & Agrawal, R. K. 2013 An Introductory Study on Time Series Modeling and Forecasting. LAP Lambert Academic Publishing, Germany.
Ahluwalia
R. S.
,
Rai
S. P.
,
Jain
S.
,
Bhisham
K.
 &
Dobhal
D. P.
2013
Assessment of snowmelt runoff modeling and isotope analysis a case study from the western Himalayas, India
.
Annals of Glaciology
54
(
62
),
299
–
304
.
Google Scholar
Crossref
Search ADS
 
Akaike
H.
1974
A new look at the statistical model identification
.
IEEE Transaction on Automatic Control
19
,
716
–
723
.
http://dx/doi.org/10.1109/TAC.1974.1100705
.
Google Scholar
Crossref
Search ADS
 
Astuti
S. W.
 &
Jamaludin
2018
Forecasting Surabaya Jakarta train passengers with SARIMA model
.
IOP Conference Series: Materials Science and Engineering
407
,
012105
.
https://iopscience.iop.org/article/10.1088/1757-899X/407/1/012105/meta
.
Google Scholar
Crossref
Search ADS
 
Box, G. E. 1976 Science and statistics. Journal of the American Statistical Association 71 (356), 791–799.
Box
G. E. P.
 &
Jenkins
G. M.
1976
Time series analysis forecasting and control – rev. ed Oakland, California, Holden-Day, 37 (2), 238-242
.
Box, G. E., Jenkins, G. M., Reinsel, G. C. & Ljung, G. M. 2015 Time Series Analysis: Forecasting and Control. John Wiley & Sons, Hoboken, NJ.
Chen
P.
,
Niu
A.
,
Liu
D.
,
Jiang
W.
 &
Ma
B.
2018
Time series forecasting of temperatures using SARIMA: an example from Nanjing
. In
IOP Conference Series: Materials Science and Engineering
.
IOP Publishing
, p.
052024
.
Google Scholar
Crossref
Search ADS
 
Colbeck
S. C.
1978
The physical aspects of water flow through the snow
.
Advances in Hydroscience
11
,
165
–
206
.
Google Scholar
Crossref
Search ADS
 
Dhanju
M. S.
1983
Studies of Himalayan snow cover area from satellites. Hydrological Applications of Remote Sensing and Remote Data Transmission, Proceedings of the Hamburg Symposium, IAHS (145), 401–409
.
Di Marco
N.
,
Avesani
D.
,
Righetti
M.
,
Zaramella
M.
,
Majone
B.
 &
Borga
M.
2021
Reducing hydrological modeling uncertainty by using MODIS snow cover data and a topography-based distribution function snowmelt model
.
Journal of Hydrology
599
,
126020
.
Google Scholar
Crossref
Search ADS
 
Dozier
J.
 &
Painter
H. T.
2004
Multispectral and hyperspectral remote sensing of alpine snow properties
.
Annual Reviews of Earth and Planetary Science
32
,
465
–
494
.
Google Scholar
Crossref
Search ADS
 
Eni
D.
 &
Adeyeye
F. J.
2015
Seasonal ARIMA modelling and forecasting of rainfall in Warri Town, Nigeria
.
Journal of Geoscience and Environment Protection
3
,
91
–
98
.
Google Scholar
Crossref
Search ADS
 
Eze
O. N.
,
Asogwa
A.
,
Obetta
K.
 &
Ojide
C. O.
2020
Time-series analysis of federal budgetary allocations to the education sector in Nigeria (1970–2018)
.
American Journal of Applied Mathematics and Statistics
8
(
1
),
1
–
8
.
Google Scholar
 
Farsi
M.
,
Hosahalli
D.
,
Manjunatha
B. R.
,
Gad
I.
,
Atlam
E. S.
,
Ahmed
A.
,
Elmarhomy
G.
,
Mahmoud
E.
 &
Ghoneim
O. A.
2020
Parallel genetic algorithms for optimizing the SARIMA model for better forecasting of the NCDC weather data
.
Alexandria Engineering Journal
60
,
1299
–
1316
.
(2021)
.
Google Scholar
Crossref
Search ADS
 
Hall
D. K.
,
Riggs
G. A.
,
Salomonson
V. V.
,
DiGirolamo
N. E.
 &
Bayr
K. A.
2002
MODIS snow-cover products
.
Remote Sensing of Environment
83
,
181
–
194
.
Google Scholar
Crossref
Search ADS
 
Hall
D. K.
,
Riggs
G. A.
,
Foster
J. L.
 &
Kumar
S. V.
2010
Development and evaluation of a cloud-gap-filled MODIS daily snow-cover product
.
Remote Sensing of Environment
114
,
496
–
503
.
Google Scholar
Crossref
Search ADS
 
Haritashya
U. K.
,
Singh
P.
,
Kumar
N.
 &
Gupta
R. P.
2006
Suspended sediment from the Gangotri Glacier: quantification, variability, and associations with discharge and air temperature
.
Journal of Hydrology
321
(
1–4
),
116
–
130
.
Google Scholar
 
Hazarika
J.
,
Goswami
K.
,
Bharali
S.
,
Bora
S. L.
 &
Patowary
A. N.
2020
Prediction of monthly mean temperature by SARIMA model
.
International Journal of Advanced Science and Technology
29
,
3937
–
3949
.
Google Scholar
 
Immerzeel
W. W.
,
Droogers
P. d. J. S.
 &
Bierkens
M. F. P.
2009
Largescale monitoring of snow covers and runoff simulation in Himalayan river basins using remote sensing
.
Remote Sensing of Environment
113
,
40
–
49
.
Google Scholar
Crossref
Search ADS
 
Jabal
Z. K.
,
Khayyun
T. S.
 &
Alwan
I. A.
2022
Impact of climate change on crop productivity using MODIS-NDVI time series
.
Civil Engineering Journal
8
(
6
).
http://dx.doi.org/10.28991/CEJ-2022-08-06-04
.
Google Scholar
 
Jain
S. K.
,
Goswami
A.
 &
Saraf
A. K.
2010
Snowmelt runoff modeling in a Himalayan basin with the aid of satellite data
.
International Journal of Remote Sensing
31
,
6603
–
6618
.
Google Scholar
Crossref
Search ADS
 
Jain
S. K.
,
Goswami
A.
 &
Sarab
A.
2018
Accuracy assessment of MODIS, NOAA, and IRS data in snow cover mapping under Himalayan conditions
.
International Journal of Remote Sensing
29
(
20
),
5863
–
5878
.
Google Scholar
Crossref
Search ADS
 
Khandelwal
I.
,
Adhikari
R.
 &
Verma
G.
2015
Time series forecasting using hybrid ARIMA and ANN models based on dwt decomposition
.
Procedia Computer Science
48
,
173
–
179
.
Google Scholar
Crossref
Search ADS
 
Kour
R.
,
Patel
N.
 &
Krishna
A. P.
2016
Effects of terrain attributes on snow-cover dynamics in parts of Chenab basin, western Himalayas
.
Hydrological Sciences Journal
61
(
10
),
1861
–
1876
.
Google Scholar
 
Kulkarni
A. V.
,
Rathore
B. P.
,
Singh
S. K.
 &
Ajai
2010
Distribution of seasonal snow cover in central and western Himalayas
.
Annals of Glaciology
51
(
54
),
123
–
128
.
Google Scholar
Crossref
Search ADS
 
Kumar
V.
,
Singh
P.
 &
Jain
S. K.
2005
Rainfall trends over Himachal Pradesh, Western Himalaya, India
. In
Proceedings Conference on Development of Hydro Power Projects – A Prospective Challenge, Shimla
,
April 20–22
.
Central Board of Irrigation and Power
,
New Delhi
,
India
.
Google Scholar
 
Luce
C. H.
 &
Tarboton
D. G.
2004
The application of depletion curves for parameterization of subgrid variability of snow
.
Hydrological Processes
18
(
8
),
1409
–
1422
.
doi:10.1002/hyp.1420
.
Google Scholar
Crossref
Search ADS
 
Luce
C. H.
,
Tarboton
D. G.
 &
Cooley
K. R.
1999
Sub-grid parameterization of snow distribution for an energy and mass balance snow cover model
.
Hydrological Processes
13
(
12-13
),
1921
–
1933
.
Google Scholar
Crossref
Search ADS
 
Lutz
S. R.
,
Mallucci
S.
,
Diamantini
E.
,
Majone
B.
,
Bellin
A.
 &
Merz
R.
2016
Hydro climatic and water quality trends across three Mediterranean river basins
.
Science of the Total Environment
571
,
1392
–
1406
.
Google Scholar
Crossref
Search ADS
PubMed
 
Mann
R.
 &
Gupta
A.
2022
Temporal trends of rainfall and temperature over two sub divisions of Western Ghats
.
High Tech and Innovative Journal
3
.
http://dx.doi.org/10.28991/HIJ- SP2022-03-03
.
Google Scholar
 
Matson
M.
,
Roeplewski
C. F.
 &
Varnadore
M. S.
1986
An Atlas of Satellite-Derived Northern Hemisphere Snow Cover Frequency
.
National Weather Service
,
Washington, DC
, p.
75
.
Google Scholar
 
Prasad
V. H.
 &
Roy
P. S.
2005
Estimation of snowmelt runoff in Beas Basin
.
India Geocarto International
20
(
2
),
41
–
47
.
Google Scholar
Crossref
Search ADS
 
Raicharoen
T.
,
Lursinsap
C.
 &
Sanguanbhokai
P.
2003
Application of critical support vector machine to time series prediction
. In
International Symposium on Circuits and Systems
.
IEEE
, Vol.
5
, pp.
V-741
–
V-744
.
Google Scholar
Crossref
Search ADS
 
Rautela, K. S., Kumar, D., Gandhi, B. G. R., Kumar, A. & Dubey, A. K. 2023 Long-term hydrological simulation for the estimation of snowmelt contribution of Alaknanda River Basin, Uttarakhand using SWAT. Journal of Water Supply: Research and Technology-Aqua 72 (2), 139–159.
Rebala, G., Ravi, A. & Churiwala, S. 2019 An Introduction to Machine Learning.Springer, Germany.
Richer
E. E.
2009
Snowmelt Runoff Analysis and Modeling for the Upper Cache La Poudre River Basin, Colorado
.
MsS Thesis
,
Colorado State University-Fort Collins
,
USA
.
Google Scholar
 
Rizul
N. S.
,
Umarie
I.
,
Munandar
K.
 &
Wardoyo
A. E.
2023
Calibration and validation of CN values for watershed hydrological response
.
Civil Engineering Journal
9
(
1
).
http://dx.doi.org/10.28991/CEJ-2023-09-01-06
.
Google Scholar
 
Shamir
E.
 &
Georgakakos
K. P.
2007
Estimating snow depletion curves for American river basins using distributed snow modeling
.
Journal of Hydrology
334
(
1–2
),
162
–
173
.
Google Scholar
 
Shukla
S.
,
Kansal
M. L.
 &
Jain
S. K.
2017
Snow cover area variability assessment in the upper part of the Satluj river basin in India
.
Geocarto International
32
(
11
),
1285
–
1306
.
Google Scholar
Crossref
Search ADS
 
Sichangi
A. W.
,
Wang
L.
,
Yang
K.
,
Chen
D.
,
Wang
Z.
,
Li
X.
,
Zhou
J.
,
Liu
W.
 &
Kuria
D.
2016
Estimating continental river basin discharges using multiple remote sensing data sets
.
Remote Sensing of Environment
179
,
36
–
53
.
Google Scholar
Crossref
Search ADS
 
Singh
P.
,
Ramasastri
K. S.
 &
Kumar
N.
1995
Topography influences precipitation distribution in different ranges of the western Himalayas
.
Nordic Hydrology
26
(
4–5
),
259
–
289
.
Google Scholar
 
Singh
P.
,
Jain
S. K.
 &
Kumar
N.
1997
Estimation of snow and glacier contribution to the Chenab River, Western Himalaya
.
Mountain Research and Development
17
(
1
),
49
–
56
.
Google Scholar
Crossref
Search ADS
 
Steele, C., Dialesandro, J., James, D., Elias, E., Rango, A. & Bleiweiss, M. 2017 Evaluating MODIS snow products for modelling snowmelt runoff: Case study of the Rio Grande headwaters. International Journal of Applied Earth Observation and Geoinformation 63, 234–243.
Zhou
X. B.
,
Xie
H. J.
 &
Hendrickx
J. M. H.
2005
Statistical evaluation of remotely sensed snow-cover products with constraints from streamflow and SNOTEL measurements
.
Remote Sensing of Environment
94
,
214
–
231
.
https://doi.org/10.1016/j.rse.2004.10.007
.
Google Scholar
Crossref
Search ADS
 
© 2023 The Authors
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).