Land Surface Temperature Relationship with the Land Use/Land Cover Indices Leading to Thermal Field Variation in the Turkish Republic of Northern Cyprus

The increase in the Land Surface Temperature (LST) caused by global warming and extreme weather events is significantly increasing Urban Hot Spots (UHS) while impacting the environmental quality within urban areas. In this research, with the help of Landsat 5, 7, and 8 satellite images, the evolution of Land Use/Cover (LULC), LST, and Surface Urban Heat Island (SUHI) in the Turkish Republic of Northern Cyprus (TRNC) between the years 1985 and 2020 is examined. The main aim is to evaluate how these variables, together with the Normalized Difference Built-up Index (NDBI), Urban Index (UI), Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI), Built Index, and Proportion Vegetation (PV) have influenced the variability of the UHS and the level of thermal comfort according to the Thermal Field Variance Index (UTFVI). The results report significant average rises of the variables: LST (6.62 °C), SUHI (0.97 °C), UHS (8.61%) of the UI and NDBI indices and class 6 of UTFVI (88.53%) that are related to the variability of the LULC corroborated by statistical analysis. Our results provide valuable information on the future development of urban areas in Northern Cyprus to make them more resilient and sustainable to rising temperatures.


Introduction
Humanity is currently facing one of its most critical and urgent challenges: global warming and its effect on extreme weather events (Kovats et al. 2005;Song et al. 2020). Rapid population growth can be seen as an essential catalyst for the expansion of urban areas into rural areas, while urbanization is understood as the main contributor to climate change ). According to a report by the United Nations Organization (UNO), currently 50% of the present population lives in urban areas, and this is expected to increase to 70% by 2050 (UNO 2018). This critical change means that in the next 30 years, the urban population will increase by 2600 million people (Mukherjee and Singh 2020), which will require an urbanization process of approximately 1,527,000 km 2 (Schneider et al. 2010).
On the one hand, urbanization can be considered a fundamental means of expanding industry and transport, generating urban and economic growth. On the other hand, it modifies the local urban climate by increasing the LST (Scolozzi and Geneletti 2012;Song et al. 2020). Numerous studies have explored positive correlations between high temperatures and urban areas. The development of green areas can be consideredone of the main measures that allow mitigating increases in LST caused by urban areas (Founda and Santamouris 2017;Hua et al. 2020;Yang et al. 2020). The increase in temperatures is joined by the phenomenon of urban climate change, Urban Heat Island (UHI), which intensifies due to pollution generated by numerous human activities (Santamouris 2020). Thus, the study of Yao et al. (2021) on mainland Chinese cities reveals that temperature 1 3 Published in partnership with CECCR at King Abdulaziz University increases result from human activities as well as impermeable urban areas.
The problems and impacts (climatic, environmental and socioeconomic) generated by the UHI on urban areas are widely discussed in scientific research, reducing the inhabitants' quality of life (Rozos et al. 2013). Among them, the decrease in air and water quality (Feizizadeh and Blaschke 2013), reduction in biodiversity (Čeplová et al. 2017), increase in the cost of energy (Santamouris 2020), increase in mortality (Arbuthnott and Hajat 2017) and even changes in the regional climate (Sarrat et al. 2006).
In this way, there is a process of decreasing people's physical and mental well-being due to the changes that the increases in coverage produce on the climate (Das and Das 2020). In this sense, it is essential to evaluate citizens' quality of life and the consequences that high temperatures can generate in their lifestyles. Within the scientific community, several indices allow for assessing the thermal well-being of urban areas: Physiological Equivalent Temperature (PET), the Urban Heat Island Intensity Index (UHIII) or the UTFVI. The latter is the most frequently used one since it allows for identifying high-temperature thermal spaces defined as UHS (Amindin et al. 2021). It relates them to the different typologies coverage by identifying the different LULC. Several studies have warned of the increase in UHS values due to the increase in LST (Sharma et al. 2021). For example, the study between the years 2011 and 2019 on the city of New Okhla (India) reported increases in UHS and LST of 33.56% and 6.4 °C, respectively (Sharma et al. 2021), the study between the years 1995 and 20,018 over the city of Ahvaz (Iran) reported increases in UHS and LST of 4% and 3.2 °C. Research on the metropolitan area of the city of Granada (Spain) between 1985 and 2020 described increases of 2.2 °C in the LST and 20.4% in the UHS (Hidalgo and Arco 2022), while a study on the Suez area (Egypt) between 1984 and 2014 revealed an increase of 16% in UHS and 4.5 °C in the LST (Ahmed 2018). Comprehensive research (Yao et al. 2019) carried out between 2001 and 2017 covering approximately 400 cities in China reported equal daytime and nighttime increases in SUHI. However, these studies are consistent in pointing out the areas of bare soil, without vegetation, or industrial zones as the places where the most significant increases in LST occur and, therefore, where the UHS are located. In addition, these have worse environmental quality indices.
In recent years, remote sensing has emerged as one of the most used methodologies for determining these phenomena, allowing studies of large urban areas of LST, UHS and LULC to be carried out (Song et al. 2018). The study of the relationship between the LST and land use and land cover is essential to understand the effects produced by variations in cover on temperature increases (Hidalgo and Arco 2022;Tepanosyan et al. 2021). Therefore, studying and identifying the relationships mentioned above to mitigate these changes in the urban environment becomes critical for maintaining the climatic balance within a city. In addition, these analyses can provide important information regarding the future development of urban areas that makes them resilient and sustainable to rising temperatures.
The reasons that justify this research are: 1. The study country does not have LST, UHS and UTFVI studies that consider temporal variability as proposed in this research. 2. The study period contemplated in this research (35 years) is considered an appropriate question since it will allow studying the variability of the indices over a long period of time. This goes against current studies (Barbieri et al. 2018;Mukherjee and Singh 2020) that use specific values to extrapolate to more extended periods and obtain global results.
Throughout this research, we intend to answer the following questions: 1. How have the variables taken into account in the area under study evolved between 1985 and 2020? 2. What are the relationships between these variables taking into account the insular situation of the country? 3. How do LULC changes influence the evolution of UHS? 4. What evolution has the area under study had in relation to environmental quality according to the UTFVI system? 5. Bearing in mind that the country has numerous development perspectives, can the results of this research be considered in future urban developments in order to minimize LST increases?

Study Area
This study was conducted in the northern part of Cyprus, an Island in the eastern basin of the Mediterranean Sea ( Fig. 1), situated between latitudes 30.33° and 35.11° N, and longitudes 32.23° and 34.55° E. The landscape is shaped by coastal areas, plain land on the inner parts of the islands, and mountainous areas. According to the TRNC Statistical Institute, with an area of approx. 3242 square kilometres, the population of North Cyprus was 382.836 in 2020 projected from the last census data obtained in 2011 (TRNC 2021). North Cyprus comprises of six districts or cities: Lefkoşa (Nicosia), Gazimağusa (Famagusta), Girne (Kyrenia), Güzelyurt (Morphou) and Lefke (TRNCPIO 2022). Except for the capital city Lefkoşa, which is an inland city, all urban centers are coastal cities.
The climate in Cyprus is considered typical Mediterranean with hot to dry summers and warm to rainy winters (TRNCPIO 2022). According to the data obtained from the TRNC Statistical Institute for 2020, the hottest months are from June to August and the coldest months from December to February. Rainfall is generally perceived between October and March, while the highest rainfall was observed in the mountainous areas in March. Temperature changes can vary related to the location as coastal, inland, or mountainous area. For example, while the average winter temperature is between 12 and 15 °C, in mountainous areas, it cooled down to approx. 5.6 °C in January 2020. Similarly, the average temperature in summer is 32 °C, while the temperature in inland areas, such as Nicosia, can reach up to 43 °C (TRNC 2021).
In the last decades, North Cyprus has faced a vast urbanization process, which resulted in a decrease of natural and farmland and increased building development in rural areas (Atak et al. 2019;Fuladlu et al. 2021). Especially, larger urban centers such as Nicosia as well as the coastal cities such as Famagusta growing outside their boundaries, increasing suburban development and sprawl (Constantinides 2001;Fuladlu et al. 2021). Figure 2 shows graphically the methodology followed in this research. The objective pursued by this research is to analyze the development of the LST and LULC variables between the years 1985 and 2020 in the Turkish Republic of Northern Cyprus (Cyprus). In turn, it is intended to analyze how these variables together with the rest of the analyzed indices: UI, NDVI, NDBI, NDWI, PV and Built up index; have allowed the development of new UHS and the deterioration of UTFVI air quality. To do this, Landsat satellite images (5, 7 and 8) have been used to obtain the evolution of the indicated indices and the LULC map using the Support Vector Machine (SVM) system. Next, using Landsat thermal images, the LST has been obtained and the evolution of the UHS has been analyzed together with the level of thermal comfort with the UTFVI methodology. With the help of statistical analysis and the Data Panel technique, the correlation between the data obtained and the relationship between the variables has been determined. This statistical analysis methodology is presented as a novel element with respect to the traditional methods of multiple correlations usually used in research, since it allows the incorporation of a larger number of variables and data by admitting the inclusion of the individual effects of a particular area and by date to get comprehensive results. Therefore, using this methodology, it is possible to obtain more precise and complex results than traditional methods.

Methodology
First, the NDVI, PV, NDWI, NDBI, Built up index and UI indices have been determined. Next, the LST has been obtained using Landsat 5 images for the years 1985, 1990 and 1995, Landsat 7 for the years 2000, 2005 and 2010 and Landsat 8 for the years 2015 and 2020 using the software QGIS. Following the SVM method, LULC maps were obtained for the indicated years. Using a precision matrix, the precision of the land cover was determined through the cross-tabulation system between the classified category and the reference category (Campbell 1996;Xu et al. 2009;Yoo et al. 2019). Next, the UHS zones were identified and the UTFVI index of the area under study was calculated. The analysis of the data obtained in this research has been carried out statistically with the help of the software STATA.

Landsat
The images for the years 1985 and 1990 were obtained from Landsat 5, the images for the years 2000, 2005 and 2010 were obtained from Landsat 7. Finally, the images for 2015 and 2020 were obtained from Landsat 8. All were downloaded from the Service United States Geological Survey (USGS) since they are subjected to an atmospheric correction process prior to discharge. Optical and thermal bands from the Landsat dataset were used to determine indices and LST. Landsat 5 has six optical bands (1-5 and 7) with a pixel size of 30 m of resolution and a thermal band (6) of 120 m of pixel size. Landsat 7 maintains the bands and sizes of its predecessor in the optical bands but decreases the pixel size of the thermal band to 60 m. Finally, Landsat 8 increases the number of optical bands to 8 (bands 1-7 and 9) and the number of thermal bands to 2 (bands 10 and 11). The first ones keep the same resolution, but the thermal ones increase up to a pixel size of 100 m. For our research and to determine the LST, only band 10 Landsat 8 was used. All the bands were applied through a pan sharpening process in order to reduce the pixel size of the bands. The area under study is under the path of the Landsat 5, 7 and 8 satellites. Although the weather in Cyprus is usually sunny throughout the year, the images were acquired for the month of June each year. In this way, the land cover classification process is improved through the methodology used, since cloud cover is close to the 0% value. 8 images have been used, 1 per year.

Radiance Spectral ( L λ )
TIRS band digital numbers (DN) have been converted to spectral radiance using the following Eq. (1) (Kafer et al. 2019): where L is the spectral radiance expressed in W/(m 2 •sr•μm); M L is 3,342 × 10 -4 W/(m 2 •sr•μm) indicated in the Landsat metadata file; Q Cal is the digital value (DN) of the bands that varies in a range from 0 to 255 and A L is a dimensionless rescaling value located in the image file.

Brightness Temperature
Next, and in order to obtain the brightness temperature (T) in º Centigrade, the spectral radiance has been converted. This is the actual temperature recorded by the satellite with the precondition that the emissivity value is one, exptressed as follows in , Eq. (2) (Weng et al. 2004): Landsat 8: K 1 = 774.8853 and K 2 = 1321.0789; Landsat 7: K 1 = 666.09 and K 2 = 1282.71 and Landsat 5: K 1 = 607.76 and K 2 = 1260.56 and L is the spectral radiance derived from Eq. (1).

Emissivity (ε)
In the present study, the ε correction based on the NDVI determination method was performed using Eq. (3) (Sharma et al. 2021) as follows: where the Pv is obtained from NDVI (Eq. 7).

LST
LST was estimated using ε values derived from Eq. (3) and using Eq. (4) (Weng et al. 2004). Equation (5) is used to determine the coefficient C 2 as follows : where LST is the Land Surface Temperature, T is the Brightness temperature of the Landsat, is the wavelength of the emitted radiation (Landsat 8: =10.8 µm and Landsat 5 and 7: = 11.457 µm), is the emissivity of the ground, h is Planck's constant with value 6.626 × 10 -34 J s, c is the speed of light with value 2.998 × 108 m/s and s is Boltzmann's constant with value 1.38 × 10 -23 J/K (Weng et al. 2004). (2)

LULC
The following were identified: vegetation, urban land, bodies of water, bare soil and agriculture. Using Sentinel 2 images from the year 2020 and a precision matrix, the precision was determined by means of a crossed system between both categories: the classified one and the reference one.

NDVI THM
It was calculated using the following Eq. (6): To calculate the PV according to Eq. (7) (Yu et al. 2014). It is necessary to take the results of Eq. (6) as follows: where NDVI is obtained from Eq. (6) and the NDVI max are the maximum values of the NDVI and NDVI min are the minimum values of the NDVI.

Normalized Difference Water Index (NDWI)
The NDWI allows us to determine the hydric stress existing in the vegetation and the humidity of the soil. It was calculated using band green and NIR according to Eq. (9) (Zha et al. 2003):

Urban Hotspots (UHS)
Urban Hotspots (UHS) are mapped using the Eq. (12) (Guha 2017;Sharma et al. 2021): where µ and are the mean values and the standard deviation of the LST of the zone in °C, respectively.

UTFVI
UTFVI is calculated using the following Eq. (13)  where LST is the temperature obtained from Landsat at each point on the surface (°C) and T mean is the average temperature determined with Landsat for the entire area (°C) (Sharma et al. 2021).

Strategy of Analysis
The Panel Data are a model that has been used in studies of this type and requires the need to perform multiple (11) Built up index = NDBI − NDVI, regressions (Alcock et al. 2015;Chen et al. 2011;Fang and Tian 2020). To choose which method is appropriate, follow these steps: (Chen et al. 2011): First, it is necessary to determine whether the effects are random or fixed using the Hausman test (Seto and Kaufmann 2003). 2. Using the Wooldridge and Wald tests, evaluate the model (Labra 2014). For the analysis, the STATA data software, version 16, was used Eq. (14) as follows: where β is the coefficient, X it is the explanatory variable, it is the error residual, and i reported the individual effects.

Evolution of the Indices (NDVI, PV, NDWI, NDBI, Built Up Index and UI)
The space-time analysis of the NDVI, PV, NDWI, NDBI, UI, and Built up index spectral indices between the years 1985 and 2020 can be found in Figs. 3, 4, 5, 6, 7, and 8, while the trend measurement's central and dispersion is represented in Tables 1 and 2. The NDVI index measures the amount and vigor of vegetation on the earth's surface. The highest mean value reported was 0.332 in 2020, while the lowest mean value was obtained in 2005 (0.191). The PV index presents the highest average value of 0.377 in the year 2020 while the lowest average value is 0.158 in the year 2010. The average increase during the years under study has been 176.8% for the NDVI and 58.09% for PV. It is observed how the lowest   in urban areas during the considered years of this study .
In this way, those places with higher values in the NDVI and PV indices present lower values of the NDBI, UI and Built up index and vice versa. It is observed how the identification and determination of these values can be considered an adequate technique to evaluate the evolution of the country's urban and rural environment.

Evaluation of the LULC
The evolution of the LULC between the years under study  can be seen in Fig. 9, while the variation of the different land uses is represented in Table 3.
The coverage identified as bodies of water presents the highest value in 2010 (24,416 ha), while the lowest value is in the year 2000 (1009 ha). The vegetation presents the      On the contrary, there are increases of 245.01%, 716.51% and 11.38% in soils identified as bodies of water, buildings and bare soil, respectively. It is evident that the increase in built-up coverage is due to the rapid urban development that the Turkish Republic of Northern Cyprus has faced in recent decades, motivated mainly by tourism and the construction of second homes. However, the increase in bare ground cover could also be due to speculation leading to an abandonment of farmland for urban use. Moreover, the increase in production costs, the variability in rainfall, and progressive industrialization are circumstances to be considered since these conditions are very common while dealing with the arable lands of the Mediterranean basin. However, according to the data presented in Table 3, these last coverages show constant increases throughout the period studied. On the contrary, the coverage farming, vegetation and water bodies present non-constant decreases with strong oscillations. This could be due to the fact that the development of these areas also depends on the climatic conditions of the area, variables depending on the rainfall of each year. In this way, it would be related that these coverages suffered decreases since the year 2000 the previously determined NDWI humidity index reports that the country is suffering from a moderate-extreme drought. These circumstances have led to the construction of reservoirs and swamps to guarantee the water supply to the population and the need to have a supply system from Turkey.
The most critical changes in the LULC occur in the center of the country, in the areas near the capital (Nicosia), and the coastal areas of the east of the country, the Famagusta region, and nearby areas. As indicated above, the rapid growth of the population and tourism are circumstances that justify these changes in coverage.
The matrix presents some results that can be found in Table 4. This was made with Sentinel 2 images to verify the LULC maps obtained in the year 2020. The accuracy was 84%, the Kappa coefficient was 0.774, and the Tau value was 0.783. However, and after the completion of the matrix, a manual correction was made with the non-coinciding points.

Evolution of LST and SUHI
The variability of the LST and SUHI between the years 1985 and 2020 is shown in Figs. 10 and 11. Table 5 shows the basic statistics of the values obtained. In general, it can be seen that rural areas have lower LST and SUHI values compared to urban areas that have higher temperatures. This implies a lower rate of growth than the average and maximum temperatures. Figure 10 shows the relationship between the PV, NDVI and LST indices. As a result, the LST is lower when the first two indices are higher and vice versa.
The lowest maximum SUHI value was 17.80 °C in the year 2015 while the highest value was 21.59 °C in the year 1995. The lowest minimum SUHI value was −10.27 °C in the year 2020 while the value highest was −8.50 °C in 1985. The lowest mean minimum LST value was 6.11 °C in the year 2010 while the highest value was 10.36 °C in the year 2000. The SUHI spatial statistics of the studied area show a continuous increase from 1985-2020 of 0.97 °C, representing an increase of 12.96%. Figure 11 shows constant increases in the SUHI variable throughout the study area. However, substantial increases are seen in the years: 1990, 2000, 2005 and 2015. It is observed that there is a relationship between the SUHI and the NDVI and PV variables. The higher the value of the first, the lower the values of the last variables and vice versa.  Figure 12 shows the variability of the UHS in the study area between the years under study. It can be seen how between the years studied there has been a significant increase in these. However, this increase acquires a higher growth rate from the readings of the year 2000. The results show ( In 1985, the areas identified as UHS (0.48% of the country's total area) were located in the south of the territory, coinciding with the town of Nicosia. As of the year 2000, these zones' growth began scattered throughout the entire extension. Since 2010 there has been a significant increase in UHS zones, reaching an extension of 6.75%. Figure 13 shows the UHS and the SUHI together between the years 1985 and 2020. A meaningful coincidence between both variables is observed since the areas with the highest SUHI values coincide with the areas classified as UHS and vice versa.

UHS
In Fig. 14, the UHS zones and the 2020 LULC map of a part of the study area have been obtained jointly. The objective of this representation is to analyze the predominance of covers with areas of higher temperatures. According to the results obtained, the UHS occupies an area of 75.66 ha (0.025%) of the coverage called vegetation, an area of 18,225 ha (6.15%) of the built-up coverage, an area of 72,446 ha (24.47%) of the coverage farmland and an area of 205,824 ha (69.35%) of bare soil coverage. Therefore, the types of land where the UHS are located are: bare soil, farmland, and built-up area. On the contrary, the cover that presents the lowest concentration of UHS is the vegetation. The study of the coverages that present    Figure 15 shows the variability experienced in the UTFVI index between the years under study. It can be seen how the red areas, classified as the worst areas for the purposes of the SUHI, have been increasing. Table 7 shows the evolution of the affected area by year under study. The predominance of two large areas is observed: classes 1 (excellent since thermal comfort is optimal) and 6 (worse since thermal comfort is minimal). The areas with the lowest values in the PV and NDVI indices are located in the thermal comfort zones classified as worst, as well as the built-up areas and with high values in the NDBI and UI indices.

Evaluation Using UTFVI
During the entire analyzed period, the best UTFVI classes have presented land decreases of −15.69%, −61.76% and −53.38%, respectively. On the contrary, the worst classes of UTFVI have presented soil increases of 157.87%, 42.37% and 88.53%, respectively. In this way, in the year 2020, 26.55% (1,131,608 ha) of the area under study presents an ecological index UTFVI of the class 6 type. Therefore, an ecological index qualified as Worst with unfavorable conditions for the development of human activities.

Relationship Between LST and Index
The data panel method has been used to establish the relationships between the LST and the NDVI, NDBI, NDWI, PV, UI, Built up index and LULC indices of the study area. In the first place, the Pearson correlation coefficient was determined and then the Data Panel was carried out. Tables 8 and 9 show the results of the analysis.
In general, the LST presents negative correlations with: NDVI (−0. From the statistical analysis, a negative relationship and greater than 99% is reported between the NDVI, NDWI, PV and LST indices and a positive relationship and greater than 99% with the LULC, UI and Built up index indices. Furthermore, a decent relationship can be seen between the data analyzed according to the values of F, R 2 and Prob > Chi 2 since the adjustment level is greater than 99% significance.

SUHI and LST
The results are shown in Tables 10 and 11. It is observed that the SUHI presents an important positive relationship with the LST (0.415). The results of the analysis report a positive relationship greater than 99% between the SUHI and the LST. A good concordance of the data is appreciated based on the variables R 2 , F and the values of Prob > chi 2 .

Relationship Between UHS and LST, SUHI and LULC
Tables 12 and 13 show the results of the data analysis obtained. In general, the UHS presents a negative correlation with the LULC (−0.130) and a positive correlation with the LST (0.277) and SUHI (0.112).
The statistical analysis results using the Data Panel technique report a statistically significant and positive relationship of 95% with the LST and SUHI variables. This becomes negative but not significant with the LULC variable. The F statistic, R 2 and Prob > chi 2 values show a good result between the dependent and independent variables with an adjustment level of 99% significance. These results are in line with the analytical values obtained previously.

Discussion
A general decrease in the mean values of the NDVI and PV indices has been demonstrated in the studied area. However, higher values are reported in rural areas with vegetation as  (Fuladlu et al. 2021), placing our results in the same line. Studies on Iran (Amindin et al. 2021) identified 30% reduction in green spaces between 1995 and 2016 and 30% increase in bare soil cover. In our research, the reduction in green areas was lower (16.32%) mainly due to the protection offered by these areas as a system for minimizing climate change. It is necessary to indicate that the decrease in green areas and the increase in urban areas does not necessarily imply a decrease in the PV and NDVI indices, but also the variability in rainfall. Its reduction as a consequence of climate change considerably affects the state of the vegetation (Li et al. 2002). This is in line with the decrease in the NDWI index values. A significant increase in LST and SUHI values has been detected throughout the years under study. The main reason for this condition is the significant growth of urban areas and its relationship with the decrease of areas for agricultural use and green areas in general. This circumstance was evidenced by the authors Kotharkar and Surawar (2016), who concluded that an increase in green areas produces a decrease in LST and SUHI and vice versa. In our research, higher temperatures have been reported in areas with higher UI, Built up index and NDBI values and lower temperatures in areas with higher PV and NDVI values, evidenced by statistical analysis. The highest correlation of the LST occurred with the NDBI variable, while the lowest correlation was with the PV variable. Other studies have shown this relationship between the LST and NDBI since the latter is not related to changes due to climate change (Anjos et al. 2020;Shahfahad et al. 2021;Tepanosyan et al. 2021). The use of construction systems using waterproof materials with high thermal absorption causes them to heat up with solar radiation and release that heat into the atmosphere at night, producing an increase in LST. As the same effect does not occur in green areas, the differences between rural and green areas are more significant, increasing the SUHI phenomenon (Saaroni et al. 2018;Wu et al. 2019;Yang et al. 2020). This is in line with findings from other studies and allows us to validate our results (Ahmed 2018; Luo and Wu 2021;Sharma et al. 2021). However, our results differ from findings of studies conducted in some regions of Iran (Amindin et al. 2021) and India (Sharma et al. 2021). The reasons could be determined by the difference in population size and surface quality between the studied areas, since our study area is smaller with less population.
A significant increase in the areas classified as UHS has been observed between the evaluated years from 1985 to 2020. This circumstance would be driven by the increase in the LST and SUHI and by the decrease in the vegetation zones. This has been corroborated by the statistical analysis carried out. In such a way that the relationship between LULC and UHS is weak while it reports a high correlation between LST and SUHI. This is motivated by the fact that the reported UHS areas are not developed in all country coverages but mainly in bare soil and farmland coverages. In the area under study, these areas have higher values of LST and SUHI and provide a significant increase in UHS. Numerous studies confirm that solar radiation in rural areas with little vegetation is higher in the morning than in urban areas. Keep in mind that the Landsat 5, 7 and 8 satellites pass through the Turkish Republic of Northern Cyprus in the morning. The leading cause of the higher solar radiation in rural areas is that the shadows generated by buildings and trees prevent the heating of the enclosures and the impermeable surfaces located on the floors of the streets and the cooling rates of the areas with vegetation (Karakuş 2019;Li and Meng 2018;Yang et al. 2020).
On the other hand, some studies carried out with satellite images show how the scarcity of vegetation has a warming effect (Estoque et al. 2017;Lin et al. 2015) in urban areas, while it has a cooling effect when the vegetation is abundant (Lin et al. 2015;Tan et al. 2017;Yu et al. 2017). These effects are not only due to the processes of shading and evapotranspiration, but also depend on the rates of cooling or heating by convection and transpiration. As our results show, this would alter the LST of an area and explain the behavior of the SUHI. Numerous studies (Guha et al. 2018;Shahfahad et al. 2021;Sharma et al. 2021) show these circumstances with results very similar to those reported here.
Lastly, a significant deterioration UTFVI of the Turkish Republic of Northern Cyprus has been evidenced. Hence, the areas with the worst thermal comfort (classes 4, 5 and 6) and the highest temperatures have increased in contrast to the best areas (classes 1, 2 and 3) and the lowest temperatures have been minimized. Therefore, in 2020, the worst category of UTFVI is developed mainly on built-up or bare soil soils, those that usually have higher temperatures. Furthermore,

3
Published in partnership with CECCR at King Abdulaziz University various studies (Guha 2017;Luo and Wu 2021;Majumder et al. 2021;Shahfahad et al. 2021) on other cities or urban areas that corroborate significant increases in worst classes motivated by a rapid increase in the urbanization of the areas and the LST, supporting the results presented here.

Conclusions
In the present study, the evolution of LULC, LST and SUHI between the years 1985 and 2020 in Turkish Republic of Northern Cyprus was analyzed by integrating UI, NDVI, NDBI, NDWI, Built up index and PV, and how it have impacted the increase in UHS and the deterioration of UTFVI. These studies are becoming elements of great importance and in the case of the Turkish Republic of Northern Cyprus, it can represent an important advance. Although the country has experienced significant economic and building development in recent years, it still has high prospects for urban development.
During the analyzed period, the area under study has experienced increases in the mean values of the NDBI, Built up and UI indices. On the contrary, reductions in the mean values of NDVI, NDWI and PV have been experienced. Furthermore, LULC coverage shows significant increases in built-up and bare soil coverage and decreases in vegetation and farmland coverage, indicating that the area under study has experienced rapid urban growth by transforming natural and agricultural land into urbanized or bare soil. Land speculation and the agricultural crisis suffered in the Mediterranean countries since the beginning of the twentieth century has motivated this circumstance. A high-level growth of the LST and SUHI has also been evidenced from the year 2005, producing a deterioration of thermal comfort and an increase in the areas classified as UHS. For the Turkish Republic of Northern Cyprus, it is essential to identify and monitor these areas to implement actions and measures to improve these unfavorable situations for the population.
In the field of practice, our results provide further scientific evidence of the relationship already identified between the variables investigated. Therefore, the findings could be helpful to future sustainable growth; public administrations and urban planners could implement new mitigation and resilience measures that allow mitigating the increases in LST, SUHI and UHS and improve the inhabitants' quality of life through a decrease in the values of UTFVI.
Author Contributions Not applicable as there is only one author.
Funding Funding for open access publishing: Universidad de Granada/ CBUA. Not applicable.

Availability of Data and Material Not applicable.
Code Availability Not applicable.

Conflict of Interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.

Ethics Approval
The author indicates that all the ethical principles governing the publication of a research article in a journal have been followed.

Consent to Participate
The author agrees to participate in the review process and subsequent publication in the event of such an event.

Consent for Publication
If the article is accepted, the author consents to the publication and transfer of the information to the journal.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.