Tourism Effect on the Spatiotemporal Pattern of Land Surface Temperature (LST): Babolsar and Fereydonkenar Cities (Cases Study in Iran)

: The purpose of this study is to investigate the effect of tourism on Land Surface Temperature (LST), an issue which has rarely been considered in the tourism development literature. In this research, remote sensing techniques have been used to analyze the changes in the LST and spectral indices including the Normalized Difference Vegetation Index (NDVI), Modiﬁed Normalized Difference Water Index (MNDWI) and Enhanced Built-Up and Bareness Index (EBBI). The data used were based on Landsat Collection 1 Surface Reﬂectance (SR) images taken in June and August. They were analyzed over 32 years in the years 1987, 1993, 1999, 2009, 2014 and 2019. The study area included the cities of Babolsar and Fereydonkenar and their suburbs in Mazandaran Province in the north of Iran and south of the Caspian Sea. First the tourism zones were separated from other land use zones and then the changes in land use and LST in each of the zones were studied for each year based on the trend of 32-year change. The results of Pearson correlation in the whole area for each main land use zone showed that there was a signiﬁcant inverse relationship between the LST and the NDVI and MNDWI indices. This relationship was direct and signiﬁcant for the EBBI index. Moreover, the results of one-way analysis of variance (ANOVA) and Tukey test showed that the LST changes in the tourism zones during the study period were signiﬁcantly different from the other zones, so that the tourism zones always experienced lower LST. The ﬁndings also showed that, in the tourism zones, the values of the NDVI and MNDWI indices showed an increasing trend compared to the urban zone. Therefore, increasing the values of these indices due to the development of green space and its regular irrigation in tourism zones has led to a signiﬁcant decrease in the LST. The applied results of this research in the urban planning and tourism literature indicate that any model of physical development such as urban development does not necessarily lead to an increase in the LST, and this is entirely dependent on the physical design strategies.


Introduction
Climate change implication has been one of the primary research topics in the field of tourism that has attracted many researchers in recent years. The focus of such studies is on climate change's significant role in the tourism sector [1][2][3][4][5]. This is important because climate change vulnerabilities have negative consequences for the future of tourism, especially in areas such as Africa, the Middle East, South Asia, and the small island developing states [6]. In this regard, the dominance of studies on the effects of climate change in tourism has led to the neglect of the other side of this relationship, namely the effects of tourism on climatic parameters. Tourism has significant influence on climate parameters [7][8][9]. It is, mainly through transportation, a major contributor to greenhouse Land 2021, 10, 945 3 of 24 neglected so far in tourism literature; second, to provide scientific evidence from a developing country from which there is little information regarding tourism effects on climate parameters. This article consists of several sections: the introduction section describes LST and its relationship to built-up environments, and is focused on filling a research gap in tourism literature. In the research method section, data collection and analysis techniques are described. In the results and discussion section, the spatiotemporal changes of LST in the study area and its relationship with tourism spatial changes are reviewed and analyzed. Finally, the results are summarized and synthesized in the conclusion section.

Study Area
The study area (Figure 1), the coastal cities of Babolsar, Fereydonkenar and their hinterlands, are located in Mazandaran Province, in the North of Iran and South of the Caspian Sea. These cities, with a population of about 100,000 people [73], have been one of the most popular domestic tourism destination in Iran for the past 50 years due to their proximity to the Caspian Sea, their well-developed tourism services and their hospitable communities. According to the results of a survey in Mazandaran Province, Babolsar is the second most attractive destination for domestic tourism after Ramsar [74,75]. As coastal cities with agricultural hinterlands, they have experienced increasing demand for tourism over recent decades. According to survey estimates, the area was the destination for nearly six million domestic tourists in 2017 [74]. The most important attractions in this area according to the survey are the sea, the climate and the variety and abundance of shopping centers. Using important climatic parameters such as temperature, relative humidity, wind, radiation and bio-climatic indicators presented by Baker et al., the climatic features of Babolsar and Fereydoonkenar show that, from May to November, Babolsar has optimal natural conditions for tourist activities in open space in terms of thermal comfort. Although in the two months of July and August due to the relative increase of temperature and high relative humidity, a sultry state prevails in these cities, the wind blowing, especially in the coastal area, make this situation tolerable and even optimal [76].  The growing trend of tourism in the study area over the past ten years has been such that the area has changed from a mostly agricultural area to a multi-functional area focusing on tourism and university education services [74].
The growth of tourism in the area in the two main forms of overnight stay and second homes has led to significant spatial changes. Tourism-induced spatial changes are mainly caused by the change of the coastal areas and bare land into second homes, shopping malls, hotels, restaurants, entertainment centers, urban public spaces and development of infrastructure such as roads, water and electricity [74].
Before addressing the effects of tourism in the study area, it is important to describe the main types of tourism development. As mentioned earlier, Babolsar and Fereydonkenar have two main types of tourism, overnight tourists and second home tourists. The main services to the overnight tourists have been developed in the coastal zone of the area, including hotels, inns, restaurants, shopping and entertainment centers. However, second home tourists are divided into two different categories, detached second homes, and attached second homes. The detached second homes have been developed mainly in the western wards of Babolsar, detached from the main urban area; their main feature is low density as well as vast yards and green space. In contrast, the attached second homes with higher density have been developed in the northern area and are attached to the urban zone.

Data and Methods
In this research, the relationship between the spatial changes in tourism and changes in LST has been investigated using remote sensing techniques ( Figure 2). It should be noted that all the analyses related to the land use/cover and LST changes were done during 32 years and in 6 periods, i.e., in the years 1987,1993,1999,2009,2014 and 2019 (Table 1).

Data and Methods
In this research, the relationship between the spatial changes in tourism and changes in LST has been investigated using remote sensing techniques ( Figure 2). It should be noted that all the analyses related to the land use/cover and LST changes were done during 32 years and in 6 periods, i.e., in the years 1987,1993,1999,2009,2014 and 2019 (Table  1). To classify the images and calculate the main spectral indices and LST (Figure 3), the data were obtained from Landsat Collection 1 SR images taken in June and August (Table  1). SR improves comparison between multiple images taken from the same area by con- To classify the images and calculate the main spectral indices and LST (Figure 3), the data were obtained from Landsat Collection 1 SR images taken in June and August (Table 1). SR improves comparison between multiple images taken from the same area by considering atmospheric effects such as aerosol scattering and thin clouds. This can significantly contribute to the detection and characterization of the land surface change [77,78]. In addition, the images have been selected carefully to observe two important points: (1) images with the least cloudiness and (2) images belonging to summer, because in this season the study area has the least cloudiness and the vegetation of the region is in its optimal condition in terms of components affecting growth (temperature and humidity). According to the physical characteristics of the area, all the land uses/cover was classified into three main categories: vegetation, built-up and bare lands. sidering atmospheric effects such as aerosol scattering and thin clouds. This can significantly contribute to the detection and characterization of the land surface change [77,78]. In addition, the images have been selected carefully to observe two important points: (1) images with the least cloudiness and (2) images belonging to summer, because in this season the study area has the least cloudiness and the vegetation of the region is in its optimal condition in terms of components affecting growth (temperature and humidity). According to the physical characteristics of the area, all the land uses/cover was classified into three main categories: vegetation, built-up and bare lands.  [80]; Enhanced Built-Up and Bareness Index (EBBI) is able to map built-up and bare land areas using a single calculation. The EBBI is the first built-up and bare land index that applies near infrared (NIR), short wave infrared (SWIR), and thermal infrared (TIR) channels simultaneously [81]; Proportional Vegetation (PV) gives the estimation of area under each land cover type. The vegetation and bare soil proportions are acquired from the NDVI of pure pixels [82]; Land Surface Emissivity (LSE) is an essential parameter for retrieving LST. The LSE indicates the actual ability of an object to emit energy relative to a black object with the same temperature [83].
Based on the reflectance characteristics and the NDVI, MNDWI, EBBI as well as the LST indices, the satellite images of the study area were finally classified ( Figure 4). In this Figure 3. Calculation of the NDVI, MNDWI and EBBI and LST indices from satellite images. Normalized Difference Vegetation Index (NDVI), also called the standardized vegetation index, is a comprehensive reflection of vegetation type, coverage form, and growth conditions in unit pixel [79]; Modified Normalized Difference Water Index (MNDWI) uses green and SWIR bands for the enhancement of open water features. It also diminishes built-up area features that are often correlated with open water in other indices [80]; Enhanced Built-Up and Bareness Index (EBBI) is able to map built-up and bare land areas using a single calculation. The EBBI is the first built-up and bare land index that applies near infrared (NIR), short wave infrared (SWIR), and thermal infrared (TIR) channels simultaneously [81]; Proportional Vegetation (PV) gives the estimation of area under each land cover type. The vegetation and bare soil proportions are acquired from the NDVI of pure pixels [82]; Land Surface Emissivity (LSE) is an essential parameter for retrieving LST. The LSE indicates the actual ability of an object to emit energy relative to a black object with the same temperature [83].
Based on the reflectance characteristics and the NDVI, MNDWI, EBBI as well as the LST indices, the satellite images of the study area were finally classified ( Figure 4). In this regard, first, the educational points with suitable distribution were selected for three main land use/cover classes (built-up, bare land, vegetation). For this purpose, as shown in Table 2, the MNDWI, NDVI and EBBI spectral indices were used to identify the training points of water bodies, vegetation, built-up and bare lands, respectively [78][79][80]. Then, the satellite images were classified using a supervised classification method and Support Vector Machine (SVM) Algorithm. The results of the classifications were evaluated using 130 well-distributed land sampling points based on the images belonging to the year 2019 ( Figure 3). In this regard, the overall accuracy and kappa coefficient were obtained and the results confirmed the acceptable accuracy of the method used (overall accuracy = 96.07 and kappa coefficient = 0.93).
the results confirmed the acceptable accuracy of the method used (overall accuracy = 96.07 and kappa coefficient = 0.93).  After identifying and classifying the various classes of land use/cover, the LST characteristics and its changes were studied separately for each land use class over the six time periods. Then, in order to investigate the probability of a significant difference between the 32-year averages of the LST of the three land use classes, vegetation, built-up and bare  Table 2. Different types of land use/cover classes and related indices in the studied area.

LULC Description Identification Reference
Vegetation All areas covered with green space including agricultural lands, parks, urban green spaces and forests NDVI > 0.35 [84] Water All water-covered areas (sea, river and dam) MNDWI > 0 [80] Built Up Man-made lands including city, village and related infrastructure 0 < EBBI < 0.35 [81] Bare Lands Bare lands EBBI ≥ 0.35 [81] After identifying and classifying the various classes of land use/cover, the LST characteristics and its changes were studied separately for each land use class over the six time periods. Then, in order to investigate the probability of a significant difference between the 32-year averages of the LST of the three land use classes, vegetation, built-up and bare lands, one way ANOVA and Tukey technique were used. Moreover, to examine the relationship between the LST and the corresponding values of the NDVI, MNDWI and EBBI indices in each cell, Pearson correlation method was used for each year. In other words, the relationship between the type of land use/cover and the LST was investigated for each cell.
In order to investigate the relationship between the spatial changes in tourism and changes in the LST, first based on the 2019 images, the area of land uses related to the tourism sector, including second villa homes, second apartment homes and tourism services, were separated from the other land uses, especially urban spaces. Due to the relatively similar physical characteristics of the tourism and urban space, in order to separate these two zones two fieldwork steps were taken.
In the first step, eight real estate dealers and investment consultants with more than 10 years of experience in the area were selected using the snowball sampling technique and they were interviewed through open interview questions. The interviewees were asked to separate tourism spaces from the urban space according to their experience and in accordance with the map (Scale = 1:10,000). In order to verify the opinion of the real estate consultants, they were also asked about the address of the streets separating the tourism blocks. As mentioned, the development pattern of second homes in the study area includes second villa homes detached from the urban area and second apartment homes attached to the urban area. Thus, during the interview process, the two spaces were separated from each other. The experience and high mastery of real estate consultants, who even participated in many tourism projects in the region, provided the research team with accurate information on identifying areas for tourism development and their historical trends.
After receiving the views of the real estate consultants and confirmatory opinions from the last interviewed consultants, in order to control and verify the interviewees' ideas, the plotted areas were re-examined by the research team. In the study area, it was relatively easy to distinguish second homes from native urban wards since the first preference of tourists is not to interfere with the native urban wards. Therefore, tourism spaces are separated from the native spaces using the streets as borders and sometimes fencing. On the other hand, the architecture of tourism spaces with the characteristic of more durable materials, modern design, use of new technology in construction and greater area, can be distinguished from the main urban wards. At the same time, in the case of blocks in which there was some interference between the urban and tourism spaces, if more than 70% of the land use in the block was dedicated to tourism activities or second homes, then that block would be included in the tourism zone. It is noteworthy that only a few blocks had this condition. However, due to the considerable physical distance from the main urban area, fencing and even the presence of security at the entrances, the identification of the detached tourism zone was made with great care, relying on the opinion of real estate consultants and field examinations. At the end of this stage, a vector border layer was created for all land uses in the form of four main zones: the urban zone, vegetation zone, detached second home zone (detached from urban area) and attached second home zone (attached to urban zone). This vector layer was used to analyze all pixels along with the LST values and land use characteristics in the next steps.
The area of different land use/cover classes (vegetation, built-up and bare lands) and LST were re-calculated separately for each of the four main zones (urban zone, vegetation zone, detached second home zone and attached second home zone). These calculations were repeated for all six time periods to obtain land use/cover changes in each zone over time. Thus, for each image pixel, the values related to the land use class, zone name, values of each of the NDVI, MNDWI and EBBI indices and finally the LST were recorded. Since Landsat images were presented in 30-m pixels (thermal images were also resized with 30-m pixels), the study area included 207,870 pixels. The values of each pixel were recorded in a table for statistical analysis. Thus, it was possible to run statistical analyses on the LST value and its relationship with the spectral indices at the level of each pixel, separately for the four main zones. The linear trend (α) formula presented below was also used to calculate the annual rate of change of the NDVI, MNDWI, EBBI and LST time series for each pixel.
where α is the slope of the line, β is a constant, n is the number of years, y indicates the numerical value of the indices and the LST and x indicates the years studied. After calculating the annual rate of change of the NDVI, MNDWI, EBBI and LST at the level of each cell separately for the four main zones, Pearson correlation coefficient was used to investigate the relationship between the LST changes and changes in the spectral indices. The purpose of this stage was to investigate in more detail the effect of spatial changes in tourism on changes in the LST. To further clarify this issue, one-way ANOVA and Tukey technique were used to investigate the significant differences in the LST changes and the NDVI, MNDWI and EBBI indices in the four main zones. The results of this section helped the researchers to investigate the effect of the spatial changes in the attached and detached tourism zones on changes in the spectral indices and the LST, and compare the results with the effect of the spatial changes in the urban zone and vegetation zone.
In addition to statistical analysis, according to the relationship between the mean and standard deviation of the annual change rate, the spatial and temporal patterns of the LST changes in the study area were divided into three categories: usual temperature changes, unusually high or low temperature changes, very unusually high or low temperature changes. In this regard, the values between X ± σ were considered as usual changes, values greater than X + σ and smaller than X + 2σ were considered as unusually high changes and values greater than X + 2σ were considered as very unusually high changes. On the other hand, values smaller than X − σ and greater than X − 2σ were considered as unusually low changes and values smaller than X − 2σ were considered as very unusually low changes.

Land Use Changes in the Study Area
In this study, in order to investigate the possible impact of tourism on changes in the LST, first land use changes in the area were analyzed over a period of 32 years (in six time periods) from 1987 to 2019. In this regard, first the land use/cover changes in the area were analyzed in three classes, i.e., built-up, bare lands and vegetation, and then, by separating the area of tourism spaces from other land uses, the effect of the spatial changes in tourism on land use changes and changes in the LST were examined.After the easing of economic and security unrest caused by the Iran-Iraq war, from 1987 onwards the entry of domestic tourists, both as overnight and second-home tourists, gradually increased in study area [74]. In the development process of tourism spaces, proximity to the coastal area and relative distance from native settlements have been the criteria for investors of second homes and tourism services [74]. On the other hand, due to the relative salinity (unsuitable for agriculture), first bare lands with a suitable location were transformed into projects for the construction of second home complexes and tourism services. Between 1987 and 2019, the area of bare land decreased from 49 to 35 km 2 with a negative annual change rate of 2.53 km 2 . Vegetation coverage also decreased with a negative annual change rate construction of second home complexes and tourism services. Between 1987 and 2019, the area of bare land decreased from 49 to 35 km 2 with a negative annual change rate of 2.53 km 2 . Vegetation coverage also decreased with a negative annual change rate of 2.41 km 2 , while the area of built-up land use increased with an annual change rate of 4.62 km 2 from 9.5 to 32.26 km 2 . This shows that the area of built-up land use has increased by about 340% from 1987 to 2019 ( Figures 5 and 6).

The Impact of Land Use Changes on LST Changes
At this stage of the research, to investigate the relationship between the la use/cover and the LST, the LST changes associated with each of the land use classes (bu up, bare land and vegetation) were examined for six time periods from 1987 to 2019 years). For this purpose, first the six-period average of the LST was calculated separat for each pixel and then the results were compared based on the three main land use c ses. The results show that, in general, the average LST of all land use classes displaye decreasing trend from 1987 to 1999. However, from 1999 to 2019, they increased. In meantime, due to shading and evapotranspiration processes, the vegetation class has b cooler in all the years under study with a difference of 5 °C between the built-up and b lands (Figures 7 and 8).
The results of ANOVA test showed that there was a statistically significant differe between the mean of the LST of the studied land use classes (1987-2019) (F= 14.033, 0.0001). In a more detailed investigation, the results of the Tukey HSD test also show that the LST values in the vegetation class significantly differed (difference ≃ 5.5 °C) fr the LST in the bare lands and built-up area during the six periods studied (p = 0.00 However, this difference was insignificant between the bare lands and the built-up a (p = 0.998) (Tables 3 and 4). In this regard, the highest values of the LST were observed the bare lands in 1987 and 1999, while in 2014 and 2019 the built-up lands had the high LST values (Figure 7). The bare lands became warmer in 1987 and 1999 due to their sa cover, which had to support construction in the following years, with the developmen tourism, especially in the coastal lands. At the same time, a large part of the bare lan identified in the images of 2014 and 2019 were related to farmlands. The crops were h vested earlier than planned due to the earlier start of the hot season.

The Impact of Land Use Changes on LST Changes
At this stage of the research, to investigate the relationship between the land use/cover and the LST, the LST changes associated with each of the land use classes (built-up, bare land and vegetation) were examined for six time periods from 1987 to 2019 (32 years). For this purpose, first the six-period average of the LST was calculated separately for each pixel and then the results were compared based on the three main land use classes. The results show that, in general, the average LST of all land use classes displayed a decreasing trend from 1987 to 1999. However, from 1999 to 2019, they increased. In the meantime, due to shading and evapotranspiration processes, the vegetation class has been cooler in all the years under study with a difference of 5 • C between the built-up and bare lands (Figures 7 and 8).     3.673 * The mean difference is significant at the 0.01 level.
As described in the research method section, in the next stage the correlation between the LST values and the corresponding values of the NDVI, MNDWI and EBBI indices were calculated at the pixel level ( Table 5). The results of this stage showed that, in all the six time periods studied, there was a direct and significant correlation between the LST and EBBI index, indicating that the LST increases with increases in built-up and bare land areas. Due to the decrease of bare land area during the studied periods, it can be found  (Figure 7). The bare lands became warmer in 1987 and 1999 due to their sand cover, which had to support construction in the following years, with the development of tourism, especially in the coastal lands. At the same time, a large part of the bare lands identified in the images of 2014 and 2019 were related to farmlands. The crops were harvested earlier than planned due to the earlier start of the hot season.  As described in the research method section, in the next stage the correlation between the LST values and the corresponding values of the NDVI, MNDWI and EBBI indices were calculated at the pixel level ( Table 5). The results of this stage showed that, in all the six time periods studied, there was a direct and significant correlation between the LST and EBBI index, indicating that the LST increases with increases in built-up and bare land areas. Due to the decrease of bare land area during the studied periods, it can be found that the positive correlation between the LST and EBBI is mainly due to the built-up area expansion. In contrast, the results showed that there is a significant inverse relationship between the LST and the NDVI and MNDWI indices, which indicate the vegetation and water bodies, respectively. This study shows that the correlation between the LST and NDVI changed from −0.4 in 1999 to −0.79 in 1987. The correlation values of these two variables during the six time periods studied show a significant inverse correlation. In contrast, the correlation between the LST and EBBI was significant and direct, so that in all the six time periods studied, the correlation between these two variables was at least from 0.69 to 0.88.

The Impact of Tourism Spatial Changes on LST
To investigate the effects of spatial development in tourism on the LST, the whole study area was divided into four zones including the urban zone, vegetation zone, the second home zone detached from the urban zone (with all surrounding tourism services), and the second home zone attached to the urban zone based on the images of 2019 (research method section). By separating the tourism areas, it was possible to study the impact of tourism on the spatial changes of the region and finally on the changes of the LST. Findings show that the built-up area in the two tourism zones increased from 1.32 km 2 to 9.32 km 2 from 1987 to 2019. Compared to the total built-up area in the region in these two years, tourism zones made up 14.57% and 30.18%, respectively. On the other hand, the built-up area in the urban zone increased from 5.92 to 14.74 km 2 from 1987 to 2019 which made up 65.59% and 47.74% of the total built-up area in the region, respectively. In other words, the built-up area in the two tourism zones increased by 708.1% from 1987 to 2019, while this increase was 248.8% in the urban zone. Therefore, the expansion of the built-up area in the tourism zones was more than 2.5 times as much as in the urban zone. This indicates the role of tourism in the occurrence of significant spatial changes in the study area (Table 6). In addition, the 32-year changes in the three classes of land use/cover, built-up, bare land and vegetation, in each of the four zones revealed more detailed information about the spatial changes within each zone (Figures 9 and 10). Findings show that the area of bare land in the attached tourism zone, with an annual change rate of −0.65, decreased from 4.62 km 2 in 1987 to 1.26 km, while in the detached tourism zone the annual change rate was −0.79. In contrast, the area of vegetation coverage in the two tourism zones increased. The vegetation area in the detached tourism zone increased by an annual growth rate of 1.98, from 0.76 km 2 in 1987 to 2.26 km 2 (Figure 9). In addition, the 32-year changes in the three classes of land use/cover, built-up, bare land and vegetation, in each of the four zones revealed more detailed information about the spatial changes within each zone (Figures 9 and 10). Findings show that the area of bare land in the attached tourism zone, with an annual change rate of −0.65, decreased from 4.62 km 2 in 1987 to 1.26 km, while in the detached tourism zone the annual change rate was −0.79. In contrast, the area of vegetation coverage in the two tourism zones increased. The vegetation area in the detached tourism zone increased by an annual growth rate of 1.98, from 0.76 km 2 in 1987 to 2.26 km 2 (Figure 9).    In addition, the 32-year changes in the three classes of land use/cover, built-up, bare land and vegetation, in each of the four zones revealed more detailed information about the spatial changes within each zone (Figures 9 and 10). Findings show that the area of bare land in the attached tourism zone, with an annual change rate of −0.65, decreased from 4.62 km 2 in 1987 to 1.26 km, while in the detached tourism zone the annual change rate was −0.79. In contrast, the area of vegetation coverage in the two tourism zones increased. The vegetation area in the detached tourism zone increased by an annual growth rate of 1.98, from 0.76 km 2 in 1987 to 2.26 km 2 (Figure 9).   In this stage of the research, in order to distinguish the effects of tourism on LST from other land uses, paired comparisons were made between the main zones. In this regard, first the correlation between the annual change rate of the LST and the annual change rate of the NDVI, MNDWI and EBBI indices was calculated in the studied zones during the six time periods (from 1987 to 2019) ( Table 7). The results of the Pearson correlation test show that there was a significant direct relationship between the trend of changes in the LST and the trend of changes in the EBBI index. However, this relationship was inverse in the trend of changes in the NDVI and MNDWI indices. This shows that the rising trend of the LST is directly related to the expansion of built-up and bare lands. In contrast, in areas with increasing vegetation and water body, the LST tended to decrease over the 32-year period ( Table 7 and Figure 11).
LST and the trend of changes in the EBBI index. However, this relationship was inverse in the trend of changes in the NDVI and MNDWI indices. This shows that the rising trend of the LST is directly related to the expansion of built-up and bare lands. In contrast, in areas with increasing vegetation and water body, the LST tended to decrease over the 32year period (Table 7 and Figure 11).  In the next step, using an ANOVA statistical test, the trend of changes in the LST and NDVI, MNDWI and EBBI indices were compared statistically between the main zones in the area. The results showed that the 32-year trend in changes in the LST was significantly different in all four main zones, with 95% confidence. This result is also true for the trend of changes in the spectral indices (Table 8 and Figure 12).  In the next step, using an ANOVA statistical test, the trend of changes in the LST and NDVI, MNDWI and EBBI indices were compared statistically between the main zones in the area. The results showed that the 32-year trend in changes in the LST was significantly different in all four main zones, with 95% confidence. This result is also true for the trend of changes in the spectral indices (Table 8 and Figure 12). development model had a higher EBBI rate on the one hand, and on the other hand the development of the urban zone on the agricultural lands and the neglect of green space development has caused a continuous decrease in the NDVI index. Therefore, during the study period, the LST of the urban zone has always been higher than that of the tourism zones.   The results show that the trend in annual changes in the LST is increasing in the whole study area, but its values show a significant difference between the main zones, so that the increasing trend of the LST in the urban and vegetation zones was higher than that of the tourism zones. This indicates that the nature of spatial development of tourism, although seemingly similar to the nature of urban development, has an inverse effect on the LST in the tourism zones. It is noteworthy that the record of higher LST in the vegetation zone, which mainly includes farmlands, is due to the bareness of these lands after harvest at the time of study. Therefore, the result of the analysis of the above statistical test on the difference in the LST of the vegetation zone and tourism zones is invalid.
At the beginning of the research, with the aim of investigating the effect of tourism on the LST in more detail, two separate tourism zones were selected, as attached to and detached from the urban area. Thus, if the results are mixed with the urban zone, it is possible to separate these two types of tourism development model. Although the results of the paired comparisons show a significant difference in the changes in the LST and the spectral indices between the two tourism zones, a comparison of the results of the Tukey test makes it clear that both tourism zones have a different effect from that of the urban zone on changes in the LST. By examining the impact of tourism on the spectral indices (Table 9), the reason for this difference becomes clear. According to what has been said, tourism has moderated the LST in tourism zones both through the development of sandy bare lands and the development of green spaces (increase of the NDVI index). In contrast, in the urban zone, the NDVI index with a negative change rate and the EBBI index with a higher positive change rate show a significant difference to the tourism zone. In contrast to tourism development, the urban development model had a higher EBBI rate on the one hand, and on the other hand the development of the urban zone on the agricultural lands and the neglect of green space development has caused a continuous decrease in the NDVI index. Therefore, during the study period, the LST of the urban zone has always been higher than that of the tourism zones.
In the following, in addition to the above-mentioned statistical tests, first, the annual trend of the LST changes in the whole area was calculated (at the pixel level) to depict the spatial pattern of 32-year changes in the LST, and then the distribution pattern of the usual/unusual and very unusual values was shown (research method section). By adapting the boundaries of the four main zones with their trend values, the majority of pixels with very unusually low and unusually low values trend belonged to the tourism zones. This result shows that, during the 32-year period, with the development of tourism spaces on bare lands and the development of green spaces in the tourism areas, the LST experienced a declining trend with a significant difference from the other zones. In contrast, the majority of pixels with a very unusually high value (warmer) were seen in the urban zone ( Figure 13).
A more detailed study of the spatiotemporal pattern of the trend in LST changes in the study area shows that the pixels with the usual change value had a LST change between 0.057 to 0.181 • C per year. Peng et al. [85] also reported similar amounts of long-term LST change as usual in a study on Chinese cities.
In the study area, most of the pixels with usual LST changes included parcels with no land use change during the 32-year period, and therefore it can be said that they have become warmer in line with global warming. In contrast, pixels with unusually high and very unusually high (warmer) changes were divided into two groups. The first group included parcels with vegetation cover that went under construction between 1987 and 2019. These areas included parcels integrated in the urban zone or parcels with changed land use as a result of the development of infrastructures, e.g., Fereydonkenar ring road. Other pixels with unusually high LST changes included infrastructure development projects such as ports, piers, and offshore access routes (in the north at the coastal line) that have created new built-up areas inside the sea. A more detailed study of the spatiotemporal pattern of the trend in LST changes in the study area shows that the pixels with the usual change value had a LST change between 0.057 to 0.181 °C per year. Peng et al. [85] also reported similar amounts of longterm LST change as usual in a study on Chinese cities.
In the study area, most of the pixels with usual LST changes included parcels with no land use change during the 32-year period, and therefore it can be said that they have become warmer in line with global warming. In contrast, pixels with unusually high and very unusually high (warmer) changes were divided into two groups. The first group included parcels with vegetation cover that went under construction between 1987 and 2019. These areas included parcels integrated in the urban zone or parcels with changed land use as a result of the development of infrastructures, e.g., Fereydonkenar ring road. Other pixels with unusually high LST changes included infrastructure development projects such as ports, piers, and offshore access routes (in the north at the coastal line) that have created new built-up areas inside the sea.
There were two areas with unusually and very unusually low (colder) LST changes. The first area includes bare lands which changed to vegetation coverage during the study period. These areas mainly included agricultural lands located in the east and southeast of Babolsar. The second area was located in the tourism zones (detached and attached zones). The bare lands in these zones were used for tourism-related infrastructure. These constructions included second homes, green space in the tourism zone, shopping and recreation centers. This finding is further evidence for the significant impact of spatial changes in tourism on changes in the LST.

Discussion
The main results of this study showed that tourism, especially second homes, over a period of almost 30 years, as a strong driving force, has caused widespread land use change in the study area. As a result, the spatial changes caused by the development of tourism, in subsequent years, have caused significant changes in the LST of the region.
Previous research has repeatedly shown that the amount of outgoing energy (Albedo and Outgoing Longwave radiation (OLR)) is strongly influenced by the physical characteristics of the land surface. Thus, any change in land surface characteristics will lead to changes in the properties of radiation and ultimately changes in temperature [35][36][37][38][39]. In this regard, the LST has been considered as an important indicator in environmental studies [93].
In this study, the results of the ANOVA test showed the significant effect of three main land use/land cover classes, including built-up, bare land and vegetation on the LST. Accordingly, due to evapotranspiration and shading processes, vegetation has a There were two areas with unusually and very unusually low (colder) LST changes. The first area includes bare lands which changed to vegetation coverage during the study period. These areas mainly included agricultural lands located in the east and southeast of Babolsar. The second area was located in the tourism zones (detached and attached zones). The bare lands in these zones were used for tourism-related infrastructure. These constructions included second homes, green space in the tourism zone, shopping and recreation centers. This finding is further evidence for the significant impact of spatial changes in tourism on changes in the LST.

Discussion
The main results of this study showed that tourism, especially second homes, over a period of almost 30 years, as a strong driving force, has caused widespread land use change in the study area. As a result, the spatial changes caused by the development of tourism, in subsequent years, have caused significant changes in the LST of the region.
Previous research has repeatedly shown that the amount of outgoing energy (Albedo and Outgoing Longwave radiation (OLR)) is strongly influenced by the physical characteristics of the land surface. Thus, any change in land surface characteristics will lead to changes in the properties of radiation and ultimately changes in temperature [35][36][37][38][39]. In this regard, the LST has been considered as an important indicator in environmental studies [93].
In this study, the results of the ANOVA test showed the significant effect of three main land use/land cover classes, including built-up, bare land and vegetation on the LST. Accordingly, due to evapotranspiration and shading processes, vegetation has a significant effect in reducing the LST by 5 • C. The results of the correlation test also confirmed a significant inverse relationship between the LST and NDVI and MNDWI indices. The correlation between the LST and NDVI during the six study periods also showed a strong relationship from −0.4 to −0.79. A review of similar research also shows that, first, the relationship between the LST and NDVI is significant and inverse. The value of this correlation has been reported to be −0.4 and more [100,101,107,110]. In contrast, the correlation between the LST and EBBI was significant and direct, so that in all six time periods studied, the correlation between these two variables was from at least 0.69 to 0.88. Peng et al. [85] have also confirmed such a strong correlation between these two variables in their research.
This result shows that the LST has increased with increases in the area of built-up and bare lands, while increasing the area of vegetation and water bodies has reduced the LST in the study area. These particular results is in line with previous research showing the significant effect of vegetation on the LST [37,[100][101][102][103][104][105][106][107][108][109][110]. In this regard, Rahman et al. [97], by modeling the spatial changes in the coastal lands and their effects on the LST characteristics, showed that physical changes due to urban development and built-up lands have a direct and significant effect on the LST, and this could pose a challenge for the residents of the studied areas in the future. The warmer built-up surfaces and the moderating role of vegetation on the LST were also confirmed in this study. Hellings et al. [99] referred to the impact of urban development and built-up space on increasing the LST in European cities, showing that green spaces and vegetation coverage reduce the LST by between 4 to 6 • C compared to built-up spaces and impermeable urban surfaces. The results of the study by Maheng et al. [37] of the city of Colombo in Sri Lanka show that a 30% increase in green space in a built-up urban environment moderates the average air temperature by 0.1 • C. In contrast, changing the use of suburban vegetation to urban spaces will increase the average temperature from 27.73 to 27.75 • C. Although this research has been carried out in different geographical conditions, its results are in line with the results of research from major European cities [99], Sri Lanka [37] and Indonesia [39] with temperate climates, showing a significant effect of vegetation in decreasing the surface temperature in different climates.
In this study, more detailed results were obtained by separating the tourism zone from the urban zone and vegetation. Unlike urban development models, the tourism development model in the study area has occurred in such a way that it did not increase the LST, but moderated it. The results of ANOVA and Tukey tests showed that there is a significant difference between tourism zones and other zones based on the LST, NDVI, MNDWI and EBBI indices, so that during the 32-year period, the tourism zones showed higher values of positive changes in the NDVI and MNDWI indices and lower values of changes in the EBBI index. This means that tourism in the study area has improved and has developed vegetation and water bodies in the tourism zones, but at the same time, the area of vegetation coverage has declined in the region. Therefore, the LST values in tourism zones have always been lower than in the other zones. This finding is fundamentally different from the findings of Chu et al. [94] in Hainan Island, China. By dividing the study area into two main zones, tourism and non-tourism, Chu et al. [94] showed that the development of tourism over a period of 20 years has a direct and significant effect on increasing the LST so that the average difference between the surface temperature of the tourism zone and other areas was more than 1.2 • C per year. The tourism development model in Hainan Island caused extensive destruction of vegetation in the region and replaced it with impermeable surfaces related to tourism structures. In contrast, the development of tourism in the Babolsar area, on the one hand, has occurred through the change in use of bare lands that have been covered with sand. Sandy surfaces, with lower specific heat capacity than built-up surfaces, have a greater effect on increasing the LST. Therefore, the development of bare lands to the built-up tourism spaces has led to a decrease in the LST. On the other hand, the development of tourism spaces in the area has led to the development of green spaces that are irrigated regularly (Figures 9 and 12). One of the most important marketing strategies of tourism investors in the study area is the development of well-designed green spaces within tourism projects. Accordingly, although the goal of investors is mainly advertising and beautification of tourism areas, this measure has affected the microclimate of the tourism zone by reducing the LST. In this regard, Goldblatt et al. [107] conducted a study of microclimate on the campus scale, citing the relationship between the NDVI, NDBI and LST. They found that, in comparison with the physical characteristics of the bare and built-up lands, developed vegetation within a campus has a significant effect on cooling the LST.
The findings of this study showed that tourism with a specific physical development pattern can have a significant effect on changes in the LST. However, it cannot be argued confidently that tourism development reduces the LST because, as Chu et al. [94] in their research on Hainan Island showed, that the effect of tourism on the LST depends on its physical development pattern. In the case investigated in this study (Babolsar and Fereydonkenar), tourism had a significant moderating effect on the LST in the tourism zones through developing green spaces that were constantly irrigated, and constructing on bare lands that were generally covered with sand. Although this study focused on the effects of the physical development of tourism instead of urban development, its results on the correlation between the land cover indices such as NDVI, MNDWI and EBBI with the LST are consistent with previous research [92][93][94][95][96][97][98][99][100][101][102][103][104][105][106][107][108][109][110].
The findings of this study are significant because, at the beginning of the research, there was no similar study in the relevant literature and, so far, only one similar piece of research has been published [94]. This research has tried to deal with the relationship between the physical changes and changes in the LST from another perspective. On the other hand, a large part of the research related to the relationship between land use/cover characteristics and LST has been conducted in temperate and humid regions with rich vegetation, and relatively less research has been done in countries with arid climates such as Iran [96]. In this regard, providing evidence from different contexts will help to enrich the relevant literature. This research faces limitations, which can be addressed in future research. Although the temperature data and the related indices were analyzed in six time periods (32 years) on a pixel level from the whole area, they all belonged to a limited period of summer when the vegetation and cloud cover of the area were in optimal condition. Therefore, it would be useful to obtain more data covering the other times of the year to enrich the analysis, though this may increase the volume of data and input of the analysis. Therefore, other researchers could contribute to the enrichment of the relevant literature by providing new evidence through conducting a similar research over longer periods and in different climates, according to different patterns of tourism development.

Conclusions
This study investigated the possible effect of spatial changes in tourism on the LST. Although several studies have been conducted to investigate the effect of land use/cover changes on the LST, only one study has been published before this article on the effect of tourism. Thus, more evidence and different research conditions and methods are required. In this research, in order to classify the satellite images, the Landsat Collection 1 SR images taken in June and August were analyzed during 32 years in six periods, i.e., the years 1987, 1993, 1999, 2009, 2014 and 2019. In line with the purpose of the research in the first stage, the spatial changes and then the LST changes were studied in the area over a period of 32 years and in the six time periods mentioned above. Then, the correlation between the LST changes and the changes in the NDVI, MNDWI and EBBI indices was examined in the area. In the second stage, with the aim of analyzing the effects of tourism on changes in the LST, based on the field method, the whole area was divided into four main zones including the urban zone, the vegetation zone, the tourism zone attached to the urban area and the tourism zone detached from the urban area. Then, using the data belonging to the period of 32 years, the LST changes and its relationship with the NDVI, MNDWI and EBBI indices were calculated separately for each zone. Finally, the LST values and their relationship with the spectral indices were statistically compared between the four main zones.
The results showed that tourism over a period of 32 years (1987-2019), along with a significant change in land use/cover, has played an important role in the spatial changes of the study area. Similarly, the results confirmed the significant correlation of LST changes in tourism zones. The results showed that the nature of the effect of tourism on the LST is different from the other types of physical change factors, such as urbanization. This means that, in the study area, tourism did not increase the LST. Instead, cooler temperatures were always recorded in the tourism zones compared to the other land use/cover zones. By examining the reasons for this in more detail, the findings showed that the development of tourism in the area with the development of green spaces has helped to increase the NDVI index values. On the other hand, since the green spaces in the tourism zones were regularly irrigated, this has helped to further decrease the temperature. The development of tourism in the bare lands in the coastal areas is also another reason for this phenomenon. The bare lands in the tourism zones were mainly covered with sand, which increased the LST due to low specific heat capacity. However, the replacement of the mentioned land use with built-up area or green spaces has led to a decrease in the LST in the tourism zones. The findings of this study showed that the factors affecting the spatial changes such as urbanization and tourism cannot necessarily have similar effects on the environmental indices such as the LST. In addition, not all tourism development patterns can have the same effects. According to the results of this research, a tourism development model that leads to the development of green space can lead to lower LST in the area.