Rehabilitated Coal Mine Tailings as Green Cooling Islands in the Metropolitan Ruhr Area (Germany)

1. Introduction

The urban heat island (UHI) effect is one of the major current and future threats to human health as cities sprawl and densify [1,2,3]. Thus, there is an extensive body of literature on available strategies to reduce the UHI effect and on improving environmental justice in city neighborhoods [4,5,6]. In this context, the cooling effects from urban greening has often been highlighted [7,8,9]. The role of urban green spaces (UGSs) in UHI mitigation has been shown in numerous thermal remote sensing analyses [7,10,11,12,13,14], successfully relating the vegetation abundance reflected by the Normalized Difference Vegetation Index (NDVI) to the land surface temperature (LST). Thus, it is now well accepted that vegetation covers decrease the land surface temperature by up to 3°C compared to the temperature in the surrounding built-up areas, whereby larger UGS areas (>50 ha) showed more substantial cooling effects than smaller green patches [15,16]. In a study by Kuang et al. [17], an LST difference of even over 6°C was found between UGS and built-up areas in some cities.

Generally, the impact of vegetation on the LST is attributed to the fact that plants drastically increase evaporative cooling processes, canopy shading as well as heat adsorption at a local scale [18,19,20]. Nevertheless, more and more research has systemized these impacts and generally differentiate between configuration and composition characteristics controlling the UGS cooling potential [21,22,23]. While configuration comprises of UGS size, shape, complexity, connectivity, and fragmentation, composition defines surface cover and relative abundances of landscape types of the UGS [24,25,26]. While it is still not possible to develop a consensus on specific ways in which varying aspects of UGS configuration interact with land surface temperature (LST) [21], compositional studies clearly showed that areas with a higher percentage of forest vegetation show a greater cooling effect than those dominated by shrub and grass vegetation forms [22,27,28]. Consequently, some studies also focused on tree characteristics such as species, age, or crown diameter and its effect on cooling [14,28]. For example, Rahman et al. [29] observed that fast-growing tree species such as Prunus umineko or P. calleryana generally show a higher stomatal conductivity, providing more cooling. Stumpe et al. (2024) observed that tree age controlled surface cooling, whereas Helletsgruber et al. [30] found that trunk circumference is a valuable indicator for estimating climate-regulating ecosystem services. Similarly, the vertical structure of trees was found to affect their cooling potential [31,32] and He et al. (2024) observed a significant negative correlation between tree canopy height and LST (r = -0.83).

Besides configuration and composition, other factors, such as water availability or terrain attributes, have been discussed in the context of the LST patterns. Concerning water availability, a strong negative correlation has often been found between vegetation moisture reflected by the Normalized Difference Moisture Index (NDMI) and LST pattern across UGS patches, indicating the high relevance of available water for evaporative cooling processes [14,33,34,35]. In this context, the Temperature Vegetation Dryness Index (TVDI) became a well-established parameter for monitoring surface soil moisture [36,37,38,39].

Concerning terrain attributes, several studies characterized the relationship between LST patterns and terrain attributes in natural environments. For example, Tan et al. [34] found a positive correlation between LST and altitude gained from a digital elevation model (DEM) in the Dongting Lake area in China. Karbalaee et al. [40] investigated the relationship between the terrain parameter and LST and observed that an increase in altitude, aspect, and slope leads to a significant decrease in LST. Also, Bai et al. [41] found a relationship between the spatial LST pattern in the Siming Mountain in China and common terrain attributes, although three-dimensional characteristics of the green space showed a higher relevance in that context. Nevertheless, in urban environments, the role of terrain characteristics in the cooling effect of greenspace patterns has generally been neglected.

The majority of these investigations focused on UGSs such as urban parks, urban forests, or urban gardens, with parks being most relevant as urban cooling islands. Consequently, the term “park cooling island (PCI)” arose [18,42,43]. More recent studies also focused on other UGS, such as peri-urban greens or cemeteries [14,44]. However, in the Ruhr Metropolitan Area (Germany), one of the biggest urban agglomerations in Germany, numerous tailing piles from the former coal-mining industry characterize the landscape. The tailing piles are man-made hills often located in residential areas systematically revegetated for local recreation several decades ago. Unlike typical park cooling islands, piles form sloped green areas with high relief energieswhich can thus provide additional cooling effects through drainage of cool air into adjacent residential areas during night time [45]. Although such green tailing piles can potentially play an important role in neighborhood cooling, these green spaces have not been previously considered in the context of UHI mitigation strategies.

Therefore, this study aims at systematically characterizing and understanding the surface cooling of the tailing piles in this region in order to assess their potential as cooling islands. In detail, the objectives of this research were:

(I)

Characterization of the summertime LST footprints of the piles compared to other common UGSs in the Ruhr Metropolitan Region.

(II)

Understanding mean summertime LST values of the piles in the context of vegetation and terrain attributes using the k-mean classification procedure.

(III)

Understanding pixel-based summertime LST values of the piles in the context of vegetation, soil and terrain attributes using random forest regression modeling.

2. Materials and Methods

2.1. Study Area

We focused on the Ruhr Metropolitan Region in North Rhine-Westphalia (Germany) since it is the largest urban economic zone in Germany with about 5.5 million inhabitants and a total area of about 4,400 km² [46,47].

The area covers the old-industrialized Ruhr area, one of Europe's most important coal fields, where millions of tons of coal were extracted through deep mining in the past century. As a by-product of these activities, large amounts of mine tailings were produced during shaft drilling and coal washing , consisting of sandstone, shale or slate rocks [48]. Most of the mining waste material was dumped in tailing piles, with only smaller amounts being used as underground stowing materials. In the early 1980`s, the regional mining company (Ruhrkohle AG) developed a tailing pile rehabilitation concept, which was integrated into the regional development plan. In this concept, the tailing piles were designed to blend optimally into the surrounding landscape, assuring the protection and restoration of ecological functions [48].

Before establishing a vegetation cover, tailing piles were generally covered with shallow soils (up to 10 cm) and, if necessary, deep soils (up to 1.8 m) for improving water retention and nutrient availability. In 1985, the greening concept [49] was implemented, where shrub and three species were selected and planted depending on soil cover characteristics. In most cases, the pile boundary areas were planted with blackthorn (Prunus spinosa), buckthorn (Hippophae rhamnoides), willow (Salix sp.), and common maple (Acer pseudoplatanus). The main areas were vegetated with red oak (Quercus rubra), wild-growing fruits, birch (Betula sp.), rowanberry (Sorbus aucuparia), and horse chestnut (Aesculus hippocastanum).

Within this study, we focused on tailing piles and their surface cooling potential, but also analyzed the thermal footprints of more common UGSs (parks, allotments, and cemeteries) as a reference. To delineate parks (P), allotments (A), or cemeteries (C) from their surroundings, we used the ‘Authorative Topographic-Cartographic Information System’ (ATKIS). It is the official German nationwide open source digital database for topographic spatial data, updated annually [50]. However, since within the ATKIS database, tailing piles (TP) are not assigned to one ATKIS object group, we used a digital map provided by Heiko Geyer of the Ruhr Regional Association (pers. comm.) for pile delineation.

2.2. Data Acquisition

2.2.1. Land Surface Temperature (LST)

To identify the surface cooling potential of all tailing piles across the Ruhr Area, Landsat 8 scenes for the Ruhr Area were selected for the periods April to September which are available from 2013 to 2023 [51]. The winter months with low sun level (October to March from 2013 to 2023) were also selected, separately. Within these selected timeslots, all relevant Landsat 8 images were filtered concerning a cloud cover of less than 5%, guaranteeing an area-wide pixel coverage. Selected dates of satellite overpassing and the corresponding images are listed and characterized also with respect to the climate conditions during satellite overpassing in Table A1 (Appendix A).

Based on the selected Landsat scenes, the LST and spectral indices were calculated according to the approach of Stumpe et al. [33]. The LST of all selected Landsat 8 OLI/TIRS images was determined based on the spectral band 10 (C2L1 dataset, TIRS 1, 10.6-11.19 μm) using the single band algorithm according to Avdan and Jovanovska [52] and Stathopoulou and Cartalis [53]. The selected Landsat images were also used to calculate spectral indices characterizing vegetation and soil properties as follows.

2.2.2. Vegetation Characteristics (NDVI, NDMI, Height)

The NDVI index was calculated as a proxy for greenness due to its strong correlation with chlorophyll content. It was calculated from the near-infrared Landsat 8 OLI band 5 (0.85-0.88 μm) and the red Landsat 8 OLI band 4 (0.64-0.67 μm) according to [54]. NDVI values range from -1 to 1, with positive values representing vegetated areas and negative values representing non-vegetated areas [54].

The NDMI index is used as proxy for vegetation moisture [55]. This index was calculated from the near-infrared 8 OLI band 5 (0.85-0.88 μm) and the short-wave infrared Landsat 8 OLI band 6 (1.57-1.65 μm). NDMI ranges between -1 and 1, with low values indicating high water stress and high values showing no water stress in the soil-leaf canopy cover level [56].

Vegetation height (VH) was determined based on a normalized digital surface model (nDSM) according to Yu et a. [57]. The corresponding DSM was derived by airborne laser scanning (ALS) with a horizontal and vertical accuracy each of 0.5*0.5 m [50]. For further analysis, the arithmetic means of the VH were calculated for each Landsat grid cell of 30*30 m.

2.2.3. Soil Characteristics (NDBaI, TVDI)

The NDBaI (Normalized Difference Bareness Index) was used in this study since the tailing piles are often characterized by a poorly vegetated and, thus, bare plateau. The index is calculated from the short-wave infrared Landsat 8 OLI band 6 (1.57-1.65 μm) and the thermal infrared Landsat 8 OLI band 10 (10.60-11.19 μm). NDBaI ranges between -1 and 1, where decreasing values indicate an increasing proportion of bare land [58].

The TVDI (Temperature Vegetation Dryness Index) can be used to assess soil moisture. It was determined based on the method proposed by Sandholt et al. [36]. This index is based on the triangular relationship between LST and NDVI values in large areas representing a wide range of surface moisture contents, from wet to dry and from bare soil to fully vegetated surfaces [36]. While the triangle's lower, nearly horizontal wet edge indicates maximal soil wetness and potential evaporation, the dry edge of the triangle represents limiting soil moisture and evaporation. Based on these relationships, the TVDI index was defined by the following equation:

TVDI = (LST - LSTmin) / (LSTmax - LSTmin)

where LST is the observed land surface temperature, and LSTmin and LSTmax are the minimum (wet edge) and maximum (dry edge) functions of NDVI, respectively. For calculation, minimum and maximum LST values of NDVI intervals with a width of 0.01 were extracted across the whole study region and separately for each Landsat summer scene. Then, separately for the extracted minimum and maximum values of the different scenes, these data were linearly fitted, and the corresponding slopes (b and d) and offsets (a and c) were used in the following equation:

LSTmin = (a + b *NDVI)

LSTmax = (c + d *NDVI)

Thus, for each pixel's observed NDVI values, a LSTmin and LSTmax value was calculated and used for TVDI calculation (equation x). TVDI values generally vary between 0 and 1, with a TVDI value of 1 indicating limited soil water content (dry edge), and a TVDI value of 0 (wet edge) indicating sufficient soil water.

To guarantee robust TVDI calculation, regression outliers were eliminated using the Cook`s Distance which summarizes how much a regression model changes when removing the i<sup>th</sup> observation. Data points from regression processing were excluded if the Cooks`s Distance was higher than 4/n. After outlier correction, we obtained validated regression models with R² generally higher than 0.96.

2.2.4. Terrain Attributes (Altitude, Slope, Aspect, Curvature)

The altitude (m ASL) of the piles was obtained from digital elevation models (DEMs) derived from airborne laser scanning (ALS) with 1*1 m horizontal resolution and vertical accuracy of 0.01 m from campaigns in 2019 and 2020 [50].

From the DEM, we determined the slope (degree), aspect (radians), and general curvature (degree) as primary terrain attributes using the SAGA Next Gen (Terrain Analysis, Morphometric) implemented in QGIS. Also, hill shading attributes were calculated for each Landsat scene according to the corresponding azimuth and altitude of light [59].

For each terrain variable, the arithmetic means within each Landsat grid cell of 30*30 m were calculated for further statistical analyses. However, since the aspect is a circular variable, the arithmetic means, and statistics that depend on the mean were inappropriate. Thus, we calculated the circular mean ($\overline{a}$) as described in the following equation:

$$\overline{a}=\arctan 2\left({\sum }_{j=1}^n\sin a_j,{\sum }_{j=1}^n\cos a_j\right)$$

(x)

where the mean direction was found by calculating the arctan2 function of the vector sum of the sinus and cosinus functions, separately. However, arctan2 values are generally varying between 0 and $\pi$ as well as between 0 and -$\pi ,$ representing, in our case east and west exposition of the tailing piles, respectively. However, since we assumed that according to the azimuth of light north and south expositions of the tailing piles were more important for their LST pattern, the absolute values (ABS) of the circular means were also calculated for Random Forest Regression (RFR) modeling. Values from 0 to ${\displaystyle \frac{\pi }{2}}$ and from ${\displaystyle \frac{\pi }{2}}$to $\pi$ then expressed the north and south exposition of the tailing piles, respectively. However, we could no longer distinguish between the east and west expositions.

2.3. Statistical Analysis

To evaluate the importance of independent factors on the LST pattern of the piles, cluster analyses and random forest regression (RFR) were used according to the approach of Stumpe et al. [14]. Thus, we followed the same statistical procedure, as briefly described below.

2.3.1. Cluster Analyses

In the first step, we typified the tailing piles based on their mean LST. According to the different selected Landsat summer scenes, we used each tailing pile's corresponding mean LST values for k-mean classification.

The conventional k-mean algorithm defines k centroids at a rate of one for each cluster so that an object is more similar to objects in its cluster than to objects in a different cluster. The algorithm aims to minimize the sum of square error in the Euclidean distance over all data points [60,61,62]. The optimal number of clusters was identified using the Silhouette method as described in Kodinariya and Makwana [61]. Generally, this method compares within-cluster distances with between-cluster distances, where the greatest difference between clusters indicates the best number of clusters. Analysis was conducted in R version 4.3.3 [63] using the Factoextra package [64].

The average silhouette score and the mean cluster silhouette scores were used to evaluate cluster accuracy. Silhouette scores generally range from -1 to 1, where a threshold is commonly set at 0.5, especially for the average silhouette score. That means scores of more than 0.5 indicate high-quality cluster assignments.

For a final TP typification, the resulting clusters were characterized according to their total tailing pile area, spectral indices, and terrain attributes, which were averaged for each cluster.

2.3.2. Random Forest Regression (RFR)

In the next step, random forest regressions were performed to identify different factors controlling the LST pattern of the tailing piles pixel-wise (30*30 m). We also tested multiple linear regression (MLR), but they failed in prediction accuracy, with R² being lower than those of the random forest models.

We used the LST distribution of the tailing piles as the dependent variable for RFR modeling, while separate RFR models were run to explain the LST distribution of the different Landsat summer scenes. The separate modeling results were finally averaged. As independent and explaining variables, we used various vegetation (NDVI, NDMI, and VH) and TP terrain characteristics (altitude, slope, curvature, aspect (ABS circular), and hill shading values, while hill shading was calculated separately for each Landsat scenes in dependence of sun angle and azimuth.

The RFR analyses generally belong to the more complex supervised learning algorithms [65,66] and are based on several random decision trees fitted to a dataset with a mean prediction of the ensemble as output. Each tree is grown from a random subset of the features, whereas large numbers of trees run in parallel with no interaction, making the RFR often powerful and accurate [14,66]. However, the number of trees must be optimized to obtain the most stable performance of the RF models. Analysis was conducted in R [63]using the RandomForest package [67].

2.3.3. Feature Importance Analysis

To identify the importance of the independent variables, variables are removed from RFR analyses, assessing how the prediction changes [68]. As indicators, on the one hand, the mean decrease accuracy (%IncMSE) was used, which is a measure of how much the model’s accuracy decreases when a particular variable is randomly permuted. A higher %IncMSE value represents higher variable importance [69]. On the other hand, the mean decrease Gini (IncNodePurity) was used, representing the total decrease in node impurity from splitting on a predictor in the tree construction process averaged over all trees [70]. Both indicators do not differentiate between positive or negative effects on the LST variance. Analysis was conducted in R [63] using the Metrics package [71].

2.3.4. Model Validation

In the first step, each dataset was randomly divided into train (70%) and test (30%) data to describe model performance. In the second step, regression models were fitted using the training datasets so that, in the third step, the dependent variable (LST) was predicted using the trained models and the corresponding test datasets. Finally, the predicted and the observed variables were used to calculate error metrics, such as the coefficient of determination R2, the root mean squared error (RMSE), and the mean absolute error (MAE).

3. Results

3.1. Tailing Pile Characteristics

3.1.1. Spatial Distribution across the Ruhr Area

In the first step, tailing piles were characterized by their spatial distribution (Figure 1) and area expansion (Table 1) across the Ruhr Area to assess their role as potential cooling islands within the urban green infrastructure.

The tailing piles are distributed across the main settlements in the Ruhr Area, as indicated by high Normalized Build-up Index (NDBI) values, so they can generally be characterized as urban greens. Nevertheless, larger tailing piles often are found more adjacent to periurban areas indicated by lower NDBI values, while the smaller ones are often surrounded by urban areas represented by high NDBI values.

With 82 tailing piles, their number is much smaller than the more common conventional urban greens such as parks, allotment gardens, or cemeteries (Table 1). Nevertheless, tailing piles comprise large areas with a mean size of 0.323 km^2, which is about ten times larger than the mean size of the other UGSs. Consequently, the piles cover 0.51 % of the total Ruhr Area, which is close to the percentages of parks, allotments, and cemeteries, with 1.46, 0.79 and 0.69%, respectively (Table 1).

3.1.2. Morphological Attributes

Generally, coal tailing piles are singular hills with elevations reaching 50 to 200 m ASL (Figure 7A) in the rather flat Ruhr area landscape of predominantly 20-60 m ASL. The single hill character is well reflected in the homogenous distribution of the pile aspect values across all directions (Figure 2B). The slope of the pile hillsides varies between 0 and about 40° with a mean of about 14° (Figure 2A). The relatively high frequency of slope values leveling around 0° reflects the hilltop plateaus and the typical step relief stemming from the pile construction. The nDSM histogram in Figure 2C reflects the distribution of plant heights, derived from the aerial laserscan measurements [50]. The predominance of low plant heights (0-0.5 m) shows that large areas of the piles are either free of vegetation or covered by grassland. Often, tailing plateaus are deliberately left bare to illustrate their industrial origin and numerous pathways were constructed for recreational access, also contributing to this nDSM class. The largest part of the piles is covered by shrub and tree vegetation reaching heigths of up to 10 m, with only few areas containing higher trees, reflecting the fact that most piles were only rehabilitated in the last 2-4 decades.

3.1.3. Thermal Footprints

Table 1 summarizes the mean LST values of all tailing piles and the other UGSs for the analysed five different satellite summer scenes. Clearly, the mean LSTs of the piles are similar to those of the other UGSs. Interestingly, the standard deviation of the LST means was generally distinctly higher for the piles than for the other UGS, indicating a higher heterogeneity within this UGS class. The LST differences between the UGS classes could not be tested by ANOVA due to the different LST variances.

When looking at the frequency distribution of mean summer LSTs on the piles at the pixel level, the histogram shows a pronounced left-skewed normal distribution (Figure 3A). This underlines the very high LST variability among and within the piles. At the same time, the histogram shows that the high LST variance is not only based on a few extreme values but forms a robust distribution of pile surface properties on a small scale. One factor responsible for this may be exposition, since at the time of satellite overpasing around 11 AM, south facing slopes should be more exposed to radiation than north facing slopes. This was confirmed, when the pixels were differentiated according to their position on differently exposed slopes (Figure 3B).

Clearly, south and south-east facing slopes were the warmest and the northern to western slopes were the coolest. However, the LST differences between the slope expositions are not significant and amount to only 1.2°C and thus only partly explain the large LST variability shown in Figure 3A.

3.2. Thermal Typification of Tailing Piles

Based on the mean LST values of the different tailing piles, a cluster analysis was performed, and the analysis included the mean LST values of each satellite summer scene separately. After clustering, we further aimed to characterize warmer or colder pile clusters, considering other properties describing pile vegetation status or topographical parameters for a final pile categorization. Therefore, we hoped to identify different characteristics between the colder and warmer pile classes and be further able to provide recommendations for pile management.

For valid cluster results, we used the within-cluster sum of squares in the first step according to the elbow methods defining the optimal number of four clusters (Figure 4A). Based on this, the k-mean algorithm was used for LST mean clustering. Figure 4 and Table 1 show the cluster assignments and their quality, respectively. Cluster A includes the tailing piles with the highest LST values, which are described by a mean LST of 28.17 °C (Table 2) followed by clusters B, C, and D with a descending order of LST means with 26.05, 24,31, and 22.88°C, respectively. Although the clusters showed a clear separation by mean LSTs, the corresponding standard deviations and mean cluster silhouette scores vary between the clusters, indicating a different quality of cluster assignment. Cluster A contains only seven tailing piles and showed the highest standard deviation and the lowest mean cluster silhouette score of 0.22. Clusters C and D enclose the most tailing piles and thus are more valid clusters with the lowest standard deviation of LST means and the highest mean cluster silhouette scores of 0.33 and 0.40, respectively. Nevertheless, the overall silhouette score of the clustering was calculated at 0.68, indicating a general high-quality cluster assignment that enabled further cluster characterization.

Mean terrain attributes such as the mean tailing pile altitude (DEM) or the mean slope and curvature showed no significant differentiation between the four clusters (ANOVA analyses, p-values>0.05). The mean tailing pile area seems to be related to the cluster LST characteristics since the mean area sizes decline with the mean cluster temperatures, but nevertheless, the mean area sizes between the clusters showed no significant differences (ANOVA analyses, p-values>0.05). However, the spectral indices (NDVI, NDMI) and the VH showed significant differences between the four clusters, with increasing mean index values being related to decreasing mean cluster LSTs. Soil moisture content (TVDI) also differed significantly between the four clusters, indicating more water limitations in the hottest clusters A and B and least in the cool cluster D. With values ranging from -0.558 to -0.536, the soil bareness index (NDBaI) was quite similar among the four clusters, although differences were partly significant. In summary, this shows that vegetation parameters and soil moisture are closely related to the differentiation between the four pile LST types.

For visualization of these results, Figure 5, shows one arbitrarily selected tailing pile from each cluster as aerial images (A1-D1), DEMs (A2-D2) and as LST patterns (A3-D3). As indicated by the aerial images and the DEMs, the tailing piles generally form hills higher than the surrounding with sloping sides. However, concerning the pile typification, the aerial images also show the distribution of green areas and bare soils. Clearly, the tailing pile “Kohlenhuck” (Figure 5A-1) belonging to the hottest cluster is characterized by a high proportion of bare soil (Figure 5A-3), while the cool piles “Tetraeder” and “Lohberg Nord” (Figure 5C,D) are dominated by tree vegetation.

3.3. Controlling Factors of the LST Distribution of the Tailing Piles

Based on the cluster results, we focused on factors controlling the LST distribution on the tailing piles using a detailed pixel-based analysis. Therefore, in the first step, we run RFR models with the LST as dependent and various independent variables shown in Fig. 6, with models being run separately for the five satellite summer scenes. The model performance was then evaluated based on the mean values from the five models (Table 3). Generally, we obtained high RFR model performances with a mean R² of 0.85 and a RMSE of 0.39 (Table 3). In order to evaluate the importance of the independent variables across all satellite summer scenes, the %IncNodeMSE values are visualized as mean values in Figure 6A. In accordance with the cluster characterization, the LST variance of the piles was found to be mainly controlled by vegetation properties, which is clearly indicated by high %IncNodeMSE values of the NDVI, the NDMI, and the VH. Although the tailing piles form hilly landscapes, terrain attributes such as pile altitude, slope, curvature, aspect or the hillshade showed relatively low %IncNodeMSE values, thus being less relevant for the LST distribution on the piles. As a control, we also analyzed these dependencies for one Landsat winter scene (Figure 6B), resulting in a less accurate RFR model with a R² of only 0.75 (Table 3). However, in winter with low sun altitude, the hillshade aspect dominates the LST distribution of the piles, especially the aspect. Due to the dominance of annual herbs and deciduous trees and shrubs, vegetation is inactive in winter and thus is of minor importance for LST variability during this season.

Since shading effects were found to be relevant for the spatial LST pattern in winter, we assumed that they may also matter in summer at certain pile altitudes. Most of the piles reach altitudes up to 125 m, with a maximum frequency between 75 and 100 m ASL (Figure 7A) and only few being higher than 125 m. In order to test for the importance of shading, RFR models for the Landsat summer scenes were run on datasets filtered depending on the maximal pile altitudes in 25 m intervals (Figure 7B). Interestingly, the importance of the aspect, as expressed in the %IncNodeMSE value increased with maximum pile elevation. Due to the decreasing sizes of the datasets, the standard deviation of the %IncNodeMSE values also increases with pile elevation (Figure 7B). Nevertheless, at a maximal pile height of 125 m ASL a threshold was observed where the %IncNodeMSE significantly increased.

Figure 7. Histogram of the tailing pile heights (A) and the variable importance (IncNodeMSE) of the aspect in RFR models differentiated according to pile heights (B). The dotted line indicates the threshold value of the pile height above which a distinct increase of the aspect importance in the RFR models occurs.

Figure 7. Histogram of the tailing pile heights (A) and the variable importance (IncNodeMSE) of the aspect in RFR models differentiated according to pile heights (B). The dotted line indicates the threshold value of the pile height above which a distinct increase of the aspect importance in the RFR models occurs.

Consequently, in the next step, the total dataset was divided into two datasets covering piles with maximal pile elevation above and below 125 m ASL. Both datasets were then subjected to the RFR analyses, with the results for both analyses shown in Table 3 and Figure 7. With respect to the Landsat summer scenes, the RFR model performance showed no considerable differences between the three data sets (all pile pixel data, pixel data of piles with DGMmax< 125m ASL or with DGMmax> 125m ASL), but with respect to the Landsat winter scenes model performance improved using pile data with DGMmax > 125m. With an R² of 0.83, the LST variance of the pixel LST values of piles with DGMmax > 125m was up to 8% better explained by the independent variables than using the unfiltered pile dataset.

Concerning the variable importance, a generally more differentiated pattern of the %IncNodeMSE values across the independent variables was observed using pixel data of piles with DGMmax< 125m ASL and with DGMmax> 125m ASL, separately. For the Landsat summer scenes, the importance of the aspect and the VH increases when pixel data of piles with DGMmax> 125m ASL was used whereas simultaneously, the importance of the NDVI, NDMI, and NDBal decreases. The effects of the dataset separation on the variable importance were similar for the Landsat winter scene but less pronounced.

In summary, it is concluded that a data set separation concerning the maximal pile height resulted in more precise and more differentiated RFR modelling results carving out the importance of the pile aspect and slightly modifying the importance of vegetation properties. Nevertheless, vegetation characteristics were found to be the dominant factors controlling the summer LST pile variability.

3.4. Impact of Soil Moisture on the Pile LST Pattern

Since vegetation characteristics dominate the pile LST pattern, soil moisture as one important controlling factor for transpirational cooling and vegetation health was analyzed using the TVDI index. Since the TVDI calculation is based on the relationship between the LST and the NDVI, this index was excluded as an independent variable in the RFR models. Instead, we analyzed the spatial distribution of the TVDI related to those of the NDMI and VH, assuming that soil moisture is an important indirect factor controlling the LST pattern of the piles.

To identify more general relationships, mean TVDI pixel values from all Landsat summer scenes were analysed. To display their spatial distribution across all piles, TVDI and VH pixel mean values were grouped according to the eight main geographic directions (Figure 8).

Across all piles, mean TVDI values varied between 0.30 and 0.46, where the lowest values (moist soils) are found on the north and northwest slopes of the pile hills (Figure 8). Complementary to this, the highest values and thus soils with the lowest moisture were found in the south and southeast exposition of the piles. Interestingly, the spatial distribution of the mean TVDI seems to be more than a temporal snapshot since the mean TVDI directional values significantly correlated with the mean VH directional values (r = -0.8), which is also visualized in in Figure 8. Here, the highest vegetation heights are found on the northwest and west while the lowest occur on the south and southeast facing slopes. Apparently, soil moisture as a growth limiting factor is another indirect but important variable controlling the LST patterns on tailing piles.

4. Discussion

4.1. The Role of Tailing Piles as Cooling Urban Greens

Generally, tailing piles showed mean LSTs in a similar range as those of parks or allotment gardens across the Ruhr Area. While the role of these UGSs as cooling islands in UHI mitigation has been demonstrated in numerous studies, our study shows that rehabilitated tailing piles can be of similar significance. This is a more or less neglected phenomenon, resulting in a lack of comparable studies. To our knowledge, only the case study of Glocke et al. [45] analyzed the cooling potential of a tailing pile, also in the Ruhr Area. They showed that compared to flat urban parks such piles have an additional cooling potential during nighttime, because then cool air can drain down the slopes into adjacent residential areas. This highlights an aspect, which we couldn’t consider but which can additionally contribute to our documented potential of tailing piles to cool their immediate surroundings and thus provide an important ecosystem service.

The unique topography of tailing piles as hills in a rather flat surrounding is also partly responsible for their rather high variability in LST values compared to other UGSs. The influence of terrain factors on LST has also been studied by Mathew [72] and Tan et al. [34] [34,72], who generally explained this by terrain-dependent variation of solar incident radiation.

4.2. Factors Controlling LST on Tailing Piles

Cluster analyses based on the mean pile LST values showed that vegetation characteristics mainly differentiate the piles into four clusters. Interestingly, terrain attributes of the piles showed no statistically significant differences between the LST clusters. Similarly, the variability of mean pixel LST values across all piles was mainly explained by vegetation vitality indices (NDVI, NDMI) and vegetation height and far less by terrain factors in the RFR models, which reached up to 86% prediction accuracy. Similar strong relationships between LST and NDVI or NDMI values have been reported in numerous other studies because vegetation abundance and its moisture condition directly regulate the evaporative cooling processes [33]. Since the complex tailing pile environment shows a high heterogeneity of grass, shrub, and tree vegetation, both indices vary greatly and thus have a high impact on LST patterns. The fact that NDMI was more important than NDVI was also previously reported in several studies [33,34,35,73] and explained by the importance of plant available water as the limiting factor for transpirative cooling processes.

The influence of vegetation height on LST has been analyzed less often [74,75], most likely due to lack of available data from the expensive aerial laser scan campaigns. In a study by Chen et al. [76], the species composition of the tree canopy showed the highest impact on LST variability, but its vertical structure was also an important impact factor. Similar to our results, Gage and Cooper [77] found that NDVI played the most crucial role in explaining LST variability, but this was directly followed by tree height. They explained the additional impact of vegetation height on the LST pattern by shading effects on adjacent areas even when plant transpiration rates are low.

Contrary to our expectations, terrain attributes had only minor impacts on the summer LST patterns, although tailing piles generally have pronounced 3D hill shapes. The exposure to solar radiation and the associated effect on LST depends on elevation, aspect and steepness of slopes [41,78,79]. Due to the relatively small elevation differences of 50-150 m in our study area, this factor has little influence on LST. In our region where Landsat overpasses around early noon time, slope steepness in the range of >0° to 45° is positively related to the amount of direct solar radiation on south facing slopes, whereas the opposite is true for north facing slopes [78]. This opposing effect therefore likely masks the potential influence of slope in our regression analyses. While the impact of the aspect on the LST distribution was also negligible among all tailing piles, it increased substantially when only piles with elevations > 125 m were considered. At the relatively high solar azimuth angels of 60 to 63° for our summer satellite images, the incidence of sunlight was apparently impeded on north-exposed hillsides only of these higher tailing. For the winter image, aspect and the hillshade parameter modeled by ArcGIS became the most important controlling variables for LST because vegetation was bare or dormant, but also because the solar zenith angle drops to about 15° so that shading effects on slopes that are not exposed to the sun become more prominent, thus becoming an important factor for LST variability.

Studies like that of Bai et al. [41] or Karbalee et al. [40] which identified terrain attributes such as elevation, slope and aspect to control LST variability were generally conducted in natural mountainous areas where elevation differences were much greater and which are then associated with changes in the composition and density of the vegetation. The tailing piles are only low artificial hills with a man-made plant cover that was varied mainly for aesthetic reasons and is therefore not related to any terrain attributes.

Although the aspect parameter generally showed a low impact on the LST patterns, it was a relevant variable for describing the relationship between soil moisture estimated by the TVDI index and vegetation height. We observed that slopes with north-northwest exposition generally showed the lowest TVDIs i.e., highest soil moisture values, and accordingly, the greatest vegetation heights. Complementary, south, and southeast exposed hillsides generally showed low soil moisture, which most likely caused the lower vegetation heights. The effect of soil moisture variation on vegetation growth is well-known and is often a result of topographic variation [80,81,82]. Thus, studies frequently used the Topographic Wetness Index (TWI) to describe and predict vegetation height. For example, in the study of Mohamedou et al. [81] the TWI was used to improve the growth-prediction accuracy of LiDAR-based vegetation heights. However, to the best of our knowledge, the TVDI has never been used in this context. Unlike the TWI, the TVDI does not consider water movement in soils, but it appears to be a promising factor in predicting vegetation height, especially on tailing piles. Since the piles were constructed from fragmented bedrock debris, seepage water will most likely not move laterally downhill but will drain into the pile body, making it unavailable for plant roots, independent of slope position.

4.3. Implications for Urban Planning

Generally, the LST footprint of the tailing piles was found to be similar to those of other urban green spaces such as parks, cemeteries or allotments. Parks have been well characterized as effective cooling islands in numerous studies, so even for city residents, an awareness of their health-related benefits has been created [7]. Although the tailing piles in the Ruhr Area are used for local recreation, city residents are not aware of their cooling potential and their potential to improve healthy living conditions, especially during summertime. Thus, we advise urban planners to promote the remarkable benefit of the tailing piles as cooling islands to increase visitor numbers and improve the quality of urban life.

Within the tailing pile areas, the presence of taller vegetation and its moisture content were found to be the most important factors controlling LST variability during summer. Along south and southeast facing slopes, soil water availability limits tree growth and transpirational cooling. Increasing soil water capacity would thus increase the cooling potential of these expositions. This could be achieved during tailing pile rehabilitation by increasing the soil cover thickness or by selecting soil substrates with high soil water retention for these sites. On established vegetated piles, soil additives that provide high water capacities, such as organic mulches or biochar [83,84] could be spread out on the surfaces or mixed into the topsoil to achieve improvements of the cooling potential.

Another cooling feature of tailing piles that is not shared by other more flat urban green spaces is the downhill drainage of cool air into the surrounding neighborhoods during night, as described by Glocke et al. [45]. According to our data, this should be most pronounced along north and northwest facing slopes where daytime temperatures are lowest. Since this flux of cool air can easily be blocked by obstacles, such as walls or buildings at the hillfoot (Glocke et al. 2023), we advise urban planners to improve the pile cooling potential by removing barriers in the drainage pathways, especially in north and northwest expositions, to alleviate night-time urban heat in the surrounding neighborhoods.

4.4. Limitations and Open Questions

Many parts of this analysis are based on mean values from dates with warm to hot air temperatures, low cloud cover and no rainfall during the previous 2 weeks. This assures that the weather conditions were comparable across the whole study area before and during satellite imaging. However, the variability of LST within each UGS class und between the four UGS classes differed from date to date, indicating that season and other dynamic parameters such as soil moisture or plant phenology, possibly also management affect the cooling potential differently in the UGSs. Working with mean values allows for a better classification of the pile types, but masks temporal effects that may be relevant for LST variability at the local scale. On the other hand, the RFR models explaining the LST variability at the pixel level which are also based on mean values from the summer dates show a very good model performance, with R² > 0.8 and low RMSE and MAE values. This indicates that temporal effects on the cooling potential of the tailing piles may be of minor importance. But this remains to be investigated.

Another flaw in the data analysis is based on the temporal incongruency between satellite scenes and the DEM and nDSM data. While the Landsat scenes stem from 2013 to 2023, the LiDAR campaigns were in 2019 and 2020. Since the rehabilitation of most tailing piles was completed before 2013 [85], the tree heights determined in 2019/20 may not be representative for the height distributions three to six years before this date. However, this was the best data available for the region and once again, the model performance based on the mean values from the satellite scenes indicate that the found relationships are solid and meaningful.

Finally, the LSTs of the tailing piles may not only originate from solar radiation but also from internal sources since several piles in the region contain patches of smoldering coal remains, leading to localized internal temperatures of 60 to 400°C (Kürten et al., 2014; Bezirksregierung Arnsberg, pers. comm.). But the effect of these smoldering patches on surface temperatures or their location and spatial extent are unknown and could only be estimated from night-time aerial or satellite imaging, which both are currently not available at the required spatial resolution. In any case, such “hot spots” can only have very minor effects on surface LST values since these are very well predicted from soil and vegetation parameters in our RFR models.