Introduction

Human activities change weather patterns at global and local scales (^{IPCC, 2018}), causing cities experience consistently higher temperatures than surrounding rural areas (^{Oke, 1982}; ^{Howard, 1988}; ^{Eniolu et al., 2013}). The Urban Heat Island (UHI) effect is a pressing policy issue, because of climate change, and that 60% of the world´s population is forecast to live in urban areas by 2030 (^{UN, 2018}) .

In Mexico, since 1990, urbanization has rapidly grown in tropical latitudes (^{Reyes and López, 2011}). These trends of tropical urbanization and rising temperatures are seen in Tuxtla Gutiérrez (TGZ), which is the second most populated city of the Mexican tropics and home to 10% of the population of the state of Chiapas.

In the last three decades, temperatures in TGZ have risen in the dry season (^{De la Mora et al., 2016}). However, there is little information about how much urbanization has contributed to rising temperatures and UHI impacts, people´s wellbeing, economy and public health (^{Coello, 2015}; ^{Díaz-Nigenda et al., 2018}; ^{López-Pérez et al., 2019}). Therefore, we examine the relationship between land-use change and extreme temperatures in TGZ.

UHI results from thermal equilibrium between atmospheric temperature (AT) and land cover with different land surface temperature (LST) (^{Oke, 1982}). UHI intensity depends on the albedo and thermal inertia of construction materials, building density, aspect ratio, sky view factor (SVF), industrial development, traffic congestion, cloud cover, wind speed, precipitation, city geometry, and relative humidity. However, the main diver of UHI is land-use change (^{Li et al., 2018}), a factor that explains about 70% of the total variance in land surface temperature (LST) (^{Oke, 1973}; ^{Imhoff et al., 2010}).

Recently, remote sensing data have helped UHI research by providing low-cost LST data at adequate temporal frequency (^{Sobrino et al., 2004}; ^{Amanollahi et al., 2016}; ^{Avdan and Jovanovska, 2016}). Sensors such as MODIS (Moderate Resolution Imaging Spectroradiometer) and AVHRR (Advanced Very High Resolution Radiometer), offer thermal infrared information suitable for regional studies due to their spatial scale and spectral resolution (^{Mas, 2011}; ^{Amanollahi et al., 2016}). Landsat imagery has been widely used at a finer scale for research (^{Zhao-Liang, Bo-Hui, et al., 2013}; ^{Chen et al., 2017}; ^{Tan et al., 2017}).

LST requires information on radiance and emissivity which can be measured in situ or by remote sensing (^{Sobrino et al., 2004}, ²⁰⁰⁸; ^{Cook et al., 2014}; ^{Jeevalakshmi et al., 2017}; ^{Tan et al., 2017}). Landsat 4, 5, and 7 images have one thermal band and Landsat 8 has two. Before Landsat 8, NASA and the United States Geological Services worked to develop LST products based on the thermal archive data of Landsat missions (^{Cook et al., 2014}). For the purpose, they automated two methods: a) Single channel technique for imagery with one thermal band, and b) Multichannel technique (Split window) for Landsat 8 (^{Weng et al., 2004}; ^{Sobrino et al., 2008}; ^{Tan et al., 2017}). The accuracy of both methods depends on the atmospheric profile information available for the scene location and scene capture time (^{Zhao-Liang, Hua, et al., 2013}; ^{Cook et al., 2014}) especially for the emissivity retrieval.

Surface emissivity is a measure of the energy emitted as thermal radiation from a given land coverage (^{Sobrino et al., 2012}). When measured with a remote sensor, it varies with geometry, surface materials, surface roughness, ground humidity, view angle, and atmospheric information such as relative humidity, wind speed, cloud cover and precipitation (^{Kahle, 1987}; ^{Suomi, 2014}; ^{Ermida et al., 2017}; ^{Gulbe et al., 2017}). In absence of atmospheric information, researchers can use other methods for surface emissivity retrieval. The most common method is the Normalized Difference Vegetation Index threshold method, NDVI^THM (^{Sobrino et al., 2008}), which assigns an emissivity value to a given land coverage, bare soil or full vegetation, according to estimated NDVI values. However, NDVI^THM has limitations; surfaces such as asphalt, or other common materials in urban environments, can cause inaccurate estimations of LST (^{Sobrino et al., 2008}). A more accurate method than NDVI^THM for urban areas is the Emissivity Classification Method (ECM), which assigns emissivity values to each land cover-land use class, based on the ASTER spectral library (ASL, https://speclib.jpl.nasa.gov/) reported by ^{Sobrino et al., 2012}.

In this research we examined the relationship between LULC and extreme temperatures in TGZ. Fist we estimated the surface UHI for the city of TGZ for three years (2001, 2011, and 2017) from historical analysis of air temperature. Later, we used Landsat 5 (TM) and 8 (OLI) imagery to define land use/land cover. Based on the cover types we retrieved LST by ECM method, as emissivity data were not available for the study area.

Characterization of the area

The city of Tuxtla Gutiérrez (TGZ) is the capital of Chiapas, a state located in southeastern Mexico between 16º 48' 46" N, 93° 13' 12" O and 16° 40' 49" N, 93 ° 02' 19" O coordinates (Figure 1). Elevation ranges from 435 m to 917 m (Figure 1). The climate is warm and sub-humid with an average annual rainfall of 915 mm the majority of which falls in the summer (^{INEGI, 2017}). The historical records registered 19 days with temperatures above 30 °C and 25.7 days of dry conditions per year with a mean annual temperature of 25.4ºC. The hottest month is April with a mean maximum historical temperature of 45.5ºC. The city of TGZ occupies a valley surrounded on the south by Mactumatzá Hill and on the east by The Sumidero Canyon National Park; both barriers hinder air flow and influence winds from the northwest with an average speed of 9.7 km/hr, reaching gusts during winter of 36 km/hr (^{Díaz-Nigenda et al., 2018}).

Figure 1 Location of weather stations in the city of Tuxtla Gutiérrez, Chiapas, Mexico.  

The Sabinal River crosses the city through The Sumidero Canyon National Park to the Grijalva River. The predominant types of soil are Leptosols, Regosols, Luvisol and Vertisols. The city is surrounded by 4.9 km of Natural Protected and Conservation Areas of low deciduous forest (^{ICIPLAM, 2013}). Almost 75% of the urban flora is exotic (^{Román, 2017}) and distributed among 193 parks and a country club (^{ICIPLAM, 2015}).

The city, as the main economic center of the state, had a population of 537,102, in 2010, which is an increase of 10% from 2005 (^{INEGI, 2013}; ^{CEIEG, 2014}) and between 2007 and 2018 it grew in area from 7.8 km² to 107 km² (^{HATG, 2007}), an increase of 11% annually. Since 1986, residential development has an annual increase of 3% in TGZ (^{Silva et al., 2015}). The Municipal Citizen Planning Institute (ICIPLAM) of TGZ has identified an irregular land use-land cover (LULC) change in the city since the 1980s, one that limits urban planning and interfere with the implementation of green areas as part of urban developments (^{ICIPLAM, 2013}) since vegetation competes with housing space.

Materials and methods

Methodology

We were interested in the impact of extreme temperatures on human wellbeing and comfort. Therefore, we used vulnerability levels due to extreme air temperatures defined by The Territory and Urban Development Ministry (^{SEDATU, 2014}) as reference values for the analysis of results (Table III).

Table III Vegetation and human vulnerability due to extreme air temperatures. 

Temperatures	Class	Vulnerability
< 28 °C	Pleasant	Comfortable wellness
28 a 31 °C	Discomfort	Evapotranspiration increases. Headaches increase in humans.
31.1 a 33 °C	Extreme discomfort	Dehydration is evident. Hoppers and heavy pollution particle pollution increases, appearing at cities.
33.1 a 35 °C	Stress condition	Plants evapotranspire excessively and wilt. Forest fire hazard increases.
> 35 °C	Upper tolerance limit	Heat strokes occur, with unconsciousness in some people. The diseases increase.

Most LST research has been done during the dry season because cloud cover is low and higher atmospheric humidity can increase the error of emissivity and temperature. For this reason, we confined our research to the dry season in TGZ.

To calculate LST, we applied conventional and thermal remote sensing corrections to the three scenes downloaded (^{Chander et al., 2009}; ^{USGS, 2013}). This correction was done in two stages. First, we considered the build-up of an emissivity layer according to the ECM, then, we retrieved LST by transforming thermal brightness (Figure 4). The combination of band thermal and radiometric correction (Figure 4), allowed to obtain LST values for each image.

Figure 2 Diagram for land surface temperature (LST) retrieve based on Landsat imagery. Dash arrows indicate the band thermal correction process up to the thermal brightness while the radiometric correction and supervised classification to assign emissivity is indicated by the black arrows.  

Figure 3 Thermal hourly distribution registered at day and night in the city of Tuxtla Gutiérrez from years 2013 to 2018. 

Figure 4 Average daily temperature and precipitation for the city of Tuxtla Gutiérrez, from 1951 to 2016 year.  

Though atmospheric conditions are the main source of error for LST retrieval (^{Cook et al., 2014}), we made no correction of this sort due to a lack of information on relative humidity, wind speed, and atmospheric pressure.

Land use-land cover image classification

We radiometrically corrected each band of the images up to reflectance (Table IV) using equations 1 and 2 (^{Chander et al., 2009}). The correction served to sort the images, into nine LULC classes and thereby assign an emissivity value to each class.

Lλ=Lminλ+Lmaxλ-LminλQcalmax-Qcalmin*(Qcal-Qcalmin) (1)

ρλ=π∙Lλ∙d2ESUNλ∙cos⁡θs (2)

Table IV Radiometric correction parameters for TM 5 and OLI 8 sensors.  

TM 5 Sensor, 2001 year
Band	QcalMAX	QcalMIN	LMAX	LMIN	d	ESUN	SZ
2	255	1	365	-2.84	1.0005897	1796	31.31
3	255	1	264	-1.17		1536
4	255	1	221	-1.51		1031
2011 year
2	255	1	365	-2.84	0.999269	1796	30.13
3	255	1	264	-1.17		1536
4	255	1	221	-1.51		1031
OLI 8 Sensor, 2017 year
	ML	AL	Mp	Ap	d		SZ
3	0.011861	-59.3059	0.00002	-0.1	0.9994378		27.84
4	0.010002	-50.0101
5	0.0061207	-30.6037

* Qcalmin = Minimum quantized value of pixel [dimensionless], Qcalmax = Maximum quantized value of pixel [dimensionless], LMIN = Spectral at-sensor radiance that is scaled to Qcalmin [W/(m2 sr μm)], LMAX = Spectral at-sensor radiance that is scaled to Qcalmax [ W/(sr m2 sr μm)], ESUN= Mean exoatmospheric solar irradiance [W/(m2 μ m)], Angle SZ= Solar zenith angle [degrees] (^{Chander et al., 2009}; ^{USGS, 2013}).

Based on a high-resolution images and the orthophoto mosaic, we generated for each image classification training polygons consisting of nine classes of LULC (Table IV). These polygons were used to train a maximum likelihood algorithm (MAXLIKELIHOOD) and classify each image (^{Melesse et al., 2007}; ^{Sobrino et al., 2012}) .

The bare soil class denotes areas where the surface consists of bedrock (mainly limestone). The bare agriculture soil class denotes agriculture soil without vegetative cover. The rest of the classes are self-explanatory (Table V).

Table V Emissivity values by type of sensor (TM 5 and OLI 8) and LULC classes. 

Sensor	TM 5	OLI 8
Band	B6	B10
Range (µm)	10.4 a 12.5	10.6 a 11.19
Water*	0.99	0.99
Trees+	0.99	0.99
Asphalt road+	0.963	0.962
Green grass*	0.981	0.98
Bare soil+	0.956	0.942
Bare agriculture soil*	0.968	0.941
Concrete ceiling+	0.957	0.943
Metal ceiling*	0.046	0.049
Roof tile ceiling+	0.942	0.939

* Source: Aster Spectral Library. ⁺Source: (^{Sobrino et al., 2012})

Emissivity Classification Method (ECM)

Emissivity Classification Method (ECM), which assigns emissivity values based on previous measurements or from libraries such as the ASTER Spectral Library (ASL, https://speclib.jpl.nasa.gov/) to generate land use and land cover (LULC) classification (^{Sobrino et al., 2012})

We applied ECM for emissivity retrieval (Table V) to avoid homogenizing the entire urban area into two or three values of emissivity, which would have happened if we relied on NDVI^THM.

Land surface temperature (LST) retrieval

Once an emissivity raster coverage was estimated, LST values were obtained, following equation 1. We applied radiometric correction of thermal bands up to radiance. Then, using equation 3 (^{Chander et al., 2009}) and conversion from Kelvin to Celsius degrees with equation 4, we estimated thermal brightness values from radiance.

TB=K2lnK1Lλ+1  (3)

TBC=TB-273.15 (4)

Finally, we calculated per-pixel LST using equation (5):

LST=TBc1+λTBcρlnE (5)

TB is the thermal brightness [K]; TB_c is the thermal brightness in Celsius degrees [C]; LST is the land surface temperature in Celsius degrees [C]; K₁ is the first calibration constant [W/ (m² sr μm)]; K₂ second calibration constant [ K]; L_λ is the spectral radiance at-sensor [W/(m² sr μm)]; ln is natural logarithm; E is the emissivity; λ is the center of the thermal band [µm]; ρ is a constant that combines the Boltzman and Plank´s Law and the speed of light with a value of 14,387.7 [µmºC] (Table VI).

Table VI Thermal correction parameters to estimate the thermal brightness values. 

Sensor TM 5
Band	K1	K2
6	607.76	1260.56
Sensor TIRS 8
10	774.89	1321.08

Results and discussion

Daily air temperature data, from 1951 to 2016, showed that after 1980 temperatures above 40 °C have become more frequent by a rate of 18% and minimum air daily temperatures above 28 °C have increased by a rate of 14%.

Hourly records from 2013 to 2018 show that maximum daily temperatures occur between 14:10 and 17:00 hrs. Daytime temperatures, between 06:00 and 19:59 hrs, range from 8.6 °C to 48.4 ºC. Nighttime temperatures, between 20:00 and 05:59 hrs, range from 9.2 °C to 32.5 °C. During the time scale study, 60% of daytime temperatures and 63% of nighttime temperatures occurred between 25 °C to 30 °C, suggesting that nights are as warmer as days in TGZ (Figure 3), but it is necessary to continue registering hourly air temperatures to confirm this trend.

As there are no weather stations located on the periphery of the city, we were unable to compare temperatures in the city center and the surrounding rural areas. In addition, there are no night Landsat images that record land surface temperature that would allow us to study diurnal variation in UHI.

On average, the hottest months, are April and May, with average daily temperature of 27.2 and 27.8 °C, respectively. During these months, the maximum temperatures occur with clear skies and low relative humidity, ideal conditions for classification of LULC and LST retrieval (Figure 4).

This information was used to determine the Landsat scenes (Table II) in order to reach an optimum accuracy of the classification processes (Table VII).

The classification processed showed that tree coverage decreased by 21% in TGZ (Table VIII). This loss of tree cover was concentrated along the Sabinal River, Botanical Garden, Zoo, some urban parks, and natural reserves areas on the south of the city (Figures 6 to 8).

Figure 5 Total area for each land surface temperature (LST) range analyzed in the years of study. 

Figure 6 Land surface temperature and land cover during April 5th, 2001 in the city of Tuxtla Gutiérrez, Chiapas. 

Figure 7. Land surface temperature and land cover during April 1st, 2011 in the city of Tuxtla Gutiérrez, Chiapas. 

Figure 8 Land surface temperature and land cover during April 1st, 2017 in the city of Tuxtla Gutiérrez, Chiapas. 

The LST distribution for the three years of the study are given in Table IX and Figure 5.

There are street trees in the area but quantifying their impact on temperature and surface variation requires finer-resolution imagery, species identification and tree health recognition (^{Román, 2017}). In the southwestern part of the city, where substantial residential development took place during the early 2000s (^{Silva et al., 2015}), vegetation coverage has decreased and BAS, concrete, and roof tile ceiling has increased. This change in land cover was accompanied by temperature changes in TGZ. Specifically, temperatures shifted from discomfort temperatures (29 °C to 31 °C) in 2001 to extreme discomfort (31 °C to 33 °C) and stress condition temperatures (33 °C to 35 °C) in 2011 and 2017, respectively.

The urban tree land-cover class had the lowest LST. In 2001, around 25% of urban tree areas were located at pleasant to discomfort temperatures zones, reaching only 11% by 2017. Large continuous areas of dense tree cover had the lowest LST, for example, parks such as Caña Hueca, Joyyo Mayu, Tuchtlán, Botanical Garden, Parque Oriente and Parque Patricia, with areas from 3 to 20 ha, registered temperatures below 28 ºC LST. By contrast, parks smaller than 3 ha or with very dispersed and low tree density, such as urban parks or wasteland (1.3 to 13.3 ha), had higher temperatures (Figure 9). It is evident at the spatial scale of the present study, that the Sabinal River has a cooling effect (Figure 6 to 8), which crosses the city in the west-east direction, and its riparian vegetation. For the three years analyzed, we observed that along this main-channel, and in its open sky sections the maximum temperatures ranged from 28 to 35 ºC. We observed, the same mitigating effect in other tributary channels to the east of the city.

Figure 9 Comparison of land surface temperature (LST) registered at two contiguous woodland areas with different tree density.  

The BAS class is made up of dark soil clays (mainly Vertisols, Regosols, and Leptosols) devoid of vegetation that is distributed at the periphery of the city. Its extent decreased 67% (3,650 to 1,214 ha) between 2001 and 2017, changing into roof tile ceiling, asphalt road, and concrete ceiling classes. For the three years analyzed, BAS mean temperature was 35 °C, but when it changes into asphalt road its mean value decreased 1 °C for and 3 °C for concrete ceiling. Thermal characteristics of soil and the day/night behavior of UHI highly influence the BAS class.

Soil physicochemical characteristics determine its capacity to absorb and retain heat. Dark and damp soils with high organic matter content tend to retain a greater amount of caloric energy and lose it more slowly. By contrast, light, stony and dry soils have a heterogeneous behavior; they can be heated considerably but also lose energy quickly or absorb it at lower rates and therefore show less caloric emission (^{Owens and Rutledge, 2005}). Climate also influences soil temperature; in tropical regions soils tend to have higher temperatures (from 2 to 4 °C above the AT), due to microbial activity, compared to soils of the same type in temperate regions (^{Buol, 2005}).

Vertisols are dark soils, with a high content of clays considered isohyperthermic, thermal and mesic (^{Virmani et al., 1982}; ^{Palma-lópez et al., 2002}; ^{FAO, 2007}). These are the type of BAS that presents the maximum retrieve LST for the three dates of the present study. They are located on the west of the city, as temporary agricultural land at rest; that is, without crops at the time the images were acquired. During field verification we observed moisture and stubble. This land use faces conversion pressure, as 11 residential developments were registered nearby.

Leptosols and Regosols are shallow, stony, grouped as Entisols, and are considered hyperthermic and thermal soils (^{Lorenz, 1995}; ^{Palma-lópez et al., 2002}; ^{FAO, 2007}). Leptosols are located on the north and east of the city on steep slopes. Between 2000 and 2015 these areas showed the greatest urban sprawl and land conversion. Now, however, they are areas of conservation. On the other hand, Regosols are located on the south of the city. Their thermal behavior is the most heterogeneous. They register temperatures lower than 28 °C to 39 °C. Possibly due to the stony conditions, Regosols´ temperature rises considerably during the day, but drops quickly in cloudy conditions, or at night. However, this behavior has not been verified since there are no Landsat scenes of the study area at night.

The physical characteristics of each land coverage are also related to the day/night behavior of UHI, a phenomenon that has been extensively discussed (^{Oke et al., 1991}; ^{Jaureguı, 1997}; ^{Prakash et al., 1999}; ^{Azevedo et al., 2016}; ^{Ando and Ueyama, 2017}; ^{Kato et al., 2018}). Bohnenstengel et al. (2011), simulated the behavior of the UHI in the city of London, evaluating the energy flow hourly. The study found that the temperatures of rural and urban areas depend on thermal inertia. Usually, thermal inertia of rural areas is low. Increasing LST in rural areas, requires less energy than doing so urban areas. Nevertheless, diurnal behavior of both areas, rural and urban, depends not only on thermal inertia but also on cloudiness, relative humidity, direction and speed of wind, and constructive materials. Similar observations in Arizona and Mexico City where peri urban areas recorded higher temperatures than cities´ downtowns during the day and the inverse behavior at nights (^{Jaureguı, 1997}; ^{Tejeda-Martínez and Jáuregui-Ostos, 2005}; ^{Imhoff et al., 2010}; ^{Ruddell and Dixon, 2014}).

Areas with both, metal ceiling and asphalt roads, such as malls like Plaza de las Americas, Plaza Cristal and Plaza Polyforum showed thermal peaks (temperature above 40 °C). These results described are consistent with those reported by ^{Coello (2015)} who used a thermographic camera to retrieve temperatures for different coverages at the main campus of The University of Science and Art of Chiapas between 13:00 and 15:00 hrs. Coello reported that maximum temperatures occurred at the metal ceiling of sport courts (60.5 °C), followed by football camp with synthetic grass (55 to 59 °C), bare soil with the large range of temperature variation (46 to 58 °C), asphalt roads (52 to 57 °C), concrete (47 to 54 °C), scattered vegetation (33 to 38 °C) and swimming pool (28 to 29 °C) (^{Coello, 2015}).

One transect profile was generated across TGZ form east to west direction to validate the changes in LST due to of land cover changes occurring from 2001 to 2017 (Figure 10). LST has increased due to transformation of trees into SDA or non-evaporating surfaces and decreased when land converted to trees (Table X).

Figure 10 Transect profile across the city of Tuxtla Gutiérrez, Chiapas used for validating the relation of land surface temperature (LST) with land use and land change cover. 

The mean temperatures reported at Table X were those occurring at transect profile (Figure 10). The highest increase in LST occurred when trees became BAS (3 °C), followed by asphalt road (1.7 °C) and bare soil/concrete ceiling (1.4-1.3 °C). On the other hand, TGZ city had reforestation programs at median strips, parks and vacant lots in 2005 and 2016, and these programs appear to have a modest cooling effect flourishing (0.1 and 1.2 °C) by 2017 year.

On the 29th and 30th of May 2019 (Table XI) between 15:00:00 and 17:17:00 hrs, temperature verification performed in 10 sites across the city of TGZ. The sites were chosen considering the areas that reported high LST during 2017. May 29th presented sunny conditions, while May 30th was cloudy with drizzle. On both dates, the maximum temperatures were recorded on site 03, in the east of the city, which had BAS Leptosol with pebbly black earth. The second highest temperatures appeared on site 02, in the west, with BAS Vertisol and black soil with residual stubble. Synthetic pit grass was the second warmest coverage, followed by asphalt. Direct temperature measurement showed the same relative patterns, although there were some differences between directly measured and remotely sensed temperatures.

Notwithstanding, this study was limited by the resolution and only day time availability of the satellite imagery, the previous studies (^{Coello, 2015}) and verification performed supported the trend of the different land uses and coverages in relation with LST. The direction of land use change determines the raised or decreased of LST, but in general, when vegetation is transformed into any other land use the LST increased at least 1 °C and decreased the same in the opposite direction. The minimum variation of LST detected in land use change from asphalt to tree (-0.1 °C) may be caused by the imagery resolution (30 m per pixel) and height of trees, since reforestation programs are recent, and according to ^{Coello (2015)}, the refreshing effect of vegetation may be more than 10 °C in comparison with asphalt and BAS cover.

The cooling effect of vegetation in tropical cities may contribute to the wellbeing of the population by mitigating UHI, but more studies are needed in economic and social terms in each city for a strategy to determine the sites and area of reforestation, the maintenance work required, the best species for reforestation and the availability of soil and water needed.

Conclusions

For the city of TGZ, the materials with the highest thermal sensitivity were BAS and asphalt roads. The change of tree cover into asphalt road, bare soil or BAS cover caused LST increased 1 °C to 3 °C. When asphalt or BAS cover classes turned into the tree class, LST decreased from 0.1 °C to 1 °C, suggesting that the cooling effect of trees is not immediately reached after reforestation. The vegetation effect in LST is appraised if woodlands are greater than 3 ha indicating a limitation of the analysis caused by the spatial resolution of the imagery used.

The ECM for the thermal assessment captured the surface heterogenicity of the TGZ city avoiding underestimations induced for NDVI methodologies. Nevertheless, the in-situ measures of emissivity values and night temperature assessment will improve results and are desirable for further investigations.

The approach implemented in this study to carry out the analysis could be used to evaluate the surface thermal dynamic of cities when monitoring weather networks are not available or are not wide enough to cover the study area properly.