Analysis of Near-Surface Temperature Lapse Rates in Mountain Ecosystems of Northern Mexico Using Landsat-8 Satellite Images and ECOSTRESS

: Mountain ecosystems provide environmental goods, which can be threatened by climate change. Near-Surface Temperature Lapse Rate (NSTLR) is an essential factor used for thermal and hydrological analysis in mountain ecosystems. The aims of the present study were to estimate NSTLR and to identify its relationship with aspect, Local solar zenith angle (LSZA) and Evaporative Stress Index (ESI) for two seasons of the year in a mountain ecosystem at the North of Mexico. Normalized Land Surface Temperature (NLST) was estimated using environmental and topographical variables. LSZA was calculated from slope to consider the effect of solar position. NSTLR was estimated through simple linear models. Observed NSTLR was 9.4 ◦ C km − 1 for the winter and 14.3 ◦ C km − 1 for the summer. Our results showed variation in NSTLR by season. In addition, aspect, LSZA and ESI also inﬂuenced NSTLR regulation. In addition, Northwest and West aspects exhibited the highest NSTLR. LSZA angles closest to 90 ◦ were related with a decrease in NSTLR for both seasons. Finally, ESI values associated with less evaporative stress were related to lower NSTLR. These results suggest potential of Landsat-8 LST and ECOSTRESS ESI to capture interactions of temperature, topography, and water stress in complex ecosystems.


Introduction
Mountain ecosystems are characterized by being located at high elevations and showing environmental heterogeneity within their altitudinal gradient [1,2]. Mountain ecosystems provide a large quantity of ecosystem services [3,4]. One of most important ecosystem services is water supply, it is estimated more than half of water resources for human consumption are directly or indirectly from these areas [5][6][7]. In addition, forested areas and grasslands within the altitudinal gradient of mountain ecosystems are important for the processes of climatic regulation, net primary productivity, and soil conservation [5,[7][8][9][10].
These ecosystems and the environmental processes occurring within them are very sensitive to climate change and other impacts from human activities, such as deforestation and overgrazing [8,11]. Globally, mountain ecosystems are threatened by human activities, including land use change degradation and human-caused climate change, which is affecting forests distribution [12,13]. For examples, changes in forest and glacier distribution related to climate change have been documented in Himalayan ecosystems [13,14]. These processes modify the humidity, precipitation and temperature patterns throughout the altitudinal gradient and affect water availability, endangering water reserves for human population, and affecting fauna and flora [14,15]. In spite of the importance and fragility of these ecosystems, research on climatic processes in mountain areas has been limited by accessibility constraints [11,16].
The analysis of mountain ecosystems trough a thermal approach could provide information about environmental dynamics of these ecosystems. An essential variable to studied thermal and climatic conditions in mountainous area is near-surface temperature lapse rate (NSTLR), defined as the rate of change in temperature by elevation [17]. NSTLR is defined as inverse relationship between elevation and temperature [17][18][19]. It is generally assumed to have a value of 6.5 • C km −1 , inferred from the NSTLR of the atmospheric column within the troposphere [20,21], and is generally the constant value of NSTLR in mountain areas. Nevertheless, several studies have concluded that the NSTLR is locally variable, depending on the region and season, with varying results of the thermal gradient value and the factors affecting it [18,[22][23][24]. There is consequently a need to further analyse the factors influencing the seasonal and spatial variability in NSTLR.
NSTLR can be affected by topographic variables, such as aspect and slope and other variables such as evapotranspiration (ET) [1,25], although many questions remain regarding the environmental drivers of NSTLR. The characterization of NSTLR in mountain ecosystems is limited by the availability and spatial distribution of weather stations, often limiting a spatially detailed characterization of climatic data [24,26,27].
Higher resolution DEM and satellite climate information offers a new opportunity to better analyze thermal and hydric processes in mountain ecosystems [13,28]. Remote sensing techniques utilizing satellite images from passive remote sensors can be a useful alternative to study mountain ecosystems, since they can provide information in areas without water stations data and areas under low accessibility conditions, providing historical data at regional and global scales [29]. The thermal dynamics can be monitored from land surface temperature (LST) data estimated from satellite images, which provide spectral information about the processes of energy exchange between solar radiation and the terrestrial surface [30].
Several studies have been utilized satellite information to estimate and analyze the NSTLR in mountain ecosystems, using the relationship between LST and a digital elevation model (DEM) [8,18,26,[31][32][33]. Nevertheless, it has been observed that some environmental variables (i.e., solar incidence, moisture, land cover characteristics) cause interference of NSTLR estimated from LST data [18,34]. Therefore, it has been suggested normalizing LST data with respect to environmental variables, potentially improving the accuracy of estimation of NSTLR from satellite information [8,35]. This approach has been less explored and requires further evaluation, particularly for the relatively less utilized medium resolution sensors. In particular, studies using Landsat Land Surface Temperature product to analyse factors driving thermal gradients in forest ecosystems are still relatively scarce [18,26,36]. In the present study, we used LST derived from Landsat-8 TIRS, which was normalized using random forests to obtain NSTLR information. This information allowed us to analyze the variability of NSTLR with respect to biophysical parameters, such as evapotranspiration (ET) and topographic characteristics.
The knowledge about NSTLR and the factors influencing it at a regional scale can allow to improve the understanding the thermal and hydric dynamics of mountain ecosystems [18] and could be used to develop climatic models, hydrological models and analysis of soil moisture [37]. For example, the effect of slope which depends on solar radiation flux incidence on the terrestrial surface; it can be quantified through the Local Solar Zenith Angle (LSZA). This parameter is calculated considering the interaction of the slope relative to sun position at time of satellite pass [38].
In addition to the previously documented contrasting effects of vegetation and topography on NSTLR [18,26], other parameters potentially affecting thermal gradients have received less attention. For example, very few studies have analysed the variation of ther-Remote Sens. 2022, 14, 162 3 of 16 mal gradients with respect to remotely sensed evapotranspiration (ET) indices, particularly from medium resolution sensors.
ET is a key hydrologic component for the energy balance of terrestrial surface, influencing water availability [39]. Until recently, there remained a large gap in our ability to monitor ET concurrently at both fine spatial and fine temporal scales globally [40]. The relatively recently launched (June 2018) ECOSTRESS sensor provides a good opportunity for monitoring daily ET at medium (70 m) resolution globally [40,41]. There remains a need to evaluate ECOSTRESS mission performance, particularly the ET estimates in humid, temperate, mountainous areas [41]. In addition, to our best knowledge, no previous study has analysed the variations Landsat-derived NSTLR against ECOSTRESS ET indices such as the Evaporative Stress Index (ESI).
Therefore, the aim of current study was to examine the variability of NSTLR for two contrasting year seasons (winter-summer) in a forest management area at the Sierra Madre Occidental (SMO) in Durango State, Mexico and to analyze the influence of topography and other biophysical factors such as evaporative stress, through Landsat-8 and ECOSTRESS images. The study had following hypotheses: (1) NSTLR in the study area is different than the constant value (6.5 • C km −1 ), (2) NSTLR is influenced by different topographic and biophysical characteristics (LSZA, aspect and ESI). Improved NSTLR characterization could contribute to ecological modeling and water management planning, including water balances assessment [31,42]. Furthermore, understanding the effects of topography on temperature gradients could contribute to ongoing fire risk analysis e.g., [43,44].

Study Area
The area of study was the management unit UMAFOR 1001, in the central part of the SMO, at the North of Durango State, Mexico ( Figure 1). The SMO is the highest mountain range in México, it has a one of the main forest reserves in Mexico, which gives it great economic and ecological relevance given the wide variety of environmental services it provides [45][46][47]. UMAFOR 1001 is bounded between 25.57 • and 26.84 • N latitude and 105.05 • and 106.55 • W longitude with an approximate 11,336 km 2 of area. The average elevation of study area is 2131 m above the mean sea level. The vegetation includes more than 75 species of pine and oak forests in the upper and middle elevations [48]. Among these species, Quercus sideroxyla Humb and Bonpl, Pinus durangensis Martínez, Quercus grisea Liebm and Pinus teocote Humb and Bonpl stand out [49], while mid-slopes and flat terrains are characterized by semiarid vegetation, such as shrublands, chaparral, natural grasslands and induced pasturelands [47]. Climate ranges from subhumid to semiarid, with an annual precipitation average of 664 mm and mean temperatures ranging from 8.7 (average of coldest month) to 16.1 • C (average of the warmest month), with an annual average of 11.1 • C [48,50]. UMAFOR 1001 has a variable topography ranging from undulating to rugged topography and altitudinal characteristics that make it valuable to analyze regional thermal processes.

Images Acquisition
The Digital Elevation Model (DEM) from the Shuttle Radar Topography Mission (SRTM), with a resolution of 30 m, was downloaded from the US Geological Service webpage (https://earthexplorer.usgs.gov) (accessed on 25 July 2020). Slope and aspect were obtained from the DEM, using QGIS 3.14 [51].
Google Earth Engine (GEE) platform was used to visualize, process, and download Landsat-8 OLI/TIRS (with a 30 m resolution) [52]. The GEE allows processing large datasets in a cloud-based platform, accessing satellite data such as Landsat collections. The images correspond to LANDSAT/LC08/C01/T1_SR collection from winter and summer of 2019 (Table 1). In addition, a cloud and cloud shadow mask from Quality assessment band (BQA) band was applied to remove pixels where attributes might be altered by climatic  [50]. Temperate forests include pine forests, oak forests, and mixed forests. The disturbance forest classification corresponds to forests that suffered some alteration and are under a natural regeneration process; (c) Location of study area.

Images Acquisition
The Digital Elevation Model (DEM) from the Shuttle Radar Topography Mission (SRTM), with a resolution of 30 m, was downloaded from the US Geological Service webpage (https://earthexplorer.usgs.gov) (accessed on 25 July 2020). Slope and aspect were obtained from the DEM, using QGIS 3.14 [51].
Google Earth Engine (GEE) platform was used to visualize, process, and download Landsat-8 OLI/TIRS (with a 30 m resolution) [52]. The GEE allows processing large datasets in a cloud-based platform, accessing satellite data such as Landsat collections. The images correspond to LANDSAT/LC08/C01/T1_SR collection from winter and summer of 2019 (Table 1). In addition, a cloud and cloud shadow mask from Quality assessment band (BQA) band was applied to remove pixels where attributes might be altered by climatic conditions; this process is a prerequisite to perform a more trustable thermal analysis; this process that must be considered every time that satellite images are processed [53,54].   [50]. Temperate forests include pine forests, oak forests, and mixed forests. The disturbance forest classification corresponds to forests that suffered some alteration and are under a natural regeneration process; (c) Location of study area.

Spectral Indices Estimation
The spectral indices had estimated employing Landsat-8 bands. These indices were estimated using the Equations (1)-(3) in the GEE platform and was used for the normalization LST.
Remote Sens. 2022, 14, 162 where NDVI is the Normalized difference vegetation index developed by Rouse, et al. [55], which is used for analyzing the photosynthetic activity in the terrestrial surface. NDVI was estimated using red band (B4) and Near infrared band (B5). Whereas NDBI is the Normalized difference built-up index (NDBI), this index allows emphasizing the impervious surfaces such as bare and built-up land [56]. NDBI was calculated with the Equation (2), using the Shortwave-infrared band (B6) and Near infrared (B5) data of Landsat-8. Finally, the Normalized difference water index (NDWI) provides wetness information, as soil moisture, vegetation, and built-up land [57]. The NSWI was calculated using B3 (green band) and B5 (Near infrared band).

Land Surface Temperature Estimation
LST was estimated using the inversion of Planck's function (Equation (4)). Where, TB is brightness temperature, defined as the temperature of blackbody. TB values come from band 10 coming from Google Earth Engine Landsat collection, λ is wavelength in radiance values, and ρ is a constant value of 14,380.
Land surface emissivity (ε) was calculated using the methodology of Sobrino and Raissouni [58], as shown in Equation (5), using the constants predefined according to the sensor (0.986 and 0.0004) and NDVI, NDVImin, NDVImax calculated in previously step.
The LST obtained was in Kelvin and it was converted to Celsius by subtracting 273.15. LST images was processed in the GEE platform using the code developed by Ermida, et al. [59] applying the methodology previously described.

Local Solar Zenith Angle
In order to analyze the effect of slope relative to sun position at the moment of satellite pass, LSZA, the angle between the zenith vertical and solar radiation (ranging from 0 • to 90 • ) was calculated. The local LSZA was calculated utilizing Equation (6), which was developed by Dozier and Frew [60].
cos Θ = cos Zs cos S + sin Zs sin S cos(As − A) where Θ is LSZA, Zs and As are the solar zenith angle and azimuth angle, respectively, and A and S are the aspect and slope of terrain, in degrees, obtained from the DEM.

Evaporative Stress Index
Evaporative stress index (ESI) is an ECOSTRESS product, which is a useful tool to analyze drought processes in areas where precipitation data are limited. It is based on evapotranspiration and potential evapotranspiration obtained from satellite thermal bands and some auxiliary products, i.e., vapor pressure data and relative humidity [61]. ESI can range between 0 and 1, where 0 is an area with high hydric stress, while 1 represents an area without hydric stress [39]. Likewise, ECOSTRESS products have a spatial resolution of 70 × 70 m [40] and its products include information for geometric correction and cloud detection [62]. Raster data was downloaded from the LP DAAC platform from US Geological Service (https://lpdaac.usgs.gov) (accessed on 19 August 2020). Images were downloaded for the UMAFOR 1001 polygon; the selected images were those which had the highest number of quality pixels for winter (24 February 2019) and summer (21 June 2019).

Near-Surface Temperature Lapse Rate Estimation
The NSTLR estimation trough the LST-DEM linear regression has a well-known error, which is due of effect of surface characteristics [8,63,64]. To solve this effect the LST obtained from Landsat-8 was normalized relative to biophysical characteristics and solar incidence angle. We used the approach developed for Firozjaei, Fathololoumi, Alavipanah, Kiavarz, Vaezi and Biswas [8], through a normalization of LST relative to the surface biophysical characteristics and the solar incidence angle (Equation (7)).
where NLST is LST modeled based on the surface biophysical characteristics and solar incidence angle, f is non-linear linking model. Random Forest regression method was used to build the non-linear linking model between LST obtained and the predictors (NDVI, NDBI, NDWI and LSZA). Finally, NSTLR was calculated by a simple linear regression between NLST (dependent variable) and DEM (independent variable), where regression slope represents the NSTLR value [8,33]. The raster layers (NLST, Aspect, LSZA and ESI) were aligned and aggregated to the same spatial context; this process was made with Bicubic interpolation raster algorithm. A database was constructed extracting information for every pixel, these processes were performed using QGIS 3.14 [51]. The NSTLR was analyzed separately to each aspect (i.e.,  (8)) for each one of the data ranges analyzed for aspect, slope and ESI ( Table 2). The goodness of fit of the 19 fitted linear models was evaluated by the coefficient of determination (R 2 ) and root square mean error (RMSE).
where b and m are coefficients of the linear regression between NLST and elevation (DEM); m is the slope value of linear regression, which represents the NSTLR. Finally, an overview of the workflow described in the previous sections is presented in Figure 2. Finally, an overview of the workflow described in the previous sections is presented in Figure 2.

Results.
The accuracy of Random Forest model to normalization LST data were to winter and summer were 0.71 R 2 and 0.86 R 2 , respectively. Meanwhile, the mean winter NLST of the study area ranged between 5.8 and 26.7 °C, with an average value of 17.4 °C. Mean summer LST ranged in the study area between 18.1 and 47.8 °C, with an average value of 34.3 °C. ESI values ranged for the study area between 0.03 and 0.7 in winter and between 0.3 and 0.9 in summer ( Table 2).
NLST maps for winter and summer show the LST distribution in the study area, where the LST decreasing in the zones with highest elevation, while in the lower areas was found higher temperature ( Figure 3). Linear regression between LST and elevation for winter had an R 2 of 0.74, with a RMSE of 1.9 °C; for summer, R 2 was 0.81 and RMSE 2.4 °C. Linear models showed that elevation influenced LST more in summer. The observed dispersion of data from the linear regression in both seasons ( Figure 4) suggests that there is a complex relationship where interacting multiple factors, in addition to elevation. The complementary factors which affect the NSTLR can be related to land cover, wetness of soil and amount of solar radiation, caused by slope and angle.

Results
The accuracy of Random Forest model to normalization LST data were to winter and summer were 0.71 R 2 and 0.86 R 2 , respectively. Meanwhile, the mean winter NLST of the study area ranged between 5.8 and 26.7 • C, with an average value of 17.4 • C. Mean summer LST ranged in the study area between 18.1 and 47.8 • C, with an average value of 34.3 • C. ESI values ranged for the study area between 0.03 and 0.7 in winter and between 0.3 and 0.9 in summer ( Table 2).
NLST maps for winter and summer show the LST distribution in the study area, where the LST decreasing in the zones with highest elevation, while in the lower areas was found higher temperature ( Figure 3). Linear regression between LST and elevation for winter had an R 2 of 0.74, with a RMSE of 1.9 • C; for summer, R 2 was 0.81 and RMSE 2.4 • C. Linear models showed that elevation influenced LST more in summer. The observed dispersion of data from the linear regression in both seasons ( Figure 4) suggests that there is a complex relationship where interacting multiple factors, in addition to elevation. The complementary factors which affect the NSTLR can be related to land cover, wetness of soil and amount of solar radiation, caused by slope and angle.
The NSTLR in the areas with different features of aspect, LSZA and ESI, was calculated the DEM-NLST lineal regression model separately and the results is showed in Table 3. The goodness of fit of the 19 fitted linear models shown R 2 above 0.5 in almost all variables. However, the differences in R 2 values suggest that ESI condition strongly affects the DEM-NLST relationship.

Temperature Lapse Rate
In our study area, observed NSTLR was 9.4 • C km −1 for winter and 14.3 • C km −1 for summer. Likewise, the NSTLR values was different in the different aspects, the summer ranged between 14.2 and 15.1, while summer shows NSTLR ranged between 9.2 and 10.8. The highest NSTLR was found for Northweast aspect and the lowest for Southeast aspect at both seasons of the year ( Figure 5). Remote Sens. 2022, 14, x FOR PEER REVIEW 8 of 17  The NSTLR in the areas with different features of aspect, LSZA and ESI, was calculated the DEM-NLST lineal regression model separately and the results is showed in Table  3. The goodness of fit of the 19 fitted linear models shown R 2 above 0.5 in almost all variables. However, the differences in R 2 values suggest that ESI condition strongly affects the DEM-NLST relationship.   The NSTLR in the areas with different features of aspect, LSZA and ESI, was calculated the DEM-NLST lineal regression model separately and the results is showed in Table  3. The goodness of fit of the 19 fitted linear models shown R 2 above 0.5 in almost all variables. However, the differences in R 2 values suggest that ESI condition strongly affects the DEM-NLST relationship.   We utilized the LSZA equation to quantify the effect of solar position on the ter slope at acquisition time, since the solar position changes the angle and intensity of solar energy above the terrain [65,66]. NSTLR decreased with increasing LSZA for b seasons, possibly caused by a less direct impact of sunlight at higher LSZA, wi difference between seasons of more than 5 °C km −1 (Figure 6). We utilized the LSZA equation to quantify the effect of solar position on the terrain slope at acquisition time, since the solar position changes the angle and intensity of the solar energy above the terrain [65,66]. NSTLR decreased with increasing LSZA for both seasons, possibly caused by a less direct impact of sunlight at higher LSZA, with a difference between seasons of more than 5 • C km −1 (Figure 6). We utilized the LSZA equation to quantify the effect of solar position on the ter slope at acquisition time, since the solar position changes the angle and intensity of solar energy above the terrain [65,66]. NSTLR decreased with increasing LSZA for b seasons, possibly caused by a less direct impact of sunlight at higher LSZA, with a di ence between seasons of more than 5 °C km −1 (Figure 6).  NSTLR against ESI ranges are shown in Figure 7. Observed ESI ranges were 0.2-1 for summer and 0-0.8 for winter (Figure 7). Observed NSTLR in our study was lower in areas with higher ESI for both seasons. These results suggest that areas with the higher water availability have lower thermal gradients and this effect is stronger in summer season.
Remote Sens. 2022, 14, x FOR PEER REVIEW 11 NSTLR against ESI ranges are shown in Figure 7. Observed ESI ranges were 0.2-1 summer and 0-0.8 for winter (Figure 7). Observed NSTLR in our study was lower in a with higher ESI for both seasons. These results suggest that areas with the higher w availability have lower thermal gradients and this effect is stronger in summer season

Discussion
In areas with a low availability of field weather stations, satellite data can provid alternative information, capturing a higher degree of variability caused by interaction climate with topography or evaporative processes, the latter influenced by land co properties. NSTLR have been shown to have a middle uncertainty when estimated f satellite data, nevertheless the methodology developed by Firozjaei et al. [8] allows mating NSTLR with more accuracy.

Discussion
In areas with a low availability of field weather stations, satellite data can provide an alternative information, capturing a higher degree of variability caused by interactions of climate with topography or evaporative processes, the latter influenced by land cover properties. NSTLR have been shown to have a middle uncertainty when estimated from satellite data, nevertheless the methodology developed by Firozjaei et al. [8] allows estimating NSTLR with more accuracy.
There is a clear difference of NLST between both seasons of year (winter and summer) explained by amount of energy received in summer, that is two times more than in winter [67]. In addition, biophysical and topographical features were key factors to NLST distribution, the effect of these variables seems to be stronger in winter. He at al. [18] noted that NLST was also influenced by aspect, slope, and solar radiation. On the order hand, the relative fine resolution scale of thermal satellite information might be evincing the topographic influences more accurately and cause the dispersion data observed in density plots of elevation-NLST.
The most widely accepted altitudinal NSTLR value is 6.5 • C km −1 . This value is given as a function of temperature in a vertical atmospheric column of wet air is frequently accepted as a NSTLR constant in mountainous environments [20,24,68]. This general value does not address the particularities of every study area and might be inexact in some specific locations [69]. For example, Soto Molina and Delgado Granados [70] found NSTLR to vary between locations in their analysis in Nevado de Toluca, Mexico, and consequently suggested that it should be quantified at each location. The seasonal and spatial variation observed in our study supports the importance to analyze NSTLR at a regional and season level. Our observed NSTLR were higher than the values previously reported by several studies [32,42,71]; this could be related to the influence of the land use types analysed on the thermal characteristics of the land Surface. Future studies could focus on analysing NSTLR variations against land use categories [72].
In the present study, the NSTLR was higher in summer, this agrees with higher NSTLR observed in summer in previous studies, such as Tang and Fang [27] study in a mountainous location in China or Navarro-Serrano, López-Moreno, Azorin-Molina, Alonso-González, Tomás-Burguera, Sanmiguel-Vallelado, Revuelto and Vicente-Serrano [71] analysis in the main mountain ranges in Spain. In contrast, Li et al. [73] observed a more marked NSTLR in winter season in Northern China mountain areas. This discrepancy could be related to our study area being dominated by semiarid ecosystems, with pasturelands in the lower slopes. It could also possibly be explained by the precipitation regimes characteristic of the area of study, receiving precipitation in summer season [74], therefore having lower evaporative stress in spite of the high temperatures. On the other hand, lower precipitation in winter resulted in lower water availability for evapotranspiration and a higher hydric stress. Conversely, studies reporting higher NSTLR in winter than in summer are generally associated with mountain ecosystems with higher humidity conditions [24,27,73,75]. These results confirm the utility of accounting for evaporative stress in understanding the seasonal and spatial variations of the thermal gradient.
More specifically, our study reveals the potential of the relatively recent ECOSTRESS sensor to characterize the impact of vegetation evaporative stress on thermal gradients from the Landsat sensor at relatively fine spatial scales. The high variation observed in summer NSTLR across ESI intervals, with a range from wetter areas with lower NSTLR, to areas with higher water stress and with higher NSTLR values, might be related to contrasting interactions of vegetation cover and water stress conditions. Areas with high vegetation cover limit sunlight incidence on terrain surface, cooling surfaces through evapotranspiration, while bare ground areas are associated to high surface temperatures; this might result in high NSTLR, mainly in summer season, where the solar radiation has higher effects [18,26]. While our findings agree with previous studies reporting effects of water availability on LST in environments with altitudinal variability (e.g., [1,27]), such previous analyses were conducted utilizing other water availability variables such as NDVI [1], unlike the ECOSTRESS ESI data utilized here.
The aspect was an important characteristic for the NSTLR values, it had higher values in the Northwest aspect of the two seasons studied (winter and summer). In the north hemisphere, southern and eastern aspects generally present a more homogeneous temperature through the altitudinal gradient. That agrees with, Liu and Li [76] in a study in China mountains where observed a higher NSTLR in North and West aspects, which are opposite to solar radiation direction. Aspect plays a key role in thermal processes in mountain areas; the higher NSTLR observed for N and NW aspects for both seasons might be related to a non-homogeneous heating process in aspects opposed to solar radiation, as opposed to a more homogeneous heating in aspects oriented towards solar radiation (i.e., S and SE) [18]. In the north hemisphere, southern and eastern aspects generally present a more homogeneous temperature through the altitudinal gradient.
Meanwhile, higher NSTLR in the 0-15 • LSZA range are very possibly a consequence of a higher solar radiation, which in turns affects moisture conditions [75], which also affects NSTLR (i.e., higher NSTLR under water stress conditions) (e.g., [77]). This supports the effect of solar angle and slope interactions on the thermal conditions in mountain ecosystems observed by several authors [18,24,69].
These results suggest sensitivity of NLST derived of Landsat-8 thermal information. Besides, ECOSTRESS product, ESI allows to analyze dynamics of vegetation water stress and temperature in our study area across an altitudinal gradient, possibly capturing variations between vegetation types and tree cover characteristics, a point that deserves further analysis. NSTLR analysis can help to understand the altitudinal migration of forest species that has been reported in several locations globally as a consequence of climate change [78][79][80][81]. Detailed NSTLR and ESI maps might aid in identifying fireprone areas [82,83]. Future studies might explore the inclusion of local relationships of topographic effects on temperature (e.g., [84]) to refine ongoing research of fire risk (e.g., [43,44,[85][86][87]) or biomass or tree growth modeling (e.g., [88,89]) in these diverse and complex ecosystems.

Conclusions
Near-surface temperature lapse rate (NSTLR) is one of the main features to understand climatic dynamics in mountain environments. Our results showed variation in NSTLR, quantified from NLST derived of Landsat-8 surface temperature and environmental variables, by season. Winter and summer NSTLR were different than the global average constant value widely utilized (6.5 • C km −1 ). NSTLR also varied widely between aspects, with the highest value observed for Northwest aspect in both seasons. In addition, Local solar zenith angle (LSZA) and evapotranspiration also influenced NSTLR regulation. The areas with greater LSZA showed lower NSTLR in both seasons. Evaporative Stress Index (ESI) showed a negative relation with NSTLR values, suggesting a role of water stress in increasing the NSTLR, this being more marked in summer season.
The results from the current study suggest potential of Landsat-8 LST and ECOSTRESS ESI to capture the interactions of topographic factors and vegetative stress on NSTLR gradients at relatively fine spatial scales, which can aid in water management planning and conservation planning. Further factors, such as vegetation types and tree cover are possibly influencing these water stress-topography-temperature gradient interactions, a point that should be further analyzed in future studies. Furthermore, future studies might focus on a multitemporal analysis of thermal gradients, including further quantifying topography and water stress interactions with vegetation and fuel characteristics, which might be beneficial for enhancing fire risk or tree growth predictions in these complex ecosystems.