Spatial Conﬁguration and Extent Explains the Urban Heat Mitigation Potential due to Green Spaces: Analysis over Addis Ababa, Ethiopia

: Urban green space (UGS) is considered a mitigative intervention for urban heat. While increasing the UGS coverage is expected to reduce the urban heat, studies on the e ﬀ ects of UGS conﬁguration have produced inconsistent results. To investigate this inconsistency further, this study conducted a multi-spatial and multi-temporal resolution analysis in the Addis Ababa city metropolitan area for assessing the relationship between UGS patterns and land surface temperature (LST). Landsat images were used to generate land cover and LST maps. Regression models were developed to investigate whether controlling for the proportion of the green area (PGS), fragmentation, shape, complexity, and proximity distance can a ﬀ ect surface temperature. Results indicated that the UGS patches with aggregated, regular and simple shapes and connectivity throughout the urban landscape were more e ﬀ ective in decreasing the LST as compared to the fragmented and complicated spatial patterns. This ﬁnding highlighted that in addition to increasing the amount of UGS, optimizing the spatial structure of UGS, could be an e ﬀ ective and useful action to mitigate the urban heat island (UHI) impacts. Changing the spatial size had a signiﬁcant inﬂuence on the interconnection between LST and UGS patterns as well. It also noted that the spatial arrangement of UGS was more sensitive to spatial scales than that of its composition. The relationship between the spatial conﬁguration of UGS and LST could be changed when applying di ﬀ erent statistical methods. This result underlined the importance of controlling the e ﬀ ects of the share of green spaces when calculating the impacts of the spatial conﬁguration of UGS on LST. Furthermore, the study highlighted that applying di ﬀ erent statistical approaches, spatial scale, and coverage of UGS can help determine the e ﬀ ectiveness of the association between LST and UGS patterns. These outcomes provided new insights regarding the inconsistent ﬁndings from earlier studies, which might be a result of the di ﬀ erent approaches considered. Indeed, these ﬁndings are expected to be of help more broadly for city planning and urban heat mitigation.


Introduction
Over half the global population dwells in cities [1]. Urbanization is typically associated with providing better socio-economic conditions. On the other end, urbanization has also been associated with environmental and ecological stresses [2][3][4]. One of the widely recorded and evident impacts methods to accomplish this aim: (1) using a single city study to explore temporal dynamics [43,57], and (2) a comparative study of multiple cities applying similar approaches [45,53]. The present study adopts the first approach, to address the following questions: (i) does the spatial pattern of UGS affect temperatures differently in the city with changes over time? and (ii) do different statistical approaches or analytical units result in different effects of the spatial configuration of UGS on LST?
To answer these questions, the study analyzed the relationships between the LST and spatial configuration of UGS by using different statistical methods, and varied spatial configuration across different periods in the metropolitan area of Addis Ababa, Ethiopia. Landsat images were used to produce the land cover and LST maps. As cities have mostly limited space for greening [16,29,43], this study provides empirical results that can improve the understanding of the seemingly contradictory impacts of the spatial configuration of UGS patterns on LST noted in earlier studies. Such empirical results can potentially help planners and managers in regulating the spatial arrangement of urban green space to mitigate UHI. Furthermore, the adopted integrated approaches, multiple statistical, and landscape metrics analysis at varied analytical units, provide comprehensive and detailed information for environmental planning issues.

Study Area
The study area, presented in Figure 1, includes the city of Addis Ababa and its surrounding regions, which cover 1711.74 km 2 . It is located between 8 • 54 00" N to 9 • 8 00" N and 38 • 36 00" E to 38 • 54 00" E. The mean altitude of the study site is 2685 m above sea level. Addis Ababa is the capital, and most important commercial, political, and industrial city of Ethiopia. The city also serves as the seat of the African Union in addition to the Oromia National, Regional State. The city has a total area of about 527 km 2 and 4.3 million residents in the year 2017 [3], which cover 16.9 percent of the country's urban inhabitants and 3 percent of the country's overall inhabitants. The city's population is expected to grow to 5.9 million in 2030 at a mean annual growth rate of 4.1 percent [61]. The climate in the study site is warm and temperate. Typically, March, April, and May are the warmest months of the year, with values of 25 • C for the mean highest temperatures. November and December are the coldest months with the lowest average temperature values of 8 • C. The rainfall is moderate to heavy during the summer (June to August) with a mean yearly rainfall of 1143 mm. The dispersed nature of Addis Ababa city's spatial growth [3] and the rapid reduction of urban green space [62,63] are likely to have considerable impacts on spatial patterns of temperature [42,51,64]. Hence, Addis Ababa metropolitan area was an interesting case study to explore the environmental effects of urbanization.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 24 statistical approaches or analytical units result in different effects of the spatial configuration of UGS on LST?
To answer these questions, the study analyzed the relationships between the LST and spatial configuration of UGS by using different statistical methods, and varied spatial configuration across different periods in the metropolitan area of Addis Ababa, Ethiopia. Landsat images were used to produce the land cover and LST maps. As cities have mostly limited space for greening [16,29,43], this study provides empirical results that can improve the understanding of the seemingly contradictory impacts of the spatial configuration of UGS patterns on LST noted in earlier studies. Such empirical results can potentially help planners and managers in regulating the spatial arrangement of urban green space to mitigate UHI. Furthermore, the adopted integrated approaches, multiple statistical, and landscape metrics analysis at varied analytical units, provide comprehensive and detailed information for environmental planning issues.

Study Area
The study area, presented in Figure 1, includes the city of Addis Ababa and its surrounding regions, which cover 1711.74 km 2 . It is located between 8°54′00″ N to 9°8′00″ N and 38°36′00″ E to 38°54′00″ E. The mean altitude of the study site is 2685 m above sea level. Addis Ababa is the capital, and most important commercial, political, and industrial city of Ethiopia. The city also serves as the seat of the African Union in addition to the Oromia National, Regional State. The city has a total area of about 527 km 2 and 4.3 million residents in the year 2017 [3], which cover 16.9 percent of the country's urban inhabitants and 3 percent of the country's overall inhabitants. The city's population is expected to grow to 5.9 million in 2030 at a mean annual growth rate of 4.1 percent [61]. The climate in the study site is warm and temperate. Typically, March, April, and May are the warmest months of the year, with values of 25 °C for the mean highest temperatures. November and December are the coldest months with the lowest average temperature values of 8 °C. The rainfall is moderate to heavy during the summer (June to August) with a mean yearly rainfall of 1143 mm. The dispersed nature of Addis Ababa city's spatial growth [3] and the rapid reduction of urban green space [62,63] are likely to have considerable impacts on spatial patterns of temperature [42,51,64]. Hence, Addis Ababa metropolitan area was an interesting case study to explore the environmental effects of urbanization.

Remote Sensing Data and Processing
Landsat is regularly adopted for computing LST in UHI investigations [43] and is openly available. Three Landsat images for 1987, 1999, and 2019 acquired during the dry season (February) with clear atmospheric conditions were used in this study. These images were chosen based on the availability of the data for the same months. This month (February) was selected based on previous studies that reported the UHI effects during the dry season were more significant [34,59]. The images were downloaded from the United States Geological Survey. The details of the data implemented for the study are shown in Table 1. The multispectral bands were applied for classification of land-use/land cover (LULC) classes, while the thermal bands were used to calculate satellite brightness temperatures.

Land Cover Mapping
The general idea of remotely sensed data and the methodological method employed for this investigation are illustrated in Figure 2

Remote Sensing Data and Processing
Landsat is regularly adopted for computing LST in UHI investigations [43] and is openly available. Three Landsat images for 1987, 1999, and 2019 acquired during the dry season (February) with clear atmospheric conditions were used in this study. These images were chosen based on the availability of the data for the same months. This month (February) was selected based on previous studies that reported the UHI effects during the dry season were more significant [34,59]. The images were downloaded from the United States Geological Survey. The details of the data implemented for the study are shown in Table 1. The multispectral bands were applied for classification of landuse/land cover (LULC) classes, while the thermal bands were used to calculate satellite brightness temperatures.

Land Cover Mapping
The general idea of remotely sensed data and the methodological method employed for this investigation are illustrated in Figure 2 flow diagram.  Satellite data pre-processing, spectral, spatial, and radiometric rectifications were done using standard remote sensing methods [65][66][67]. The study then mapped the LULC classes based on 30 m resolution of Landsat images, bands 5-4-3 for Landsat 8 OLI, and bands 4-3-2 for Landsat-5. Depending on its spectral signature, each pixel in the image was allocated to the category that fits best. Ancillary data, such as Google Earth images, were used to help in classification. The five defined categories were: built-up area, bare land, water body, vegetation cover, and grassland. A supervised classification approach with the support vector machine (SVM) classifier algorithm was employed. The correctness of land cover classified results was validated by the accuracy assessment approaches [68][69][70]. The study further merged the vegetation cover, grasslands, and water body as UGS to analyze the relationships between the spatial patterns of UGS and land LST.

Land Surface Temperature Retrieval
Several factors affect LST, including surface emissivity [71,72], elevation [19,73], wind velocity, time of day, and the phenomena of precipitation [74], clouds [4,43,72], atmospheric situations, satellite sensor features, and the mathematical approaches used to retrieve LST [75,76]. It is not feasible to explicitly account for all these components in this study due to the following reasons: (1) the atmospheric data are not available before January 2000. (2) Though the atmospheric rectification based on satellite imagery using generic parameters is important, it has been revealed to the result in only a minor error for LST retrieval [77]. (3) This study implemented LST for the exploration of UGS patterns' effect on urban land temperatures.
where, TB values in Kelvin (K); K1, K2 = Thermal conversion constants for the band.

Estimation of the Land Surface Emissivity
Land surface emissivity (ε) is one of the essential parameters in retrieving LST. Surface materials have different emissivity due to their different viewing angles, materials, and roughness. Normalized difference vegetation index (NDVI) is one of the most commonly applied vegetation indexes. NDVI is calculated from the red (R) and near-infrared (NIR) bands-for Landsat 5, and 8, using Equation (4).
where, RED is a red spectral band, and NIR is a near-infrared spectral band of Landsat images. Hence, to compute the LSE, the NDVI-based emissivity approach proposed by Van de Griend and Owe [78] was applied in this study (Equation (5)).
where, ε s and ε v denote the emissivity of bare soil pixels and vegetation, respectively. Where P v indicates the fraction of vegetation [79], which was calculated using Equation (6): In the above, NDVI min and NDVI max are 0.2 and 0.5, respectively in a universal condition [80]. The term d ε includes geometrical distribution influence of the internal reflections and natural surface. As suggested by Valor and Caselles [81], the values of ε s and ε v are 0.960 and 0.985, respectively, for vegetation structures and unknown emissivity. The study also applied these emissivity values in the calculation. Moreover, for rough and heterogeneous surfaces, d ε can increase to 2% of emissivity [80] and the values of d ε are calculated as proposed by Sobrino et al. [82], using Equation (7): where, F is a shape factor with a mean value of 0.55, assuming different types of geometrical distributions [82]. Finally, the rectified land surface temperature maps ( Figure 3) were produced using the Equation (8) as:

The Spatial Pattern of Urban Green Space Analysis
Four landscape metrics: (1) Percent cover of green space (PGS) (among composition metrics), and three configuration metrics: (2) patch density (PD), (3) edge density (ED), and (4) mean nearestneighbor distance (MNN) ( Table 2) were applied to compute the spatial patterns of UGS. These metrics denote the main characteristics that can describe the spatial patterns of UGS, including composition, complexity, configuration, and fragmentation of patches. These metrics were chosen based on their prevalence in the scientific literature and environmental significance [2,16,39,49,53]. The spatial metrics are sensitive to scales in the spatial pattern explorations [16,45,57]. The dependency on the scale also impacts on the results of UGS and LST analysis [16,45]. Thus, the study

The Spatial Pattern of Urban Green Space Analysis
Four landscape metrics: (1) Percent cover of green space (PGS) (among composition metrics), and three configuration metrics: (2) patch density (PD), (3) edge density (ED), and (4) mean nearest-neighbor distance (MNN) ( Table 2) were applied to compute the spatial patterns of UGS. These metrics denote the main characteristics that can describe the spatial patterns of UGS, including composition, complexity, configuration, and fragmentation of patches. These metrics were chosen based on their prevalence in the scientific literature and environmental significance [2,16,39,49,53]. The spatial metrics are sensitive to scales in the spatial pattern explorations [16,45,57]. The dependency on the scale also impacts on the results of UGS and LST analysis [16,45]. Thus, the study used 9 varied sizes of analytical units starting from 120 m (which is equal to the pixel size of the TM thermal band) to

Correlation Analysis
A Pearson correlation matrix was built based on the correlation between LST and the spatial patterns of the UGS. A partial correlation analysis was performed to assess the configuration metrics (PD, ED, and MNN) and the LST by controlling for the influence of the PGS. Controlling for the influence of PGS is important, as it had the highest correlation with the configuration metrics as well as the LST across the study period in almost all tested spatial scales. Thus, the Pearson correlation evaluation is reliable to get a spurious connection between configuration metrics and LST.

Regression Analysis
For each of the studied period, multiple regression models were built to select the best model that describes the highest variance in LST. To avoid redundancy (inter-correlation) within the chosen spatial metrics [43,84], a multi-collinearity condition number (MCN) coupled with the variance inflation factor (VIF), for each metric was calculated. The values of MCN and VIF for each

Correlation Analysis
A Pearson correlation matrix was built based on the correlation between LST and the spatial patterns of the UGS. A partial correlation analysis was performed to assess the configuration metrics (PD, ED, and MNN) and the LST by controlling for the influence of the PGS. Controlling for the influence of PGS is important, as it had the highest correlation with the configuration metrics as well as the LST across the study period in almost all tested spatial scales. Thus, the Pearson correlation evaluation is reliable to get a spurious connection between configuration metrics and LST.

Regression Analysis
For each of the studied period, multiple regression models were built to select the best model that describes the highest variance in LST. To avoid redundancy (inter-correlation) within the chosen spatial metrics [43,84], a multi-collinearity condition number (MCN) coupled with the variance inflation factor (VIF), for each metric was calculated. The values of MCN and VIF for each independent variable were less than 30 and 7.5, respectively, which is acceptable for modeling [43].
The study then employed the ordinary least squares (OLS) for every built model, with the assumption of independent error terms. OLS is usually the initial step for diagnostic and to use as a benchmark for comparison. This method neglected the issue of spatial autocorrelation [45,53]. Spatial data contains spatial dependence and heterogeneity from spatial autocorrelation, which can amplify confounding in the analysis. Spatial clustering in error may come from spatial heterogeneity or spatial dependence. To address this, spatial auto-regression (SAR) models that incorporated spatial autocorrelation into the modeling are more suitable to explore the association between spatial patterns of UGS and LST [43,45,85]. Regarding SAR, the neighborhood interconnection of the response variable was estimated by a (n × n) framework of spatial weights and was incorporated in the ordinary multiple linear regression to consider spatial autocorrelation [86]. Accordingly, two spatial regression approaches: spatial lag and error models were applied to investigate the impacts of the spatial pattern of UGS on LST. The spatial lag model expects that the spatial autoregressive occurs only in the response variable [86]. The spatial lag model is expressed in Equation (9).
where Wy is the spatial lag variable, ρ is a spatial autoregressive coefficient, X is the independent regression variable, β is the regression coefficients, and ε is the error term vector. Conversely, the spatial error model considers the spatial impacts that are not wholly defined by the factors within the error terms and is presented in Equation (10).
where Wµ is a vector of spatially lagged errors, and λ is a coefficient of spatial autoregressive The study implemented the Lagrange Multiplier statistics to match the applied modeling methods for each year. The standard coefficients (beta weights) were used to examine the relative significance of composition and configuration metrics on estimating LST [43,45]. The regression model was run using GEODA version 1.8.16 [87]. confounding in the analysis. Spatial clustering in error may come from spatial heterogeneity or spatial dependence. To address this, spatial auto-regression (SAR) models that incorporated spatial autocorrelation into the modeling are more suitable to explore the association between spatial patterns of UGS and LST [43,45,85]. Regarding SAR, the neighborhood interconnection of the response variable was estimated by a (n × n) framework of spatial weights and was incorporated in the ordinary multiple linear regression to consider spatial autocorrelation [86]. Accordingly, two spatial regression approaches: spatial lag and error models were applied to investigate the impacts of the spatial pattern of UGS on LST. The spatial lag model expects that the spatial autoregressive occurs only in the response variable [86]. The spatial lag model is expressed in Equation (9).

Summary of Experimental Steps
where Wy is the spatial lag variable, ρ is a spatial autoregressive coefficient, X is the independent regression variable, β is the regression coefficients, and ε is the error term vector. Conversely, the spatial error model considers the spatial impacts that are not wholly defined by the factors within the error terms and is presented in Equation (10).
where Wμ is a vector of spatially lagged errors, and λ is a coefficient of spatial autoregressive The study implemented the Lagrange Multiplier statistics to match the applied modeling methods for each year. The standard coefficients (beta weights) were used to examine the relative significance of composition and configuration metrics on estimating LST [43,45]. The regression model was run using GEODA version 1.8.16 [87].    Figure 6 shows the changes in the land use pattern of the Addis Ababa metropolitan regions from the 1987 to 2019 study period. The area of non-urban areas (vegetation, grass, agriculture, and water) has rapidly been reduced across the study periods. In contrast, the area of the urban region has significantly increased over the past 32 years (Figure 6). The overall and kappa coefficients of LULC classes were 87.8, 88.2, and 89.4 in 1987, 1999, and 2019, respectively. The kappa values were 0.85 for 1987 and 1999, respectively, and 0.88 for 2019. These values indicated that the generated data had responsibly high accuracy.

LULC Changes in Addis Ababa Metropolitan Regions
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 24

LULC Changes in Addis Ababa Metropolitan Regions
Figures 6 shows the changes in the land use pattern of the Addis Ababa metropolitan regions from the 1987 to 2019 study period. The area of non-urban areas (vegetation, grass, agriculture, and water) has rapidly been reduced across the study periods. In contrast, the area of the urban region has significantly increased over the past 32 years (Figure 6). The overall and kappa coefficients of LULC classes were 87.8, 88.2, and 89.4 in 1987, 1999, and 2019, respectively. The kappa values were 0.85 for 1987 and 1999, respectively, and 0.88 for 2019. These values indicated that the generated data had responsibly high accuracy.  Table 3 shows the relative change patterns of UGS and the corresponding landscape metrics. The proportion of urban green space (PGS) in Addis Ababa metropolitan decreased gradually from 44.3% in 1987 to 33.4% in 1999, followed by a considerable decrease to 26.3% in 2019 (Table 3). In  Table 3 shows the relative change patterns of UGS and the corresponding landscape metrics. The proportion of urban green space (PGS) in Addis Ababa metropolitan decreased gradually from 44.3% in 1987 to 33.4% in 1999, followed by a considerable decrease to 26.3% in 2019 (Table 3). In terms of fragmentation and aggregation, the PD values of the UGS gradually increased between 1987 and 1999, and then sharply increased further in 2019. The result indicated that the UGS was disaggregated across the study periods, but UGS patch fragmentation was more significant in 2019, as evidenced by a sharp rise in the PD values: 7 in 1999 to 27 in 2019. The shape complexity of UGS, which was calculated by ED, revealed that changes in the overall complexity of the patches became more complex and irregular patterns from 1987 to 2019. This trend was high during the 1999-2019 periods, as evidenced by a rapid increase in the ED values, from 111.4 m/ha in 1999 to 2015.1 m/ha in 2019 (Table 3). Whereas, the connectivity among UGS patches increased notably from 1987 to 2019 and was generally increased by 41 m during this study period (from 42 m of MNN in 1987 to 91.36 m in 2019) ( Table 3). In general, it was inferred that Addis Ababa lost its green space from 1987 to 2019, at about 12% between 1999 and 2019. Concerning the spatial patterns, UGS patches have become smaller, more complex in shape, and more fragmented overall.

Effects of Spatial Patterns of UGS on LST
Across the study periods, the Pearson correlation analysis revealed that most metrics were significant at all analytical units, except MNN at large spatial scales (Table 4). Composition metrics (PGS) showed a negative correlation with LST, suggesting that the more UGS patches give a higher cooling effect. Among the configuration metrics, PD revealed a positive correlation with LST, indicating that the low fragmented UGS patches show a more effective cooling signature. In contrast, the ED had a negative correlation with LST across all analytical units. In contrast, MNN revealed a positive, significant correlation with LST only at analytical units of 120 m, 240 m, and 360 m. Similar to PGS, the intensity of the correlations between the three configuration metrics (PD, ED, and MNN) and LST also generally decreased with the increase of the analysis scale.
After controlling for the impact of PGS, the correlations (i.e., partial correlations) between LST and configuration metrics were significantly changed and are shown in Table 4. From the table, four features were notable: (i) the intensity of the partial correlations, calculated by the correlation coefficients significantly declined relative to their equivalent Pearson correlation coefficients. (ii) Some configuration metrics and LST became insignificant (e.g., MNN). (iii) Significantly correlated at only smaller scales (e.g., ED was significantly correlated to LST at less than or equal to 360 m, and reduced as the spatial magnitude increased), and (iv) More importantly, the association between some configuration metrics and LST shifted their directions. For instance, the correlation between MNN and LST were shifted from positive to negative when the analytical unit was between 360 and 960 m, 360 to 840 m, and 360 to 840 m in the year 1987, 1999, and 2019, respectively. In addition, the correlation of ED with LST was shifted from negative to positive across all analytical units except at 1080 m in the year 1987. Table 4. Correlation coefficients between Landscape metrics and LST. The bold rows indicate the partial correlation analysis, where for proportion of urban green space (PGS), the control variables were the configuration metrics, and for configuration metrics, the control variable was PGS.

Influences of UGS Spatial Size Pattern on LST
The correlation between the implemented landscape metrics and LST appears to be regularly decreasing with spatial size (Table 4). Whereas, some configuration metrics can influence on LST separately of the composition metrics at smaller scales. For example, the correlation of PD continued to be significant with LST at analytical units between 120 m to 600 m. This tendency is also noted when considering the influence of PGS across the study periods. The results imply that the patch density has a more prominent effect on LST than the other configuration metrics at smaller scales. It is noted that the optimum spatial scale of UGS for exploring the relationship between LST and spatial metrics could be between 120-600 m.
Indeed, this study did not find nor does it seem feasible to precisely identify a scale at which the correlation between LST and metrics converge. However, most metrics were generally more predictable at 240 m. For instance, after controlling for the impact of composition metrics (PGS), connectivity between the green spaces revealed a significant positive effect on LST only on the spatial scale of 240 m, as evidenced by MNN values in Table 4. Thus, 240 m × 240 m was applied as the optimum size for further analysis. Table 5 shows that the results from spatial error, Lag, and OLS linear regression models, the PGS had significant negative impacts on LST, across all the study periods. Additionally, PGS was the most significant indicator of LST, playing a considerably more significant role in forecasting LST relative to the other spatial configuration factors (Table 5). Among the three configuration metrics, PD played a much more significant and positive effect in estimating LST (Table 5). Other metrics, ED and MNN had positively significant effects on LST in the years 1985 and 1999. The results indicated that composition, fragmentation, and shape complexity characteristics of UGS spatial patterns affected LST. Moreover, the results showed that the spatial error model outperformed its correspondent Lag, and OLS models, for each of the periods. This was evidenced by the higher R 2 values and the lower AIC of the spatial error model compared to the other two models (Table 5).

Discussion
The relative importance of both spatial composition and configuration of UGS on LST was almost similar across the study periods. PGS was the critical factor in estimating LST. This is consistent with prior studies that suggested that the amount of green space plays a more important role in the cooling effects of the city than their configurations [42,45,46,85]. The overall results suggest that the fraction of green patches, their geometrics shape, degree of aggregation and fragmentation, and the proximity distance of patches to each other, are the main determinants of LST in urban regions. In other words, for the greater cooling effect, more patches of UGS distribution and a reasonable spatial configuration are important considerations to impact LST and reduce the urban heat.

Scale-Dependence of UGS Spatial Patterns and LST Relationships
The association between PGS and LST was negative at all spatial scales, across the study periods. This is consistent with the past results that showed that the amount of vegetation cover is a significant feature in mitigating UHI effects [32,45,53]. The study also found that increasing the spatial scale, the relationship between landscape metrics and LST, including both composition and configuration, becomes weaker, highlighting a scale effect. The results are consistent with the recent studies found elsewhere [46]. Furthermore, a nonlinear relation is noted between LST and vegetation cover-the LST reduction by enhancing the amount of vegetation cover was relatively lower at larger scales ( Table 6). This is likely due to higher diversity in green features of land covers and complex combinations typically occurring with the higher green cover than for small scales. Hence, smaller green fractions are likely less sensitive to the effect of changes in the proportion of vegetation cover on LST.  As noted, while the study did not find the exact spatial scale at which the greatest correlation has occurred, but it appeared that most metrics became more predictable at 240 m. This could be a robust feature considering the following three reasons, as demonstrated in this study. Primarily, after controlling for the effects of PGS, patch density of UGS revealed significant positive impacts on LST at the spatial extents between 120 m and 600 m in each studied period (see Table 4). Secondly, after controlling for the impact of composition metrics, connectivity between the green spaces revealed a significant positive effect on LST only on the spatial scale of 240 m in the year 1999, as evidenced by MNN values in Table 4. Third, such spatial scales were also suggested in earlier studies [28,43,57,58], though varied data used and approaches applied. Together this suggests that a spatial scale around 240 m could be an optimal spatial size to investigate the connection between UGS spatial patterns and LST.
The study results also revealed that spatial autocorrelation affected the association between LST and landscape metrics. Results from the spatial error and spatial lag models were higher than OLS in terms of R 2 , regression coefficients, and AIC values. This indicated that temperature in a specific region was affected by and affected the temperature in its surrounding samples (via advection). A similar finding was also suggested by [45,85]. The findings indicated that the LST at a location was not only dependent on the spatial arrangement of UGS in that region but also broader distribution in nearby areas.

The Spatial Composition and Configuration of UGS Effects on LST
While the amount of urban green space was the critical reason for cooling urban LST, urban greens arrangement was also associated with LST distributions. Arrangement of green space, including patch density, patches connectivity, shape complexity, and edge density, revealed positive effects on LST.
These spatial configurations of urban green spaces' effects on LST; however, varied with the earlier studies in terms of significance, direction, and magnitude. The positive association between PD and LST in this study, for instance, is consistent with studies [42,53,85], suggesting aggregated greenspace is better than a fragmented greenspace in reducing LST, while others reported the opposite [16,49,51]. This study found that all analytical units revealed a positive correlation of the green space patch density with LST. In addition, after controlling for the effects of a green space percent, the PD was positively significant to LST within 120-600 m across the study periods. Similarly, the contradictory results of UGS configuration effects on the LST in different cities were documented by several studies. Edge density of UGS, for instance, was found to be negatively associated with LST in Kansas, USA [42], Baltimore, USA [45], and Berlin, Germany [52], but positive in Singapore [53]. Rising of the edge density may potentially result in a rise of shade offered by vegetation to adjoining surfaces [51,85]. It may also increase energy flow and coupling between vegetation and their neighboring areas [56]. Thereby, taking into account only the process of shading, rising edge density will result in lower LST. Conversely, high edge density is mostly a result of more disaggregated UGS. As small and dispersed UGS has a higher temperature than that of aggregated and larger patches [43,53], which can decrease the effectiveness of surface albedo and evapotranspiration. This study also highlighted that the edge density of UGS had positive effects on LST across the study periods, suggesting that the complex shape and irregular patches of green spaces, the higher in LST. In other words, the UGS with simpler and regular patches leads to reduced LST.
The inconsistent results of configuration metrics reported by different researchers are likely due to the amount of UGS considered in the studies. For instance, in 1987, the amount of UGS was 44.3%. This figure reduced to 38.4% in 1999 and finally reached 26.3% in 2019 (see Table 3). At the same time, the effects of configuration metrics were changed in terms of direction, significance, and magnitude. For instance, the average nearest distances among each patch of UGS had a positive impact on LST in 1987 (see Table 4). However, in the case of data for 2019, it revealed insignificant, even changed its direction to negative at larger scales (see Table 4). This suggests that the amounts of green spaces can decide the effectiveness of configuration metrics to predict the relationship between LST and UGS patterns. Simultaneously, high fragmentation, complex shape pattern, and decreasing the connectivity among each patch observed in 2019 than the other periods. This may suggest that when patches are highly complicated and disconnected, the configuration's role could be declined. These outcomes can improve the understanding of the irregularity impacts of the spatial configuration of UGS on LST from earlier studies.
Furthermore, the irregularity impacts of configuration metrics on LST could be due to the spatial scale effects. The effects of configuration metrics were variable across the study periods. The average nearest distance between patches (MNN) was, for instance, notable in scale dependence with LST, where the relationships of MNN shifted from positive at a smaller spatial size to negative at the larger spatial size (see Table 4). That is, the correlation between configuration patterns of UGS and LST was sensitive to the spatial scales. The multi-scale analysis can assist resolve these inconsistencies. By exploring the relationship between UGS and LST pattern metrics under various analytical units, the insight into the impacts of UGS patterns on LST can be enhanced.
In general, considering the significant impacts of spatial size [16,28,43,45], the satellite image resolution [42,55], and the choice of statistical approach [51,57], comparing the UGS and LST relations across the studies could be challenging. This multi-scale and multi-year analysis with medium-resolution thermal data assists the interpretation and extension of prior that considered a single snapshot data in multiple cities, e.g., [39,45]. This could have the potential to identify the spatial composition and configuration relationships with LST due to its multiple periods in the same area to confirm its consistency.

Methodological Implications for Urban Greenspace Spatial Pattern Analysis
The study findings emphasize the need for controlling the influences of the percentage of green spaces when analyzing the spatial arrangement of UGS impacts on LST. For each studied period, after controlling for the influences of PGS, the association between LST and configuration metrics radically shifted, related to outcomes from the Pearson correlation analysis. The correlation between edge density and LTS, for instance, changed from negative to positive at all scales across the study periods. Additionally, the correlation between LTS and average nearest neighbor patch distance has shifted from positive to negative after or equal to 480 m. This is due to most of the configuration factors naturally correlated to PGS (Table A1; [45,85]). For instance, edge density had a significant negative connection with LST as per Pearson correlation analysis (r = −0.54, −0.49, and −0.67, p < 0.01 Table 3) in 1987, 1999, and 2019, respectively at the scale of 240 m. The correlation was due to the highly significant positive association between edge density and PGS (r = 0.83, 0.77, and 0.74, p < 0.01 Table A1) in 1987, 1999, and 2019, respectively. After controlling for the influences of PGS, edge density indeed had a significant positive connection with LST. Thus, it is important to apply a statistical approach such as multiple regression models and partial correlation, instead of Pearson correlation, to assess the relative contributions configuration and percent cover of green spaces to the LST. Applying Pearson correlation analysis alone could produce distorted results.
This research also suggested that the spatial pattern of UGS in Addis Ababa can be efficiently computed using only two metrics: PD and PGS, which represent configuration and composition features of spatial patterns. Each of these two metrics is the most important metrics related to their matching element, as distinguished by OLS, Lag, and Spatial error models. These metrics could be steadily applied to describe the temporal dynamics in LST distributions. These findings can help by minimizing the need to measure the redundant metrics as well as interpreting multiple metrics.
Furthermore, statistical methods, such as structural equation modeling and patch analysis, have been increasingly applied to detect the complexity and nested associations among social situations, land cover, landscape patterns, and surface temperatures [42,88]. Such assessments potentially assist in detecting the direct or indirect influences of the spatial structure of UGS on LST. Among the demonstrated spatial regression approaches, the spatial error model was more convenient than OLS in predicting LST, based on the pattern of UGS, as evidenced by their regression coefficient, R 2 , and AIC, across the study periods. This is especially relevant considering that the spatial patterns of UGS and LST are more likely to invalidate the assumption of independence among each patch, which did not occur in the frequently used statistical approaches such as OLS and Pearson correlation analysis. This study, therefore, again highlights that for spatial datasets, including the exploration of UGS spatial pattern effects on LST, the spatial regression approaches offer more robust evidence by incorporating spatial autocorrelation into the modeling to overcome the spatial clustering error from spatial dependence. These findings thus underlined the results of earlier studies [43,85] that the importance of the spatial regression model can overcome the concerns of spatial autocorrelation when exploring LST-UGS pattern relationships, which cannot be considered by OLS.

Conclusions
Urban green space has considerable cooling potential in mitigating urban heat. It is well-understood that increasing the amount of UGS can reduce the temperature of the city, but studies on the effects of UGS configuration have revealed inconsistent results. This study conducted multi-temporal, multiple statistical assessments, and varied scales to examine the connection between LST and green space patterns in Addis Ababa metropolitan area and further investigate the reasons for the inconsistency in prior studies. Landsat images were used to generate the LST and land cover maps. The main findings are summarized here.
(1) The UGS patches with aggregated, regular and simple shape and connectivity across the urban landscape are more effective in decreasing the LST than that of fragmented and complicated shape patterns. Thus, in addition to simply increasing the amount of UGS, optimizing the spatial structure of UGS, can be an effective and useful action to mitigate the UHI impacts. (2) The study also confirmed that the spatial pattern of UGS (in Addis Ababa) could be efficiently computed using only two metrics: PGS and PD, which represent, composition and configuration characteristics of spatial patterns. (3) Changing the spatial size had a considerable influence on the association between LST and UGS patterns. As the spatial scale increased, the influences of landscape metrics on LST, including both composition and configuration, become weaker, indicating spatial scale effects. (4) It is also noted that the configuration metrics (PD, ED, and MNN) are more sensitive to spatial scale than composition metrics (PGS). (5) The relationships between configuration metrics and LST could be changed when applying different statistical methods. This result underlines the importance of controlling the effects of PGS when calculating the impacts of the spatial configuration of UGS on LST. These outcomes provide additional understanding of the varying findings from earlier studies, likely confounded by the different approaches used (e.g., Partial correlation analysis versus Pearson correlation). (6) The relationships between LST and landscape metrics could be affected by spatial autocorrelation. In general, the study highlighted that applying different statistical approaches, spatial scale, and abundance of UGS can appropriately determine the effectiveness of configuration metrics to predict the association between LST and UGS patterns.
Urban heat island reduction is a potential mitigative strategy to decrease heat stress. Urban planning can increase the resilience of urban regions through enhancements to the green infrastructure and built environments. The findings of this study have potential implications for urban planners, suggesting more green space distribution with careful consideration of spatial configuration is important for effectively achieving the reduction in urban heating. Green infrastructure planning needs to consider both composition and configuration analyses for scenario planning. However, this study did not calculate the intensity of UHI, and future studies can incorporate it.
Appendix A Table A1. Expressive statistics of landscape metrics of UGS and LST.

Scale
Variables