2.1 The study area
The seven sub-watersheds are located between the Iraq – Iran and have a total area of about 3770 km2 (Fig. 1). The main watershed occupies vast area of the western Iran and northeastern Iraq. The study area is characterized by several, relatively large, peripheral valleys along with numerous smaller ones. Most of these valleys originate within western Iran and flow toward eastern Iraq. The Iraqi part of the catchment is mainly flat land. Its elevation ranges from 400 m above mean sea level (a.m.s.l.) at the international border line to about 25 m (a.m.s.l.) at the southeast end of the region. The Iranian part consists of high lands with an elevation of more than 2500 m (a.m.s.l.). It includes a series of hills and mountain ranges with a basin divide runs along the peaks within Iran (Rahi et al. 2019). The dominant climate in the Iraqi part of the watershed is a subtropical desert climate (BWh type of Köppen-Geiger climate classification system), and of type BSh (hot semi-arid) in the Iranian part of the sub-watersheds. Ali Al-Gharbi has a yearly temperature of 30.79°C, which is 4.02% higher than Iraq's average. The average annual precipitation is 16.85 mm and there are 32.41 rainy days. In Dehloran in Iran, the northern part of the sub-watershed, temperatures hover around 19°C during the day and 11°C at night. The average annual rainfall is 20.48 mm, and the humidity is close to 41%. According to the digital soil map of the world (Sanchez et al. 2009), the study area has three main soil textures: clay loam, loam, and sand, (Fig. 2). The three textures occupy 1257 km2 (33%), 2090 km2 (55%), and 423 km2 (11%), for clay loam, loam, and sand, respectively. The distribution of loam and sand textures which belong to the B and A hydrological groups, indicates that the soil is highly permeable, resulting in high infiltration rates and low runoff. Loam and sandy soils cover the northern and eastern portions of the study area, while clayey loam covers the middle and large portion of the southern part, indicating that these parts are likely to generate runoff and floods. The land use/land cover (LULC) map for the study area is shown in Fig. 3. This map represents the LULC for the year 2020 and is gotten from the ESRI web site (Sentinel-2 10m Land Use/Land Cover Timeseries Downloader (arcgis.com). The original global map, developed from Europe Space Agency (ESA) Sentinel-2 imagery at 10 m by Impact Observatory's deep learning AI land classification model, is based on a massive database of National Geographic Society-labeled image pixels. Six LULC classes, namely, water, trees, crops, built area, shrub, and bare area are found in the study area. The most dominant class is bare area with 2792 km2 (about 74% of the total area), followed by crops with 486 km2 (13%), and shrub class with 467 km2 (12%). The other classes only cover a small part (1%) of the study area.
From the geological point of view, the Mesopotamian Foredeep basin and a small part of Zagros Fold belt occupy the Iraqi part of the study area, whereas the Zagros Fold Belt encompasses the Iranian Territory. The Mesopotamian Foredeep basin is an elongated epicontinental basin that formed over an earlier continental or migration basin (Fouad 2010). It contains several buried structures, including folds, faults, and diapiric structures. On the other hand, the Zagros fold is a zone of deformed crustal rocks, arising from the collision of the Arabian the Arabian margin and the Eurasian plate following the closure of the Neo-Tethys ocean during the Tertiary (Falcon 1974; Le Garzic et al. 2019).
2.2 Work framework
To assess the watershed health in the considered 7 sub-watersheds, five steps were adopted (Fig. 4): (i) the monthly average of grided rainfall data of type PERSIANN (Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks) system developed by the Center for Hydrometeorology and Remote Sensing (CHRS) at the University of California, Irvine (UCI) was gathered firstly for the period 1980–2020 and for 9 locations inside and around the studied sub-watershed (Fig. 1). The PERSIANN system uses neural network function classification/ approximation procedures to compute an estimate of rainfall rate at each 0.25° x 0.25° pixel of the infrared brightness temperature image provided by geostationary satellites (Sorooshian et al. 2000; Nguyen et al. 2019). (ii) The non-parametric Sen’s slope and Mann-Kendall tests were used to assess the rate of change in rainfall data and the significance of rainfall time series trends, respectively. (iii) The SPI drought index was estimated using the rainfall data on monthly time scale using DrinC 1.7 software (Tigkas et al. 2015). (iv) the RRV indicators were calculated based on the estimated SPI values. (v) the watershed health index (WHI) was calculated based on the RRV indices and the spatial and temporal variation of WHI were presented as maps and discussed to show the main trend of it.
2.3 Rainfall trend analysis
In statistics, trend analysis determines whether or not a time series of a random variable (in this case, monthly rainfall) is increasing or decreasing over time or whether its probability distribution has shifted over time (Yenigün et al. 2008). Parametric, non-parametric, and mixed methods are often used to study the trend of a time series of variable. In hydrology, non-parametric analysis is preferred over parametric analysis due to the robustness of non-parametric techniques in dealing with missing or tied data, seasonality, non-normality, non-linearity, and serial correlation (McLeod et al. 1991). For the purpose of this study, the non-parametric Mann-Kendall and Sen’s slope estimator were used in trend detection of rainfall dataset. More Mathematical detail about these tests can be found in Gilbert (1987).
2.4 SPI calculations
The Standardized Precipitation Index (SPI) is the most commonly used indicator worldwide for detecting and characterizing meteorological droughts. The SPI index, which was developed by McKee et al. (1993), and described in detail by Edwards (1997), measures precipitation anomalies at a given location, based on a comparison of observed total precipitation amounts for an accumulation period of interest (e.g. 1, 3, 12, 48 months), with the long-term historic rainfall record for that period. The historic precipitation data is fitted to a gamma probability distribution, which is then transformed into a normal distribution using an equal probability transformation. (Guttman 1999). Positive SPI numbers suggest wet circumstances with more precipitation than the median, while negative SPI values indicate dry conditions with less precipitation than the median (Belayneh and Adamowski 2012).
In most cases, the probability distribution that fit climatological precipitation time series is gamma distribution (Thom 1966). The gamma distribution is defined by its frequency or probability density function as:
\(g\left(x\right)=\frac{1}{{\beta }^{\alpha }{\Gamma }\left(\alpha \right)}{x}^{\alpha -1}{e}^{-x/\beta }\) , for x > 0, (1)
where \(\alpha >0, \beta >0, and x>0\) are the shape parameter, the scale parameter, and the amount of precipitation, respectively. \({\Gamma }\left(\alpha \right)\) is the Gamma function and defined as (Cacciamani et al. 2007):
$${\Gamma }\left(\alpha \right)={\int }_{0}^{\infty }{y}^{\propto -1}{e}^{-y}dy$$
2
To model the observed precipitation data with a gamma probability function, it is necessary to estimate \(\alpha\) and \(\beta\) parameters. The maximum likelihood solutions are frequency used to optimally evaluate these parameters according to the following formulas (Thom 1966):
$$\widehat{\propto }=\frac{1}{4A}\left(1+\sqrt{1+\frac{4A}{3}}\right)$$
$$\widehat{\beta }=\frac{\stackrel{-}{x}}{\widehat{\propto }}$$
3
$$A=ln\left(\stackrel{-}{x}\right)-\frac{\sum ln\left(x\right)}{n}$$
4
where n is the number of precipitation observations
The estimated \(\alpha\) and \(\beta\) parameters are utilized to find the cumulative probability of an observed precipitation event for the given month and time scale by integration the \(g\left(x\right)\) with respect to x (Cacciamani et al. 2007):
$$G\left(x\right)={\int }_{0}^{\infty }g\left(x\right)dx=\frac{1}{\widehat{\beta }{\Gamma }\left(\widehat{\propto }\right)}={\int }_{0}^{x}{x}^{\widehat{\propto }-1}{e}^{-x/\beta }dx$$
5
Equation 5 becomes the incomplete gamma function by letting\(t=x/\widehat{\beta }\)
$$G\left(x\right)=\frac{1}{{\Gamma }\left(\widehat{\propto }\right)}{\int }_{0}^{x}{t}^{\widehat{\propto }-1}{e}^{-t}dt$$
6
Because the gamma function is undefined at x = 0 and a precipitation distribution can have zeros, the cumulative probability becomes: (Edward and McKee 1997)
$$H\left(x\right)=q+\left(1-q\right)G\left(x\right)$$
7
where q is the probability of a zero.
The \(H\left(x\right)\) is transformed to the standard normal distribution with mean zero and variance of one (Z score) to get the value of SPI.
Due to the difficulty involved with the approach described above and the need to produce figures for all stations at all time scales and for each year (Edward and McKee 1997), the Z or SPI value is more easily obtained computationally using an approximation provided by Abramowitz et al. (1988) that convert cumulative probability to the standard normal random variable Z:
$$Z=SPI=-\left(t-\frac{{c}_{0}+{c}_{1}t+{c}_{2}{t}^{2}}{1+{d}_{1}t+{d}_{2}{t}^{2}+{d}_{3}{t}^{3}}\right)$$
for \(0<H\left(x\right)\le 0.5\) (8)
$$Z=SPI=+\left(t-\frac{{c}_{0}+{c}_{1}t+{c}_{2}{t}^{2}}{1+{d}_{1}t+{d}_{2}{t}^{2}+{d}_{3}{t}^{3}}\right)$$
for \(0.5<H\left(x\right)\le 1.0\) (8)
where
$$t=\sqrt{ln\left(\frac{1}{{\left(H\left(x\right)\right)}^{2}}\right)}$$
for \(0<H\left(x\right)\le 0.5\) (9)
$$t=\sqrt{ln\left(\frac{1}{{\left(1-H\left(x\right)\right)}^{2}}\right)}$$
for \(0.5<H\left(x\right)\le 1.0\) (8)
and \({c}_{0}\)=2.515517, \({c}_{1}\)= 0.802853, \({c}_{2}\)=0.010328, \({d}_{1}\)=1.432788, \({d}_{2}\)=0.189268, \({d}_{3}\)=0.001308.
In this study, the SPI values were calculated using the DrinC software package (Tigkas et al. 2015) using a monthly time scale.
2.5 The SPI critical value in the RRV calculation
In this study, SPI thresholds for RRV indicators are mostly determined by hazard intensity contained within the SPI index. The SPI hazard classes become persistently negative when the SPI, whichever timescale is investigated, reaches a value of -1 (McKee et al. 1993). Despite this, there is not a standard, as some researchers chose a threshold that is less than zero, but not quite − 1, while others classified drought as less than − 1 (Sadeghi and Hazbavi 2017). For the present research, the value of 0.1 is adopted from the previous studies ( Hazbavi et al. 2018; Hazbavi and Sadeghi 2017; Sadeghi and Hazbavi 2017), and for the areas that are closest to the study area and have similar climatic conditions (arid to semi-arid).
2.6 RRV indicators calculation
In the context of meteorological drought, the concept of RRV is applied through analyzing the temporal variation of the SPI. As previously stated, the system is in a satisfactory condition when the SPI is above the threshold (0.1 here); and when the SPI is below the threshold, the system is in unsatisfactory condition (failure event). Reliability (Rel) is defined by the probability that a system is in a satisfactory state (Hashimoto et al. 1982). Mathematically, it can be defined as:
$${R}_{el}=1-\frac{\sum _{j=1}^{M}d\left(j\right)}{T}$$
(9)
The failure event duration is given by \(d\left(j\right)\), M is the number of failure events, and T is the total number of time intervals (12 here).
The second indicator, Res is a measure of how quickly a watershed can return to a satisfactory level after falling below the threshold. It can calculate using the following equation (Al-Mohammdawi et al. 2022):
$${R}_{es}={\left\{\frac{1}{M}\sum _{j=1}^{M}d\left(j\right)\right\}}^{-1}$$
(10)
Finally, the Vul refers to the likely magnitude of a failure, if one occurs (Hashimoto et al. 1982). It can be computed as:
$${V}_{ul}=\frac{1}{M}\sum _{j=1}^{T}\left\{\left[\frac{{L}_{obs}\left(i\right)-{L}_{std}\left(i\right)}{{L}_{std}\left(i\right)}\right]\right\}H\left[{L}_{obs}\left(i\right)-{L}_{std}\left(i\right)\right]$$
(11)
The \({L}_{obs}\left(i\right)\) is the calculated SPI at the ith time step, \({L}_{std}\left(i\right)\) is the critical threshold of SPI (0.1 in our case), and H( ) is the Heaviside function which ensures that only failure events are involved in the Vul calculation (Hazbavi and Sadeghi 2017).
Both reliability and resiliency indicators have values between 0 and 1 (probability values), whereas the vulnerability could take any value. The larger the reliability and resiliency values, the better the watershed health is. On the contrary, the larger the vulnerability, the worse the health of the watershed. However, before we used these indicators to assess the watershed health, these indicators should be standardized to reflect the better or worst situations. In this study, the three indicators were standardized to refer to the better situation. The following equations were utilized to standardize the three indicators: (Al-Abadi et al. 2017)
For large the better (the cost type): for reliability and resiliency indicators
$${y}_{i}=\frac{{x}_{i}-{x}_{i\left(min\right)}}{{x}_{i\left(max\right)}-{x}_{i\left(min\right)}}$$
(12)
For the smaller the better (the efficiency type): for vulnerability indicator
$${y}_{i}=\frac{{x}_{i\left(max\right)}-{x}_{i}}{{x}_{i\left(max\right)}-{x}_{i\left(min\right)}}$$
(13)
where \({y}_{i}\) is the standardized indicator value, \({x}_{i}\) is the original indicator value, \({x}_{i\left(min\right)}\) and \({x}_{i\left(max\right)}\) are the non-standardized minimum and maximum values of the indicator, respectively.
2.7 Watershed health index
The watershed health index (WHI) for each station was calculated based on the RRV for each sub-watershed using the geometric mean according the following equation: (Sadeghi and Hazbavi 2017)
$$WHI={\left({R}_{el}\times {R}_{es}\times {V}_{ul}\right)}^{\frac{1}{3}}$$
(14)
The WHI is classified into five equal classes to reveal the healthy status of the meteorological point (Sadeghi and Hazbavi 2017): 0-0.2 (very unhealthy), 0.21–0.4 (unhealthy), 0.41–0.60 (moderate healthy), 0.61–0.80 (healthy), and 0.80-1.00 (very healthy).
To reveal the spatial distribution of WHI across the whole study area and for the studied period, the calculated WHI was interpolated using well known deterministic inverse distance weighted (IDW) interpolation technique in the ESRI ArcGIS map 10.8.1. IDW is one of the simplest and most popular interpolation techniques (Babak and Deutsch 2009). IDW predicts values for any unmeasured location based on measured values surrounding that location. Each measured point is assumed to have a local influence that diminishes with distance. The IDW model gives greater weight to points closest to the prediction location and diminishes as the distance increases (Shepard 1968; Borges et al. 2016).