Time Series Analyses and Forecasting of Surface Urban Heat Island Intensity Using ARIMA Model in Punjab, Pakistan

Abstract

In the context of rapid urbanization, Urban Heat Island (UHI) is considered as a major anthropogenic alteration in Earth environments, and its temporal trends and future forecasts for large areas did not receive much attention. Using land surface temperature (LST) data from MODIS (Moderate Resolution Imaging Spectro-radiometer) for years 2006 to 2020, we quantified the temporal trends of daytime and nighttime surface UHI intensity (SUHII, difference of urban temperature to rural temperature) using the Mann-Kendall (MK) trend test in six major cities of the Punjab province of Pakistan and estimated the future SUHII for the year 2030 using the ARIMA model. Results from the study revealed that the average mean SUHII for daytime was noted as 2.221 °C and the average mean nighttime SUHII was noted as 2.82 °C for the years 2006 to 2020. The average mean SUHII for daytime and nighttime exhibited increasing trends for all seasons and annually, and for the daytime spring season it showed a maximum upward trend of 0.486 °C/year (p < 0.05) and for the nighttime annual SUHII with an increasing rate of 0.485 °C/year (p < 0.05) which exhibited a maximum upward trend. The ARIMA model forecast suggested an increase of 0.04 °C in the average daytime SUHII and an increase of 0.1 °C in the average nighttime SUHII until 2030. The results from this study highlight the increasing trends of daytime and nighttime SUHII, ARIMA also forecasted an increase in daytime and nighttime SUHII, suggesting various strategies are needed for an effective mitigation of the UHI effect.

1. Introduction

Urbanization is a complex process that changes the cover of land from rural areas to industrial and urban areas. It is a vital component of land surface conversion and one of the most important demographic and spatial trends across the globe [1]. It produces distinguished spatial patterns that are influenced by different factors such as local physical characteristics and the transport network [2]. Currently, the urban population of the world is about 4.22 billion (55.3%) and by 2050 it will increase to 6.68 billion (68.37%) [3]. Rapid urbanization gradually replaces natural land cover by artificial features such as roads, roofs and hard structures, leading to an increase in nontranspiring and nonevaporating surface in cities [4]. This can have a direct effect on the energy of the surface of the land and have an impact on the biophysical properties of the surface of the land [5]. In particular, urbanization causes a decrease in the greenery required for humans in cities, which increases the heat that badly affects the hydrology of the environment and climate change [6]. One of the most prominent effects on urbanization is generating the urban heat island (UHI) effect. In 1833, Howard [7] initially presented the idea of UHI; UHI is a phenomenon in which the temperatures of urban areas are more than the temperature of its associated rural areas [8]. In recent decades, a large number of researchers have studied the negative effects of UHI on human life, ecosystem, reduction of biodiversity [9], climate change [10], destruction of vegetation [11], increase in the rate of disease and mortality [12], and reduction in water and air quality [13].

Five different methods, including fixed stations, remote sensing, mobile traverses [12,13], vertical sensing, and energy balances, can be used to measure the two types of UHI: (1) atmospheric UHI (AUHI) and (2) surface UHI (SUHI) [14]. AUHI utilized temperature data from ground weather stations [15]; all over the world, weather stations are unevenly distributed and they are very few in numbers, which causes the low observation density of results [16]. Therefore, it is very difficult to use data from the weather station for a large study area. On the opposite side, remote sensing technology (RS) can obtain continuous data for land surface temperature (LST) over a large study area. In 1972, Rao [17] performed the first analysis of SUHI and after that many researchers adopted RS data to perform SUHI analyses [18,19,20,21]. Generally, UHI showed its peak at nighttime, while the SUHI intensity (SUHII) showed greater values during the daytime [22].

Advanced techniques are necessary to comprehensively monitor the spatial distribution of UHI. These techniques are also used to measure periodic and dynamic changes in urban thermal environments [23]. LST products from Moderate Resolution Imaging Spectroradiometer (MODIS) were adopted by numerous researchers to study SUHI on a large scale due to the high temporal resolution and large area coverage of MODIS [24]. Studies from the past showed that SUHI have different temporal and spatial trends across the globe. For example, nine cities in central Europe showed high SUHII values on a monthly basis, and daytime SUHII values for the summer season showed the highest values [24]. The study by Roth et al. [25] showed that the SUHII at daytime was greater than the SUHI at nighttime for the cities of Los Angeles, Seattle and Vancouver. Zhou et al. [26] studied the SUHII for southern China, and the results showed that the daytime SUHII was greater than in other areas. The study by Yao et al. [27] focused on the temporal trends of SUHII in 31 cities in China for the study period 2001–2015, and the results showed the expansion and increase in SUHII. In addition, many studies focused on analyzing the fundamental factors influencing SUHII, for example, climatic conditions, heat emission by human activities, vegetation cover and urban areas [28]. A study showed that daytime SUHII for 419 cities globally had a significantly negative correlation with vegetation cover, and nighttime SUHII is insignificantly correlated with vegetation cover [20]. The contributions of these studies to understanding the primary causes of the SUHII effect are remarkable.

Forecasting models that are data-driven and can predict UHI in the future are necessary due to the growing impacts of UHI in major urban centers. Regression analysis [29], weighted regression analysis [30], principal component analysis (PCA), factor analysis (FA), canonical correlation analysis (CCA), cluster analysis (CCA), and spectral analysis [19] are all statistical methods frequently used in the development of UHI forecast models. However, most regression analysis is not good enough to capture nonlinear relationships [31,32]. Additionally, an artificial neural network model was used to predict weather; while highly accurate, the model can only make small adjustments to a series of data values [33]. Consequently, the Autoregressive Integrated Moving Average (ARIMA) model is used to predict the SUHI of the study area, as it improves the accuracy of the prediction while requiring fewer data inputs.

The diurnal and seasonal variations of the SUHII temporal variations were thoroughly studied. It was generally agreed that SUHII was at its peak during the summer. In general, daytime SUHII was stronger in warm cities compared to cold cities and in summer compared to winter. The SUHII did not change significantly at night as the city or season changed [34,35,36,37,38,39,40]. The interannual trend in SUHII has recently received more attention [18,19,41,42,43,44,45,46,47,48]. In some cities in China [8,14,41,43,48,49,50,51,52,53], India [54], Seoul Korea [31], Sri Lanka [55], Nigeria [56] Tehran City [57] and the Mediterranean [58], SUHII indicated a significant increasing trend. Similarly, very few researchers focused on studies related to SUHI in Pakistan [59,60]. A study of the expansion of SUHI in coastal areas of Pakistan has recently been conducted [60] in addition to this, only the cities of Karachi [61] and Lahore [62] have attracted attention to study SUHI patterns. Rizvi et al. [60] explored SUHII of coastal areas in Pakistan for the year 2017 and the result revealed that of five selected coastal cities Karachi showed the highest mean SUHII (2.60 °C), following by Ormara having the mean SUHII of 2.35 °C, the mean of Pasni city was recorded as 1.40 °C, Gwadar has the mean SUHII of 1.35 °C and Jiwani has the minimum mean SUHII (0.462 °C). A study analyzed the LST variations in the city of Lahore from 1996 to 2016 using Landsat data and according to the results of this study, the LST of Lahore increased by 4.8 °C and according to the prediction it will increase by 2 °C by the year 2035 [61]. Bilal et al. studied the UHI of Karachi city using MODIS data for years 2000–2020, the results showed that the mean SUHII for Karachi in the last 20 years (for months of January and May) is 0.15 °C [62]. Although previous studies concluded that SUHII variations typically occur on a variety of temporal time scales, such as diurnal variations, seasonal variations, and the interannual trend (i.e., temporal trends of SUHII are composed of seasonal variation and interannual trend, presenting the difference in days and nights in different cities). Using day- and night-time LST data from various cities, it is possible to analyze the diurnal variations of SUHII separately. However, the analysis of temporal variations is affected by the coupling of seasonal variation and annual trend in time series. So, to describe the temporal variations of SUHII more precisely, the seasonal variation and interannual trend must be separated. To create sustainable and resilient cities for SUHI, daytime and nighttime variations of SUHI must be considered. Understanding the SUHII trends can help urban planners avoid SUHI and urbanization-related issues when planning/designing urban growth. Therefore, we examined the temporal trends of daytime and nighttime SUHII and its forecast until 2030 in different cities of Punjab province of Pakistan by using the autoregressive integrated moving average (ARIMA) model. This study may help policymakers, urban planners, and other decision makers in their efforts to address the growing problem of SUHI in cities and develop sustainable urban designs. Furthermore, this research will be more helpful to those cities that may require additional consideration when planning and managing urbanization and SUHI mitigation strategies to reduce the rising urban surface heat island for a livable environment and sustainable cities.

This study aims to use the MODIS data for a detailed study of the SUHII phenomenon in the most populous province of Pakistan. To address the primary goal, precise objectives of this study are: (1) to evaluate the spatial distribution of SUHII for each city, (2) to measure the temporal trends of daytime and nighttime SUHII effect for each city, and (3) time series analyses and forecasting the urban LST, rural LST, and SUHII using ARIMA model.

2. Materials and Methods

2.1. Study Area

Punjab is the largest province of Pakistan, with a population of 109.98 million, of which 63.14% is the rural population and 36.86% is the urban population [63] and Punjab is the second largest province by area, occupying a total area of 205,344 km² [56] (Figure 1). Punjab showed a maximum population growth of 49.40% in Pakistan during the last Census-1998 [57]. Punjab has a subtropical climate with a humid summer and cold winter, Punjab received rainfall between 100 mm in the south and 600 mm in the northwest region [58]. The variation in temperature in Punjab is between −2 °C and 45 °C, but in some extreme conditions the mercury can reach 47 °C in summer and can drop to −5 °C in winter [60]. Punjab is the leading province in the economic and social development of the country and, with that, it is also the fastest growing province of the country. Due to this growth and urbanization, Punjab is facing various urban thermal issues. It is essential to address these thermal issues for sustainable urban management. The six main cities in Punjab were chosen to examine the spatial patterns and future trends of LST in this study, which are Lahore, Faisalabad, Gujranwala, Multan, Sialkot, and twin cities of Rawalpindi/Islamabad (RWP/ISB),

Table 1. Area, population, and population density of selected cities (according to the 2017 census).

City Name	Total Area (km2)	Total Population	Urban Percentage	Population Density (km−2)
Lahore	1772	11,119,985	100	6275.39
Faisalabad	5857	7,882,444	47.79	1345.82
ISB/RWP	6191	7,405,748	53.00	1616.715
Gujranwala	3622	5,011,066	58.85	1383.51
Multan	3720	4,746,166	43.38	1275.85
Sialkot	3016	3,894,938	29.39	1291.43

showed the details of these cities.

2.2. Data

The LST products of the MODIS satellite provided the ability to obtain continuous data with appropriate spatial resolution to differentiate between rural and urban areas. MODIS LST products are widely used to study SUHII on a large scale [63,64,65]. In this study, MODIS MOD11A2 with 1 km spatial resolution and 8 days temporal resolution [66,67] was used to evaluate LST trends and the effect of SUHII. To delineate the urban and rural area for each city annually, the MODIS MCD12Q1 composite land cover product with a spatial resolution of 500 m was chosen, the International Geosphere Biosphere Programme (IGBP) MCD12Q1 classification scheme provides 17 categories [68] To match the resolution of LST data, the MCD12Q1 product was rescaled to 1 km.

2.3. Methods

This study can be divided into three parts: (1) delineation of rural and urban areas by MCD12Q1 product, (2) calculation of SUHII using MOD11A2 LST data, (3) analysis of the trend of daytime and nighttime SUHII using the Mann-Kendall (MK) trend test (4) and prediction of daytime and nighttime SUHII for each city in Punjab using the autoregressive integrated moving average model (ARIMA). The complete methodology flow chart is illustrated in Figure 2.

2.3.1. Delineation of Rural and Urban Areas

For the delineation of rural and urban areas for each city in this study, first water pixels with a value of 17 pixels have been removed, as it is necessary to remove these pixels to eliminate the possible effect of these pixels on LST calculations [68,69]. Subsequently, urban pixels with the 13-pixel value were extracted for each city and the remaining pixels after the water and urban pixels in each city were considered rural areas [70].

2.3.2. SUHII Calculation

After defining the extent of the urban and rural area for every city, we use the following formula to calculate SUHII:

SUHII = Urban LST − Rural LST

(1)

where SUHII denotes the intensity of the urban heat island of the city, Urban LST represents the average LST from urban areas and Rural LST denotes the LST average of rural areas.

2.3.3. Mann-Kendall Test for Trend

Mann-Kendall (MK) is a nonparametric technique for identifying trend in time series data [71] MK used classifications of time-based data and treats each entry as a reference for comparison with all entries in time series [72]. Generally, MK is widely adopted to check trends in environmental time series data [73]. The statistics to calculate MK are given below:

$$S={{\displaystyle \sum }}_{p=1}^{n-1}{{\displaystyle \sum }}_{q=p+1}^{n}sgn\left(T_q-T_p\right)$$

(2)

$$sgn\left(T_q-T_p\right)=\{\begin{array}{c}if \left(T_q-T_p\right)<0; then -1\\ if\left(T_q-T_p\right)=0; then 0\\ if\left(T_q-T_p\right)>0;then 1\end{array}$$

(3)

where n is the number of data entries, T_p and T_q denote the successive entries in the time series for time p and q, sgn is the representation of the function that takes −1 if T_q − T_p < 0, takes 0 if T_q − T_p = 0, and takes 1 if T_q − T_p > 0. In the case of n > 10, the normal estimation is used (Z value of Kendall). To calculate the Z value, we need to calculate the variance of S VAR (S) [74].

$$VAR\left(S\right)=\frac{1}{18}\left[n\left(n-1\right)\left(2n+5\right)-{{\displaystyle \sum }}_{j=1}^i{y}_j\left(y_j-1\right)\left(2y_j=5\right)\right]$$

(4)

The equation considers the numerous tied entries (equal entries); where i represents the number of these equal trend values or groups, y_j is the representation of the number of entries presented in the jth group. After computing S and VAR (S) we can calculate the value of Z.

$$Z=\{\begin{array}{c}if S>0;then \frac{S-1}{\sqrt{VAR\left(S\right)}}\\ if S=0;then 0\\ if S<0;then \frac{S+1}{\sqrt{VAR\left(S\right)}}\end{array}$$

(5)

The positive value Z indicates the increasing trend and the negative value of Z indicates the decreasing trend in time series. Z values were tested at significance levels of 1% and 5%.

2.3.4. Sen’s Slope Estimator

The Sen slope estimate derived by Sen [75] has been a widely adopted method for estimating the slope of trends in time series [76] The representation of a linear model is as follows:

$$f\left(t\right)=St+C$$

(6)

where S is the slope, the function f(t) is the representation of time series, t is the time, and C is constant. S can be calculated by the following equation:

$$S_i=\frac{x_p-x_q}{p-q}$$

(7)

Here i = 1,2,3, …, n, x_p and x_q are the entries for times p and q and p > q. Calculating the median values of n of S_i can be done by the following equation [71].

$$S=\{\begin{array}{c}if n is odd then S_{\frac{n+1}{2}} \\ if n is even then \frac{S_{\frac{n}{2}}+S_{\frac{n+2}{2}}}{2}\end{array}$$

(8)

A positive value of S_i denotes the upward trend, a negative value of S_i denotes the downward trend, and zero value of S_i denotes that there is no trend in time series.

2.4. ARIMA Modeling

ARIMA (p, d, q) consists of two integration models; the first is the autoregressive (AR) model and the second is the moving average (MA) model. ARIMA is a widely adopted model to study time-series data. Meanwhile, LST data are time series data and in general time series data have trends in it that are non-stationary [77]. Removal of nonstationarity is performed by finite differencing of data points (d) in the ARIMA model [78]. In ARIMA, p represents the order of autoregression and q is the representation of the order of moving average [79]. To check the stationarity of the original data, the Dickey-Fuller unit root (ADF) test is a common method; if the test results suggest that the data are nonstationary, then differencing is applied to data to make in stationary [80]. When variance and mean of time series data are time-independent, then a time series is considered as stationary series [81]. The ARIMA (p, d, q) model can be presented by following mathematical equation:

$$\phi \left(L\right){\left(1-L\right)}^{d}yt=\theta \left(LT\right)\epsilon t$$

(9)

$$\left(1-{{\displaystyle \sum }}_{i=1}^p{\phi }_i{L}^i\right){\left(1-L\right)}^d{y}_t=\left(1+{{\displaystyle \sum }}_{j=1}^q{\theta }_j{L}^j\right)\epsilon t$$

(10)

where the order of the autoregressive, integrated, and moving average parts is represented as p, d, and q, respectively. However, ARIMA is not always useable for time series forecasting, in case of seasonal variations due to climatic data; seasonal ARIMA (SARIMA (p, d, q) is used instead of normal ARIMA, where index refers to the observations per season. Identifying the correct values for p and q is a challenging task in ARIMA. The autocorrelation function (ACF) and the partial autocorrelation function (PACF) facilitate this task, since these functions evaluate the relationship between the entries of time series. It is useful to plot ACF and PACF against their lags to identify the values of p and q [82].

The root mean square error (RMSE) and the mean absolute percentage error (MAPE) are the most common methods for measuring forecast accuracy in time series data. RMSE and MAPE can be calculated from the following formulas:

$$RMSE=\sqrt{\frac{1}{n}}{{\displaystyle \sum }}_{t=1}^n{\left(y_t+z^t\right)}^2$$

(11)

$$MAPE=\frac{1}{n}{{\displaystyle \sum }}_{t=1}^n\left|\frac{y_t-z_t}{x_t}\right|\times 100$$

(12)

where t represents the time, $y_t$ is representation of original values on time t, $z_t$ denotes the predicted value sat time t, sequence of time points denoted by n. The performance of the model can be measured by the values of RMSE and MAPE, a lower value for RMSE and MAPE shows better performance [83,84]. Validation is an important step to test the robustness of the ARIMA model. We also used the one-year (12 months) method to validate the ARIMA model. Figure 3 illustrates the complete flow for the ARIMA model.

3. Results

3.1. Distribution of the Average SUHII for the Last 15 Years

All selected cities showed positive values for nighttime SUHII from 2006 to 2020 in all seasons and annually, for daytime SUHII most cities showed positive values except for few cities (

Table 2. Averaged annual and seasonal SUHII for each city from 2006 to 2020.

City	Daytime SUHII (°C)					Nighttime SUHII (°C)
	Annual	Winter	Spring	Summer	Autumn	Annual	Winter	Spring	Summer	Autumn
Lahore	2.924	1.837	3.902	3.542	2.393	3.232	3.382	4.639	2.078	2.841
Faisalabad	3.139	1.609	3.915	3.854	3.147	3.213	3.540	3.809	2.139	3.377
ISB/RWP	−0.219	−0.923	−0.160	0.582	−0.416	1.042	0.750	1.268	1.121	1.024
Gujranwala	2.108	0.717	2.324	2.266	3.071	3.336	3.436	4.980	1.514	3.449
Multan	3.327	1.823	3.947	4.118	3.417	3.069	2.822	4.007	2.102	3.357
Sialkot	2.045	0.065	2.017	2.944	3.146	3.045	3.269	4.425	1.542	2.982
Average	2.221	0.855	2.657	2.877	2.460	2.823	2.866	3.855	1.749	2.838

). The maximum value of SUHII was observed at nighttime in the spring season (3.855 °C) and the minimum value of SUHHI was observed at daytime in the winter season (0.855 °C) for the averaged six cities. In case of individual city, Gujranwala showed the highest SUHII for nighttime spring (4.980 °C) and ISB/RWP exhibited the lowest value for SUHII for daytime winter (−0.923 °C).

3.2. Statistical Summary of SUHII

The statistical summary of the daytime and nighttime SUHII for each city is presented in

Table 3. Statistical summary of daytime and nighttime SUHII for Punjab cities.

City	Daytime SUHII (°C)		Nighttime SUHII (°C)
	Mean	SD	Mean	SD
Lahore	2.924	1.289	3.229	1.288
Faisalabad	3.139	1.304	3.214	0.916
ISB/RWP	−0.219	1.295	1.042	0.452
Gujranwala	2.108	2.476	3.332	1.541
Multan	3.327	1.633	3.068	1.023
Sialkot	2.045	2.154	3.041	1.427
Average	2.221	1.370	2.821	0.975

. Statistical results revealed that for the daytime SUHII Multan city has the highest mean value (3.327 °C) and ISB/RWP showed the lowest mean value of −0.219 °C. In case of nighttime SUHII Gujranwala showed maximum mean value of 3.332 °C and like daytime SUHII ISB/RWP showed minimum mean value of 0.452 °C.

3.3. Fifteen-Year Temporal Trends of SUHII for Punjab

Figure 4 showed the temporal trends of averaged daytime SUHII for selected six cities from 2006 to 2020; all results including annually, spring, summer, autumn and winter showed upward trends for SUHII, the highest significant upward trend is recorded in spring season (0.486 °C/year, p < 0.05) and the lowest significant increasing trend shown by annual data (0.410 °C/year, p < 0.05). Temporal trends for nighttime were shown in Figure 5; nighttime trends for SUHII are similar to daytime SUHII and trends from all four seasons and annually are recorded as upward. The highest significant increasing trend was recorded annually (0.485 °C/year, p < 0.05) and the minimum significant increasing trend was shown by the summer season (0.352 °C/year, p < 0.05). For daytime SUHII, the summer and autumn seasons exhibited insignificant increasing trends and for nighttime SUHII, two seasons (winter and spring) exhibited insignificant increasing trends.

Table 4. The last 15 years of temporal trends of SUHII for Punjab (Z represents the MK value and S represents the Sen’s slope value).

City	Daytime SUHII (° C/Year)						Nighttime SUHII (°C/Year)
		Annual	Winter	Spring	Summer	Autumn	Annual	Winter	Spring	Summer	Autumn
Lahore	Z	0.048	0.029	0.219	−0.105	−0.448 *	0.295	0.048	0.162	0.181	0.371 *
	S	0.002	0.007	0.041	−0.016	−0.040	0.032	0.016	0.036	0.028	0.050
Faisalabad	Z	0.124	0.333	0.124	0.067	0.124	0.543 **	0.181	0.390 *	0.371 *	0.486
	S	0.016	0.056	0.023	0.008	0.007	0.039	0.037	0.039	0.039	0.034
ISB/RWP	Z	0.581 **	0.314	0.543 **	0.410 *	0.562 **	0.486 *	0.219	0.333	0.162	0.314
	S	0.095	0.037	0.178	0.114	0.092	0.025	0.016	0.034	0.013	0.017
Gujranwala	Z	0.276	−0.124	0.295	0.124	0.124	0.314	−0.048	0.390 *	0.371 *	0.314
	S	0.031	−0.015	0.057	0.042	0.015	0.030	−0.003	0.077	0.036	0.049
Multan	Z	−0.105	0.200	0.086	−0.219	−0.371 *	0.314	0.067	0.143	0.162	0.352 *
	S	−0.014	0.054	0.016	−0.060	−0.094	0.025	0.008	0.024	0.023	0.042
Sialkot	Z	0.371 *	0.371 *	0.257	0.314	−0.086	0.333 *	−0.067	0.257	0.333	0.371 *
	S	0.040	0.062	0.054	0.063	−0.006	0.030	−0.021	0.073	0.056	0.041
Average	Z	0.410 *	0.448 *	0.486 *	0.276	0.086	0.485 *	0.067	0.314	0.352 *	0.390 *
	S	0.031	0.036	0.061	0.027	0.006	0.032	0.006	0.043	0.029	0.034

** p < 0.01 and * p < 0.05.

presented details of the temporal trends of daytime and nighttime SUHII for each city; for daytime annual SUHII, eight of the nine cities showed an increasing trend and only Multan exhibited an insignificant decreasing trend. ISB/RWP showed the highest significant increasing trend at a rate of 0.583 °C/year (p < 0.01) and Sialkot exhibited the lowest significant increasing trend of 0.371 °C/year (p < 0.05). For winter daytime SUHII only Gujranwala exhibited a downward trend and all other cities exhibited upward trends, the highest upward trend was recorded in Sialkot (0.371 °C/year, p < 0.05). Meanwhile, increasing trends in SUHII during the daytime were observed in all cities for the spring season and the maximum significant increase trend was observed in ISB/RWP (0.543 °C/year, p < 0.05). Additionally, two cities (Lahore and Multan) exhibited downward trends and all other cities showed upward trends for the summer season and a maximum significant upward trend was observed in ISB/RWP (0.410 °C/year, p < 0.05). Lastly, for the autumn season, six cities showed increasing trends for daytime SUHII, the highest significant increasing trend was shown in ISB/RWP with an increasing rate of 0.562 °C/year (p < 0.01).

All cities showed upward trends for the annual, spring, summer and autumn season, and only two cities (Gujranwala and Sialkot) in the winter season showed downward trends for the nighttime SUHII. For the annual nighttime SUHII the highest significant upward was recorded in Faisalabad with an increasing rate of 0.543 °C/year (p < 0.01). The winter season showed all insignificant trends for the nighttime SUHII. Meanwhile, a significant increase trend was observed for the spring season in Faisalabad with a rate of 0.390 °C/year (p < 0.05). Additionally, Gujranwala showed a significant increasing trend (0.371°C/year, p < 0.05) for the summer season. In the last, Lahore showed a maximum significant upward trend with a rate of 0.371 °C/year (p < 0.05) and Multan showed a minimum significant increasing trend (0.352 °C/year, p < 0.05).

3.4. ARIMA Model for Daytime and Nighttime SUHII

As SUHII was calculated from LST, which is a seasonal phenomenon, and plots of the significant ACF and PACF confirms at multiple lags of 12, so SARIMA with seasonal differencing of 12 was adopted to perform forecasting of daytime and nighttime SUHII. Initially, the ADF unit root test was applied to check the stationarity of the data and then the values for p and q along with the seasonal order P and Q were recognized as SARIMA (0,0,1) (1,1,1)12 to estimate the SUHII daytime and nighttime for each city. Values for RMSE and MAPE were calculated by testing data for the individual model.

Table 5. RMSE, MAPE and p(ADF) for the fitted ARIMA Model.

City	Daytime SUHII (°C)			Nighttime SUHII (°C)
	p-Value (ADF)	RMSE	MAPE	p-Value (ADF)	RMSE	MAPE
Lahore	<0.0001	0.74	0.24	<0.0001	0.66	0.19
Faisalabad	<0.0001	0.77	0.3	<0.0001	0.53	0.14
ISB/RWP	<0.0001	0.75	2.19	<0.0001	0.38	0.63
Gujranwala	<0.0001	0.82	0.76	<0.0001	0.66	0.56
Multan	<0.0001	1	0.46	<0.0001	0.64	0.22
Sialkot	<0.0001	0.91	1.39	<0.0001	0.68	0.32
Average	<0.0001	0.43	2.98	<0.0001	0.43	0.12

showed the values of RMSE, MAPE and p for the ADF test.

Forecasting was estimated for monthly daytime and nighttime SUHII averages after selection of the best-fit SARIMA model. Monthly average data of daytime and nighttime SUHII for each city from January 2006 to December 2019 were taken as validation data for the model, and then forecasting was done on monthly average of SUHII from January 2020 to December 2030. Figure 6 showed the forecasting values for the average SUHII during the daytime for each city and Figure 7 showed the forecasting values for the average SUHII during the nighttime for each city with 95% upper and lower confidence level. All forecasted values seem very close to the observed values, which confirmed the good forecasting of the selected SARIMA model. The forecast of the averaged daytime and nighttime SUHII endorsed the increasing trends of SUHII, as the averaged daytime SUHII of six cities for year 2030 was forecasted at 2.37 °C, which is higher than the average daytime SUHII of six cities for year 2020 (2.33 °C) and the averaged nighttime SUHII of six cities for year 2030 was forecasted at 3.00 °C, which was higher than the SUHII of 2020 (2.99 °C).

4. Discussion

Statistical and spatial distributions of daytime and nighttime SUHII for six cities in Punjab were analyzed on a seasonal and annual scale for years 2006–2020. The average daytime SUHII for all six cities varied from daytime SUHII 0.855 °C in winter to nighttime SUHII 3.855 °C in spring. The average mean for the daytime SUHII for the six cities was recorded as 2.221 °C and the average nighttime SUHII for the six cities was recorded as 2.82 °C.

Temporal trends for daytime and nighttime SUHII were evaluated using MK trend analyses for every city from 2006 to 2020. Temporal trends for the average daytime and nighttime SUHII of six cities exhibited upward trends annually and seasonally. These results of increasing SUHII are similar to the results of recent study in Islamabad [52] that showed an increase in the effect of UHI due to the change in LULC, increasing SUHII was also observed in recent study in metropolitan cities of Pakistan [58,85]. Furthermore, the increasing trends of SUHII are similar to the results of Kuwait [86], Penang (Malaysia) [87], the mega-cities of Bangladesh [88,89,90,91,92] and Wuhan (China) [93]. The highest significant upward trend for daytime SUHII was observed in the spring season (0.486 °C/year, p < 0.05), for nighttime the highest significant upward trend was shown by annual SUHII (0.485 °C/year, p < 0.05). Almost every city showed an increasing trend for daytime and nighttime SUHII for all seasons and annually; maximum increasing trend for annual daytime SUHII was observed in ISB/RWP at a rate of 0.583 °C/year (p < 0.01), Sialkot with an increasing rate of 0.371 °C/year (p < 0.05) exhibited the highest upward trend for daytime SUHII of winter season, the largest significant rising trend was observed in ISB/RWP (0.543 °C/year, p < 0.05) daytime SUHII for spring season, for summer season daytime SUHII maximum significant upward trend was observed in ISB/RWP (0.410 °C/year, p < 0.05) and for autumn season daytime SUHII the highest significant rising trend was shown by ISB/RWP with an increasing rate of 0.562 °C/year (p < 0.01). For the annual SUHII nighttime Faisalabad showed a maximum significant upward trend with an increase rate of 0.543 °C/year (p < 0.01), no city showed a significant trend for nighttime SUHII in the winter season, for the spring season the highest significant increasing trend for nighttime SUHII was shown by Faisalabad (0.390 °C/year p < 0.05), Gujranwala presented a maximum significant rising trend with the rate of 0.371°C/year (p < 0.05) for nighttime SUHII in the summer season, and the largest upward trend for nighttime SUHII in winter season was recorded in Lahore (0.371 °C/year (p < 0.05). Increasing trends in SUHII indicate that man-made structures are causing the increase in urban LST. SARIMA (0,0,1) (1,1,1)12 was fitted for the monthly daytime and nighttime SUHII of six cities in Punjab from 2006 to 2020 after checking the stationarity of the data. The ADF test was adopted [86] to check the stationarity of the data giving p values <0.0001 for every city, ACF and PACF were applied to achieve the values for p, d, and q to set the SARIMA model [84]. The RMSE and MAPE values confirmed the fitness of the SARIMA model. The SARIMA model confirmed the increasing trends of SUHII at daytime and nighttime for six cities of Punjab and forecast higher SUHII values in 2030.

The increasing trends of daytime and nighttime SUHII by the MK trend test is pointing to an increase in urban LST in cities of Punjab which can be related to the decrease in vegetation due to urbanization [19,63]. Forecasting of the increase in SUHII and the trends of increasing SUHII in all cities of Punjab need the attention of policy makers for future planning of cities and management of cities with higher values of SUHII such as ISB/RWP.

4.1. Implications

In recent years, a noticeable increase in daytime and nighttime SUHII has been observed in Pakistan, and this trend is expected to continue. Although the effects of urbanization on UHI are inevitably, there are mitigation measures that can be taken in planning, management, and policy making to reduce surface temperature and make the urban thermal environment more livable. Pavement is a significant element of the city and occupies a significant portion of the urban area; as a result, it contributes significantly to the formation of the urban heat island effect (UHI) [94]. The utilization of cool pavements, which have a lower surface temperature and a lower rate of sensible heat flux to the atmosphere, is one of the most effective mitigation strategies that can be used to combat the UHI effect. It is more effective to fully implement the strategy of cool and evaporative roofs in high-density urbanized areas, while in low-density urbanized areas, there is considered more space for planting, which makes street trees an effective strategy for mitigating the effects of the urban heat island effect [62]. In addition, urban reflective materials that prevent the direct transformation of incident solar radiation into sensible heat have the potential to reduce the impact of UHI [65,95]. In general, the method proposed in this study can help us evaluate the efficacy of mitigation measures by characterizing the temporal variations of historical and future UHI. This allows us to determine whether the rising trend of UHI is changing when mitigation measures are implemented, which in turn provides decision makers with technical support.

4.2. Limitations and Future Work

MODIS data (MOD11A2 and MCD12Q1) were used in the study to examine the area of urbanization and its effects on SUHI for cities of Punjab. Researchers widely used MODIS data to examine spatiotemporal trends of SUHII [21,72,96], especially for a large study area, processing of MODIS data is convenient, which makes it a good choice. Analyzing the effects of urbanization on SUHII for Punjab using higher spatial resolution data (e.g., Landsat) can be very time consuming. Furthermore, data with higher spatial resolution have low temporal resolution. Future studies should adopt higher spatial resolution, which can ensure more accuracy than MODIS data. Higher spatial resolution data and multisensor fusion (MODIS and Landsat) can be used to investigate urbanization effects on VC in future studies [97,98]. Additionally, this study adopted the ARIMA model for forecasting which performed well; however, future studies should focus on machine learning models for forecasting. The present study revealed the significant effects of uncontrolled urbanization on SUHII in Punjab cities, however future studies should focus on the consequences of this effect on human health and quality of life in these cities. Future studies also can focus on each city separately and consider more variables like land cover, land use and vegetation index.

5. Conclusions

In this study, trends and forecasting for daytime and nighttime SUHII for six cities in Punjab were evaluated for the study period 2001 to 2020 using the MK trend test and the ARIMA model, respectively. The following conclusions were drawn from the results of this study:

The average mean for the daytime SUHII for the six cities was recorded as 2.221 °C and the average nighttime SUHII for the six cities was recorded as 2.82 °C.

The averaged daytime SUHII showed upward trends for all seasons and annually, spring season with an increasing rate of 0.486 °C/year (p < 0.05) showed a maximum significant rising trend and annual SUHII with an increasing rate of 0.410 °C/year (p < 0.05) showed a minimum significant rising trend.

The average nighttime SUHII for all six cities exhibited an increase for all four seasons and annually. Annual SUHII with an increasing rate of 0.485 °C/year (p < 0.05) exhibited highest significant increasing trend and summer SUHII with an increasing rate of 0.352 °C/year (p < 0.05) exhibited lowest significant increasing trend.

The ARIMA model forecast an increase of 0.04 °C in the average SUHII at daytime of six cities of Punjab and an increase of 0.01 °C in the nighttime SUHII for all six cities of Punjab.

The increase in SUHII can cause a decrease in vegetation in urban areas of Punjab. Furthermore, increasing trends for SUHII in most cities of Punjab indicate a rapid and unplanned urbanization that requires immediate attention of policy makers.