Effect of Land Use/Cover Changes on Urban Cool Island Phenomenon in Seville, Spain

Abstract

This paper analyses Seville’s surface urban heat island (SUHI) phenomenon, comparing spatial and temporal patterns of land surface temperature (LST) during July 1987, 2000 and 2017. Landsat data captured throughout three July months were analyzed for the different years, techniques of geographic information systems, ecological variables and geospatial approaches and used to carry out the analysis. The results indicate that from 1987 to 2017, the averaged LST has increased by 9.1 °C in the studied area. The urban areas are colder than their surroundings, suggesting the role of baresoil and cultivated land in the reversal of the SUHI phenomenon. The results show that a fraction of green space has a high unstandardized coefficient (β) through the three time periods. A decreasing trend is also observed in the standardized β “fraction of impervious surface” in the three time periods. The linear regression analysis shows a negative relationship of mean LST with impervious surface fraction due to the presence of shadows projected by buildings, and a positive relationship with green space fraction caused by the influence of baresoil and cultivated land that inverts the LST behavior pattern. The study concludes that there is a need to implement SUHI mitigation strategies during the initial phases of engineering projects where the origin of this problem can be acted upon, since the process of creating streets and public space offers a valuable opportunity to restore the environmental quality and diminish the effects generated by climate change.

1. Introduction

Climate change is the greatest environmental challenge facing humanity today due to its global scale and profound social and economic implications. In 2015, during the United Nations Climate Change Conference in Paris, the first balanced global agreement was adopted to address global warming and the objective of limiting temperature increase to 2 °C by 2100 was agreed. In the words of Laurent Fabius, “a global climate agreement is a universal necessity that must be undertaken by all countries thus promoting climate solidarity and supporting the mobilization of financing and technological development” [1].

Cities only represent the 2% of the earth´s surface [2]; however, 60% of the word´s energy consumption and more than the 70% of carbon dioxide emissions take place in cities [2,3,4]. Therefore, cities contribute significantly to climate change.

Furthermore, there is a global trend to migrate from rural to urban areas. This migratory flow, much more notorious in developed countries, represents a social, economic and technical challenge. According the United Nations, by 2050, the population in cities will increase by 2.5 billion [5].

In recent decades because of rapid urbanization, biodiversity and several natural habitats [6,7], as well as soils, destined to agriculture and forests have been lost. There has also been an increase in the temperature in cities closely related to urban air quality, energy consumption and the health of at-risk populations [8,9,10,11,12,13,14,15]. All this shows the conclusion of the important protagonist that human activities play in the alteration and modification of the climate system [16,17,18,19,20,21,22]. In this situation, it is essential to know and understand the role and impact of different land uses on the environmental system [23,24].

Urban heat island (UHI) is a thermal phenomenon in which the temperature of urban spaces is higher than nearby rural zones [12,20,22,25,26,27,28,29,30,31]. In 1818, Luke Howard was the first researcher to discover London’s thermal variation, pointing out that the temperature of the city centre was 3.7 °F higher than the surrounding fields [25]. Subsequently, in 1958 Gordon Manley was the first to define this thermal variation as an urban heat island [32]. The interest in the study of the urban heat island is well justified, since urban areas constitute unique sectors within the climate of the region where they are located. This interest responds not only to the need for knowledge to achieve a more pleasant environment for city dwellers, but also to the need to analyze and anticipate the changes these entail and the repercussions they may have on the climate.

There are two types of UHI: the first is based on the fact that air temperature is higher during the night, while the second is based on the fact that land surface temperature (LST) tends to be higher during the day due to solar radiation [27].

The formation of a surface urban heat island (SUHI) is principally caused by changes in the landscape due to urban development, mainly in the reduction of vegetation in urban areas, the properties of urban materials used, urban geometry, anthropogenic heat, and climate and location, all of which generate an increase in LST [26,33,34,35,36].

LST is a key variable in SUHI generation and constitutes a determining factor in surface radiation and energy exchange [37], as well as the control of heat distribution between the surface and atmosphere [38,39]. Therefore, it is necessary to model and predict environmental changes, as well as to analyze and comprehend the dynamics of the LST and its connection with changes of anthropogenic origin [40].

The composition of land use and land cover (LULC) is one of the key factors influencing LST [41,42,43]. LST generally has a positive correlation with an impervious surface and a negative correlation with areas of vegetation [39,44].

Currently, it is still difficult to identify the distribution and spatial pattern of SUHI with temperature data observed in situ because of the lack of total coverage of the studied zone and the limited spatial resolution of the data [45,46,47].

Remote sensing provides a very promising approach for the generalized study of LST and therefore makes it possible to monitor the spatial outline of the SUHI. The high spatial resolution thermal satellite data obtained during the daytime have been extensively used to identify and evaluate SUHI at the mesoscale, i.e., in the space occupied by the entire city that usually extends between 10 and 200 km [29,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66].

Among the available satellite data, Landsat data have been widely used in numerous SUHI case studies around the world thanks to their precision, such as the detection capability in terms of the spatial and temporal resolutions with the uncertainties and due to the fact that data acquisition is freely accessible with spatial and temporal coverage for most of the areas susceptible to develop SUHI [60,63,67,68]. On the other hand, the integration of remote sensors and geographic information systems (GIS) has been recognized as a powerful and effective tool to detect the different uses of urban land and changes in its coverage which are decisive in the understanding of the relationship between human activities and the SUHI [69,70,71,72].

The aim of this research is to establish the relationship between the composition and pattern of the urban landscape and the formation and evolution of the SUHI in the Seville City. Various landscape variables are explored to analyze the spatial and temporal variations of the land surface temperature in the analyzed area. The land cover groups considered in this analysis were: impervious surface (including surfaces found in urban and suburban landscapes such roads, parking lots, drive ways, sidewalks, roofs and industrial areas), green space (GS) (including land that is covered with grass, trees, shrubs, or other vegetation), water (W) (including all bodies of water) and other (BSC) (including all the lands not classified as green space, impervious surfaces, and water. This category mainly includes baresoil and cultivated land).

The methodology and the results obtained from this research will constitute a powerful tool that will allow researchers to know which zones are more vulnerable to the formation of the SUHI and why. This research establishes the starting point of a wider and more ambitious study, in which the factors that influence the formation of the SUHI for different cities in the world will be analyzed. By means of the comparative analysis of the results obtained from each of the analyzed cities, it will be determined to what extent the strategies carried out for the reduction and mitigation of the SUHI have been effective, and to conclude that it is essential to implement these strategies in the design and initial phases of the engineering projects to improve the adaptation of cities to the climate change and increase their resilience.

2. Methodology

2.1. Study Area

Seville is the largest and most populated province in Andalusia, with an approximate area of 14,000 km² and 1.9 million inhabitants [73]. It is bordered to the north by Badajoz, to the east by Cordoba, to the south by Cadiz and Malaga, and to the west by Huelva. It is surrounded by the Sierra Morena mountain range to the north, by the Sierras Subbéticas mountain range to the south-east and between them is the Guadalquivir river in the middle and lower portion of the basin, where the valley opens up towards the Atlantic forming a wide area of countryside and marshes [74]. The climatic spectrum of Seville City is very broad, ranging from the subtropical climate in the Guadalquivir basin to the temperate Mediterranean climate of humid winters and long hot summers, and to the areas of the Sierra Norte to the south of the province with harsher climatic characteristics. The average annual temperature is 18.6 °C; with July being the warmest month with an average of 27.8 °C and January being the coldest month with an average of 10.3 °C. The average rainfall is 576 mm; July is the driest month with an average of 1 mm and November is the wettest month with an average of 87 mm [75]. The study area encompasses 40,000 ha, spanning the city of Seville and part of its adjacent municipalities. The elevation of the study area varies between −9.97 m and 163.95 m, with an average elevation of 30.96 m above sea level. The slope varies from 0 to 64.88 degrees, with an average slope of 2.51 degrees (Figure 1). This study area includes built-up areas, arable land, grassland areas and shrubland.

2.2. Satellite Data Used and Pre-Processing

For this study, Landsat satellite images captured in 1987; July 01; GMT 10:27:27 (Landsat 5 TM), 2002; July 02; GMT 10:50:59 (Landsat 7 ETM+) and 2017; July 03; GMT 11:02:20 (Landsat 8 OLI/TIRS) have been used (http:/earthexplorer.usgs.gov/). For the selection of the satellite images, those that were cloud-free or with minimum cloud coverage (less than 10%) were considered. The study area is contained entirely within path 202 and row 34. All images were acquired in the same month, during the dry season (Figure 2).

In order to carry out land cover mapping, image classification, index derivation and LST retrieval, all satellite images were subjected to two pre-processing procedures consisting of radiometric calibration and atmospheric correction using ArcMap software. The atmospheric correction was performed using the dark-object subtraction (DOS) model proposed by [76]. This model postulates that atmospheric mist increases the digital number (DN) value in areas of clean, deep and calm water, where the physical characteristics must have zero reflectance. The representative value of that difference is subtracted in each band, in all the pixels of the scene. To estimate LST it is necessary to convert the DN from each of the thermal bands to radiance values as a measure of the amount of energy that reaches the satellite [37,77,78]. These radiance values are then used to perform the conversion to surface brightness temperature (expressed in Kelvin degrees), considering emissivity, vegetation fraction, the normalized vegetation index, and calibration constants.

2.3. Land Cover Mapping

For the years 1987, 2002 and 2017, the Landsat imageries were obtained in order to obtain the classification of the land cover maps. For this purpose, the maximum likelihood supervised classification process was used (Figure 3) [36,79]. The following land cover groups were considered in this analysis:

• Impervious surface (IS), including surfaces found in urban and suburban landscapes such roads, parking lots, driveways, sidewalks, roofs and industrial areas.
• Green space (GS), including land that is covered with grass, trees, shrubs, or other vegetation.
• Water (W), including all bodies of water.
• Other (BSC), including all the lands not classified as green space, impervious surfaces, and water. This category mainly includes baresoil and cultivated soil.

The bodies of water and the impervious surfaces were extracted from the images using the modified normalized difference water index (MNDWI) (Equation (1)) and the visible red and NIR-based built-up index (VrNIR-BI) (Equation (2)), respectively [80,81,82,83,84].

The green spaces were extracted from the images via the normalized difference vegetation index (NDVI) (Equation (3)) [85].

$$\mathrm{MNDWI}=\frac{{\mathsf{\rho }}_{\mathrm{Green}}-{\mathsf{\rho }}_{\mathrm{SWIR}1}}{{\mathsf{\rho }}_{\mathrm{Green}}+{\mathsf{\rho }}_{\mathrm{SWIR}1}}$$

(1)

$$\mathrm{VrNIR}-\mathrm{BI}=\frac{{\mathsf{\rho }}_{\mathrm{Red}}-{\mathsf{\rho }}_{\mathrm{NIR}}}{{\mathsf{\rho }}_{\mathrm{Red}}+{\mathsf{\rho }}_{\mathrm{NIR}}}$$

(2)

$$\mathrm{NDVI}=\frac{{\mathsf{\rho }}_{\mathrm{NIR}}-{\mathsf{\rho }}_{\mathrm{Red}}}{{\mathsf{\rho }}_{\mathrm{NIR}}+{\mathsf{\rho }}_{\mathrm{Red}}}$$

(3)

where ρ_Green, ρ_Red, ρ_NIR, and ρ_SWIR1 refer to the surface reflectance values of bands 2, 3, 4 and 5 for Landsat 5 TM and Landsat 7 ETM+ and bands 3, 4, 5, 6 for Landsat 8 OLI/TIRS.

In order to verify the accuracy obtained in the classification of each of the land uses and land cover, different bands were combined (

Table 1. Composition of Bands for Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 OLI/TIRS satellites.

Composite Bands	Landsat 5 TM	Landsat 7 ETM+	Landsat 8 OLI/TIRS
Vegetation	4, 3, 2	4, 3, 2	5, 4, 3
Urban areas	7, 5, 3	7, 5, 3	7, 6, 4
Agriculture	5, 4, 1	5, 4, 1	6, 5, 2
Water	4, 5, 3	4, 5, 3	5, 6, 4
Natural colour	3, 2, 1	3, 2, 1	4, 3, 2

) and the information obtained from the MNDWI, VrNIR-BI and NDVI indices was used as a reference. A total of 1000 reference points generated by the stratified random sampling technique were used [86]. Finally, an overall accuracy of 86% was obtained.

2.4. LST Retrieval

To retrieve the LST from remote sensing satellite thermal bands, the NDVI were used to derive the emissivity of the land surface. Subsequently, using these values, the at-satellite brightness temperatures were scaled [87]. The pre-processed thermal bands created in Section 2.2 were used as follows in Equation (4) [88,89]:

$$\mathrm{LST}=\frac{{\mathrm{T}}_{\mathrm{B}}}{1+\left(\mathsf{\lambda }\times \frac{{\mathrm{T}}_{\mathrm{B}}}{\mathsf{\rho }}\right)\times \ln \mathsf{\epsilon }}$$

(4)

where T_B = at-satellite brightness temperature in degrees Kelvin; λ = wavelength of emitted radiance (λ = 11.5 µm for Landsat 5 TM band 6 and Landsat 7 ETM+ band 6 [89,90], and 10.8 µm for Landsat 8 OLI/TIRS band 10 [77]; ρ = h $\times c/\sigma$ (1.438 $\times$ 10⁻² mK), $\sigma$= Boltz-mann constant (1.38 $\times$ 10⁻³⁴ Js), and c = velocity of light (2.998 $\times$ 108 m/s); and $\epsilon$ is the land surface emissivity estimated using the NDVI method [87]. The resulting LST values were then changed from degrees Kelvin to degrees Celsius (°C).

2.5. Spatial Analysis

2.5.1. Urban-Rural Gradient Analysis

This analysis determines the spatial variability of LST, the spatial distribution of impervious surface, green space and baresoil and cultivated soil across the urban-rural gradient of the study area. Firstly, Seville’s city centre was located and called kilometre 0. Multiple ring buffer zones were then created around the centre of each study area defined by the time period analyzed, with an interval of 300 m distance. Mean LST and the density of impervious surface, green space, baresoil and cultivated land were obtained for each of the rings. For each of the rings, the number of pixels of each of the land uses was counted, and each pixel had dimensions of 30 × 30. Multiplying the number of pixels by the area of each pixel, we obtained the total area of each land uses in each of the rings. The fraction of each of the uses was obtained by establishing the relationship between the total area of each of the land uses in each ring and the total area of the ring. Finally, mean LST and the density of impervious surface, green space, baresoil and cultivated land were obtained for each of the rings.

2.5.2. Multiresolution Grid-Based Analysis

This study focuses on determining the effect of impervious surface, green space and baresoil and cultivated land on LST. For this purpose, a set of polygonal grids adjusted to LST raster maps was generated, considering Seville’s city centre calculated in Section 2.5.1. The size of each grid was 210 m × 210 m [33,77,91]. Finally, the mean LST and the density of impervious surface, green space and baresoil and cultivated land were obtained for each of the grids.

2.5.3. Statistical Analysis

By means of bivariate correlation analysis and scatter diagrams, the correlation between mean LST and the density of impervious surface, green space, baresoil and cultivated areas measured in each buffer zone (urban-rural gradient analysis) and grid (analysis based on multiple resolution grid) were determined.

To regress the correlation between mean LST and the independent variables considered, in each of the 9216 polygonal grids a multiple linear regression (MLR) by ordinary least squares (OLS) method was used [92,93,94,95].

The OLS MLR model is given as Equation (5):

$$Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\dots +{\beta }_k{X}_k+\epsilon$$

(5)

where Y is the dependent variable (mean LST), ${\beta }_0$ is the intercept coefficient, ${\beta }_1+{\beta }_2,\dots ,{\beta }_k$ are the regression coefficients or slopes, $X_1+X_2,\dots ,X_k$ are the explanatory variables, and $\epsilon$ is the standard error.

Seven explanatory parameters were selected based on the hypothesis of the significant influence they have on the variations of LST in the study area. These variables include the fractions of IS, GS, and BSC; mean elevation, mean slope, mean aspect and mean hill shade. Both the dependent variable (mean LST) and the seven explanatory parameters were estimated on the set of polygonal grids used for the analysis described in Section 2.5.2. The fractions of impervious surface, green space and baresoil and cultivated areas that were considered were obtained from the land cover maps developed in Section 2.3. The remaining four explanatory topographic variables were obtained from the digital elevation model using ArcMap software. It has been obtained by the interpolation of the land class obtained from LIDAR flights 1st Coverage of the National Plan for Aerial Orthophotography (PNOA).

3. Results

3.1. Impervious Surface, Green Space and Baresoil and Cultivated vs. LST

Land cover classification maps indicate that the study area has undergone rapid urbanization (Figure 3). Between 1987 and 2002, the IS area remained practically constant, even declining by 130 ha. It was between 2002 and 2017 that the IS area saw the largest increase of 3727 ha, equivalent to a total increase of 3597 ha in 30 years. The GS area has decreased by 2260 ha in 30 years. For BSC, between 1987 and 2002, the area decreased by 4575 ha, whilst between 2001 and 2017 this area increased by 3057 ha. In the 30 years analyzed, the BSC area decreased by 1519 ha.

The LST maps of the study area for 1987, 2002 and 2017 are shown in Figure 4. The average temperature for July 1987 was 29.92 °C, rising to 31.86 °C in 2002 and 39.02 °C in 2017. It can be observed that from 1987 to 2017 the land surface temperature has increased by 9.1 °C, the most notable increase being 7.16 °C in the time interval between 2002 and 2017; this coincides with the period in which the urbanization process intensified, whereas the increase for the interval between 1987 and 2002 was 1.94 °C. In general, the highest values were found in the city center, in the areas destined for the characteristic cultivation and baresoil. Figure 4 shows that urban areas are colder than the surroundings, suggesting the role of baresoil and cultivated land in the reversal of the SUHI phenomenon.

Figure 5 shows mean LST of the LULC classes for 1987, 2002 and 2017. In the three years analyzed, the urban class (impervious surface) together with baresoil and cultivated land had the highest mean LST. The temperature difference between both (impervious surface and baresoil and cultivated land) was 1.13 °C for 1987, 1.22 °C for 2002 and 0.63 °C for 2017, being the highest mean LST in the three years studied for baresoil and cultivated land than for the urban class.

As for green space and water classes, there were no changes in LST during the first study period (1987 to 2002). Nevertheless, mean LST increased considerably in the period of 2002–2017; 6.65 °C for the green space class and 9.37 °C for the water class. Therefore, the results indicate that urban areas along with baresoil and cultivated land have influenced the spatial pattern of LST.

3.2. Impervious Surface, Green Space and Baresoil and Cultivated vs. LST along Urban-Rural Gradient

Regarding the urban-rural gradient, mean LST remains practically constant (Figure 6). The results indicate that there was a downturn of the SUHI singularity (i.e., surface urban cool island). The agricultural land used for the cultivation of flowers, cereals and fruits does not fit with the normal spatial-temporal patterns of cultivation, but rather the requirements of the market. The differing cultivation on agricultural area and the varying phenology of vegetation in the area of urban green spaces causes seasonal variations in the spatial layout of vegetation that influence the spatial distribution of LST [96,97]. This is why there is no similar pattern of behavior between fraction of IS and mean LST alongside the rural urban gradient.

Across the three time periods analyzed, fraction of IS decreased while fraction of GS increased and then remained practically constant along the rural urban gradient. As for fraction of BSC, this followed a similar evolution to the fraction of GS with the difference in that it increases as the distance to the city center gets higher. In the studied time, the linear correlation analysis presents a negative relationship of mean LST with fraction of IS caused by the shadows projected by buildings, and a positive correlation with the fraction of GS caused by the influence of barren and cultivated lands that inverts the LST behavior pattern. In all cases, the correlation between mean LST and fraction of IS and GS is weak (Figure 7). The sample number varies according to the year and the land use/cover. Thus, for the year 1987 the sample size was 131,824, 105,495 and 194,022 for GS, IS and BSC, respectively. For the year 2002 it was 188,195 for GS, 104,149 for IS and 143,138 for BSC. Finally, for the year 2017 the sample size was 106,710, 145,565 and 177,150 for GS, IS and BSC, respectively.

On the other hand, the linear regression analysis between mean LST and baresoil and cultivated land is shown in Figure 8, where a positive relationship between them is observed with a stronger correlation with mean LST.

3.3. Impervious Surface, Green Space and Baresoil and Cultivated vs. LST at Multiple Resolution

Across grid sizes in Seville City, the correlation is stronger between mean LST, green space and baresoil and cultivated density than between mean LST and impervious surface density over the three years. From this three-year analysis, the strongest correlation between mean LST and green space density occurred in 1987, while 2002 was the year with the strongest correlation between mean LST and baresoil and cultivated land. In the case of impervious surface density, this correlation occurred in 2002 and 2017, with both being very similar (Figure 9 and Figure 10).

3.4. Landscape Variables Influencing Surface Temperature Spatial Variations

Table 2. Result of the OLS MLR analysis (Dependent Variable: Mean LST; N = 9216).

1987	OLS MLR Analysis
Variables	Unstandardized Coefficient		Standardized Coefficients (β)	Sig.	VIF
	β	Std. Error
(Constant)	33.005	1.048		0.000
Fraction of GS	−0.065	0.001	−0.762	0.000	1.423
Fraction of IS	−0.016	0.001	−0.176	0.000	1.257
Fraction of BSC	0.042	0.001	0.523	0.000	1.323
Mean Elevation	0.007	0.001	0.066	0.000	1.247
Mean Slope	0.040	0.012	0.028	0.001	1.288
Mean Hillshade	−0.008	0.006	−0.011	0.001	1.459
Mean Aspect	0.001	0.001	0.018	0.000	1.377
R2 = 0.545; Adjusted R2 = 0.545
2002	OLS MLR Analysis
Variables	Unstandardized Coefficient		Standardized Coefficients (β)	Sig.	VIF
	β	Std. Error
(Constant)	34.802	1.108		0.000
Fraction of GS	−0.052	0.001	−0.644	0.000	1.372
Fraction of IS	−0.007	0.001	−0.074	0.000	1.269
Fraction of BSC	0.032	0.001	0.342	0.000	1.240
Mean Elevation	0.011	0.001	0.109	0.000	1.233
Mean Slope	0.100	0.012	0.071	0.000	1.268
Mean Hillshade	−0.010	0.006	−0.015	0.001	1.457
Mean Aspect	0.003	0.001	0.047	0.000	1.370
R2 = 0.457; Adjusted R2 = 0.456
2017	OLS MLR Analysis
Variables	Unstandardized Coefficient		Standardized Coefficients (β)	Sig.	VIF
	β	Std. Error
(Constant)	45,289	1.008		0.000
Fraction of GS	−0.055	0.001	−0.642	0.000	1.284
Fraction of IS	−0.004	0.001	−0.046	0.000	1.236
Fraction of BSC	0.038	0.001	0.456	0.000	1.280
Mean Elevation	0.007	0.001	0.081	0.000	1.176
Mean Slope	0.016	0.011	0.013	0.001	1.257
Mean Hillshade	−0.035	0.006	−0.058	0.000	1.454
Mean Aspect	0.007	0.001	0.109	0.000	1.371
R2 = 0.429; Adjusted R2 = 0.429

reviews the results of the OLS MLR analysis. In respect to the variance inflation factor (VIF) values for all variables used in this study at the three periods, there was a low multicollinearity between explanatory or independent variables. Additionally, the results of the OLS MLR analysis present that, together, the explanatory variables considered in the analysis were significant in explaining a relevant number of spatial variations in mean LST at the three time periods (R² = 0.55 (1987); R² = 0.46 (2002); R² = 0.43 (2017)) (ρ < 0.001). The individual regression coefficients β of the explanatory variables were also statistically significant.

As for the standardized regression coefficients, the results present that a fraction of GS has a high β through the three time periods, while mean hillshade and slope had the lowest β for 1987, 2002 and 2017. A decreasing trend is also observed in the standardized β “Fraction of IS” in the three time periods (Table 2). This is consistent with the influence of baresoil and cultivated land on mean LST, which shows an increasing trend.

4. Discussion and Conclusions

Seville is the political and economic capital of the autonomous region of Andalusia which has played a very important role in its economic growth. In 2017, Seville generated a wealth of 39,500 million euros, 4% more than the previous year and above the Andalusian average (3%). The gross domestic product (GDP) generated by Seville was 24.50% of the Andalusian total, estimated at 161,111 million euros. In the province of Seville, 336 companies were founded in 2017, representing 25.03% of all the companies founded in Andalusia and almost 4.26% of those created in Spain [73]. Its population has also experienced considerable growth, rising from 1.4 million in 1987 to 1.9 million in 2017.

As shown in Figure 3, Seville has experienced rapid urban growth, as indicated by the substantial increase in the area classified as impervious surface (urban class). In 2017, the urban class increased by 37.87% with respect to the year 1987. In the first period analyzed (1987–2002), Seville underwent various changes such as the conversion into the administrative and political capital of the new autonomous region of Andalusia, and the celebration of the 1992 Universal Exposition (EXPO 92). During this period the construction of communication infrastructures, bridges, the elimination of railway barriers, as well as the urbanization of La Cartuja, was completed. There was a major urban expansion, especially towards the north (Pino Montano) and southeast (East Seville and Bermejales) of the city. Urban regeneration was also undertaken in historical sectors, which gave rise to a clear gentrification, causing a population migration and certain economic activities to move towards more suburban areas [98]. The growth observed in the second period analyzed (2002–2017) was to a certain extent caused by the urban development of Seville in the previous period. The area of impervious surface increased by 39.80% in this period of time.

In relation to green spaces, these are mainly distributed between the northern and southern zones of the Guadalquivir, corresponding to the fertile lowlands and marshes, respectively (Figure 3). In 1987 the percentage of green spaces in the analyzed area was 30.55%, while in 2017 this percentage decreased by 24.85%.

In the province of Seville, 864 km² of the total area is farmland; 134 km² of which is made up of meadows and pastures and 265 km² is forest land. In the area studied, the percentage of cultivated land in 1987 was 44.97%. In 2002 this percentage fell to 32.88% as a consequence of the urban development of the area and increased to 41.25% in 2017 at the expense of green space area.

As for LST, the highest values were found in the city center as well as in baresoil and cultivated land (Figure 4). This is due to the fact that surface characteristics, such as constructions, streets and other impervious surfaces, absorb more heat (solar radiation) than surfaces where the cover is vegetal (Figure 5) [33,77,79,99,100]. The fact that baresoil may have an LST similar to or even higher than that of the urban class (constructed or impervious surfaces) has been reported by other researchers [101,102]. In bare and dry soil, cities may have more vegetation than their environs, reversing the most common urban scenario in which relative lack of breathable vegetation and available moisture, as well as soil “waterproofing,” are among the reasons of UHI formation. The phenomenon of a “cool island,” which can result from the aforementioned, has been reported by several authors as being a predominantly daily phenomenon [79,103,104,105,106]. In this study, Landsat data capture three July months analyzed for different years. The satellite images were obtained between 10:27 and 11:02 AM, which may justify the fact that the city center had a lower LST than the nearby zones. Similarly, it has also been found that the characteristics of dry soil [107] and bare soil [108] produce high thermal values which, if part of the non-urban environment, might help an urban heat dissipation effect [109].

The date and time of acquisition of the satellite images would justify the existence of this heat dissipating effect found in the study area.

On one hand, the satellite data were obtained in the morning, which increases the possibility of detecting a weakened or dissipated urban heat island. The reversal (an urban heat sink) that happened in the day in the zone under study depends on certain surface situations. For any superficial material, specific internal properties, such as heat capacity, thermal conductivity and inertia, have a big importance in the control of the temperature of a body in balance with its environment [110]. These properties change according to the type of soil and its moisture content [111]. Dry, bare and low-density soils have been associated with high LST as a result of fairly low thermal inertia [112,113]. Soil emissivity depends on soil moisture situations and soil density [113]. Thus, for areas characterized by partial vegetation cover, the thermal surface properties can have a great influence in the measurement of LST through the thermal courses of conduction, convection and radiation.

In addition, satellite images were captured in early July, revealing additional surface features that would help in the creation of an urban heat sink. The crops were largely pre-emergent, and the satellite images used resulted in baresoil. As mentioned above, this baresoil state can strongly contribute to an increase in surface temperatures [108].

Therefore, there seems to have occurred a series of temporal and superficial characteristics favorable to the development of the heat sink within the study area and that would tend, consequently, to increase the temperature of the surface of the predominant agricultural land in the non-urban environment.

As can be seen in Figure 8, for the three years analyzed there is a positive relationship between the baresoil and cultivated land and the mean LST, similar to that obtained by other authors for impervious surface and mean LST [33,77]. The mixture of urban expansion, crop rotation, poorly managed cropland and vegetation degradation might have caused increases in bare/semi-baresoil, especially in the latter. The highest temperature found on desolate and cultivated land in the city contributes to the generation of a cool, urban island surface. Dry, barren soil has a low heat conductivity capacity and warms up quickly at dawn, while urban areas save solar energy. This phenomenon justifies the fact that the LST remains practically constant along the rural urban gradient and that the typical correlations obtained by other authors between IS and GS with mean LST are reversed [89,114]. This is demonstrated in Figure 6 where there is no decrease in LST as the distance to the city center gets higher. Although a temperature difference is observed between the three years analyzed, it remains practically constant along the urban gradient and even increases after 15 km which coincides with the beginning of the baresoil and cultivated land category [115]. The temperature of the urban class (impervious surface) follows a pattern consistent with that obtained by other authors [77]. It decreases as the distance to the city center increases, i.e., as a fraction of IS decreases, mean LST decreases. However, this is not followed by a decrease in mean LST as a consequence of the temperature reached by baresoil and cultivated land. This would justify the negative correlation of the fraction of IS and mean LST. The same happens with green space, although in Figure 6 there is an increase in the GS along the urban gradient, which translates into a decrease in temperature. If we observe Figure 7, we see that the correlation of fraction of GS with mean LST is positive due to the increase in baresoil and cultivated land.

From Figure 10, it is detected that impervious surface had the greatest influence on the average LST for the year 2002 with a slope value of 0.0188, followed by the year 2017 (slope = 0.0177) and the year 1987 (slope = 0.0143). It is also observed that green space had the greatest impact on mean LST for 1987 (slope = −0.0608), followed by 2017 (slope = −0.0546) and 2002 (slope = −0.0529). This confirms the urban growth experienced in the period 2002–2017 as a consequence of the changes implemented in the period 1987–2002.

The results show that in the case of fraction of IS, the correlations with mean LST are noticeably positive and increase as time goes by, taking a value from 0.0256 in the year 1987 to 0.0527 for 2017, which is consistent with the increase in baresoil and cultivated land. Conversely, the fraction of GS maintains a negative correlation with mean LST which decreases over time, from 0.5132 in 1987 to 0.4114 in the year 2017. This reinforces the fundamental role played by green space and baresoil and cultivated land in the formation of the UHI.

Cities are adopting an urban disperse model in that they tend to occupy increasingly larger areas with the removal of certain sectors outside the city limits (office parks, industrial activities, low-density residences, university institutions, etc.) for the creation of dormitories, etc. This is partly motivated by the emergence of a series of factors such as the increase in land prices, changing perceptions on the quality of life that influence the construction of new housing (the building of closed residential complexes with private gardens and the high value of being in close contact with nature, among others), the dominance of the car over the city, and so on. All of this means that cities need an ever-increasing consumption of energy and materials; therefore, making them less sustainable.

The solution to this problem is closely linked to efficient urban planning, in which measures that are based on an exhaustive territorial, economic and sociological analysis are adopted, aimed at restoring the cities’ environmental quality and reducing the effects generated by climate change.

In summary, it is vital to introduce an assessment planning culture in the framework of climate change [115,116,117,118]. The control of urban expansion, the increase of green areas (including roofs and building façades) as well as the percentage of permeable soil, the modification of the albedo of materials and pavements (increasing the degree of reflection of incoming solar radiation), the integration of artificial water bodies, the promotion of urban ventilation, the layout of buildings and, in general, the composition of urban morphology in order to facilitate air circulation, generates urban canyons and eases temperatures [119]. These are all elements that must be included in the daily practices of urban and land planning [120].

It concludes that there is a need to implement UHI mitigation strategies during the design and initial phases of the engineering project, from where the origin of this problem can be acted upon, since the process of creating streets and public space offers a valuable opportunity to restore the environmental quality of our cities and to diminish the effects generated by climate change.