Evaluation of the effects of land cover changes and urbanization on land surface temperature A remot

2.1. Study area

The study area is located in Mascara city, mainly, in sub-watershed of Oued Fekan(fig.1.b), located in north-western of Algeria between latitude 35° 7'N to 35° 31' N and longitude 0° 25' W to 0° 4.143' E, thereby covering a total geographical area of approximately 1195 km². It includes the Ghriss plain which covers an agricultural area of more than 72000 hectares. It is limited from the north by the Beni-Chougrane Mountains, from the south by the Saida Mountains, from the west by Bouhanifia Mountains and in the east, we find the trays of Tighennif.

Figure 1. Location map of study area (a) Sub-watersheds of the Macta, (b) sub-watershed of Oued Fekan. 2.1.1. Topography The plain corresponds to a flat stretch with an average altitude of 450 m, overhanging by high edge ridges that exceed the 1000 m of altitude in the south (Nesmoth Mountain). Non-permanent watercourses drain this topographic catchment area which is upstream in the mountains of BeniChougrane and Nesmoth, while its downstream is located in the western part. The western part is the outlet which feeds the Bouhanifia dam.

Figure 2. Digital Elevation Model of Oued Fekan. 2.1.2. Climatic and rainfall context This region is characterized by a semi-arid climate where the annual average of precipitations is around 313 mm/year for the period of 1976-2001 [26]. For the same period, the actual evapotranspiration was estimated at 289 mm/year while the average annual temperature is 16.5 °C. The accumulated rainfall between December and February corresponds on average to 37% of the total annual rainfall. The spring and autumn rains correspond respectively to 30% and 28% of this total, while the summer rains, which are stormy and torrential, represent only 5% of the total. The area has no significant spatial rainfall variability although the landforms are more watered than the central plain. Very high temperatures are measured during the period from June to August with maximum values exceeding 40 °C during the day. It favours the evaporation of the waters of the superficial aquifer. The wadis of the plain have non-perennial flows with very low or null annual average flows.

Figure 3. Evolution of rainfall in Mascara, Matemoure, Tizi, Ghriss and Froha meteorological station. 911

2.1.3. Soil representation The study of pedology is essential to understand the nature and types of the soils that exist in Oued Fekan region. This part has been permanently taken from works made at the ANRH (National Agency of Hydraulic Resources) level where the plain has been the subject of many soil surveys. Lands which outcropping the plain are from the sedimentary formation with variable texture, composed predominately of recent and older alluvium. Generally, soils are of calcimagnesic type, but occasionally we meet isohumic and poorly evolved soils [18]. All these types of soil in the study area are favourable to the infiltration of precipitation water for recharging plio-quaternary groundwater and represent the nourishing ground, naturally fertile. Its significant agricultural potential characterizes the plain since more than 90% of its total area is currently used for agricultural activity.

2.2. Data sources

Our study area is fully included in the Landsat 8, Path/Row 197/36. As shown on the table 1, the images used in this analysis were acquired in 2018, respectively in February, May, August and November so that it covers four different seasons. Images of Landsat 8 OLI-TIRS (Operational Land Imager Thermal Infrared Sensor) are downloaded from the Earth Explorer US Geological Survey image database (Landsat Imagery Archive). Table 1. Specification of the used satellites images. Satellite

OLI (Operational Land Imager) Landsat 8and TIRS 197/36 (Thermal Infrared Sensor)

Band 1: 0.435-0.451 Band 2: 0.452–0.512 Band 3: 0.533–0.590 Band 4: 0.636–0.673 Band 5: 0.851–0.879 February 30 Band 6: 1.566–1.651 May (For TIRS bands resolution 2018 10:32 Band 7: 2.107–2.294 August is 100 m) Band 8: 0.50–0.68 November (For band 8 resolution is 15 Band 9: 1.363–1.384 m) Band 10: 10.60–11.19 (Thermal band) Band 11: 11.50–12.51 (Thermal band)

2.2.1. Preprocessing of OLI Data The USGS EROS Center provides standard Landsat 8 data products. It consists of quantized and calibrated scaled Digital Numbers (DN) representing multispectral image data acquired by both the Thermal Infrared Sensor (TIRS) and Operational Land Imager (OLI). The products are delivered in 16-bits; we rescaled them to the Top of Atmosphere (TOA) radiance using radiometric rescaling coefficients provided in the product metadata file (MTL file). 𝐿λ = 𝑀𝐿 𝑄𝑐𝑎𝑙 + 𝐴𝐿

Where Lλ is TOA spectral radiance (Watts/ (m2 * srad * μm)), ML is Band-specific multiplicative rescaling factor from the metadata (RADIANCE_MULT_BAND_x, where x is the band number), AL is Band-specific additive rescaling factor from the metadata

(RADIANCE_ADD_BAND_x, where x is the band number), Qcal is Quantized and calibrated standard product pixel values (DN). The spectral radiance in the optical range is converted into reflectance after correction for atmospheric effects using commercially-available MODTRAN Software radiative transfer code [27]. These reflectances are subsequently exploited to calculate the normalized difference vegetation index (NDVI) and the normalized difference built-up index (NDBI).NDVI optimize contrasts between the infrared and the visible bands. This index makes it possible to visualize, on a single channel, the dynamic responses related to the density of vegetative cover. NDVI is included between -1 and +1.When the area corresponding to the soil has a strong chlorophyllin activity; it marks a high value [28]. NDVI is then calculated as the combination between near infrared and red band: (𝑅

NDBI based on the analysis of the unique spectral responses of built-up areas in seven Landsat TM bands. The original NDBI approach developed by Zha is implemented based on Landsat OLI bands covering NIR and SWIR wavelength [29], using the following equations: (𝑅

2.2.2. Pre-processing of TIRS Data The following relationship express the spectral space-reaching radiance (Lsat(λ)) measured by the sensor in the thermal infrared: 𝐿 𝑠𝑎𝑡(λ) = [𝜀0 𝐿λ (T0 ) + (1 − ε0 )L↓ atm(λ)] ∗ τλ + L↑ atm(λ)]

Where: Lλ(T0) is the radiance of a blackbody target of kinetic temperature T0. λ is atmospheric transmission. Latm(λ) is the down welling or sky radiance. Latm(λ) is the upwelling or atmospheric path radiance. Both atmospheric parameters (λ ,Latm(λ) , Latm(λ)) are estimated, by the web atmospheric correction parameters calculator, at time of satellite overpass [30]. 0 is the surface emissivity which is estimated from the vegetation index (NDVI) [31]: 𝜀0 = 1.0094 + 0.047 × ln(𝑁𝐷𝑉𝐼)

Spectral radiances Lλ(T0) reflected from the Earth surface can be deducted by inversion of Equation 4. After converting TIRS band data from DN to spectral radiance, we convert them to the top of atmosphere brightness temperature using the thermal constants provided in the metadata file: 𝑇=

Where: T Top of atmosphere brightness temperature (K) Lλ TOA spectral radiance (Watts/( m2 * srad * μm)) K1 Band-specific thermal conversion constant from the metadata (K1_CONSTANT_BAND_x, where x is the thermal band number 10 or 11) K2 Band-specific thermal conversion was constant from the metadata (K2_CONSTANT_BAND_x, where x is the thermal band number 10 or 11).

Table 2. Band-specific thermal conversion constant from the metadata K1 and K2 values. Thermal constant K1 K2

TIRS bands were designed to allow the calculating surface temperature. For our analysis and calculation of LST, we used only bands 10 from the Thermal Infrared Sensor (TIRS) of the Landsat 8 satellite. According to Wang, since TIRS band 11 data from the Landsat 8 have considerable uncertainty, he also recommended using only TIRS band 10 data as a single spectral band for any LST estimation [32]. Finally, the land surface temperature is calculated as follows: 𝐿𝑆𝑇 =

Where: Tb = At-satellite brightness temperature (in Kelvin) ε = emissivity λ =.wavelength of emitted radiance (11.5 μm) ρ = h × c/σ (1,438 × 10-2 m K) σ = Boltzmann constant (1,38 × 10-23 J K-1) h = Planck's constant (6,626 × 10-34Js) c = Velocity of light (2,998 × 108 m sec-1)

Figure 4. Flowchart of retrieving LST values from Landsat 8 image.

3.1. Spatiotemporal variation of LST

In order to study the behaviour of LST in our study area, the LST image was generated, as shown in figure 6. The area under study is characterized by recurrent drought, where the semi-arid climate is predominant. Two significant periods describe this area, dark and rainy period from November to March and another hot and dry, from April to October [18]. Winter (from December to February) is usually cold enough as the absolute minimum of the air temperature descends to -2 °C. Summer (from June to August) is usually hot and dry since the absolute maximum of the air temperature is equal to +50 °C. A hot drying wind from the South called "sirocco" accentuate the thermal maxima. The remote sensing approach allowed measuring the LST of the entire study area, the interpretation of the spatial distribution of LST with the type of land use becomes possible. As it can be seen from Figure 6, especially fig 6.a, boundaries in south-eastern part are the coldest; on the other hand, the northwest part is warmer. It represents the Bouhanifia Mountains, an area known for its wealth of thermal water. It is found that the warmer area was also marked in the centre of the plain. Particularly where areas were filed agricultural and now marked as built-up areas. Temperature map shows that LST marks important values mainly at the built-up area and urban fallow land, which is going to be built up. During the winter season, the minimum and maximum of land surface temperature are 6.93 and 33.74 °C, respectively; for the summer season, the most basic marked value of LST is 25.67 °C while the maximum value is 51.21 °C.

An average of 43.06 °C during the hot summer season against an average of 18.524 °C during the cold winter, which clearly explains why this study area is classified as semi-arid in its type of climate. The absence of a strong relationship between LST and UHI during the four seasons; especially in summer where bare land is predominant; returns to urban sprawl type exerted on the non-compact plain (except for major cities that are not metropolitan); otherwise it is discontinuous, scattered or diffuse and generally linear along the roads. This type of urban sprawl explains the homogeneity of LST marked in almost all the plain. Despite being an agricultural area, marked values are essential.

Figure 6. Spatial distribution of LST for four seasons in degrees Celsius (°C). Validation of the spatial distribution of surface temperature and its relation to the type of land occupation, including UHI effects, is made as it can be seen in figure 7. The resulting LST class was located on the imagery. It was observed that Oued Fergoug show very low LST, also parcels with a high density of vegetation. The urban areas of Mascara city show LST higher than the vegetated areas and water-body but lower than the bare soil. Urban areas of Mascara, Tighenif and Ghriss cities have marked a difference in the surface temperature higher than adjacent agricultural parcels of the plain about 2 to 3 °C in winter. This difference in LST changed along with the variation of season; these results involve studying the relationship between Normalized Difference Built Index and Land Surface Temperature and analyzing implications of urban growth on thermal characteristics.

Figure 7. Visual interpretation of LST distribution using Reference Maplink

3.2. Spatiotemporal variation of NDVI and NDBI

Spatial representation of NDVI and NDBI distribution over the four seasons are shown respectively in figure 8 and figure 9. Firstly, relationships between NDVI and NDBI was analysed; a weak and negative correlation was marked between the two variables. An inverse relationship between NDVI and NDBI was marked. Height values of NDVI indicate vegetated areas; while height values of NDBI indicate built-up areas and barren land. Both indices were used to differentiate the land use/cover types of our study area. The most significant correlation occurred in summer. It suggests that NDVI can respond to NDBI change and reflect the evolution of the urban construction land (Figure 10).

Figure 10. Scatter plots of NDVI (x-axes) and NDBI (y-axes) for four seasons 919

3.3. Relationships between the LST and NDVI

In order to study the relation of LST with the respective indices, scatter plots were obtained for LST/NDVI, LST/NDBI, as shown in Figure 11 – 12, respectively. Areas with high temperature generally tend to have lower NDVI values. This aspect has mainly corroborated from the analysis of the relationship between LST and NDVI. As shown in the LST images, the LST values of non-vegetated areas (dried river beds, barren land, built-up) are higher than those marked for areas of water bodies and vegetation, including agricultural land. The results obtained from the correlation between LST and NDVI showed non-linear behaviour between them. In urban regions, NDVI cannot explain LST due to more significant heterogeneity in land cover, especially in urban areas that include soil, buildings, vegetation and water. NDVI cannot well describe the LST of the study area since the trend of LST and NDVI relationship during seasons is different. Affected by the seasons, the relationship between NDVI and LST cannot achieve a stable level which indicates that the correlation between the two indices is influenced by season.

Figure 11. Scatter plots of LST (x-axes) and NDVI (y-axis) for four seasons As shown by Figure 12, higher LST values are usually found in the areas of higher NDBI values because of the predominance of high-density urban areas, including commercial, industrial, and residential developments. In contrast, analysis of the above data shows that LST correlates closely with NDBI. There is a strong positive correlation between the values of both indices, indicating a direct relationship between them. As NDBI values increase LST values also

increase; generally in the same areas. Instead NDVI, NDBI can be used to study LST in the region containing different land covers.

Figure 12. Scatter plots of LST (x-axes) and NDBI (y-axis) for four seasons.

3.5. A further investigation between LST, NDBI and NDVI

A deeper interpretation between the mean of LST, NDBI and NDVI was applied over all the heterogeneous area that encompasses the three cities of Mascara, Maoussa and Tighennif. Thus a large area of agricultural land showed a very strong positive relationship between the mean of LST and NDBI for all seasons (r2> 0.86). The supported method to set the pace of NDBI/NDVI is explained in the following relation: PaceNDBI= (MAXNDBI-MINNDBI)/n PaceNDVI= (MAXNDVI-MINNDVI)/n Where: n = number of points retained; (Mostly we fixed n at 50 points). The obtained results are presented in Fig. 13 and Fig. 14: Contrariwise, the relationship between the mean of LST and NDVI was not accurate, in all graphs (Figure 13), we note that there is a stabilisation indicating a strong negative correlation between NDVI and LST(r2 of NDVI vs. LST ranges from 0.001 to 0.049). These inverse linear correlations begin when NDVI>0.2 for spring and summer. This same negative correlation starts

instead of when NDVI reaches the value of 0.18 for February and November. Results obtained show a negative correlation which varies by season; for values lower than the one mentioned before (NDVI<0.2 for May and August, NDVI<0.18 for February and November) including negative values which represent mainly non-vegetated area, urban surfaces and bare soil, no stable trend was identified. NDVI has an evident seasonal change which implies using this index is complicated to understand UHI.

Conversely, as shown in figure 14, the co-linear correlation between NDBI and the mean of LST (r2 of NDBI vs. LST ranges from 0.861 to 0.925) was marked. LST is almost co-linear with air temperature as indicated by significant correlation value (0.861 for February, 0.915 for May, 0.925 for August and 0.91 for November 2018) at 0.02 level of significance. Its less variation with the season (not subject to seasonal effects) makes it an apparent index for future analysis of spatiotemporal UHI effects. The relationship between NDBI and LST was more stable comparing to NDVI and LST. NDBI can adequately describe the LST since they were positively correlated in four seasons (without seasonal effects.). The full index can be used as a supplemental indicator to the traditionally applied NDVI for further analysing seasonal changes in UHI effects.

3.6. Spatial dependence and correlation among NDBI, NDVI and LST from two profiles

To study a deeper Spatio-temporal distribution of LST, NDBI and NDVI, which is profiles of each index are drawn in two different areas. The first, located in Mascara city, covers buildings, roads, settlement and vegetation. While the second, located in the western part, covers water of Bouhanifia dam, sparse settlement and grassland that around it. Horizontal and vertical profiles have been made across the study area to represent the LULC wise LST (Figure 15). Profiles lines are shown in high-resolution Google Earth images acquired in 05/04/2018.

Figure 15. Selected LST, NDBI and NDVI Profiles. By analysing obtained maps and profiles, it is identified that built-up areas like buildings, roads are warmer than vegetation-covered areas and grassland, but water bodies remain coolest of all. The urbanization process in sub-watershed of Oued Fekan had a significant impact on its thermal characteristics since the mean temperatures were higher in urban areas than in the surrounding suburban/rural areas. Heavily industrialized districts and built-up land, densely populated, are closely associated with high-temperature anomalies. Except in wetlands and waters where the NDVI profile shows a better correlation, both profiles of LST and NDBI varies in parallel with a strong correlation. All three indices (LST, NDBI and NDVI) used in our study, can

describe, in a qualitative manner, the spatial distribution and temporal variation in the plain thermal environment. Extracting LST from thermal remote sensing data is useful for studying other subjects such as evapotranspiration. Results obtained showed that industrialized districts and built-up land, densely populated, are closely associated with high-temperature anomalies. We recommend controlling the urbanization exercised in Ghriss plain, avoiding compacted urban sprawl, with thoughtful planning of the distribution of industry, which will mitigate the UHI effect for a healthy urban environment and sustainable development.