Impacts of Vegetation and Topography on Land Surface Temperature Variability over the Semi-Arid Mountain Cities of Saudi Arabia

Land surface temperature (LST) can fully reflect the water–heat exchange cycle of the earth surface that is important for the study of environmental change. There is little research on LST in the semi-arid region of Abha-Khamis-Mushyet, which has a complex topography. The study used LST data, retrieved from ASTER data in semi-arid mountain areas and discussed its relationship with land use/land cover (LULC), topography and the normalized difference vegetation index (NDVI). The results showed that the LST was significantly influenced by altitude and corresponding LULC type. In the study area, during the summer season, extreme high-temperature zones were observed, possibly due to dense concrete surfaces. LST among different types of land use differed significantly, being the highest in exposed rocky areas and built-up land, and the lowest in dense vegetation. NDVI and LST spatial distributions showed opposite trends. The LST–NDVI feature space showed a unique ABC obtuse-angled triangle shape and showed an overall negative linear correlation. In brief, the LST could be retrieved well by the emissivity derived NDVI TES method, which relied on upwelling, downwelling, and transmittance. In addition, the LST of the semi-arid mountain areas was influenced by elevation, slope zenith angle, aspect and LULC, among which vegetation and elevation played a key role in the overall LST. This research provides a roadmap for land-use planning and environmental conservation in mountainous urban areas.


Introduction
Mountains are well-known drivers of climate change and growth in global mass and energy balance. Extreme topographical variability introduced in mountain areas creates complex spatial-temporal patterns of meteorological and hydrological process, and ecological parameters with high gradients of attitude at short distances [1][2][3][4]. Mountain ecosystems are also very vulnerable to natural and anthropogenic changes, compared to lowland ecosystems. Mountain ecosystem restoration from disruptions is usually slow or failing. Analysis of atmospheric surface processes in mountain regions heterogeneous landscape in terms of the terrain complexity. It was situated in an essential zone of Afromontane, where cold and semi-arid climates characterize the city's environment [61]. The mean annual rainfall for the last 55 years (1965-2019) was 245 mm, among which most of the precipitation occurred during the months of February to June, while the mean minimum and maximum air temperature were 9.4 °C and 30.8 °C, respectively. The watershed is vulnerable to high-intensity rainfall, and several rural areas undergo flash flooding in or during the winter season [62].

Methodology
In this study, ASTER was used and radiometric calibration and geometric co-registration was carried out and used for the preparation of land use/land cover, vegetation density map (NDVI), and estimating the LST. Topographic elevation map was developed for the study area using the Digital Elevation Model (DEM) obtained from ALOS PALSAR RTC data [63]. The raw DEM data was processed by filling the sinks to cover local depressions. Slope, aspect, and drainage network maps were derived directly from the elevation data (i.e., DEM). A field survey was conducted from 24-27 June 2019, comparing the estimated satellite surface temperature values with ground-based measurements (in-situ). The ArcGIS and image processing software, i.e., ArcGIS 10.3, ENVI 4.8, ERDAS imagine 9.2, FLAASH 4.8, and Microsoft Office 2010 software and statistics computation, SPSS 24 software were used to carry out this study. For field observations of surface temperature, infrared thermometer and Garmin-GPS were used. Figure 2 demonstrates the flow diagram for the adopted methodology.

Methodology
In this study, ASTER was used and radiometric calibration and geometric co-registration was carried out and used for the preparation of land use/land cover, vegetation density map (NDVI), and estimating the LST. Topographic elevation map was developed for the study area using the Digital Elevation Model (DEM) obtained from ALOS PALSAR RTC data [63]. The raw DEM data was processed by filling the sinks to cover local depressions. Slope, aspect, and drainage network maps were derived directly from the elevation data (i.e., DEM). A field survey was conducted from 24-27 June 2019, comparing the estimated satellite surface temperature values with ground-based measurements (in-situ). The ArcGIS and image processing software, i.e., ArcGIS 10.3, ENVI 4.8, ERDAS imagine 9.2, FLAASH 4.8, and Microsoft Office 2010 software and statistics computation, SPSS 24 software were used to carry out this study. For field observations of surface temperature, infrared thermometer and Garmin-GPS were used. Figure 2 demonstrates the flow diagram for the adopted methodology.

ASTER Sensor
ASTER is a multi-spectral high-spatial-resolution image that was launched on NASA's Terra spacecraft in 1999, which transmits data from a visible to the thermal infrared region, covering an area of 3600 km 2 in 14 spectral bands. ASTER comprises three different spectral band categories; in the first category there are four bands with a spatial resolution of 15 m in the region of visible and near-infrared (VNIR); in the second category there are 6 bands with a spatial resolution of 30 m in the region of shortwave infrared (SWIR); and in the third category, there are five bands with a spatial resolution of 90 m in the region of thermal infrared (TIR) [64].

Radiometric Calibration
Radiometric calibration is crucial to the development of high-quality scientific data. The ability to detect and quantify changes in the Earth's environment depends on the sensors, which can quantify the surface of the earth over time. A comprehensive description of scientific data from a global range of remote sensor data requires the identification of product artefacts and the study of earth process changes [65]. Determining the spectral emission of the sensor is the crucial step for the transfer of image data from the number of sensors to a standard radiometric scale. "The radiometric calibration of the ASTER sensors involve the rescaling of the raw digital numbers (Q) transmitted from the satellite to the calibrated digital numbers (Qcal), using the same radiometric scale for all scenes processed on the ground for a specific period".

ASTER Sensor
ASTER is a multi-spectral high-spatial-resolution image that was launched on NASA's Terra spacecraft in 1999, which transmits data from a visible to the thermal infrared region, covering an area of 3600 km 2 in 14 spectral bands. ASTER comprises three different spectral band categories; in the first category there are four bands with a spatial resolution of 15 m in the region of visible and near-infrared (VNIR); in the second category there are 6 bands with a spatial resolution of 30 m in the region of shortwave infrared (SWIR); and in the third category, there are five bands with a spatial resolution of 90 m in the region of thermal infrared (TIR) [64].

Radiometric Calibration
Radiometric calibration is crucial to the development of high-quality scientific data. The ability to detect and quantify changes in the Earth's environment depends on the sensors, which can quantify the surface of the earth over time. A comprehensive description of scientific data from a global range of remote sensor data requires the identification of product artefacts and the study of earth process changes [65]. Determining the spectral emission of the sensor is the crucial step for the transfer of image data from the number of sensors to a standard radiometric scale. "The radiometric calibration of the ASTER sensors involve the rescaling of the raw digital numbers (Q) transmitted from the satellite Atmosphere 2020, 11, 762 6 of 28 to the calibrated digital numbers (Qcal), using the same radiometric scale for all scenes processed on the ground for a specific period".
Digital number (DN) values conversion to radiance (at-sensor) using Equation (1) is expressed in "where, UCC is unit conversion coefficients from the HDF file" (Source: ASTER User Handbook [42]).

Atmospheric Correction for the ASTER Channel
The electromagnetic energy that goes through the atmosphere attenuates the radiation before reaching the sensor. Therefore the elimination of the atmospheric influences is important. ASTER's atmospheric thermal data correction uses a clear-sky approach that cannot be applied near the clouds. This is based on the equation of radiative transfer and uses the "MODerate resolution atmospheric TRANsmission (MODTRAN)", Version 3.5 radiative transfer code. The system includes all the required standard atmospheric and geometric parameters, including temperature, channel reflection, water vapor, aerosol, CO 2 mixing ratio (ppm), ozone, number of disort stream, radiance image scale factor, etc. This atmospheric module used an algorithm that was based on observations by Kaufman et al. [6], and the parameters were associated with each image's geometry as input for atmospheric correction. For more information, please see Palluconi et al. [66].

Retrieval of Emissivity and LST
In terms of emissivity variation, natural surfaces at a pixel scale are not homogeneous. The surface emissivity value of Earth shows a strong seasonal presence and dependence on the landscape. "Emissivity depends in particular on the type of land cover, soil moisture, organic soil composition, vegetation density and structure i.e., broadband emissivity for densely vegetated areas, which are usually around 0.96-0.98, but for bare soils, it can be lower than 0.90 [67]" because of the variations in TES's scaling of high-low MMD emissivity spectra. It is appropriate, irrespective of the TIR data, to use ASTER's VNIR data to estimate vegetation proportion cover and surface emissivity. Various methods for determining emissivity from NDVI were formulated [58,68,69].
ASTER satellite data were used to derive the land surface emissivity (LSE) by taking the vegetation cover ratio per pixel in combination with NDVI. The "proportion of vegetation cover" was calculated as given in Equation (2) [70]: "where, i g = NDVI value of pure soil; i v = NDVI value of pure vegetation pixel; i is the NDVI value of mixed pixels".
"where ρNIR v − ρRED v are reflectance in the NIR (band5) and red (band4) for pure pixels of vegetation; ρNIR g − ρRED g are reflectance in the NIR (band5) and red (band4) for pure pixels of the soil, respectively". For this, the model had two phases-in the first phase, the pure pixel vegetation and the bare soil properties were derived from the NDVI, and in the second phase, the emissivity was determined by integrating the vegetation cover proportion per pixel with the NDVI and the emissivity value error estimate. The proportion of mixed pixel NDVI (consisting of mixed pixels of vegetation and bare soil) was calculated in Equation (4) below Atmosphere 2020, 11, 762 7 of 28 "where, i v Pure pixel of vegetation NDVI; i g Pure pixel of soil NDVI; di NDVI error" The "emissivity and cavity effects" was estimated as mentioned below, using Equations (5) and (6) where ∆ε was calculated from Equation (6) or, on a simplified form: "where, P v vegetation cover proportion, P s fraction of soil cover, ε v emissivity vegetation canopy, ε s emissivity bare soil, F, G,F shape factors considered the ratio of radiation from the vegetation to the ground, the ratio from the ground to the side and the ratio from the other side [70], respectively".
The dε term refers to the radiation that enters the sensor indirectly through internal reflections between the walls and the rugged terrain because the natural surfaces are not black.
Such an impact meant that the radiation that the sensor detected was higher than that emitted directly by the roughness components called the cavity effect of a rough surface [70]. The cavity effects value estimation of the term "dε is the mean weighted value" was considered to be the emissivity value of different land covers. "Weighted value dε of 0.04 was adopted in this study". The emissivity values calculated the term. "ε v and ε s , i.e., pure vegetation emissivity values and soil pixels were considered to be 0.978 and 0.914", respectively, from the literature referred to in "Valor and Caselles [15]".
to achieve accurate surface temperatures, the brightness temperatures with emissivity and other atmospheric parameters need to be corrected. The spectral radiance of thermal bands at the sensor is stated in Equation (8) where, LS j = "spectral radiance observed by the sensor", ε j = "surface emissivity at wavelength j", L BB j (T) = "spectral radiance from a blackbody at surface temperature T", L j sky = "spectral radiance incident upon the surface from the atmosphere (downwelling), from MODTRAN", L j atm = "spectral radiance emitted by the atmosphere (upwelling), from MODTRAN" τ j = "spectral atmospheric transmission, from MODTRAN".
Using the MODTRAN radiative transfer model, the radiance data for atmospheric effects was adjusted to achieve the radiance produced by the surface (L j ).
The basic data files were provided as inputs to perform the MODTRAN 4, which were calculated for the parameters required to further the investigation. The standard atmospheric parameter (tropical Atmosphere 2020, 11, 762 8 of 28 climate) was taken from this analysis. The resulting parameters were used to measure the radiance using the equation below: where, L BB j (T) = Hence, where, C 1 = "First radiation constant = 3.74151 × 10 −16 (W m 2 )" C 2 = "Second radiation constant = 1.44 × 10 4 (µm K)" λ j = "wavelength of channel j, (m)" If the surface emissivity was known, the reflected sky radiation could be corrected using the aforementioned equation. The surface temperature could then be determined using Equation (11):

LST Result from Emissivity Derived from the Proportion of Vegetation Cover in Conjunction with NDVI
It is reasonable to use ASTER's VNIR data to estimate vegetation cover and surface emissivity, regardless of TIR bands, due to differences in the scaling of low and high MMD emissivity spectra in TES. Figure 3 shows the LST spatial distribution retrieval from emissivity derived from the proportion of vegetation cover, in conjunction with NDVI. The result obtained from the satellite-derived LST of 2019, ranged from 22.00 • C to 55.16 • C, with an average of 44.01 • C and a standard deviation of 3.37. A field survey was conducted from 24-27 June 2019, comparing the estimated satellite surface temperature values with ground-based measurements (in-situ). All in-situ temperature measurements were taken for various surface areas and were measured at a typical distance of 1-2 m, based on the ease of access of the feature in a given place. Temperature measurements were carried out for dense vegetation (tree cover), exposed rocky area, and the concrete material, that was more dominant in the area in Table 1, illustrates vegetation cover temperature of 29.46 • C for field observation and 31.26 • C for satellite measurement. In contrast, concrete (asphalt: parking area) and concrete built-up field observation indicated a temperature of 50.42 • C and 41.12 • C and satellite measurement indicated a temperature of 48.77 • C and 43.29 • C, respectively. The error in these two surface temperature measurement methods was within the range of 1-3 • C. The coefficient of correlation between the derived satellite and the observed surface temperature field was very high, at 0.98 (significant at 0.01 level). Considering the high atmospheric attenuation in semi-arid mountainous urban climate, it was a reasonable value. Therefore, the comparison between satellite-driven LST and ground-based observed LST was in good agreement with the relative surface temperature analysis of the thermal environmental studies [75,76].  The LST was strong in north and east and low in west and south, showing a relatively steady decline from northeast to southwest. Additionally, the LST created a high-value area in downtown Abha-Khamis-Mushyet cities and other land-intensive building areas, showing a daily high-temperature (heat island) geographical phenomenon. In the Eastern part of the study area, i.e., the Khamis-Mushyet town, the LST formed a high-value area whereas, in the western part, i.e., Abha City, LST was low-value. The combination of elevation, NDVI, and LULC data were used to overlay analysis, which showed that the trend was mainly attributable to the cumulative effects of LULC, NDVI, and elevation. Moreover, the integration of high-resolution thermal maps revealed high-temperature zones distinct from the city's typical urban thermal environment. Apart from this, the highest LST was noted in the northern study area, due to exposed rocky area.

Relationship between LST and LULC
LST was ranked as follows, among different types of LULC-exposed rock > built-up land > bare soil > scrubland > agricultural cropland > sparse vegetation > dense water bodies ( Table 2). LST of sparse vegetation and built-up land was high, including its standard deviation. Generally, the LST of built-up, bare soil, and agricultural cropland with high a anthropogenic influence was higher than  The LST was strong in north and east and low in west and south, showing a relatively steady decline from northeast to southwest. Additionally, the LST created a high-value area in downtown Abha-Khamis-Mushyet cities and other land-intensive building areas, showing a daily high-temperature (heat island) geographical phenomenon. In the Eastern part of the study area, i.e., the Khamis-Mushyet town, the LST formed a high-value area whereas, in the western part, i.e., Abha City, LST was low-value. The combination of elevation, NDVI, and LULC data were used to overlay analysis, which showed that the trend was mainly attributable to the cumulative effects of LULC, NDVI, and elevation. Moreover, the integration of high-resolution thermal maps revealed high-temperature zones distinct from the city's typical urban thermal environment. Apart from this, the highest LST was noted in the northern study area, due to exposed rocky area.

Relationship between LST and LULC
LST was ranked as follows, among different types of LULC-exposed rock > built-up land > bare soil > scrubland > agricultural cropland > sparse vegetation > dense water bodies ( Table 2). LST of sparse vegetation and built-up land was high, including its standard deviation. Generally, the LST of built-up, bare soil, and agricultural cropland with high a anthropogenic influence was higher than other LULC classes; the LST of the water body and, dense vegetation that had low anthropogenic influence was significantly low. Hence, LULC classes should be planned realistically, and the effect of coolness should be encouraged by water and vegetation. Table 2 presents the results of the post-hoc test related to the LST and NDVI variations between various classes of LULC in the study area. Agricultural cropland was predominantly dry land with limited plant density, while sparse vegetation was mostly dry soil, Acacia woodlands, i.e., Acacia origena and A. gerrardii and low vegetation coverage. Agricultural cropland and sparse vegetation coverage showed only a substantial difference in surface emissivity, with no significant difference in LSTs; these influences could also be seen in the mean NDVI value of agricultural land and sparse vegetation. Influenced by differing heat flux, thermal conductivity and thermal radiation properties, the surface temperature of dense vegetation, water bodies, exposed and built-up rock showed significant variations, indicating that the various forms of LULC contributed significantly to the local thermal environment (Figure 3). Analysis of surface temperature distribution over built-up land revealed that high impermeable surface densities and low vegetation cover led to adverse shifts in the study area's thermal environment.

Characteristics LST and NDVI Distribution
The statistics of surface temperature and NDVI by LU/LC types were obtained by superimposing LULC with surface temperature, and NDVI. Figure 4a,c infer that the NDVI and LST displayed patterns of inverse spatial distribution. The LST showed opposite patterns in two separate ways, with the corresponding NDVI. The histogram profile of NDVI and LST over different LULC is shown in Figure 5. The surface temperature in Abha-Khamis-Mushyet twin cities developed a high-value area, showing an urban thermal pattern phenomenon that is exceptional for exposed rocky areas; further, the low-value NDVI zones in this region exhibited a small impact. The high LST peaks were related to the low NDVI value; meanwhile, the high NDVI peak areas, corresponded exactly to the low LST value.  On the contrary, the low LST and NDVI water body indicated its effective cooling effect. Compared to the built-up areas, the vegetation has "shown a significantly lower surface temperature, as dense vegetation can reduce the heat stored and surface structures by transpiration and minimize the soil" background. The LST of water bodies appeared to be low because the water slowly warmed up during the day, "possessing high thermal inertia and convection and turbulence" (for example, wave action). The water temperature was generally consistent with a vegetable coverage of about zero, and the water LST correlated positively with NDVI [77]. Meanwhile, the main part of the bare soil and the exposed rocky area was the absence of vegetation coverage.
Consequently, in the discussion of the linear coefficient relationship between NDVI and LST, between various LULC, the bare soil, exposed rock area, and water bodies were not taken into account. Figure 6 shows the linear regression coefficient between LST and NDVI for the LULC class as a whole, dense vegetation, agricultural land, and sparse vegetation. The most significant "coefficient of regression between LST and NDVI values found in dense vegetation (R 2 = −0.534) tailed agricultural cropland vegetation (R 2 = −0.494) and sparse vegetation (R 2 = −0.452)". In dense vegetation, agricultural croplands and sparse vegetation, with high to moderate correlation was observed between NDVI with surface temperature, which could use linear regression to simulate summer surface temperature, if the NDVI values were known ( Figure 6). The surface temperature could, therefore, be precise and reasonably simulated using the NDVI value during the summer or lower precipitation periods. Figure 4 illustrates that the LSTs of the various LULC classes had a negative linear regression coefficient with the NDVI value. In the linear regression coefficient, the values were ranked as follows-dense vegetation > agricultural cropland > sparse vegetation > scrubland, representing that the surface temperature of the dense vegetation decreased rapidly with On the contrary, the low LST and NDVI water body indicated its effective cooling effect. Compared to the built-up areas, the vegetation has "shown a significantly lower surface temperature, as dense vegetation can reduce the heat stored and surface structures by transpiration and minimize the soil" background. The LST of water bodies appeared to be low because the water slowly warmed up during the day, "possessing high thermal inertia and convection and turbulence" (for example, wave action). The water temperature was generally consistent with a vegetable coverage of about zero, and the water LST correlated positively with NDVI [77]. Meanwhile, the main part of the bare soil and the exposed rocky area was the absence of vegetation coverage.
Consequently, in the discussion of the linear coefficient relationship between NDVI and LST, between various LULC, the bare soil, exposed rock area, and water bodies were not taken into account. Figure 6 shows the linear regression coefficient between LST and NDVI for the LULC class as a whole, dense vegetation, agricultural land, and sparse vegetation. The most significant "coefficient of regression between LST and NDVI values found in dense vegetation (R 2 = −0.534) tailed agricultural cropland vegetation (R 2 = −0.494) and sparse vegetation (R 2 = −0.452)". In dense vegetation, agricultural croplands and sparse vegetation, with high to moderate correlation was observed between NDVI with surface temperature, which could use linear regression to simulate summer surface temperature, if the NDVI values were known ( Figure 6). The surface temperature could, therefore, be precise and reasonably simulated using the NDVI value during the summer or lower precipitation periods. Figure 4 illustrates that the LSTs of the various LULC classes had a negative linear regression coefficient with the NDVI value. In the linear regression coefficient, the values were ranked as follows-dense vegetation > agricultural cropland > sparse vegetation > scrubland, representing that the surface temperature of the dense vegetation decreased rapidly with increased NDVI value, followed by agricultural cropland and sparse vegetation, whereas the scrubland surface temperature declined slower rate. Thus, the surface temperature of dense and sparse vegetation and agricultural cropland was more influenced by vegetation cover than the scrubland and built-up land, and its cooling impact was more apparent than that of other types of LULC. Consequently, an adequate increase in green areas in urban lands improved the thermal environment in the city.
Atmosphere 2020, 11, x FOR PEER REVIEW 13 of 26 increased NDVI value, followed by agricultural cropland and sparse vegetation, whereas the scrubland surface temperature declined slower rate. Thus, the surface temperature of dense and sparse vegetation and agricultural cropland was more influenced by vegetation cover than the scrubland and built-up land, and its cooling impact was more apparent than that of other types of LULC. Consequently, an adequate increase in green areas in urban lands improved the thermal environment in the city.

Spatial Structure of LST and LULC Relationship
Due to the unusual occurrence of single land use, the LST-NDVI relationship should be studied in a spatial combination of various forms of LULC. In figure 6, the LST-NDVI feature space represents a unique triangle shape with an obtuse-angled ABC. A represents a low LST water body and a relatively negative NDVI value; B represents bare soils, exposed rock and built-up with high-LST and low-NDVI value; and C represents dense vegetation with the lowest-LST value and highest-NDVI value. In Figure 6, the different types coincided with each other in the LST-NDVI feature space; built-up, mostly mixed with exposed rock and bare soil had high evaporation and low moisture content; sparse vegetation was spread in the transitional zone between the agricultural cropland and dense vegetation, with high vegetation cover, relatively abundant soil moisture and low evapotranspiration. The dense vegetation's main distribution was localized mainly and followed Wadi from north-east to south-west (Figure 7).

Spatial Structure of LST and LULC Relationship
Due to the unusual occurrence of single land use, the LST-NDVI relationship should be studied in a spatial combination of various forms of LULC. In Figure 6, the LST-NDVI feature space represents a unique triangle shape with an obtuse-angled ABC. A represents a low LST water body and a relatively negative NDVI value; B represents bare soils, exposed rock and built-up with high-LST and low-NDVI value; and C represents dense vegetation with the lowest-LST value and highest-NDVI value. In Figure 6, the different types coincided with each other in the LST-NDVI feature space; built-up, mostly mixed with exposed rock and bare soil had high evaporation and low moisture content; sparse vegetation was spread in the transitional zone between the agricultural cropland and dense vegetation, with high vegetation cover, relatively abundant soil moisture and low evapotranspiration. The dense vegetation's main distribution was localized mainly and followed Wadi from north-east to south-west (Figure 7).

Relationship between LST and Topographical Parameters
The topographical impacts on surface temperature were examined in-depth, in terms of elevation, slope, and aspect with a basic understanding of surface temperature change patterns, as shown in Figures 4 and 6. The surface temperature and altitude spatial distribution exhibited an inverse pattern. With increasing altitude, the surface temperature gradually decreased. Topographical distribution alternated between the valleys and hills in east-west and northeast-southwest directions. Accordingly, the high-value peaks and low-value valleys of the surface temperature were distributed alternately. Moreover, the surface temperature had a notably vertical variation rule in the region, with a relatively high altitude difference. Taking Alsouda mountain (located in western part; elevation 2730 m MSL) and Khamish-Mushayet city (located in the eastern part) as examples, Figure 4d shows a low altitude in the Khamish-Mushayet area (Elevation 2000 m MSL), where there is a high surface temperature value. Figure 7 illustrates that the surface temperature of high hills gradually declined with the increase in elevation.

Effect of Elevation on the LST
To determine the effect of elevation on the surface temperature over different LULC, the relationship between the LST and the elevation was investigated through zonal statistical analysis, using ArcGIS software (Table 3) and scatter plot (Figure 7c), which showed the effect of elevation on the LST, similar to the effect of elevation in the air temperature [78]. Table 3 shows the statistical analysis LST over the various LULC, at the classified elevation zone. The elevation effect on the LST was clearly expressed over various LULC. The built-up land, the mean NDVI and LST at the elevation of 1564-2000 m was 0.061 and 46.14 • C, respectively, whereas at the elevation of 2501-2736 m, NDVI increase gradually and LST decrease gradually to 0.123 and 39.29 • C. A similar pattern of a gradual decrease in LST was found in the overall LULC classes because of the holistic impacts from surface type, solar radiation, meteorological conditions, the thermal property of the underlying materials, etc. The major gradual decrease in inclination was found over the densely vegetated area. However, statistical analysis of LST and elevation revealed that with a rise in elevation, the overall trend of LST declined. Due to the strong interactive effects from different factors, the relationship was thus more complex.

Effect of Aspect on the LST
In mountainous regions, the slope and aspect are known to affect radiation conditions. Different aspects alter the conditions of sun shading and cause the thermal interactions of the land surfaces. To strengthen the knowledge of the topographic effect on LST, the changing LST trends were assessed in eight cardinal directions "(North  (Figure 4e). Zonal statistics of the elevation zone was performed on each direction using NDVI and LST values to obtain the corresponding statistics. As depicted in Table 4, analysis findings from the variances in the cardinal directions was the highest in the east (67.5-112.5)  The statistical analysis of aspect with NDVI also drew some valuable finding. The lower NDVI value was found in the east (67.5-112.5) and the southeast (112. 5-157.5), at an elevation of 2501-2736 m. This could be attributed to higher temperatures and shifts in precipitation trends over the past three decades, which could account for the observed rise in junipers' drought stress at these altitudes [79]. These situations seem to suggest that the influence of aspect difference on LST should be a strong feature for various locations with different seasons. As for the process, the difference in surface solar radiation in each cardinal direction would cause directional changes. Slopes facing the sun receive more solar radiation, allowing faster heating than the slopes away from sunlight. Additionally, the LST of slopes facing the sun was higher than the LST of slopes faced away from the sun.

Effect of Slope on the LST
Besides the aspect impact, the slope factor also significantly regulated the quantity of direct radiation exposure that is commonly more global radiation percentage than the sky diffuse radiation, through both the angle between the normal slope and the solar beam. This might be why incoming solar radiation on earth's surface was quite susceptible to the slope and influenced the heating process of the surface. The image acquisition time for the Abha-Khamish-Mushyet area was 10:49:09 a.m., and the research area's mountain path was southeast for the sunlit slope and the shady north-west slope (Figure 4f). In the analysis, the average sunlit slope LST was 2.32 • C above the shady slope. To analyze the slope impact separately, we categorized the local solar zenith angle into five zones, as per the natural break method, and the zonal statistics of the elevation zone was performed at each slope range, using the NDVI and LST values, to obtain the corresponding statistics (Table 5). A specific phenomenon could be identified from the LST value presented in Table 5, which showed that there was a significant declining trend with a higher slope angle. Variations in LST values at various angles could range from 1564-2000 to 2501-2736, approximately 3-4 • C, indicating an important effect of the surface slope on LST. Though under different slope conditions, the frequency of the influence changed and it was linked to the atmospheric impacts and types of LULC.

Effect of Vegetation on the LST
It was evident that changes in the LST values were highly linked with changes in the topographical conditions. However, the relationship was not unique, and the vegetation often played a crucial role in the changes to the LST. Sandholt et al. [80] reported that the vegetation not only regulated the LST through the surface energy balance, but also affected the sensor's visible bare soil and vegetation coverage and changes in the radiative temperature between the vegetation canopy and the soil influence on the LST. In addition, the thermal characteristics of the surface were also partially affected by the vegetation conditions. For example, the LST assessment of the study area was chosen to explore the impacts to examine the LST variations under different vegetation conditions. As seen in Tables 3 and 4, NDVI was used as an indicator of vegetation and categorized into different ranges. In most circumstances, the areas with lower NDVI value could be considered as the bare surfaces (Figure 4c). Though a declining trend could be observed with the increase in the NDVI value, and this trend indicated the vegetation's thermal impacts, including evaporative control, thermal properties, and surface radiative balance. With vegetation canopy increases, the LST decreased within a relatively low-temperature range. This decrease might be reflected in part from the triangle space investigated in many studies between LST and NDVI.

Spatial Characteristics of NDVI and LST Distribution
Examining surface temperature patterns, transect (profile) was drawn across the image, i.e., north-east and south-west, along the Abha Wadi shown in Figure 7. This transect was made to observe a relationship in LST-NDVI-Elevation. Such transects enabled us to understand the factors that influence the city's thermal environment, as well as the role of vegetation with topographical variation. From the spatial profile of transects, it could be seen that there were high peaks LST value (box I) related to a low 'peaks' of the NDVI. Meanwhile, the zone highlighted with boxes (II and III) depicted the high peaks of NDVI value, which corresponded exactly to the low LST 'valley'. Two boxes in the Abha Wadi profile, i.e., I and IV, were characterized by comparatively low NDVI and LST. The low LST and NDVI suggest its strong cooling effect, as seen in box IV. In addition, in the spatial profile of transects, LST correlated negatively with elevation and moderately with the low linear coefficient of R 2 = −0.323. In contrast, LST was negatively correlated with the NDVI and showed a high linear coefficient of R 2 = −0.534.

Discussion
In this study, LSE retrieved from NDVI-derived emissivity method and its integration into the TES algorithm was applied to ASTER multi-spectral thermal data to accurately estimate LST. In the TES algorithm, surface emissivities were determined based on TIR image data alone, which led to scaling errors in the estimation of emissivity in heterogeneous surfaces, especially in the vegetated areas. The TES algorithm problem MMD threshold was where the target emissivities were close to the threshold value of 0.03. It related to the scaling errors that ultimately resulted in a scaling-value jump that induced discontinuity in the image of a transition from bare soil to the vegetation. Due to variations in the scaling of high-low MMD emissivity spectra in TES, it was appropriate to use ASTER's VNIR data to estimate vegetation proportion cover and surface emissivity, irrespective of TIR data. The results showed that emissivity derived from the NDVI method was more reliant on the upwelling, downwelling, and transmittance and would provide the best results [74]. The main strength of emissivity derived from the NDVI algorithms was the potential within the algorithms to incorporate atmospheric correction [81].
There are various advantages and disadvantages in contrasting the emissivity approach derived from NDVI with the TES methods. The TES algorithm provided a significant accuracy for LSE. However, the TES accuracy for surfaces with low MMD values was expected to decrease, as was the case with vegetated areas [74]. The results from the NDVI approaches were predicted to be marginally better in this case of LSE. Nevertheless, the emissivity-derived NDVI approach did not extend to surfaces with negative NDVI values, i.e., water, the emissivity spectra were also well-known for all surfaces. The efficiency of the TES method over low MMD surfaces need further studying and validation, besides adapting the NDVI method to rock covers. The ability to retrieve LST, together with LSE, was one advantage of the emissivity derived from the NDVI approach [74,81]. The emissivity method derived from NDVI offered improvements in surface emissivity prediction from visible and NIR data. This point should be considered when evaluating the NDVI-derived emissivity method; it referred to the estimation of the vegetation proportion cover from NDVI values (Equation (2)), according to Carlson and Ripley [82]. This was a simple method of estimating the vegetation proportion cover and was sufficiently accurate for using the NDVI method. Furthermore, a multi-spectral analysis by endmembers might result in an improvement in the estimation of vegetation proportion; although the principal source of error in the NDVI methodology is the selection of soil emissivity.
An increase in vegetation areas is the main option discussed for local climate mitigation [53]. Vegetation in urban areas governs climate mainly through shading [54], carbon dioxide sequestration [55], and evapotranspiration [56]. The mitigation potential of urban vegetation requires additional research because vegetation and local climate are strongly associated, and local plant species should be preferred as a strong mitigation option; however, exotic trees are common in the urban context. When vegetation becomes denser, there is a general decreasing trend in the LST (exhibited by the increase in NDVI). The dense vegetation canopy effectively prevents incident solar radiation and improves the evaporative cooling effect of vegetation, which prevents a rise in surface temperature [23,[57][58][59]. Therefore, the surface temperature decreases with the increase in vegetation cover. The terrain impact evaluation on the LST indicates that the change in temperature pattern shows significant variations with respect to the main factors associated with elevation, vegetation cover, slope and aspect, consistent with He et al. [57] findings. The general geographic phenomenon was observed for LST spatio-temporal variation at different seasons or different latitudinal locations. The seasonal MODIS LST correlation analysis [51,78] not only confirmed the important impact of topography and land cover on LST variation but also offered insights into the topography effect on LST by quantitatively examining temperature drop variability concerning various factors, slope, aspect, vegetation, and solar radiation.
The urban thermal effect of the built-up and densely populated area was notable in Abha-Khamis-Mushayet cities during the summer season (June). As such, apart from those induced by the classical urban thermal effect, several zones of high LST with various sizes and shapes were noticed in the study area's northern and central portion. From the LST distribution, the highest LST of more than 53 • C was observed where no construction was found. Abnormal regions of high LST were located mostly in the northern and central sections of the study region, with no vegetation cover and exposed rocky areas (Figure 4b). Bare soil contributed to low vegetation cover [83], surface water infiltration, loss of water [84], human intrusion, and degradation causes extreme soil erosion [20], resulting in exposed rock exposure and rocky desertification [85]. Therefore, surface substances' specific heat capacity and thermal inertia were low. Consequently, the LST area grew rapidly during the day and was higher in the surrounding area, producing an abnormal zone of high LST. Thereby, the LST of the exposed rocky areas characterized not only by the classical urban thermal effect but also the abnormal zone of high LST, exacerbated by the desertification in the rocky area. In future studies, LST could be used as a critical element and index to examine and assess the degree of desertification in rocky areas.
Twin cities Abha-Khamis-Mushayet have prominent mountain features, with mountains representing a more significant portion of the study area. Figure 5D indicates a tilting of the terrain from southwest to northeast, while a decrease in LST from northeast to southwest stretched along the mountains, suggesting a dominant role in the LST of the region by altitude, and there was a negative correlation between LST and elevation in the quantitative relation. Two principal factors were responsible for this pattern. (1) Surface temperature influenced by air-temperature, was also decreased with height. (2) In semi-arid mountainous regions, Abha-Khamis-Mushayet twin cities, the vegetation cover on the southwest side was relatively good compared to the north and the northeast. So the solar radiation that the ground surface received was often spread in latent heat, and low LST. In the meantime, the intensity and length of solar radiation varied in various slope directions. In this analysis, the LST also differed across multiple slope directions, consistent with the Wen et al. [86] findings. Different types of LULC show different thermal conductivity, thermal capacity, roughness, and albedo, resulting in LST variations. In the analysis, the LSTs of different LULC types differed in the same direction at the same elevation. In overview, the LST in this present study area was affected by altitude, slope aspect, and LULC type, in which elevation and vegetation played a key role in the overall LST pattern.
Current literature on the LST-NDVI relationship varied substantially in combinations of different LULC classes (Table 6). In this research, the LST and NDVI scatter plot diagram showed an obtuse-angled triangle distribution in combinations of different LULC classes, consistent with the Zhou et al. [87] findings. LST-NDVI showed a negative correlation and NDVI scatter plots showed an "obtuse-angled triangle" distribution Still, it varied from what was found by other researchers, due to variations between quantitative and qualitative analyses or variations in soil backgrounds and the sensor used in the studied area. The LST-NDVI feature space showed a unique ABC obtuse-angled triangle shape an overall showed a negative linear correlation. Clearly, A characterized a water body with low-LST and somewhat negative NDVI value; B characterized bare soils, exposed rock, and built-up with high-LST and low NDVI value; and C characterized dense vegetation with the lowest LST value and the highest NDVI value, which was compatible with the findings of Zhou et al. [87]. Nevertheless, it should be noted that the analysis still had limitations. Soil moisture should be a very significant factor affecting LST's evolution. Several studies showed their indirect relationship and established methods to infer soil surface moisture, based on LST variation [28,44,76]. Consequently, in future research, the surface moisture effect should be known to some degree.

Conclusions
In the present study, LSE was retrieved from the proportion of vegetation cover in conjunction with NDVI and integrated into the TES algorithm to derive LST from ASTER data. The spatial distribution patterns of LST and factors of influence in the semi-arid mountain cities of Abha-Khamish-Mushayet were summarized by combining LULC classes and NDVI with DEM data. The LST variations between various LULC classes were then investigated and the characteristics of the distribution of NDVI and LST in the study area and their quantitative relationship were discussed. The main contributions of this research were described (1) retrieving LST from the proportion of vegetation cover in conjunction with NDVI integrated into the TES algorithm and comparing the estimated satellite surface temperature values with ground-based measurements (in-situ), which was consistent with the actual value. Influenced by topographical parameters, i.e., elevation, slope, aspect, and underlying types of LULC, LST exhibited an overall decreasing trend from northeast to southwest and stretched across the extent of the mountain, suggesting a role of LST by altitude. The imaging time, i.e., summer, abnormal regions of high-temperature were located mainly in the northern and central sections of the study region, with no vegetation cover and exposed rocky areas. Bare soil contributed to low vegetation cover, surface water infiltration, loss of water, human intrusion, and degradation caused extreme soil erosion, resulting in exposed rock exposure and rocky desertification. Therefore, surface substances' specific heat capacity and thermal inertia were low. Consequently, the LST area grew rapidly during the day and was higher in the surrounding area, producing an abnormal zone of high LST. (2) Several comparative analysis indicated that the difference in LST was important among most types of LULC; the exposed rock area exhabited the highest mean LST of 45.67 • C, followed by built-up land of 44.44 • C, and water bodies and dense vegetation exhibited the lowest LST of 25.75 • C and 36.88 • C, respectively. LSTs were ranked as follows, among different types of land use/land cover: exposed rock > built-up > bare soil > scrubland > agricultural > sparse vegetation > dense water bodies. Generally, the LST of built-up, bare soil, and agricultural cropland with high anthropogenic influence intensity was higher than other LULC classes; the LST of the water body and, dense vegetation, which had low anthropogenic influence intensity, and was significantly low. Hence, LULC classes should be planned realistically, and the effect of coolness should be encouraged by water and vegetation. (3) For the LULC type, water bodies, bare soil, and exposed rocky area were not considered in the discussion on the linear coefficient relationship between NDVI and LST, as it had a nonlinear relationship. The most significant "coefficient of regression between the LST and NDVI values were found in dense vegetation (R 2 = −0.534), tailed agricultural cropland vegetation (R 2 = −0.494), sparse vegetation (R 2 = −0.452)". In dense vegetation, agricultural croplands and sparse vegetation, with a high to moderate correlation was observed between NDVI with surface temperature, which could use linear regression to simulate summer surface temperature, if the NDVI values were known. The surface temperature could, therefore, be precise and reasonably simulated using the NDVI value during the summer or lower precipitation periods. (4) The LST-NDVI feature space showed a unique ABC obtuse-angled triangle shape. Conclusively, A represented a water body with low LST and relatively negative NDVI value; B represented bare soils, exposed rock, and built-up with high LST and low NDVI value; and C represented dense vegetation with the lowest LST value and the highest corresponding NDVI value. (5) The vegetation cover on the southwest side was relatively good compared to the north and the northeast. So the solar radiation that the ground surface received often spreads in latent heat, and low LST. In the meantime, the solar radiation length and intensity varied in various directions of the slope. The LST also varied in different slope directions. Different types of LULC had varied thermal conductivity, thermal capacity, roughness, and albedo, resulting in LST variations. In the analysis, the LSTs of different LULC types differed in the same direction at the same elevation. The LST was affected by elevation, slope aspect, and LULC class, in which elevation and vegetation played a key role in the overall LST pattern. This study offers a framework for land-use planning and environmental protection in urban mountainous areas.