Land Surface Temperature Dynamics in Response to Changes in Land Cover in An-Najaf Province, Iraq

Abstract

Land surface temperature (LST) is a critical environmental indicator affected by land cover (LC) changes. Currently, the most convenient and fastest way to retrieve LST is to use remote sensing images due to their continuous monitoring of the Earth's surface. The work intended to investigate land cover change and temperature response inAn-Najaf province. Landsat multispectral imageries acquired inAugust 1989, 2004, and 2021 were employed to estimate land cover change and LST responses. The findings exhibited an increase in water bodies, built-up areas, plantations, and croplands by 7.78%, 7.27%, 6.98%, 3.24%, and 7.78%, respectively, while bare soil decreased by 25.27% for the period (1989-2021). This indicates a transition from barren lands to different land cover types. The contribution index (CI) was employed to depict how changes in land cover categories altered mean region surface temperatures. The highest LSTs recorded were in bare lands (42.2℃, 44.25℃, and 46.9℃), followed by built-up zones (41.6℃, 43.96℃, and 44.89℃), cropland (30.9℃, 32.96℃, and 34.76℃), plantations (35.4℃, 36.97℃, and 38.92℃), and water bodies (27.3℃, 29.35℃, and 29.68℃) respectively, in 1989, 2004, and 2021. Consequently, these changes resulted in significant variances in LST between different LC types.

1. Introduction

The term land coverrefersto the physicalstate ofthe surface, which is often classified into several categories: forests, farmlands, arid lands, water bodies, and built-up areas (Khwarahm, 2021). The changes in land use/land cover (LULC) are important environmental indicators that reveal human-environment interactions which are detrimental to the climate system. Because of its role in the interchange of earth’s surface energy, surface matter, and physical and chemical processes with the atmosphere (Imran et al., 2021), thus LST is one ofthe environmental criteria affected by land cover changes (Deng et al., 2018). Among the most notable changes are the differences in thermal characteristics of the radioactive surface. More solar energy can be stored and transformed into sensible heat due to the thermal characteristics of soil, built-up land, and impervioussurfaces, also removing the trees and shrubs decreasesthe natural cooling influences of evaporation and shading (Doomi et al., 2016).

Remote sensing (RS) monitoring of LST offers a broad range of information and better spatial consistency when compared with the traditional observation manners employed at meteorological centers (AL–Ruwashdi and AL–Khakani, 2022). The spatiotemporal LULCchange with physical parameters such as surface radiation and emissivity is typically comprehended utilizing RS in conjunction with the geographical information system (GIS) method (Patil et al., 2018). As a result, this technique is a lot used in the field of thermal environmental studies (Khafid et al., 2020).

Many researchers used remote sensing techniquesto examine the impact of land cover changes on LST in different parts ofthe world. For example, Faqe Ibrahim (2017) evaluated and analyzed the impact of land cover changes on LST in Duhok city between 1990 and 2016. The study findings demonstrated that rising land surface temperatures are significantly influenced by changes in land use and cover. Where built-up zones and arid land typically experience the highest temperatures. Obiefuna et al. (2018) examined the impact of land cover alteration on LST in Lagos city, Nigeria between 1984 and 2015. Their findings demonstrated that the growing urbanization of Lagos city changed the ground thermal environment as evidenced by the LST increase. Also, Ogunjobi et al. (2018) investigated the effects of urban expansion on temporal changes of LST in Sokoto Metropolis, Nigeria from 1986 to 2016. The study revealed a positive link between the built-up region and LST.This evolution may be the result of human activitiesthrough urban growth along with its potential impacts on the urban climate.

Although there is a considerable lot of literature on the effects ofland cover change on land surface thermal characteristics in many places in the world, there are very few of them in An-Najaf Province which has seen an increase in surface temperature due to urban sprawl and irregular human activities.Thisstudy aims to: (1) assess the changes in land cover for the periods (1989–2004) and (2004–2021), (2) estimate the change in land surface temperature, and (3) evaluate the influence of changes in LULC on LST values for the period mentioned above. The findings of this study will interest urban planners, environmentalists, and decision-makers in developing ways to mitigate uncontrolled urbanization and inconvenience thermal.

2. Materials and Methods

2.1. Study Area

This paper investigates the changes in LST and LULC in An-Najaf province, which is positioned in the Middle Euphrates zone of Iraq, 160 km southwest of the capital, Baghdad, with an overall area of around 28,824 km² . It extends from 42°50′–44°44′ E (longitudes) and 29°50′–32°15′ N (latitudes). The average altitude of the city is about 70 m above sea level. A portion of An-Najaf province was selected as a study zone with an area of 1376.75 km² (Fig. 1) within the longitudes of 44°8′–44°34′ E and the latitudes of 31°40′–32°19′ N.The climate ofthe studied region ranges from semi-humid to semi-dry.

Fig. 1. Study site, from left to right: An-Najaf Province site in Iraq, study position from An-Najaf province, and RGB image (band combination 3, 5, 7) for study site from Landsat-8 OLI (2021).

The annual mean temperature reaches 32°C, the average temperature for August is 44°C, and its monthly means are decreasing to 14°C in January. The total annual precipitation ranges from (36 mm) to (190 mm), which often falls during the autumn, winter, and spring, while it interrupts in the summer.

2.2. Data Source

This study employed the geometrically rectified Landsat imagesfrom the official website of the United States Geological Survey (USGS) for the path/row: 168/036 ofAn- Najaf province.The image scenes were dated 9August 1989 and 2August 2004 from Landsat TM 5, and 1 August 2021 from Landsat 8 TIRS/OLI. The images are captured in the dry season to ensure that it isfree from clouds,so no atmospheric correction is required. For facilitating the image handling in the GIS environment, all data were geo-corrected and projected to UniversalTransverse Mercator, Zone 38N, based on the World Geodetic System (WGS) of 1984. Due to its widespread availability, Landsat imagery is commonly utilized for land cover and LST research. All multispectral and thermal bands have a medium spatial resolution, which is necessary for accurate land cover classification and LST extraction from thermal bands (Nse et al., 2020).

2.3. Classification of Land Cover

The supervised classification was performed using the maximum likelihood classifier (MLC) by the ArcGIS 10.8 program to classify the LULC types in the study zone. MLC performs better than other specific classification algorithms where it considers variability within category distributions (Guechi et al., 2021), and it utilizes a statistical method for pattern identification (Nse et al., 2020).

Using field survey data as well as prior knowledge of the area, the supervised classification was created. The region has been classified into five categories namely, cropland, bare soil, plantation, built-up, and water bodies. About 50 to 80 training samples were allocated for each of the five land cover classes, depending on the area of each category. MLC creates the probability function using data collected from training locations. The technique then compares each pixelfromthe image to known pixels(training locations) and classifies unknown pixels to one of the land cover pixels from the image to known pixels (training locations) and classifies unknown pixels to one of the land cover categories based on likeness and the highest likelihood of belonging to a previously known class (Jensen, 2005). The classified images for 1989, 2004, and 2021 were then used to determine the rate and amount of LULC patterns change across the study period.

2.4. Accuracy Assessment

The classification accuracy assessment is an essential stage in determining the validity of the results. Any image classification process’s objective is to get the maximum level of accuracy feasible (Sharma et al., 2011).Themost commonmethod to assess classification accuracy is to construct a collection of randomly categorized points of the classified image and then compare them to the ground truth data in the confusion matrix. A confusion matrix considers a common approach for validation, it is a two-dimensional table used to implement accuracy measures. Producer accuracy (PA), user accuracy (UA), overall accuracy (OA), and kappa coefficient (KC) are accuracy assessment criteria derived from the confusion matrix (Thakkar et al., 2017). For each land cover category, the accuracy of both the producer and the user are computed individually, based on the row allocation and matrix column as in Eq. (1) and (2) (Patel and Kaushal, 2010):

\(\begin{aligned}P A=\frac{\text { Number of correct samples in a type }}{\text { Sum of samples in the row }} \times 100\end{aligned}\)       (1)

\(\begin{aligned}U A=\frac{\text { Number of correct samples in a type }}{\text { Sum of samples in the colum }} \times 100\end{aligned}\)       (2)

The percentage of correctly identified samples ofthe total number of samples over the map is referred to as overall accuracy, and is calculated as in the Eq. (3) (Sreelekha and Reddy, 2019):

\(\begin{aligned}O A=\frac{\text { Total Number of correc samples }}{\text { Total Number of samples }} \times 100\end{aligned}\)       (3)

The Kappa coefficient uses to evaluate the image classification accuracy (Cohen, 1960). It is a measure ofthe agreement that depends on the variance between the true agreement in the confusion matrix and the agreement of chance. The Kappa coefficient value ranges between 0 and 1, and it can be calculated as follows (Thakkar et al., 2017):

\(\begin{aligned}K=\frac{\text { The total sum of correct }- \text { Sum of all the }(\text { row total } \times \text { column total })}{\text { Total squard }- \text { Sum of all the }(\text { row total } \times \text { column total })}\end{aligned}\)       (4)

2.5. Land Surface Temperature (LST) Retrieval

This study adopted the thermal bands of Landsat 5 (band 6) and Landsat 8 (band 10) to retrieve the LST ofAn-Najaf province. Wang et al.(2015)found that the Landsat 8 data product for band 11 of the TIRS sensor has great uncertainty. Bendib et al.(2017) proposed that Landsat 8’s TIRS band 10 is employed as one spectral band for land surface temperature estimation. For this reason, band 10 of Landsat 8 was chosen to retrieve the LST. Four stages have been applied to calculate LST, the first stage differs in Landsat 5 than in Landsat 8 images, while the other three steps are the same for the two Landsat images. The LST can be retrieved using metadata files from the Landsat series. The steps for retrieving LST are as follows:

2.5.1. Digital Number Conversion to Spectral Radiance (L_λ)

For Landsat 5,spectralradiance was calculated from the digital values of the thermal band 6, applying the following formula (Zareie et al., 2016):

\(\begin{aligned}L_{\lambda}=\left(\frac{L_{\text {MAX }}-L_{\text {MIN }}}{Q_{\text {CALMAX }}-Q_{\text {CALMIN }}}\right) \times\left(Q_{\text {CAL }}-Q_{\text {CALMIN }}\right)+L_{\text {MIN }}\end{aligned}\)       (5)

Where L_λ representsthe spectral radiance measured by (Watts/m²·sr·μm), Q_CALMIN represents the minimum value of Digital Number (DN) (1), Q_CALMAX represents the maximum value of Digital Number (DN) (255), Q_CAL is the digital number of the band 6, L_MIN is the minimum radiance (1.238), L_MAX is the maximum radiance (15.600). While formula (6) has been utilized to calculate the spectralradiance forthe Landsat 8 taken from the USGS website:

L_λ = M_L × Q_CAL + A_L       (6)

M_L is the radiance multiplicative scaling factor of the band (10 or 11) and, A_L is the band-specific additive rescaling factor. L_MIN, L_MAX, Q_CALMIN, Q_CALMAX, M_L, and A_L values were obtained fromthe attachedmetadata file with Landsat images.

2.5.2. Spectral Radiance (L_λ) Conversion to AT-Satellite Brightness Temperature (T_B)

After calculating the spectralradiation L_λ, the second step was to convert L_λ to (T_B), both Landsat 5 and Landsat 8 images utilized the Eq. (7) (Zareie et al., 2016):

\(\begin{aligned}T_{B}=\frac{K_{2}}{\ln \left(\frac{K_{1}}{L_{\lambda}}+1\right)}\end{aligned}\)       (7)

Where T_B indicates at-satellite brightness temperature in (°C), K₁ and K₂ are the calibration constants for thermal bands for Landsat 5 and Landsat 8, which can also be obtained from the metadata file (Table 1).

Table 1. Calibration constant of thermal bands

2.5.3. Land Surface Emissivity Estimation (ε)

The surface emissivity of a body is defined as its ability to convert heat energy into quantifiable radiant energy. The ability of the ground surface to turn heat energy into radiating energy is critical for the retrieval of LST. The emissivity (ε) was computed according to the proposed equation by Sobrino et al. (2004) as shown below:

ε =0.004 × P_v + 0.986       (8)

P_v representsthe vegetation proportion estimated based on Weng et al., (2004), as in Eq. (9):

\(\begin{aligned}P_{v}=\left(\frac{N D V I-N D V I_{\min }}{N D V I_{\max }-N D V I_{\min }}\right)^{2}\end{aligned}\)       (9)

Where: NDVI is the normalized difference vegetation index, NDVI_min is the minimum value of NDVI and NDVI_max isthemaximumvalue ofNDVI.The following relation was used to express the NDVI:

\(\begin{aligned}N D V I=\frac{\rho(\text { Band5 })-\rho(\text { Band } 4)}{\rho(\text { Band })+\rho(\text { Band } 4)}\end{aligned}\)       (10)

Where ρ (band4) and ρ (band5) represent the spectral reflectance of bands 4 and 5, respectively, of Landsat 8 images. The associated metadata file with Landsat images containsthe parameters employed to transform the DN to values of spectral reflectance.

2.5.4. LST Calculation

In the last step, the brightness temperature was converted to the LST in degrees Celsius as indicated in Eq. (11) (Mujabar, 2019):

\(\begin{aligned}L S T=\frac{T_{B}}{1+\left(\frac{\lambda \times T_{B}}{\rho}\right) \times \ln \varepsilon}-273.15\end{aligned}\)       (11)

Where λ = emitted radiance’s wavelength (λ = 11.5 μm), ρ = hc / σ (1.438×10^–2 m·K), σ = Bolzmann’s constant (1.38×10^–23 J/K), h = Planck’s constant (6.26×10^–34 J·s), c = light’s velocity (2.998×10^–8 m/s).

2.6. Land Cover Relation to LST

The impact ofland cover on the cooling or warming of any region relies on the kind of land cover and the proportion of the total area occupied by each category. For instance, vegetation and water/humid areas have a cooling effect on the surface because of latent heat transfer. However, although there is a cooling effect, the total value depends on the ratio of the total area it occupies (Anderson et al., 1976). The contribution index (CI) measures the extent to which a land cover type is warming or cooling, taking into account the ratio of the total area it covers. The CI utilizes to associate long-term land cover changes with intensities of LST. For each land cover type, the CI was calculated for all the periods studied using Eq. (12) (Odindi et al., 2017).

CI = Dt × S       (12)

Where CI denotes the relative LST contribution of the entire region, Dt is the difference between the cover type’s average temperature and the entire region’s average temperature, and S is the land cover type’s proportional area.

3. Results and Discussion

3.1. Classification Accuracy Evaluation

The maps ofland use land coverfor 1989, 2004, and 2021 were created using the maximum likelihood classifier, where the study area was divided into five distinct categories. Due to classification techniques and image acquisition methods, land cover maps typically have some errors (Ogunjobi et al., 2018).

Arandom sample from 250 pixels was used to verify the accuracy of classified images, then compared to remote sensing images obtained by Google Earth. The confusion matrix with relevant statistics (overall accuracy and kappa coefficient) was applied to assess the classification accuracy findings. Table 2 provides the accurate evaluation results for the classified images forthree years.Overall classification accuraciesfor 1989, 2004, and 2021 were 93.33%, 92.12%, and 93.84%, respectively. These outcomes are consistent with the required minimum accuracy (85%) recommended by Anderson (Anderson, 1971).

Table 2. Classified images accuracy results for 1989, 2004, and 2021

For the three years (1989, 2004, and 2021), kappa coefficients were 0.91, 0.90, and 0.92, respectively. Since the land covers obtained by Google Earth are fully compatible with the land covers identified at different random locations, the accuracy of the classified categories was high.According to Congalton and Green (2019), the assessment results were highly accepted.

3.2. Changes Analysis for Land Cover

The final categorized maps were created by viewing various land cover classes in the studied region for 1989, 2004, and 2021, as in Fig. 2.

Fig. 2. Land cover maps of the study region for (A) 1989, (B) 2004, and (C) 2021.

Table 3 displays the areal distribution of different categories and their percentages from 1989 to 2021, and Fig. 3 shows this graphically.

Table 3. Areal distribution and percentage of land cover in 1989, 2004, and 2021

Fig. 3. The percentages scheme of land cover categories for the three years.

Throughout the research period, bare soil was the prevailing kind ofland cover.As a result ofthe drought of Bahr An-Najaf, where vast portions of it converted into sandy areas, and the decline of agricultural zones, the exposed soil increased from 729.41 km² (52.98%) to 745.92 km² (54.18%) between 1989 and 2004. During the same period, plantations decreased from 267.36 km² (19.42%) to 249.6 km² (18.13%) due to poor land management policy. Also, water bodies decreased from 43.79 km² (3.18%) to 19.16 km² (1.39%). While built-up areas increased from 40.35 km² (2.93%) to 74.06 km² (5.38%). Finally, croplands slightly declined from 295.84 km² (21.49 %)to 288.01 km² (20.92%).

Between 2004 and 2021, there was a considerable decrease in a barren land, which fell to 398.07 km² (28.91%) of the total area, due to the conversion of a large portion of it to residential areas for suiting the expanding population’s housing needs. Also, due to large parts of plantations being transformed into residential areas, the plantation area hasfallen to 345.73 km² (25.11%), increasing the built-up area to 174.1 km² (12.65%) in 2021. On the other hand, the croplands (most of which represent rice fields) increased to 332.65 km² (24.16%) to meet the needs of the local population for the rice crop, which is the primary food people. Finally, the water bodies have risen significantly to 126.2 km² (9.17%), represented by the Bahr An-Najaf depression, which is one of the characteristic water features in the region, due to increases in rainfall and thus increased torrential water, which is the primary source of Bahr An-Najaf water.

3.3. Spatial Distribution of LST between 1989 and 2021

Due to the different reflections of various land covers, LST differsfor different land cover categories. The findings demonstrate thatforthe years 1989, 2004, and 2021, the LSTofthe examined area varied between around (26.62–48.96°C) (average 38.46°C), (27.22–51.29°C) (average 40.61°C), and (28.24–53.36°C) (average 42.59°C), respectively.

According to the spatial distribution of LSTthrough the studied periods, the maximum LST was at the bare soil zone, whereasthe minimum LST was at the water bodies areas(Fig. 4). The LST at the bare soil area was higherthan that ofthe built-up region, this might be due to the study area’s climatic condition where the bare soil zone of An-Najaf province represents a part of the desert area. This result is consistent with several studiesthatshowed exposed soils had the highest LST level. Faqe Ibrahim (2017) and Ogunjobi et al. (2018), mentioned that LST in various land covers might vary depending on the data acquisition time since built-up surfaces acquire heat more slowly than bare soil, whereas, bare soil radiates the heat more rapidly than the built-up surfaces. Also, it noted that as the built-up area extended, there was a corresponding increase in the LST rates around the urban area, which could be due to radiating the heat created by the built-up region, as confirmed by the priors studies (Orimoloye et al., 2018;Babalola et al., 2016;Rotem-Mindali et al., 2015).

Fig. 4. LST Maps of An-Najaf city in August for the three years, (A) 1989, (B) 2004, and (C) 2021.

Bare soil areas showed the highest rates of LSTs reached (42.20°C in 1989, 44.25°C in 2004, and 46.90°C in 2021), followed by built-up areas(41.60°C in 1989, 43.96°C in 2004, and 44.89°C in 2021). The rates of LSTfor plantation land were (35.40°Cin 1989, 36.97°C in 2004, and 38.92°C in 2021). Following water bodies, croplands exhibited the lower LSTs of (30.90°C in 1989, 32.96°C in 2004, and 34.76°C in 2021), which represents rice crop fields that are often submerged in water, limiting the amount of heat stored in soil and surface structures via evaporation and transpiration, this reduces the temperature of this type of land cover. While, the water bodies showed the lowest rates of LST values reached (27.30°C in 1985, 29.35°C in 2004, and 29.68°C in 2021). Fig. 5 illustratesthemean LSTof variousland coversin 1989, 2004, and 2021.

Fig. 5. The differences in mean LST with various land coverings.

3.4. The Relation between Land Cover Changes and LST Dynamics

The land cover change map demonstrated that in locations where the land cover classes transformed into buildings or bare land, the surface temperature had increased at a significant rate. Where the most significantrise in the surface temperature was observed when water bodies converted into the barren area (19.6°C),followed by croplandsinto bare soil(14.54°C), croplands into the built-up area (12.29°C), and water bodies into croplands (5.52°C), whereas a decrease of LSTnoticed when bare soil transformed into both water bodies (–15.63°C), croplands (–12.06°C), plantation (–7.72°C) and built-up areas (–2.28°C), as well as converted plantations into croplands (–5.21°C) (as illustrated in Table 4).

Table 4. The difference in mean LST with various land covering​​​​​​​

The contribution index for each land cover category was calculated for the three years as in Table 5. The contribution of water bodies in cooling, as evidenced by CI, ranged from about -0.28 in 1989, -0.11 in 2004, to -1.16 in 2021 as a result of the alteration of water bodies areas during the study years. Dense vegetation areas have been characterized by high negative CIs (1.44,-1.61,-1.87for1989,2004, and2021,respectively), indicating that plant density reducesthe warming ofthe region during the dry season, with the highest negative contributions in 2021, as a result of the expansion of rice crop fields. There was also a negative CI, but they were lower than in dense vegetation (-0.20, -0.66, to -0.92 for 1989, 2004, and 2021, respectively) in the low vegetation zones that reduce the cooling effect of these areas. In general, bare soil areas exhibited higher positive CIs; however, these values decreased over for three years, from 1.98 to 1.80 to 1.25, respectively, as a result of a decline in bare soilregions and the presence of scattered farms, which lessen the warming effect.

Table 5. Thermal contribution for each land cover category based on a contribution index (CI)​​​​​​​

Due to the expansion of built-up zones during the period time, the built-up zone’s positive warming impact grew from 0.09 to 0.17 to 0.30 in 1989, 2004, and 2021, respectively. In general, the cooling impact in the two areas of dense/low vegetation coverincreased over the study period, besides the decreased warming effect in the bare lands, because of the growth and decline in the two zones, respectively.

4. Conclusions

Multispectral remote sensing data was employed to determine land cover types and the changes in LST for the period between 1989 and 2021. Based on the results, it was concluded that multispectral Landsat images and the maximum likelihood approach accurately extract spatial patterns of land cover. The currentstudy revealed the status ofland cover dynamics and their potential repercussions on the thermal characteristics ofthe land surface in the studied region. The land cover classification result illustrated that both waterbodies, cropland, plantation, and built-up areas experienced an increase while, bare land decreased, throughout the study, as a result of economic, social, and political factors.

The results illustrated that the value of LST differs by the variousland cover classesforinstance, bare land and built-up areas increased radiant temperature. The study indicated that LST in the barren land was higher than for the built-up areas, unlike many previous studies, because of the semi-arid environment of the region. The increasing water level and Bahr An-Najaf area led to a marked decline in LSTin the western parts ofthe study area. LSThas generally declined along the eastern and southeastern strip of the region due to the increasing cultivation ofrice fieldsto fill the local need for this crop, which is the primary food source of the population.As a result ofthe mentioned facts, we could recommend to urban designersthat they improve green space both within and around the urban centers to reduce the rates of LST.