Spatiotemporal variability of land surface temperature in north-western Ethiopia

The aggravating deforestation, industrialization, and urbanization are becoming the principal causes for environmental challenges worldwide. As a result, satellite-based remote sensing helps to explore the environmental challenges spatially and temporally. This investigation analyzed the spatiotemporal variability in land surface temperature (LST) and its link with elevation in the Amhara region, Ethiopia. The Moderate Resolution Imaging Spectroradiometer (MODIS) LST data (2001–2020) were used. The pixel-based linear regression model was used to explore the spatiotemporal variability of LST changes. Furthermore, Sen’s slope and Mann-Kendall trend test were used to determine the magnitude of temporal shifts of the areal average LST and evaluate trends in areal average LST, respectively. Coefficient of variation (CV) was also used to analyze spatial and temporal variability in seasonal and annual LST. The seasonal LST CV varied from 1.096–10.72%, 0.7–11.06%, 1.29–14.76%, and 2.19–10.35% for average autumn (September to November), summer (June to August), spring (March to May), and winter (December to February) seasons, respectively. The highest inter-annual variability was observed in the eastern, northern, and south-western districts than that in the other parts. The seasonal spatial LST trend varied from −0.7–0.16, −0.4–0.224, 0.6–0.19, and −0.6–0.32 for average autumn, summer, spring, and winter seasons, respectively. Besides, the annual spatial LST slope varied from −0.58 to 0.17. Negative slopes were found in the central, mid-western, and mid-northern districts in annual LST, unlike the other parts. The annual variations of mean areal LST decreased insignificantly at the rate of 0.046°C year−1 (P<0.05). However, the inter-annual variability trend of annual LST increased significantly. Generally, the LST is tremendously variable in space and time and negatively correlated with elevation.


Introduction
Anthropogenic activities are becoming the principal cause for increasing greenhouse gas emissions which results in devastating environmental problems. The growing increment in population leads to increasing deforestation, forest degradation, agriculture, urbanization, and industrialization. This event has an extensive effect on ozone depletion which in turn boosts land surface temperature (LST). Besides, the rapid changes in global climatic conditions, biological diversity, and energy consumption are becoming common (Güneralp et al. 2017;Yang et al. 2020aYang et al. , 2020b, because a significant amount of natural and plantation forests, shrublands, and croplands are destroyed and replaced by artificial impervious surfaces which could be the source of enormous environmental influences on this planet (Miles and Esau 2020;Zhou et al. 2016;Fabeku et al. 2018;Jiang et al. 2015). Specifically, this situation causes a fall in evapotranspiration through photosynthesis, increase in runoff, and LSTs which in turn increase the occurrence of urban heat islands (UHIs) (Farina 2012;Ibitoye et al. 2017;Singh et al. 2014;Ye et al. 2017). UHI refers to an intensification in LST and has a great effect on urban climate, human well-being, biodiversity, and sustainable development (Chen et al. 2017;Qiao et al. 2020;Yang et al. 2020aYang et al. , 2020b. In short, UHI could be expressed as the situation where an urban environment has a higher temperature than a rural environment (Aina et al. 2017). Global warming happened because of the absorption of incoming radiation from the sun in a huge amount and re-radiating less while its cooling occurs when the outgoing solar radiation is significant. No country could be immune from the adverse global warming impacts since it has no geographical boundary. So, it requires the collective efforts of every nation.
LST is the crucial element for the investigation of climate in an urban area. Because it is very pertinent for the assessment of climate change and UHIs (Amiri et al. 2009;Khorchani et al. 2018;Mirzaei et al. 2020). It helps to evaluate the temperature dynamics of the earth's surface that are intensively correlated with climate change based on long-term data (Jiang et al. 2015;Peng et al. 2018;Zhu et al. 2018). Several studies stated that LST change is highly linked with land use/ land cover (LULC) change. When a significant amount of GHGs are released, naturally available GHGs would be denser so that an abundant amount of outgoing solar radiation would be blocked and will cause LST increment. Though benefits like getting warmth and obtaining a long time growing season for crops could be obtained from increasing temperature, it causes great socioeconomic and environmental challenges. Consistent assessments of LST discrepancies are vital for recognizing land surface energy budget and its interactions with the atmosphere properly (Maffei et al. 2018;Mildrexler et al. 2018;Thorne et al. 2016). LST could be expressed as an indicator of interactions among the components of the climate system. It further shows the thermal reaction to the urban coverage, height of construction, land surface coverage, and energy consumption (Yang et al. 2020a(Yang et al. , 2020bLi et al. 2013;Khandelwal et al. 2017;Miles and Esau 2020).
The temperature could be measured using traditional observation and remote-sensing techniques. Remote sensing is among the indirect approaches which are helpful to evaluate the LST (Shwetha and Kumar 2016). In recent times, it has become an effective method to get data at a large-scale spatial and temporal coverage for investigating the spatiotemporal changes in vegetation, land use/land cover, and LST (Berger et al. 2017;Li et al. 2013;Tomlinson et al. 2011). Moreover, it increasingly helps to assess biological, physical, and other features in this planet (Li et al. 2020;Qureshi et al. 2020). LST is an essential constituent of the climate system that is influenced prominently by interactions of the atmosphere and the earth's surface. The spatiotemporal variability in LST and its consequences are wide ranging across urban and rural environments. In equatorial and temperate regions, urban areas usually experienced higher surface temperatures than rural areas, unlike arid areas. This is mainly because of the intervention by humans and relatively low natural vegetation (Parvez et al. 2019). Several studies showed that remote sensing has the benefit of delivering uninterrupted, greatest quality, and real data of large areas as compared to meteorological stations data (Khanal et al. 2020;Li et al. 2013;Ngie et al. 2014). Ground-based meteorological station's LST data are scarce at a worldwide scale. Due to this, it is generally considered inadequate to examine its variability at a higher spatial scale (NourEldeen et al. 2020).
Numerous investigations have been conducted about the LST spatiotemporal variability using MODIS data in Asian and West African countries (Fabeku et al. 2018;Fathizad et al. 2017;Fonseka et al. 2019;Ibitoye et al. 2017;Phan et al. 2018;Qiao et al. 2020). However, limitations of such studies are available in Ethiopia, particularly over the Amhara region. Some researchers (Degefu and Bewket 2014;Dawit et al. 2019;Gebrehiwot and Veen 2013) investigated climate change using only time series data of air temperature and rainfall in the region so far although LST is a key indicator for environmental and climate change. Besides, the progress of conducting researches on LST variability is still very low in Ethiopia, specifically in the Amhara region. Nevertheless, deforestation is aggravating to expand agricultural lands because of the growing population in this region. This circumstance influences global warming. Moreover, a severe drought that affects crop and livestock production is being observed in the region. This investigation, therefore, would offer vital information about the LST dynamics and its link with elevation. It further helps for drought assessment. This investigation analyzed the spatiotemporal variability of LST and its link with elevation.

Materials and methods
The study area Amhara region is situated between 8°45′N to 13°45′N latitude and 35°46′E to 40°25′E longitudes (Fig. 1). The area coverage was expected to be 156,960 km 2 , and its altitudes were found between 513 and 4462 m a.s.l. Most districts are found above 1500 m a.s.l. According to Ayalew (2012), the yearly mean air temperature ranges from 15 to 21°C. Rain-fed farming is the predominant livelihood source of the communities. About 89% of the people are involved principally in mixed agriculture (Ayalew 2012).

Data and preprocessing
The MODIS (Moderate Resolution Imaging Spectroradiometer) data were used. According to Justice et al. (1998), MODIS refers to a multispectral device on the Aqua and Terra satellites that quantity high spatial resolution constituents of the climate system at 36 visible and ultraviolet frequencies with 0.4-14.4-mm spectral range. Nowadays, MODIS images are increasingly employed to develop numerous data for environmental investigations and monitoring, including LST products. A Terra product was used, because it contains longer time series data.
The monthly products on a 0.05°geographical grid (MOD11C3) were used although a 1-km daily spatial resolution data was available. This was done because the accessibility and trustworthiness of LST data increase with spatial and temporal aggregation (Bosilovich 2006;Li et al. 2018;Luintel et al. 2019). Version 6 MOD11C3 LST data (2001-2020) were downloaded from USGS LPDAAC FTP server (https://lpdaac. usgs.gov/data_access/data_pool). MOD11C3 MODIS LST data comprises both daytime and nighttime LST data. However, the monthly daytime MODIS LST data was used. Moreover, the quality flag in MOD11C3 was employed to screen the observations with a reduced quality. Preprocessing comprised subsetting, reprojecting, and eliminating false data points. These false data points incorporate pixels influenced by atmospheric disturbances including clouds. Several researchers pointed out the significance of screening cloud-contaminated data points in time series environmental and climatic data investigation (Julien et al. 2006;Julien and Sobrino 2010). Besides, the temperature data of Bahir Dar, Debre Berhan, Lalibela, Metema, and Chagni weather stations (2001-2014) were acquired from n is the length of the time series, and i is the number of years. Tmaxi and LSTi are the maximum temperature of weather stations and MODIS land surface temperature in the year i, respectively, and Tmax and LST are the mean maximum temperature of weather stations and the mean MODIS land surface temperature data, respectively the National Meteorological Service Agency (NMSA) of Ethiopia. The weather stations measure maximum temperature (Tmax) and minimum temperature (Tmin). However, only the maximum temperature (Tmax) was used due to its compatibility with daytime MODIS LST data (Phan et al. 2019). The weather stations' monthly maximum temperature (Tmax) data was used to evaluate the monthly daytime MODIS LST using some evaluation metrics like the Pearson correlation coefficient (r), mean absolute error (MAE), and root mean square error (RMSE) (Benali et al. 2012). The equations for MAE, RMSE, and r are presented in Table 1. Moreover, Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) was used and acquired from the United States Geological Survey (USGS) (https:/earthexplorer.usgs.gov/).

Trend analysis
The pixel-based linear regression model was used to evaluate the trends of LST changes pixel-wise in space and time using the seasonal and annual observations (2001-2020)    respectively (Qian et al. 2016;NourEldeen et al. 2020). In this investigation, the dependent and explanatory variables were LST and year, respectively. The pixel-wise analysis in the LST trend was calculated using Eq. 1 ( NourEldeen et al. 2020;Zhang et al. 2011).
where LSTi refers to land surface temperature in year i, n refers to length of time (n=20), and ti refers to an index number for 2001 to 2020 (1-20).
In addition, Sen's slope was computed to determine the magnitude of temporal shifts of the areal average LST. Several studies indicated that this approach is influenced to less extent by missing data and outliers unlike linear regression (Fernandes and Leblanc 2005;Porter et al. 2002). When a linear trend is available, the degree of the monotonous trend could be calculated and determined using Sen's slope estimator shown by Eq. 2 (Sen 1968): where β is the median of slope values between the yi and yj data measurements in phases i and j (i < j), respectively. The negative and positive values of β show downward and upward trends, respectively. Besides, the sign and value of β indicate the course and steepness of the trend, respectively. Furthermore, the non-parametric Mann-Kendall trend test method was executed to evaluate the trends in areal average LST values over the Amhara region. It is usually employed to assess monotonic trends of large time series climatic, environmental, and hydrological data. It is commonly affected to some extent when data is missing and distributed unevenly. Moreover, it is less vulnerable to outliers. This is for the reason that the ranks of observations are considered unlike their actual values (Libiseller and Grimvall 2002). The null hypothesis (H0) was there is no trend. Meaning, the Yi data are randomly ordered and tested against the alternative hypothesis (H1). The Mann-Kendall statistics (S) was calculated using the formulae developed by Mann (1945) and shown below by Eq. 3: where Yi and Yj are consecutive data values for n-length data and When the data are distributed identically and independently, Sen's slope estimator mean is zero, and Sen's slope estimator variance was computed using the formulae shown below by Eq. 5: where n refers to the length of data, m refers to tied group number in time series (a tied group refers to a group of sample data with similar value), and ti refers to the number of data points in the ith group. The Z statistics were computed using Eq. 6: To test the monotone trends, a significance level (α=0.05) was used. The decision for the two-tail hypothesis test was made by comparing the calculated Z with critical values. The H0 is accepted if the critical value is greater than the absolute value of calculated Z or if the P-value is higher than the selected significance level. On the contrary, the H0 is rejected if and only if the absolute value of the calculated Z is greater than the critical value. Once the H0 is rejected, the value of Z is negative and positive for declining and increasing trends, respectively (Hamlaoui-Moulai et al. 2013). Furthermore, it was considered statistically significant if the H1 is accepted.

Analysis of coefficient of variation (CV)
The LST variability for each pixel was analyzed seasonally and annually in space and time by computing the CV (Gidey et al. 2018;Muthoni et al. 2019). The analysis of CV was undertaken to evaluate the annual and seasonal LST variability relative to the mean percentage (2001-2020). The CV was calculated using the formulae shown below: where CV refers to the coefficient of variation value of LST, σ refers to the standard deviation of LST, and x is the long-term mean of LST.

Distribution of LST
The spatial distribution in average seasonal LST (2001LST ( -2020 spatially ranged from 43. 45-16.62°C, 39.89-14.59°C, 50.39-21.102°C, and 43.164-20.39°C for autumn, summer, spring, and winter, respectively (Fig. 2). During the spring and winter seasons, the highest LST was experienced in the western, north-western, and north-eastern districts. In opposition, the lowest LST was experienced in the southern and eastern districts. Moreover, the eastern, north-eastern, and north-western districts of the region showed the highest LST in the autumn and summer seasons. In contrast, the central, southern, and northern districts experienced the lowest LST (Fig. 2). The reason might be associated with most of the eastern, northeastern, and north-western districts had lower elevation (Fig.  3). However, most of the central, southern, and northern districts had higher elevations. This result is in line with the findings of several studies (Qiao et al. 2020;He et al. 2018;Phan et al. 2018;Peng et al. 2020) that revealed areas with higher elevation showed lower LST. LST for the spring season was greater than that for the other seasons. Besides, the average annual LST spatial distribution ranged from 41.976 to 18.565°C (Fig. 4). Figure 5 shows the seasonal LST distribution CV (%) in space. The seasonal LST CV ranges from 1.096-10.72%, 2.19-10.35%, 1.29-14.76%, and 0.7-11.06% for average autumn, summer, spring, and winter seasons, respectively (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). For the winter season, the maximum interannual variability in LST was observed in the eastern district.

Spatiotemporal variability in LST
In opposition, the northern, central, and western parts showed less inter-annual variability. Moreover, the maximum interannual variability was experienced in the eastern, south-western, and northern districts in the spring season. In contrast, less inter-annual variability was experienced in the north-eastern, north-western, and western districts. Besides, the eastern, southern, and north-eastern districtsshowed the maximum inter-annual variability. However, the north-western and north-eastern districts showed thelowest inter-annual variability in the autumn season. In opposition, the maximum interannual variability was experienced in the southern, central, and western parts in the summer season. However, less inter-annual variability was experienced in the eastern, north-western, and south-western districts. The CV of LST Fig. 6 Spatial distribution of annual LST CV (%) (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020) for the spring season (1.29 %< CV<14.76%) was greater than that for the other seasons. Besides, the maximum inter-annual variability of LST was recorded in the spring season than other seasons. Furthermore, Fig. 6 illustrates the annual LST CV (%) distribution in space. The spatial annual LST CV (2001-2020) ranged from 0.97 to 10.37%. The maximum interannual variability was experienced in the eastern, northern, and south-western districts. In contrast, less inter-annual variability was detected in the western, north-western, and northeastern districts. The changes in vegetation cover of land surfaces might be the reason for the maximum inter-annual variability in the eastern, northern, and south-western districts. Because, the earth surface's thermal properties with high vegetation are different from earth surfaces with low vegetation cover. The rise of vegetation cover is favorable to increasing soil moisture content, which can reduce soil erosion and desertification. In the equatorial and temperate regions, the areas with low vegetation cover are generally known for higher surface temperatures than the areas with high vegetation cover (Cai et al. 2016;Liang et al. 2017;Xiao et al. 2018;Ye et al. 2017). Figure 7 depicts the seasonal spatial LST trend of the Amhara region. The seasonal spatial LST trend varied from −0.6-0.19, −0.4-0.224, −0.7-0.16, and −0.6-0.32 for average autumn, summer, spring, and winter seasons, respectively (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). During the spring and winter seasons, negative LST slope values were experienced in the western, mid-northern, and central districts. However, positive slope values were Fig. 7 Spatial LST trend of winter season (a), spring season (b), summer season (c), and autumn season (d) (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020) experienced in the eastern, southern, south-western, and north-western districts. Besides, the negative LST slope values were also experienced in the southern, eastern, central, and mid-northern districts, while positive slope values were experienced in the western, south-western, north-western, and north-eastern districts in autumn and summer seasons. The annual spatial LST slope varied from −0.58 to 0.17 (Fig. 8). Negative slopes were found in the central, mid-western, and mid-northern districts. However, positive slopes were found in the north-western, north-eastern, southern, eastern, and south-western districts in annual LST. An increasing trend of annual LST in the north-western, north-eastern, southern, eastern, and south-western districts could be associated with the influences of anthropogenic activities like agricultural expansion, deforestation, and variation of elevation. Several studies (Berger et al. 2017;Jiang et al. 2015;Yang et al. 2020aYang et al. , 2020bFabeku et al. 2018;Qiao et al. 2020;Phan et al. 2018) stated a strong negative association is an available elevation and vegetation with LST.

Spatiotemporal trend
On the monthly basis, the trends of spatial average LST decreased insignificantly (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020) in all months except in January, March, August, and September (P<0.05) (Table 3). Similarly, the seasonal variability of areal average LST is presented in Fig. 9. LST decreased insignificantly at 0.036°C year −1 , 0.041°C year −1 , 0.074°C year −1 , and 0.005°C year −1 in autumn, summer, spring, and winter seasons at 5% significance level, respectively (Table 3 and Fig.  9). Moreover, the annual average spatial LST variability decreased insignificantly at 0.046°C year -1 . But, the inter-annual variability of autumn, summer, spring, and winter seasons LST increased insignificantly from 2001 to 2020. However, the inter-annual variability of annual LST increased significantly ( Fig. 10 and Table 4). Generally, the results revealed that there is a slight variation in LST as shown in Tables 3 and  4 and Figs. 9 and 10. But, this variation is not significant. The reason for exhibiting a nonsignificant LST trend could be associated with the spatial coverage and not considering most of the drought-affected lowland areas, out of the Amhara region in Ethiopia.     (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020) respectively. Furthermore, the annual spatial LST slope varies from −0.58 to 0.17. The western, mid-northern, and central areas show a general warming trend in the winter and spring seasons. Besides, the southern, eastern, central, and midnorthern parts show a warming trend during the summer and autumn seasons. In contrast, a general decreasing trend was found in the eastern, southern, south-western, and northwestern districts during the winter and spring seasons. A decreasing trend was also experienced in western, south-western, north-western, and north-eastern districts in the summer and autumn seasons. The central, mid-western, and midnorthern districts show a decreasing trend while the northwestern, north-eastern, southern, eastern, and south-western districts show an increasing trend in LST annually. The rate of decrement trend in spatial mean LST varies between 0.005 and 0.074°C year −1 . The spring season LST trend is insignificantly decreasing greater than other seasons at 0.074°C year −1 (P<0.05). The annual variations of mean LST decreased insignificantly at the rate of 0.046°C year −1 . Moreover, a negative correlation was observed between elevation and LST. The findings of this study would provide invaluable information to plan and implement appropriate adaptation and mitigation strategies by the experts of forestry, environment, and climate change. It is expected that different land use/land cover (LULC) types influence the LST. Therefore, an empirical study on the spatiotemporal dynamics of LST and its association with land use/land cover (LULC) change should be conducted. Furthermore, researches should be conducted by considering the spatial coverage of the whole regions of the country to provide more reliable information on the spatiotemporal LST dynamics in Ethiopia.