Using the MODIS Sensor for Snow Cover Modeling and the Assessment of Drought Effects on Snow Cover in a Mountainous Area

Abstract

Snow is one of the essential factors in hydrology, freshwater resources, irrigation, travel, pastimes, floods, avalanches, and vegetation. In this study, the snow cover of the northern and southern slopes of Alborz Mountains in Iran was investigated by considering two issues: (1) Estimating the snow cover area and the (2) effects of droughts on snow cover. The snow cover data were monitored by images obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor. The meteorological data (including the precipitation, minimum and maximum temperature, global solar radiation, relative humidity, and wind velocity) were prepared by a combination of National Centers for Environmental Prediction-Climate Forecast System Reanalysis (NCEP-CFSR) points and meteorological stations. The data scale was monthly and belonged to the 2000–2014 period. In the first part of the study, snow cover estimation was conducted by Multiple Linear Regression (MLR), Least Square Support Vector Machine (LSSVM), Group Method of Data Handling (GMDH), Multilayer Perceptron (MLP), and MLP with Grey Wolf Optimization (MLP-GWO) models. The most accurate estimations were produced by the MLP-GWO and GMDH models. The models produced better snow cover estimations for the northern slope compared to the southern slope. The GWO improved the MLP’s accuracy by 10.7%. In the second part, seven drought indices, including the Palmer Drought Severity Index (PDSI), Bahlme–Mooley Drought Index (BMDI), Standardized Precipitation Index (SPI), Multivariate Standardized Precipitation Index (MSPI), Modified Standardized Precipitation Index (SPImod), Joint Deficit Index (JDI), and Standardized Precipitation-Evapotranspiration Index (SPEI) were calculated for both slopes. The results showed that the effects of a drought event on the snow cover area would remain up to 5 (or 6) months in the region. The highest impact of drought appears after two months in the snow cover area, and the drought index most related to snow cover variations is the 2–month time window of SPI (SPI₂). The results of both subjects were promising and the methods can be examined in other snowy areas of the world.

1. Introduction

Snow cover in many regions around the world has a direct impact on human life, including engineering, irrigation, travel, recreation, and hydrology. Additionally, floods and avalanches are important natural hazards, which are strongly associated with snow [1]. Another major aspect of snow is its interaction with vegetation. The melting of this water source contributes to vegetation survival and fertility, and has a profound effect on the phenomenon of soil frost and plant chilling injury [2,3,4].

The snow cover area is one of the most important parameters of the hydrological and climate cycles. Field measurements can provide clues for observable changes in much of our knowledge on the cryosphere. Over the last 20 years, due to the lack of appropriate accessibility, the influence of various topographic and physiographic features on the estimation of cryosphere hydrological parameters, and the insufficiency of meteorological stations at a high altitude have made the necessity of using indirect methods more prominent at regional and global scales. Moreover, the snowmelt process has a high energy requirement that results in atmospheric warming. Therefore, any change in snow cover has a direct impact on the variability of the water and energy cycle [3,4].

Snow cover can be estimated using modeling, measurement stations, and remote sensing applications [5]. Various satellite data have been used to identify the extent of snow cover and, in recent years, the Moderate Resolution Imaging Spectroradiometer (MODIS) has been one of the most commonly used approaches for snow cover monitoring [6,7,8,9,10]. The disadvantage of this method is that the MODIS sensor was established in 2000, and before that, there were no images available for snow cover monitoring. An alternative method is the use of estimation models, which can use meteorological datasets (which are available for periods before 2000) for snow cover estimation.

Statistical models and artificial intelligence-based methods have been widely used in hydrological modeling studies [6]. However, there are relatively few studies on the application of these methods for modeling snow parameters [7]. Tabari et al. (2010) compared snow depth and snow water equivalent (SWE) estimations of multivariate linear regression (MLR), discriminant function analysis, ordinary kriging, ordinary kriging-multivariate linear regression, the artificial neural network (ANN), and the neural network-genetics algorithm (NNGA). The results indicated that the NNGA, ANN, and MLR methods were able to predict SWE at the desirable level of accuracy [8]. Lee et al. (2019) used ANN to estimate the maximum daily fresh snow accumulation (MDFSA) in South Korea. The results showed that the ANN model has an acceptable performance [9]. Khadka et al. (2014) studied the impact of climate change on snow and glacier melting and its subsequent runoff in the Tamakoshi Basin located in the Hindu Kush Himalaya region. They used HADCM3, CGCM3, SDSM, and SERES to predict snow cover changes for the future, based on temperature and precipitation data and MLR [10]. From the literature review, it could be observed that there is no published study related to estimating the snow cover area using meteorological variables. Therefore, in the current study, machine learning and regression models were used for this issue for the first time.

Furthermore, the surface features of snowfall from a hydroclimatic perspective play a major role in various areas, such as flood control, watershed management, the water supply, soil erosion, and droughts. A drought affects almost all of the hydrological and agricultural sections and their own related index clarifies each of them. For example, a meteorological drought is shown by the Percent of Normal Precipitation Index (PNPI) and Rainfall Anomaly Index (RAI). The Surface Water Supply Index (SWSI) is commonly used for hydrological drought monitoring. Furthermore, the Standardized Precipitation Index (SPI) is an index that can be related to all of the drought types, in its different time windows. For example, SPI’s 1-month time window is related to a meteorological drought and its 3– and 6–month time windows are related to drought conditions from the perspective of soil moisture and surface streamflow, respectively [11]. Moreover, the 9–month to 12–month time windows of SPI show agricultural droughts and those of 12–month to 24–month relate to droughts from the perspective of the groundwater level [11,12,13]. As another investigation of the current study, the researchers will evaluate the effect of drought on snow cover variations and try to find the index most related to the snow cover area.

2. Materials and Methods

2.1. Area of Study

Central Alborz is part of the interconnected Alborz range. It has the highest convexity and transverse thickness and height and extends from the Karaj valley to Semnan fault. The Central Alborz range is divided into the northern and southern slopes and has three walls. At the northern wall, there are limited altitudes in Tehran and Semnan provinces and the rest are situated in Mazandaran province. Secondly, the middle wall forms the northern boundary of Tehran province and is the highest part of the central Alborz range. The highest point of Alborz is Mount Damavand, with a 5671 m height. This massive mountain wall continues to form the Kandovan Mountain and then the Taleghan Mountain in the northwest, to the point where the Alamut river connects to the Taleghan river. To the northeast, this wall extends to Firuzkuh and Savadkuh altitudes to the Firuzkuh river valley (the main branch of Hablehroud), which runs south of its eastern slopes. It continues to the Shahmirzad altitude and Semnan city. Thirdly, the south wall is the third part of the central altitude, with rivers such as Jajrood, Karaj, Hablehrud, and Golrudbar cutting through it and creating gates along the north and south. Topographically, it includes mountains and plains. In this region, snow plays a key role in the hydrological and hydroclimate cycles, and a significant portion of the total annual runoff in the area results from snowmelt. The geographical location of the northern, central, and southern Alborz range and the distribution of the stations under study are shown in Figure 1.

2.2. Data

2.2.1. Snow Cover Data

To monitor snow cover, MODIS model and Digital Elevation Model (DEM) images were used. To that end, the MOD10A1 and MOD10A2 of MODIS data were first obtained from the NASA National Snow and Ice Data Center (NSIDC) database with a spatial resolution of 500 × 500 m in HDF format. Heavy snowfall in central Alborz usually occurs from early November to May and the snow cover persists until June. For this area, images were available for the 15 years from 2000 to 2014. To process the images, first, a preprocessing operation including geometric, radiometric, and atmospheric corrections was applied to them in ENVI 5.3 software. This operation included calibrating images, converting image coordinates to real terrain coordinates (WGS84 UTM 39N), and removing clouds from images. The following provides a description of the MODIS sensor’s usage applied in the current study to monitor the snow cover area.

One of the five sensors on the Terra satellite is the MODIS sensor, which was launched on December 18, 1999. This sensor’s coverage significantly increased with the launch of the Aqua satellite on May 4, 2002. The MODIS sensor has a spatial resolution (250, 500, and 1000 m), 12-bit radiometric resolution, high temporal resolution (1 to 2 days), and 36-band spectral separation of bands 20 to 23, excluding radiation bands and bands 1 to 7, which are reflective. It covers wavelengths from 0.4 microns (visible light) to 14.4 microns (thermal infrared). The MODIS sensor provides snow cover detection maps using the Normalized Difference Snow Index (NDSI) on a macro-scale and provides a rapid application in regional studies [13]. NDSI is a spectral measure of the relative magnitude of the characteristic reflectance difference between the visible and short-wave infrared reflectance of snow in the MODIS sensor, employed to detect changes in snow cover using the Normalized Difference Vegetation Index (NDVI) algorithm [14]. NDSI uses snow spectral reflectance advantages that exist in the visible band with a high reflectance and infrared spectrum with a low reflectance, and NDSI as an algorithm for detecting snow from clouds and snow-covered areas with a set of thresholds applied and calculated pixel to pixel [14]. In the first method, the dual reflectance criteria (i.e., the pixel reflectance values in band 6 > 11% and the pixel reflectance values in band 4 ≥ 10%) are applied if NDSI > 0.4, based on Equation (1). In the next step, three conditional tests are performed to extract the NDSI values, according to Equation (1). Similar to many spectral rationing methods, this index reduces atmospheric effects [15]. Applying the NDSI will result in the formation of pixels with a value of −1 to +1, for which values ranging from −1 to 0 indicate areas where there is no snow and those from 0 to +1 indicate areas where snow has positive coefficients due to its lightness and heaviness (depending on its depth). When the snow depth is higher, the number is closer to one, and for a lower depth of snow, the number tends toward zero. The MODIS snow map algorithm from bands 4 and 6 of this sensor is automatically implemented to extract the Normalized Difference Snow Index (NDSI) and is calculated based on Equation (1) [16].

$$NDSI=\frac{MODI{S}_{Band4}-MODI{S}_{Band6} }{MODI{S}_{Band4}+MODI{S}_{Band6}}=\frac{green-SWIR}{green+SWIR}$$

(1)

In this equation, $green$ is the green coverage, which is obtained from band 4 (0.545–0.665 µm) and has a reflectance ≥ 10%. Band 6 (0.841–0.876 µm) has a reflectance > 11% and relates to Short-Wave Infrared (SWIR).

2.2.2. Meteorological Data

In this study, the combined use of the National Centers for Environmental Prediction-Climate Forecast System Reanalysis (NCEP-CFSR) database and observations due to insufficient meteorological stations located at altitude (from the Iranian Meteorological Organization) is provided. The NCEP-CFSR points and the ground stations are specified in Figure 1. The meteorological variables used include the minimum temperature (T_min), maximum temperature (T_max), global solar radiation (GSR), relative humidity (RH), precipitation (P), and wind velocity (W), and their units are $°\mathrm{C}$, $°\mathrm{C}$, $\frac{\mathrm{M}\mathrm{J}}{{\mathrm{m}}^2·\mathrm{day}}$,$\%$, $\mathrm{mm},$ and $\frac{\mathrm{m}}{\mathrm{s}}$, respectively. These meteorological data, prepared on monthly scales, belong to the period of 1979–2014 (36 years). These datasets were integrated to be used as inputs for snow cover estimation. As there was one data series for snow cover for each of the northern and southern slopes, the meteorological data had to be integrated into one equivalent series for the estimation and investigation. Therefore, the Thiessen polygon method was used to merge the data into one series, for each of the two slopes of Alborz mountain (according to Guhathakurta et al.’s [17] study). The statistical characteristics of

Table 1. Statistical characteristics of the monthly data used in this study.

Variable	Slope	Mean	StDev	C.V.	Minimum	Maximum	Skewness	Kurtosis
Snow Cover (SC) $(\mathrm{k}{\mathrm{m}}^2$)	Northern Slope	2307	3153	136.65	0.00	14066	1.43	1.40
	Southern Slope	1956	3290	168.23	0.00	13538	1.89	2.68
Minimum Temperature (Tmin) $(°\mathrm{C}$)	Northern Slope	2.19	7.63	347.90	−19.61	14.09	−0.37	−0.87
	Southern Slope	3.32	8.69	261.81	−20.12	16.72	−0.31	−0.98
Maximum Temperature (Tmax) $(°\mathrm{C}$)	Northern Slope	14.52	9.93	68.39	−4.43	30.09	−0.11	−1.39
	Southern Slope	15.87	10.78	67.90	−4.48	32.22	−0.13	−1.43
Global Solar Radiation (GSR) $(\frac{\mathrm{M}\mathrm{J}}{{\mathrm{m}}^2·\mathrm{day}}$)	Northern Slope	19.37	6.93	35.81	8.35	32.11	−0.01	−1.29
	Southern Slope	20.91	7.12	34.06	9.25	33.05	−0.09	−1.34
Relative Humidity (RH) $(\%$)	Northern Slope	0.59	0.16	26.48	0.27	0.92	0.33	−1.05
	Southern Slope	0.48	0.18	38.31	0.22	0.87	0.65	−0.82
Precipitation (P) $(\mathrm{mm}$)	Northern Slope	46.92	32.64	69.55	0.84	160.35	0.74	0.24
	Southern Slope	24.00	19.28	80.36	0.10	81.35	0.97	0.60
Wind Velocity (W) $(\frac{\mathrm{m}}{\mathrm{s}}$)	Northern Slope	2.64	0.38	14.24	1.62	3.66	0.16	0.00
	Southern Slope	2.86	0.51	17.77	1.64	4.11	−0.01	−0.47

refer to the equivalent series.

In this investigation, the drought indices were calculated by precipitation and temperature, recorded during the 1979–2014 period. However, according to the limitation of the snow cover data (because of the nonexistence of snow cover images before 2000), investigating the relationship between drought indices and snow cover was conducted for the period of 2000–2014. Furthermore, the estimation of snow cover was implemented for the 2000 to 2014 period. All of the variables, including the snow cover area and meteorological variables, were provided at a monthly scale for all 12 months of the years under study (Table 1 relates to the monthly scale data).

2.3. Drought Indices

Drought indices are used to monitor drought events. Droughts usually begin by precipitation deficit and affect other water resources, appearing as different types, such as meteorological droughts, agricultural droughts, etc. Different types of droughts can be monitored by special individual indices. For example, the Bhalme and Mooley Drought Index (BMDI) is a meteorological drought index and the Palmer Drought Severity Index (PDSI) is an agricultural drought index. On the other hand, there is the Standardized Precipitation Index (SPI), which refers to different types of drought in its different time windows [11]. There are also multivariate drought indices, such as the Joint Deficit Index (JDI) and Multivariate SPI (MSPI), which can monitor different types of drought simultaneously [12,18].

In the current research, seven different drought indices were employed to investigate the effect of drought on snow cover.

Table 2. Brief information about the drought indices employed in this study.

Name	Developed by	Input Variable(s)	Scale	More Information and Details about the Calculation Steps
Palmer Drought Severity Index (PDSI)	Palmer, 1965 [19]	Precipitation; temperature	Monthly	Palmer, 1965 [19]; Van Der Schrier et al., 2011 [20]
Bahlme and Mooley Drought Index (BMDI)	Bahlme and Mooley, 1980 [21]	Precipitation	Monthly	Bahlme and Mooley, 1980 [21]
Standardized Precipitation Index (SPI)	McKee et al., 1993 [22]	Precipitation	Monthly	McKee et al., 1993 [22]; Svoboda et al., 2012 [11]
Multivariate Standardized Precipitation Index (MSPI)	Bazrafshan et al., 2014 [23]	Precipitation	Monthly	Bazrafshan et al., 2014 [23]; Bazrafshan et al., 2015 [18]; Aghelpour et al., 2020 [12]
Modified Standardized Precipitation Index (SPImod)	Kao and Govindaraju, 2010 [24]	Precipitation	Monthly	Kao and Govindaraju, 2010 [24]
Joint Deficit Index (JDI)	Kao and Govindaraju, 2010 [24]	Precipitation	Monthly	Kao and Govindaraju, 2010 [24]; Bazrafshan et al., 2015 [18]
Standardized Precipitation-Evapotranspiration Index (SPEI)	Vicente-Serrano et al., 2010 [25]	Precipitation; temperature	Monthly	Vicente-Serrano et al., 2010 [25]

shows the drought indices used. Additionally, references for their calculation methods and supplementary information are provided.

2.4. Applied Models

2.4.1. Multiple Linear Regression (MLR)

Multiple linear regression is a linear model implemented between one or more independent variables (input variables) and one dependent variable (target variable). This model is based on the regression coefficients for input variables, which are optimized by the least squares algorithm. In the current research, the inputs were the meteorological variables and the target was the snow cover area. The mathematical form of this model is as follows:

$$y=a+b_1{x}_1+b_2{x}_2+\cdots +b_n{x}_n+e,$$

(2)

where y is the target (dependent) variable, $a$ is the model’s constant coefficient or intercept of the regression line, $b_i$ represents the regression coefficients, $x$ represents the input (independent) variables, and $e$ is the model’s error.

2.4.2. Least Square Support Vector Machine (LSSVM)

Originating from support vector machines (SVMs), LSSVM is a powerful approach for solving nonlinear classification problems, function estimation, and regression. The SVM, developed by Cortes and Vapnik [26], is based on the principle of construction risk minimization (SRM), which reduces the overestimation of the model [27]. The parameters of this model are γ, representing the adjustment constant, and σ2, representing the radial basis function (RBF) kernel width optimized by the least squares algorithm in this study [6].

2.4.3. Group Method of Data Handling (GMDH)

Ivakhnenko presented GMDH in 1970 at the Ukrainian Institute of Technology. This is a nonlinear multivariate statistical analysis method that can be used to determine the complex structure of input–output relationships in complex systems [28]. GMDH is based on the separation of second- or third-order nested polynomials to simulate regression and classification topics [28]. GMDH is a multilayer network in which each layer is composed of multiple nodes. The results from each layer are transferred to the next layer; in fact, the results of one layer are used as inputs for the next layer [29,30].

2.4.4. Multilayer Perceptron (MLP)

The concept of the perceptron was first introduced by McCulloch and Pitts [31] as an artificial neuron. An MLP network provides a nonlinear relationship between the input and output vectors. This is accomplished by connecting neurons from one layer to another one (previous or next layer). The output of each neuron is multiplied by weight coefficients and given as the input for a nonlinear excitation function [32]. In the training phase, the training data are given to the perceptron, and the grid weights are then adjusted to minimize the error between the target and output of the model, or to reach the default number of training times. Then, like all modeling processes, different inputs (not present in the training phase) are used for model validation. The training of neural networks is generally very complex and can be stated to be an optimization problem with a large number of variables [33].

2.4.5. Multilayer Perceptron-Grey Wolf Optimization (MLP-GWO)

This model results from the integration of the MLP model and the GWO algorithm. The GWO algorithm is a population-based meta-heuristic algorithm that is inspired by the behavior of gray wolves during hunting [34]. This algorithm consists of three main steps: (1) Tracking and approaching; (2) pursing and encircling; and (3) attacking. There are four types of wolf in this algorithm. The alpha wolf is the main driver of the algorithm, the beta and delta wolves act as assistants to the alpha wolf, and the omega wolf follows both the alpha and beta wolves. Gray wolves hunt their prey through surrounding them. At the hunting stage, the gray wolf attacks the surrounded prey. The alpha wolf follows the hunting process, and the beta and delta wolves sometimes participate. For mathematically simulating wolf hunting behavior, it is usually assumed that alpha, beta, and delta wolves have better information on the location of prey [35]. Gray wolves attack and hunt when the target stops moving (see [34] for more information on this optimization algorithm). A schematic diagram of the models used to estimate the snow cover area is shown in Figure 2.

2.5. Model Performance Criteria

Several statistical criteria are used in each study to evaluate the performance of the models. These criteria for current research include the Root Mean Square Error (RMSE), the Normalized Root Mean Square Error (NRMSE), Nash–Sutcliff (NS), and the Willmott Index (WI). The mathematical equations of these indices are as follows:

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{i}}-{\mathrm{E}}_{\mathrm{i}}\right)}^2},$$

(3)

$$\mathrm{NRMSE}=\frac{\sqrt{\frac{1}{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{i}}-{\mathrm{E}}_{\mathrm{i}}\right)}^2}}{{\mathrm{O}}_{\max }-{\mathrm{O}}_{\min }},$$

(4)

$$\mathrm{NS}=1-\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{O}}_{\mathrm{i}}-{\mathrm{E}}_{\mathrm{i}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right)}^2},$$

(5)

$$\mathrm{WI}=1-\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{O}}_{\mathrm{i}}-{\mathrm{E}}_{\mathrm{i}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\left(\left|{\mathrm{E}}_{\mathrm{i}}-\overline{\mathrm{O}}\right|+\left|{\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right|\right)}^2}.$$

(6)

In these equations, O_i represents the observational data, E_i represents the estimated data, Ō represents the average observational data, O_max represents the maximum observational data, and O_min represents the minimum observational data. The closer values of WI and NS are to 1, and RMSE and NRMSE are to 0, the better the performance of the models.

3. Results

3.1. Investigating the Relationship between Changes in the Snow Cover Area and Drought Indices

The seven drought indices of PDSI, BMDI, SPI, MSPI, SPImod, JDI, and SPEI were initially calculated. Then, the relationship between these indices and the snow cover area for both the northern and southern slopes of the Alborz Mountains was evaluated by the Spearman correlation test. The results of the Spearman test are presented in

Table 3. Spearmen correlation test results on the relationship between the indexes and snow cover for the northern slope.

Index

SPI1

SPI2

SPI3

SPI4

SPI5

SPI6

SPI7

SPI8

SPI9

SPI10

SPI11

SPI12

Correlation coefficient

0.622 **

0.723 **

0.720 **

0.641 **

0.503 **

0.327 **

0.157

0.001

−0.100

−0.128

−0.111

−0.024

Index

SPI13

SPI14

SPI15

SPI16

SPI17

SPI18

SPI19

SPI20

SPI21

SPI22

SPI23

SPI24

Correlation coefficient

0.112

0.241 **

0.319 **

0.332 **

0.280 **

0.194 *

0.075

−0.045

−0.137

−0.177 *

−0.138

−0.053

Index

SPEI1

SPEI2

SPEI3

SPEI4

SPEI5

SPEI6

SPEI7

SPEI8

SPEI9

SPEI10

SPEI11

SPEI12

Correlation coefficient

0.094

0.097

0.062

0.036

0.010

−0.015

−0.008

−0.020

−0.008

−0.007

−0.006

0.007

Index

SPEI13

SPEI14

SPEI15

SPEI16

SPEI17

SPEI18

SPEI19

SPEI20

SPEI21

SPEI22

SPEI23

SPEI24

Correlation coefficient

0.017

0.011

0.020

0.005

0.005

−0.006

−0.015

−0.014

0.013

0.014

0.038

0.043

Index

SPImod1

SPImod2

SPImod3

SPImod4

SPImod5

SPImod6

SPImod7

SPImod8

SPImod9

SPImod10

SPImod11

SPImod12

Correlation coefficient

0.129

0.131

0.070

0.017

−0.019

−0.028

−0.015

−0.022

−0.010

−0.014

−0.015

−0.011

Index

SPImod13

SPImod14

SPImod15

SPImod16

SPImod17

SPImod18

SPImod19

SPImod20

SPImod21

SPImod22

SPImod23

SPImod24

Correlation coefficient

0.007

0.002

0.008

−0.009

−0.015

−0.026

−0.036

−0.023

−0.007

0.003

0.022

0.031

Index

MSPI1-3

MSPI1-6

MSPI1-9

MSPI1-12

MSPI3-6

MSPI3-12

MSPI6-12

MSPI12-24

MSPI24-48

JDI

BMDI

PDSI

Correlation coefficient

0.085

0.041

0.039

0.018

−0.005

−0.002

−0.021

−0.028

0.037

0.008

0.024

−0.034

* Correlation is significant at the level of 0.05. ** Correlation is significant at the level of 0.01. The bold numbers are the significant correlations.

and

Table 4. Spearmen correlation test results on the relationship between the indexes and snow cover for the southern slope.

Index

SPI1

SPI2

SPI3

SPI4

SPI5

SPI6

SPI7

SPI8

SPI9

SPI10

SPI11

SPI12

Correlation coefficient

0.617 **

0.638 **

0.577 **

0.441 **

0.272 **

0.058

−0.122

−0.267 **

−0.328 **

−0.305 **

−0.198 *

−0.041

Index

SPI13

SPI14

SPI15

SPI16

SPI17

SPI18

SPI19

SPI20

SPI21

SPI22

SPI23

SPI24

Correlation coefficient

0.120

0.222 **

0.252 **

0.209 *

0.115

−0.010

−0.136

−0.246 **

−0.300 **

−0.274 **

−0.188 *

−0.078

Index

SPEI1

SPEI2

SPEI3

SPEI4

SPEI5

SPEI6

SPEI7

SPEI8

SPEI9

SPEI10

SPEI11

SPEI12

Correlation coefficient

0.019

−0.079

−0.186 *

−0.190 *

−0.165 *

−0.124

−0.085

−0.044

−0.033

−0.024

−0.046

−0.032

Index

SPEI13

SPEI14

SPEI15

SPEI16

SPEI17

SPEI18

SPEI19

SPEI20

SPEI21

SPEI22

SPEI23

SPEI24

Correlation coefficient

−0.026

−0.044

−0.051

−0.065

−0.101

−0.100

−0.105

−0.087

−0.033

0.074

0.020

0.017

Index

SPImod1

SPImod2

SPImod3

SPImod4

SPImod5

SPImod6

SPImod7

SPImod8

SPImod9

SPImod10

SPImod11

SPImod12

Correlation coefficient

0.081

0.047

−0.001

−0.065

−0.097

−0.128

−0.100

−0.091

−0.063

−0.039

−0.030

−0.021

Index

SPImod13

SPImod14

SPImod15

SPImod16

SPImod17

SPImod18

SPImod19

SPImod20

SPImod21

SPImod22

SPImod23

SPImod24

Correlation coefficient

−0.020

−0.023

−0.041

−0.070

−0.090

−0.101

−0.096

−0.069

−0.049

−0.035

−0.013

−0.005

Index

MSPI1-3

MSPI1-6

MSPI1-9

MSPI1-12

MSPI3-6

MSPI3-12

MSPI6-12

MSPI12-24

MSPI24-48

JDI

BMDI

PDSI

Correlation coefficient

0.029

–0.027

–0.031

–0.016

–0.061

–0.023

–0.007

–0.017

0.010

–0.040

0.022

–0.064

* Correlation is significant at the level of 0.05. ** Correlation is significant at the level of 0.01. The bold numbers are the significant correlations.

.

Here, 24 time windows of the indices SPI, SPEI, and SPImod were calculated for each slope. The results indicate a significant relationship (p = 0.01) between the snow cover area and SPI for both the northern and southern slopes of Alborz and some SPEI time windows for the southern slope. The highest correlation coefficients were observed in SPI 1-, 2-, 3-, 4-, 5-, and 6-month time windows and accordingly, MSPI was calculated in time windows of 1–3 months (MSPI_1–3), 1–6 months (MSPI_1–6), 1–9 months (MSPI_1–9), and 1–12 months (MSPI_1–12). The five main time windows mentioned in the introduction of MSPI [23] were also calculated and investigated (including MSPI_3–6, MSPI_6–12, MSPI_3–12, MSPI_12–24, and MSPI_24–48). The drought indices, including SPImod (in all time windows), MSPI (in all time windows), JDI, BMDI, and PDSI, were not significantly related to snow cover variations. The highest correlation coefficients for the northern slope were 0.723 and 0.720, which were reported for SPI₂ and SPI₃, respectively. For the southern slope, the highest correlation coefficients were 0.617 and 0.638, which respectively belong to SPI₁ and SPI₂. The most important reasons for this are the higher humidity and cloudiness and consequently lower evaporation rates of the northern slopes, and conversely, the higher received radiation rates for the southern slopes. This causes the snow to melt later on the northern slopes. Therefore, the effects of snow cover on droughts will remain longer and consequently, larger time windows of SPI can be related to the snow cover area for the northern slopes. As a result, SPI₂ can be identified as the best time window of SPI and the drought variable (among these indices) most related to the snow cover variations in this region. In general, it can be reported that the drought status will be most effective during the increasing and decreasing of the snow cover area with a two-month time lag, remaining for up to five (or six) months.

3.2. Results of Snow Cover Estimation

3.2.1. Input Selection

Prior to modeling, the correlation matrix was used at the input selection stage. The correlation results of the variables are shown in Figure 3.

There is a rational and strong negative correlation between the maximum and minimum temperatures, as well as between radiation and snow cover areas, for both slopes. Moreover, in terms of the positive correlations, the relative humidity has a strong correlation with the snow cover and precipitation has a weaker positive correlation coefficient, but both of them are significant at the 0.01 level. The variable wind speed had the weakest correlation for both slopes, being significant at the 0.05 level for the northern slope and insignificant for the southern slope. The opposite direction of the wind speed’s correlation for these two slopes can be related to the influential synoptic systems. The northern slope of Alborz Mountains is affected by the fronts from the Mediterranean Sea, Black Sea, and Caspian Sea. The absolute majority of rainfall due to these systems falls at the northern slope and a much smaller amount falls in the southern zones, which causes a big difference between these two areas’ climates (north of these mountains, there are humid and super-humid climates and, to the south, there is an arid climate). The wind speed close to the earth’s surface, especially at altitude and in mountainous regions, is significantly related to the atmospheric fronts’ speed. Therefore, when the wind speed increases, a front transfers the moisture to altitude and causes snowfall, so the snow cover area will increase. However, the front loses its humidity and after descending the altitudes in the southern part, it usually produces chinook wind. This phenomenon causes snowmelt and decreases the snow cover area, so there will be a negative relation between snow cover and wind speed in this region.

In these datasets, the precipitation refers to rainfall data. Also at high altitudes, the snowfall events are transformed into water equivalents of snow and count as rainfall. It is obvious that snowfall can increase the snow cover and causes a positive correlation with the snow cover area. In contrast, warm precipitation causes a negative correlation. As we can see for the correlations, there are positive correlations (+0.333 for the southern slope and +0.479 for the northern slope) between the snow cover area and precipitation. To extend this issue, it can be stated that, at lower altitudes, warmer seasons have warm precipitation, which melts the snow and causes a negative correlation with the snow cover area, and colder seasons have both cold and warm precipitation. However, almost all of the precipitation at high altitudes occurs in the form of snowfall, which increases the snowy area and has a positive correlation with the snow cover area. Since the snow cover area exists at higher altitudes, cold precipitation and snowfall will be more effective in terms of the snow cover. Therefore, in the general form of this region, precipitation will be positively correlated with the snow cover area.

The correlation of snow cover with meteorological variables is somewhat stronger for the northern slope, because the standard deviation, skewness, and coefficient of variation of this phenomenon (snow cover area) are lower for the northern slope, so it is somewhat closer to exhibiting a normal distribution. This phenomenon can lead to less variability and, consequently, greater consistency with other meteorological variables. To select the inputs, the meteorological variables for each slope were sorted separately by the correlation coefficients, and the input scenarios were selected as given in

Table 5. Input scenarios of the models for each slope.

Zone	Name of the Input Scenario	Input Variables	Scenarios for the Models
			MLR	SVM	GMDH	MLP	MLP-GWO
Northern Slope	Scenario 1	Tmin	MLR1	SVM1	GMDH1	MLP1	MLP-GWO1
	Scenario 2	Tmin; Tmax	MLR2	SVM2	GMDH2	MLP2	MLP-GWO2
	Scenario 3	Tmin; Tmax; RH	MLR3	SVM3	GMDH3	MLP3	MLP-GWO3
	Scenario 4	Tmin; Tmax; RH; GSR	MLR4	SVM4	GMDH4	MLP4	MLP-GWO4
	Scenario 5	Tmin; Tmax; RH; GSR; P	MLR5	SVM5	GMDH5	MLP5	MLP-GWO5
	Scenario 6	Tmin; Tmax; RH; GSR; P; W	MLR6	SVM6	GMDH6	MLP6	MLP-GWO6
Southern Slope	Scenario 1	Tmin	MLR1	SVM1	GMDH1	MLP1	MLP-GWO1
	Scenario 2	Tmin; RH	MLR2	SVM2	GMDH2	MLP2	MLP-GWO2
	Scenario 3	Tmin; RH; Tmax	MLR3	SVM3	GMDH3	MLP3	MLP-GWO3
	Scenario 4	Tmin; RH; Tmax; GSR	MLR4	SVM4	GMDH4	MLP4	MLP-GWO4
	Scenario 5	Tmin; RH; Tmax; GSR; P	MLR5	SVM5	GMDH5	MLP5	MLP-GWO5
	Scenario 6	Tmin; RH; Tmax; GSR; P; W	MLR6	SVM6	GMDH6	MLP6	MLP-GWO6

.

Therefore, input scenarios, inputs (meteorological variables), and the target (snow cover area) were arranged, and in each scenario, 75% of the number of months was selected for training and 25% for the test. Then, five methods, including MLR, LSSVM, GMDH, MLP, and MLP-GWO, were used for estimating the snow cover.

3.2.2. Models’ Performances

This section only covers the test part, and the results are listed in

Table 6. Results of snow cover estimation by the models.

Slope	Input Scenario	Train				Test
		NRMSE	RMSE (km2)	WI	NS	NRMSE	RMSE (km2)	WI	NS
Northern Slope	MLR1	0.116	1631.166	0.914	0.746	0.152	1386.662	0.939	0.762
	LSSVM1	0.097	1366.625	0.948	0.822	0.175	1601.406	0.937	0.683
	GMDH1	0.105	1475.582	0.931	0.792	0.133	1216.993	0.956	0.817
	MLP1	0.100	1402.681	0.948	0.812	0.158	1440.066	0.948	0.744
	MLP-GWO1	0.100	1403.624	0.943	0.812	0.128	1170.852	0.961	0.831
	MLR2	0.115	1621.970	0.916	0.749	0.154	1405.929	0.939	0.756
	LSSVM2	0.083	1169.151	0.964	0.870	0.170	1551.450	0.939	0.702
	GMDH2	0.095	1338.373	0.948	0.829	0.127	1161.180	0.961	0.833
	MLP2	0.088	1239.241	0.958	0.853	0.144	1314.323	0.953	0.786
	MLP-GWO2	0.091	1285.300	0.953	0.842	0.138	1264.108	0.956	0.802
	MLR3	0.110	1544.673	0.927	0.772	0.145	1329.095	0.944	0.782
	LSSVM3	0.081	1144.194	0.965	0.875	0.139	1271.612	0.955	0.800
	GMDH3	0.084	1179.983	0.963	0.867	0.107	980.987	0.971	0.881
	MLP3	0.093	1314.088	0.949	0.835	0.111	1017.925	0.968	0.872
	MLP-GWO3	0.092	1296.527	0.948	0.840	0.094	859.117	0.973	0.909
	MLR4	0.110	1543.266	0.927	0.773	0.145	1323.238	0.945	0.784
	LSSVM4	0.077	1077.685	0.969	0.889	0.126	1151.492	0.959	0.836
	GMDH4	0.080	1120.458	0.967	0.880	0.089	816.011	0.979	0.918
	MLP4	0.075	1051.216	0.971	0.895	0.118	1077.687	0.961	0.856
	MLP-GWO4	0.068	956.988	0.977	0.913	0.092	839.345	0.977	0.913
	MLR5	0.110	1541.353	0.927	0.773	0.147	1338.466	0.944	0.779
	LSSVM5	0.075	1058.106	0.971	0.893	0.133	1214.064	0.953	0.818
	GMDH5	0.077	1087.333	0.969	0.887	0.097	887.575	0.976	0.903
	MLP5	0.078	1096.451	0.966	0.885	0.118	1078.002	0.958	0.856
	MLP-GWO5	0.071	996.632	0.974	0.905	0.116	1059.867	0.963	0.861
	MLR6	0.107	1501.507	0.932	0.785	0.103	937.472	0.971	0.891
	LSSVM6	0.088	1236.120	0.958	0.854	0.106	972.237	0.972	0.883
	GMDH6	0.080	1119.229	0.966	0.880	0.091	829.099	0.978	0.915
	MLP6	0.083	1171.790	0.962	0.869	0.096	881.555	0.975	0.904
	MLP-GWO6	0.086	1204.090	0.959	0.862	0.083	754.711	0.980	0.930
Southern Slope	MLR1	0.148	1999.178	0.876	0.637	0.147	1885.137	0.882	0.646
	LSSVM1	0.124	1678.581	0.921	0.744	0.109	1392.661	0.941	0.807
	GMDH1	0.129	1748.936	0.923	0.722	0.104	1328.284	0.950	0.824
	MLP1	0.128	1726.925	0.923	0.729	0.111	1429.751	0.943	0.797
	MLP-GWO1	0.124	1679.358	0.927	0.744	0.105	1351.662	0.949	0.818
	MLR2	0.146	1975.748	0.882	0.645	0.145	1861.597	0.879	0.655
	LSSVM2	0.125	1688.168	0.919	0.741	0.116	1486.640	0.923	0.780
	GMDH2	0.128	1726.927	0.921	0.729	0.103	1321.075	0.946	0.826
	MLP2	0.126	1709.780	0.923	0.734	0.102	1311.089	0.943	0.829
	MLP-GWO2	0.126	1709.434	0.913	0.734	0.098	1261.877	0.953	0.842
	MLR3	0.146	1974.002	0.890	0.646	0.121	1557.307	0.917	0.759
	LSSVM3	0.110	1484.501	0.941	0.800	0.130	1661.526	0.913	0.725
	GMDH3	0.112	1512.714	0.937	0.792	0.110	1406.948	0.935	0.803
	MLP3	0.116	1574.361	0.933	0.775	0.114	1460.509	0.940	0.788
	MLP-GWO3	0.118	1595.235	0.939	0.769	0.107	1377.436	0.945	0.811
	MLR4	0.146	1982.769	0.889	0.643	0.123	1581.594	0.914	0.751
	LSSVM4	0.102	1379.232	0.950	0.827	0.125	1607.478	0.916	0.743
	GMDH4	0.095	1286.933	0.960	0.849	0.105	1341.722	0.947	0.821
	MLP4	0.103	1393.533	0.952	0.823	0.111	1422.194	0.948	0.799
	MLP-GWO4	0.105	1414.971	0.949	0.818	0.102	1302.612	0.953	0.831
	MLR5	0.147	1987.693	0.888	0.641	0.129	1649.801	0.905	0.729
	LSSVM5	0.100	1358.052	0.952	0.832	0.131	1676.468	0.916	0.720
	GMDH5	0.101	1369.514	0.950	0.830	0.114	1462.780	0.930	0.787
	MLP5	0.117	1587.929	0.923	0.771	0.112	1442.699	0.932	0.793
	MLP-GWO5	0.120	1629.483	0.924	0.759	0.107	1368.283	0.939	0.814
	MLR6	0.146	1978.638	0.889	0.644	0.127	1634.376	0.901	0.734
	LSSVM6	0.108	1468.020	0.941	0.804	0.134	1719.101	0.882	0.706
	GMDH6	0.094	1266.889	0.960	0.854	0.087	1111.549	0.965	0.877
	MLP6	0.113	1535.174	0.940	0.786	0.122	1558.616	0.910	0.758
	MLP-GWO6	0.109	1472.911	0.941	0.803	0.115	1477.708	0.922	0.783

Bold numbers are the best estimations.

. It can be observed that the MLR model is the weakest model in most cases.

This indicates that the nonlinear relationships between the meteorological variables and the snow cover area are stronger than the linear relationships. The best performance in the northern and southern slopes belongs to the MLR6 and MLR3 scenarios, respectively. The equations derived from MLR models are given below:

$$SC=-741.58T_{min}+533.87T_{max}+8990.42RH-152.77GSR+8.09P+1462.76W-10412.99,$$

(7)

$$SC=-841.78T_{min}+538.59T_{max}+7151.76RH-7154.55.$$

(8)

These equations are applicable for these two slopes for estimating the snow cover if no snow cover data are available. Equation (8) is the best MLR model for southern slopes, as it has the least error and can also be used with just three variables, including $T_{min}$, $T_{max}$, and $RH$.

Subsequently, the SVM and MLP artificial intelligence models performed relatively well in terms of the accuracy. Both of these models displayed their best performance with scenario 6 and the model performance can be considered very desirable, with NS > 0.75 and NRMSE < 0.10. However, the NRMSE value for these two models is more than 0.10 for the southern slope.

Among the AI-based models, the GMDH and/or MLP-GWO models produced the best results in all scenarios. GMDH provided its best estimation with the fourth input scenario (GMDH4), exhibiting results of RMSE = 816.011 km², NRMSE = 0.089, NS = 0.918, and WI = 0.979 for the northern slope, and with the sixth input scenario (GMDH6), displaying results of RMSE = 1111.549 km², NRMSE = 0.087, NS = 0.877, and WI = 0.965, for the southern slope. The most accurate performance of MLP-GWO for the northern slope resulted from the MLP-GWO6 input scenario, with RMSE = 754.711 km², NRMSE = 0.083, NS = 0.930, and WI = 0.980. For the southern slope, the MLP-GWO2 input scenario produced the best estimation, with RMSE = 1261.877 km², NRMSE = 0.098, NS = 0.842, and WI = 0.953. The NRMSE < 0.10 and NS > 0.80 for the southern slope and NS > 0.90 for the northern slope show that the performances of GMDH and MLP-GWO can be highly desirable and superior to the other models in the estimation of snow cover.

The GMDH model, with many optimization parameters, such as the number of layers and the number of each layer’s neurons, can provide several make ups to adjust its network structure, whereas the LSSVM model only uses two parameters for optimization (Table 6). The MLP model is similar to the GMDH in terms of having several optimization parameters. However, it uses the GWO algorithm, which can increase the accuracy of MLP (with an average of 10.7% improvement) by selecting the optimal parameters. The optimal parameters of the machine learning approaches are provided in

Table 7. Parameters of the machine learning models.

Zone	Inputs	Models
		LSSVM		GMDH					MLP				MLP-GWO
		γ	σ2	Number of Neurons in Layers					Number of HIDDEN layers	Number of Neurons	Transfer Function (Input to Hidden Layer)	Transfer Function (Hidden Layer to Output)	Quantities and Values
				L1 *	L2	L3	L4	L5
Northern slope	Scenario 1	19.4	7.9	1	-	-	-	-	1	3	Log-sigmoid	Linear	Maximum Number of Iterations = 500 Number of agents = 30 Best search agent = 0.2–1.4
	Scenario 2	109.5	1.3	1	-	-	-	-	1	1	Log-sigmoid	Linear
	Scenario 3	99.2	3.0	3	3	3	3	1	1	3	Log-sigmoid	Linear
	Scenario 4	234.7	9.8	6	15	35	1	-	1	13	Log-sigmoid	Linear
	Scenario 5	905.0	22.1	15	15	15	1	-	1	15	Log-sigmoid	Linear
	Scenario 6	602.8	208.9	15	15	15	1	-	1	23	Log-sigmoid	Linear
Southern slope	Scenario 1	38.6	6.5	1	-	-	-	-	1	8	Log-sigmoid	Linear
	Scenario 2	117.1	35.8	1	-	-	-	-	1	17	Log-sigmoid	Linear
	Scenario 3	96.3	3.5	3	3	3	3	1	1	10	Log-sigmoid	Linear
	Scenario 4	198.7	7.1	6	15	35	35	1	1	12	Log-sigmoid	Linear
	Scenario 5	613.7	15.3	10	35	35	35	1	1	6	Log-sigmoid	Linear
	Scenario 6	110.8	29.9	15	35	35	35	1	1	5	Log-sigmoid	Linear

* L: Layer, i.e., L1 shows the number of neurons in the first layer.

.

Correlation graphs were used to compare the correlation between the snow cover estimation and the actual values (Figure 4 and Figure 5). Here, the weakest coefficient of determination (for both slopes) is also reported in MLR estimation. For the northern slope (Figure 4), the highest R² is observed in the GMDH outputs (R² = 0.92 in GMDH4 and GMDH5), and then in the MLP-GWO model estimation (where R² reaches 0.91 in MLP-GWO3 and MLP-GWO4). The mean coefficients of determination are lower for the southern slope (Figure 5). The weakest coefficient of determination belongs to MLR1 and MLR2, with R² = 0.68. The strongest result of this slope belongs to GMDH6, with R² = 0.88. These values indicate the favorable estimates of the snow cover provided by the GMDH and MLP-GWO models.

Taylor diagrams were plotted for the models’ outputs, as shown in Figure 6 and Figure 7. As there are a large number of models, only two Taylor diagrams were plotted for each slope: The first diagram displays the MLR, SVM, and GMDH outputs, and the second diagram displays the MLP and MLP-GWO outputs. The weakest estimates (points farthest from Obs) are points G and H, which belong to LSSVM1 and LSSVM2, respectively, for the northern slope (Figure 6). The most accurate estimates are the points P, Q, u, and x, which were obtained for GMDH4, GMDH5, MLP-GWO3, and MLP-GWO6, respectively. These points are located closest to the Obs point between the two dashed lines of RMSE = 1000 km² and RMSE = 500 km² and below the R = 0.95 line. For the southern slope (Figure 7), the weakest estimates were obtained for points A and B, which belong to MLR1 and MLR2, and they are close to the R = 0.80 line. The most accurate estimation for this slope is the R point of GMDH6. This point is the closest mode to the R = 0.95 line and RMSE = 1000 km² dashed line. The effect of using the GWO algorithm on MLP optimization can be observed for both slopes, where, in all cases (with similar inputs), the MLP-GWO points are in a better position than the MLP points, indicating the relatively better performance of this model in estimating snow cover areas.

In another graphical form, the frequency distribution of the observed and estimated snow cover data was plotted by a violin plot (Figure 8). Here, the frequency distribution of the most accurate performances of each model is examined in the form of violins [36]. The similarity of the violins to the actual snow cover data violins may indicate that the frequency distribution of the estimated data is close to the actual data. For the northern slope, the GMDH outputs, followed by the MLP-GWO outputs, have the most curvature similar to the observational data violin (particularly GMDH, which performed very desirably in showing data curvature greater than the third quartile). The violin of the southern slope observational data is more complex. This complexity is due to the presence of two months of snow cover at the upper tail end of the distribution. Among the models, the GMDH model was able to estimate one of these two exceptional months, which is beyond the capability of the other four models, particularly MLR and LSSVM. In estimating the violin curvature of the observational data, GMDH generally has the closest possible curvature to the observational data.

4. Discussion

Due to the lack of studies on snow cover area estimation, we can follow studies on the estimation of other snow variables [9], as discussed. Estimates of snow depth in the German Black Forest mountains were obtained using inputs similar to this study (except for the fact that the average temperature was used instead of the maximum and minimum temperatures) [37]. In the study, the MODIS sensor was used to monitor the snow depth and the estimator models were similar to the artificial intelligence-based models. The results of this study with 0.83 < R < 0.93, like the current study, indicated the accuracy of artificial intelligence models for estimating the desired snow variables. The MLP model with precipitation, minimum, and maximum temperature inputs was used to estimate the maximum daily fresh snow accumulation (MDFSA) in South Korea [9]. This model presented a desirable result, with a correlation coefficient of R = 0.90, which is consistent with the results of the present study. Binaghi et al. [38] used the RBF neural network to estimate the thickness (depth) of snow in the central Alpine area in Italy. Meteorological inputs (temperature and precipitation) were used in the study and NRMSE ranged from 0.04 to 0.07, indicating a better accuracy than the current study. This difference may be due to the nature of the variable (snow depth) itself, which can be explained by the fact that its estimation (based on the normalized RMSE) is more accurate than estimating the snow cover area.

5. Conclusions

This study investigated the relationships between several drought indexes, including PDSI, BMDI, SPI, MSPI, SPImod, JDI, and SPEI, and the snow cover area, and compared the accuracies of the artificial intelligence-based models, MLR, LSSVM, GMDH, MLP, MLP-GWO, and MLR in estimating the snow cover area using meteorological inputs, precipitation, the minimum and maximum temperature, the global solar radiation, the relative humidity, and the wind velocity. The evaluation of the correlations between drought indices and the snow cover area showed that the effect of drought on snow cover would remain up to 5 (or 6) months in the region. The most related index is the time window 2 of SPI (SPI₂), which shows that the effectiveness of a drought event will appear after two months based on snow cover variations. The outcomes of this study obtained in the modeling part confirmed that MLP-GWO and GMDH models can be considered acceptable for estimating the snow cover area of a region. The application of these data-driven models is significant in the following cases: (1) Difficulty in field measurements of the snow cover area; (2) a lack of accessibility to satellite images for monitoring the snow cover area; (3) estimation of the snow cover area in missing months; and (4) estimation of the snow cover area in years without satellite images. This study was conducted on a monthly scale, but a limitation for snow studies, at a daily scale, is the separation of snow from cloud. It can cause errors in monitoring the snow cover area and produces missing daily data. For this subject, the 8–day images of MODIS, or images of the other sensors such as Landsat, SPOT, ASTER, etc., can be used. The results of this study in both sections are acceptable and promising, and are well worth researching elsewhere in the world. They suggest using and testing other types of drought indices to further investigate snow cover variations’ effects on droughts. Furthermore, other artificial intelligence-based models, such as ANFIS, RBF, and GRNN, and different optimization algorithms, such as the Genetic Algorithm, Particle Swarm Optimization, Firefly Algorithm, etc., are suggested for snow cover estimation.