Spatiotemporal dynamic of environmental indices of watershed sustainability in connection with land-use change

ABSTRACT

Exploring spatial and temporal land-use changes is valuable for local governments to address issues of sustainability and planning policy where urbanization and industrialization are taking place. Besides anthropogenic effects, natural driving forces like climate change may also affect sustainability. However, such relationships have not been studied minutely. Hence, this study first investigates the land-use changes and their relationship with land surface temperature (LST) for the Shazand Watershed, Iran, in 1986, 1998, 2008, and 2016 coincided with supplementary industrialization stages. Furthermore, the relations among LST and other biophysical parameters, including Normalized Difference Vegetation Index (NDVI), Normalized Difference Buildup Index (NDBI), and Normalized Difference Water Index (NDWI), were analyzed, and corresponding variations were explored. The results indicated that the mean LST of the study watershed has an increasing trend from 1986 to 2008 due to land-use change and drought intensification. Later, LST decreased in 2016. Lower LST was associated with irrigation farming and orchard, and higher LST was related to sparse oak forest areas. There was also a negative correlation between LST and NDVI. As a result, it was inferred that greenery declined LST. Conversely, a positive correlation was found between LST and NDBI resulting from the built-up areas. Since LST could influence biological, physical, chemical processes, it can therefore be supported as an effective index for environmental sustainability assessment.

KEYWORDS:

• Temperature–emissivity separation algorithm
• land surface temperature
• Shazand Watershed
• sustainability index
• sustainability trend
• Watershed monitoring

Introduction

The global ecosystem has been changed by increasing industrialization and urbanization (Peng et al. 2017), leading to land degradation. Facing the effects of land degradation has increased the attention on monitoring, evaluating, and regulating the state and sources of risk and the sustainability degree of ecosystem health (Dernbach and Mintz 2011).

Watershed sustainability refers to a condition in which the use of an ecological system in a manner that satisfies current needs without compromising the needs or options of the future generation (Mirchooli et al. 2021; Teklitz et al. 2021). To achieve the long-term sustainability of human-induced watersheds, interdisciplinary approaches are necessary to assess the long-term social and economic pressure (Fathizad et al. 2017). Different watersheds with various spatial scales must have the same particular framework of sustainability, with their own set of indicators. Different methods for assessing ecosystem sustainability, individual indicator or compound indices, product-related assessment, and integrated assessment tools (Ness et al. 2007). Among these, sustainability assessment using an indicator and index method is a helpful tool to track sustainability trends. Understanding sustainability trends helps prepare short-term projections for the future (Ness et al. 2007). In this regard, the selection of the proper indicator is vital in sustainability assessment. In the same vein, satellite remote sensing has the potential to provide accurate, timely, and lower-cost information about phenomena and natural resources (Zareie et al. 2016). Since land-use maps can explain human interferences such as urbanization and industrial development on the environment, the information on spatiotemporal land-use/cover change is therefore essential, which can help policy/decision-makers apprehend the sustainability-related problems (Ibrahim, Samah, and Fauzi 2012). August and implement revitalization measures especially in degraded areas (Spalevic et al. 2017). Land-use/cover results from climate change and human exploitation (Zhao et al. 2020). They are essential driving forces of socio-economic impact on watershed ecosystems (Asadolahi et al. 2018) that can affect soil erosion (Spalevic et al. 2016), ecosystem services (Peng et al. 2020), water quality (Karcher, VanBriesen, and Nietch 2013), and even soil quality (Hinge, Surampalli, and Goyal 2019). The important biophysical indices such as Normalized Difference Vegetation Index (NDVI), Normalized Difference Built-Up Index (NDBI), and Normalized Difference Water Index (NDWI) are mutated owing to land-use/cover changes (Pal and Ziaul 2017), leading to an alteration in sustainability. It also influences the energy swapping between the land surface and atmosphere (Li et al. 2017; Wang et al. 2016). Accordingly, Land Surface Temperature (LST) is one of the most critical factors excessively controlled land-use/cover change (Wang et al. 2016; Fathizad et al. 2017).

LST is an important factor that controls and affects biological, physical, and chemical processes on the earth’s surface (Shi and Zhang 2017). It affects the biophysical factors such as vegetation, impervious surface, albedo, and, accordingly, the net radiation and water balance (El-Zeiny and Effat 2017). It can also provide important information about the physical characteristic of surface and climate change. Therefore, studying the spatiotemporal dynamic of LST can be of great use as an indicator for sustainability evaluation and land-use planning (Rajendran and Mani 2015).

Remote Sensing (RS) and Geographic Information System (GIS) are also extensively used in LST analysis and land surface (Bokaie et al. 2016). Remote-sensing satellite data, for instance, Landsat (Grover and Singh 2016), ASTER (Wang and Liang 2009), MODIS (Duan et al. 2014), and NOAA AVHRR (Pu et al. 2006), have been used to study thermal environment and to derive biophysical parameters such as vegetation and built up indices. Landsat data have the highest resolution among these images and are easy to obtain (Zhang et al. 2017).

Various facets of LST have been studied to analyze the relationship between land-use/cover and LST by Grover and Singh (2016) in India, Zhang et al. (2017) in China, Mushore et al. (2017) in Zimbabwe, Jafari and Hasheminasab (2017) in Iran and Gohain, Mohammad, and Goswami (2021). The relation of some particular land-use/cover such as industrial land (Rao et al. 2018), urbanized land (Sheng et al. 2017), and type of vegetation cover (Adams and Smith 2014) with LST was also studied. Some researchers (e.g., Son et al. 2012; Sruthi and Aslam 2015) used LST to monitor drought in different areas. The relations among LST and some indices were further investigated by AAli, Marsh, and Smith (2017), Roşca et al. (2017), and Fathizad et al. (2017). LST, Urban Heat Island (UHI), and their relations with land-use/cover changes were analyzed in different urban areas (Uysal and Polat 2015; Li et al. 2016; Zhang et al. 2017; Tarawally et al. 2018). The impacts of landscape composition and pattern on LST and UHI were examined by Zhou, Huang, and Cadenasso (2011), Song et al. (2014), and Estoque, Murayama, and Myint (2017). These studies help better understand the concepts and foundations of LST and its impact across time and space. However, the relations among different biophysical parameters and LST are known, and the spatiotemporal interaction of LST and biophysical parameters (Tayyebi, Shafizadeh-Moghadam, and Tayyebi 2018; Zhang et al. 2017) have rarely been addressed. In contrast, it is vital to manage nexus among different segments of the earth system adaptively.

Therefore, the current study has been formulated to investigate the dynamic of spatiotemporal variation of LST variation in connection with land-use change, using NDVI, NDBI, and NDWI over the last 30 years. The results would be helpful for local managers to quantify human-induced effects on the biophysical condition of the study watershed. Besides that, a general and applicable approach is also developed, which can be applied for drawing adaptable strategies for managing watersheds in developing countries.

Materials and methods

Study area

Geographically, the Shazand Watershed, as one of the sub-watershed of the Namak Lake, is located in Markazi Province, Iran, elongated from 33^◦ 44’ to 34^◦ 12’ N latitude and 49^◦ 4’ to 49^◦ 52’ E longitude with a total area of 1740 km². The climate of the watershed is moderate semi-arid to cold semi-arid according to the Embereger classification, with an average annual rainfall of 420 mm and an average temperature of 12 ^◦C (Davudirad, Sadeghi, and Sadoddin 2016; Sadeghi and Hazbavi 2017; Mirchooli et al. 2021).

The only forest area of the watershed (ca. 0.01 km²) with a very sparse cover of Quercus persica L., a light and symmetry canopy, is located in the southwest. Moreover, some of the primary vegetation covers of the watershed rangeland include Astragalus verus, Stipa barbata, Acanthophyllum microcephalum, and Artemisia sieberi (Talab et al. 2010).

The important industrial areas of Imam Khomeini Refinery, petrochemicals, and thermal power plant and their developmental and supplementary stages were performed in the study watershed over the last 30 years. The study was dissected into four-node years of 1986, 1998, 2008, and 2016 due to data availability and synchronization with industrial development in the area. A general view of the study watershed with supplementary information has been depicted in Figure 1.

Figure 1. General view, elevation distribution (m), and location of the Shazand Watershed in Markazi Province, Iran.

The rapid urban and industrial expansion has been observed in the Shazand Watershed led to a declining land capability (Davudirad et al. 2016; Mirchooli et al. 2021). The exerted problems and pressure by urbanization and industrialization motivated researchers to evaluate the land-use change and other related problems in the watershed (Davudirad, Sadeghi, and Sadoddin 2016; Hazbavi and Sadeghi 2017).

Satellite data processing

In this study, Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI) were collected and used for 1986, 1998, 2008, and 2016 to investigate the spatiotemporal variation of LST and to calculate the remote sensing-based indices as shown in Table. 1. These values have been calculated in ENVI 5.3, ArcGIS 10.3, and Excel 2013.

Table 1. The characteristics of the Landsat image used in the present study

Landsat images	Date	Path/Row	The wavelength of emitted radiance (μm)
Landsat 5 Thematic Mapper (TM)	14. 07. 1986	165/36 165/37	11.45 μm
Landsat 5 Thematic Mapper (TM)	12. 07. 1998
Landsat 7 Enhanced Thematic Mapper Plus (ETM+)	26. 07. 2008
Landsat 8 Operational Land Imager (OLI)	13. 07. 2016		10.895 μm, band 10 12 μm, band 11

Research method

Land-use classification

The land-use map of the Shazand Watershed for 1986, 1998, and 2008 was obtained from existing records (e.g., Davudirad, Sadeghi, and Sadoddin 2016), and the land-use map of 2016 was generated by supervised classification and maximum likelihood algorithm (Davudirad, Sadeghi, and Sadoddin 2016). Some eight land-uses/covers viz. bareland, dry farming, forest, irrigation farming, orchard, rangeland, residential, and outcrops were identified. Nearly 1340 samples were collected during the field survey and divided at 30% and 70% for validation and training. Afterward, necessary geo-referenced, radiometric and atmospheric corrections were performed. The error matrix assessed the overall accuracy of land-use, and the producer’s and user’s accuracy were also calculated.

LST retrieval

The TM, ETM+, and OLI images for July (1986, 1998, 2008, and 2016) were used to generate the LST maps using temperature–emissivity separation algorithm through four following main steps (Sahani, 2021). Considering the Sahzand watershed located in the center of Iran with an arid and semi-arid climate and the date of image in summer (July), it experiences the lowest value of rainfall and wind and consequent ill-effect on image analysis.

Step 1 – Conversion of digital number into spectral radiance

To estimate the LST in Landsat 5 and 7, digital number (DN) was converted to spectral radiance using Equation (1)(1) ${\mathrm{L}}_{\lambda }=\frac{{\mathrm{L}}_{\mathrm{M}\mathrm{a}\mathrm{x},\lambda }-{\mathrm{L}}_{\mathrm{M}\mathrm{i}\mathrm{n},\lambda }}{{\mathrm{Q}}_{\mathrm{C}\mathrm{a}\mathrm{l},\mathrm{M}\mathrm{a}\mathrm{x}}-{\mathrm{Q}}_{\mathrm{C}\mathrm{a}\mathrm{l},\mathrm{M}\mathrm{i}\mathrm{n}}}.({\mathrm{Q}}_{\mathrm{C}\mathrm{a}\mathrm{l}}-{\mathrm{Q}}_{\mathrm{C}\mathrm{a}\mathrm{l},\mathrm{M}\mathrm{i}\mathrm{n})}+{\mathrm{L}}_{\mathrm{M}\mathrm{i}\mathrm{n},\lambda }$(1) :

(1) ${\mathrm{L}}_{\lambda }=\frac{{\mathrm{L}}_{\mathrm{M}\mathrm{a}\mathrm{x},\lambda }-{\mathrm{L}}_{\mathrm{M}\mathrm{i}\mathrm{n},\lambda }}{{\mathrm{Q}}_{\mathrm{C}\mathrm{a}\mathrm{l},\mathrm{M}\mathrm{a}\mathrm{x}}-{\mathrm{Q}}_{\mathrm{C}\mathrm{a}\mathrm{l},\mathrm{M}\mathrm{i}\mathrm{n}}}.({\mathrm{Q}}_{\mathrm{C}\mathrm{a}\mathrm{l}}-{\mathrm{Q}}_{\mathrm{C}\mathrm{a}\mathrm{l},\mathrm{M}\mathrm{i}\mathrm{n})}+{\mathrm{L}}_{\mathrm{M}\mathrm{i}\mathrm{n},\lambda }$(1)

where L_λ is spectral at-sensor radiance (Wm⁻² sr⁻¹ μm⁻¹), L_Min,λ, and L_Max,λ are spectral radiance scaled to Q_Cal.Min and Q_{Cal, Max}, respectively. Q_{Cal, Min} is the minimum quantized calibrated pixel value (corresponding to L_Min), Q_{Cal, Max} is the maximum quantized calibrated pixel value (corresponding to L_Max), and Q_Cal is the quantized calibrated pixel value in digital number (Zareie et al. 2016). The “Spectral Subset” in the “Radiometric Correction” tools was used to select the band. All calculations were performed in the ENVI environment.

Equation (2)(2) ${\mathrm{L}}_{\lambda }={\mathrm{M}}_{\mathrm{L}}\ast {\mathrm{Q}}_{\mathrm{C}\mathrm{a}\mathrm{l}}+{\mathrm{A}}_{\mathrm{L}}-{\mathrm{O}}_{\mathrm{i}}$(2) was also used for calculation of spectral radiance (L_λ):

(2) ${\mathrm{L}}_{\lambda }={\mathrm{M}}_{\mathrm{L}}\ast {\mathrm{Q}}_{\mathrm{C}\mathrm{a}\mathrm{l}}+{\mathrm{A}}_{\mathrm{L}}-{\mathrm{O}}_{\mathrm{i}}$(2)

where M_L represents the band-specific multiplicative rescaling factor, Q_cal is the Band 10 image, A_L is the band-specific additive rescaling factor, and O_i is the correction value for band 10 (El-Zeiny and Effat 2017).

Step 2 – Conversion of spectral radiance to brightness temperature

Spectral radiance was then converted to brightness temperature (T_B, in Kelvin) using a thermal constant in metadata with the help of Equation (3)(3) ${\mathrm{T}}_{\mathrm{B}}=\frac{{\mathrm{K}}_2}{ln\left(\frac{{\mathrm{K}}_1}{{\mathrm{L}}_{\lambda }}\right)+1}$(3) (Artis and Carnahan 1982; Pal and Ziaul 2017):

(3) ${\mathrm{T}}_{\mathrm{B}}=\frac{{\mathrm{K}}_2}{ln\left(\frac{{\mathrm{K}}_1}{{\mathrm{L}}_{\lambda }}\right)+1}$(3)

where K₁ and K₂ are the constant coefficients obtained from metadata (Wm⁻² sr⁻¹ μm⁻¹).

Step 3 – Estimation of land surface temperature (LST)

The brightness temperature was ultimately converted to LST using Equation 4(4) $\mathrm{L}\mathrm{S}\mathrm{T}=\frac{{\mathrm{T}}_{\mathrm{B}}}{1+\left(\lambda \ast \frac{{\mathrm{T}}_{\mathrm{B}}}{\rho }\right)\mathrm{l}\mathrm{n}\epsilon }$(4) (Artis and Carnahan 1982):

(4) $\mathrm{L}\mathrm{S}\mathrm{T}=\frac{{\mathrm{T}}_{\mathrm{B}}}{1+\left(\lambda \ast \frac{{\mathrm{T}}_{\mathrm{B}}}{\rho }\right)\mathrm{l}\mathrm{n}\epsilon }$(4)

where λ is the wavelength of emitted radiance (11.45 μm for Landsat 4, 5 and 7; 10.895 μm for Landsat 8 band 10; 12 μm for Landsat 8 band 11), ɛ is land surface emissivity, ρ = h.c/σ, in which σ is Stefan Boltzmann’s constant (5.67 × 10⁻⁸Wm⁻²K⁻⁴), h is Plank’s constant (6.626 × 10⁻³⁴J S), and c is the light velocity (m s⁻¹) ɛ is land surface emissivity; if NDVI <0.2 then the pixel is considered as bare soil, a mean ɛ is 0.97 if NDVI >0.5; the pixels correspond to dense vegetation areas and ɛ is 0.99 and if 0.2 ≤ NDVI ≤0.5. Each pixel corresponds to the mixing of bare soil and vegetation. ɛ was then computed using Equation (5)(5) $\epsilon ={\epsilon }_{\mathrm{v}}\ast {\mathrm{P}}_{\mathrm{v}}+{\epsilon }_{\mathrm{s}}\ast \left(1-{\mathrm{P}}_{\mathrm{v}}\right)+\mathrm{d}\epsilon$(5) (Sobrino, Raissouni, and Li 2001):

(5) $\epsilon ={\epsilon }_{\mathrm{v}}\ast {\mathrm{P}}_{\mathrm{v}}+{\epsilon }_{\mathrm{s}}\ast \left(1-{\mathrm{P}}_{\mathrm{v}}\right)+\mathrm{d}\epsilon$(5)

which ɛ_v is the emissivity of vegetation and ɛ_s is bare soil emissivity. P_v and dɛ were also calculated using Equations (6)(6) $\mathrm{P}\mathrm{v}=\frac{\left(\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)}{\left(\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)}$(6) and (7(7) $\mathrm{d}\epsilon =\left(1-{\epsilon }_{\mathrm{s}}\right)\ast \left(1-{\mathrm{P}}_{\mathrm{v}}\right)\ast \mathrm{F}\ast {\epsilon }_{\mathrm{v}}$(7) ) (Carlson and Ripley 1997):

(6) $\mathrm{P}\mathrm{v}=\frac{\left(\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)}{\left(\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)}$(6)

(7) $\mathrm{d}\epsilon =\left(1-{\epsilon }_{\mathrm{s}}\right)\ast \left(1-{\mathrm{P}}_{\mathrm{v}}\right)\ast \mathrm{F}\ast {\epsilon }_{\mathrm{v}}$(7)

where F is the form factor averaging 0.55 according to the different geometric distribution.

Equation (5)(5) $\epsilon ={\epsilon }_{\mathrm{v}}\ast {\mathrm{P}}_{\mathrm{v}}+{\epsilon }_{\mathrm{s}}\ast \left(1-{\mathrm{P}}_{\mathrm{v}}\right)+\mathrm{d}\epsilon$(5) was eventually briefed to Equation (8)(8) $\epsilon =0.004\ast {\mathrm{P}}_{\mathrm{v}}+0.986$(8) (Sobrino, Jiménez-Muñoz, and Paolini 2004):

(8) $\epsilon =0.004\ast {\mathrm{P}}_{\mathrm{v}}+0.986$(8)

Step 4 – Unit conversion of LST

The unit of calculated LST was converted from Kelvin (^◦K) to degree Celsius (^◦C) by subtracting a constant value of 273.15.

Validation of LST

To validate estimated LST, the recorded data for two synoptic stations of Arak and Shazand located in and at the vicinity of the study Watershed (Figure 1) were used. Although the Shazand station was established in 2007, the only available and synchronized data of 2016 could be applied for the study. The ground-based data of LST were received from the Markazi Meteorological Bureau.

Development of NDVI, NDBI, and NDWI maps

Landsat 5, 7, and 8 images were employed to calculate NDVI, NDBI, and NDWI. NDVI values vary between −1 and 1, in which negative values were linked to water, low positive values demonstrated bareland and soil, and high positive values pertained to vegetation and greenness (Zhang et al. 2017; Lloyd and Dennison 2018).

The NDBI was sensitive to built-up (residential) areas. As a result, positive values represented built-up area, and negative values exhibited other land-uses. NDBI was expressed using Equation (9)(9) $NDBI=\frac{\left(MIR-NIR\right)}{\left(MIR+NIR\right)}$(9) (Zhang et al. 2017):

(9) $NDBI=\frac{\left(MIR-NIR\right)}{\left(MIR+NIR\right)}$(9)

where MIR is Band 5 (for Landsat TM), and Band 6 (for Landsat 8), and NIR is Band 4 (for Landsat TM), and Band 5 (for Landsat 8).

The NDWI was also used to distinguish the water content of the vegetation canopy, calculated using Equation (10)(10) $NDWI=\frac{\left(NIR-MIR\right)}{\left(NIR+MIR\right)}$(10) (Liu, Yao, and Wang 2016):

(10) $NDWI=\frac{\left(NIR-MIR\right)}{\left(NIR+MIR\right)}$(10)

where NIR is Band 4 (for Landsat TM), and Band 5 (for Landsat 8), and MIR is Band 5 (for Landsat TM), and Band 6 (for Landsat 8).

To examine the relationship between LST and study indices (i.e., NDVI, NDBI, NDWI), some 10,000 random points were designated across the study area, and appropriate modeling was made using a scatter plot. Furthermore, to better understand land-use influences on biophysical parameters and LST, three transect lines were made across the study area to investigate the fluctuations of NDVI, NDBI, NDWI, and LST in different land-uses/covers. These transects were selected to have the most land-use/cover changes along the lines.

Results and discussion

Land-use/cover maps

In the current study, land-uses were monitored from 1986 to 2016. Eight land-use classes (including bareland, dry farming, forest, irrigation farming, orchard, rangeland, residential, and outcrops) were detected in Landsat images using maximum likelihood classification in 1986, 1998, 2008 (Davudirad, Sadeghi, and Sadoddin 2016) and 2016 and mapped in Figure 2 and corresponding data have been summarized in Table 2. In this way, the interpolate accuracy for the land-use map in 2016 was 87.58%. The land-use maps indicated that the largest area of the watershed is covered by rangeland in all study periods.

Table 2. Land-use distribution (km²) from 1986 to 2016 for the Shazand Watershed, Iran

Node year	Rangeland	Dry farming	Irrigation farming	Orchard	Residential	Outcrops	Bareland	Sparse forest
1986	51.29	31.44	5.13	0.11	1.40	8.10	1.89	0.01
1998	40.37	31.74	12.23	2.61	2.85	8.10	2.09	0.01
2008	39.32	27.70	11.81	5.34	4.88	8.10	2.84	0.01
2016	38.31	25.38	11.81	7.79	5.16	8.10	3.44	0.01

Figure 2. Land-use maps in study years in the Shazand Watershed, Iran.

The results (Table 2) showed that settlement and industrial areas increased from 1.40% in 1986 to 5.1% in 2016, verifying human interference in the environment. Regarding the trend of industrial development, there was a conversion and degradation of rangeland in all periods. According to land-use maps, the area of orchards increased simultaneously with industrial development. At the same time, with the increase in the orchard area, bareland was also expanded. Between 1986 and 1998, a remarkable increase was observed in irrigation farming, i.e., some doubles. Furthermore, in 1998–2008, the establishment and related developmental projects of Mohajeran City led to urban extension and encroachment to other land-uses. Overall, the results further showed that the green spaces and orchard areas increased simultaneously as industrial areas were developed.

The spatial pattern of NDVI, NDBI, NDWI, and LST

Given that these indicators are closely related, the spatial value of LST was validated by NDVI and NDBI (Grover and Singh 2016). The NDVI maps of 1986, 1998, 2008, and 2016 are shown in Figure 3. The mean NDVI values for the whole watershed are also illustrated in Table 3. The weighted mean NDVI of the Shazand Watershed was about 0.18. By considering that this index proves the vegetation health, stress, and greenness or biomass (Grover and Singh 2016) and the values more than 0.3 illustrated healthy vegetation (Sahana, Ahmed, and Sajjad 2016), the situation of vegetation cover, especially rangeland (covering 42.32% of the watershed) was not found satisfactory. It indicated a decreasing trend of mean value from 1986 to 2008, corresponding with declining precipitation and drought occurrence. After that, there was a slight increase in 2016, which inferred little improvement in green cover. According to the results, sustainability based on the NDVI indicator had a decreasing trend by 2008. Results also realized that the central and southwest of the watershed had the highest mean value of NDVI in the studied years associated with orchard and irrigation lands. In contrast, outcrops, bareland, and poor rangelands had the lowest mean value of NDVI in different years. The former land-uses with the highest NDVI and the latter with the least NDVI caused the lowest and highest LST, respectively.

Table 3. Mean values of study indices of NI, NDWI, NDBI, and LST in node years in the Shazand Watershed, Iran

Year	NDVI	NDWI	NDBI	LST (^◦C)
1986	0.187	0.005	−0.005	38.82
1998	0.176	0.001	−0.001	40.41
2008	0.170	0.002	−0.002	42.21
2016	0.195	0.030	−0.030	38.41

Figure 3. Spatiotemporal distribution of NDVI in study years in the Shazand Watershed, Iran.

Figure 4 shows that the NDBI had positive value in built-up and bareland areas (no vegetation). The mean values of NDBI were −0.005, −0.001, −0.002, and −0.03 in 1986, 1998, 2008, and 2016, respectively. Changes in NDBI will help decision-makers for management and sustainable usage of the lands, as suggested by Sinha, Verma, and Ayele (2016).

Figure 4. Spatiotemporal distribution of NDBI in study years in the Shazand Watershed, Iran.

The spatial distribution of NDWI was also similar to NDVI in the Shazand Watershed (Figure 5). NDWI and NDBI took similar bands, but the main difference is what pattern they highlighted (Pal and Ziaul 2017), so the distributions of the two indices were similar. The mean values of NDWI are given in Table 3.

Figure 5. Spatiotemporal distribution of NDWI in study years in the Shazand Watershed, Iran.

The LST maps of the study area in 1986, 1998, 2008, and 2016 were ultimately developed, as shown in Figure 6. The mean LST was 38.82, 40.41, 42.21, and 38.41 ^◦C in 1986, 1998, 2008, and 2016. The figures show an increasing trend from 1986 to 2008 due to residential and impervious areas. Furthermore, drought as a natural phenomenon was also a vital driving force that exacerbated the positive trend in LST from 1998 to 2008. However, the effect of other climatic variables such as wind speed on LST has been reported by Tayyebi, Shafizadeh-Moghadam, and Tayyebi (2018) for the metropolitan of Tehran, capital of Iran. The highest value of LST occurred in the western part of the watershed in 1998, and the lowest one occurred in the central part of the watershed in 2008.

Figure 6. Spatiotemporal distribution of LST (^◦C) in study years in the Shazand Watershed, Iran.

To validate the LST retrieval from Landsat images, the ground-based measurements of LST for Arak and Shazand meteorological stations of Markazi Province were collected. The difference between retrieval LST and measured data was 2.81 to 3.84, with the corresponding accuracy beyond 90.98%. The detailed comparisons have been summarized in Table 4. The comparison should not be materialized for 1986, 1998, and 2008 in the Shazand station due to the unavailability of measured data.

Table 4. Estimated LST (°C) and corresponding ground-based data for study stations in the Shazand Watershed, Iran

Station	Variable	1986	1998	2008	2016
Arak	Measured LST Retrieved LST Difference Accuracy (%)	41.10 37.81 3.29 91.99	43.00 40.19 2.81 93.46	42.60 38.76 3.84 90.98	42.50 39.68 2.82 93.36
Shazand	Measured LST Retrieved LST Difference Accuracy (%)	----	----	----	35.20 31.79 3.41 90.31

Relation of NDVI, NDBI, NDWI, and LST

These remote sensing-based indices were helpful to interpret land cover more accurately. The variations of NDVI, NDBI, and NDWI were regressed with LST variations whose relationships have been shown in Table 5. Relations among remote sensing-based indices and mean LSTs were analyzed using linear regression. The results indicated a negative and significant relationship (R2 = 0.25–0.31, P < 0.01) between LST and NDVI.

Table 5. Relationships between LST and NDVI, NDBI, and NDWI in 1986, 1998, 2008, and 2016 based on 10,000 random points

2016	2008	1998	1986
LST = −18.5NDVI+40.9	LST = −19.3NDVI+41.2	LST = −17.5NDVI+40.8	LST = −18.02NDVI+39.4	NDVI
R^2 = 0.85	R^2 = 0.93	R^2 = 0.89	R^2 = 0.91
LST = 25.0NDBI+39.2	LST = 27.0NDBI+42.3	LST = 26.3NDBI+40.4	LST = 24.0NDBI+38.9	NDBI
R^2 = 0.93	R^2 = 0.92	R^2 = 0.82	R^2 = 0.62
LST = −25.0NDWI+39.2	LST = −27.0NDWI+42.3	LST = −26.3NDWI+40.4	LST = −24.0NDWI+38.9	NDWI
R^2 = 0.93	R^2 = 0.92	R^2 = 0. 82	R^2 = 0. 62

NDBI was a sensitive index to the built-up area. It was, hence, positively correlated to LST. Though, Bouhennache et al. (2016) reported that other indices such as urban index (UI) could distinguish the urban and soil features better than NDBI. It is also exhibited a negative correlation between LST and NDWI (R2 = 0.25–0.31, P < 0.01) in different years. It was therefore inferred that the lower LSTs, the more water or moisture. These negative and positive relations among indices might be due to the difference in remote sensing-based indices’ capability in recognizing land-uses. For instance, there was only 5% of watershed covered by urban and industrial areas, even in 2016, with the largest comparative area of urbanization. NDBI was defined even for other land-use classes. Hence, the correlation between LST and NDBI was low. This result is similar to that reported by Ali, Marsh, and Smith (2017). Chen and Zhang (2017) also reported the highest correlation coefficient between LST, NDVI, and NDBI. They further indicated that NDBI and NDVI-MNDWI were good indicators for LST estimation. Guha et al. (2018) found a strong relationship between LST and NDVI in Florence and Naples, Italy. Guha and Govil (2021) indicated a strong relationship between NDVI and LST, insignificant for low and high LST zones.

Overall, the results reflected significant consistency with other studies. Roşca et al. (2017) concluded a strong correlation between LST and NDVI in areas with a lack of vegetation cover. These results also agreed with Pal and Ziaul (2017), Zhang et al. (2017), and Madanian et al. (2018). Moreover, the results were consistent with Ali, Marsh, and Smith (2017) except for NDWI, which showed a positive relation with LST and NDWI.

Land-use/cover and LST relationship

Landsat images are helpful tools for investigating the land-use changes and understanding their effect on LST values in different categories. This study concluded that LST in high vegetation-cover areas was lower than the other land-uses. The mean LST was calculated in various land-use classes in studied years (Table 6). The LST pattern showed the highest value of the mean LST was mainly distributed in the sparse forest and dry farming. Besides, low values of LST were located mainly in irrigation farming and orchard areas. This was due to the influence of vegetation on LST, which was influenced by evapotranspiration and emissivity. These areas had the highest mean value of NDVI. The sparse residential areas with more vegetation and evaporation cooling effect rather than dense built-up areas (Mushore et al. 2017) were found in the study watershed.

Table 6. LSTs (^◦C) in different land-uses and study years in the Shazand Watershed, Iran

Land-use/cover	1986	1998	2008	2016
Rangeland	39.78	41.09	42.68	39.73
Dry farming	39.68	42.61	44.40	40.19
Irrigation farming	30.55	35.55	39.30	33.11
Orchard	28.25	36.26	41.69	37.14
Outcrops	35.18	37.71	37.89	36.62
Bareland	37.85	40.04	40.45	38.24
Residential	37.20	38.82	41.73	36.83
Sparse forest	40.50	42.48	43.02	39.65

Overall, the mean LST increased from 1986 to 2016 in most land-use classes. The minimum value of LST occurred in the orchard in 1986, and the maximum value was located in the sparse forest in 2008. The industrial and urban areas had a lower value of LST. Similar results had been reported by Zareie et al. (2016) for Yazd Region, Iran. Furthermore, the forest areas had the highest value of LST in the study period. A sparse oak forest reserve is located on the southern slope of the watershed. The mean value of NDVI for the oak forest was 0.14, which is almost equal to dry farming and rangeland. The other land-uses with high values of LST were rangeland and dry farming. The mean value of LST in 1998 was more than the values in 1986. All land-use classes had the maximum values of LST in 2008 compared to other study years. Most land-use changes occurred in 1998–2008 (Davudirad, Sadeghi, and Sadoddin 2016). These results agreed with Jafari and Hasheminasab (2017) that approved the impact role of land-use changes on LST. In addition, Tan et al. (2020) showed that the LST values significantly changed across different land-use types. Mukherjee and Singh (2020) indicated an increase in LST of the study areas, as they experienced an increase in built-up areas. Recently, Mustafa et al. (2021) concluded an increasing trend for LST due to land-use change and population density in Beijing, China.

It should be noted that the areas with no changes in land-uses also had different LSTs verifying other affecting factors. In other words, climate change was effective on LST changes in the region located in the arid and semi-arid belt in the northern hemisphere. This could be attributed to drought, which occurred in 2008, so the annual rainfall was the lowest in studied years, as Zaree et al. (2016) mentioned. Besides, the annual rainfall in 2008 was lower than in other years in most meteorological stations located in the study watershed, as shown in Figure 7. The same land-use classes had the lowest and the highest values of LST due to the reasons mentioned above in 2016.

Figure 7. The average rainfall of the stations located in the watershed for the period 1986–2016.

As depicted in Figure 8, three transects were arbitrarily considered across the watershed to study the effects of land-uses on biophysical parameters and LST. To this aim, the profile graph was created using an interpolated line in the ArcGIS environment. In general, it is inferred from these transects that the values of LST have increased as a result of land-use changes and indices in different parts of the profile. The southern profile line (Figure 6, Bottom-right, CD) crossed rangeland in the center and east, and then dry farming and forest areas in the west. The LST in irrigation farming in the west of the watershed was the lowest value. Transect AB corresponded to the highest NDVI and NDWI. The northern profile line (AB) passed through dry, industrial, irrigation, and dry farming from east to west. The vertical profile line (EF) crossed through rangeland, industrial, dry farming, and rangeland areas from north to south. The detailed results of LST values along the profile lines in 2016 have been shown in Figures 8–10.

Figure 8. Variations of LST, NDVI, NDBI, and NDWI along with the transect AB (Figure 6, bottom-right) for 2016 in the Shazand Watershed, Iran.

Figure 9. Variations of LST, NDVI, NDBI, and NDWI along the transect CD (Figure 6, bottom-right) for 2016 in the Shazand Watershed, Iran.

Figure 10. Variations of LST, NDVI, NDBI, and NDWI along the transect EF (Figure 6, bottom-right) for 2016 in the Shazand Watershed, Iran.

Conclusion

The dynamics of LST was examined in the current study in connection with land-use/cover and some biophysical parameters using multi-temporal Landsat images. The results showed a negative correlation between NDVI and LST and a positive correlation between NDBI, NDWI, and LST. Besides, most changes were related to the conversion of rangeland to irrigation farming due to urbanization and industrialization in the study period of some 30 years. The land-use conversions resulted in a progressing mean LST in the studied node years from 1986 to 2016. The spatiotemporal dynamic of LST showed that hotspot areas with the highest LST were associated with sparse oak forest and rangeland, resulting from the effects of land-use change and climate change. It can be inferred from the results of correlations among LST and biophysical parameters that the protection and maintenance of greenery spaces as an essential variable in moderate climate are necessary. So politicians and decision-makers should adopt greenery strategies and protect the oak forest against degradation as a practical approach for reducing hotspot effects. The present study results on land-use changes will provide information about the urbanization development process and help policymakers and regional managers adoptively manage the watershed ecosystem. Moreover, identifying the areas with the highest LST helps decision-makers for proper environmental monitoring and sustainability studies and could lead to management decisions to protect natural resources. Of course, more insight studies with a more extended period and high-resolution information and data are needed to draw an integrated conclusion.

Data availability

Some or all data, models, or code generated or used during the study are available in a repository or online in accordance with funder data retention policies. (https://earthexplorer.usgs.gov/)

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was partially supported by the Tarbiat Modares University Agrohydrology Research Group [Grant No. IG-39713], Iran, concerning the corresponding author.