2.1 Study area
The Zagros climatic zone covers Iran from the northwest to the southwest (Fig. 1). The region has a temperate and humid climate characterized by seasonal precipitation. According to the climatic classification of (Alijani 1999), the regions consist of three climate classes: sub-humid subtropical, cold, and semi-mountainous areas. The average annual precipitation (P) varies between 250 and 800 mm in different parts of the region. Rainfall historically occurs during winter and spring. Summers are hot and dry (mean temperature: 35 ºC) and winters are cool (mean temperature: 7 ºC) (Fathizadeh et al. 2013). Past studies show that only six synoptic stations in the area have data with adequate quality for this research. Table 1 shows the characteristics of these stations.
Table 1
The synoptic stations used in the study (data from Iranian Meteorological Organization)
Station
|
Longitude
|
Latitude
|
Elevation (m)
|
Sanandaj
|
47.01
|
35.25
|
1373.4
|
Kermanshah
|
47.15
|
34.35
|
1318.5
|
Shahrekord
|
50.84
|
32.29
|
2048.9
|
Shiraz
|
52.63
|
29.56
|
1488
|
Fasa
|
53.72
|
28.90
|
1268
|
Khorramabad
|
48.40
|
33.43
|
1147.8
|
2.2 Dataset
Minimum temperature, maximum temperature, irradiation, and precipitation data were obtained daily from six synoptic stations in the Zagros region. The common statistical basis was considered from 1984 to 2013, and pre-processing of data took place.
2.3 Downscaling statistical data
For this purpose, past and future climates at the six stations were simulated. The basic requirement of the model at the calibration stage is a file that identifies past climate at the stations ( Semenov and Stratonovitch 2010). This file was compiled using daily precipitation data, minimum temperature, maximum temperature, and irradiance at synoptic stations, taking 1984-2013 as the base period. According to the latest IPCC recommendation, the use of multiple models in climatic simulations, instead of running the models individually, can be effective in reducing the uncertainties in the model. To that end, we used four GCMs (BCC-CSM1, CANESM2, HADGEM2-ES, NORESM1-M, Table 2), and three representative concentration pathways (RCP 2.6, 4.5, and 8.5), and calculated the arithmetic mean for each group in R. Finally, precipitation, minimum temperature, maximum temperature, and irradiance were projected for future periods (2010-2039) and (2040-2069) and (2070-2099).
Table 2
List of GCMs in IPCC AR5 (Miao et al, 2014)
GCM
|
Model
|
Resolution
|
Source
|
BCC-CSM1
|
It is a fully coupled global climate-carbon model including interactive vegetation and the global carbon cycle in which the atmospheric component, ocean component, land component, and sea ice component are fully coupled and interact with each other through fluxes of momentum, energy, water, and carbon at their interfaces. The exchange of atmospheric carbon with the land biosphere is calculated at each model time step (20 min).
|
160 × 320
|
Beijing Climate Center, China Meteorological Administration, China
|
CanESM2
|
This is the fourth generation of the atmospheric-ocean public circulation model. There are 40 vertical levels with a distance of 10 meters near the surface (there are 16 levels at a height of 200 meters) to nearly 400 meters deep in the ocean. Spherical horizontal coordinates are used with lattice distance ~ 1.41 ° latitude and 0.94 ° latitude.
|
64 × 128
|
National Center for Atmospheric Research (NCAR), United States
|
HadGEM2-ES
|
The HadGEM2-ES is a coupled ground system model used by the Hadley Met Office to simulate CMIP5. The HadGEM2-ES climatic model includes an atmospheric GCM with a horizontal and vertical resolution of N96 and L38 and an oceanic GCM with a horizontal resolution of 1 degree (rising to 3.1 degrees at the equator) and 40 vertical levels.
|
145 × 192
|
Met Office Hadley Center, Unites Kingdom
|
NorESM
|
The Norwegian Ground System Model (NorESM) is one of 20 climate models producing output for the NorESM1-M CMIP5, with a horizontal resolution of ~ 2 ° for atmospheric and terrestrial components and 1 ° for the ocean and ice components.
|
96 × 144
|
Norwegian Climate Centre, Norway
|
2.4 Standardized Precipitation Index (SPI)
The standardized precipitation index is one of the most important indices in drought monitoring. This index is one of the few time-sensitive drought indices. SPI allows for determining the time scale, depending on which aspects of the impact of droughts are more important (e.g., agricultural resources, hydrological features, etc.). The time scale can range from one month to several years. Also, there are more indicators available for analyzing the severity, duration, and frequency of drought using SPI compared to other indexes (Huang et al. 2016; Lin et al. 2020; Nabaei et al. 2019; Gebrewahid et al. 2017; Teuling et al. 2013; Won et al. 2018). In this study, the time scale was set to one year. Only the precipitation parameter is used in the calculation of SPI. The precipitation of each station is calculated on the timescale. Table 3 shows the classification of SPI and SPEI indices according to (McKee 1995).
Table 3
The Classification of SPI and SPEI indices (McKee, 1995)
Extremely Wet
|
Very wet
|
Moderately Wet
|
Normal Wet
|
Normal Dry
|
Moderately Dry
|
Severely Dry
|
Extremely Dry
|
+2.0
|
1.5-1.99
|
1.0-.149
|
-0.99-0
|
0-0.99
|
-1.0- -1.4
|
-1.5–1.99
|
-2.0 or less
|
2.5 Standardized Precipitation Evapotranspiration Index (SPEI)
SPEI was first introduced by (Serrano et al. 2010). This index is a multivariate index in which precipitation data and evapotranspiration data are combined. This index is calculated similarly to the standardized precipitation index. However, SPI utilizes monthly precipitation, while SPEI utilizes monthly precipitation difference and potential evapotranspiration (PET), representing climate water balance. SPEI has been widely used to study meteorological droughts (Chen et al. 2015; Shao et al. 2018; Paulo et al. 2012; Ren et al. 2012; Wang et al. 2015; Wolf et al. 2011; Yao et al. 2018; Zhao et al. 2015). Table 3 presents how this index is classified. In this study, SPEI was calculated on a one-year time scale in R.
SPEI is derived from equations (1) to (7):
$$m=6.75\times {10}^{-7}{I}^{3}-7.71\times {10}^{-5}{I}^{2}+1.79\times {10}^{-2}I$$
1
$$PET=16K{\left(\frac{10T}{I}\right)}^{m}$$
2
\(i={\left(\frac{T}{5}\right)}^{1.514} (\)3)
$$K=\left(\frac{N}{12}\right)\left(\frac{NDM}{30}\right) \left(4\right)$$
Where T is the average monthly temperature in Celsius, m is the coefficient of dependence on I, I is the sum of 12-month heat indexes, k is the correction factor in terms of month and latitude, NDM is the number of days in a month, and N is the maximum number of hours of sunlight. Thus, with a value for PET, the difference between the precipitation (P) and PET for the month i is calculated:
$${\text{D}}_{\text{i}}={\text{P}}_{\text{i}}-{\text{PET}}_{\text{i}}$$
5
The calculation of SPEI requires a three-parameter distribution. The probability distribution function of series D is based on the following equation:
$$F\left(x\right)=\left[1+{\left(\frac{\alpha }{x-y}\right)}^{\beta }\right]{ }^{-1}$$
6
Where, α, β and 𝑥 are the parameters of scale, shape and position for D values (Singh et al. 1993). The values of F(x) for series D in different time scales are well suited to the observed experimental values, so SPEI can be simply calculated from the standard values of F (x).
\(SPEI=W-\)(C0+C1W+C2W2/1+d1w+d2w2+d3w2) (7)
Where 𝑊 =\(\sqrt{2\text{l}\text{n} \left(\text{P}\right)}\) for 𝑃≤0.5, P being the probability of higher values of D is determined. The values of 𝐶0, 𝐶1 and 𝐶2 and also 𝑑1, 𝑑2 and 𝑑3 are fixed. The mean value of SPEI is zero and the standard deviation is equal to one. SPEI is a standardized variable, so it can be compared with other SPEI values from different places and times. A SPEI value of zero is equivalent to the 50% cumulative probability D (Mavromatis 2007).
2.6 Kernel smoothing method
The kernel smoothing method provides promising results for two simulated and one real-world data set. At the same time, it offers a great alternative for Kriging approaches, in particular in cases where the estimation process has failed. This can be the case when either there are few samples or the observed data are not suitably designed for computer experiments. In this condition, the working unit is a polygon created through Delaunay triangulation. The approach is flexible enough to embed different types of kernel and weight functions. As the method can be easily executed, it can be used as a criterion to validate other interpolation methods (Muhlenst 2009). This method has been used in past studies such as (Hessl et al. 2018; Xu et al. 2020).
Kernel smoothing can be a good option given noisy measurements x(ti) of processes at irregularly spaced times ti. The smoother is a weighted average (Eqs 8 and 9).
$$\text{x}\left(\text{t}\right)=\frac{\sum _{\text{i}=1}^{\text{N}}\text{K}\left(\text{t}-\text{t}\text{i}\right)\text{x}\left(\text{t}\text{i}\right)}{\sum _{\text{i}=1}^{\text{N}}\text{K}\left(\text{t}-\text{t}\text{i}\right)}$$
8
$$K\left(T\right)=\text{exp}\left(-\frac{{T}^{2}}{2{\tau }^{2}}\right)$$
9