Analyzing the correlation of forest and wetland with land surface temperature by using geospatial technology: a case of Yayo district, Southwestern Ethiopia

Abstract

Forests and wetlands played a crucial role in minimizing the increasing trends of land surface temperature (LST). This study aimed to analyze the contribution of forests and wetlands on LST in Yayo district using multi-spectral and multi-temporal satellite. Landsat images of 1986, 2003 and 2021 were downloaded from United State Geological Survey website with 30 m spatial resolution. Results show that the mean LST was increased from 22.8 °C in 1986 to 27 °C in 2003 over cultivated land. As the results of forest conservation and wetland protection, the LST in Yayo district is experencing a declining trend in the recent year. Result reveal that the normalized difference vegetation index and Normalized Difference Water Index show a strong negative relationship with LST. Over the last two decades, the LST was declined by 1.5 °C due to the protection of forest and wetland ecosystem. Thus, governmental and non-governmental organization should encourage the conservation and protection of forest and wetland ecosystems.

Keywords:

• Correlation
• conservation
• uSGS
• normalized difference vegetation index
• normalized deference water index

Introduction

The role of forests and wetland in adapting climate change is very high. Studies clearly documented the role of forest and wetland ecosystems in climate change mitigation and biodiversity conservation worldwide (Erwin 2009; Junk et al. 2013; Seddon et al. 2020; Huang et al. 2021; Yan et al. 2022). Both forest and wetland have huge potential in carbon sequestration, balancing ecosystem, restoring ground water, purifying toxins, stabilize microclimates and minimize the potential occurrence of climate extremes such as drought and floods (Gebre and Gebremedhin 2019; He et al. 2021; Kattel 2022; Zhou et al. 2022). However, deforestation and forest degradation were the main cause for climate change (Van Khuc et al. 2018; Moisa, Dejene and Gemeda 2022; Yang et al. 2023). Previous studies confirmed that there is an inverse relationship between land surface temperature (LST) and forest and wetland ecosystems (Singh et al. 2017; Muro et al. 2018; Kayumba et al. 2021; Li, Wang, et al. 2022; Pei et al. 2023). To understand the role of vegetation cover and wetland ecosystem Merga et al. (2022) analyzed the relationship between normalized difference vegetation index (NDVI), and normalized difference water index (NDWI) with LST and confirm a negative correlation in Didessa River sub-basin of Western Ethiopia. In addition, assessment of land use land cover (LULC) impact on LST is very decisive for monitoring spatial distribution of surface temperature (Moisa, Dejene, Hinkosa, et al. 2022; Moisa et al. 2022b).

Several studies have been conducted on the impacts of LULC change on LST in different parts of Ethiopia (Moisa, Dejene and Gemeda 2022; Moisa, Dejene, Hinkosa, et al. 2022). According to Tafesse and Suryabhagavan (2019), LST in Adama has been increased due to uncontrolled LULC changes which accelerate loss of vegetation covers. Balew and Korme (2020), recognize that LULC change brought higher accumulation of LST in and around Bahir Dar city. Merga et al. (2022) also reported that over expansion of agriculture and illegal settlement causes dramatically decline of vegetation which increase LST in Didessa River Sub-basin in Western Ethiopia. Although substantial studies have been conducted in different parts of the country, no studies reported the role of safeguarding forests and wetlands in mitigating LST, they were dealing with LULC changes, attitudes on some economical plants like coffee Arabica and woody plants (Moat et al. 2017; Getahun and Keno 2019; Oostdijk 2021). Forest conservation is one of the key techniques to minimize the increasing trends of LST. The southwestern part of Ethiopia is relatively better in terms of forest cover as compared to other parts of Ethiopia. Study by Tadese et al. (2021) reported that four biosphere reserves, including Yayo, Majang, Sheka, and Kafa were established between 2010 and 2017 in southwestern part of Ethiopia. Forest and wetland have been a crucial role in minimizing LST, but little research attention has been given until to date in the study area. Geospatial technology was used to analyze the relationship betwwn forest, wetland and LST to offer evidence based information for decision makers and natural resource manager for sustainable use of forest and wetland to mitigate climate change. Therefore, the present study aimed to fill the existing research gap by analyzing the correlation between forest and wetland with LST using geospatial technology in Yayo district, southwestern Ethiopia.

Materials and methods

Description of the study area

This study was conducted in Yayo district which located in the Ilu Ababor Zone, of Oromia National Regional State, SW Ethiopia (Figure 1). Geographically, the Yayo district lies between 8⁰7′00″ and 8⁰21′00″ N and 35⁰46′30″ and 36⁰7″E and the topography of the study area varies between 1,226 and 2,534.5 above sea level. It covers an area of 807.7 km². The mean annual temperature is about 23.74°C, and the district is known by its rainy condition with precipitation of about 1,625 mm/year (Tadese et al. 2021). Yayo is located within the Eastern Afromontane Biodiversity Hotspot, and montane rainforest fragments with wild Coffee arabica. Three management zones namely; core, buffer and transition zones were established next to its designation as a United Nations Educational, Scientific and Cultural Organization (UNESCO) biosphere reserve. Yayo Biosphere Reserve has a surface are of 167,021 ha, of which 27,733 ha core area, 21,552 ha buffer zone, and 117,736 ha transition zone. The core zone is free from human intervention and is fully protected; the buffer zone is an important ecological corridors, connecting the core zone to the transitional zones, while the transition zone is free for wise use of natural resources (Tadese et al. 2021). The transitional zone is found adjacent to the buffer zone and includes agricultural lands, wetlands, grasslands, settlements, and forest fragments. These collectively made a conducive battle for conservation of the forest in present study area (Abera et al. 2021).

Figure 1. Location map of the study area.

Data sources and descriptions

Landsat images of Thematic Mapper (TM) in 1986, Enhanced Thematic Mapper Plus (ETM+) of 2003 and Operational Land Imager (OLI) of 2021 were downloaded from United States Geological Survey (USGS) website (http://eathexplorer.usgs.gov/) with 30 m spatial resolution by using path 170 and row 54 which are cloud free images during the dry season (January and February) (Table 1). Earth Resources Data Analysis System (ERDAS) Imagine 2015 software was used for LULC classification. ArcGIS 10.3 was performed for analyzing and visualizing spatial data. Landsat images of 1986, 2003 and 2021 were used in the present study due to free from strip and cloud cover. From the downloaded Landsat images, multispectral bands were used for classification of LULC and calculation of NDVI and NDWI from 1986 to 2021 in the Yayo district.

Table 1. Data sources and its descriptions.

Data types	Sensor	Acquisition year	Path and row	Resolution (m)	Source
Landsat 5	TM	January 1, 1986	170 & 054	30	(http://eathexplorer.usgs.gov/)
Landsat 7	ETM+	January 1, 2003	170 & 054	30	(http://eathexplorer.usgs.gov/)
Landsat 8	OLI/TIRS	February 1, 2021	170 & 054	30	(http://eathexplorer.usgs.gov/)
ASTER DEM	–	–	–	30	(http://eathexplorer.usgs.gov/)

Methods of data analysis

Analysis of LULC classes

The LULC in Yayo district was classified from Landsat images of TM of 1986, ETM + of 2003 and OLI/TIRS of 2021. The LULC of the study area were classified into five classes: i.e. bare land, cultivated land, forest land, grassland and wetland by using supervised classification method with maximum likelihood algorithm.

Accuracy assessment

For accuracy assessment of LULC classes, 100 ground control points were collected using GPS, while Google Earth was used inaccessible areas in the study area. All LULC classes on the map were validated with reality on the ground using the collected ground points. 100 ground control points collected by randomly sampled 20 ground control points for forest land, cultivated land, bare land, grassland and wetland. The overall accuracy was calculated by dividing all the pixels properly classified by total number of pixels in the matrix Equation (1)(1) $$\mathit{Overall}\mathrm{ }\mathit{accuracy}=\frac{⅀X_{ij}}{N}\mathrm{*}100$$(1) . The user accuracy is calculated by multiplying the total number of correctly classified pixels in a class by the sum of the values in the rows of the same class Equation (2)(2) $$\mathit{User}\mathrm{ }\mathit{accuracy}=X_{ii}/X_j\mathrm{ }\mathrm{*}100$$(2) . The producer accuracy index is calculated by dividing the total of the values in the same class’s column by the number of correctly identified pixels in that class as shown in Equation (3)(3) $$\mathit{Producer}\mathrm{ }\mathit{accuracy}=X_{ii}/X_k\mathrm{*}100$$(3) . (1) $$\mathit{Overall}\mathrm{ }\mathit{accuracy}=\frac{⅀X_{ij}}{N}\mathrm{*}100$$(1) where N is total number of accuracy sites, $X_{ij}$ is sum of diagonal values. (2) $$\mathit{User}\mathrm{ }\mathit{accuracy}=X_{ii}/X_j\mathrm{ }\mathrm{*}100$$(2) where $X_{ii}$ is total number of corrected pixels and $X_j$ is column total. (3) $$\mathit{Producer}\mathrm{ }\mathit{accuracy}=X_{ii}/X_k\mathrm{*}100$$(3) where X_ii is the total number of corrected pixels and X_k is row total.

Land use land cover change detection

Rate of LULC change

The extent of forest cover change was calculated following recent studies (Moisa and Gemeda 2021; Moisa, Dejene, Hinkosa, et al. 2022) by using Equation (4)(4) $$\mathit{Rate}\mathrm{ }of\mathrm{ }\mathit{change}\mathrm{ }\left(\frac{{km}^2}{\mathit{year}}\right)=\frac{{\mathit{X}}^2-{\mathit{X}}^1}{\mathit{Z}}$$(4) . (4) $$\mathit{Rate}\mathrm{ }of\mathrm{ }\mathit{change}\mathrm{ }\left(\frac{{km}^2}{\mathit{year}}\right)=\frac{{\mathit{X}}^2-{\mathit{X}}^1}{\mathit{Z}}$$(4) where, X² is an area of LULC in square kilometer in time 2, X¹ is an area of LULC in square kilometer in time 1; Z is time interval between X2 and X1 in years.

Trends of land use land cover change

Trends of LULC used to show negative numbers in the trend indicate a decrease (losses), whereas positive values indicate an increase (gain) in that LULC as indicated Equation (5(5) $$\mathit{Trend}\mathrm{ }of\mathrm{ }\mathit{LULC}=Y2-Y1$$(5) ) (Moisa and Gemeda 2021; Wang, Wang, et al. 2022). (5) $$\mathit{Trend}\mathrm{ }of\mathrm{ }\mathit{LULC}=Y2-Y1$$(5) where: Y2 is the current year, Y1 is the previous year

Normalized difference vegetation index estimation

The NDVI used to calculate the quantity of above ground green vegetation cover (Moisa, Dejene and Gemeda 2022; Moisa, Dejene, Hinkosa, et al. 2022) as indicated in Equation (6)(6) $$\mathit{NDVI}=\frac{\mathit{NIR}-R}{\mathit{NIR}+R}$$(6) . (6) $$\mathit{NDVI}=\frac{\mathit{NIR}-R}{\mathit{NIR}+R}$$(6) where NDVI is Normalized Difference Vegetation Index,

NIR is the near infrared band and R is the red band

Landsat 8 data used Bands 5 (Infrared) and 4 (red), whereas Landsat 5 and 7 used Bands 4 and 3, to calculate the NDVI values.

Normalized difference water index (NDWI)

The NDWI is used to simulate water areas, which often have significant thermal changes (Zhou and Wang 2011; Siqi and Yuhong 2020). Green band (band 2 for Landsat 5 and 7, band 3 for Landsat 8) and near infrared (band 4 for Landsat 5 and 7, band 5 for Landsat 8) reflectance measurements were used to build the formula Equation (7)(7) $$\mathit{NDWI}=\frac{\mathit{Green}-\mathit{NIR}}{\mathit{Green}+\mathit{NIR}}$$(7) . The Near-Infrared (NIR) and Green (G) channels are used to calculate the NDWI (Jackson et al. 2004; Wang, Liu, et al. 2022; Fang et al. 2023). (7) $$\mathit{NDWI}=\frac{\mathit{Green}-\mathit{NIR}}{\mathit{Green}+\mathit{NIR}}$$(7)

Retrieval of land surface temperature

The thermal band of Landsat TM 1986, Landsat ETM+ 2003, and Landsat OLI/TIRS 2021 were used to determine the LST in the study area (Kumari et al. 2018; Zhang et al. 2021; Wolteji et al. 2022). Thermal bands are band 6 for Landsat 5 TM 1986 and Landsat ETM+ 2003, and band 10 for Landsat OLI/TIRS 2021. The Mono-window approach was used to calculate LST (Chibuike et al. 2018; Li, Charles, et al. 2022).

Step 1: Conversion of digital number (DN) into radiance

The digital numbers of TM and ETM + were converted to radiance using Equation (8)(8) $$\text{L}\lambda =\frac{\mathit{LMAX}\lambda -\mathit{LMIN}\lambda }{\mathit{QCALMAX}-\mathit{QCALMIN}}\mathrm{*}\left(DN-\mathit{QCALMIN}\right)+\mathit{LMIN}\lambda$$(8) (Chander et al. 2009). ETM + has a digital number between 0 and 255. (8) $$\text{L}\lambda =\frac{\mathit{LMAX}\lambda -\mathit{LMIN}\lambda }{\mathit{QCALMAX}-\mathit{QCALMIN}}\mathrm{*}\left(DN-\mathit{QCALMIN}\right)+\mathit{LMIN}\lambda$$(8) where;

QCAL is the quantized calibrated pixel value in DN

LMINλ is the spectral radiance that is scaled to QCALMIN in Watts/(meter squared*ster*µm

MAXλ is the spectral radiance that is scaled to QCALMAX in Watts/(meter squared*ster*µm

QCALMIN is the minimum quantized calibrated pixel value (corresponding to LMINλ) in DN

QCALMAX is the maximum quantized calibrated pixel value (corresponding to LMAXλ) in DN = 255

The digital number of band 10 from Landsat 8 TIRS was converted into radiance values using Equation (9)(9) $$L_{\lambda }=\left(M_L\mathrm{*}{Q}_{\mathit{Cal}}\right)+A_L$$(9) . (9) $$L_{\lambda }=\left(M_L\mathrm{*}{Q}_{\mathit{Cal}}\right)+A_L$$(9) where;

L_λ is the Top of Atmosphere (TOA) spectral radiance Watts/(meter squared*ster*µm)

M_L is Band specific multiplicative rescaling factor from the meta data (Radiance multi band x)

A_L is Band specific additive rescaling factor from the meta data (Radiance Add Band x, where x is the band number)

Qcal is quantized and calibrated standard product pixel values of DN

Step 2: Conversion to brightness temperature

LST was calculated using brightness temperature based on land surface emissivity (Yang et al. 2017; Xu et al. 2022). TM and ETM + values from band 6 were transformed from spectral radiance to brightness temperature (Yue et al. 2007) Equation (10)(10) $$\mathrm{TB}=\frac{K2}{ln\left(\frac{K1}{L\lambda }+1\right)}$$(10) . (10) $$\mathrm{TB}=\frac{K2}{ln\left(\frac{K1}{L\lambda }+1\right)}$$(10) where; TB is effective at satellite temperature in Kelvin, K₂ is Calibration constant 2.

K₁ is Calibration constant 1, Lλ is Spectral radiance in Wm² ster^-1 µm^-1

Step 3: Estimation of land surface emissivity using NDVI

The NDVI has been used to calculate the abundance of vegetation cover on the earth’s surface (Elmore et al. 2000; Li et al. 2021; Moisa, Tiye, et al. 2022). The NDVI value was derived using multispectral bands from Landsat images taken in 1986, 2003 and 2021. Bands 4 and 3 of Landsat 5 and 7 were utilized to measure near infrared (NIR) and red, respectively. Furthermore, in Landsat 8, bands 5 and 4 were employed to assess NIR and red, respectively.

A proportional vegetation ($P_V)$ calculation was done based on the NDVI measurements to calculate land surface emissivity (Neinavaz et al. 2020; Zhang et al. 2021). The vegetation proportion was calculated (Carlson and Ripley, 1997) using Equation 11(11) $${P_V\mathrm{ }=\mathrm{ }\left[\frac{\mathrm{NDVI}\text{-}\mathrm{ }\mathrm{NDVImin}}{\mathrm{NDVImax}\text{-}\mathrm{NDVImin}}\right]}^2$$(11) . (11) $${P_V\mathrm{ }=\mathrm{ }\left[\frac{\mathrm{NDVI}\text{-}\mathrm{ }\mathrm{NDVImin}}{\mathrm{NDVImax}\text{-}\mathrm{NDVImin}}\right]}^2$$(11) where; PV is proportion vegetation, NDVI is Normalized difference vegetation index,

NDVImin is minimum value of NDVI, NDVImax is maximum value of NDVI

An estimation of LST was performed (Chibuike et al. 2018; Liu et al. 2022) and calculated using Equation (12)(12) $$\mathrm{Ƹ}=0.004\mathrm{*}{P}_V+0.986$$(12) . (12) $$\mathrm{Ƹ}=0.004\mathrm{*}{P}_V+0.986$$(12)

Where; $\text{Ƹ is emissivity and }{P}_V\text{ is vegetation proportion}$

The calculated land surface temperature was corrected for emissivity (Srivastava et al. 2010; Zhu et al. 2023) using Equation (13)(13) $$\mathrm{LST}\mathrm{ }=\mathrm{ }\frac{\mathrm{TB}}{1\mathrm{ }\left(\mathrm{\lambda }\frac{\mathrm{TB}}{\mathrm{P}}\right)\text{InƸ}}$$(13) . (13) $$\mathrm{LST}\mathrm{ }=\mathrm{ }\frac{\mathrm{TB}}{1\mathrm{ }\left(\mathrm{\lambda }\frac{\mathrm{TB}}{\mathrm{P}}\right)\text{InƸ}}$$(13) where;

LST is land surface temperature (in Kelvin), TB is radiant surface temperature (in Kelvin)

$\text{λ is the wavelegnth of emitted radiance }(10.8\mathrm{\mu }\mathrm{m}),$ $\text{P is h*c}/\mathrm{\sigma }(1.438\mathrm{*}{10}^{-1}\mathrm{ }\mathrm{mK}),$ $\text{h is Planc}{\mathrm{k}}^{\prime }\text{s constant }(6.26\mathrm{*}{10}^{10-34}\text{ J s}),$ $\text{c is the velocity of light }\left(2.998\mathrm{*}{10}^8\frac{\mathrm{m}}{\mathrm{s}}\right),$ $\text{σ is Stefan Boltzman’n}\text{s constant }(1.38\mathrm{*}{10}^{-23\mathrm{ }}\mathrm{J}\mathrm{ }{\mathrm{K}}^{-1}),$ $\text{Ƹ is Land surface emissivity}.$

To calculate LST in degree Celsius, 273.15 should be subtracted from LST in kelvin (Equation (14)(14) $$°\mathrm{C}=K-273.15$$(14) ). (14) $$°\mathrm{C}=K-273.15$$(14)

Pearson correlation coefficient

Pearson correlation coefficient (Pearson’s r) is the most widely used correlation coefficient and measures of the statistical relationship between two quantitative variables (Liu et al. 2022; Šverko et al. 2022; Zhou et al. 2022; Burgund et al. 2023). Pearson’s r values range from −1 to +1. The Pearson’s r between 0 and 3 is weak, between 3 and 5 is moderate and greater than 5 is strong. A negative Pearson’s r value indicates a negative relationship, and a positive value indicates a positive relationship between the two variables.

Results and discussions

Accuracy assessment

The accuracy of the classified LULC maps was assessed to ensure their trustworthiness. The reference data was compared to the classified LULC types. The overall accuracy results of LULC classification are 86.8, 83.4 and 85.7% for the year 1986, 2003, and 2021. As a result, for the study periods 1991, 2003, and 2021, the kappa coefficients were 0.85, 0.82, and 0.83, respectively. In addition, the result of user accuracy and producer accuracy of 1986, 2003 and 2021 were stated in (Table 2).

Table 2. Accuracy assessment of LULC for 1986, 2003 and 2021.

LU/LC Types	1986		2003		2021
	Producers accuracy (%)	Users accuracy (%)	Producers accuracy (%)	Users accuracy (%)	Producers accuracy (%)	Users accuracy (%)
Cultivated	89	82.6	87.5	95	88.6	92.8
Grassland	90.1	88.2	92.3	91.8	93.3	95.2
Forest land	92.3	93.5	94.3	95.8	96.4	94.7
Bare land	84	86	92	84.9	92.2	90.9
Wetland	82.6	81.7	90.5	85	91.7	92
Overall Accuracy	86.8%		83.4%		85.7%
Kappa coefficient	0.85		0.82		0.83

Analysis of LULC change

As agricultural communities, there are rapid conversions of LULC for different economic activities. As compared to other region, the southwestern parts of Ethiopia are covered by dense and open forest. For instance, Yayo is one of the few remaining area of intact Afromontane forests in southwestern part of Ethiopia. As the demand for agricultural expansion caused by rapid population growth, the forest in this particular area is under threats. To minimize the threat of forest, UNESCO gives high attention for this particular area and established the Yayo Biosphere Reserve in 2010. Southern part and northern parts of the study area was conserved in 2021 than 2003. As the result, wetland area and forest land were increased in final year of the study period. To examining the extent of change over different land cover classes, we classified into cultivated land, bare land, forest land, grass land and water body in all study period (Figure 2). The results revealed that the share of cultivated land in 1986 was highest, which covers about 165 km². In contrst, cultivated land was declined due to forest conservation and protection. Similarly, bare land decreased from 30.4 to 1.3 km². Moreover, forest cover shows substantial increment from 319.4 km² in 1986 to 488.6 km² in 2021 (Table 3). The increment of forest cover and wetlands is associated with the establishment of Yayo biosphere reserve (Tadese et al. 2021). The result of this study was consistence’s with the previous studies (Liu et al. 2020; Negash et al. 2021; Moisa, Dejene and Gemeda 2022; Moisa et al. 2022a).

Figure 2. LULC map of the study area.

Table 3. Calculated LULC of the study area.

LULC Types	1986		2003		2020
	Area (km^2)	Area (%)	Area (km^2)	Area (%)	Area (km^2)	Area (%)
Bare land	30.4	3.8	47.9	5.9	1.3	0.2
Cultivated land	165.0	20.4	170.0	21.0	53.2	6.6
Forest land	319.4	39.5	412.8	51.1	488.6	60.5
Grassland	152.7	18.9	98.6	12.2	115.5	14.3
Wetland	140.2	17.4	78.4	9.7	149.0	18.4
Total	807.7	100.0	807.7	100.0	807.7	100.0

The results indicate that, forest land; grassland and wetland land cover classes experienced a positive trend with an increment by 169.2, 8.8 km² between 1986 and 2021. This increment implies the existence of management practices for restoration and rehabilitation of degraded natural forests, grassland and wetland ecosystem. In contrast, cultivated land, grassland, and bare experienced a negative trend by an area of 111.7, 37.2, and 29.1 km², respectively. The declining of cultivated land, grassland and bare land indicates the presence conservation actions in the study area (Figure 3). Similar findings have been reported in southwestern parts of Ethiopia as the results of conservation measures (Tadese et al. 2021; Moisa, Dejene, Hinkosa, et al. 2022).

Figure 3. LULC trends.

Results show that cultivated land and bare land were decreased by the rate of 3.3 km²/year and 0.9 km²/year from 1986 to 2021, respectively. However, forest land, wetland and grassland were increased by the rate of 2.1, 1.7 and 0.3 km²/year over the study period, respectively (Table 4). The rate by which the cultivated land and bare land decreased due to the recovery of forest land, grassland and wetland in the study area. This decreasing rate of cultivated land doesn’t mean the total collapses of cultivation activities rather than applying improved and sustainable agricultural practices on a right place at time. Local community is planting trees on their land in the form of agro forestry for different purposes. These findings are consistent with previous study (Moisa and Gemeda 2021). Tadese et al. (2021) also reported the increasing trend of forest cover as the results of Reducing Emissions from Deforestation and Forest Degradation (REDD+) since 2010 in similar geographical region.

Table 4. Rate of LULC types.

LULC classes	1986–2003	2003–2021	1986–2021
	km^2/year	km^2/year	km^2/year
Bare land	−1.8	0.1	−0.9
Cultivated land	−1.3	−5.2	−3.3
Forest land	12.8	−8.6	2.1
Grassland	−4.6	5.3	0.3
Wetland	−5.1	8.5	1.7

Post classification and Land Use Transfer Matrix (LUTM) methods were used for LULC change detection. The results revealed that about 9.7 km² of cultivated land and 0.4 km² of bare land converted to forest land between 1986 and 2021. Similarly, about 46 km² of grass land area also shifted to forest land. However, about 489.5 km² (60.6%) of the study area was remained unchanged (Table 5; Figure 4). These findings are consistency with previous studies (Moisa and Gemeda 2021; Moisa et al. 2022b).

Figure 4. LULC matrix of the study area.

Table 5. LULC matrix from 1986 to 2021.

LULC classes		2021(km^2)
		Bare land	Cultivated land	Forest land	Grassland	Wetland	Total
1986 (km^2)	Bare land	0.4	12.2	0.4	13.8	3.6	30.4
	Cultivated land	0.8	34.2	9.7	93.2	27.1	165.0
	Forest land	0.0	0.3	309.6	4.8	4.7	319.4
	Grassland	0.1	5.6	46.0	42.1	59.1	152.8
	Wetland	0.0	0.9	25.7	10.3	103.2	140.1
	Total	1.3	53.2	391.3	164.2	197.7	807.7

Spatiotemporal distributions of NDVI and NDWI

The NDVI was used to evaluate the status of forest cover in the study area. The maximum value of the NDVI implies higher vegetation and vice versa. The result revealed that maximum NDVI values declined from 0.63 in 1986 to 0.54 in 2003 as conservation actions do not get attention by governments until the establishment of Yayo biosphere reserve. Declining of NDVI values are associated with agricultural expansion to forest areas to meet the high demand for agricultural production in the region. After the establishment of biosphere reserve and protection of the forest by the government in the study area, the NDVI value was increased to 0.55 (Figure 5). The result of this study was more agreement with the previous studies (Li et al. 2020; Merga et al. 2022; Moisa, Dejene and Gemeda 2022; Moisa, Merga, et al. 2022). In addition NDVI values, which can provide information on vegetation status of the study area, interviews were conducted with key informants to understand the status of forest cover after the establishment of the Yayo Biosphere Reserve. Thus, key informants interview minimizes uncertainty about the actual status of forest cover on the ground. This is because the NDVI values represent only a few extreme pixels, making it very difficult to draw conclusions about the status of forests in a particular region with a high degree of uncertainty. Results show that the central and southern parts of the study area has higher NDVI values whereas, northern parts has lower NDVI values from 1986 to 2021.

Figure 5. NDVI map of the study area.

The status of water bodies were determined by NDWI. The results of this study revealed that the maximum values of NDWI were declined from 0.13 to 0.05 in 1986 to 2003 before conservation initiatives. However, the values of NDWI increased by 0.36 after conservation measures and later establishment of Yayo biosphere reserve. This implies the higher deterioration of wetland over the past decades was reversed after the establishment of Yayo Biosphere Reserve. These findings are validated by key informant interviews. Ethiopia government recognized the importance of biosphere reserve in the protection of natural resources in southwestern parts of Ethiopia including Yayo biosphere reserve (Tadese et al. 2021). The southern parts of the study area occupied with lower NDWI value while the northern and northwestern parts are characterized with higher NDWI values (Figure 6).

Figure 6. NDWI map of the study area.

LST in response to LULC change

The year selection for this study was depended on the establishment of Yayo biosphere reserve by UNESCO which located in Yayo district. The mean LST directly inter related to status of LULC change in the study area. Results show that the LST of Yayo district was increased from 19 to 21 °C in 1986 to 2003. Before the Yayo Biosphere Reserve is established in 2010, the LST increased by 2 °C from 1986 to 2003. This could be occurred due conversion of forest and wetland to cultivated land which exhibit the status of carbon stock and rises the temperature in the study area. Every country is experiencing the effects of climate change, but depending adaptive capacity, the severity may vary. Temperatures in Ethiopia warmed-up in every region. Climate change has been seen in wettest parts particularly Southwestern Ethiopia (Gemeda et al. 2021; Wang et al. 2023). Climate change causes for declining of rainfall during the main crop growing season and have an impact on agricultural productions and results in food insecurity. However, between 2003 and 2021 the mean LST decreased from 21 to 17.5 °C due to forest cover enhancement and wetland conservation practices. This conservation measure contributes to the decline of LST by 1.5 °C from 1986 to 2021. The declining of LST over the study period might be due to establishment of Yayo Biosphere Reserve as one of forest conservation mechanisms (Tadese et al. 2021). Spatially, central part of the study area was experienced with higher LST while northern and southern parts were mostly encompassing with lower LST due to conservation of forest land and wetland (Figure 7). The area covered with high and low LST was validated by using ground control points which collected by Google earth and handled GPS. The results of this study are inconsistency with previous studies (Merga et al. 2022; Moisa et al. 2022b). Contrasting results are due to variations in geographical locations, including the extent of forest and wetland ecosystems. Merga et al. (2022) conducted around Didessa River Sub-basin, while the work of Moisa et al. (2022b) was around the Anger River Sub-basin. Both sites have less forests and wetlands than Yayo.

Figure 7. LST map of the study area.

Correlation between LULC and LST

The relationship between LULC and LST was analyzed over study period. The results show that the mean LST over cultivated land was increased from 22.8 to 27 °C between 1986 and 2003. Similarly, the mean LST over forest declined from 21.1 °C in 2003 to 19.6 °C after forest conservation strategies (Table 6). Between 1986 and 2003 the LST over bare land, cultivated land, forest, grassland and wetland were raised by 6.8, 4.2, 3.6, 5.7 and 5.8 °C, respectively. However, from 2003 to 2021, substantial decline of the mean LST has been observed in all LULC classes. The status of LST in bare land, cultivated land, forest, grassland and wetland were decreased by 3, 0.6, 1.5, 1.3 and 2.6 °C, respectively. The declining of LST from 2003 onwards indicates the presence of conservation practies. This result is inconsistency with previous studies (Moisa, Dejene and Gemeda 2022; Merga et al. 2022) who reported the increasing trend of LST over different LULC classes (Moisa, Dejene and Gemeda 2022; Merga et al. 2022).

Table 6. Relationship of LULC and LST.

LULC types	1986	2003	2021
	Mean LST (°C)	Mean LST (°C)	Mean LST (°C)
Bare land	22.9	29.7	26.7
Cultivated land	22.8	27.0	26.4
Forest land	17.5	21.1	19.6
Grassland	19.2	24.9	23.6
Wetland	17.1	22.9	20.3

The correlation between LST and NDVI

Results show that LST and NDVI have inverse relationship (Figure 8). Because an increase in NDVI value leads to a decrease in LST. The Pearson’s r values indicates that LST and NDVI has strong negative correlation (r= −0.99). The adjusted R-square value also indicates strong relation between the two variables (0.978). This implies that higher vegetation density related to lower LST (Figure 8). The decline of vegetation cover due to expansion of agricultural land was the main causes for increments LST (Moisa, Dejene and Gemeda 2022; Merga et al. 2022; Yang, Li, et al. 2022). In addition, the study more consistence with the previous studies (Yang et al. 2021; Zhao et al. 2021), LST was increased during vegetation cover was declined and vice versa.

Figure 8. Correlation between LST and NDVI.

The correlation between LST and NDWI

The results show a strong negative association between LST and NDWI. The Pearson’s r values indicate a strong negative correlation between LST and NDWI (r= −0.93). The adjusted R-Square value also shows a strong relation between the two variables (0.86). Degradation of NDWI causes increments of LST and vice versa (Figure 9). The results of this study are consistent with studies reporting an inverse relationship between LST and NDWI (Zhou and Wang 2011; Yang, Dou, et al. 2022). In addition, the result of the present study in line with which conducted in mountain in the southeastern Tibetan Plateau and indicated that, LST has strong influence on water resources and glacial changes (Zhao et al. 2021).

Figure 9. Correlation between LST and NDWI.

Regression analysis of LST

In this study, inferential statistics were used in the form of a one-way analysis of variance (ANOVA) to examine the effect of LULC on LST (Table 7). This regression analysis was performed between independent variable (NDWI and NDVI) and dependent variable (LST). Bala et al. (2018) and Mensah et al. (2020) reported a non-linear relationship between LST with NDVI and NDWI. The LST of the study area was decreased due to increasing values of NDVI and NDWI. Yang et al. (2021) stated that, LST has strong relationship with vegetation cover (NDVI) and water resources in the Tibetan Plateau.

Table 7. Regression analysis of LST with NDVI and NDWI.

Variable	Coefficients	Standard error	t Stat	p-Value	Lower 95%	Upper 95%
LST	53.44467	3.162513197	16.89943	6.22E-07**	45.966515	60.922825
NDWI	24.00133	5.047826436	4.754786	0.002072*	12.06512	35.937546
NDVI	−78.4	8.34543903	−9.39435	3.22E-05**	−98.1338	−58.66614

Dependent variable: LST. Coefficient is statistically significant at (*p < 0.01 and ** p<0.001).

Conclusion

Wetlands and forest played a key role in minimizing LST and climate change by absorbing and storing carbon to mitigate atmospheric greenhouse gases effects on ecological balance. In the study area, forest land was increased by an area of 319.4 km² (39.5%) in 1986 to 488.6 km2 (60.5%) in 2021. Wetland land cover is also increased by 8.8 km². The mean LST was decreased by 1.5 °C due to conservation initiatives around Yayo district over the past 35 year. The results of this study can be useful for policy makers to invest on forest and wetland conservation to minimize the problem of climate change in general and LST is particularly.

Consequently, the mean LST was increased from 1986 to 2003 by 22.8–27 °C in cultivated land. However, in 2021 the mean of LST decreased to 26.4 °C after conservation measure was taken for Yayo Forest. In addition, the LST was decreased from 2003 to 2021 with the mean of 21.1–19.6 °C after wise use of natural resources. The mean LST in cultivated land increased from 22.8 °C in 1986 to 27 °C in 2003, while 26.4 °C in 2021 due to forest and conservation in the study area. Similarly, mean LST was increased in 2003, but decreased in 2021 due initiation of forest and wetland conservation. The result was validated by using Google earth and GPS ground truth data the area captured with high LST and low LST. In addition, LULC classes were validated by collected ground truth data. To minimize LST, natural resource managers and environmental analyst should increase public awareness on the importance of forest and wetland conservation. Further studies which considers socio economic data about the impacts of wetland and forest degradation on the rising of LST should be conducted in the study area. Moreover, high resolution satellite images can be used to classify LULC classes at five year interval by minimizing uncertainties and soil dryness influence on LST estimation.

Ethical approval

The authors agreed to publish the manuscript on Scandinavian Journal Forest Research.

Author contributions

MBM involved in research design, data collection, data analysis, and draft manuscript writing. YNW, GGN, MDN and BTG involved in data analysis MMG, KTD and DOG participated in methodology, data analysis and manuscript edition. All authors read and approved the final manuscript.

Acknowledgments

The authors would like to acknowledge Wollega University Shambu Campus Faculty of Technology, Faculty of Agriculture, Faculty of Resource Management and Economics, Wollega University Nekemte Campus, Raya University and Jimma University College of Agriculture and Veterinary Medicine for the existing facilities to conduct this study. No funding received for this research.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Available in this manuscript.